VolumeTest Class Reference

Inheritance diagram for VolumeTest:

Public Member Functions
	LIBMESH_CPPUNIT_TEST_SUITE (VolumeTest)
	CPPUNIT_TEST (testTwistedVolume)
	CPPUNIT_TEST (testEdge3Volume)
	CPPUNIT_TEST (testEdge3Invertible)
	CPPUNIT_TEST (testEdge4Invertible)
	CPPUNIT_TEST (testQuad4Invertible)
	CPPUNIT_TEST (testTri3TrueCentroid)
	CPPUNIT_TEST (testQuad4TrueCentroid)
	CPPUNIT_TEST (testPyramid5TrueCentroid)
	CPPUNIT_TEST (testHex8TrueCentroid)
	CPPUNIT_TEST (testPrism6TrueCentroid)
	CPPUNIT_TEST (testHex20PLevelTrueCentroid)
	CPPUNIT_TEST (testQuad4AspectRatio)
	CPPUNIT_TEST (testQuad4Warpage)
	CPPUNIT_TEST (testQuad4MinMaxAngle)
	CPPUNIT_TEST (testQuad4Jacobian)
	CPPUNIT_TEST (testTri3AspectRatio)
	CPPUNIT_TEST (testTet4DihedralAngle)
	CPPUNIT_TEST (testTet4Jacobian)
	CPPUNIT_TEST (testC0PolygonSquare)
	CPPUNIT_TEST (testC0PolygonQuad)
	CPPUNIT_TEST (testC0PolygonPentagon)
	CPPUNIT_TEST (testC0PolygonHexagon)
	CPPUNIT_TEST (testC0PolyhedronCube)
	CPPUNIT_TEST (testC0PolyhedronHexagonalPrism)
	CPPUNIT_TEST_SUITE_END ()
void	setUp ()
void	tearDown ()
void	testTri3TrueCentroid ()
void	testQuad4TrueCentroid ()
void	testPyramid5TrueCentroid ()
void	testHex8TrueCentroid ()
void	testPrism6TrueCentroid ()
void	testHex20PLevelTrueCentroid ()
void	testTwistedVolume ()
void	testEdge3Volume ()
void	testEdge3Invertible ()
void	testEdge4Invertible ()
void	testQuad4Invertible ()
void	testQuad4AspectRatio ()
void	testTri3AspectRatio ()
void	testQuad4Warpage ()
void	testQuad4MinMaxAngle ()
void	testQuad4Jacobian ()
void	testTet4DihedralAngle ()
void	testTet4Jacobian ()
void	testC0Polygon (const std::vector< Point > &points, Real expected_volume)
void	testC0PolygonSquare ()
void	testC0PolygonQuad ()
void	testC0PolygonPentagon ()
void	testC0PolygonHexagon ()

Protected Member Functions
std::pair< std::unique_ptr< Elem >, std::vector< std::unique_ptr< Node > > >	construct_elem (const std::vector< Point > &pts, ElemType elem_type)
Elem *	buildC0Polyhedron (std::vector< std::shared_ptr< Polygon > > sides, MeshBase &mesh)
void	testC0PolygonMethods (MeshBase &mesh, unsigned int n_sides)
void	testC0PolyhedronMethods (MeshBase &mesh, bool has_midnode)
void	testElemVolume (const Elem *elem, Real expected_volume)
void	testC0PolyhedronCube ()
void	testC0PolyhedronHexagonalPrism ()
bool	test_elem_invertible (const std::vector< Point > &pts, ElemType elem_type)
void	test_true_centroid_and_volume (ElemType elem_type)

Detailed Description

Definition at line 25 of file volume_test.C.

Member Function Documentation

buildC0Polyhedron()

Elem * VolumeTest::buildC0Polyhedron ( std::vector< std::shared_ptr< Polygon > >  sides, MeshBase &  mesh  )

inlineprotected

Definition at line 997 of file volume_test.C.

  {
    std::unique_ptr<libMesh::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));
    libmesh_assert(dynamic_cast<C0Polyhedron *>(elem));
    mesh.prepare_for_use();
 
    return elem;
  }

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_node(), libMesh::libmesh_assert(), mesh, libMesh::MeshBase::prepare_for_use(), and libMesh::DofObject::set_id().

Referenced by testC0PolyhedronCube(), and testC0PolyhedronHexagonalPrism().

construct_elem()

std::pair< std::unique_ptr< Elem >, std::vector< std::unique_ptr< Node > > > VolumeTest::construct_elem ( const std::vector< Point > &  pts, ElemType  elem_type  )

inlineprotected

Definition at line 968 of file volume_test.C.

  {
    const unsigned int n_points = pts.size();
 
    // Create Nodes
    std::vector<std::unique_ptr<Node>> nodes(n_points);
    for (unsigned int i=0; i<n_points; i++)
      nodes[i] = Node::build(pts[i], /*id*/ i);
 
    // Create Elem, assign nodes
    std::unique_ptr<Elem> elem = Elem::build(elem_type, /*parent*/ nullptr);
 
    // Make sure we were passed consistent input to build this type of Elem
    libmesh_error_msg_if(elem->n_nodes() != n_points,
                         "Wrong number of points "
                         << n_points
                         << " provided to build a "
                         << Utility::enum_to_string(elem_type));
 
    for (unsigned int i=0; i<n_points; i++)
      elem->set_node(i, nodes[i].get());
 
    // Return Elem and Nodes we created
    return std::make_pair(std::move(elem), std::move(nodes));
  }

References libMesh::Elem::build(), libMesh::Node::build(), and libMesh::Utility::enum_to_string().

Referenced by test_elem_invertible(), testQuad4AspectRatio(), testQuad4Jacobian(), testQuad4MinMaxAngle(), testQuad4Warpage(), testTet4DihedralAngle(), testTet4Jacobian(), and testTri3AspectRatio().

CPPUNIT_TEST() [1/24]

VolumeTest::CPPUNIT_TEST

(

testC0PolygonHexagon

)

CPPUNIT_TEST() [2/24]

VolumeTest::CPPUNIT_TEST

(

testC0PolygonPentagon

)

CPPUNIT_TEST() [3/24]

VolumeTest::CPPUNIT_TEST

(

testC0PolygonQuad

)

CPPUNIT_TEST() [4/24]

VolumeTest::CPPUNIT_TEST

(

testC0PolygonSquare

)

CPPUNIT_TEST() [5/24]

VolumeTest::CPPUNIT_TEST

(

testC0PolyhedronCube

)

CPPUNIT_TEST() [6/24]

VolumeTest::CPPUNIT_TEST

(

testC0PolyhedronHexagonalPrism

)

CPPUNIT_TEST() [7/24]

VolumeTest::CPPUNIT_TEST

(

testEdge3Invertible

)

CPPUNIT_TEST() [8/24]

VolumeTest::CPPUNIT_TEST

(

testEdge3Volume

)

CPPUNIT_TEST() [9/24]

VolumeTest::CPPUNIT_TEST

(

testEdge4Invertible

)

CPPUNIT_TEST() [10/24]

VolumeTest::CPPUNIT_TEST

(

testHex20PLevelTrueCentroid

)

CPPUNIT_TEST() [11/24]

VolumeTest::CPPUNIT_TEST

(

testHex8TrueCentroid

)

CPPUNIT_TEST() [12/24]

VolumeTest::CPPUNIT_TEST

(

testPrism6TrueCentroid

)

CPPUNIT_TEST() [13/24]

VolumeTest::CPPUNIT_TEST

(

testPyramid5TrueCentroid

)

CPPUNIT_TEST() [14/24]

VolumeTest::CPPUNIT_TEST

(

testQuad4AspectRatio

)

CPPUNIT_TEST() [15/24]

VolumeTest::CPPUNIT_TEST

(

testQuad4Invertible

)

CPPUNIT_TEST() [16/24]

VolumeTest::CPPUNIT_TEST

(

testQuad4Jacobian

)

CPPUNIT_TEST() [17/24]

VolumeTest::CPPUNIT_TEST

(

testQuad4MinMaxAngle

)

CPPUNIT_TEST() [18/24]

VolumeTest::CPPUNIT_TEST

(

testQuad4TrueCentroid

)

