Triangular Metrics

automesh metrics tri --help
Quality metrics for an all-triangular finite element mesh

Usage: automesh metrics tri [OPTIONS] --input <FILE> --output <FILE>

Options:
  -i, --input <FILE>   Mesh input file (exo | inp | stl)
  -o, --output <FILE>  Quality metrics output file (csv | npy)
  -q, --quiet          Pass to quiet the terminal output
  -h, --help           Print help

automesh implements the following triangular element quality metrics:

• Maximum edge ratio ${\mathrm{ER}}_{\max }$
• Minium scaled Jacobian ${\mathrm{SJ}}_{\min }$
• Maximum skew
• Element area
• Minimum angle ${\theta }_{\min }$

A brief description of each metric follows.

Maximum Edge Ratio

• ${\mathrm{ER}}_{\max }$ measures the ratio of the longest edge to the shortest edge in a mesh element.
• A ratio of 1.0 indicates perfect element quality, whereas a very large ratio indicates bad element quality.
• Knupp et al.¹ (page 26) indicate an acceptable range of [1.0, 1.3].

Minimum Scaled Jacobian

• ${\mathrm{SJ}}_{\min }$ evaluates the determinant of the Jacobian matrix at each of the corners nodes, normalized by the corresponding edge lengths, and returns the minimum value of those evaluations.
• Knupp et al.¹ (page 29) indicate an acceptable range of [0.5, 2*sqrt(3)/3] $\approx$ [0.5, 1.2].
• A scaled Jacobian close to 0 indicates that the triangle is poorly shaped (e.g., very thin or degenerate), which can lead to numerical instability.
• A scaled Jacobian of 1 indicates that the triangle is equilateral, which is the ideal shape for numerical methods. This is achieved through scaling as follows:

$$\mathrm{S}{\mathrm{J}}_{\min }=\frac{\sin ({\theta }_{\min })}{\sin (60^{\circ })}$$

• In the preceeding equation,
  • $60^{\circ }$ is the smallest angle of an equilateral triangle, and
  • ${\theta }_{\min }$ is the smallest angle of the subject triangle.

Maximum Skew

• Skew measures how much an element deviates from being a regular shape (e.g., in 3D a cube; in 2D a square or equilateral triangle). A skew value of 0 indicates a perfectly regular shape, while higher values indicate increasing levels of distortion.
• Knupp et al.¹ does not give a definition of skew for triangles, so we provide our definition below. For a triangle where ${\theta }_{\min }$ is the smallest angle of the triangle,

$$\mathrm{ske}{\mathrm{w}}_{\max }=\frac{60^{\circ }-{\theta }_{\min }}{60^{\circ }}$$

• For an equilateral triangle, ${\theta }_{\min }=60^{\circ }$ and $\mathrm{ske}{\mathrm{w}}_{\max }=0$.
• In the limit as ${\theta }_{\min }\to 0^{\circ }$ $\mathrm{ske}{\mathrm{w}}_{\max }\to 1$.

Element Area

• Measures the area of the element.

Minimum Angle ${\theta }_{\min }$

• The smallest the three angles of a triangle.
• Ideal value: $60^{\circ }$ (for an equilateral triangle).
• The ${\theta }_{\min }$ should be maximized (kept as close to $60^{\circ }$ as possible) and above a minimum threshold (e.g., $15^{\circ }$ to $30^{\circ }$).
• Acute angles (close to $0^{\circ }$) cause high errors and numerical instability.

Unit Tests

We use the ABAQUS input file single_valence_04_noise2.inp. We import the file into Cubit and create a triangular surface mesh:

import abaqus mesh geometry  "/Users/chovey/autotwin/automesh/tests/input/single_valence_04_noise2.inp" feature_angle 135.00
surface 1 scheme trimesh minimum size 100
surface 2 scheme trimesh minimum size 100
surface 3 scheme trimesh minimum size 100
surface 4 scheme trimesh minimum size 100
surface 5 scheme trimesh minimum size 100
# surface 6 scheme trimesh minimum size 100 # there is no side 6, two sides were merged
delete mesh surface all propagate
surface all scheme trimesh
mesh surface all
quality tri all aspect ratio global draw mesh list detail
quality tri all scaled jacobian global draw mesh list detail
quality tri all element area global draw mesh list detail
export stl ascii "/Users/chovey/autotwin/automesh/tests/input/single_valence_04_noise2.stl" mesh  overwrite

We also examine several .stl files and process them, for example,

