From fd0e56554c4d2f8ba00a87faa3b23c250b4e2c13 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 08:09:11 -0600
Subject: [PATCH 01/16] Add Verdict-style SKEW_ANGLE quality metric for quads

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>
---
 include/enums/enum_elem_quality.h |  3 ++-
 src/geom/elem_quality.C           | 16 +++++++++++++++
 src/geom/face_quad.C              | 33 +++++++++++++++++++++++++++++++
 src/utils/string_to_enum.C        |  1 +
 4 files changed, 52 insertions(+), 1 deletion(-)

diff --git a/include/enums/enum_elem_quality.h b/include/enums/enum_elem_quality.h
index abccd998785..1b7c30412a7 100644
--- a/include/enums/enum_elem_quality.h
+++ b/include/enums/enum_elem_quality.h
@@ -52,7 +52,8 @@ enum ElemQuality : int {
                   EDGE_LENGTH_RATIO,
                   MAX_DIHEDRAL_ANGLE,
                   MIN_DIHEDRAL_ANGLE,
-                  SCALED_JACOBIAN};
+                  SCALED_JACOBIAN,
+                  SKEW_ANGLE};
 }
 
 #endif
diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index 2353fb72d6e..b4fe55e5817 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -58,6 +58,10 @@ std::string Quality::name (const ElemQuality q)
       its_name = "Skew";
       break;
 
+    case SKEW_ANGLE:
+      its_name = "Skew Angle";
+      break;
+
     case SHEAR:
       its_name = "Shear";
       break;
@@ -162,6 +166,17 @@ std::string Quality::describe (const ElemQuality q)
       break;
 
     case SKEW:
+      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:
       desc << "Maximum |cos A|, where A\n"
            << "is the angle between edges\n"
            << "at element center.\n"
@@ -410,6 +425,7 @@ std::vector<ElemQuality> Quality::valid(const ElemType t)
           SHEAR,
           SIZE,
           SKEW,
+          SKEW_ANGLE,
           STRETCH,
           TAPER,
           WARP
diff --git a/src/geom/face_quad.C b/src/geom/face_quad.C
index c4ca8a57a11..cd721548e46 100644
--- a/src/geom/face_quad.C
+++ b/src/geom/face_quad.C
@@ -400,6 +400,34 @@ Real Quad::quality (const ElemQuality q) const
           }
       }
 
+      // Verdict/CUBIT "skew" metric: the maximum |cos A|, where A is
+      // the angle between the principal axes of the element. The
+      // principal axes are the vectors connecting the midpoints of
+      // opposite edges. A value of 0 indicates a perfectly orthogonal
+      // (unskewed) element; larger values (up to 1) indicate
+      // increasing skew. This differs from the SKEW metric above,
+      // which is Knupp's algebraic skew (1 is ideal).
+      // See: C. J. Stimpson et al., "The Verdict Geometric Quality
+      // Library," Sandia report SAND2007-1751, 2007.
+    case SKEW_ANGLE:
+      {
+        const Point x0 = this->point(0), x1 = this->point(1),
+                    x2 = this->point(2), x3 = this->point(3);
+
+        // Principal axes: midpoint-to-midpoint of opposite edges.
+        const Point X1 = (x1 - x0) + (x2 - x3);
+        const Point X2 = (x3 - x0) + (x2 - x1);
+
+        const Real n1 = X1.norm(), n2 = X2.norm();
+
+        // Degenerate element: return 0 (the Verdict convention) if
+        // either principal axis has zero length.
+        if (n1 == 0. || n2 == 0.)
+          return 0.;
+
+        return std::abs((X1 * X2) / (n1 * n2));
+      }
+
       // This test returns 0 if a Quad:
       // 1) Is "twisted" (i.e. has an invalid numbering)
       // 2) Has nearly parallel adjacent edges
@@ -541,6 +569,11 @@ std::pair<Real, Real> Quad::qual_bounds (const ElemQuality q) const
       bounds.second = 1.;
       break;
 
+    case SKEW_ANGLE:
+      bounds.first  = 0.;
+      bounds.second = 0.5;
+      break;
+
     case DISTORTION:
       bounds.first  = 0.6;
       bounds.second = 1.;
diff --git a/src/utils/string_to_enum.C b/src/utils/string_to_enum.C
index 18e86cd1c6e..7f911244a07 100644
--- a/src/utils/string_to_enum.C
+++ b/src/utils/string_to_enum.C
@@ -410,6 +410,7 @@ std::map<SolverType, std::string> enum_to_solvertype =
 std::map<std::string, ElemQuality> elemquality_to_enum {
    {"ASPECT_RATIO"       , ASPECT_RATIO},
    {"SKEW"               , SKEW},
+   {"SKEW_ANGLE"         , SKEW_ANGLE},
    {"SHEAR"              , SHEAR},
    {"SHAPE"              , SHAPE},
    {"MAX_ANGLE"          , MAX_ANGLE},

From 49b953ccd775052ac319c130d712294f8fa37a09 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 08:09:21 -0600
Subject: [PATCH 02/16] Add SKEW_ANGLE quality metric for hexes

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>
---
 src/geom/cell_hex.C     | 38 ++++++++++++++++++++++++++++++++++++++
 src/geom/elem_quality.C |  1 +
 2 files changed, 39 insertions(+)

diff --git a/src/geom/cell_hex.C b/src/geom/cell_hex.C
index e1c2bab0198..2d8f0c85330 100644
--- a/src/geom/cell_hex.C
+++ b/src/geom/cell_hex.C
@@ -491,6 +491,43 @@ Real Hex::quality (const ElemQuality q) const
             return (den == 0.) ? 0 : (8. / den);
           }
       }
+
+      // Verdict/CUBIT "skew" metric: the maximum |cos A| over the
+      // three pairs of principal axes, where A is the angle between a
+      // pair of axes. Each principal axis is the sum of the vectors
+      // connecting the midpoints of opposite faces along one logical
+      // direction. A value of 0 indicates a perfectly orthogonal
+      // (unskewed) element; larger values (up to 1) indicate
+      // increasing skew. This differs from the SKEW metric above,
+      // which is Knupp's algebraic skew (1 is ideal).
+      // See: C. J. Stimpson et al., "The Verdict Geometric Quality
+      // Library," Sandia report SAND2007-1751, 2007.
+    case SKEW_ANGLE:
+      {
+        const Point
+          x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3),
+          x4 = point(4), x5 = point(5), x6 = point(6), x7 = point(7);
+
+        // Principal axes, one per logical (xi, eta, zeta) direction.
+        const Point
+          X1 = (x1 - x0) + (x2 - x3) + (x5 - x4) + (x6 - x7),
+          X2 = (x3 - x0) + (x2 - x1) + (x7 - x4) + (x6 - x5),
+          X3 = (x4 - x0) + (x5 - x1) + (x6 - x2) + (x7 - x3);
+
+        const Real n1 = X1.norm(), n2 = X2.norm(), n3 = X3.norm();
+
+        // Degenerate element: return 0 (the Verdict convention) if any
+        // principal axis has zero length.
+        if (n1 == 0. || n2 == 0. || n3 == 0.)
+          return 0.;
+
+        // Normalize, then take the largest |cos| among the three
+        // pairs of principal axes.
+        const Point X1h = X1 / n1, X2h = X2 / n2, X3h = X3 / n3;
+        return std::max({std::abs(X1h * X2h),
+                         std::abs(X1h * X3h),
+                         std::abs(X2h * X3h)});
+      }
 #endif // LIBMESH_DIM >= 3
 
       /**
@@ -517,6 +554,7 @@ std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
       break;
 
     case SKEW:
+    case SKEW_ANGLE:
       bounds.first  = 0.;
       bounds.second = 0.5;
       break;
diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index b4fe55e5817..a32272a8f16 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -475,6 +475,7 @@ std::vector<ElemQuality> Quality::valid(const ElemType t)
           SHEAR,
           SIZE,
           SKEW,
+          SKEW_ANGLE,
           STRETCH,
           TAPER
         };

From 39d652bbe27abcbdcaba8f1d13f612ce88e3f305 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 08:09:26 -0600
Subject: [PATCH 03/16] Add unit tests for the SKEW_ANGLE quality metric

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>
---
 tests/geom/volume_test.C | 110 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 110 insertions(+)

diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index fa7e9ce995d..18ed9d140c3 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -42,6 +42,8 @@ public:
   CPPUNIT_TEST( testQuad4Warpage );
   CPPUNIT_TEST( testQuad4MinMaxAngle );
   CPPUNIT_TEST( testQuad4Jacobian );
+  CPPUNIT_TEST( testQuad4SkewAngle );
+  CPPUNIT_TEST( testHex8SkewAngle );
   CPPUNIT_TEST( testTri3AspectRatio );
   CPPUNIT_TEST( testTet4DihedralAngle );
   CPPUNIT_TEST( testTet4Jacobian );
@@ -823,6 +825,114 @@ public:
     }
   }
 
+  void testQuad4SkewAngle()
+  {
+    LOG_UNIT_TEST;
+
+    // The SKEW_ANGLE metric is the Verdict "skew": the maximum |cos A|
+    // between the element's principal axes. 0 means perfectly
+    // orthogonal (unskewed), larger (up to 1) means more skewed.
+
+    // Case 1: A unit square has orthogonal principal axes, so its skew
+    // angle is exactly 0.
+    {
+      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0)};
+      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
+      libmesh_ignore(nodes);
+
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.0, /*actual=*/elem->quality(SKEW_ANGLE), TOLERANCE);
+    }
+
+    // Case 2: For a rhombus with interior angle theta, the two
+    // principal axes are separated by theta, so the skew metric is
+    // |cos(theta)|. This is also invariant to rigid body rotation, so
+    // we rotate the rhombus about the z-axis before checking.
+    {
+      auto test_rhombus_quad = [this](Real theta)
+      {
+        const Real ct = std::cos(theta);
+        const Real st = std::sin(theta);
+        std::vector<Point> pts = {
+          Point(0, 0, 0),
+          Point(1, 0, 0),
+          Point(1. + ct, st, 0),
+          Point(     ct, st, 0)};
+
+        // Rotate all points about the z-axis by 30 degrees to confirm
+        // the metric is rotation invariant.
+        const Real cr = std::cos(libMesh::pi / 6);
+        const Real sr = std::sin(libMesh::pi / 6);
+        RealTensorValue Rz(cr, -sr, 0,
+                           sr,  cr, 0,
+                           0,    0, 1);
+        for (auto & pt : pts)
+          pt = Rz * pt;
+
+        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
+        libmesh_ignore(nodes);
+
+        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::abs(ct), /*actual=*/elem->quality(SKEW_ANGLE), TOLERANCE);
+      };
+
+      // theta = pi/2 -> |cos| = 0 (orthogonal)
+      test_rhombus_quad(libMesh::pi / 2);
+      // theta = pi/3 -> |cos| = 0.5
+      test_rhombus_quad(libMesh::pi / 3);
+      // theta = pi/6 -> |cos| = sqrt(3)/2
+      test_rhombus_quad(libMesh::pi / 6);
+    }
+
+    // Case 3: A degenerate quad with a zero-length principal axis
+    // returns 0 by the Verdict convention.
+    {
+      std::vector<Point> pts = {Point(0, 0, 0), Point(0, 0, 0), Point(1, 1, 0), Point(1, 1, 0)};
+      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
+      libmesh_ignore(nodes);
+
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.0, /*actual=*/elem->quality(SKEW_ANGLE), TOLERANCE);
+    }
+  }
+
+  void testHex8SkewAngle()
+  {
+    LOG_UNIT_TEST;
+
+    // Case 1: A unit cube has mutually orthogonal principal axes, so
+    // its skew angle is exactly 0.
+    {
+      std::vector<Point> pts = {
+        Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
+        Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)};
+      auto [elem, nodes] = this->construct_elem(pts, HEX8);
+      libmesh_ignore(nodes);
+
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.0, /*actual=*/elem->quality(SKEW_ANGLE), TOLERANCE);
+    }
+
+    // Case 2: Shear the top face of a unit cube by k in the
+    // x-direction. Only the zeta principal axis tilts, becoming
+    // (k, 0, 1); the xi and eta axes stay orthogonal. The largest
+    // |cos| between any pair of axes is then k / sqrt(k^2 + 1).
+    {
+      auto test_sheared_hex = [this](Real k)
+      {
+        std::vector<Point> pts = {
+          Point(0,   0, 0), Point(1,   0, 0), Point(1,   1, 0), Point(0,   1, 0),
+          Point(k,   0, 1), Point(1+k, 0, 1), Point(1+k, 1, 1), Point(k,   1, 1)};
+        auto [elem, nodes] = this->construct_elem(pts, HEX8);
+        libmesh_ignore(nodes);
+
+        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::abs(k) / std::sqrt(k*k + 1),
+                                /*actual=*/elem->quality(SKEW_ANGLE), TOLERANCE);
+      };
+
+      // k = 1 -> 1/sqrt(2) ~ 0.7071
+      test_sheared_hex(1.0);
+      // k = 0.5 -> 0.5/sqrt(1.25) ~ 0.4472
+      test_sheared_hex(0.5);
+    }
+  }
+
   void testTet4DihedralAngle()
   {
     LOG_UNIT_TEST;

From d48a69f06ea17830adf797b709b6273f27e65726 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 09:57:54 -0600
Subject: [PATCH 04/16] Implement the SIZE (relative size) quality metric

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>
---
 src/geom/elem.C         | 73 +++++++++++++++++++++++++++++++++++++++++
 src/geom/elem_quality.C |  7 ++--
 2 files changed, 77 insertions(+), 3 deletions(-)

diff --git a/src/geom/elem.C b/src/geom/elem.C
index 9440ddd9307..1f61c0dd422 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -1996,6 +1996,79 @@ Real Elem::quality (const ElemQuality q) const
         return min_node_area;
       }
 
