All tri for polygonal and polyhedra elements

src/mesh/mesh_modification.C

@@ -30,7 +30,11 @@
 #include "libmesh/function_base.h"
 #include "libmesh/cell_tet4.h"
 #include "libmesh/cell_tet10.h"
+#include "libmesh/cell_c0polyhedron.h"
+#include "libmesh/cell_polyhedron.h"
 #include "libmesh/elem_range.h"
+#include "libmesh/face_c0polygon.h"
+#include "libmesh/face_polygon.h"
 #include "libmesh/face_tri3.h"
 #include "libmesh/face_tri6.h"
 #include "libmesh/libmesh_logging.h"
@@ -449,6 +453,18 @@ void MeshTools::Modification::all_tri (MeshBase & mesh)
   if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
     max_subelems = 6;
 
+  // 2D polygons and 3D polyhedra can be split into an arbitrary
+  // number of triangles/tetrahedra depending on their topology, so we
+  // have to scan the mesh to find the largest split we will need.
+  for (const Elem * elem : mesh.element_ptr_range())
+    {
+      if (const Polygon * poly = dynamic_cast<const Polygon *>(elem))
+        max_subelems = std::max(max_subelems, poly->n_subtriangles());
+      else if (const Polyhedron * polyhedron = dynamic_cast<const Polyhedron *>(elem))
+        max_subelems = std::max(max_subelems, polyhedron->n_subelements());
+    }
+  mesh.comm().max(max_subelems);
+
   new_elements.reserve (max_subelems*n_orig_elem);
 
   // If the original mesh has *side* boundary data, we carry that over
@@ -495,8 +511,9 @@ void MeshTools::Modification::all_tri (MeshBase & mesh)
           libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
 
         // The new elements we will split the quad into.
-        // In 3D we may need as many as 6 tets per hex
-        std::array<std::unique_ptr<Elem>, 6> subelem {};
+        // In 3D we may need as many as 6 tets per hex, and in 2D
+        // we may need as many subtriangles as a polygon has.
+        std::vector<std::unique_ptr<Elem>> subelem(max_subelems);
 
         switch (etype)
           {
@@ -1246,6 +1263,58 @@ void MeshTools::Modification::all_tri (MeshBase & mesh)
               break;
             }
 
+          case C0POLYGON:
+            {
+              // Split a C0Polygon into the triangles defined by its
+              // current triangulation.  This relies on the user having
+              // a valid triangulation (the constructor sets a default
+              // one, and the user can refresh it via retriangulate()
+              // after moving nodes).
+              const C0Polygon * polygon = cast_ptr<const C0Polygon *>(elem);
+              const unsigned int n_subtri = polygon->n_subtriangles();
+              for (unsigned int t = 0; t != n_subtri; ++t)
+                {
+                  const std::array<int, 3> tri = polygon->subtriangle(t);
+                  if (tri[0] < 0 || tri[1] < 0 || tri[2] < 0)
+                    libmesh_not_implemented_msg
+                      ("Cannot convert a C0Polygon whose triangulation\n"
+                       "introduces special (non-vertex) points");
+                  subelem[t] = Elem::build(TRI3);
+                  subelem[t]->set_node(0, elem->node_ptr(tri[0]));
+                  subelem[t]->set_node(1, elem->node_ptr(tri[1]));
+                  subelem[t]->set_node(2, elem->node_ptr(tri[2]));
+                }
+
+              break;
+            }
+
+          case C0POLYHEDRON:
+            {
+              // Split a C0Polyhedron into the tetrahedra defined by its
+              // current tetrahedralization.  If the polyhedron required
+              // a mid-element node, the user is expected to have added
+              // that node to the mesh during construction; we just
+              // reference it via the polyhedron's node pointers.
+              const C0Polyhedron * polyhedron =
+                cast_ptr<const C0Polyhedron *>(elem);
+              const unsigned int n_sub = polyhedron->n_subelements();
+              for (unsigned int t = 0; t != n_sub; ++t)
+                {
+                  const std::array<int, 4> tet = polyhedron->subelement(t);
+                  if (tet[0] < 0 || tet[1] < 0 || tet[2] < 0 || tet[3] < 0)
+                    libmesh_not_implemented_msg
+                      ("Cannot convert a C0Polyhedron whose triangulation\n"
+                       "introduces special (non-vertex) points");
+                  subelem[t] = Elem::build(TET4);
+                  subelem[t]->set_node(0, elem->node_ptr(tet[0]));
+                  subelem[t]->set_node(1, elem->node_ptr(tet[1]));
+                  subelem[t]->set_node(2, elem->node_ptr(tet[2]));
+                  subelem[t]->set_node(3, elem->node_ptr(tet[3]));
+                }
+
+              break;
+            }
+
             // No need to split elements that are already simplicial:
           case EDGE2:
           case EDGE3:
@@ -1382,7 +1451,7 @@ void MeshTools::Modification::all_tri (MeshBase & mesh)
               // the same on all processors, therefore keeping the Mesh
               // in sync.  Note: we offset the new IDs by the max of the
               // pre-existing ids to avoid conflicting with originals.