CPPUNIT_TEST() [19/24]

VolumeTest::CPPUNIT_TEST

(

testQuad4Warpage

)

CPPUNIT_TEST() [20/24]

VolumeTest::CPPUNIT_TEST

(

testTet4DihedralAngle

)

CPPUNIT_TEST() [21/24]

VolumeTest::CPPUNIT_TEST

(

testTet4Jacobian

)

CPPUNIT_TEST() [22/24]

VolumeTest::CPPUNIT_TEST

(

testTri3AspectRatio

)

CPPUNIT_TEST() [23/24]

VolumeTest::CPPUNIT_TEST

(

testTri3TrueCentroid

)

CPPUNIT_TEST() [24/24]

VolumeTest::CPPUNIT_TEST

(

testTwistedVolume

)

CPPUNIT_TEST_SUITE_END()

VolumeTest::CPPUNIT_TEST_SUITE_END

(

)

LIBMESH_CPPUNIT_TEST_SUITE()

VolumeTest::LIBMESH_CPPUNIT_TEST_SUITE

(

VolumeTest

)

setUp()

void VolumeTest::setUp ( )

inline

Definition at line 57 of file volume_test.C.

  {
  }

tearDown()

void VolumeTest::tearDown ( )

inline

Definition at line 61 of file volume_test.C.

  {
  }

test_elem_invertible()

bool VolumeTest::test_elem_invertible ( const std::vector< Point > &  pts, ElemType  elem_type  )

inlineprotected

Definition at line 1271 of file volume_test.C.

  {
    // Construct Elem of desired type
    auto [elem, nodes] = this->construct_elem(pts, elem_type);
    libmesh_ignore(nodes);
 
    // Return whether or not this Elem has an invertible map
    return elem->has_invertible_map();
  }

References construct_elem(), and libMesh::libmesh_ignore().

Referenced by testEdge3Invertible(), testEdge4Invertible(), and testQuad4Invertible().

test_true_centroid_and_volume()

void VolumeTest::test_true_centroid_and_volume ( ElemType  elem_type )

inlineprotected

Definition at line 1284 of file volume_test.C.

  {
    // Check that derived class true_centroid() gives same result as
    // the base class Elem::true_centroid() on a 2x2x2 mesh of 10%
    // distorted elements.
    {
      ReplicatedMesh mesh(*TestCommWorld);
 
      MeshTools::Generation::build_cube(mesh,
                                        /*nelem=*/2, /*nelem=*/2, /*nelem=*/2,
                                        /*xmin=*/-1, /*xmax=*/1,
                                        /*ymin=*/-1, /*ymax=*/1,
                                        /*zmin=*/-1, /*zmax=*/1,
                                        elem_type);
 
      MeshTools::Modification::distort(mesh,
                                       /*factor=*/0.1,
                                       /*perturb_boundary=*/false);
 
      for (const auto & elem : mesh.element_ptr_range())
        {
          Point derived_centroid = elem->true_centroid();
          Point base_centroid = elem->Elem::true_centroid();
 
          // Debugging: check results in detail
          // auto flags = libMesh::out.flags();
          // libMesh::out << std::scientific << std::setprecision(16);
          // libMesh::out << "derived_centroid = " << derived_centroid << std::endl;
          // libMesh::out << "base_centroid = " << base_centroid << std::endl;
          // libMesh::out.flags(flags);
 
          CPPUNIT_ASSERT(derived_centroid.absolute_fuzzy_equals(base_centroid, 50*TOLERANCE*TOLERANCE));
 
          // Make sure that base class and "optimized" routines for computing the cell volume agree
          Real derived_volume = elem->volume();
          Real base_volume = elem->Elem::volume();
          LIBMESH_ASSERT_FP_EQUAL(base_volume, derived_volume, 50*TOLERANCE*TOLERANCE);
        }
    }
  }

References libMesh::TypeVector< T >::absolute_fuzzy_equals(), libMesh::MeshTools::Generation::build_cube(), libMesh::MeshTools::Modification::distort(), mesh, libMesh::Real, TestCommWorld, and libMesh::TOLERANCE.

Referenced by testHex8TrueCentroid(), testPrism6TrueCentroid(), and testPyramid5TrueCentroid().

testC0Polygon()

void VolumeTest::testC0Polygon ( const std::vector< Point > &  points, Real  expected_volume  )

inline

Definition at line 913 of file volume_test.C.

  {
    Mesh mesh(*TestCommWorld);
 
    const auto np = points.size();
    std::unique_ptr<Elem> polygon = std::make_unique<C0Polygon>(np);
 
    for (auto p : make_range(np))
      polygon->set_node(p, mesh.add_point(points[p], p));
 
    polygon->set_id() = 0;
    Elem * elem = mesh.add_elem(std::move(polygon));
 
    const Real derived_volume = elem->volume();
    const Real base_volume = elem->Elem::volume();
    LIBMESH_ASSERT_FP_EQUAL(base_volume, derived_volume, TOLERANCE*TOLERANCE);
    LIBMESH_ASSERT_FP_EQUAL(derived_volume, expected_volume, TOLERANCE*TOLERANCE);
 
    this->testC0PolygonMethods(mesh, np);
  }

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::make_range(), mesh, libMesh::Real, testC0PolygonMethods(), TestCommWorld, libMesh::TOLERANCE, and libMesh::Elem::volume().

Referenced by testC0PolygonHexagon(), testC0PolygonPentagon(), testC0PolygonQuad(), and testC0PolygonSquare().

testC0PolygonHexagon()

void VolumeTest::testC0PolygonHexagon ( )

inline

Definition at line 955 of file volume_test.C.

  {
    LOG_UNIT_TEST;
    testC0Polygon({{0,0}, {1,0}, {1.5,0.5}, {1,1}, {0,1}, {-0.5, 0.5}}, 1.5);
  }

References testC0Polygon().

testC0PolygonMethods()

void VolumeTest::testC0PolygonMethods ( MeshBase &  mesh, unsigned int  n_sides  )

inlineprotected

Definition at line 1012 of file volume_test.C.

  {
    Elem * elem = mesh.query_elem_ptr(0);
    bool found_elem = elem;
    mesh.comm().max(found_elem);
    CPPUNIT_ASSERT(found_elem);
 
    if (!elem)
      return;
 
    CPPUNIT_ASSERT_EQUAL(elem->type(), C0POLYGON);
    CPPUNIT_ASSERT_EQUAL(elem->n_nodes(), n_sides);
    CPPUNIT_ASSERT_EQUAL(elem->n_sub_elem(), n_sides);
    CPPUNIT_ASSERT_EQUAL(elem->n_sides(), n_sides);
    CPPUNIT_ASSERT_EQUAL(elem->n_vertices(), n_sides);
    CPPUNIT_ASSERT_EQUAL(elem->n_edges(), n_sides);
 
    // Even number of sides
    if (!(n_sides%2))
      {
        const unsigned int nsover2 = n_sides/2;
        for (auto i : make_range(nsover2))
          {
            CPPUNIT_ASSERT_EQUAL(elem->opposite_side(i), i+nsover2);
            CPPUNIT_ASSERT_EQUAL(elem->opposite_side(i+nsover2), i);
          }
      }
 
    for (unsigned int i : make_range(n_sides))
      {
        CPPUNIT_ASSERT(elem->is_vertex(i));
        CPPUNIT_ASSERT(!elem->is_edge(i));
        CPPUNIT_ASSERT(!elem->is_face(i));
        CPPUNIT_ASSERT(elem->is_node_on_side(i,i));
        CPPUNIT_ASSERT(elem->is_node_on_edge(i,i));
        CPPUNIT_ASSERT(elem->is_node_on_side((i+1)%n_sides,i));
        CPPUNIT_ASSERT(elem->is_node_on_edge((i+1)%n_sides,i));
        std::vector<unsigned int> nodes = elem->nodes_on_side(i);
        CPPUNIT_ASSERT_EQUAL(nodes.size(), std::size_t(2));
        CPPUNIT_ASSERT(nodes[0] == i ||
                       nodes[1] == i);
        CPPUNIT_ASSERT(nodes[0] == (i+1)%n_sides ||
                       nodes[1] == (i+1)%n_sides);
        std::vector<unsigned int> edge_nodes = elem->nodes_on_edge(i);
        CPPUNIT_ASSERT(nodes == edge_nodes);
 
 
        CPPUNIT_ASSERT_EQUAL(elem->local_side_node(i,0), i);
        CPPUNIT_ASSERT_EQUAL(elem->local_side_node(i,1), (i+1)%n_sides);
        CPPUNIT_ASSERT_EQUAL(elem->local_edge_node(i,0), i);
        CPPUNIT_ASSERT_EQUAL(elem->local_edge_node(i,1), (i+1)%n_sides);
 
        auto edge = elem->side_ptr(i);
        CPPUNIT_ASSERT_EQUAL(edge->type(), EDGE2);
        CPPUNIT_ASSERT_EQUAL(edge->node_ptr(0), mesh.node_ptr(i));
        CPPUNIT_ASSERT_EQUAL(edge->node_ptr(1), mesh.node_ptr((i+1)%n_sides));
      }
 
    CPPUNIT_ASSERT(!elem->is_flipped());
    elem->flip(&mesh.get_boundary_info());
    CPPUNIT_ASSERT(elem->is_flipped());
    elem->flip(&mesh.get_boundary_info());
  }

