all_tri.C

Go to the documentation of this file.

#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/mesh_tools.h>
#include <libmesh/boundary_info.h>
 
#include "test_comm.h"
#include "libmesh_cppunit.h"
 
#include <cmath>
 
 
using namespace libMesh;
 
class AllTriTest : public CppUnit::TestCase
{
public:
  LIBMESH_CPPUNIT_TEST_SUITE( AllTriTest );
 
  // 2D tests
#if LIBMESH_DIM > 1
  CPPUNIT_TEST( testAllTriTri );
  CPPUNIT_TEST( testAllTriQuad );
  CPPUNIT_TEST( testAllTriQuad8 );
  CPPUNIT_TEST( testAllTriQuad9 );
  CPPUNIT_TEST( testAllTriC0Polygon );
  CPPUNIT_TEST( testAllTriC0PolygonOctagon );
#endif
 
  // 3D tests
#if LIBMESH_DIM > 2
  CPPUNIT_TEST( testAllTriPrism6 );
  CPPUNIT_TEST( testAllTriPrism18 );
  CPPUNIT_TEST( testAllTriPrism20 );
  CPPUNIT_TEST( testAllTriPrism21 );
  CPPUNIT_TEST( testAllTriPyramid5 );
  CPPUNIT_TEST( testAllTriPyramid14 );
  CPPUNIT_TEST( testAllTriC0PolyhedronCube );
  CPPUNIT_TEST( testAllTriC0PolyhedronHexagonalPrism );
#endif
 
  CPPUNIT_TEST_SUITE_END();
 
protected:
  // Helper function called by the test implementations, saves a few lines of code.
  void test_helper(ElemType elem_type,
                   dof_id_type n_elem_expected,
                   std::size_t n_boundary_conds_expected)
  {
    ReplicatedMesh mesh(*TestCommWorld);
 
    // Build a 2x1 2D or 1x1x1 3D mesh, ask to split it into simplices
    const unsigned int dim = Elem::type_to_dim_map[elem_type];
    MeshTools::Generation::build_cube(mesh,
                                      /*nx=*/dim < 3 ? 2 : 1,
                                      /*ny=*/1,
                                      /*nz=*/dim < 3 ? 0 : 1,
                                      /*xmin=*/0., /*xmax=*/1.,
                                      /*ymin=*/0., /*ymax=*/1.,
                                      /*zmin=*/0., /*zmax=*/1.,
                                      elem_type);
 
    MeshTools::Modification::all_tri(mesh);
 
    // Make sure that the expected number of elements is found.
    CPPUNIT_ASSERT_EQUAL(n_elem_expected, mesh.n_elem());
 
    const BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
    // Make sure the expected number of BCs is found.
    CPPUNIT_ASSERT_EQUAL(n_boundary_conds_expected, boundary_info.n_boundary_conds());
 
    // In these tests we expect all our elements to have the proper
    // orientation, and we set up the inputs such that our outputs
    // should be non-curved and should have no boundary sides that
    // aren't on the previous external boundaries
    for (const Elem * elem : mesh.element_ptr_range())
      {
        CPPUNIT_ASSERT(!elem->is_flipped());
        CPPUNIT_ASSERT(elem->has_affine_map());
        for (auto s : elem->side_index_range())
          if (!elem->neighbor_ptr(s))
            CPPUNIT_ASSERT_EQUAL(boundary_info.n_boundary_ids(elem, s), 1u);
      }
 
    // We set up inputs with measure 1
    LIBMESH_ASSERT_NUMBERS_EQUAL(1, MeshTools::volume(mesh), TOLERANCE);
  }
 
public:
  void setUp() {}
 
  void tearDown() {}
 
  // 4 TRIs no-op
  void testAllTriTri() { LOG_UNIT_TEST; test_helper(TRI3, /*nelem=*/4, /*nbcs=*/6); }
 
  // 2 quads split into 4 TRIs.
  void testAllTriQuad() { LOG_UNIT_TEST; test_helper(QUAD4, /*nelem=*/4, /*nbcs=*/6); }
 
  // 2 QUAD8s split into 4 TRIs.
  void testAllTriQuad8() { LOG_UNIT_TEST; test_helper(QUAD8, /*nelem=*/4, /*nbcs=*/6); }
 
  // 2 QUAD9s split into 4 TRIs.
  void testAllTriQuad9() { LOG_UNIT_TEST; test_helper(QUAD9, /*nelem=*/4, /*nbcs=*/6); }
 
  // 2 PRISMs split into 6 TETs with 2 boundary faces per side.
  void testAllTriPrism6() { LOG_UNIT_TEST; test_helper(PRISM6, /*nelem=*/6, /*nbcs=*/12); }
  void testAllTriPrism18() { LOG_UNIT_TEST; test_helper(PRISM18, /*nelem=*/6, /*nbcs=*/12); }
  void testAllTriPrism20() { LOG_UNIT_TEST; test_helper(PRISM20, /*nelem=*/6, /*nbcs=*/12); }
  void testAllTriPrism21() { LOG_UNIT_TEST; test_helper(PRISM21, /*nelem=*/6, /*nbcs=*/12); }
 
  // 6 PYRAMIDs split into 12 TETs with 2 boundary faces per side
  void testAllTriPyramid5() { LOG_UNIT_TEST; test_helper(PYRAMID5, /*nelem=*/12, /*nbcs=*/12); }
  void testAllTriPyramid14() { LOG_UNIT_TEST; test_helper(PYRAMID14, /*nelem=*/12, /*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;
    Real poly_volume = octagon->volume();
    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());
 
    LIBMESH_ASSERT_NUMBERS_EQUAL(poly_volume, MeshTools::volume(mesh), TOLERANCE);
  }
 
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=*/s);
 
    // 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();
 
    // Keep this map for future checks
    std::vector<std::array<int, 4>> sub_elem_sides_to_parent_side(elem->n_sub_elem());
    for (unsigned int b = 0; b < elem->n_sub_elem(); ++b)
      sub_elem_sides_to_parent_side[b] = poly->subelement_sides_to_poly_sides(b);
 
    Real poly_volume = poly->volume();
    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());
 
    // Check that the boundary info and numbering is as expected
    for (const Elem * e : mesh.element_ptr_range())
      for (const auto s : make_range(e->n_sides()))
      {
        // Get parent side
        const auto side = sub_elem_sides_to_parent_side[e->id()][s];
        // Check boundary exists on the tet side if the tet side is on a poly side
        if (side != invalid_int)
          CPPUNIT_ASSERT_EQUAL(mesh.get_boundary_info().has_boundary_id(e, s, side), true);
      }
    LIBMESH_ASSERT_NUMBERS_EQUAL(poly_volume, MeshTools::volume(mesh), TOLERANCE);
  }
 
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);
  }
};
 
 
CPPUNIT_TEST_SUITE_REGISTRATION( AllTriTest );