-              subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
+              subelem[i]->set_id( max_orig_id + max_subelems*elem->id() + i );
 
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
               subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);

tests/mesh/all_tri.C

@@ -1,13 +1,19 @@
 #include <libmesh/libmesh.h>
 #include <libmesh/replicated_mesh.h>
 #include <libmesh/elem.h>
+#include <libmesh/cell_c0polyhedron.h>
+#include <libmesh/cell_polyhedron.h>
+#include <libmesh/face_c0polygon.h>
+#include <libmesh/face_polygon.h>
 #include <libmesh/mesh_generation.h>
 #include <libmesh/mesh_modification.h>
 #include <libmesh/boundary_info.h>
 
 #include "test_comm.h"
 #include "libmesh_cppunit.h"
 
+#include <cmath>
+
 
 using namespace libMesh;
 
@@ -28,6 +34,8 @@ public:
   CPPUNIT_TEST( testAllTriQuad );
   CPPUNIT_TEST( testAllTriQuad8 );
   CPPUNIT_TEST( testAllTriQuad9 );
+  CPPUNIT_TEST( testAllTriC0Polygon );
+  CPPUNIT_TEST( testAllTriC0PolygonOctagon );
 #endif
 
   // 3D tests
@@ -36,6 +44,8 @@ public:
   CPPUNIT_TEST( testAllTriPrism18 );
   CPPUNIT_TEST( testAllTriPrism20 );
   CPPUNIT_TEST( testAllTriPrism21 );
+  CPPUNIT_TEST( testAllTriC0PolyhedronCube );
+  CPPUNIT_TEST( testAllTriC0PolyhedronHexagonalPrism );
 #endif
 
   CPPUNIT_TEST_SUITE_END();
@@ -112,6 +122,190 @@ public:
   void testAllTriPrism18() { LOG_UNIT_TEST; test_helper_3D(PRISM18, /*nelem=*/6, /*nbcs=*/12); }
   void testAllTriPrism20() { LOG_UNIT_TEST; test_helper_3D(PRISM20, /*nelem=*/6, /*nbcs=*/12); }
   void testAllTriPrism21() { LOG_UNIT_TEST; test_helper_3D(PRISM21, /*nelem=*/6, /*nbcs=*/12); }
