From e3e378429ab4088cfabfa5b9a4d9f0eb7146cd2b Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Thu, 21 May 2026 10:58:16 -0600 Subject: [PATCH 1/4] Implement all_tri support for C0Polygon elements Splits each C0Polygon into the TRI3s defined by its existing triangulation. The per-element subelement container is changed from a fixed-size array to a vector since polygons can produce more than the 6 subelements assumed for 2D/3D primitives, and max_subelems is now computed by also scanning polygons in 2D meshes. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/mesh/mesh_modification.C | 43 +++++++++++++++++++++++++++++++++--- 1 file changed, 40 insertions(+), 3 deletions(-) diff --git a/src/mesh/mesh_modification.C b/src/mesh/mesh_modification.C index 2512ac914ab..674efc1dbc5 100644 --- a/src/mesh/mesh_modification.C +++ b/src/mesh/mesh_modification.C @@ -31,6 +31,8 @@ #include "libmesh/cell_tet4.h" #include "libmesh/cell_tet10.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 +451,15 @@ 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 can be split into an arbitrary number of triangles + // depending on their number of sides, so we have to scan the mesh + // to find the largest split we will need. + if (mesh.mesh_dimension() == 2) + for (const Elem * elem : mesh.element_ptr_range()) + if (const Polygon * poly = dynamic_cast(elem)) + max_subelems = std::max(max_subelems, poly->n_subtriangles()); + 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 +506,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, 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> subelem(max_subelems); switch (etype) { @@ -1246,6 +1258,31 @@ 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(elem); + const unsigned int n_subtri = polygon->n_subtriangles(); + for (unsigned int t = 0; t != n_subtri; ++t) + { + const std::array 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; + } + // No need to split elements that are already simplicial: case EDGE2: case EDGE3: @@ -1382,7 +1419,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); From e2b8d41dd2a0a4934f13580076643318bee6b485 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Thu, 21 May 2026 11:27:37 -0600 Subject: [PATCH 2/4] Add unit test for polygon conversion to tet Co-Authored-By: Claude Opus 4.7 (1M context) --- tests/mesh/all_tri.C | 84 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/tests/mesh/all_tri.C b/tests/mesh/all_tri.C index 18d10d8a86a..2a8c789433e 100644 --- a/tests/mesh/all_tri.C +++ b/tests/mesh/all_tri.C @@ -1,6 +1,8 @@ #include #include #include +#include +#include #include #include #include @@ -8,6 +10,8 @@ #include "test_comm.h" #include "libmesh_cppunit.h" +#include + using namespace libMesh; @@ -28,6 +32,8 @@ public: CPPUNIT_TEST( testAllTriQuad ); CPPUNIT_TEST( testAllTriQuad8 ); CPPUNIT_TEST( testAllTriQuad9 ); + CPPUNIT_TEST( testAllTriC0Polygon ); + CPPUNIT_TEST( testAllTriC0PolygonOctagon ); #endif // 3D tests @@ -112,6 +118,84 @@ 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(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 octagon = std::make_unique(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()); + } }; From 75a5427e2c483b730a2375aab49deb533e8733bf Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Thu, 21 May 2026 12:54:11 -0600 Subject: [PATCH 3/4] Implement all_tri support for C0Polyhedron elements Splits each C0Polyhedron into the TET4s defined by its existing tetrahedralization. The max_subelems scan is extended to also walk polyhedra in 3D so the per-element subelement vector and global ID offsets accommodate polyhedra that produce more than the 6 tets a HEX8 would. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/mesh/mesh_modification.C | 42 +++++++++++++++++++++++++++++++----- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/src/mesh/mesh_modification.C b/src/mesh/mesh_modification.C index 674efc1dbc5..09bae6a8da7 100644 --- a/src/mesh/mesh_modification.C +++ b/src/mesh/mesh_modification.C @@ -30,6 +30,8 @@ #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" @@ -451,13 +453,16 @@ 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 can be split into an arbitrary number of triangles - // depending on their number of sides, so we have to scan the mesh - // to find the largest split we will need. - if (mesh.mesh_dimension() == 2) - for (const Elem * elem : mesh.element_ptr_range()) + // 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(elem)) max_subelems = std::max(max_subelems, poly->n_subtriangles()); + else if (const Polyhedron * polyhedron = dynamic_cast(elem)) + max_subelems = std::max(max_subelems, polyhedron->n_subelements()); + } mesh.comm().max(max_subelems); new_elements.reserve (max_subelems*n_orig_elem); @@ -1283,6 +1288,33 @@ void MeshTools::Modification::all_tri (MeshBase & mesh) 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(elem); + const unsigned int n_sub = polyhedron->n_subelements(); + for (unsigned int t = 0; t != n_sub; ++t) + { + const std::array 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: From 1adba02a26d72ed7289b70b122032eae69eced4d Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Thu, 21 May 2026 13:39:35 -0600 Subject: [PATCH 4/4] Add unit tests for C0Polyhedron all_tri conversion Covers both tetrahedralization paths: a cube (no mid-element node) and a hexagonal prism (mid-element node added to make the split work). Co-Authored-By: Claude Opus 4.7 (1M context) --- tests/mesh/all_tri.C | 110 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/tests/mesh/all_tri.C b/tests/mesh/all_tri.C index 2a8c789433e..5249693114b 100644 --- a/tests/mesh/all_tri.C +++ b/tests/mesh/all_tri.C @@ -1,6 +1,8 @@ #include #include #include +#include +#include #include #include #include @@ -42,6 +44,8 @@ public: CPPUNIT_TEST( testAllTriPrism18 ); CPPUNIT_TEST( testAllTriPrism20 ); CPPUNIT_TEST( testAllTriPrism21 ); + CPPUNIT_TEST( testAllTriC0PolyhedronCube ); + CPPUNIT_TEST( testAllTriC0PolyhedronHexagonalPrism ); #endif CPPUNIT_TEST_SUITE_END(); @@ -196,6 +200,112 @@ public: 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 & points, + const std::vector> & nodes_on_side) + { + ReplicatedMesh mesh(*TestCommWorld, /*dim=*/3); + + for (auto p : index_range(points)) + mesh.add_point(points[p], /*id=*/p); + + std::vector> 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(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 mid_elem_node; + std::unique_ptr polyhedron = + std::make_unique(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(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 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> 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 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> 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); + } };