Add quality metrics for quads (and all elements for some)

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

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
 
 

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.

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:

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:
       {

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
         };