References libMesh::C0POLYGON, libMesh::ParallelObject::comm(), libMesh::EDGE2, libMesh::Elem::flip(), libMesh::MeshBase::get_boundary_info(), libMesh::Elem::is_edge(), libMesh::Elem::is_face(), libMesh::Elem::is_flipped(), libMesh::Elem::is_node_on_edge(), libMesh::Elem::is_node_on_side(), libMesh::Elem::is_vertex(), libMesh::Elem::local_edge_node(), libMesh::Elem::local_side_node(), libMesh::make_range(), libMesh::Parallel::Communicator::max(), mesh, libMesh::Elem::n_edges(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::n_sub_elem(), libMesh::Elem::n_vertices(), libMesh::MeshBase::node_ptr(), libMesh::Elem::nodes_on_edge(), libMesh::Elem::nodes_on_side(), libMesh::Elem::opposite_side(), libMesh::MeshBase::query_elem_ptr(), libMesh::Elem::side_ptr(), and libMesh::Elem::type().

Referenced by testC0Polygon().

testC0PolygonPentagon()

void VolumeTest::testC0PolygonPentagon ( )

inline

Definition at line 949 of file volume_test.C.

  {
    LOG_UNIT_TEST;
    testC0Polygon({{0,0}, {1,0}, {1.5,0.5}, {1,1}, {0,1}}, 1.25);
  }

References testC0Polygon().

testC0PolygonQuad()

void VolumeTest::testC0PolygonQuad ( )

inline

Definition at line 943 of file volume_test.C.

  {
    LOG_UNIT_TEST;
    testC0Polygon({{0,0}, {1,0}, {1,2}, {-1,1}}, 2.5);
  }

References testC0Polygon().

testC0PolygonSquare()

void VolumeTest::testC0PolygonSquare ( )

inline

Definition at line 937 of file volume_test.C.

  {
    LOG_UNIT_TEST;
    testC0Polygon({{0,0}, {1,0}, {1,1}, {0,1}}, 1);
  }

References testC0Polygon().

testC0PolyhedronCube()

void VolumeTest::testC0PolyhedronCube ( )

inlineprotected

Definition at line 1161 of file volume_test.C.

  {
    ReplicatedMesh mesh(*TestCommWorld);
 
    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(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);
 
    // See notes in elem_test.h
    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
 
    // Build all the sides.
    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]));
      }
 
    const auto poly = dynamic_cast<C0Polyhedron *>(buildC0Polyhedron(sides, mesh));
    testElemVolume(poly, 1);
  #ifdef LIBMESH_ENABLE_EXCEPTIONS
    testC0PolyhedronMethods(mesh, /*midnode*/false);
 
    // Check routine for subtet side to poly side mapping
    // NOTE: this test is hard coded to the elements used right now (to be reworked)
    const auto subtet0_sides_to_poly_sides = poly->subelement_sides_to_poly_sides(0);
    CPPUNIT_ASSERT_EQUAL(libMesh::invalid_int, subtet0_sides_to_poly_sides[0]);
    CPPUNIT_ASSERT_EQUAL(1, subtet0_sides_to_poly_sides[1]);
    CPPUNIT_ASSERT_EQUAL(0, subtet0_sides_to_poly_sides[2]);
    CPPUNIT_ASSERT_EQUAL(2, subtet0_sides_to_poly_sides[3]);
  #else
    testC0PolyhedronMethods(mesh, /*midnode*/true);
  #endif
  }

References libMesh::MeshBase::add_point(), buildC0Polyhedron(), libMesh::index_range(), libMesh::invalid_int, mesh, libMesh::MeshBase::node_ptr(), testC0PolyhedronMethods(), TestCommWorld, and testElemVolume().

testC0PolyhedronHexagonalPrism()

void VolumeTest::testC0PolyhedronHexagonalPrism ( )

inlineprotected

Definition at line 1213 of file volume_test.C.

  {
    ReplicatedMesh mesh(*TestCommWorld);
 
    mesh.add_point(Point(0, -2, 0), 0);
    mesh.add_point(Point(-1, -1, 0), 1);
    mesh.add_point(Point(-1, 1, 0), 2);
    mesh.add_point(Point(0, 2, 0), 3);
    mesh.add_point(Point(1, 1, 0), 4);
    mesh.add_point(Point(1, -1, 0), 5);
    mesh.add_point(Point(0, -2, 1), 6);
    mesh.add_point(Point(-1, -1, 1), 7);
    mesh.add_point(Point(-1, 1, 1), 8);
    mesh.add_point(Point(0, 2, 1), 9);
    mesh.add_point(Point(1, 1, 1), 10);
    mesh.add_point(Point(1, -1, 1), 11);
 
    // See notes in elem_test.h
    // With this ordering, we'll need to use a mid-node to tetrahedralize it
    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} };
 
    // Build all the sides.
    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]));
      }
 
    const auto poly = dynamic_cast<C0Polyhedron *>(buildC0Polyhedron(sides, mesh));
    CPPUNIT_ASSERT_EQUAL((unsigned int)13, poly->n_nodes());  // we should have a mid node
    testElemVolume(poly, 6);
    testC0PolyhedronMethods(mesh, /*midnode*/true);
 
    // Check routine for subtet side to poly side mapping
    const auto subtet0_sides_to_poly_sides = poly->subelement_sides_to_poly_sides(0);
    CPPUNIT_ASSERT_EQUAL(0, subtet0_sides_to_poly_sides[0]);
    CPPUNIT_ASSERT_EQUAL(libMesh::invalid_int, subtet0_sides_to_poly_sides[1]);
    CPPUNIT_ASSERT_EQUAL(libMesh::invalid_int, subtet0_sides_to_poly_sides[2]);
    CPPUNIT_ASSERT_EQUAL(libMesh::invalid_int, subtet0_sides_to_poly_sides[3]);
  }

References libMesh::MeshBase::add_point(), buildC0Polyhedron(), libMesh::index_range(), libMesh::invalid_int, mesh, libMesh::MeshBase::node_ptr(), testC0PolyhedronMethods(), TestCommWorld, and testElemVolume().

testC0PolyhedronMethods()

void VolumeTest::testC0PolyhedronMethods ( MeshBase &  mesh, bool  has_midnode  )

inlineprotected

