C0polygon/Polyhedron Element Quality Metrics

include/geom/elem.h

@@ -1087,11 +1087,14 @@ class Elem : public ReferenceCountedObject<Elem>,
    * .) EDGE_LENGTH_RATIO - ratio of maximum to minimum edge (in 2D,
    *    side) length, where the min/max is taken over all Elem edges.
    * .) MIN,MAX_ANGLE - The minimum (respectively maximum) angle
-   *    between all pairs of adjacent Elem edges, in degrees. Note
-   *    that, in 3D, these are *not* the dihedral angles (angle
-   *    between adjacent planar faces of the element), which we plan to
-   *    add support for in the future. In 2D, we compute the angle
-   *    between adjacent sides for this metric.
+   *    between all pairs of adjacent Elem edges, in degrees. In 3D,
+   *    these are *not* the dihedral angles between adjacent planar
+   *    faces of the element. In 2D, we compute the angle between
+   *    adjacent sides for this metric.
+   * .) MIN,MAX_DIHEDRAL_ANGLE - In 3D, the minimum (respectively
+   *    maximum) unoriented angle between adjacent side planes, folded
+   *    into the range [0, 90] degrees. In 2D, these are equivalent to
+   *    MIN,MAX_ANGLE.
    */
   virtual Real quality (const ElemQuality q) const;
 

src/geom/cell_polyhedron.C

@@ -494,6 +494,19 @@ std::pair<Real, Real> Polyhedron::qual_bounds (const ElemQuality q) const
       bounds.second = 180.;
       break;
 
+    // The generic dihedral metrics measure the unoriented angle
+    // between adjacent face planes, so their values are in [0, 90].
+    // Use the same recommended lower bounds as MIN,MAX_ANGLE.
+    case MIN_DIHEDRAL_ANGLE:
+      bounds.first  = 30.;
+      bounds.second = 90.;
+      break;
+
+    case MAX_DIHEDRAL_ANGLE:
+      bounds.first  = 60.;
+      bounds.second = 90.;
+      break;
+
     case JACOBIAN:
     case SCALED_JACOBIAN:
       bounds.first  = 0.5;

src/geom/elem_quality.C

@@ -46,6 +46,10 @@ std::string Quality::name (const ElemQuality q)
   switch (q)
     {
 
+    case EDGE_LENGTH_RATIO:
+      its_name = "Edge Length Ratio";
+      break;
+
     case ASPECT_RATIO:
       its_name = "Aspect Ratio";
       break;
@@ -211,18 +215,22 @@ std::string Quality::describe (const ElemQuality q)
 
     case MAX_DIHEDRAL_ANGLE:
       desc << "Largest angle between all adjacent pairs of sides (in 2D, equivalent to MAX_ANGLE).\n"
+           << "In 3D, this is the largest unoriented angle between adjacent side planes, in the range [0, 90].\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (90 -> 135)\n"
-           << "Triangles: (60 -> 90)";
+           << "Triangles: (60 -> 90)\n"
+           << "C0Polyhedra: (60 -> 90)";
       break;
 
     case MIN_DIHEDRAL_ANGLE:
       desc << "Smallest angle between all adjacent pairs of sides (in 2D, equivalent to MIN_ANGLE).\n"
+           << "In 3D, this is the smallest unoriented angle between adjacent side planes, in the range [0, 90].\n"
            << '\n'
            << "Suggested ranges:\n"
            << "Quads: (45 -> 90)\n"
-           << "Triangles: (30 -> 60)";
+           << "Triangles: (30 -> 60)\n"
+           << "C0Polyhedra: (30 -> 90)";
       break;
 
     case CONDITION:
@@ -490,7 +498,34 @@ std::vector<ElemQuality> Quality::valid(const ElemType t)
         break;
       }
 
