diff --git a/src/mesh/mesh_modification.C b/src/mesh/mesh_modification.C index 2512ac914ab..09bae6a8da7 100644 --- a/src/mesh/mesh_modification.C +++ b/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(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); // 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, 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 +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(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; + } + + 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: @@ -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); diff --git a/tests/mesh/all_tri.C b/tests/mesh/all_tri.C index 18d10d8a86a..5249693114b 100644 --- a/tests/mesh/all_tri.C +++ b/tests/mesh/all_tri.C @@ -1,6 +1,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -8,6 +12,8 @@ #include "test_comm.h" #include "libmesh_cppunit.h" +#include + 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(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()); + } + +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); + } };