+      // Relative size metric: min(J, 1/J), where J is the determinant
+      // of the "weighted" nodal Jacobian, A * W^{-1}. Following the
+      // other algebraic metrics (SHAPE, SKEW, JACOBIAN), the reference
+      // (weight) matrix W is the identity, i.e. the canonical unit
+      // reference element (unit-length edges meeting at right angles),
+      // for which det(W) = 1. J is therefore the element's nodal
+      // Jacobian determinant (area in 2D, volume in 3D spanned by the
+      // edges meeting at a node), averaged over the corner nodes.
+      //
+      // Both undersized (J < 1) and oversized (J > 1) elements are
+      // penalized, and an element the size of the unit reference
+      // 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.
+    case SIZE:
+      {
+        // 1D elements don't have interior corners, so this metric does
+        // not really apply to them.
+        const auto N = this->dim();
+        if (N < 2)
+          return 1.;
+
+        // Average the nodal Jacobian determinant over the corner
+        // nodes. This uses the same nodal Jacobian construction as the
+        // JACOBIAN metric above.
+        Real sum_node_area = 0.;
+        unsigned int n_corners = 0;
+
+        for (auto n : this->node_index_range())
+          {
+            // Get list of edge ids adjacent to this node.
+            auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
+
+            // Skip any nodes that don't have dim() adjacent edges (see
+            // the JACOBIAN metric above for the Pyramid apex caveat).
+            if (adjacent_edge_ids.size() != N)
+              continue;
+
+            // Construct oriented edges pointing away from node n.
+            std::vector<Point> oriented_edges(N);
+            for (auto i : make_range(N))
+              {
+                auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
+                auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
+                if (node_0 != n)
+                  std::swap(node_0, node_1);
+                oriented_edges[i] = this->point(node_1) - this->point(node_0);
+              }
+
+            // Unscaled nodal area (2D) or volume (3D).
+            Real node_area = (N == 2) ?
+              cross_norm(oriented_edges[0], oriented_edges[1]) :
+              std::abs(triple_product(oriented_edges[0], oriented_edges[1], oriented_edges[2]));
+
+            sum_node_area += node_area;
+            ++n_corners;
+          }
+
+        // No usable corners, or a degenerate (zero-size) element: return
+        // 0 (the lowest quality).
+        if (n_corners == 0)
+          return 0.;
+
+        const Real J = sum_node_area / n_corners;
+        if (J == 0.)
+          return 0.;
+
+        return std::min(J, 1. / J);
+      }
+
       // Return 1 if we made it here
     default:
       {
diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index a32272a8f16..f4f06454286 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -335,9 +335,10 @@ std::string Quality::describe (const ElemQuality q)
       break;
 
     case SIZE:
-      desc << "min (|J|, |1/J|)\n"
-           << '\n'
-           << "|J| = norm of Jacobian matrix.\n"
+      desc << "Relative size: min(J, 1/J),\n"
+           << "where J is the determinant\n"
+           << "of the nodal Jacobian relative\n"
+           << "to the unit reference element.\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (0.3 -> 1)\n"

From b9f9e070aa9cdfa0fb261d9cd17e95f072eeffca Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 09:57:59 -0600
Subject: [PATCH 05/16] Add unit tests for the SIZE quality metric

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>
---
 tests/geom/volume_test.C | 65 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 65 insertions(+)

diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index 18ed9d140c3..ea32f420cc5 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -44,6 +44,8 @@ public:
   CPPUNIT_TEST( testQuad4Jacobian );
   CPPUNIT_TEST( testQuad4SkewAngle );
   CPPUNIT_TEST( testHex8SkewAngle );
+  CPPUNIT_TEST( testQuad4Size );
+  CPPUNIT_TEST( testHex8Size );
   CPPUNIT_TEST( testTri3AspectRatio );
   CPPUNIT_TEST( testTet4DihedralAngle );
   CPPUNIT_TEST( testTet4Jacobian );
@@ -933,6 +935,69 @@ public:
     }
   }
 
+  void testQuad4Size()
+  {
+    LOG_UNIT_TEST;
+
+    // Relative SIZE = min(J, 1/J), where J is the element's nodal
+    // Jacobian determinant measured against the unit reference
+    // element. A unit square is ideal (1); larger and smaller elements
+    // score below 1 symmetrically.
+    auto size_of = [this](const std::vector<Point> & pts)
+    {
+      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
+      libmesh_ignore(nodes);
+      return elem->quality(SIZE);
+    };
+
+    // Unit square -> J = 1 -> 1
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+      /*actual=*/size_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0)}), TOLERANCE);
+
+    // 2x2 square -> nodal area J = 4 -> min(4, 1/4) = 0.25
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.25,
+      /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,2,0), Point(0,2,0)}), TOLERANCE);
+
+    // 2x1 rectangle -> J = 2 -> 0.5
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.5,
+      /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0)}), TOLERANCE);
+
+    // Unit-edge rhombus with interior angle pi/6 -> J = sin(pi/6) = 0.5
+    // -> min(0.5, 2) = 0.5
+    {
+      const Real c = std::cos(libMesh::pi/6), s = std::sin(libMesh::pi/6);
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::sin(libMesh::pi/6),
+        /*actual=*/size_of({Point(0,0,0), Point(1,0,0), Point(1.+c,s,0), Point(c,s,0)}), TOLERANCE);
+    }
+  }
+
+  void testHex8Size()
+  {
+    LOG_UNIT_TEST;
+
+    auto size_of = [this](const std::vector<Point> & pts)
+    {
+      auto [elem, nodes] = this->construct_elem(pts, HEX8);
+      libmesh_ignore(nodes);
+      return elem->quality(SIZE);
+    };
+
+    // Unit cube -> J = 1 -> 1
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+      /*actual=*/size_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0),
+                          Point(0,0,1), Point(1,0,1), Point(1,1,1), Point(0,1,1)}), TOLERANCE);
+
+    // 2x2x2 cube -> nodal volume J = 8 -> min(8, 1/8) = 0.125
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.125,
+      /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,2,0), Point(0,2,0),
+                          Point(0,0,2), Point(2,0,2), Point(2,2,2), Point(0,2,2)}), TOLERANCE);
+
+    // 2x1x1 box -> J = 2 -> 0.5
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.5,
+      /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0),
+                          Point(0,0,1), Point(2,0,1), Point(2,1,1), Point(0,1,1)}), TOLERANCE);
+  }
+
   void testTet4DihedralAngle()
   {
     LOG_UNIT_TEST;

From 98b747d42c71f0e47990e6c3e95c6a4d21a3629c Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 10:04:17 -0600
Subject: [PATCH 06/16] Implement the TAPER quality metric for quads

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>
---
 src/geom/face_quad.C     | 27 +++++++++++++++++++++++++++
 tests/geom/volume_test.C | 35 +++++++++++++++++++++++++++++++++++
 2 files changed, 62 insertions(+)

diff --git a/src/geom/face_quad.C b/src/geom/face_quad.C
index cd721548e46..bbd084e7581 100644
--- a/src/geom/face_quad.C
+++ b/src/geom/face_quad.C
@@ -324,6 +324,33 @@ Real Quad::quality (const ElemQuality q) const
           return std::sqrt(2) * min_edge / d_max;
       }
 
+      // Maximum ratio of lengths derived from opposite edges. This uses
+      // the same convention as the Hex TAPER metric (of which the Quad
+      // is the single-face case): for each of the two pairs of opposite
+      // edges we 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 indicating no taper, i.e. both pairs of opposite
+      // edges are equal in length (as for any parallelogram).
+    case TAPER:
+      {
+        const Real d01 = this->length(0,1);
+        const Real d12 = this->length(1,2);
+        const Real d23 = this->length(2,3);
+        const Real d03 = this->length(0,3);
+
+        // Longer length of each opposite-edge pair.
+        const Real max0 = std::max(d01, d23);
+        const Real max1 = std::max(d12, d03);
+
+        // Degenerate element with a zero-length pair of opposite edges:
+        // return 0 (the lowest quality).
+        if (max0 == 0. || max1 == 0.)
+          return 0.;
+
+        return std::min(std::min(d01, d23) / max0,
+                        std::min(d12, d03) / max1);
+      }
+
     case SHAPE:
     case SKEW:
       {
diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index ea32f420cc5..588dc7e9bb7 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -46,6 +46,7 @@ public:
   CPPUNIT_TEST( testHex8SkewAngle );
   CPPUNIT_TEST( testQuad4Size );
   CPPUNIT_TEST( testHex8Size );
+  CPPUNIT_TEST( testQuad4Taper );
   CPPUNIT_TEST( testTri3AspectRatio );
   CPPUNIT_TEST( testTet4DihedralAngle );
   CPPUNIT_TEST( testTet4Jacobian );
@@ -998,6 +999,40 @@ public:
                           Point(0,0,1), Point(2,0,1), Point(2,1,1), Point(0,1,1)}), TOLERANCE);
   }
 