Definition at line 1079 of file volume_test.C.

  {
    Elem * elem = mesh.query_elem_ptr(0);
    bool found_elem = elem;
    mesh.comm().max(found_elem);
    CPPUNIT_ASSERT(found_elem);
 
    if (!elem)
      return;
 
    CPPUNIT_ASSERT_EQUAL(elem->type(), C0POLYHEDRON);
 
    const int V = elem->n_vertices();
    const int E = elem->n_edges();
    const int F = elem->n_faces();
 
    if (!has_midnode)
      CPPUNIT_ASSERT_EQUAL(int(elem->n_nodes()), V);
    else
      CPPUNIT_ASSERT_EQUAL(int(elem->n_nodes()), V + 1);
    CPPUNIT_ASSERT_EQUAL(int(elem->n_sides()), F);
 
    // Euler characteristic for polygons homeomorphic to spheres is 2
    CPPUNIT_ASSERT_EQUAL(V-E+F, 2);
 
    int FE = 0;
    int FV = 0;
 
    std::vector<int> sides_on_vertex(V, 0);
    for (auto s : make_range(F))
      {
        auto side = elem->build_side_ptr(s);
        CPPUNIT_ASSERT_EQUAL(side->type(), C0POLYGON);
        FE += side->n_edges();
        const int SV = side->n_vertices();
        FV += SV;
 
        auto nos = elem->nodes_on_side(s);
        CPPUNIT_ASSERT_EQUAL(nos.size(), std::size_t(SV));
 
        for (auto v : make_range(SV))
          {
            Node * sidenode = side->node_ptr(v);
            const unsigned int n = elem->get_node_index(sidenode);
            CPPUNIT_ASSERT(n != invalid_uint);
            ++sides_on_vertex[n];
 
            // We're just looking at first order so far
            CPPUNIT_ASSERT(side->is_vertex(v));
            CPPUNIT_ASSERT(!side->is_edge(v));
            CPPUNIT_ASSERT(elem->is_vertex(n));
            CPPUNIT_ASSERT(!elem->is_edge(n));
            CPPUNIT_ASSERT(!elem->is_face(n));
            CPPUNIT_ASSERT(elem->is_node_on_side(n, s));
            CPPUNIT_ASSERT(std::find(nos.begin(), nos.end(), n) != nos.end());
          }
      }
 
    // We hit every edge from 2 faces
    CPPUNIT_ASSERT_EQUAL(E*2, FE);
 
    // We hit every vertex from at least 3 faces
    CPPUNIT_ASSERT_LESSEQUAL(FV, V*3);  // V*3 <= FV
    for (auto sov : sides_on_vertex)
      CPPUNIT_ASSERT_LESSEQUAL(sov, 3); // 3 <= sov
 
    CPPUNIT_ASSERT(!elem->is_flipped());
  }

References libMesh::Elem::build_side_ptr(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::ParallelObject::comm(), libMesh::Elem::get_node_index(), libMesh::invalid_uint, libMesh::Elem::is_edge(), libMesh::Elem::is_face(), libMesh::Elem::is_flipped(), libMesh::Elem::is_node_on_side(), libMesh::Elem::is_vertex(), libMesh::make_range(), libMesh::Parallel::Communicator::max(), mesh, libMesh::Elem::n_edges(), libMesh::Elem::n_faces(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::n_vertices(), libMesh::Elem::nodes_on_side(), libMesh::MeshBase::query_elem_ptr(), and libMesh::Elem::type().

Referenced by testC0PolyhedronCube(), and testC0PolyhedronHexagonalPrism().

testEdge3Invertible()

void VolumeTest::testEdge3Invertible ( )

inline

Definition at line 245 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // 1.) This is the original test which started the investigation
    // of determining invertibility.  In this test, the actual
    // midpoint of nodes 0 and 1 is 0.5*(1.100328e2 + 1.176528e2) =
    // 113.8428, so we can see that the middle node is closer to the
    // left endpoint. In this case, it is too close and the element is
    // not invertible.
    bool invertible = test_elem_invertible({
      Point(-3.566160e1, -6.690970e-1, 1.100328e2),
      Point(-3.566160e1, -6.690970e-1, 1.176528e2),
      Point(-3.566160e1, -6.690970e-1, 1.115568e2)}, EDGE3);
    CPPUNIT_ASSERT(!invertible);
 
    // 2.) Just like case 1, but now node 2 is at the midpoint, so
    // this case is invertible.
    invertible = test_elem_invertible({
      Point(-3.566160e1, -6.690970e-1, 1.100328e2),
      Point(-3.566160e1, -6.690970e-1, 1.176528e2),
      Point(-3.566160e1, -6.690970e-1, 113.8428)}, EDGE3);
    CPPUNIT_ASSERT(invertible);
 
    // 3.) Non-collinear case where the mid-edge node is "above" and "way
    // past" the right endpoint. This case is not invertible
    invertible = test_elem_invertible({Point(0, 0, 0), Point(1, 0, 0), Point(3.5, 1.5, 0)}, EDGE3);
    CPPUNIT_ASSERT(!invertible);
  }

References libMesh::EDGE3, and test_elem_invertible().

testEdge3Volume()

void VolumeTest::testEdge3Volume ( )

inline

Definition at line 201 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    Mesh mesh(*TestCommWorld);
    MeshTools::Generation::build_line (mesh, /*nelem=*/1, 0., 1., EDGE3);
    CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(3), mesh.n_nodes());
 
    auto edge3 = mesh.query_elem_ptr(0);
    if (!edge3) // We may be on a distributed mesh
      return;
 
    // Check unperturbed, straight edge case
    LIBMESH_ASSERT_FP_EQUAL(1.0, edge3->volume(), TOLERANCE*TOLERANCE);
 
    // Get references to the individual Edge3 nodes
    auto & middle_node = edge3->node_ref(2);
    auto & right_node = edge3->node_ref(1);
 
    // Check middle node perturbed in +x direction case. This should
    // not change the volume because it's still a straight line
    // element.
    middle_node = Point(0.5 + 1.e-3, 0., 0.);
    LIBMESH_ASSERT_FP_EQUAL(1.0, edge3->volume(), TOLERANCE*TOLERANCE);
 
    // Check middle node perturbed in -x direction case. This should
    // not change the volume because it's still a straight line
    // element.
    middle_node = Point(0.5 - 1.e-3, 0., 0.);
    LIBMESH_ASSERT_FP_EQUAL(1.0, edge3->volume(), TOLERANCE*TOLERANCE);
 
    // Check volume of actual curved element against pre-computed value.
    middle_node = Point(0.5, 0.25, 0.);
    right_node = Point(1., 1., 0.);
    LIBMESH_ASSERT_FP_EQUAL(1.4789428575446, edge3->volume(), TOLERANCE);
 
    // Compare with volume computed by base class Elem::volume() call
    // which uses quadrature.  We don't expect this to have full
    // floating point accuracy.
    middle_node = Point(0.5, 0.1, 0.);
    right_node = Point(1., 0., 0.);
    LIBMESH_ASSERT_FP_EQUAL(edge3->Elem::volume(), edge3->volume(), 1e-4);
  }

References libMesh::MeshTools::Generation::build_line(), libMesh::EDGE3, mesh, libMesh::MeshBase::n_nodes(), libMesh::MeshBase::query_elem_ptr(), TestCommWorld, and libMesh::TOLERANCE.

testEdge4Invertible()

void VolumeTest::testEdge4Invertible ( )

inline

Definition at line 275 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Reference Elem should be invertible
    {
      const Elem & edge4 = ReferenceElem::get(EDGE4);
      CPPUNIT_ASSERT(edge4.has_invertible_map());
    }
 
    // If node 2 goes to the left past -5/9 = -.555, the element becomes non-invertible
    {
      // x2 > -5/9, the map is still invertible
      bool invertible =
        test_elem_invertible({Point(-1, 0, 0), Point(1, 0, 0), Point(-0.5, 0, 0), Point(Real(1)/3, 0, 0)},
                  EDGE4);
      CPPUNIT_ASSERT(invertible);
 
      // x2 < -5/9, it is too close to x0 now
      invertible =
        test_elem_invertible({Point(-1, 0, 0), Point(1, 0, 0), Point(-0.57, 0, 0), Point(Real(1)/3, 0, 0)},
                  EDGE4);
      CPPUNIT_ASSERT(!invertible);
    }
 
    // If node 2 goes to the right past 5/21 ~ 0.2381, the element becomes non-invertible
    {
      // x2 < 5/21, the map should still be invertible
      bool invertible =
        test_elem_invertible({Point(-1, 0, 0), Point(1, 0, 0), Point(Real(3)/21, 0, 0), Point(Real(1)/3, 0, 0)},
                  EDGE4);
      CPPUNIT_ASSERT(invertible);
 
      // x2 > 5/21, x2 is too close to x3 now
      invertible =
        test_elem_invertible({Point(-1, 0, 0), Point(1, 0, 0), Point(Real(6)/21, 0, 0), Point(Real(1)/3, 0, 0)},
                  EDGE4);
      CPPUNIT_ASSERT(!invertible);
    }
  }

