Add quality metrics for quads (and all elements for some) by GiudGiud · Pull Request #4551 · libMesh/libmesh

GiudGiud · 2026-09-14T17:49:43Z

not rebased yet, too early for review

a different skew metric
condition number
size
taper for quads

aside from different skew metric, mostly filling unimplemented options

moosebuild · 2026-09-14T19:43:16Z

Job Coverage, step Generate coverage on a41acb4 wanted to post the following:

Coverage

	54e0d5	#4551 a41acb
	Total	Total	+/-	New
Rate	65.99%	66.30%	+0.31%	100.00%
Hits	79562	80010	+448	243
Misses	40996	40664	-332	0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

The existing SKEW metric implements Knupp's algebraic skew (1 is ideal), but its help string and CUBIT-derived suggested ranges describe the Verdict "skew": the maximum |cos A|, where A is the angle between the element's principal axes (0 is ideal). Rather than change SKEW's behavior, add a distinct SKEW_ANGLE metric implementing the Verdict definition, and correct SKEW's help string to describe what it actually computes. For a QUAD, the two principal axes are the midpoint-to-midpoint vectors of opposite edges; SKEW_ANGLE is |cos| of the angle between them. Co-Authored-By: Claude <noreply@anthropic.com>

Implement the Verdict "skew" metric for HEX elements, matching the SKEW_ANGLE metric added for quads. A hex has three principal axes, each the sum of the midpoint-to-midpoint vectors of opposite faces along one logical direction; SKEW_ANGLE is the maximum |cos| over the three pairs of axes (0 is ideal, 1 is fully skewed). Co-Authored-By: Claude <noreply@anthropic.com>

Test QUAD4 and HEX8 SKEW_ANGLE against closed-form values: a unit square/cube is 0 (orthogonal); a rhombus with interior angle theta is |cos(theta)| (also checked to be rotation invariant); a unit cube sheared by k in x is k/sqrt(k^2+1); and degenerate elements return 0. Co-Authored-By: Claude <noreply@anthropic.com>

SIZE was advertised as a valid metric with a description and suggested ranges, but no element implemented it, so quality(SIZE) fell through to the Elem::quality default and silently returned 1. Implement it generically in Elem::quality (covering quads, hexes, tris, tets) as the relative size min(J, 1/J), where J is the determinant of the nodal Jacobian. Following the other algebraic metrics (SHAPE, SKEW, JACOBIAN), the reference/weight matrix is the identity, i.e. the unit reference element, so J is the element's nodal Jacobian determinant (area in 2D, volume in 3D) averaged over the corner nodes. Both over- and undersized elements are penalized and a unit-sized element scores 1. This differs from the Verdict/CUBIT relative size, which normalizes by the mesh-average element size (not available to a per-element method), and it does not square J. Co-Authored-By: Claude <noreply@anthropic.com>

Test QUAD4 and HEX8 SIZE against closed-form values: a unit square/cube is 1; a 2x2 square is 0.25 and a 2x2x2 cube is 0.125; 2:1 rectangles and boxes are 0.5; and a unit-edge rhombus with interior angle pi/6 is sin(pi/6) = 0.5. Co-Authored-By: Claude <noreply@anthropic.com>

TAPER was advertised as a valid Quad metric with a description and suggested ranges, but Quad::quality had no TAPER case, so it fell through to the Elem::quality default and silently returned 1. Implement it as the maximum ratio of lengths derived from opposite edges, using the same convention as the Hex TAPER metric (the Quad is its single-face case): for each of the two pairs of opposite edges, form the ratio of the shorter to the longer length, and return the smallest (worst) such ratio. The value lies in (0, 1], with 1 meaning no taper (both pairs of opposite edges equal, as for any parallelogram). Co-Authored-By: Claude <noreply@anthropic.com>

Two problems found while reviewing the SKEW_ANGLE/SIZE/TAPER commits: 1. Hex qual_bounds: SKEW was left grouped with SKEW_ANGLE at (0 -> 0.5), but the corrected SKEW description (and the Quad bounds) put Knupp's algebraic skew at (0.3 -> 1). Give SKEW the (0.3 -> 1) range with the other Knupp metrics and leave SKEW_ANGLE at the Verdict (0 -> 0.5). 2. SIZE: "std::min(J, 1. / J)" mixes Real and double, which fails to compile under a single-precision (Real == float) configuration. Use Real(1) / J so both arguments are Real. Co-Authored-By: Claude <noreply@anthropic.com>

CONDITION was advertised as a valid metric with a description and suggested ranges, but no element implemented it, so quality(CONDITION) fell through to the Elem::quality default and silently returned 1. Implement it generically in Elem::quality (covering quads, hexes, tris, tets) as the maximum over the corner nodes of the nodal Jacobian condition number kappa = |A|_F * |A^{-1}|_F / N. Following the other algebraic metrics (SHAPE, SKEW), the reference matrix is the identity, so kappa = 1 for an orthogonal, equal-length (ideal) corner and grows with stretch or skew. For a 2x2 Jacobian |A^{-1}|_F = |A|_F / |det|; for 3x3 the inverse rows are the pairwise edge cross products over det. A degenerate corner (zero determinant) has an infinite condition number, reported as 0 following the convention that 0 stands in for infinity. Co-Authored-By: Claude <noreply@anthropic.com>