+    case C0POLYGON:
+      {
+        v = {
+          EDGE_LENGTH_RATIO,
+          JACOBIAN,
+          SCALED_JACOBIAN,
+          MAX_ANGLE,
+          MIN_ANGLE
+        };
 
+        break;
+      }
+
+    case C0POLYHEDRON:
+      {
+        // The generic Jacobian metrics only inspect vertices with
+        // exactly three adjacent edges, but arbitrary polyhedra may
+        // have vertices of higher valence.
+        v = {
+          EDGE_LENGTH_RATIO,
+          MAX_ANGLE,
+          MIN_ANGLE,
+          MAX_DIHEDRAL_ANGLE,
+          MIN_DIHEDRAL_ANGLE
+        };
+
+        break;
+      }
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 

tests/geom/elem_test.C

@@ -7,6 +7,7 @@
 #include <libmesh/mesh_refinement.h>
 #include <libmesh/parallel_implementation.h>
 #include <libmesh/enum_to_string.h>
+#include <libmesh/elem_quality.h>
 
 using namespace libMesh;
 
@@ -65,8 +66,9 @@ public:
         // We're building isotropic meshes, where even elements
         // dissected from cubes ought to have tolerable quality.
         //
-        // Worst I see is 2 on tets, but let's add a little tolerance
-        // in case we decide to play with rotated meshes here later
+        // Worst I see is 2 on tets and the skewed C0Polyhedron, but
+        // let's add a little tolerance in case we decide to play with
+        // rotated meshes here later.
         CPPUNIT_ASSERT_LESSEQUAL(Real(2+TOLERANCE), edge_length_ratio); // edge_length_ratio <= 2
 
         // The MIN_ANGLE and MAX_ANGLE quality metrics are also defined on all elements
@@ -84,6 +86,64 @@ public:
         // that has a 135 degree angle
         CPPUNIT_ASSERT_LESSEQUAL(135 + TOLERANCE, max_angle);
 
+        const auto assert_quality =
+          [elem](const ElemQuality quality,
+                 const Real expected_value,
+                 const Real expected_lower_bound,
+                 const Real expected_upper_bound)
+          {
+            const auto bounds = elem->qual_bounds(quality);
+            LIBMESH_ASSERT_FP_EQUAL(expected_lower_bound, bounds.first, TOLERANCE);
+            LIBMESH_ASSERT_FP_EQUAL(expected_upper_bound, bounds.second, TOLERANCE);
+            LIBMESH_ASSERT_FP_EQUAL(expected_value, elem->quality(quality), TOLERANCE);
+          };
+
+        if (elem->type() == C0POLYGON)
+          {
+            const std::vector<ElemQuality> expected = {
+              EDGE_LENGTH_RATIO,
+              JACOBIAN,
+              SCALED_JACOBIAN,
+              MAX_ANGLE,
+              MIN_ANGLE
+            };
+
+            CPPUNIT_ASSERT(expected == Quality::valid(elem->type()));
+            CPPUNIT_ASSERT_EQUAL(
+              std::string("Edge Length Ratio"),
+              Quality::name(EDGE_LENGTH_RATIO));
+
+            assert_quality(EDGE_LENGTH_RATIO, std::sqrt(Real(2)), 1., 4.);
+            assert_quality(JACOBIAN, 0.5, 0.5, 1.);
+            assert_quality(
+              SCALED_JACOBIAN, Real(1) / std::sqrt(Real(2)), 0.5, 1.);
+            assert_quality(MAX_ANGLE, 135., 108., 144.);
+            assert_quality(MIN_ANGLE, 90., 54., 108.);
+          }
+        else if (elem->type() == C0POLYHEDRON)
+          {
+            const std::vector<ElemQuality> expected = {
+              EDGE_LENGTH_RATIO,
+              MAX_ANGLE,
+              MIN_ANGLE,
+              MAX_DIHEDRAL_ANGLE,
+              MIN_DIHEDRAL_ANGLE
+            };
+
+            CPPUNIT_ASSERT(expected == Quality::valid(elem->type()));
+
+            const Real min_polyhedron_angle =
+              std::acos(Real(1) / std::sqrt(Real(17))) *
+              Real(180) / libMesh::pi;
+            assert_quality(EDGE_LENGTH_RATIO, 2., 1., 4.);
+            assert_quality(
+              MAX_ANGLE, Real(180) - min_polyhedron_angle, 60., 180.);
+            assert_quality(MIN_ANGLE, min_polyhedron_angle, 30., 180.);
+            assert_quality(MAX_DIHEDRAL_ANGLE, 90., 60., 90.);
+            assert_quality(
+              MIN_DIHEDRAL_ANGLE, min_polyhedron_angle, 30., 90.);
+          }
+
         // Notes on minimum angle we expect to see:
         // 1.) 1D Elements don't have interior angles, so the base
         //     class implementation currently returns 0 for those