References libMesh::EDGE4, libMesh::ReferenceElem::get(), libMesh::Elem::has_invertible_map(), libMesh::Real, and test_elem_invertible().

testElemVolume()

void VolumeTest::testElemVolume ( const Elem *  elem, Real  expected_volume  )

inlineprotected

Definition at line 1150 of file volume_test.C.

  {
    const Real derived_volume = elem->volume();
    const Real base_volume = elem->Elem::volume();
    LIBMESH_ASSERT_FP_EQUAL(base_volume, derived_volume, TOLERANCE*TOLERANCE);
    LIBMESH_ASSERT_FP_EQUAL(derived_volume, expected_volume, TOLERANCE*TOLERANCE);
  }

References libMesh::Real, libMesh::TOLERANCE, and libMesh::Elem::volume().

Referenced by testC0PolyhedronCube(), and testC0PolyhedronHexagonalPrism().

testHex20PLevelTrueCentroid()

void VolumeTest::testHex20PLevelTrueCentroid ( )

inline

Definition at line 150 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Test that Elem base class true_centroid() implementation works
    // for an elevated p_level HEX20
    {
#ifdef LIBMESH_ENABLE_AMR
      ReplicatedMesh mesh(*TestCommWorld);
      MeshTools::Generation::build_cube(mesh,
                                        /*nelem=*/1, /*nelem=*/1, /*nelem=*/1,
                                        /*xmin=*/-1, /*xmax=*/1,
                                        /*ymin=*/-1, /*ymax=*/1,
                                        /*zmin=*/-1, /*zmax=*/1,
                                        HEX20);
      Elem * hex20 = mesh.elem_ptr(0);
      hex20->set_p_level(1);
      Point true_centroid = hex20->true_centroid();
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(0), TOLERANCE*TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(1), TOLERANCE*TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(2), TOLERANCE*TOLERANCE);
#endif // LIBMESH_ENABLE_AMR
    }
  }

References libMesh::MeshTools::Generation::build_cube(), libMesh::MeshBase::elem_ptr(), libMesh::HEX20, mesh, libMesh::Elem::set_p_level(), TestCommWorld, libMesh::TOLERANCE, and libMesh::Elem::true_centroid().

testHex8TrueCentroid()

void VolumeTest::testHex8TrueCentroid ( )

inline

Definition at line 147 of file volume_test.C.

{ LOG_UNIT_TEST; test_true_centroid_and_volume(HEX8); }

References libMesh::HEX8, and test_true_centroid_and_volume().

testPrism6TrueCentroid()

void VolumeTest::testPrism6TrueCentroid ( )

inline

Definition at line 148 of file volume_test.C.

{ LOG_UNIT_TEST; test_true_centroid_and_volume(PRISM6); }

References libMesh::PRISM6, and test_true_centroid_and_volume().

testPyramid5TrueCentroid()

void VolumeTest::testPyramid5TrueCentroid ( )

inline

Definition at line 130 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Test Pyramid5::true_centroid() gives the correct result for a reference element
    {
      const Elem & pyr5 = ReferenceElem::get(PYRAMID5);
      Point true_centroid = pyr5.true_centroid();
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(0), TOLERANCE*TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(1), TOLERANCE*TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(0.25, true_centroid(2), TOLERANCE*TOLERANCE);
    }
 
    // Currently there is not an optimized Pyramid5::true_centroid() to compare against
    test_true_centroid_and_volume(PYRAMID5);
  }

References libMesh::ReferenceElem::get(), libMesh::PYRAMID5, test_true_centroid_and_volume(), libMesh::TOLERANCE, and libMesh::Elem::true_centroid().

testQuad4AspectRatio()

void VolumeTest::testQuad4AspectRatio ( )

inline

Definition at line 407 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Case 1: Test that rigid body rotations of a unit square
    // quadrilateral that have no effect on the quality of the
    // element.
    {
      // Construct unit square QUAD4
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0)};
      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
      libmesh_ignore(nodes);
 
      // 1a) Unit square aspect ratio should be == 1
      Real aspect_ratio = elem->quality(ASPECT_RATIO);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/aspect_ratio, TOLERANCE);
 
      // 1b) Rotate all points about x-axis by 90 degrees
      Real cost = std::cos(.5*libMesh::pi);
      Real sint = std::sin(.5*libMesh::pi);
      RealTensorValue Rx(1, 0, 0,
                         0, cost, sint,
                         0, sint, cost);
 
      for (auto & pt : pts)
        pt = Rx * pt;
 
      aspect_ratio = elem->quality(ASPECT_RATIO);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/aspect_ratio, TOLERANCE);
 
      // 1c) Rotate all points about z-axis by 90 degrees
      RealTensorValue Rz(cost, -sint, 0,
                         sint,  cost, 0,
                         0,        0, 1);
 
      for (auto & pt : pts)
        pt = Rz * pt;
 
      aspect_ratio = elem->quality(ASPECT_RATIO);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/aspect_ratio, TOLERANCE);
 
      // 1d) Rotate all points about y-axis by 270 degrees
      RealTensorValue Ry(cost,  0, sint,
                         0,     1, 0,
                         -sint, 0, cost);
 
      for (int cnt=0; cnt<3; ++cnt)
        for (auto & pt : pts)
          pt = Ry * pt;
 
      aspect_ratio = elem->quality(ASPECT_RATIO);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/aspect_ratio, TOLERANCE);
    }
 
    // Case 2: Rhombus QUAD4. This case should have an aspect ratio of
    // 1/sin(theta), where theta is the acute interior angle of the
    // rhombus.
    {
      // Helper lambda function that constructs a rhombus quad with
      // interior acute angle theta.
      auto test_rhombus_quad = [this](Real theta)
      {
        Real ct = std::cos(theta);
        Real st = std::sin(theta);
        std::vector<Point> pts = {
          Point(0, 0, 0),
          Point(1, 0, 0),
          Point(1. + ct, st, 0),
          Point(     ct, st, 0)};
        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
        libmesh_ignore(nodes);
 
        // The expected aspect ratio for the rhombus is 1/sin(theta)
        Real aspect_ratio = elem->quality(ASPECT_RATIO);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0/st, /*actual=*/aspect_ratio, TOLERANCE);
      };
 
      // 2a) Rhombus with interior angle theta=pi/6. The expected
      // aspect ratio in this case is 1/std::sin(pi/6) = 2
      test_rhombus_quad(libMesh::pi / 6);
 
      // 2b) Rhombus with interior angle theta=pi/3. The expected
      // aspect ratio in this case is 1/std::sin(pi/3) = 2/sqrt(3) = 1.155
      test_rhombus_quad(libMesh::pi / 3);
    }
 
    // Case 3) Rectangle QUAD4. The "old" and "new" aspect ratio metrics
    // return the same result for this case, whih is simply the ratio of
    // the longest to shortest side.
    {
      auto test_rectangle_quad = [this](Real a)
      {
        std::vector<Point> pts = {Point(0, 0, 0), Point(a, 0, 0), Point(a, 1, 0), Point(0, 1, 0)};
        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
        libmesh_ignore(nodes);
 
        // The expected aspect ratio for the rectangle is "a"
        Real aspect_ratio = elem->quality(ASPECT_RATIO);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/a, /*actual=*/aspect_ratio, TOLERANCE);
      };
 
      // 3a) Test case twice as long as it is tall
      test_rectangle_quad(2.0);
 
      // 3b) Test no funny numerical stuff with extremely stretched case
      test_rectangle_quad(1e6);
    }
 
    // Case 4) Degenerate QUAD with zero length side. In the "old"
    // aspect ratio metric we'd just get zero (as a stand-in for
    // infinity) for this case since the minimum side length is zero,
    // but using the new metric we get a non-zero value of 2.5. Thus,
    // in the new metric it is possible to compare the aspect ratios
    // of different degenerate QUAD elements rather than simply
    // assigning all such elements a quality value of 0. Degenerate
    // quadrilaterals are sometimes used as a "hack" to avoid having a
    // separate subdomain of TRIs, and (surprisingly) many finite
    // element operations just "work" on such elements.
    {
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 0, 0), Point(0, 1, 0)};
      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
      libmesh_ignore(nodes);
 
      Real aspect_ratio = elem->quality(ASPECT_RATIO);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/2.5, /*actual=*/aspect_ratio, TOLERANCE);
    }
 
    // Case 5) Trapezoid QUAD
    {
      auto test_trapezoid_quad = [this](Real a)
      {
        // 0 <= a <= 1/2
        std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1-a, 1, 0), Point(a, 1, 0)};
        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
        libmesh_ignore(nodes);
 
        Real aspect_ratio = elem->quality(ASPECT_RATIO);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1./(1. - a), /*actual=*/aspect_ratio, TOLERANCE);
      };
 
      // 3a) Test "simple" trapezoid with expected aspect ratio 1.5
      test_trapezoid_quad(1./3);
 
      // 3b) Test "degenerate" trapezoid with one zero length base
      test_trapezoid_quad(0.5);
    }
  }