import stl "/Users/chovey/autotwin/automesh/tests/input/one_facet.stl" feature_angle 135.00 merge make_elements
surface 1 scheme trimesh minimum size 100
delete mesh surface 1  propagate
surface 1  scheme trimesh
mesh surface 1
quality tri all aspect ratio global draw mesh list detail
quality tri all scaled jacobian global draw mesh list detail
quality tri all element area global draw mesh list detail

We verify the following element qualities:

Figure: Triangle metrics. Leading values are from automesh. Values in [brackets] are from an independent Python calculation, (see metrics_triangle.py) and agree with automesh in double precision with a tolerance of less than 2.22e-15. Except for edge ratio, all values were also verified with Cubit. Cubit uses the term Aspect Ratio; it is not the same as Edge Ratio.

• File A is single_valence_04_noise2.inp.
• File B is one_facet.stl.
• File C is an equilateral triangle with nodal coordinates at (-2, 0, 0), (2, 0, 0), and (0, 2*sqrt(3), 0) and has side length 4.0, saved to tests/input/equilateral_4.stl.
• File D is an equilateral triangle with nodal coordinates at (-0.5, 0, 0), (0.5, 0, 0), and (0, sqrt(3) / 2, 0) and has side length 1.0, saved to tests/input/equilateral_1.stl.
• File E is an off axis triangle with approximate (30, 60, 90) degree inner angles, with coordinates at (0.0, 1.0, 3.0), (2.0, 0.0, 2.0), and (1.0, sqrt(3.0) + 1.0, 1.0), saved to tests/input/off_axis.stl.
• e is the element number in the mesh.

Local Numbering Scheme

Nodes

The local numbering scheme for nodes of a triangular element:

        2
       / \
      /   \
     /     \
    /       \
   0---------1

Faces

The local numbering scheme for faces of a triangular element:

Source

metrics_triangle.py

"""This script is a quality control tool for the metrics of a tri mesh."""

from typing import Final

import numpy as np

DEG_TO_RAD: Final[float] = np.pi / 180.0
J_EQUILATERAL: Final[float] = np.sqrt(3.0) / 2.0  # sin(60 deg)


def nf(x):
    """Does a list comprehension, casting each item from a np.float64
    to a float."""
    return [float(y) for y in x]


nodal_coordinates = (
    (-0.2, 1.2, -0.1),  # single_valence_04_noise2.inp begin
    (1.180501, 0.39199, 0.3254445),
    (0.1, 0.2, 0.3),
    (-0.001, -0.021, 1.002),
    (1.2, -0.1, 1.1),
    (1.03, 1.102, -0.25),
    (0.0, 1.0, 1.0),
    (1.01, 1.02, 1.03),  # single_valence_04_noise2.inp end
    (0.0, 0.0, 1.0),  # one_facet.stl begin
    (0.0, 0.0, 0.0),
    (1.0, 0.0, 0.0),  # one_facet.stl end
    (-2.0, 0.0, 0.0),  # equilateral with edge length 4.0 start
    (2.0, 0.0, 0.0),
    (0.0, 2.0 * np.sqrt(3.0), 0.0),  # equilateral with edge length 4.0 end
    (-0.5, 0.0, 0.0),  # equilateral with edge length 1.0 start
    (0.5, 0.0, 0.0),
    (0.0, np.sqrt(3.0) / 2.0, 0.0),  # equilateral with edge length 1.0 end
    (0.0, 1.0, 3.0),  # off_axis.stl begin
    (2.0, 0.0, 2.0),
    (1.0, np.sqrt(3.0) + 1.0, 1.0),  # off_axis.stl end
)

element_node_connectivity = (
    (1, 2, 3),  # single_valence_04_noise2.inp begin
    (4, 2, 5),
    (1, 6, 2),
    (4, 3, 2),
    (4, 1, 3),
    (4, 7, 1),
    (2, 8, 5),
    (6, 8, 2),
    (7, 8, 6),
    (1, 7, 6),
    (4, 5, 8),
    (7, 4, 8),  # single_valence_04_noise2.inp end
    (9, 10, 11),  # one_facet.stl
    (12, 13, 14),  # equilateral triangle with side length 4.0
    (15, 16, 17),  # equilateral triangle with side length 1.0
    (18, 19, 20),  # off_axis.stl
)

NODE_NUMBERING_OFFSET: Final[int] = 1

