Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

Metrics

automesh metrics --help
Quality metrics for an existing finite element mesh

Usage: automesh metrics <COMMAND>

Commands:
  hex   Quality metrics for an all-hexahedral finite element mesh
  tet   Quality metrics for an all-tetrahedral finite element mesh
  tri   Quality metrics for an all-triangular finite element mesh
  help  Print this message or the help of the given subcommand(s)

Options:
  -h, --help  Print help

Hexahedral Metrics

automesh metrics hex --help
Quality metrics for an all-hexahedral finite element mesh

Usage: automesh metrics hex [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 hexahedral element quality metrics defined in the Verdict report.[^Knupp_2006]

* Maximum edge ratio ${\rm ER}_{\max}$
* Minium scaled Jacobian ${\rm SJ}_{\min}$
* Maximum skew
* Element volume

A brief description of each metric follows.

### Maximum Edge Ratio

* ${\rm 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.*[^Knupp_2006] (page 87) indicate an acceptable range of `[1.0, 1.3]`.

### Minimum Scaled Jacobian

* ${\rm 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.*[^Knupp_2006] (page 92) indicate an acceptable range of `[0.5, 1.0]`, though in practice, minimum values as low as `0.2` and `0.3` are often used.

![](img/metrics_msj.png)

Figure. Illustrate of minimum scaled Jacobian[^Hovey_2023] with acceptable range for quality occurring in `[0.3, 1.0]`.

### 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.*[^Knupp_2006] (page 97) indicate an acceptable range of `[0.0, 0.5]`.

### Element Volume

* Measures the volume of the element.

## Hexahedral Unit Tests

Inspired by Figure 2 of Livesu *et al.*[^Livesu_2021] reproduced here below

![](img/Livesu_Fig_2.png)

we examine several unit test singleton elements and their metrics.

valence | singleton | ${\rm ER}_{\max}$ | ${\rm SJ}_{\min}$ | ${\rm skew_{\max}}$ | volume
:---: | :---: | :---: | :---: | :---: | :---:
3           | ![](img/single_valence_03.png)        | 1.000000e0 (1.000)    | 8.660253e-1 (0.866)   | 5.000002e-1 (0.500)   | 8.660250e-1 (0.866)
3' (noised) | ![](img/single_valence_03_noise1.png) | 1.292260e0 (2.325) ** *Cubit (aspect ratio): 1.292* | 1.917367e-1 (0.192)   | 6.797483e-1 (0.680)   | 1.247800e0  (1.248)
4           | ![](img/single_valence_04.png)        | 1.000000e0 (1.000)    | 1.000000e0  (1.000)   | 0.000000e0  (0.000)   | 1.000000e0  (1.000)
4' (noised) | ![](img/single_valence_04_noise2.png) | 1.167884e0 (1.727) ** *Cubit (aspect ratio): 1.168* | 3.743932e-1 (0.374)   | 4.864936e-1 (0.486)   | 9.844008e-1 (0.984)
5           | ![](img/single_valence_05.png)        | 1.000000e0 (1.000)    | 9.510566e-1 (0.951)   | 3.090169e-1 (0.309)   | 9.510570e-1 (0.951)
6           | ![](img/single_valence_06.png)        | 1.000000e0 (1.000)    | 8.660253e-1 (0.866)   | 5.000002e-1 (0.500)   | 8.660250e-1 (0.866)
...         | ...                                   | ...                   | ...                   | ...                   | ...
10          | ![](img/single_valence_10.png)        | 1.000000e0 (1.000)    | 5.877851e-1 (0.588)   | 8.090171e-1 (0.809)   |  5.877850e-1 (0.588)

Figure: Maximum edge ratio, minimum scaled Jacobian, maximum skew, and volume.
Leading values are from `automesh`.
Values in parenthesis are results from [HexaLab](https://www.hexalab.net).[^Hexalab_2023]
Items with ** indicate where `automesh` and Cubit agree, but HexaLab disagrees.
  Cubit uses the term *Aspect Ratio* for Edge Ratio.

The connectivity for all elements:

```sh
1,    2,    4,    3,    5,    6,    8,    7

with prototype:

The element coordinates follow:

# 3
    1,      0.000000e0,      0.000000e0,      0.000000e0
    2,      1.000000e0,      0.000000e0,      0.000000e0
    3,     -0.500000e0,      0.866025e0,      0.000000e0
    4,      0.500000e0,      0.866025e0,      0.000000e0
    5,      0.000000e0,      0.000000e0,      1.000000e0
    6,      1.000000e0,      0.000000e0,      1.000000e0
    7,     -0.500000e0,      0.866025e0,      1.000000e0
    8,      0.500000e0,      0.866025e0,      1.000000e0

# 3'
    1,      0.110000e0,      0.120000e0,     -0.130000e0
    2,      1.200000e0,     -0.200000e0,      0.000000e0
    3,     -0.500000e0,      1.866025e0,     -0.200000e0
    4,      0.500000e0,      0.866025e0,     -0.400000e0
    5,      0.000000e0,      0.000000e0,      1.000000e0
    6,      1.000000e0,      0.000000e0,      1.000000e0
    7,     -0.500000e0,      0.600000e0,      1.400000e0
    8,      0.500000e0,      0.866025e0,      1.200000e0

# 4
    1,      0.000000e0,      0.000000e0,      0.000000e0
    2,      1.000000e0,      0.000000e0,      0.000000e0
    3,      0.000000e0,      1.000000e0,      0.000000e0
    4,      1.000000e0,      1.000000e0,      0.000000e0
    5,      0.000000e0,      0.000000e0,      1.000000e0
    6,      1.000000e0,      0.000000e0,      1.000000e0
    7,      0.000000e0,      1.000000e0,      1.000000e0
    8,      1.000000e0,      1.000000e0,      1.000000e0

# 4'
    1,      0.100000e0,      0.200000e0,      0.300000e0
    2,      1.200000e0,      0.300000e0,      0.400000e0
    3,     -0.200000e0,      1.200000e0,     -0.100000e0
    4,      1.030000e0,      1.102000e0,     -0.250000e0
    5,     -0.001000e0,     -0.021000e0,      1.002000e0
    6,      1.200000e0,     -0.100000e0,      1.100000e0
    7,      0.000000e0,      1.000000e0,      1.000000e0
    8,      1.010000e0,      1.020000e0,      1.030000e0

# 5
    1,      0.000000e0,      0.000000e0,      0.000000e0
    2,      1.000000e0,      0.000000e0,      0.000000e0
    3,      0.309017e0,      0.951057e0,      0.000000e0
    4,      1.309017e0,      0.951057e0,      0.000000e0
    5,      0.000000e0,      0.000000e0,      1.000000e0
    6,      1.000000e0,      0.000000e0,      1.000000e0
    7,      0.309017e0,      0.951057e0,      1.000000e0
    8,      1.309017e0,      0.951057e0,      1.000000e0

# 6
    1,      0.000000e0,      0.000000e0,      0.000000e0
    2,      1.000000e0,      0.000000e0,      0.000000e0
    3,      0.500000e0,      0.866025e0,      0.000000e0
    4,      1.500000e0,      0.866025e0,      0.000000e0
    5,      0.000000e0,      0.000000e0,      1.000000e0
    6,      1.000000e0,      0.000000e0,      1.000000e0
    7,      0.500000e0,      0.866025e0,      1.000000e0
    8,      1.500000e0,      0.866025e0,      1.000000e0

# 10
    1,      0.000000e0,      0.000000e0,      0.000000e0
    2,      1.000000e0,      0.000000e0,      0.000000e0
    3,      0.809017e0,      0.587785e0,      0.000000e0
    4,      1.809017e0,      0.587785e0,      0.000000e0
    5,      0.000000e0,      0.000000e0,      1.000000e0
    6,      1.000000e0,      0.000000e0,      1.000000e0
    7,      0.809017e0,      0.587785e0,      1.000000e0
    8,      1.809017e0,      0.587785e0,      1.000000e0

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 <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8333em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord"><span class="mord mathrm">ER</span></span></span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">a</span><span class="mtight">x</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span>
* Minium scaled Jacobian <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8333em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord"><span class="mord mathrm">SJ</span></span></span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span>
* Maximum skew
* Element area
* Minimum angle

A brief description of each metric follows.

### Maximum Edge Ratio

* <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8333em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord"><span class="mord mathrm">ER</span></span></span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">a</span><span class="mtight">x</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span> 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.*[^Knupp_2006] (page 26) indicate an acceptable range of `[1.0, 1.3]`.

### Minimum Scaled Jacobian

* <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8333em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord"><span class="mord mathrm">SJ</span></span></span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span> 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.*[^Knupp_2006] (page 29) indicate an acceptable range of `[0.5, 2*sqrt(3)/3]` <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.4831em;"></span><span class="mrel">≈</span></span></span></span> `[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:

<span class="katex-display"><span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:2.363em;vertical-align:-0.936em;"></span><span class="mord"><span class="mord"><span class="mord mathrm">S</span><span class="mord"><span class="mord mathrm">J</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">=</span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mord"><span class="mopen nulldelimiter"></span><span class="mfrac"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:1.427em;"><span style="top:-2.314em;"><span class="pstrut" style="height:3em;"></span><span class="mord"><span class="mop">sin</span><span class="mopen">(</span><span class="mord mathrm">6</span><span class="mord"><span class="mord mathrm">0</span><span class="msupsub"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.6001em;"><span style="top:-2.989em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mord mtight">∘</span></span></span></span></span></span></span></span></span><span class="mclose">)</span></span></span><span style="top:-3.23em;"><span class="pstrut" style="height:3em;"></span><span class="frac-line" style="border-bottom-width:0.04em;"></span></span><span style="top:-3.677em;"><span class="pstrut" style="height:3em;"></span><span class="mord"><span class="mop">sin</span><span class="mopen">(</span><span class="mord"><span class="mord mathnormal" style="margin-right:0.02778em;">θ</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:-0.0278em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mclose">)</span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.936em;"><span></span></span></span></span></span><span class="mclose nulldelimiter"></span></span></span></span></span></span></span></span>

* In the preceeding equation,
  * <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6741em;"></span><span class="mord">6</span><span class="mord"><span class="mord">0</span><span class="msupsub"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.6741em;"><span style="top:-3.063em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mord mtight">∘</span></span></span></span></span></span></span></span></span></span></span></span> is the smallest angle of an equilateral triangle, and
  * <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord mathnormal" style="margin-right:0.02778em;">θ</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:-0.0278em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span> 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.*[^Knupp_2006] does not give a definition of skew for triangles, so we provide our definition below.
For a triangle where <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord mathnormal" style="margin-right:0.02778em;">θ</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:-0.0278em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span> is the smallest angle of the triangle,

<span class="katex-display"><span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord mathrm">ske</span><span class="mord"><span class="mord mathrm" style="margin-right:0.01389em;">w</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:-0.0139em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">a</span><span class="mtight">x</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">=</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:2.0574em;vertical-align:-0.686em;"></span><span class="mord"><span class="mopen nulldelimiter"></span><span class="mfrac"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:1.3714em;"><span style="top:-2.314em;"><span class="pstrut" style="height:3em;"></span><span class="mord"><span class="mord">6</span><span class="mord"><span class="mord">0</span><span class="msupsub"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.6001em;"><span style="top:-2.989em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mord mtight">∘</span></span></span></span></span></span></span></span></span></span></span><span style="top:-3.23em;"><span class="pstrut" style="height:3em;"></span><span class="frac-line" style="border-bottom-width:0.04em;"></span></span><span style="top:-3.677em;"><span class="pstrut" style="height:3em;"></span><span class="mord"><span class="mord">6</span><span class="mord"><span class="mord">0</span><span class="msupsub"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.6741em;"><span style="top:-3.063em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mord mtight">∘</span></span></span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mord"><span class="mord mathnormal" style="margin-right:0.02778em;">θ</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:-0.0278em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.686em;"><span></span></span></span></span></span><span class="mclose nulldelimiter"></span></span></span></span></span></span>

* For an equilateral triangle, <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord mathnormal" style="margin-right:0.02778em;">θ</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:-0.0278em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">=</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:0.6741em;"></span><span class="mord">6</span><span class="mord"><span class="mord">0</span><span class="msupsub"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.6741em;"><span style="top:-3.063em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mord mtight">∘</span></span></span></span></span></span></span></span></span></span></span></span> and <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord mathrm">ske</span><span class="mord"><span class="mord mathrm" style="margin-right:0.01389em;">w</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:-0.0139em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">a</span><span class="mtight">x</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">=</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">0</span></span></span></span>.
* In the limit as <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord mathnormal" style="margin-right:0.02778em;">θ</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3175em;"><span style="top:-2.55em;margin-left:-0.0278em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">i</span><span class="mtight">n</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">→</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:0.6741em;"></span><span class="mord"><span class="mord">0</span><span class="msupsub"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.6741em;"><span style="top:-3.063em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mord mtight">∘</span></span></span></span></span></span></span></span></span></span></span></span> <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8444em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord"><span class="mord mathrm">ske</span><span class="mord"><span class="mord mathrm" style="margin-right:0.01389em;">w</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:-0.0139em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mtight"><span class="mop mtight"><span class="mtight">m</span><span class="mtight">a</span><span class="mtight">x</span></span></span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">→</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">1</span></span></span></span>.

### Element Area

* Measures the area of the element.

### Minimum Angle

* The smallest the three angles of a triangle.

## Triangular 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:

```sh
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:

fileearea (deg)
A11.508 [1.508]0.761 [0.761]0.313 [0.331]0.610 [0.610]41.2 [41.2]
A21.550 [1.550]0.739 [0.739]0.337 [0.337]0.550 [0.550]39.8 [39.8]
A31.787 [1.787]0.639 [0.639]0.440 [0.440]0.569 [0.569]33.6 [33.6]
A41.915 [1.915]0.595 [0.595]0.483 [0.483]0.402 [0.402]31.0 [31.0]
A52.230 [2.230]0.426 [0.426]0.639 [0.639]0.342 [0.342]21.7 [21.7]
A61.623 [1.623]0.700 [0.700]0.378 [0.378]0.571 [0.571]37.3 [37.1]
A71.240 [1.240]0.898 [0.898]0.149 [0.149]0.424 [0.424]51.0 [51.0]
A81.385 [1.385]0.831 [0.831]0.233 [0.233]0.443 [0.443]46.1 [46.1]
A91.606 [1.606]0.719 [0.719]0.358 [0.358]0.648 [0.648]38.5 [38.5]
A101.429 [1.429]0.806 [0.806]0.262 [0.262]0.704 [0.704]44.3 [44.3]
A111.275 [1.275]0.880 [0.880]0.172 [0.172]0.668 [0.668]49.7 [49.7]
A121.436 [1.436]0.804 [0.804]0.264 [0.264]0.516 [0.516]44.1 [44.1]
B11.414 [1.414]0.816 [0.816]0.250 [0.250]0.500 [0.500]45.0 [45.0]
C11.000 [1.000]1.000 [1.000]0.000 [0.000]6.928 [6.928]60.0 [60.0]
D11.000 [1.000]1.000 [1.000]0.000 [0.000]0.433 [0.433]60.0 [60.0]
E11.256 [1.256]0.869 [0.869]0.187 [0.187]3.273 [3.273]48.8 [48.8]

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.

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