+  void testQuad4Taper()
+  {
+    LOG_UNIT_TEST;
+
+    // TAPER = min over the two opposite-edge pairs of (shorter/longer
+    // length), in (0, 1], with 1 meaning no taper.
+    auto taper_of = [this](const std::vector<Point> & pts)
+    {
+      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
+      libmesh_ignore(nodes);
+      return elem->quality(TAPER);
+    };
+
+    // Unit square -> both pairs equal -> 1
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+      /*actual=*/taper_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0)}), TOLERANCE);
+
+    // 2x1 rectangle -> opposite edges equal (parallelogram) -> no taper -> 1
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+      /*actual=*/taper_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0)}), TOLERANCE);
+
+    // Rhombus (unit edges) -> all edges equal -> no taper -> 1
+    {
+      const Real c = std::cos(libMesh::pi/3), s = std::sin(libMesh::pi/3);
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+        /*actual=*/taper_of({Point(0,0,0), Point(1,0,0), Point(1.+c,s,0), Point(c,s,0)}), TOLERANCE);
+    }
+
+    // Symmetric trapezoid: bottom edge length 4, top edge length 2, the
+    // two slanted edges equal -> worst pair ratio = 2/4 = 0.5
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.5,
+      /*actual=*/taper_of({Point(0,0,0), Point(4,0,0), Point(3,1,0), Point(1,1,0)}), TOLERANCE);
+  }
+
   void testTet4DihedralAngle()
   {
     LOG_UNIT_TEST;

From 70eafc940c0a42fbc189cf9107cffd6dce5d10b5 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 10:08:22 -0600
Subject: [PATCH 07/16] Fix review issues in the quality metric commits

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>
---
 src/geom/cell_hex.C | 2 +-
 src/geom/elem.C     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geom/cell_hex.C b/src/geom/cell_hex.C
index 2d8f0c85330..10c27636b25 100644
--- a/src/geom/cell_hex.C
+++ b/src/geom/cell_hex.C
@@ -553,12 +553,12 @@ std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
       bounds.second = 4.;
       break;
 
-    case SKEW:
     case SKEW_ANGLE:
       bounds.first  = 0.;
       bounds.second = 0.5;
       break;
 
+    case SKEW:
     case SHEAR:
     case SHAPE:
       bounds.first  = 0.3;
diff --git a/src/geom/elem.C b/src/geom/elem.C
index 1f61c0dd422..b0f0bac05db 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -2066,7 +2066,7 @@ Real Elem::quality (const ElemQuality q) const
         if (J == 0.)
           return 0.;
 
-        return std::min(J, 1. / J);
+        return std::min(J, Real(1) / J);
       }
 
       // Return 1 if we made it here

From 872d76d91304d6b5a36ac5ee73fd4ed40b25b383 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 10:13:18 -0600
Subject: [PATCH 08/16] Implement the CONDITION quality metric

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>
---
 src/geom/elem.C         | 85 +++++++++++++++++++++++++++++++++++++++++
 src/geom/elem_quality.C |  6 ++-
 2 files changed, 89 insertions(+), 2 deletions(-)

diff --git a/src/geom/elem.C b/src/geom/elem.C
index b0f0bac05db..b50f917d6ed 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -2069,6 +2069,91 @@ Real Elem::quality (const ElemQuality q) const
         return std::min(J, Real(1) / J);
       }
 
+      // Maximum condition number of the nodal Jacobian matrix over the
+      // corner nodes. At each corner the Jacobian A has the adjacent
+      // edge vectors as its columns; its (Frobenius-norm) condition
+      // number is kappa = |A|_F * |A^{-1}|_F / N. Following the other
+      // algebraic metrics (SHAPE, SKEW), the reference (weight) matrix
+      // is the identity, so kappa = 1 for an orthogonal, equal-length
+      // (ideal) corner and grows without bound as the corner is
+      // stretched or skewed. A degenerate corner (zero Jacobian
+      // 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:
+      {
+        // 1D elements don't have interior corners, so this metric does
+        // not really apply to them.
+        const auto N = this->dim();
+        if (N < 2)
+          return 1.;
+
+        // kappa >= 1 for every matrix, so 1 is both the ideal value and
+        // a safe floor for the running maximum.
+        Real max_cond = 1.;
+
+        for (auto n : this->node_index_range())
+          {
+            // Get list of edge ids adjacent to this node.
+            auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
+
+            // Skip any nodes that don't have dim() adjacent edges (see
+            // the JACOBIAN metric above for the Pyramid apex caveat).
+            if (adjacent_edge_ids.size() != N)
+              continue;
+
+            // Construct oriented edges pointing away from node n; these
+            // are the columns of the nodal Jacobian A.
+            std::vector<Point> e(N);
+            for (auto i : make_range(N))
+              {
+                auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
+                auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
+                if (node_0 != n)
+                  std::swap(node_0, node_1);
+                e[i] = this->point(node_1) - this->point(node_0);
+              }
+
+            // Squared Frobenius norm of A.
+            Real frob_A_sq = 0.;
+            for (auto i : make_range(N))
+              frob_A_sq += e[i].norm_sq();
+
+            // |det(A)| and the squared Frobenius norm of A^{-1}.
+            Real abs_det, frob_Ainv_sq;
+            if (N == 2)
+              {
+                abs_det = cross_norm(e[0], e[1]);
+
+                // Degenerate corner: infinite condition number.
+                if (abs_det == 0.)
+                  return 0.;
+
+                // For a 2x2 matrix, |A^{-1}|_F = |A|_F / |det|.
+                frob_Ainv_sq = frob_A_sq / (abs_det * abs_det);
+              }
+            else
+              {
+                abs_det = std::abs(triple_product(e[0], e[1], e[2]));
+
+                // Degenerate corner: infinite condition number.
+                if (abs_det == 0.)
+                  return 0.;
+
+                // The rows of A^{-1} are (e1 x e2), (e2 x e0), (e0 x e1),
+                // each divided by det(A).
+                frob_Ainv_sq = (e[1].cross(e[2]).norm_sq() +
+                                e[2].cross(e[0]).norm_sq() +
+                                e[0].cross(e[1]).norm_sq()) / (abs_det * abs_det);
+              }
+
+            const Real kappa = std::sqrt(frob_A_sq * frob_Ainv_sq) / N;
+            max_cond = std::max(max_cond, kappa);
+          }
+
+        return max_cond;
+      }
+
       // Return 1 if we made it here
     default:
       {
diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index f4f06454286..629166295b7 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -249,8 +249,10 @@ std::string Quality::describe (const ElemQuality q)
       break;
 
     case CONDITION:
-      desc << "Condition number of the\n"
-           << "Jacobian matrix.\n"
+      desc << "Maximum condition number of\n"
+           << "the Jacobian matrix at each\n"
+           << "corner. 1 is ideal, larger\n"
+           << "values are worse.\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (1 -> 4)\n"

From 36f69e8a2f06fa242a2d4590ba2846c5ace6b6d3 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Mon, 14 Sep 2026 10:13:18 -0600
Subject: [PATCH 09/16] Add unit tests for the CONDITION quality metric

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>
---
 tests/geom/volume_test.C | 64 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 64 insertions(+)

diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index 588dc7e9bb7..6d1eca5ec45 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -47,6 +47,8 @@ public:
   CPPUNIT_TEST( testQuad4Size );
   CPPUNIT_TEST( testHex8Size );
   CPPUNIT_TEST( testQuad4Taper );
+  CPPUNIT_TEST( testQuad4Condition );
+  CPPUNIT_TEST( testHex8Condition );
   CPPUNIT_TEST( testTri3AspectRatio );
   CPPUNIT_TEST( testTet4DihedralAngle );
   CPPUNIT_TEST( testTet4Jacobian );
@@ -1033,6 +1035,68 @@ public:
       /*actual=*/taper_of({Point(0,0,0), Point(4,0,0), Point(3,1,0), Point(1,1,0)}), TOLERANCE);
   }
 
+  void testQuad4Condition()
+  {
+    LOG_UNIT_TEST;
+
+    // CONDITION = max over corners of the Jacobian condition number
+    // kappa = |A|_F |A^-1|_F / 2 = (|e0|^2 + |e1|^2) / (2 |e0 x e1|),
+    // which is 1 for a square corner and grows with stretch or skew.
+    auto cond_of = [this](const std::vector<Point> & pts)
+    {
+      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
+      libmesh_ignore(nodes);
+      return elem->quality(CONDITION);
+    };
+
+    // Unit square -> 1
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+      /*actual=*/cond_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0)}), TOLERANCE);
+
+    // 2x1 rectangle -> (4 + 1) / (2 * 2) = 1.25 at every corner
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.25,
+      /*actual=*/cond_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0)}), TOLERANCE);
+
+    // Unit-edge rhombus with interior angle theta -> kappa = 1/sin(theta)
+    {
+      auto rhombus_condition = [&cond_of](Real theta)
+      {
+        const Real c = std::cos(theta), s = std::sin(theta);
+        return cond_of({Point(0,0,0), Point(1,0,0), Point(1.+c,s,0), Point(c,s,0)});
+      };
+
+      // theta = pi/6 -> 1/sin(pi/6) = 2
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1./std::sin(libMesh::pi/6),
+        /*actual=*/rhombus_condition(libMesh::pi/6), TOLERANCE);
+      // theta = pi/3 -> 1/sin(pi/3) = 2/sqrt(3)
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1./std::sin(libMesh::pi/3),
+        /*actual=*/rhombus_condition(libMesh::pi/3), TOLERANCE);
+    }
+  }
+
+  void testHex8Condition()
+  {
+    LOG_UNIT_TEST;
+
+    auto cond_of = [this](const std::vector<Point> & pts)
+    {
+      auto [elem, nodes] = this->construct_elem(pts, HEX8);
+      libmesh_ignore(nodes);
+      return elem->quality(CONDITION);
+    };
+
+    // Unit cube -> 1
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
+      /*actual=*/cond_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0),
+                          Point(0,0,1), Point(1,0,1), Point(1,1,1), Point(0,1,1)}), TOLERANCE);
+
+    // 2x1x1 box: at each corner |A|_F^2 = 6, |det| = 2,
+    // |A^-1|_F^2 = 9/4, so kappa = sqrt(6 * 9/4)/3 = sqrt(6)/2
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::sqrt(Real(6))/2.,
+      /*actual=*/cond_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0),
+                          Point(0,0,1), Point(2,0,1), Point(2,1,1), Point(0,1,1)}), TOLERANCE);
+  }
+
   void testTet4DihedralAngle()
   {
     LOG_UNIT_TEST;

From d374f95834ba127bf1add9661bc43f01e7762267 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Wed, 16 Sep 2026 15:50:53 -0600
Subject: [PATCH 10/16] Fix backwards TAPER qual_bounds for quads and hexes

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>
---
 src/geom/cell_hex.C  | 7 +++++--
 src/geom/face_quad.C | 7 +++++--
 2 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/geom/cell_hex.C b/src/geom/cell_hex.C
index 10c27636b25..1348b84ebb1 100644
--- a/src/geom/cell_hex.C
+++ b/src/geom/cell_hex.C
@@ -582,8 +582,11 @@ std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
       break;
 
     case TAPER:
-      bounds.first  = 0.;
-      bounds.second = 0.4;
+      // TAPER is 1 for an untapered element and decreases toward 0 with
+      // increasing taper (see Hex::quality), so the good range runs up
+      // to 1, not down from 0.
+      bounds.first  = 0.4;
+      bounds.second = 1.;
       break;
 
     case STRETCH:
diff --git a/src/geom/face_quad.C b/src/geom/face_quad.C
index bbd084e7581..c1d0fe9cbc9 100644
--- a/src/geom/face_quad.C
+++ b/src/geom/face_quad.C
@@ -553,8 +553,11 @@ std::pair<Real, Real> Quad::qual_bounds (const ElemQuality q) const
       break;
 
     case TAPER:
-      bounds.first  = 0.;
-      bounds.second = 0.7;
+      // TAPER is 1 for an untapered element and decreases toward 0 with
+      // increasing taper (see Quad::quality), so the good range runs up
+      // to 1, not down from 0.
+      bounds.first  = 0.7;
+      bounds.second = 1.;
       break;
 
     case WARP:

From 0f909d7b18ef33c65be28518bea7e944a74a9c63 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Wed, 16 Sep 2026 15:50:53 -0600
Subject: [PATCH 11/16] Add elem_test check that valid metrics stay within
 qual_bounds

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>
---
 tests/geom/elem_test.C | 82 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 82 insertions(+)

diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C
index 943e4f8554e..8b5db38c35a 100644
--- a/tests/geom/elem_test.C
+++ b/tests/geom/elem_test.C
@@ -8,6 +8,7 @@
 #include <libmesh/parallel_implementation.h>
 #include <libmesh/enum_to_string.h>
 #include <libmesh/elem_quality.h>
+#include <libmesh/node.h>
 
 using namespace libMesh;
 
@@ -253,6 +254,86 @@ public:
       }
   }
 