@@ -96,6 +156,7 @@ public:
           CPPUNIT_ASSERT_GREATEREQUAL((std::acos(Real(2)/std::sqrt(Real(6))) * 180 / libMesh::pi) - TOLERANCE, min_angle);
 
         // MIN,MAX_DIHEDRAL_ANGLE are implemented for all 3D elements
+        // tested here.
         if (elem->dim() > 2)
           {
             const Real min_dihedral_angle = elem->quality(MIN_DIHEDRAL_ANGLE);
@@ -120,46 +181,50 @@ public:
             CPPUNIT_ASSERT_GREATEREQUAL(45 - TOLERANCE, min_dihedral_angle);
           }
 
-        // The JACOBIAN and SCALED_JACOBIAN quality metrics are also defined on all elements
+        // The generic JACOBIAN metrics require exactly dim() adjacent
+        // edges at a vertex, which is not guaranteed for arbitrary
+        // C0Polyhedron topology.
+        if (elem->type() != C0POLYHEDRON)
+          {
+            // The largest non-infinite elements in this test (Hexes and
+            // Prisms) have a nodal Jacobian (volume) of 8, e.g. the 2x2x2
+            // reference Hex. The infinite elements that we construct for
+            // this testing have a max nodal area of 64 since they are
+            // created (see setUp()) with max side lengths of 4 (4*4*4=64).
+            const Real jac = elem->quality(JACOBIAN);
+            if (elem->infinite())
+              CPPUNIT_ASSERT_LESSEQUAL   (64 + TOLERANCE, jac);
+            else
+              CPPUNIT_ASSERT_LESSEQUAL   (8 + TOLERANCE, jac);
+
+            // The smallest 2D/3D nodal areas for regular elements in this
+            // test are found in tetrahedra, which have a minimum value of
+            // 2.0. However, we return a default value of 1.0 for 0D and
+            // 1D elements here, and we have a custom distorted-pentagon
+            // C0Polygon with a 0.5 at 3 nodes, so we handle those cases
+            // too.
+            if (elem->dim() < 2)
+              CPPUNIT_ASSERT_GREATEREQUAL(1 - TOLERANCE, jac);
+            else if (elem->type() == C0POLYGON)
+              CPPUNIT_ASSERT_GREATEREQUAL(0.5 - TOLERANCE, jac);
+            else
+              CPPUNIT_ASSERT_GREATEREQUAL(2 - TOLERANCE, jac);
+
+            // The scaled Jacobian should always be <= 1. The minimum
+            // scaled Jacobian value I observed was ~0.408248 for the
+            // tetrahedral meshes. This is consistent with the way that
+            // we generate Tet meshes with build_cube(), as we don't
+            // simply refine the reference tetrahedron in that case.
+            const Real scaled_jac = elem->quality(SCALED_JACOBIAN);
+            CPPUNIT_ASSERT_LESSEQUAL   (1 + TOLERANCE, scaled_jac);
+            CPPUNIT_ASSERT_GREATEREQUAL(    Real(0.4), scaled_jac);
 
-        // The largest non-infinite elements in this test (Hexes and
-        // Prisms) have a nodal Jacobian (volume) of 8, e.g. the 2x2x2
-        // reference Hex. The infinite elements that we construct for
-        // this testing have a max nodal area of 64 since they are
-        // created (see setUp()) with max side lengths of 4 (4*4*4=64).
-        const Real jac = elem->quality(JACOBIAN);
-        if (elem->infinite())
-          CPPUNIT_ASSERT_LESSEQUAL   (64 + TOLERANCE, jac);
-        else
-          CPPUNIT_ASSERT_LESSEQUAL   (8 + TOLERANCE, jac);
-
-        // The smallest 2D/3D nodal areas for regular elements in this
-        // test are found in tetrahedra, which have a minimum value of
-        // 2.0. However, we return a default value of 1.0 for 0D and
-        // 1D elements here, and we have a custom distorted-pentagon
-        // C0Polygon with a 0.5 at 2 nodes and a custom C0Polyhedron
-        // with a 1, so we handle those cases too.
-        if (elem->dim() < 2 || elem->type() == C0POLYHEDRON)
-          CPPUNIT_ASSERT_GREATEREQUAL(1 - TOLERANCE, jac);
-        else if (elem->type() == C0POLYGON)
-          CPPUNIT_ASSERT_GREATEREQUAL(0.5 - TOLERANCE, jac);
-        else
-          CPPUNIT_ASSERT_GREATEREQUAL(2 - TOLERANCE, jac);
-
-        // The scaled Jacobian should always be <= 1. The minimum
-        // scaled Jacobian value I observed was ~0.408248 for the
-        // tetrahedral meshes. This is consistent with the way that
-        // we generate Tet meshes with build_cube(), as we don't
-        // simply refine the reference tetrahedron in that case.
-        const Real scaled_jac = elem->quality(SCALED_JACOBIAN);
-        CPPUNIT_ASSERT_LESSEQUAL   (1 + TOLERANCE, scaled_jac);
-        CPPUNIT_ASSERT_GREATEREQUAL(    Real(0.4), scaled_jac);
-
-        // Debugging
-        // libMesh::out << "elem->type() = " << Utility::enum_to_string(elem->type())
-        //              << ", jac = " << jac
-        //              << ", scaled_jac = " << scaled_jac
-        //              << std::endl;
+            // Debugging
+            // libMesh::out << "elem->type() = " << Utility::enum_to_string(elem->type())
+            //              << ", jac = " << jac
+            //              << ", scaled_jac = " << scaled_jac
+            //              << std::endl;
+          }
       }
   }
 
