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/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/cell_hex.C b/src/geom/cell_hex.C
index e1c2bab0198..1348b84ebb1 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
 
       /**
@@ -516,11 +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;
@@ -544,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/elem.C b/src/geom/elem.C
index 9440ddd9307..2fc3dc1b1e0 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,6 +1997,168 @@ Real Elem::quality (const ElemQuality q) const
         return min_node_area;
       }
 
+      // 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 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
+        // not really apply to them.
+        const auto N = this->dim();
+        if (N < 2)
+          return 1.;
+
+        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())
+          {
+            // 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;
+
+            const Real a  = nodal_det(*this, n, adjacent_edge_ids, N);
+            const Real aw = nodal_det(ideal, n, adjacent_edge_ids, N);
+
+            // Degenerate corner: worst quality.
+            if (a == 0. || aw == 0.)
+              return 0.;
+
+            const Real tau = (a / aw) * vol_ratio;
+            size = std::min(size, std::min(tau, Real(1) / tau));
+          }
+
+        return have_corner ? size : 0.;
+      }
+
+      // Maximum condition number of the nodal Jacobian over the corner
+      // 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 (0 stands in for infinity, cf. EDGE_LENGTH_RATIO).
+    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.;
+
+        // 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.
+        Real max_cond = 1.;
+
+        for (auto n : this->node_index_range())
+          {
+            // 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;
+
+            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
+            // 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);
+                }
+
+            max_cond = std::max(max_cond, std::sqrt(num1 * num2) / N);
+          }
+
+        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 2353fb72d6e..d3e3ac699ed 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,10 +166,22 @@ 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"
-           << '\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)";
@@ -234,8 +250,11 @@ 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, relative to an ideal\n"
+           << "(regular) element. 1 is ideal,\n"
+           << "larger values are worse.\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (1 -> 4)\n"
@@ -320,9 +339,12 @@ 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 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"
@@ -410,6 +432,7 @@ std::vector<ElemQuality> Quality::valid(const ElemType t)
           SHEAR,
           SIZE,
           SKEW,
+          SKEW_ANGLE,
           STRETCH,
           TAPER,
           WARP
@@ -459,6 +482,7 @@ std::vector<ElemQuality> Quality::valid(const ElemType t)
           SHEAR,
           SIZE,
           SKEW,
+          SKEW_ANGLE,
           STRETCH,
           TAPER
         };
diff --git a/src/geom/face_quad.C b/src/geom/face_quad.C
index c4ca8a57a11..c1d0fe9cbc9 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:
       {
@@ -400,6 +427,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
@@ -498,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:
@@ -541,6 +599,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/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,
diff --git a/src/utils/string_to_enum.C b/src/utils/string_to_enum.C
index 18e86cd1c6e..f0a03bb63df 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},
@@ -426,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 =
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 );             \
diff --git a/tests/geom/volume_test.C b/tests/geom/volume_test.C
index fa7e9ce995d..0605b5aa405 100644
--- a/tests/geom/volume_test.C
+++ b/tests/geom/volume_test.C
@@ -42,6 +42,13 @@ public:
   CPPUNIT_TEST( testQuad4Warpage );
   CPPUNIT_TEST( testQuad4MinMaxAngle );
   CPPUNIT_TEST( testQuad4Jacobian );
+  CPPUNIT_TEST( testQuad4SkewAngle );
+  CPPUNIT_TEST( testHex8SkewAngle );
+  CPPUNIT_TEST( testQuad4Size );
+  CPPUNIT_TEST( testHex8Size );
+  CPPUNIT_TEST( testQuad4Taper );
+  CPPUNIT_TEST( testQuad4Condition );
+  CPPUNIT_TEST( testHex8Condition );
   CPPUNIT_TEST( testTri3AspectRatio );
   CPPUNIT_TEST( testTet4DihedralAngle );
   CPPUNIT_TEST( testTet4Jacobian );
@@ -823,6 +830,277 @@ 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 testQuad4Size()
+  {
+    LOG_UNIT_TEST;
+
+    // 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);
+      libmesh_ignore(nodes);
+      return elem->quality(SIZE);
+    };
+
+    // 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);
+    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);
+    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);
+    {
+      const Real c = std::cos(libMesh::pi/6), s = 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()
+  {
+    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);
+    };
+
+    // 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);
+    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. 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);
+  }
+
+  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 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;