+  void test_quality_bounds()
+  {
+    LOG_UNIT_TEST;
+
+    // For a well-shaped, unit-scale "ideal" element of this type, every
+    // quality metric that libMesh considers valid for the type should
+    // evaluate to a value inside its suggested qual_bounds(). We build
+    // ideal (regular) shapes only for the linear element types where
+    // that is straightforward; other types are skipped.
+    std::vector<Point> pts;
+    switch (elem_type)
+      {
+      case TRI3:
+        // Equilateral triangle, unit edge length.
+        pts = {Point(0, 0, 0),
+               Point(1, 0, 0),
+               Point(0.5, std::sqrt(Real(3))/2., 0)};
+        break;
+
+      case QUAD4:
+        // Unit square.
+        pts = {Point(0, 0, 0), Point(1, 0, 0),
+               Point(1, 1, 0), Point(0, 1, 0)};
+        break;
+
+      case TET4:
+        // Regular tetrahedron, unit edge length.
+        pts = {Point(0, 0, 0),
+               Point(1, 0, 0),
+               Point(0.5, std::sqrt(Real(3))/2., 0),
+               Point(0.5, std::sqrt(Real(3))/6., std::sqrt(Real(2)/3.))};
+        break;
+
+      case HEX8:
+        // Unit cube.
+        pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
+               Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)};
+        break;
+
+      default:
+        // No ideal element constructed for this type; nothing to check.
+        return;
+      }
+
+    // Build the element from freshly created nodes; evaluating quality
+    // metrics does not require the element to belong to a mesh.
+    std::vector<std::unique_ptr<Node>> nodes(pts.size());
+    for (unsigned int i = 0; i < pts.size(); ++i)
+      nodes[i] = Node::build(pts[i], i);
+
+    std::unique_ptr<Elem> elem = Elem::build(elem_type);
+    CPPUNIT_ASSERT(elem->n_nodes() == pts.size());
+    for (unsigned int i = 0; i < pts.size(); ++i)
+      elem->set_node(i, nodes[i].get());
+
+    for (const ElemQuality q : Quality::valid(elem_type))
+      {
+        const std::pair<Real, Real> bounds = elem->qual_bounds(q);
+
+        // A metric may be listed as valid for a type without having
+        // suggested bounds defined; qual_bounds() returns (-1, -1) as a
+        // sentinel in that case, which we skip.
+        if (bounds.first == -1. && bounds.second == -1.)
+          continue;
+
+        const Real value = elem->quality(q);
+
+        std::ostringstream msg;
+        msg << "Quality metric " << Utility::enum_to_string(q)
+            << " on an ideal " << Utility::enum_to_string(elem_type)
+            << " evaluated to " << value
+            << ", outside its suggested bounds ["
+            << bounds.first << ", " << bounds.second << "]";
+
+        CPPUNIT_ASSERT_MESSAGE(msg.str(),
+                               value >= bounds.first - TOLERANCE &&
+                               value <= bounds.second + TOLERANCE);
+      }
+  }
+
   void test_maps()
   {
     LOG_UNIT_TEST;
@@ -971,6 +1052,7 @@ public:
   CPPUNIT_TEST( test_bounding_box );            \
   CPPUNIT_TEST( test_ref_elem );                \
   CPPUNIT_TEST( test_quality );                 \
+  CPPUNIT_TEST( test_quality_bounds );          \
   CPPUNIT_TEST( test_node_edge_map_consistency ); \
   CPPUNIT_TEST( test_maps );                    \
   CPPUNIT_TEST( test_static_data );             \

From 84348704ebb6b14a43a05a66eb4eb42cea74b212 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Wed, 16 Sep 2026 20:19:09 -0600
Subject: [PATCH 12/16] Reference an ideal element (not the master) in
 CONDITION and SIZE

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>
---
 src/geom/elem.C          | 153 +++++++++++++++++++++------------------
 src/geom/elem_quality.C  |  13 ++--
 tests/geom/volume_test.C |  50 +++++++------
 3 files changed, 118 insertions(+), 98 deletions(-)

diff --git a/src/geom/elem.C b/src/geom/elem.C
index b50f917d6ed..5290ca14613 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -22,6 +22,7 @@
 #include "libmesh/boundary_info.h"
 #include "libmesh/fe_type.h"
 #include "libmesh/fe_interface.h"
+#include "libmesh/tensor_value.h"
 #include "libmesh/node_elem.h"
 #include "libmesh/edge_edge2.h"
 #include "libmesh/edge_edge3.h"
@@ -1996,23 +1997,17 @@ Real Elem::quality (const ElemQuality q) const
         return min_node_area;
       }
 
-      // Relative size metric: min(J, 1/J), where J is the determinant
-      // of the "weighted" nodal Jacobian, A * W^{-1}. Following the
-      // other algebraic metrics (SHAPE, SKEW, JACOBIAN), the reference
-      // (weight) matrix W is the identity, i.e. the canonical unit
-      // reference element (unit-length edges meeting at right angles),
-      // for which det(W) = 1. J is therefore the element's nodal
-      // Jacobian determinant (area in 2D, volume in 3D spanned by the
-      // edges meeting at a node), averaged over the corner nodes.
-      //
-      // Both undersized (J < 1) and oversized (J > 1) elements are
-      // penalized, and an element the size of the unit reference
-      // 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.
+      // Relative size metric: min over the corner nodes of min(tau,
+      // 1/tau), where tau is the ratio of the corner's nodal Jacobian
+      // determinant to the determinant of the nodal Jacobian of an
+      // ideal (regular) element of the *same volume*. We take that
+      // same-volume reference to be the mean of the element's own
+      // corner nodal Jacobian determinants, so tau = 1 at every corner
+      // of an element whose Jacobian is uniform -- i.e. any affine
+      // element (parallelogram, box, or regular simplex), at any scale
+      // -- and the metric is 1. Non-uniform (non-affine) elements, e.g.
+      // tapered or sheared shapes, score below 1, and a degenerate
+      // corner (zero determinant) drives it to 0.
     case SIZE:
       {
         // 1D elements don't have interior corners, so this metric does
@@ -2021,11 +2016,10 @@ Real Elem::quality (const ElemQuality q) const
         if (N < 2)
           return 1.;
 
-        // Average the nodal Jacobian determinant over the corner
-        // nodes. This uses the same nodal Jacobian construction as the
-        // JACOBIAN metric above.
+        // Collect the nodal Jacobian determinant at each corner, using
+        // the same construction as the JACOBIAN metric above.
+        std::vector<Real> node_areas;
         Real sum_node_area = 0.;
-        unsigned int n_corners = 0;
 
         for (auto n : this->node_index_range())
           {
@@ -2049,37 +2043,55 @@ Real Elem::quality (const ElemQuality q) const
               }
 
             // Unscaled nodal area (2D) or volume (3D).
-            Real node_area = (N == 2) ?
+            const Real node_area = (N == 2) ?
               cross_norm(oriented_edges[0], oriented_edges[1]) :
               std::abs(triple_product(oriented_edges[0], oriented_edges[1], oriented_edges[2]));
 
+            node_areas.push_back(node_area);
             sum_node_area += node_area;
-            ++n_corners;
           }
 
-        // No usable corners, or a degenerate (zero-size) element: return
-        // 0 (the lowest quality).
-        if (n_corners == 0)
+        // No usable corners: return 0 (the lowest quality).
+        if (node_areas.empty())
           return 0.;
 
-        const Real J = sum_node_area / n_corners;
-        if (J == 0.)
+        // Nodal Jacobian determinant of the ideal element of the same
+        // volume (its Jacobian is uniform, so this is the mean).
+        const Real mean_node_area = sum_node_area / node_areas.size();
+        if (mean_node_area == 0.)
           return 0.;
 
-        return std::min(J, Real(1) / J);
+        Real size = 1.;
+        for (const Real a : node_areas)
+          {
+            // A zero-area corner is degenerate: worst quality.
+            if (a == 0.)
+              return 0.;
+
+            const Real tau = a / mean_node_area;
+            size = std::min(size, std::min(tau, Real(1) / tau));
+          }
+
+        return size;
       }
 
-      // Maximum condition number of the nodal Jacobian matrix over the
-      // corner nodes. At each corner the Jacobian A has the adjacent
-      // edge vectors as its columns; its (Frobenius-norm) condition
-      // number is kappa = |A|_F * |A^{-1}|_F / N. Following the other
-      // algebraic metrics (SHAPE, SKEW), the reference (weight) matrix
-      // is the identity, so kappa = 1 for an orthogonal, equal-length
-      // (ideal) corner and grows without bound as the corner is
-      // stretched or skewed. A degenerate corner (zero Jacobian
-      // 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.
+      // Maximum condition number of the nodal Jacobian over the corner
+      // nodes, measured relative to an ideal (regular) element rather
+      // than to the reference element. At each corner the physical
+      // nodal Jacobian A has the adjacent edge vectors as its columns,
+      // and W is the nodal Jacobian of the ideal corner: unit-length
+      // edges meeting at 60 degrees for a simplex, 90 degrees
+      // otherwise. Working with the corner metric tensors T_A = A^T A
+      // and T_W = W^T W (which handles a lower-dimensional element
+      // living in a higher-dimensional space), the Frobenius condition
+      // number of the weighted Jacobian A W^{-1} is
+      //   kappa = sqrt(tr(T_A T_W^{-1}) * tr(T_W T_A^{-1})) / N.
+      // This is 1 for a corner similar to the ideal one -- so an
+      // equilateral triangle or regular tetrahedron scores 1, not just
+      // a right-angled corner -- and grows with distortion. A
+      // degenerate corner has an infinite condition number, reported as
+      // 0 following the convention that 0 stands in for infinity (cf.
+      // EDGE_LENGTH_RATIO).
     case CONDITION:
       {
         // 1D elements don't have interior corners, so this metric does
@@ -2088,6 +2100,18 @@ Real Elem::quality (const ElemQuality q) const
         if (N < 2)
           return 1.;
 
+        // Ideal corner metric tensor T_W: unit-length edges meeting at
+        // the regular angle (cos = 0.5 for simplices, 0 otherwise).
+        // Unused rows/columns are left as the identity so that
+        // RealTensor's 3x3 inverse yields the correct NxN inverse.
+        const Real cos_ideal = (this->n_vertices() == N + 1) ? 0.5 : 0.;
+        RealTensor Tw(1, 0, 0,  0, 1, 0,  0, 0, 1);
+        for (auto i : make_range(N))
+          for (auto j : make_range(N))
+            if (i != j)
+              Tw(i, j) = cos_ideal;
+        const RealTensor Tw_inv = Tw.inverse();
+
         // kappa >= 1 for every matrix, so 1 is both the ideal value and
         // a safe floor for the running maximum.
         Real max_cond = 1.;
@@ -2114,40 +2138,31 @@ Real Elem::quality (const ElemQuality q) const
                 e[i] = this->point(node_1) - this->point(node_0);
               }
 