References libMesh::ASPECT_RATIO, construct_elem(), libMesh::libmesh_ignore(), libMesh::pi, libMesh::QUAD4, libMesh::Real, and libMesh::TOLERANCE.

testQuad4Invertible()

void VolumeTest::testQuad4Invertible ( )

inline

Definition at line 316 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Case 1: Test that rigid body rotations have no effect on the
    // invertibility of the reference element
    {
      // 1a) The reference element rotated into various different different planes.
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0)};
      bool invertible = test_elem_invertible(pts, QUAD4);
      CPPUNIT_ASSERT(invertible);
 
      // 1b) Rotate all points about x-axis by 90 degrees
      Real cost = std::cos(.5*libMesh::pi);
      Real sint = std::sin(.5*libMesh::pi);
      RealTensorValue Rx(1, 0, 0,
                         0, cost, sint,
                         0, sint, cost);
 
      for (auto & pt : pts)
        pt = Rx * pt;
 
      invertible = test_elem_invertible(pts, QUAD4);
      CPPUNIT_ASSERT(invertible);
 
      // 1c) Rotate all points about z-axis by 90 degrees
      RealTensorValue Rz(cost, -sint, 0,
                         sint,  cost, 0,
                         0,        0, 1);
 
      for (auto & pt : pts)
        pt = Rz * pt;
 
      invertible = test_elem_invertible(pts, QUAD4);
      CPPUNIT_ASSERT(invertible);
 
      // 1d) Rotate all points about y-axis by 270 degrees
      RealTensorValue Ry(cost,  0, sint,
                         0,     1, 0,
                         -sint, 0, cost);
 
      for (int cnt=0; cnt<3; ++cnt)
        for (auto & pt : pts)
          pt = Ry * pt;
 
      invertible = test_elem_invertible(pts, QUAD4);
      CPPUNIT_ASSERT(invertible);
    }
 
    // Case 2: Planar quad with top right vertex displaced to the position
    // (alpha, alpha). Some different cases are described below.
    // .) alpha==1: affine case, always invertible
    // .) 1/2 < alpha < 1: planar case, invertible
    // .) alpha<=1/2: planar case but node is now at center of the
    //    element, should give a zero/negative Jacobian on the displaced
    //    Node -> not invertible.
    {
      const Real alpha = .5;
 
      bool invertible =
        test_elem_invertible({Point(0, 0, 0), Point(1, 0, 0), Point(alpha, alpha, 0), Point(0, 1, 0)}, QUAD4);
 
      CPPUNIT_ASSERT(!invertible);
    }
 
    // Case 3) Top right corner is moved to (alpha, 1, 0). Element
    // becomes non-invertible when alpha < 0.
    {
      const Real alpha = -0.25;
 
      bool invertible =
        test_elem_invertible({Point(0, 0, 0), Point(1, 0, 0), Point(alpha, 1, 0), Point(0, 1, 0)}, QUAD4);
 
      CPPUNIT_ASSERT(!invertible);
    }
 
    // Case 4) Degenerate case - all 4 points at same location. This
    // zero-volume element does not have an invertible map.
    {
      const Real alpha = std::log(2);
 
      bool invertible =
        test_elem_invertible({Point(alpha, alpha, alpha),
                   Point(alpha, alpha, alpha),
                   Point(alpha, alpha, alpha),
                   Point(alpha, alpha, alpha)}, QUAD4);
 
      CPPUNIT_ASSERT(!invertible);
    }
  }

References libMesh::pi, libMesh::QUAD4, libMesh::Real, and test_elem_invertible().

testQuad4Jacobian()

void VolumeTest::testQuad4Jacobian ( )

inline

Definition at line 783 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Rhombus QUAD4. This case should have a JACOBIAN and
    // SCALED_JACOBIAN value of sin(theta), where "theta" is the
    // specified amount that we "sheared" the element by on creation.
    // The JACOBIAN and SCALED_JACOBIAN are the same for this element
    // because the edge lengths are all = 1.
    {
      // Helper lambda function that constructs a rhombus quad with
      // interior acute angle theta.
      auto test_rhombus_quad = [this](Real theta)
      {
        Real ct = std::cos(theta);
        Real st = std::sin(theta);
        std::vector<Point> pts = {
          Point(0, 0, 0),
          Point(1, 0, 0),
          Point(1. + ct, st, 0),
          Point(     ct, st, 0)};
        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
        libmesh_ignore(nodes);
 
        Real jac = elem->quality(JACOBIAN);
        Real scaled_jac = elem->quality(SCALED_JACOBIAN);
 
        // Debugging
        // libMesh::out << "jac = " << jac << std::endl;
        // libMesh::out << "scaled_jac = " << scaled_jac << std::endl;
 
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::abs(std::sin(theta)), /*actual=*/jac, TOLERANCE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/std::abs(std::sin(theta)), /*actual=*/scaled_jac, TOLERANCE);
      };
 
      // 2a) Rhombus with interior angle theta=pi/6.
      test_rhombus_quad(libMesh::pi / 6);
 
      // 2b) Rhombus with interior angle theta=pi/3.
      test_rhombus_quad(libMesh::pi / 3);
    }
  }

References construct_elem(), libMesh::JACOBIAN, libMesh::libmesh_ignore(), libMesh::pi, libMesh::QUAD4, libMesh::Real, libMesh::SCALED_JACOBIAN, and libMesh::TOLERANCE.

testQuad4MinMaxAngle()

void VolumeTest::testQuad4MinMaxAngle ( )

inline

Definition at line 680 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Case 1: Test that rigid body rotations of a unit square
    // quadrilateral that have no effect on the quality of the
    // element.
    {
      // Construct unit square QUAD4
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0)};
      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
      libmesh_ignore(nodes);
 
      // 1a) Reference Elem should have min angle == max angle == pi/2
      {
        Real min_angle = elem->quality(MIN_ANGLE);
        Real max_angle = elem->quality(MAX_ANGLE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/min_angle, TOLERANCE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/max_angle, TOLERANCE);
      }
 
      // 1b) Rotate all points about x-axis by 90 degrees
      Real cost = std::cos(.5*libMesh::pi);
      Real sint = std::sin(.5*libMesh::pi);
      RealTensorValue Rx(1, 0, 0,
                         0, cost, sint,
                         0, sint, cost);
 
      for (auto & pt : pts)
        pt = Rx * pt;
 
      {
        Real min_angle = elem->quality(MIN_ANGLE);
        Real max_angle = elem->quality(MAX_ANGLE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/min_angle, TOLERANCE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/max_angle, TOLERANCE);
      }
 
      // 1c) Rotate all points about z-axis by 90 degrees
      RealTensorValue Rz(cost, -sint, 0,
                         sint,  cost, 0,
                         0,        0, 1);
 
      for (auto & pt : pts)
        pt = Rz * pt;
 
      {
        Real min_angle = elem->quality(MIN_ANGLE);
        Real max_angle = elem->quality(MAX_ANGLE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/min_angle, TOLERANCE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/max_angle, TOLERANCE);
      }
 
      // 1d) Rotate all points about y-axis by 270 degrees
      RealTensorValue Ry(cost,  0, sint,
                         0,     1, 0,
                         -sint, 0, cost);
 
      for (int cnt=0; cnt<3; ++cnt)
        for (auto & pt : pts)
          pt = Ry * pt;
 
      {
        Real min_angle = elem->quality(MIN_ANGLE);
        Real max_angle = elem->quality(MAX_ANGLE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/min_angle, TOLERANCE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/max_angle, TOLERANCE);
      }
    }
 
    // Case 2: Rhombus QUAD4. This case should have an interior min
    // angle of theta and an interior max angle of pi-theta, where
    // "theta" is the specified amount that we "sheared" the element
    // by on creation
    {
      // Helper lambda function that constructs a rhombus quad with
      // interior acute angle theta.
      auto test_rhombus_quad = [this](Real theta)
      {
        Real ct = std::cos(theta);
        Real st = std::sin(theta);
        std::vector<Point> pts = {
          Point(0, 0, 0),
          Point(1, 0, 0),
          Point(1. + ct, st, 0),
          Point(     ct, st, 0)};
        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
        libmesh_ignore(nodes);
 
        Real min_angle = elem->quality(MIN_ANGLE);
        Real max_angle = elem->quality(MAX_ANGLE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(180) / libMesh::pi * theta,                 /*actual=*/min_angle, TOLERANCE);
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(180) / libMesh::pi * (libMesh::pi - theta), /*actual=*/max_angle, TOLERANCE);
      };
 
      // 2a) Rhombus with interior angle theta=pi/6.
      test_rhombus_quad(libMesh::pi / 6);
 
      // 2b) Rhombus with interior angle theta=pi/3.
      test_rhombus_quad(libMesh::pi / 3);
    }
  }

