diff --git a/include/geom/elem.h b/include/geom/elem.h index 5ce4e229ae..75611008ce 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -1087,11 +1087,14 @@ class Elem : public ReferenceCountedObject, * .) 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; diff --git a/src/geom/cell_polyhedron.C b/src/geom/cell_polyhedron.C index decfea6f2d..9f1bde025f 100644 --- a/src/geom/cell_polyhedron.C +++ b/src/geom/cell_polyhedron.C @@ -494,6 +494,19 @@ std::pair 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; diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C index 5975eb3045..2353fb72d6 100644 --- a/src/geom/elem_quality.C +++ b/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 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 diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C index d9458c2ba2..41ced6f6b7 100644 --- a/tests/geom/elem_test.C +++ b/tests/geom/elem_test.C @@ -7,6 +7,7 @@ #include #include #include +#include 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 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 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 diff --git a/tests/geom/elem_test.h b/tests/geom/elem_test.h index e72d7c5745..e0f227ebfd 100644 --- a/tests/geom/elem_test.h +++ b/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