+
+  // Build a C0Polygon paving (triangles, quads, hexagons) via
+  // build_square and split it into a pure TRI3 mesh.
+  void testAllTriC0Polygon()
+  {
+    LOG_UNIT_TEST;
+
+    ReplicatedMesh mesh(*TestCommWorld, /*dim=*/2);
+
+    MeshTools::Generation::build_square(mesh,
+                                        /*nx=*/2, /*ny=*/2,
+                                        /*xmin=*/0., /*xmax=*/1.,
+                                        /*ymin=*/0., /*ymax=*/1.,
+                                        C0POLYGON);
+
+    // The post-all_tri element count is the sum of the per-polygon
+    // subtriangle counts, and external sides should be preserved one
+    // for one as TRI3 sides.
+    dof_id_type n_elem_expected = 0;
+    for (const Elem * elem : mesh.element_ptr_range())
+      {
+        const Polygon * poly = dynamic_cast<const Polygon *>(elem);
+        CPPUNIT_ASSERT(poly != nullptr);
+        n_elem_expected += poly->n_subtriangles();
+      }
+
+    const std::size_t n_bcs_before =
+      mesh.get_boundary_info().n_boundary_conds();
+
+    MeshTools::Modification::all_tri(mesh);
+
+    CPPUNIT_ASSERT_EQUAL(n_elem_expected, mesh.n_elem());
+    CPPUNIT_ASSERT_EQUAL(n_bcs_before,
+                         mesh.get_boundary_info().n_boundary_conds());
+
+    for (const Elem * elem : mesh.element_ptr_range())
+      CPPUNIT_ASSERT_EQUAL(ElemType(TRI3), elem->type());
+  }
+
+  // A single regular octagon (8 sides, 6 subtriangles) exercises the
+  // path where a polygon requires more subelements than any non-polygon
+  // 2D element type would.
+  void testAllTriC0PolygonOctagon()
+  {
+    LOG_UNIT_TEST;
+
+    constexpr unsigned int n_sides = 8;
+
+    ReplicatedMesh mesh(*TestCommWorld, /*dim=*/2);
+
+    std::unique_ptr<Elem> octagon = std::make_unique<C0Polygon>(n_sides);
+    for (unsigned int i = 0; i < n_sides; ++i)
+      {
+        const Real angle = 2 * libMesh::pi * i / n_sides;
+        Node * node = mesh.add_point(Point(std::cos(angle), std::sin(angle), 0.),
+                                     /*id=*/i);
+        octagon->set_node(i, node);
+      }
+    octagon->set_id() = 0;
+    Elem * elem = mesh.add_elem(std::move(octagon));
+
+    // Mark every external side with a boundary id so we can verify
+    // boundary information is transferred to the new triangles.
+    for (unsigned int s = 0; s < n_sides; ++s)
+      mesh.get_boundary_info().add_side(elem, s, /*bnd_id=*/0);
+
+    mesh.prepare_for_use();
+
+    MeshTools::Modification::all_tri(mesh);
+
+    // n_sides - 2 = 6 subtriangles
+    CPPUNIT_ASSERT_EQUAL(dof_id_type(n_sides - 2), mesh.n_elem());
+    CPPUNIT_ASSERT_EQUAL(std::size_t(n_sides),
+                         mesh.get_boundary_info().n_boundary_conds());
+
+    for (const Elem * e : mesh.element_ptr_range())
+      CPPUNIT_ASSERT_EQUAL(ElemType(TRI3), e->type());
+  }
+
+protected:
+
+  // Builds a C0Polyhedron in mesh whose faces are described by
+  // nodes_on_side (indices into the existing mesh node list), runs
+  // all_tri, and verifies that the result is a pure TET4 mesh with the
+  // expected sub-element count and preserved boundary data.
+  void test_helper_c0polyhedron
+    (const std::vector<Point> & points,
+     const std::vector<std::vector<unsigned int>> & nodes_on_side)
+  {
+    ReplicatedMesh mesh(*TestCommWorld, /*dim=*/3);
+
+    for (auto p : index_range(points))
+      mesh.add_point(points[p], /*id=*/p);
+
+    std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
+    for (auto s : index_range(nodes_on_side))
+      {
+        const auto & nodes_on_s = nodes_on_side[s];
+        sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
+        for (auto i : index_range(nodes_on_s))
+          sides[s]->set_node(i, mesh.node_ptr(nodes_on_s[i]));
+      }
+
+    std::unique_ptr<Node> mid_elem_node;
+    std::unique_ptr<Elem> polyhedron =
+      std::make_unique<C0Polyhedron>(sides, mid_elem_node);
+    if (mid_elem_node)
+      mesh.add_node(std::move(mid_elem_node));
+    polyhedron->set_id() = 0;
+    Elem * elem = mesh.add_elem(std::move(polyhedron));
+
+    const auto * poly = cast_ptr<const C0Polyhedron *>(elem);
+    const dof_id_type n_elem_expected = poly->n_subelements();
+
+    // Mark every external face with a boundary id so we can verify
+    // boundary information is transferred to the new tets.
+    for (unsigned int s = 0; s < elem->n_sides(); ++s)
+      mesh.get_boundary_info().add_side(elem, s, /*bnd_id=*/0);
+
+    // The number of boundary triangles produced is the total number of
+    // subtriangles across the polyhedron's polygonal faces.
+    std::size_t n_bcs_expected = 0;
+    for (unsigned int s = 0; s < elem->n_sides(); ++s)
+      n_bcs_expected += sides[s]->n_subtriangles();
+
+    mesh.prepare_for_use();
+
+    MeshTools::Modification::all_tri(mesh);
+
+    CPPUNIT_ASSERT_EQUAL(n_elem_expected, mesh.n_elem());
+    CPPUNIT_ASSERT_EQUAL(n_bcs_expected,
+                         mesh.get_boundary_info().n_boundary_conds());
+
+    for (const Elem * e : mesh.element_ptr_range())
+      CPPUNIT_ASSERT_EQUAL(ElemType(TET4), e->type());
+  }
+
+public:
+
+  // A cube built as a C0Polyhedron exercises the path where the
+  // optimal tetrahedralization succeeds (no mid-element node needed).
+  void testAllTriC0PolyhedronCube()
+  {
+    LOG_UNIT_TEST;
+
+    const std::vector<Point> points =
+      { {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
+        {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1} };
+
+    const std::vector<std::vector<unsigned int>> nodes_on_side =
+      { {0, 1, 2, 3},   // min z
+        {0, 1, 5, 4},   // min y
+        {2, 6, 5, 1},   // max x
+        {2, 3, 7, 6},   // max y
+        {0, 4, 7, 3},   // min x
+        {5, 6, 7, 4} }; // max z
+
+    test_helper_c0polyhedron(points, nodes_on_side);
+  }
+
+  // A hexagonal prism exercises the fallback path where a mid-element
+  // node is added to tetrahedralize the polyhedron.
+  void testAllTriC0PolyhedronHexagonalPrism()
+  {
+    LOG_UNIT_TEST;
+
+    const std::vector<Point> points =
+      { { 0, -2, 0}, {-1, -1, 0}, {-1, 1, 0},
+        { 0,  2, 0}, { 1,  1, 0}, { 1, -1, 0},
+        { 0, -2, 1}, {-1, -1, 1}, {-1, 1, 1},
+        { 0,  2, 1}, { 1,  1, 1}, { 1, -1, 1} };
+
+    const std::vector<std::vector<unsigned int>> nodes_on_side =
+      { {0, 1, 2, 3, 4, 5},
+        {0, 1,  7,  6},
+        {1, 2,  8,  7},
+        {2, 3,  9,  8},
+        {3, 4, 10,  9},
+        {4, 5, 11, 10},
+        {5, 0,  6, 11},
+        {6, 7,  8,  9, 10, 11} };
+
+    test_helper_c0polyhedron(points, nodes_on_side);
+  }
 };