mesh_element_max_edge_lengths = []
mesh_element_edge_ratios = []
mesh_element_minimum_angles = []
mesh_element_maximum_skews = []
mesh_element_areas = []
mesh_element_jacobians = []
mesh_element_jacobians_scaled = []


def angle(a: np.ndarray, b: np.ndarray) -> float:
    """Given two vectors, find the angle between them."""
    dot_product = np.dot(a, b)
    norm_a = np.linalg.norm(a)
    norm_b = np.linalg.norm(b)

    cos_theta = dot_product / (norm_a * norm_b)

    angle_radians = np.arccos(cos_theta)
    angle_degees = np.degrees(angle_radians)

    return angle_degees


for element in element_node_connectivity:
    print(f"element with nodes: {element}")
    path = element + (element[0],)
    # print(f"  node path {path}")
    pairs = tuple(zip(element, element[1:] + (element[0],)))
    print(f"  node pairs {pairs}")
    element_edge_ratios = []
    element_minimum_angles = []
    element_minimum_jacobians = []
    edge_vectors = ()
    # edge ratios
    for pair in pairs:
        print(f"    pair {pair}")
        aa, bb = pair
        edge = np.array(nodal_coordinates[bb - NODE_NUMBERING_OFFSET]) - np.array(
            nodal_coordinates[aa - NODE_NUMBERING_OFFSET]
        )
        edge_vectors = edge_vectors + (edge,)
        edge_length = np.linalg.norm(edge)
        # print(f"    lens {edge_length}")
        element_edge_ratios.append(edge_length)

    # print(f"  edge ratios {element_edge_ratios}")

    # edge ratios
    len_max = max(element_edge_ratios)
    # print(f"  max edge ratio {len_max}")
    mesh_element_max_edge_lengths.append(len_max)

    len_min = min(element_edge_ratios)
    # print(f"  min edge ratio {len_min}")
    ratio = len_max / len_min
    mesh_element_edge_ratios.append(ratio)

    # edge vectors and then angles
    edge_vectors_pairs = tuple(zip(edge_vectors, edge_vectors[1:] + (edge_vectors[0],)))
    # print(f"  edge vectors pairs {edge_vectors_pairs}")

    for item in edge_vectors_pairs:
        # print(f"    edge vectors pair {item}")
        # flip the direction of the first vector so that it shares an origin
        # with the secon vector
        angle_degrees = angle(-1.0 * item[0], item[1])
        # print(f"    angle {angle_degrees}")
        element_minimum_angles.append(angle_degrees)

    print(f"  element angles (deg) {nf(element_minimum_angles)}")
    angle_min = min(element_minimum_angles)
    print(f"  min angle (deg) {angle_min}")
    mesh_element_minimum_angles.append(angle_min)

    jacobian_e = np.sin(angle_min * DEG_TO_RAD)
    jacobian_scaled_e = jacobian_e / J_EQUILATERAL
    print(f"  min Jacobian (sin(angle_min)) {jacobian_e}")
    print(f"  min scaled Jacobian  {jacobian_scaled_e}")
    mesh_element_jacobians.append(jacobian_e)
    mesh_element_jacobians_scaled.append(jacobian_scaled_e)

    skew_max = (60.0 - angle_min) / 60.0
    mesh_element_maximum_skews.append(skew_max)

    # Compute areas only for triangles for now.
    if len(element) == 3:
        # area of a triangle
        aa = np.linalg.norm(edge_vectors[0])
        bb = np.linalg.norm(edge_vectors[1])
        cc = np.linalg.norm(edge_vectors[2])
        # Calculate the semi-perimeter
        ss = (aa + bb + cc) / 2.0
        # Use Heron's formula to calcuate the area
        area = np.sqrt(ss * (ss - aa) * (ss - bb) * (ss - cc))
        mesh_element_areas.append(area)

print(f"\nmesh element max edge lengths: {nf(mesh_element_max_edge_lengths)}")
print(f"\nmesh element edge ratios: {nf(mesh_element_edge_ratios)}")
print(f"\nmesh element minimum angles: {nf(mesh_element_minimum_angles)}")
print(f"\nmesh element maximum skews: {nf(mesh_element_maximum_skews)}")
print(f"\nmesh element areas: {nf(mesh_element_areas)}")
print(f"\nmesh minimum scaled Jacobians: {nf(mesh_element_jacobians_scaled)}")

References

1. Knupp PM, Ernst CD, Thompson DC, Stimpson CJ, Pebay PP. The verdict geometric quality library. SAND2007-1751. Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA (United States); 2006 Mar 1. link ↩ ↩2 ↩3