Test QUAD4 and HEX8 CONDITION against closed-form values: a unit square/cube is 1; a 2x1 rectangle is 1.25; a unit-edge rhombus with interior angle theta is 1/sin(theta); and a 2x1x1 box is sqrt(6)/2. Co-Authored-By: Claude <noreply@anthropic.com>

GiudGiud · 2026-09-16T15:47:42Z

                  MIN_DIHEDRAL_ANGLE,
-                  SCALED_JACOBIAN};
+                  SCALED_JACOBIAN,
+                  SKEW_ANGLE};


this because our current SKEW actually is "SHEAR" there:
https://coreform.com/cubit_help/mesh_generation/mesh_quality_assessment/quadrilateral_metrics.htm

GiudGiud · 2026-09-16T15:48:30Z

-    case SKEW:
+    case SKEW_ANGLE:
      bounds.first  = 0.;
      bounds.second = 0.5;


the ranges for our skew were wrong, 1 was best, 0 was bad

GiudGiud · 2026-09-16T15:49:41Z

+      // element scores the ideal value of 1. This differs from the
+      // Verdict/CUBIT relative size, which normalizes J by the
+      // mesh-average element size; that requires mesh-wide context not
+      // available to this per-element method, so we use the reference
+      // element instead. Unlike the standard Verdict metric, J is not
+      // squared here.


this makes it less useful. maybe to the point i should drop the change?

Oh, I forgot to respond to this bit, sorry.

IMHO normalizing by the reference element size isn't very useful, but normalizing by the "ideal" Jacobian, the |J| for a perfectly-shaped element of the same volume(), would give something valuable.

Oh, no. Oh no, no, no.

A "perfectly-shaped element" is not, in general, proportional to our master element. Our edge/square/cube elements are fine as exemplars, but our prisms+pyramids are debateable and our tris+tets are outright wrong. For this and for CONDITION we need to be getting the Jacobian for the transformation from an ideal element, not from a master element, just like we had to do in the variational mesh smoother.

ok using the ideal element from the variable smoother now

roystgnr

I'd love another unit test too, but beggars can't be choosers, and the describe() clarification is the important thing IMHO.

roystgnr · 2026-09-16T16:44:33Z

+        // Degenerate element: return 0 (the Verdict convention) if any
+        // principal axis has zero length.
+        if (n1 == 0. || n2 == 0. || n3 == 0.)
+          return 0.;


Not a fan of this. I guess I'm okay matching their convention, but let's make sure to put it in the describe() string too.

roystgnr · 2026-09-16T16:57:11Z

+      // determinant) has an infinite condition number, reported as 0
+      // following the convention used elsewhere (e.g. EDGE_LENGTH_RATIO)
+      // that 0 stands in for infinity.
+    case CONDITION:


We thoughtfully added qual_bounds(CONDITION) recommendations for 4 different categories of elements, but we didn't previously have an implementation for quality(CONDITION)?

Oh, wow. Looks like those bounds were added by Ben, in 2002 or 2003, in libMesh's 3rd non-trivial CVS commit, the 1st commit with files for some Elem subclasses, presumably as placeholders in the hope that somebody would add the implementation someday.

This PR is awesome. "Lafayette, we are here!"

roystgnr · 2026-09-16T17:03:12Z

+      desc << "Knupp's algebraic skew metric,\n"
+           << "based on the nodal Jacobian\n"
+           << "skew matrices. 1 is ideal,\n"
+           << "smaller values are worse.\n"
+           << '\n'
+           << "Suggested ranges:\n"
+           << "Hexes: (0.3 -> 1)\n"
+           << "Quads: (0.3 -> 1)";
+      break;
+
+    case SKEW_ANGLE:


Really good fix here, but pretty damning of our test coverage that we had bounds for one metric but were computing another.