-            // Squared Frobenius norm of A.
-            Real frob_A_sq = 0.;
+            // Physical corner metric tensor T_A = A^T A, padded with the
+            // identity in unused dimensions (as for T_W above).
+            RealTensor Ta(1, 0, 0,  0, 1, 0,  0, 0, 1);
             for (auto i : make_range(N))
-              frob_A_sq += e[i].norm_sq();
-
-            // |det(A)| and the squared Frobenius norm of A^{-1}.
-            Real abs_det, frob_Ainv_sq;
-            if (N == 2)
-              {
-                abs_det = cross_norm(e[0], e[1]);
-
-                // Degenerate corner: infinite condition number.
-                if (abs_det == 0.)
-                  return 0.;
+              for (auto j : make_range(N))
+                Ta(i, j) = e[i] * e[j];
 
-                // For a 2x2 matrix, |A^{-1}|_F = |A|_F / |det|.
-                frob_Ainv_sq = frob_A_sq / (abs_det * abs_det);
-              }
-            else
-              {
-                abs_det = std::abs(triple_product(e[0], e[1], e[2]));
+            // Degenerate corner: infinite condition number.
+            if (Ta.det() == 0.)
+              return 0.;
 
-                // Degenerate corner: infinite condition number.
-                if (abs_det == 0.)
-                  return 0.;
+            const RealTensor Ta_inv = Ta.inverse();
 
-                // The rows of A^{-1} are (e1 x e2), (e2 x e0), (e0 x e1),
-                // each divided by det(A).
-                frob_Ainv_sq = (e[1].cross(e[2]).norm_sq() +
-                                e[2].cross(e[0]).norm_sq() +
-                                e[0].cross(e[1]).norm_sq()) / (abs_det * abs_det);
-              }
+            // num1 = tr(T_A T_W^{-1}), num2 = tr(T_W T_A^{-1}) over the
+            // NxN blocks (both symmetric, so summed as elementwise dot
+            // products).
+            Real num1 = 0., num2 = 0.;
+            for (auto i : make_range(N))
+              for (auto j : make_range(N))
+                {
+                  num1 += Ta(i, j) * Tw_inv(i, j);
+                  num2 += Tw(i, j) * Ta_inv(i, j);
+                }
 
-            const Real kappa = std::sqrt(frob_A_sq * frob_Ainv_sq) / N;
+            const Real kappa = std::sqrt(num1 * num2) / N;
             max_cond = std::max(max_cond, kappa);
           }
 
diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index 629166295b7..662b1bd53a1 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -251,8 +251,9 @@ std::string Quality::describe (const ElemQuality q)
     case CONDITION:
       desc << "Maximum condition number of\n"
            << "the Jacobian matrix at each\n"
-           << "corner. 1 is ideal, larger\n"
-           << "values are worse.\n"
+           << "corner, relative to an ideal\n"
+           << "(regular) element. 1 is ideal,\n"
+           << "larger values are worse.\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (1 -> 4)\n"
@@ -338,9 +339,11 @@ std::string Quality::describe (const ElemQuality q)
 
     case SIZE:
       desc << "Relative size: min(J, 1/J),\n"
-           << "where J is the determinant\n"
-           << "of the nodal Jacobian relative\n"
-           << "to the unit reference element.\n"
+           << "where J is the determinant of\n"
+           << "the nodal Jacobian relative to\n"
+           << "an ideal element of the same\n"
+           << "volume. 1 for a uniform\n"
+           << "(affine) element.\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (0.3 -> 1)\n"
diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index 6d1eca5ec45..3170d53302e 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -942,10 +942,12 @@ public:
   {
     LOG_UNIT_TEST;
 
-    // Relative SIZE = min(J, 1/J), where J is the element's nodal
-    // Jacobian determinant measured against the unit reference
-    // element. A unit square is ideal (1); larger and smaller elements
-    // score below 1 symmetrically.
+    // Relative SIZE = min over corners of min(tau, 1/tau), where tau is
+    // the corner nodal Jacobian determinant divided by that of an ideal
+    // element of the same volume (the mean nodal determinant). It is 1
+    // for any element with a uniform Jacobian (any affine element, at
+    // any scale) and drops below 1 as the Jacobian varies across the
+    // element.
     auto size_of = [this](const std::vector<Point> & pts)
     {
       auto [elem, nodes] = this->construct_elem(pts, QUAD4);
@@ -953,25 +955,25 @@ public:
       return elem->quality(SIZE);
     };
 
-    // Unit square -> J = 1 -> 1
+    // Affine elements (uniform Jacobian) all score 1, regardless of
+    // scale or shape: unit square, a larger square, a stretched
+    // rectangle, and a sheared rhombus.
     LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
       /*actual=*/size_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0)}), TOLERANCE);
-
-    // 2x2 square -> nodal area J = 4 -> min(4, 1/4) = 0.25
-    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.25,
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
       /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,2,0), Point(0,2,0)}), TOLERANCE);
-
-    // 2x1 rectangle -> J = 2 -> 0.5
-    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.5,
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
       /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0)}), TOLERANCE);
-
-    // Unit-edge rhombus with interior angle pi/6 -> J = sin(pi/6) = 0.5
-    // -> min(0.5, 2) = 0.5
     {
       const Real c = std::cos(libMesh::pi/6), s = std::sin(libMesh::pi/6);
-      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::sin(libMesh::pi/6),
+      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
         /*actual=*/size_of({Point(0,0,0), Point(1,0,0), Point(1.+c,s,0), Point(c,s,0)}), TOLERANCE);
     }
+
+    // Trapezoid with corner nodal areas {6, 6, 4, 4} (mean 5): the
+    // worst corner ratio is 4/5 = 0.8.
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.8,
+      /*actual=*/size_of({Point(0,0,0), Point(3,0,0), Point(2,2,0), Point(0,2,0)}), TOLERANCE);
   }
 
   void testHex8Size()
@@ -985,20 +987,20 @@ public:
       return elem->quality(SIZE);
     };
 
-    // Unit cube -> J = 1 -> 1
+    // Affine boxes (uniform Jacobian) score 1 at any scale.
     LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
       /*actual=*/size_of({Point(0,0,0), Point(1,0,0), Point(1,1,0), Point(0,1,0),
                           Point(0,0,1), Point(1,0,1), Point(1,1,1), Point(0,1,1)}), TOLERANCE);
-
-    // 2x2x2 cube -> nodal volume J = 8 -> min(8, 1/8) = 0.125
-    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.125,
-      /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,2,0), Point(0,2,0),
-                          Point(0,0,2), Point(2,0,2), Point(2,2,2), Point(0,2,2)}), TOLERANCE);
-
-    // 2x1x1 box -> J = 2 -> 0.5
-    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.5,
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0,
       /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0),
                           Point(0,0,1), Point(2,0,1), Point(2,1,1), Point(0,1,1)}), TOLERANCE);
+
+    // A frustum: 2x2 base, unit top shrunk toward the axis. The four
+    // bottom corners have nodal volume 4 and the four top corners 1
+    // (mean 2.5), so the worst corner ratio is 1/2.5 = 0.4.
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.4,
+      /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,2,0), Point(0,2,0),
+                          Point(0.5,0.5,1), Point(1.5,0.5,1), Point(1.5,1.5,1), Point(0.5,1.5,1)}), TOLERANCE);
   }
 
   void testQuad4Taper()

From 40dc7fedffc087c89d94adce3fff4ad66629aca4 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Thu, 17 Sep 2026 08:23:44 -0600
Subject: [PATCH 13/16] Move get_target_elem to ReferenceElem::ideal_target

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>
---
 include/geom/reference_elem.h                 |  21 ++
 include/systems/variational_smoother_system.h |  10 -
 src/geom/reference_elem.C                     | 276 ++++++++++++++++++
 src/systems/variational_smoother_system.C     | 274 +----------------
 4 files changed, 298 insertions(+), 283 deletions(-)

diff --git a/include/geom/reference_elem.h b/include/geom/reference_elem.h
index f7ac9635e6d..09ee05ea061 100644
--- a/include/geom/reference_elem.h
+++ b/include/geom/reference_elem.h
@@ -23,11 +23,17 @@
 // Local includes
 #include "libmesh/libmesh_common.h"
 
+// C++ includes
+#include <memory>
+#include <utility>
+#include <vector>
+
 namespace libMesh
 {
 
 // forward declarations
 class Elem;
+class Node;
 enum ElemType : int;
 
 /**
@@ -46,6 +52,21 @@ namespace ReferenceElem
  */
 const Elem & get (const ElemType type_in);
 
+/**
+ * \returns A freshly built "ideal" (regular) element of the given type,
+ * i.e. the optimally-shaped element that a mesh optimizer targets: an
+ * equilateral triangle, regular tetrahedron, etc., sized to the volume
+ * of the reference element. For element types that have no distinct
+ * ideal shape (e.g. quads and hexes, whose reference element is already
+ * regular), this returns a copy of the reference element.
+ *
+ * The returned Elem holds pointers into the returned Nodes, so the
+ * caller must keep the Node vector alive for at least as long as the
+ * Elem.
+ */
+std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
+ideal_target (const ElemType type);
+
 } // namespace ReferenceElem
 
 
diff --git a/include/systems/variational_smoother_system.h b/include/systems/variational_smoother_system.h
index abdc3253b00..d8218d7378d 100644
--- a/include/systems/variational_smoother_system.h
+++ b/include/systems/variational_smoother_system.h
@@ -140,16 +140,6 @@ class VariationalSmootherSystem : public libMesh::FEMSystem
    */
   virtual void solve() override;
 