@@ -935,7 +1000,7 @@ INSTANTIATE_ELEMTEST(QUADSHELL9);
 // This just tests with one pentagon, but better than nothing
 INSTANTIATE_ELEMTEST(C0POLYGON);
 
-// And this is just one truncated cube; also better than nothing
+// And this is just one skewed box; also better than nothing
 INSTANTIATE_ELEMTEST(C0POLYHEDRON);
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

tests/geom/elem_test.h

@@ -214,18 +214,20 @@ class PerElemTest : public CppUnit::TestCase
 #endif
 
 #if 1
-        // Or try a plain cube, for simpler debugging
+        // Or try an affine-skewed box.  Unlike a unit cube, its
+        // non-unit quality values cannot be confused with the generic
+        // fallback value of 1.
         _mesh->add_point(Point(0, 0, 0), 0);
-        _mesh->add_point(Point(1, 0, 0), 1);
-        _mesh->add_point(Point(1, 1, 0), 2);
-        _mesh->add_point(Point(0, 1, 0), 3);
+        _mesh->add_point(Point(2, 0, 0), 1);
+        _mesh->add_point(Point(2.25, 1, 0), 2);
+        _mesh->add_point(Point(0.25, 1, 0), 3);
         _mesh->add_point(Point(0, 0, 1), 4);
-        _mesh->add_point(Point(1, 0, 1), 5);
-        _mesh->add_point(Point(1, 1, 1), 6);
-        _mesh->add_point(Point(0, 1, 1), 7);
+        _mesh->add_point(Point(2, 0, 1), 5);
+        _mesh->add_point(Point(2.25, 1, 1), 6);
+        _mesh->add_point(Point(0.25, 1, 1), 7);
 
         // With some combinations of face triangulations, even a
-        // simple cube has no tetrahedralization that doesn't have
+        // simple box has no tetrahedralization that doesn't have
         // either interior discontinuities or a 0-volume pancake tet!
         //
         // The initial "natural" way to orient our sides is commented