In elem_test.C, should we add a new test that loops q over Quality::valid(elem->type()), evaluates elem->quality(a) and asserts (since we're on simple default meshes there) that the result is within elem->qual_bounds(q)?

added that test

The TAPER metric returns 1 for an untapered (ideal) element and decreases toward 0 as taper increases (see Quad::quality / Hex::quality), and the metric description gives the good ranges as Quads (0.7 -> 1) and Hexes (0.4 -> 1). But qual_bounds() had these reversed as (0, 0.7) and (0, 0.4), so an ideal element's TAPER value of 1 fell outside its own suggested bounds. Store the ranges the right way round. Co-Authored-By: Claude <noreply@anthropic.com>

For a well-shaped, unit-scale ideal element (equilateral triangle, unit square, regular tetrahedron, unit cube), loop over Quality::valid() for the type and assert each metric evaluates within its qual_bounds(). Metrics that are listed as valid but have no bounds defined (qual_bounds returns the (-1,-1) sentinel) are skipped. Other element types, for which a clean ideal shape is not readily constructed here, are skipped. This guards against metrics or bounds drifting out of sync, which is how the reversed TAPER bounds were found. Co-Authored-By: Claude <noreply@anthropic.com>

CONDITION and SIZE previously used the identity/master corner as the reference: a unit right-angle corner. That makes an equilateral triangle or regular tetrahedron score CONDITION 1.15 / 1.22 instead of 1, and ties SIZE to the master element's scale. Instead, measure both against the ideal (regular) element of the same volume, as the variational mesh smoother does with its target element. CONDITION: at each corner use the weighted Jacobian A * W^{-1}, where W is the ideal corner (unit edges at 60 degrees for a simplex, 90 otherwise). Via the corner metric tensors T_A, T_W (which also handles a lower-dimensional element embedded in 3D), kappa = sqrt(tr(T_A T_W^{-1}) * tr(T_W T_A^{-1})) / N. Equilateral triangles and regular tets now score 1; quads/hexes are unchanged (their ideal corner already is the right angle). SIZE: tau is the corner nodal Jacobian determinant divided by that of an ideal element of the same volume (the mean of the corner determinants), and SIZE = min over corners of min(tau, 1/tau). It is 1 for any affine element (uniform Jacobian) at any scale, and drops below 1 for non-affine (tapered/sheared) elements. This is nodal-Jacobian based, not elem->volume() based. Update the descriptions and the QUAD4/HEX8 SIZE tests (affine shapes now score 1; a trapezoid scores 0.8 and a frustum hex 0.4). Co-Authored-By: Claude <noreply@anthropic.com>

The variational smoother's get_target_elem builds the ideal (regular) element for a type -- equilateral triangle, regular tet, etc. -- sized to the reference element's volume. It is pure geometry (Elem::build, Node::build, reference_elem()), so move it into the ReferenceElem namespace in src/geom as ideal_target(), where the quality metrics can reuse it too. The smoother now calls ReferenceElem::ideal_target(); its get_target_to_reference_jacobian (which needs FEMContext) stays put. Pure refactor; no behavior change. Co-Authored-By: Claude <noreply@anthropic.com>

Replace the hand-rolled analytic ideal-corner tensors with the shared ReferenceElem::ideal_target element, so the quality metrics and the variational smoother derive "the ideal element" from one place. CONDITION: build the ideal element and take the ideal corner metric tensor T_W from its corresponding corner (scale-invariant, so the reference-volume sizing is fine). Equilateral triangles / regular tets still score 1; quads/hexes unchanged. SIZE: tau_k = (nodal_det_phys / nodal_det_ideal) * ideal_vol / this_vol, i.e. the ideal element rescaled to this element's volume, matching "same volume as the element". Affine elements score 1 at any scale; the non-affine hex frustum test now uses the true same-volume ideal value 3/7 (was 0.4 under the mean-of-nodal-dets proxy). Quad values are unchanged (mean nodal area equals the area there). Co-Authored-By: Claude <noreply@anthropic.com>

GiudGiud · 2026-09-17T14:38:35Z

+ * Elem.
+ */
+std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
+ideal_target (const ElemType type);


I put it here since reference not too far from ideal ideologically

could also just go in elem.h

GiudGiud · 2026-09-17T15:06:40Z

tests need fixing, on it

…m map These two ElemQuality values were missing from elemquality_to_enum, so Utility::enum_to_string() (and string_to_enum()) threw "No ElemQuality with enumeration N found" for them. Add both so every ElemQuality value round-trips to/from its string name. Co-Authored-By: Claude <noreply@anthropic.com>

GiudGiud and others added 9 commits September 16, 2026 09:46

GiudGiud force-pushed the quad-hex-skew-angle-metric branch from 5012aa8 to 36f69e8 Compare September 16, 2026 15:46

GiudGiud commented Sep 16, 2026

View reviewed changes

GiudGiud marked this pull request as ready for review September 16, 2026 15:50

roystgnr requested changes Sep 16, 2026

View reviewed changes

GiudGiud and others added 5 commits September 16, 2026 15:50

GiudGiud commented Sep 17, 2026

View reviewed changes

GiudGiud requested a review from roystgnr September 17, 2026 14:51

Add warning to metric description

414c695

GiudGiud force-pushed the quad-hex-skew-angle-metric branch from 2a03e73 to 414c695 Compare September 17, 2026 14:52

GiudGiud force-pushed the quad-hex-skew-angle-metric branch from 2832bc4 to a41acb4 Compare September 17, 2026 22:59

Conversation

GiudGiud commented Sep 14, 2026

Uh oh!

moosebuild commented Sep 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Coverage

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

roystgnr left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

roystgnr Sep 16, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

GiudGiud commented Sep 17, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

moosebuild commented Sep 14, 2026 •

edited

Loading

roystgnr Sep 16, 2026 •

edited

Loading