References construct_elem(), libMesh::libmesh_ignore(), libMesh::MAX_ANGLE, libMesh::MIN_ANGLE, libMesh::pi, libMesh::QUAD4, libMesh::Real, and libMesh::TOLERANCE.

testQuad4TrueCentroid()

void VolumeTest::testQuad4TrueCentroid ( )

inline

Definition at line 78 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Test Quad4::true_centroid() override
    {
      const Elem & quad4 = ReferenceElem::get(QUAD4);
      Point true_centroid = quad4.true_centroid();
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(0), TOLERANCE*TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(0, true_centroid(1), TOLERANCE*TOLERANCE);
 
      // Compare to centroid computed via generic base class implementation
      Point base_class_centroid = quad4.Elem::true_centroid();
      CPPUNIT_ASSERT(true_centroid.absolute_fuzzy_equals(base_class_centroid, TOLERANCE*TOLERANCE));
    }
 
    // Check that "optimized" Quad4::true_centroid() gives same result
    // as "generic" Elem::true_centroid() on a mesh of 10% distorted elements.
    {
      ReplicatedMesh mesh(*TestCommWorld);
 
      MeshTools::Generation::build_square(mesh,
                                          /*nx=*/3, /*ny=*/3,
                                          /*xmin=*/0., /*xmax=*/1.,
                                          /*ymin=*/0., /*ymax=*/1.,
                                          QUAD4);
 
      MeshTools::Modification::distort(mesh,
                                       /*factor=*/0.1,
                                       /*perturb_boundary=*/false);
 
      for (const auto & elem : mesh.element_ptr_range())
        {
          Point derived_centroid = elem->true_centroid();
          Point base_centroid = elem->Elem::true_centroid();
 
          // Debugging: check results in detail
          // auto flags = libMesh::out.flags();
          // libMesh::out << std::scientific << std::setprecision(16);
          // libMesh::out << "derived_centroid = " << derived_centroid << std::endl;
          // libMesh::out << "base_centroid = " << base_centroid << std::endl;
          // libMesh::out.flags(flags);
 
          CPPUNIT_ASSERT(derived_centroid.absolute_fuzzy_equals(base_centroid, TOLERANCE*TOLERANCE));
 
          Real derived_volume = elem->volume();
          Real base_volume = elem->Elem::volume();
          LIBMESH_ASSERT_FP_EQUAL(base_volume, derived_volume, TOLERANCE*TOLERANCE);
        }
    }
  }

References libMesh::TypeVector< T >::absolute_fuzzy_equals(), libMesh::MeshTools::Generation::build_square(), libMesh::MeshTools::Modification::distort(), libMesh::ReferenceElem::get(), mesh, libMesh::QUAD4, libMesh::Real, TestCommWorld, libMesh::TOLERANCE, and libMesh::Elem::true_centroid().

testQuad4Warpage()

void VolumeTest::testQuad4Warpage ( )

inline

Definition at line 600 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Case 1: Test that rigid body rotations of a unit square
    // quadrilateral that have no effect on the quality of the
    // element.
    {
      // Construct unit square QUAD4
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0)};
      auto [elem, nodes] = this->construct_elem(pts, QUAD4);
      libmesh_ignore(nodes);
 
      // 1a) Any flat element should have warp == 1
      Real warpage = elem->quality(WARP);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/warpage, TOLERANCE);
 
      // 1b) Rotate all points about x-axis by 90 degrees
      Real cost = std::cos(.5*libMesh::pi);
      Real sint = std::sin(.5*libMesh::pi);
      RealTensorValue Rx(1, 0, 0,
                         0, cost, sint,
                         0, sint, cost);
 
      for (auto & pt : pts)
        pt = Rx * pt;
 
      warpage = elem->quality(WARP);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/warpage, TOLERANCE);
 
      // 1c) Rotate all points about z-axis by 90 degrees
      RealTensorValue Rz(cost, -sint, 0,
                         sint,  cost, 0,
                         0,        0, 1);
 
      for (auto & pt : pts)
        pt = Rz * pt;
 
      warpage = elem->quality(WARP);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/warpage, TOLERANCE);
 
      // 1d) Rotate all points about y-axis by 270 degrees
      RealTensorValue Ry(cost,  0, sint,
                         0,     1, 0,
                         -sint, 0, cost);
 
      for (int cnt=0; cnt<3; ++cnt)
        for (auto & pt : pts)
          pt = Ry * pt;
 
      warpage = elem->quality(WARP);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/1.0, /*actual=*/warpage, TOLERANCE);
    }
 
    // Case 2: Unit square quadrilateral with Node 2 displaced by a distance h in the z-direction.
    {
      auto test_warped_quad = [this](Real h)
      {
        std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, h), Point(0, 1, 0)};
        auto [elem, nodes] = this->construct_elem(pts, QUAD4);
        libmesh_ignore(nodes);
 
        Real warpage = elem->quality(WARP);
        // libMesh::out << "QUAD with node 3 displaced by h = "
        //              << h << " in the z-direction , warpage = " << warpage
        //              << std::endl;
        LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(1) / (h*h + 1), /*actual=*/warpage, TOLERANCE);
      };
 
      // For h = 0.1, the expected warpage value is ~ 0.991, so not
      // much different than the value of 1.0 for a flat element.
      test_warped_quad(0.1);
 
      // For h = 0.3, the expected warpage value is ~ 0.917431, which
      // is pretty close to the lower bound (0.9) of "acceptable"
      // warpage, as suggested in the Cubit manual.
      test_warped_quad(0.3);
    }
  }

References construct_elem(), libMesh::libmesh_ignore(), libMesh::pi, libMesh::QUAD4, libMesh::Real, libMesh::TOLERANCE, and libMesh::WARP.

testTet4DihedralAngle()

void VolumeTest::testTet4DihedralAngle ( )

inline