-  /**
-   * Get the target element for a given element type.
-   * @param type Element type
-   * @return a std::pair containing the target element for type and the
-   * corresponding nodes that must be kept in scope while the target element is
-   * used.
-   */
-  static std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
-  get_target_elem(const ElemType & type);
-
   /**
    * Get the jacobians (and determinants) of the target-to-reference element mapping.
    * @param target_elem Target element.
diff --git a/src/geom/reference_elem.C b/src/geom/reference_elem.C
index 1f88d7718bf..37d5be530b6 100644
--- a/src/geom/reference_elem.C
+++ b/src/geom/reference_elem.C
@@ -25,11 +25,15 @@
 #include "libmesh/threads.h"
 #include "libmesh/enum_to_string.h"
 #include "libmesh/enum_elem_type.h"
+#include "libmesh/fuzzy_equals.h"
 
 // C++ includes
 #include <map>
 #include <sstream>
 #include <memory> // std::unique_ptr
+#include <cmath>  // std::sqrt, std::cbrt
+#include <utility> // std::pair
+#include <vector>
 
 
 //-----------------------------------------------
@@ -278,5 +282,277 @@ const Elem & get (const ElemType type_in)
 
   return *ref_elem_map[base_type];
 }
+
+std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
+ideal_target (const ElemType type)
+{
+  // Build target element
+  auto target_elem = Elem::build(type);
+
+  // Volume of reference element
+  const auto ref_vol = target_elem->reference_elem()->volume();
+
+  // Update the nodes of the target element, depending on type
+  const Real sqrt_2 = std::sqrt(Real(2));
+  const Real sqrt_3 = std::sqrt(Real(3));
+  std::vector<std::unique_ptr<Node>> owned_nodes;
+
+  const auto type_str = Utility::enum_to_string(type);
+
+  // Elems deriving from Tri
+  if (type_str.compare(0, 3, "TRI") == 0)
+    {
+
+      // The target element will be an equilateral triangle with area equal to
+      // the area of the reference element.
+
+      // Equilateral triangle side length preserving area of the reference element
+      const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
+
+      // Define the nodal locations of the vertices
+      const auto & s = side_length;
+      //                                         x        y                  node_id
+      owned_nodes.emplace_back(Node::build(Point(0.,      0.),               0));
+      owned_nodes.emplace_back(Node::build(Point(s,       0.),               1));
+      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
+
+      switch (type)
+        {
+            case TRI3: {
+              // Nothing to do here, vertices already added above
+              break;
+            }
+
+            case TRI6: {
+              // Define the midpoint nodes of the equilateral triangle
+              //                                         x         y                   node_id
+              owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00),              3));
+              owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
+              owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
+
+              break;
+            }
+
+          default:
+            libmesh_error_msg("Unsupported triangular element: " << type_str);
+            break;
+        }
+    } // if Tri
+
+  // Elems deriving from Prism
+  else if (type_str.compare(0, 5, "PRISM") == 0)
+    {
+
+      // The target element will be a prism with an equilateral triangular
+      // base with volume equal to the volume of the reference element.
+
+      // For an equilateral triangular base with side length s, the
+      // base area is s^2 * sqrt(3) / 4.
+      // The prism height that will result in equal face areas is
+      // s * sqrt(3) / 4. We choose s such that the target element has
+      // the same volume as the reference element:
+      // v = (s^2 * sqrt(3) / 4) * (s * sqrt(3) / 4) = 3 * s^3 / 4
+      // --> s = (16 * v / 3)^(1/3)
+      // I have no particular motivation for imposing equal face areas,
+      // so this can be updated if a more `optimal` target prism is
+      // identified.
+
+      // Side length that preserves the volume of the reference element
+      const auto side_length = std::cbrt(16. * ref_vol / 3.);
+      // Prism height with the property that all faces have equal area
+      const auto target_height = 0.25 * side_length * sqrt_3;
+
+      const auto & s = side_length;
+      const auto & h = target_height;
+      //                                         x        y                 z    node_id
+      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
+      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
+      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
+      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               h),  3));
+      owned_nodes.emplace_back(Node::build(Point(s,       0.,               h),  4));
+      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, h),  5));
+
+      if (type == PRISM15 || type == PRISM18 || type == PRISM20 || type == PRISM21)
+        {
+          // Define the edge midpoint nodes of the prism
+          const auto & on = owned_nodes;
+          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 6));
+          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 7));
+          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 8));
+          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 9));
+          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
+          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[5]) / 2.), 11));
+          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
+          owned_nodes.emplace_back(Node::build(Point((*on[4] + *on[5]) / 2.), 13));
+          owned_nodes.emplace_back(Node::build(Point((*on[5] + *on[3]) / 2.), 14));
+
+          if (type == PRISM18 || type == PRISM20 || type == PRISM21)
+            {
+              // Define the rectangular face midpoint nodes of the prism
+              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3] + *on[4]) / 4.), 15));
+              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4] + *on[5]) / 4.), 16));
+              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3] + *on[5]) / 4.), 17));
+
+              if (type == PRISM20 || type == PRISM21)
+                {
+                  // Define the triangular face midpoint nodes of the prism
+                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 18));
+                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4] + *on[5]) / 3.), 19));
+
+                  if (type == PRISM21)
+                    // Define the interior point of the prism
+                    owned_nodes.emplace_back(Node::build(Point((*on[9] + *on[10] + *on[11]) / 3.), 20));
+
+                }
+            }
+        }
+
+      else if (type != PRISM6)
+        libmesh_error_msg("Unsupported prism element: " << type_str);
+
+    } // if Prism
+
+  // Elems deriving from Pyramid
+  else if (type_str.compare(0, 7, "PYRAMID") == 0)
+    {
+
+      // The target element is a pyramid with an square base and
+      // equilateral triangular sides with volume equal to the volume of the
+      // reference element.
+
+      // A pyramid with square base sidelength s and equilateral triangular
+      // sides has height h = s / sqrt(2).
+      // The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
+      // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
+      // non-optimal reference element.
+
+      // Side length that preserves the volume of the reference element
+      const auto side_length = std::cbrt(3. * sqrt_2 * ref_vol);
+      // Pyramid height with the property that all faces are equilateral triangles
+      const auto target_height = side_length / sqrt_2;
+
+      const auto & s = side_length;
+      const auto & h = target_height;
+
+      //                                         x        y        z    node_id
+      owned_nodes.emplace_back(Node::build(Point(0.,      0.,      0.), 0));
+      owned_nodes.emplace_back(Node::build(Point(s,       0.,      0.), 1));
+      owned_nodes.emplace_back(Node::build(Point(s,       s,       0.), 2));
+      owned_nodes.emplace_back(Node::build(Point(0.,      s,       0.), 3));
+      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h),  4));
+
+      if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
+        {
+          const auto & on = owned_nodes;
+          // Define the edge midpoint nodes of the pyramid
+
+          // Base node to base node midpoint nodes
+          owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
+          owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
+          owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
+          owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));
+
+          // Base node to apex node midpoint nodes
+          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
+          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
+          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
+          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
+
+          if (type == PYRAMID14 || type == PYRAMID18)
+            {
+              // Define the square face midpoint node of the pyramid
+              owned_nodes.emplace_back(
+                  Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));
+
+              if (type == PYRAMID18)
+                {
+                  // Define the triangular face nodes
+                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
+                  owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
+                  owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
+                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
+                }
+            }
+        }
+
+      else if (type != PYRAMID5)
+        libmesh_error_msg("Unsupported pyramid element: " << type_str);
+
+    } // if Pyramid
+
+  // Elems deriving from Tet
+  else if (type_str.compare(0, 3, "TET") == 0)
+    {
+
+      // The ideal target element is a a regular tet with equilateral
+      // triangles for all faces, with volume equal to the volume of the
+      // reference element.
+
+      // The volume of a tet is given by v = b * h / 3, where b is the area of
+      // the base face and h is the height of the apex node. The area of an
+      // equilateral triangle with side length s is b = sqrt(3) s^2 / 4.
+      // For all faces to have side length s, the height of the apex node is
+      // h = sqrt(2/3) * s. Then the volume is v = sqrt(2) * s^3 / 12.
+      // Solving for s, the side length that will preserve the volume of the
+      // reference element is s = (6 * sqrt(2) * v)^(1/3), where v is the volume
+      // of the non-optimal reference element (i.e., a right tet).
+
+      // Side length that preserves the volume of the reference element
+      const auto side_length = std::cbrt(6. * sqrt_2 * ref_vol);
+      // tet height with the property that all faces are equilateral triangles
+      const auto target_height = sqrt_2 / sqrt_3 * side_length;
+
+      const auto & s = side_length;
+      const auto & h = target_height;
+
+      // For regular tet
+      //                                         x        y                z     node_id
+      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
+      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
+      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
+      owned_nodes.emplace_back(Node::build(Point(0.5 * s, sqrt_3 / 6. * s,  h),  3));
+
+      if (type == TET10 || type == TET14)
+        {
+          const auto & on = owned_nodes;
+          // Define the edge midpoint nodes of the tet
+
+          // Base node to base node midpoint nodes
+          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 4));
+          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 5));
+          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 6));
+          // Base node to apex node midpoint nodes
+          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 7));
+          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[3]) / 2.), 8));
+          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3]) / 2.), 9));
+
+          if (type == TET14)
+            {
+              // Define the face midpoint nodes of the tet
+              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 10));
+              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3]) / 3.), 11));
+              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[3]) / 3.), 12));
+              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3]) / 3.), 13));
+            }
+        }
+
+      else if (type != TET4)
+        libmesh_error_msg("Unsupported tet element: " << type_str);
+
+    } // if Tet
+
+  // Set the target_elem equal to the reference elem
+  else
+    for (const auto & node : target_elem->reference_elem()->node_ref_range())
+      owned_nodes.emplace_back(Node::build(node, node.id()));
+
+  // Set nodes of target element
+  for (const auto & node_ptr : owned_nodes)
+    target_elem->set_node(node_ptr->id(), node_ptr.get());
+
+  libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));
+
+  return std::make_pair(std::move(target_elem), std::move(owned_nodes));
+}
 } // namespace ReferenceElem
 } // namespace libMesh
diff --git a/src/systems/variational_smoother_system.C b/src/systems/variational_smoother_system.C
index 00cd3e63f55..624eb40c6bf 100644
--- a/src/systems/variational_smoother_system.C
+++ b/src/systems/variational_smoother_system.C
@@ -254,7 +254,7 @@ void VariationalSmootherSystem::prepare_for_smoothing()
       // Add target element info, if applicable
       if (_target_jacobians.find(elem->type()) == _target_jacobians.end())
         {
-          const auto [target_elem, target_nodes] = get_target_elem(elem->type());
+          const auto [target_elem, target_nodes] = ReferenceElem::ideal_target(elem->type());
           get_target_to_reference_jacobian(target_elem.get(),
                                            femcontext,
                                            _target_jacobians[elem->type()],
@@ -897,278 +897,6 @@ void VariationalSmootherSystem::compute_mesh_quality_info()
     libMesh::out << info;
 }
 
-std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
-VariationalSmootherSystem::get_target_elem(const ElemType & type)
-{
-  // Build target element
-  auto target_elem = Elem::build(type);
-
-  // Volume of reference element
-  const auto ref_vol = target_elem->reference_elem()->volume();
-
-  // Update the nodes of the target element, depending on type
-  const Real sqrt_2 = std::sqrt(Real(2));
-  const Real sqrt_3 = std::sqrt(Real(3));
-  std::vector<std::unique_ptr<Node>> owned_nodes;
-
-  const auto type_str = Utility::enum_to_string(type);
-
-  // Elems deriving from Tri
-  if (type_str.compare(0, 3, "TRI") == 0)
-    {
-
-      // The target element will be an equilateral triangle with area equal to
-      // the area of the reference element.
-
-      // Equilateral triangle side length preserving area of the reference element
-      const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
-
-      // Define the nodal locations of the vertices
-      const auto & s = side_length;
-      //                                         x        y                  node_id
-      owned_nodes.emplace_back(Node::build(Point(0.,      0.),               0));
-      owned_nodes.emplace_back(Node::build(Point(s,       0.),               1));
-      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
-
-      switch (type)
-        {
-            case TRI3: {
-              // Nothing to do here, vertices already added above
-              break;
-            }
-
-            case TRI6: {
-              // Define the midpoint nodes of the equilateral triangle
-              //                                         x         y                   node_id
-              owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00),              3));
-              owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
-              owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
-
-              break;
-            }
-
-          default:
-            libmesh_error_msg("Unsupported triangular element: " << type_str);
-            break;
-        }
-    } // if Tri
-
-  // Elems deriving from Prism
-  else if (type_str.compare(0, 5, "PRISM") == 0)
-    {
-
-      // The target element will be a prism with an equilateral triangular
-      // base with volume equal to the volume of the reference element.
-
-      // For an equilateral triangular base with side length s, the
-      // base area is s^2 * sqrt(3) / 4.
-      // The prism height that will result in equal face areas is
-      // s * sqrt(3) / 4. We choose s such that the target element has
-      // the same volume as the reference element:
-      // v = (s^2 * sqrt(3) / 4) * (s * sqrt(3) / 4) = 3 * s^3 / 4
-      // --> s = (16 * v / 3)^(1/3)
-      // I have no particular motivation for imposing equal face areas,
-      // so this can be updated if a more `optimal` target prism is
-      // identified.
-
-      // Side length that preserves the volume of the reference element
-      const auto side_length = std::cbrt(16. * ref_vol / 3.);
-      // Prism height with the property that all faces have equal area
-      const auto target_height = 0.25 * side_length * sqrt_3;
-
-      const auto & s = side_length;
-      const auto & h = target_height;
-      //                                         x        y                 z    node_id
-      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
-      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
-      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
-      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               h),  3));
-      owned_nodes.emplace_back(Node::build(Point(s,       0.,               h),  4));
-      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, h),  5));
-
-      if (type == PRISM15 || type == PRISM18 || type == PRISM20 || type == PRISM21)
-        {
-          // Define the edge midpoint nodes of the prism
-          const auto & on = owned_nodes;
-          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 6));
-          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 7));
-          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 8));
-          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 9));
-          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
-          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[5]) / 2.), 11));
-          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
-          owned_nodes.emplace_back(Node::build(Point((*on[4] + *on[5]) / 2.), 13));
-          owned_nodes.emplace_back(Node::build(Point((*on[5] + *on[3]) / 2.), 14));
-
-          if (type == PRISM18 || type == PRISM20 || type == PRISM21)
-            {
-              // Define the rectangular face midpoint nodes of the prism
-              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3] + *on[4]) / 4.), 15));
-              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4] + *on[5]) / 4.), 16));
-              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3] + *on[5]) / 4.), 17));
-
-              if (type == PRISM20 || type == PRISM21)
-                {
-                  // Define the triangular face midpoint nodes of the prism
-                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 18));
-                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4] + *on[5]) / 3.), 19));
-
-                  if (type == PRISM21)
-                    // Define the interior point of the prism
-                    owned_nodes.emplace_back(Node::build(Point((*on[9] + *on[10] + *on[11]) / 3.), 20));
-
-                }
-            }
-        }
-
-      else if (type != PRISM6)
-        libmesh_error_msg("Unsupported prism element: " << type_str);
-
-    } // if Prism
-
-  // Elems deriving from Pyramid
-  else if (type_str.compare(0, 7, "PYRAMID") == 0)
-    {
-
-      // The target element is a pyramid with an square base and
-      // equilateral triangular sides with volume equal to the volume of the
-      // reference element.
-
-      // A pyramid with square base sidelength s and equilateral triangular
-      // sides has height h = s / sqrt(2).
-      // The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
-      // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
-      // non-optimal reference element.
-
-      // Side length that preserves the volume of the reference element
-      const auto side_length = std::cbrt(3. * sqrt_2 * ref_vol);
-      // Pyramid height with the property that all faces are equilateral triangles
-      const auto target_height = side_length / sqrt_2;
-
-      const auto & s = side_length;
-      const auto & h = target_height;
-
-      //                                         x        y        z    node_id
-      owned_nodes.emplace_back(Node::build(Point(0.,      0.,      0.), 0));
-      owned_nodes.emplace_back(Node::build(Point(s,       0.,      0.), 1));
-      owned_nodes.emplace_back(Node::build(Point(s,       s,       0.), 2));
-      owned_nodes.emplace_back(Node::build(Point(0.,      s,       0.), 3));
-      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h),  4));
-
-      if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
-        {
-          const auto & on = owned_nodes;
-          // Define the edge midpoint nodes of the pyramid
-
-          // Base node to base node midpoint nodes
-          owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
-          owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
-          owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
-          owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));
-
-          // Base node to apex node midpoint nodes
-          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
-          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
-          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
-          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
-
-          if (type == PYRAMID14 || type == PYRAMID18)
-            {
-              // Define the square face midpoint node of the pyramid
-              owned_nodes.emplace_back(
-                  Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));
-
-              if (type == PYRAMID18)
-                {
-                  // Define the triangular face nodes
-                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
-                  owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
-                  owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
-                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
-                }
-            }
-        }
-
-      else if (type != PYRAMID5)
-        libmesh_error_msg("Unsupported pyramid element: " << type_str);
-
-    } // if Pyramid
-
-  // Elems deriving from Tet
-  else if (type_str.compare(0, 3, "TET") == 0)
-    {
-
-      // The ideal target element is a a regular tet with equilateral
-      // triangles for all faces, with volume equal to the volume of the
-      // reference element.
-
-      // The volume of a tet is given by v = b * h / 3, where b is the area of
-      // the base face and h is the height of the apex node. The area of an
-      // equilateral triangle with side length s is b = sqrt(3) s^2 / 4.
-      // For all faces to have side length s, the height of the apex node is
-      // h = sqrt(2/3) * s. Then the volume is v = sqrt(2) * s^3 / 12.
-      // Solving for s, the side length that will preserve the volume of the
-      // reference element is s = (6 * sqrt(2) * v)^(1/3), where v is the volume
-      // of the non-optimal reference element (i.e., a right tet).
-
-      // Side length that preserves the volume of the reference element
-      const auto side_length = std::cbrt(6. * sqrt_2 * ref_vol);
-      // tet height with the property that all faces are equilateral triangles
-      const auto target_height = sqrt_2 / sqrt_3 * side_length;
-
-      const auto & s = side_length;
-      const auto & h = target_height;
-
-      // For regular tet
-      //                                         x        y                z     node_id
-      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
-      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
-      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
-      owned_nodes.emplace_back(Node::build(Point(0.5 * s, sqrt_3 / 6. * s,  h),  3));
-
-      if (type == TET10 || type == TET14)
-        {
-          const auto & on = owned_nodes;
-          // Define the edge midpoint nodes of the tet
-
-          // Base node to base node midpoint nodes
-          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 4));
-          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 5));
-          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 6));
-          // Base node to apex node midpoint nodes
-          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 7));
-          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[3]) / 2.), 8));
-          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3]) / 2.), 9));
-
-          if (type == TET14)
-            {
-              // Define the face midpoint nodes of the tet
-              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 10));
-              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3]) / 3.), 11));
-              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[3]) / 3.), 12));
-              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3]) / 3.), 13));
-            }
-        }
-
-      else if (type != TET4)
-        libmesh_error_msg("Unsupported tet element: " << type_str);
-
-    } // if Tet
-
-  // Set the target_elem equal to the reference elem
-  else
-    for (const auto & node : target_elem->reference_elem()->node_ref_range())
-      owned_nodes.emplace_back(Node::build(node, node.id()));
-
-  // Set nodes of target element
-  for (const auto & node_ptr : owned_nodes)
-    target_elem->set_node(node_ptr->id(), node_ptr.get());
-
-  libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));
-
-  return std::make_pair(std::move(target_elem), std::move(owned_nodes));
-}
-
 void VariationalSmootherSystem::get_target_to_reference_jacobian(
     const Elem * const target_elem,
     const FEMContext & femcontext,

From 72a3e53fa74483243e15d11e4dce92664e9ed441 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Thu, 17 Sep 2026 08:29:53 -0600
Subject: [PATCH 14/16] Use ReferenceElem::ideal_target in CONDITION and SIZE

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>
---
 src/geom/elem.C          | 184 ++++++++++++++++++---------------------
 tests/geom/volume_test.C |  10 ++-
 2 files changed, 93 insertions(+), 101 deletions(-)

diff --git a/src/geom/elem.C b/src/geom/elem.C
index 5290ca14613..2fc3dc1b1e0 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -1999,15 +1999,13 @@ Real Elem::quality (const ElemQuality q) const
 
       // Relative size metric: min over the corner nodes of min(tau,
       // 1/tau), where tau is the ratio of the corner's nodal Jacobian
-      // determinant to the determinant of the nodal Jacobian of an
-      // ideal (regular) element of the *same volume*. We take that
-      // same-volume reference to be the mean of the element's own
-      // corner nodal Jacobian determinants, so tau = 1 at every corner
-      // of an element whose Jacobian is uniform -- i.e. any affine
-      // element (parallelogram, box, or regular simplex), at any scale
-      // -- and the metric is 1. Non-uniform (non-affine) elements, e.g.
-      // tapered or sheared shapes, score below 1, and a degenerate
-      // corner (zero determinant) drives it to 0.
+      // determinant to that of the ideal (regular) element of the same
+      // volume -- ReferenceElem::ideal_target(), rescaled to this
+      // element's volume via the volume ratio. tau = 1 at every corner
+      // of an element whose Jacobian is uniform (any affine element:
+      // parallelogram, box, or regular simplex, at any scale), so the
+      // metric is 1; non-uniform (tapered/sheared) elements score below
+      // 1, and a degenerate corner drives it to 0.
     case SIZE:
       {
         // 1D elements don't have interior corners, so this metric does
@@ -2016,82 +2014,75 @@ Real Elem::quality (const ElemQuality q) const
         if (N < 2)
           return 1.;
 
-        // Collect the nodal Jacobian determinant at each corner, using
-        // the same construction as the JACOBIAN metric above.
-        std::vector<Real> node_areas;
-        Real sum_node_area = 0.;
+        const Real vol = this->volume();
+        if (vol == 0.)
+          return 0.;
+
+        // Ideal (regular) element of the same type; its Jacobian is
+        // uniform. We compare nodal determinants after rescaling it to
+        // this element's volume, i.e. multiply by ideal_vol / this_vol.
+        const auto ideal_pair = ReferenceElem::ideal_target(this->type());
+        const Elem & ideal = *ideal_pair.first;
+        const Real vol_ratio = ideal.volume() / vol;
+
+        // Nodal Jacobian determinant at node n of element el (the same
+        // construction as the JACOBIAN metric above).
+        auto nodal_det = [](const Elem & el, const unsigned int n,
+                            const std::vector<unsigned int> & edge_ids,
+                            const unsigned int dim)
+        {
+          std::vector<Point> e(dim);
+          for (unsigned int i = 0; i != dim; ++i)
+            {
+              auto n0 = el.local_edge_node(edge_ids[i], 0);
+              auto n1 = el.local_edge_node(edge_ids[i], 1);
+              if (n0 != n)
+                std::swap(n0, n1);
+              e[i] = el.point(n1) - el.point(n0);
+            }
+          return (dim == 2) ? cross_norm(e[0], e[1])
+                            : std::abs(triple_product(e[0], e[1], e[2]));
+        };
 
+        Real size = 1.;
+        bool have_corner = false;
         for (auto n : this->node_index_range())
           {
-            // Get list of edge ids adjacent to this node.
-            auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
-
             // Skip any nodes that don't have dim() adjacent edges (see
             // the JACOBIAN metric above for the Pyramid apex caveat).
+            const auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
             if (adjacent_edge_ids.size() != N)
               continue;
+            have_corner = true;
 
-            // Construct oriented edges pointing away from node n.
-            std::vector<Point> oriented_edges(N);
-            for (auto i : make_range(N))
-              {
-                auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
-                auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
-                if (node_0 != n)
-                  std::swap(node_0, node_1);
-                oriented_edges[i] = this->point(node_1) - this->point(node_0);
-              }
-
-            // Unscaled nodal area (2D) or volume (3D).
-            const Real node_area = (N == 2) ?
-              cross_norm(oriented_edges[0], oriented_edges[1]) :
-              std::abs(triple_product(oriented_edges[0], oriented_edges[1], oriented_edges[2]));
-
-            node_areas.push_back(node_area);
-            sum_node_area += node_area;
-          }
+            const Real a  = nodal_det(*this, n, adjacent_edge_ids, N);
+            const Real aw = nodal_det(ideal, n, adjacent_edge_ids, N);
 
-        // No usable corners: return 0 (the lowest quality).
-        if (node_areas.empty())
-          return 0.;
-
-        // Nodal Jacobian determinant of the ideal element of the same
-        // volume (its Jacobian is uniform, so this is the mean).
-        const Real mean_node_area = sum_node_area / node_areas.size();
-        if (mean_node_area == 0.)
-          return 0.;
-
-        Real size = 1.;
-        for (const Real a : node_areas)
-          {
-            // A zero-area corner is degenerate: worst quality.
-            if (a == 0.)
+            // Degenerate corner: worst quality.
+            if (a == 0. || aw == 0.)
               return 0.;
 
-            const Real tau = a / mean_node_area;
+            const Real tau = (a / aw) * vol_ratio;
             size = std::min(size, std::min(tau, Real(1) / tau));
           }
 
-        return size;
+        return have_corner ? size : 0.;
       }
 
       // Maximum condition number of the nodal Jacobian over the corner
-      // nodes, measured relative to an ideal (regular) element rather
-      // than to the reference element. At each corner the physical
-      // nodal Jacobian A has the adjacent edge vectors as its columns,
-      // and W is the nodal Jacobian of the ideal corner: unit-length
-      // edges meeting at 60 degrees for a simplex, 90 degrees
-      // otherwise. Working with the corner metric tensors T_A = A^T A
-      // and T_W = W^T W (which handles a lower-dimensional element
-      // living in a higher-dimensional space), the Frobenius condition
-      // number of the weighted Jacobian A W^{-1} is
+      // nodes, measured against the ideal (regular) element rather than
+      // the reference element. At each corner the physical nodal
+      // Jacobian A and the ideal nodal Jacobian W (taken from the same
+      // corner of ReferenceElem::ideal_target) give the weighted
+      // Jacobian A W^{-1}, whose Frobenius condition number, via the
+      // corner metric tensors T_A = A^T A and T_W = W^T W (which also
+      // handles a lower-dimensional element embedded in 3D), is
       //   kappa = sqrt(tr(T_A T_W^{-1}) * tr(T_W T_A^{-1})) / N.
       // This is 1 for a corner similar to the ideal one -- so an
       // equilateral triangle or regular tetrahedron scores 1, not just
       // a right-angled corner -- and grows with distortion. A
       // degenerate corner has an infinite condition number, reported as
-      // 0 following the convention that 0 stands in for infinity (cf.
-      // EDGE_LENGTH_RATIO).
+      // 0 (0 stands in for infinity, cf. EDGE_LENGTH_RATIO).
     case CONDITION:
       {
         // 1D elements don't have interior corners, so this metric does
@@ -2100,17 +2091,34 @@ Real Elem::quality (const ElemQuality q) const
         if (N < 2)
           return 1.;
 
-        // Ideal corner metric tensor T_W: unit-length edges meeting at
-        // the regular angle (cos = 0.5 for simplices, 0 otherwise).
-        // Unused rows/columns are left as the identity so that
-        // RealTensor's 3x3 inverse yields the correct NxN inverse.
-        const Real cos_ideal = (this->n_vertices() == N + 1) ? 0.5 : 0.;
-        RealTensor Tw(1, 0, 0,  0, 1, 0,  0, 0, 1);
-        for (auto i : make_range(N))
-          for (auto j : make_range(N))
-            if (i != j)
-              Tw(i, j) = cos_ideal;
-        const RealTensor Tw_inv = Tw.inverse();
+        // Ideal (regular) element of the same type; W is its nodal
+        // Jacobian. Its scale is irrelevant here (the condition number
+        // is scale invariant), so the reference-volume sizing is fine.
+        const auto ideal_pair = ReferenceElem::ideal_target(this->type());
+        const Elem & ideal = *ideal_pair.first;
+
+        // Corner metric tensor T = A^T A at node n of element el, padded
+        // with the identity in unused dimensions so that RealTensor's
+        // 3x3 inverse yields the correct NxN inverse.
+        auto metric_tensor = [](const Elem & el, const unsigned int n,
+                                const std::vector<unsigned int> & edge_ids,
+                                const unsigned int dim)
+        {
+          std::vector<Point> e(dim);
+          for (unsigned int i = 0; i != dim; ++i)
+            {
+              auto n0 = el.local_edge_node(edge_ids[i], 0);
+              auto n1 = el.local_edge_node(edge_ids[i], 1);
+              if (n0 != n)
+                std::swap(n0, n1);
+              e[i] = el.point(n1) - el.point(n0);
+            }
+          RealTensor T(1, 0, 0,  0, 1, 0,  0, 0, 1);
+          for (unsigned int i = 0; i != dim; ++i)
+            for (unsigned int j = 0; j != dim; ++j)
+              T(i, j) = e[i] * e[j];
+          return T;
+        };
 
         // kappa >= 1 for every matrix, so 1 is both the ideal value and
         // a safe floor for the running maximum.
@@ -2118,38 +2126,21 @@ Real Elem::quality (const ElemQuality q) const
 
         for (auto n : this->node_index_range())
           {
-            // Get list of edge ids adjacent to this node.
-            auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
-
             // Skip any nodes that don't have dim() adjacent edges (see
             // the JACOBIAN metric above for the Pyramid apex caveat).
+            const auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
             if (adjacent_edge_ids.size() != N)
               continue;
 
-            // Construct oriented edges pointing away from node n; these
-            // are the columns of the nodal Jacobian A.
-            std::vector<Point> e(N);
-            for (auto i : make_range(N))
-              {
-                auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
-                auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
-                if (node_0 != n)
-                  std::swap(node_0, node_1);
-                e[i] = this->point(node_1) - this->point(node_0);
-              }
-
-            // Physical corner metric tensor T_A = A^T A, padded with the
-            // identity in unused dimensions (as for T_W above).
-            RealTensor Ta(1, 0, 0,  0, 1, 0,  0, 0, 1);
-            for (auto i : make_range(N))
-              for (auto j : make_range(N))
-                Ta(i, j) = e[i] * e[j];
+            const RealTensor Ta = metric_tensor(*this, n, adjacent_edge_ids, N);
 
             // Degenerate corner: infinite condition number.
             if (Ta.det() == 0.)
               return 0.;
 
+            const RealTensor Tw = metric_tensor(ideal, n, adjacent_edge_ids, N);
             const RealTensor Ta_inv = Ta.inverse();
+            const RealTensor Tw_inv = Tw.inverse();
 
             // num1 = tr(T_A T_W^{-1}), num2 = tr(T_W T_A^{-1}) over the
             // NxN blocks (both symmetric, so summed as elementwise dot
@@ -2162,8 +2153,7 @@ Real Elem::quality (const ElemQuality q) const
                   num2 += Tw(i, j) * Ta_inv(i, j);
                 }
 
-            const Real kappa = std::sqrt(num1 * num2) / N;
-            max_cond = std::max(max_cond, kappa);
+            max_cond = std::max(max_cond, std::sqrt(num1 * num2) / N);
           }
 
         return max_cond;
diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index 3170d53302e..0605b5aa405 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -995,10 +995,12 @@ public:
       /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,1,0), Point(0,1,0),
                           Point(0,0,1), Point(2,0,1), Point(2,1,1), Point(0,1,1)}), TOLERANCE);
 
-    // A frustum: 2x2 base, unit top shrunk toward the axis. The four
-    // bottom corners have nodal volume 4 and the four top corners 1
-    // (mean 2.5), so the worst corner ratio is 1/2.5 = 0.4.
-    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/0.4,
+    // A frustum: 2x2 base, unit top shrunk toward the axis. Bottom
+    // corners have nodal volume 4, top corners 1; the element volume is
+    // 7/3 and the ideal cube's nodal-det/volume ratio is 1, so the ideal
+    // nodal volume at this element's volume is 7/3. The worst corner
+    // ratio is (1)/(7/3) = 3/7.
+    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(3)/7,
       /*actual=*/size_of({Point(0,0,0), Point(2,0,0), Point(2,2,0), Point(0,2,0),
                           Point(0.5,0.5,1), Point(1.5,0.5,1), Point(1.5,1.5,1), Point(0.5,1.5,1)}), TOLERANCE);
   }