Definition at line 826 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Construct a tetrahedron whose projection into the x-y plane looks like the unit square,
    // and where the apex node is just slightly out of the x-y plane. This element should have
    // reasonable MIN,MAX_ANGLE values but MIN,MAX_DIHEDRAL angles of nearly 0. This is an
    // example that demonstrates one should not solely use edge angles to determine Elem quality
    // in 3D.
    std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0), Point(1, 1, 0.01)};
    auto [elem, nodes] = this->construct_elem(pts, TET4);
    libmesh_ignore(nodes);
 
    Real min_angle = elem->quality(MIN_ANGLE);
    Real max_angle = elem->quality(MAX_ANGLE);
    Real min_dihedral_angle = elem->quality(MIN_DIHEDRAL_ANGLE);
    Real max_dihedral_angle = elem->quality(MAX_DIHEDRAL_ANGLE);
 
    // Debugging
    // libMesh::out << "Squashed Tet4 min_angle = " << min_angle << " degrees" << std::endl;
    // libMesh::out << "Squashed Tet4 max_angle = " << max_angle << " degrees"  << std::endl;
    // libMesh::out << "Squashed Tet4 min_dihedral_angle = " << min_dihedral_angle << " degrees" << std::endl;
    // libMesh::out << "Squashed Tet4 max_dihedral_angle = " << max_dihedral_angle << " degrees" << std::endl;
 
    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/44.9985676771277, /*actual=*/min_angle, TOLERANCE);
    LIBMESH_ASSERT_FP_EQUAL(/*expected=*/90, /*actual=*/max_angle, TOLERANCE);
 
    // Assert that both min and max dihedral angles are less than 1 degree (a very low quality element)
    CPPUNIT_ASSERT_LESS(Real(1), min_dihedral_angle);
    CPPUNIT_ASSERT_LESS(Real(1), max_dihedral_angle);
  }

References construct_elem(), libMesh::libmesh_ignore(), libMesh::MAX_ANGLE, libMesh::MAX_DIHEDRAL_ANGLE, libMesh::MIN_ANGLE, libMesh::MIN_DIHEDRAL_ANGLE, libMesh::Real, libMesh::TET4, and libMesh::TOLERANCE.

testTet4Jacobian()

void VolumeTest::testTet4Jacobian ( )

inline

Definition at line 858 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    {
      // Same element as in testTet4DihedralAngle(). Here we verify that
      // the element is low quality since the SCALED_JACOBIAN is < O(h)
      // as h -> 0, and h can be arbitrarily small in the squashed element.
      Real h = Real(1)/100;
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0), Point(1, 1, h)};
      auto [elem, nodes] = this->construct_elem(pts, TET4);
      libmesh_ignore(nodes);
 
      Real jac = elem->quality(JACOBIAN);
      Real scaled_jac = elem->quality(SCALED_JACOBIAN);
 
      // Debugging
      // libMesh::out << "Squashed Tet4 jac = " << jac << std::endl;
      // libMesh::out << "Squashed Tet4 scaled_jac = " << scaled_jac << std::endl;
 
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/h, /*actual=*/jac, TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/h/(h*h+1)/std::sqrt(h*h+2),
                              /*actual=*/scaled_jac, TOLERANCE * TOLERANCE);
    }
 
    {
      // A representative Tet generated by calling build_cube(). This
      // element has minimum interior edge angle of approximately
      // 35.26 degrees, so it it is not a simple refinement of the
      // reference element. It is not located at the origin, but still
      // has a simple to compute value for the JACOBIAN and
      // SCALED_JACOBIAN quality metrics.
      std::vector<Point> pts = {
        Point(2.5, 2.5, 4.5),
        Point(3.5, 1.5, 5.5),
        Point(1.5, 1.5, 5.5),
        Point(2.5, 2.5, 5.5)
      };
      auto [elem, nodes] = this->construct_elem(pts, TET4);
      libmesh_ignore(nodes);
 
      Real jac = elem->quality(JACOBIAN);
      Real scaled_jac = elem->quality(SCALED_JACOBIAN);
 
      // Debugging
      // libMesh::out << "build_cube() Tet4 jac = " << jac << std::endl;
      // libMesh::out << "build_cube() Tet4 scaled_jac = " << scaled_jac << std::endl;
 
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/2,                          /*actual=*/jac, TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(1)/std::sqrt(Real(6)), /*actual=*/scaled_jac, TOLERANCE);
    }
  }

References construct_elem(), libMesh::JACOBIAN, libMesh::libmesh_ignore(), libMesh::Real, libMesh::SCALED_JACOBIAN, libMesh::TET4, and libMesh::TOLERANCE.

testTri3AspectRatio()

void VolumeTest::testTri3AspectRatio ( )

inline

Definition at line 555 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // Case 1) Reference TRI3
    {
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0)};
      auto [elem, nodes] = this->construct_elem(pts, TRI3);
      libmesh_ignore(nodes);
 
      // Compute the aspect ratio for the reference Tri
      // The expected value is ~ 1.44338
      Real aspect_ratio = elem->quality(ASPECT_RATIO);
      // libMesh::out << "Unit TRI3 aspect ratio = " << aspect_ratio << std::endl;
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(5)/2/std::sqrt(Real(3)), /*actual=*/aspect_ratio, TOLERANCE);
    }
 
    // Case 2) Equilateral TRI3
    {
      std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(0.5, 0.5*std::sqrt(3), 0)};
      auto [elem, nodes] = this->construct_elem(pts, TRI3);
      libmesh_ignore(nodes);
 
      // Compute the aspect ratio for the reference Tri
      // The expected value is 1.0
      Real aspect_ratio = elem->quality(ASPECT_RATIO);
      // libMesh::out << "Equilateral TRI3 aspect ratio = " << aspect_ratio << std::endl;
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(1), /*actual=*/aspect_ratio, TOLERANCE);
    }
 
    // Case 3) Reference TRI3 with one leg length = L >> 1
    {
      Real L = 10.;
      std::vector<Point> pts = {Point(0, 0, 0), Point(L, 0, 0), Point(0, 1, 0)};
      auto [elem, nodes] = this->construct_elem(pts, TRI3);
      libmesh_ignore(nodes);
 
      // Compute the aspect ratio for the reference Tri
      // The expected value is ~ 11.5759
      Real aspect_ratio = elem->quality(ASPECT_RATIO);
      // libMesh::out << "TRI3 with leg length L = " << L << ", aspect ratio = " << aspect_ratio << std::endl;
      LIBMESH_ASSERT_FP_EQUAL(/*expected=*/Real(1)/std::sqrt(Real(3)) * (Real(2)*L + Real(1)/(2*L)), /*actual=*/aspect_ratio, TOLERANCE);
    }
  }

References libMesh::ASPECT_RATIO, construct_elem(), libMesh::libmesh_ignore(), libMesh::Real, libMesh::TOLERANCE, and libMesh::TRI3.

testTri3TrueCentroid()

void VolumeTest::testTri3TrueCentroid ( )

inline

Definition at line 65 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    // The true_centroid() == vertex_average() == (1/3, 1/3) for reference Tri3
    {
      const Elem & tri3 = ReferenceElem::get(TRI3);
      Point true_centroid = tri3.true_centroid();
      LIBMESH_ASSERT_FP_EQUAL(Real(1)/3, true_centroid(0), TOLERANCE*TOLERANCE);
      LIBMESH_ASSERT_FP_EQUAL(Real(1)/3, true_centroid(1), TOLERANCE*TOLERANCE);
    }
  }

References libMesh::ReferenceElem::get(), libMesh::Real, libMesh::TOLERANCE, libMesh::TRI3, and libMesh::Elem::true_centroid().

testTwistedVolume()

void VolumeTest::testTwistedVolume ( )

inline

Definition at line 175 of file volume_test.C.

  {
    LOG_UNIT_TEST;
 
    ReplicatedMesh mesh(*TestCommWorld);
 
    // Build an element type that will fall back on our generic
    // quadrature-based Elem::volume()
    MeshTools::Generation::build_cube(mesh,
                                      /*nelem=*/1, /*nelem=*/1, /*nelem=*/1,
                                      /*xmin=*/-1, /*xmax=*/1,
                                      /*ymin=*/-1, /*ymax=*/1,
                                      /*zmin=*/-1, /*zmax=*/1,
                                      PRISM21);
 
    // Pick an element and twist it
    Elem * prism6 = mesh.elem_ptr(0);
    prism6->point(1) *= -1;
    prism6->point(1) += 2*prism6->point(0);
 
    // The real test here is that volume() doesn't throw
    const Real vol = prism6->volume();
    const Real gold_vol = 3+Real(5)/9;
    CPPUNIT_ASSERT_LESS(TOLERANCE, std::abs(vol-gold_vol));
  }

References libMesh::MeshTools::Generation::build_cube(), libMesh::MeshBase::elem_ptr(), mesh, libMesh::Elem::point(), libMesh::PRISM21, libMesh::Real, TestCommWorld, libMesh::TOLERANCE, and libMesh::Elem::volume().

---

The documentation for this class was generated from the following file:

• tests/geom/volume_test.C