From 414c695be2efa6d47126ba32d77be701e7ba2fd1 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Thu, 17 Sep 2026 08:45:07 -0600
Subject: [PATCH 15/16] Add warning to metric description

---
 src/geom/elem_quality.C | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index 662b1bd53a1..d3e3ac699ed 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -180,7 +180,8 @@ std::string Quality::describe (const ElemQuality q)
       desc << "Maximum |cos A|, where A\n"
            << "is the angle between edges\n"
            << "at element center.\n"
-           << '\n'
+           << "NOTE: some degenerate elements\n"
+           << "score 0 if zero-length along principal axis.\n"
            << "Suggested ranges:\n"
            << "Hexes: (0 -> 0.5)\n"
            << "Quads: (0 -> 0.5)";

From a41acb4ebf52cbc7401328f0421d4301fc42a7c7 Mon Sep 17 00:00:00 2001
From: Guillaume Giudicelli <guillaume.giudicelli@inl.gov>
Date: Thu, 17 Sep 2026 12:06:31 -0600
Subject: [PATCH 16/16] Register EDGE_LENGTH_RATIO and SCALED_JACOBIAN in the
 ElemQuality enum 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>
---
 src/utils/string_to_enum.C | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/utils/string_to_enum.C b/src/utils/string_to_enum.C
index 7f911244a07..f0a03bb63df 100644
--- a/src/utils/string_to_enum.C
+++ b/src/utils/string_to_enum.C
@@ -427,7 +427,9 @@ std::map<std::string, ElemQuality> elemquality_to_enum {
    {"ASPECT_RATIO_GAMMA" , ASPECT_RATIO_GAMMA},
    {"SIZE"               , SIZE},
    {"JACOBIAN"           , JACOBIAN},
+   {"SCALED_JACOBIAN"    , SCALED_JACOBIAN},
    {"TWIST"              , TWIST},
+   {"EDGE_LENGTH_RATIO"  , EDGE_LENGTH_RATIO},
   };
 
 std::map<ElemQuality, std::string> enum_to_elemquality =