libMesh::ReferenceElem Namespace Reference

This namespace implements singleton reference elements for each fundamental element type supported by libMesh. More...

Functions
const Elem &	get (const ElemType type_in)
std::pair< std::unique_ptr< Elem >, std::vector< std::unique_ptr< Node > > >	ideal_target (const ElemType type)

Detailed Description

This namespace implements singleton reference elements for each fundamental element type supported by libMesh.

Author

Benjamin S. Kirk

Date

2013

Namespace providing access to reference geometric element types.

Function Documentation

get()

const Elem & libMesh::ReferenceElem::get

(

const ElemType

type_in

)

Returns

A constant reference to the reference element of the user-requested type.

Definition at line 258 of file reference_elem.C.

{
  ElemType base_type = type_in;
 
  // For shell elements, use non shell type as the base type
  if (type_in == TRISHELL3)
    base_type = TRI3;
 
  if (type_in == QUADSHELL4)
    base_type = QUAD4;
 
  if (type_in == QUADSHELL8)
    base_type = QUAD8;
 
  if (type_in == QUADSHELL9)
    base_type = QUAD9;
 
  init_ref_elem_table();
 
  // Throw an error if the user asked for an ElemType that we don't
  // have a reference element for.
  libmesh_error_msg_if(ref_elem_map[base_type] == nullptr || type_in == INVALID_ELEM,
                       "No reference elem data available for ElemType " << type_in
                       << " = " << Utility::enum_to_string(type_in) << ".");
 
  return *ref_elem_map[base_type];
}

References libMesh::Utility::enum_to_string(), libMesh::INVALID_ELEM, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::QUADSHELL8, libMesh::QUADSHELL9, libMesh::TRI3, and libMesh::TRISHELL3.

Referenced by NodalNeighborsTest::do_test(), main(), libMesh::Elem::reference_elem(), SideVertexAverageNormalTest::testEdge2(), SideVertexAverageNormalTest::testEdge3(), VolumeTest::testEdge4Invertible(), SideVertexAverageNormalTest::testHexes(), WhichNodeAmITest::testHexes(), SideVertexAverageNormalTest::testPrisms(), WhichNodeAmITest::testPrisms(), VolumeTest::testPyramid5TrueCentroid(), SideVertexAverageNormalTest::testPyramids(), WhichNodeAmITest::testPyramids(), VolumeTest::testQuad4TrueCentroid(), SideVertexAverageNormalTest::testQuads(), SideVertexAverageNormalTest::testTets(), WhichNodeAmITest::testTets(), VolumeTest::testTri3TrueCentroid(), SideVertexAverageNormalTest::testTris(), and MeshSmootherTest::testVariationalSmoother().

ideal_target()

std::pair< std::unique_ptr< Elem >, std::vector< std::unique_ptr< Node > > > libMesh::ReferenceElem::ideal_target

(

const ElemType

type

)

Returns

A freshly built "ideal" (regular) element of the given type, i.e. the optimally-shaped element that a mesh optimizer targets: an equilateral triangle, regular tetrahedron, etc., sized to the volume of the reference element. For element types that have no distinct ideal shape (e.g. quads and hexes, whose reference element is already regular), this returns a copy of the reference element.

The returned Elem holds pointers into the returned Nodes, so the caller must keep the Node vector alive for at least as long as the Elem.

Definition at line 287 of file reference_elem.C.

{
  // Build target element
  auto target_elem = Elem::build(type);
 
  // Volume of reference element
  const auto ref_vol = target_elem->reference_elem()->volume();
 
  // Update the nodes of the target element, depending on type
  const Real sqrt_2 = std::sqrt(Real(2));
  const Real sqrt_3 = std::sqrt(Real(3));
  std::vector<std::unique_ptr<Node>> owned_nodes;
 
  const auto type_str = Utility::enum_to_string(type);
 
  // Elems deriving from Tri
  if (type_str.compare(0, 3, "TRI") == 0)
    {
 
      // The target element will be an equilateral triangle with area equal to
      // the area of the reference element.
 
      // Equilateral triangle side length preserving area of the reference element
      const auto side_length = std::sqrt(4. / sqrt_3 * ref_vol);
 
      // Define the nodal locations of the vertices
      const auto & s = side_length;
      //                                         x        y                  node_id
      owned_nodes.emplace_back(Node::build(Point(0.,      0.),               0));
      owned_nodes.emplace_back(Node::build(Point(s,       0.),               1));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s), 2));
 
      switch (type)
        {
            case TRI3:
            case TRISHELL3: {
              // Nothing to do here, vertices already added above
              break;
            }
 
            case TRI6:
            case TRI7: {
              // Define the midpoint nodes of the equilateral triangle
              //                                         x         y                   node_id
              owned_nodes.emplace_back(Node::build(Point(0.50 * s, 0.00),              3));
              owned_nodes.emplace_back(Node::build(Point(0.75 * s, 0.25 * sqrt_3 * s), 4));
              owned_nodes.emplace_back(Node::build(Point(0.25 * s, 0.25 * sqrt_3 * s), 5));
 
              if (type == TRI7)
                {
                  // Define the interior "bubble" node at the centroid
                  const auto & on = owned_nodes;
                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 6));
                }
 
              break;
            }
 
          default:
            libmesh_error_msg("Unsupported triangular element: " << type_str);
            break;
        }
    } // if Tri
 
  // Elems deriving from Prism
  else if (type_str.compare(0, 5, "PRISM") == 0)
    {
 
      // The target element will be a prism with an equilateral triangular
      // base with volume equal to the volume of the reference element.
 
      // For an equilateral triangular base with side length s, the
      // base area is s^2 * sqrt(3) / 4.
      // The prism height that will result in equal face areas is
      // s * sqrt(3) / 4. We choose s such that the target element has
      // the same volume as the reference element:
      // v = (s^2 * sqrt(3) / 4) * (s * sqrt(3) / 4) = 3 * s^3 / 4
      // --> s = (16 * v / 3)^(1/3)
      // I have no particular motivation for imposing equal face areas,
      // so this can be updated if a more `optimal` target prism is
      // identified.
 
      // Side length that preserves the volume of the reference element
      const auto side_length = std::cbrt(16. * ref_vol / 3.);
      // Prism height with the property that all faces have equal area
      const auto target_height = 0.25 * side_length * sqrt_3;
 
      const auto & s = side_length;
      const auto & h = target_height;
      //                                         x        y                 z    node_id
      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               h),  3));
      owned_nodes.emplace_back(Node::build(Point(s,       0.,               h),  4));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, h),  5));
 
      if (type == PRISM15 || type == PRISM18 || type == PRISM20 || type == PRISM21)
        {
          // Define the edge midpoint nodes of the prism
          const auto & on = owned_nodes;
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 6));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 7));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 8));
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 9));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[5]) / 2.), 11));
          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
          owned_nodes.emplace_back(Node::build(Point((*on[4] + *on[5]) / 2.), 13));
          owned_nodes.emplace_back(Node::build(Point((*on[5] + *on[3]) / 2.), 14));
 
          if (type == PRISM18 || type == PRISM20 || type == PRISM21)
            {
              // Define the rectangular face midpoint nodes of the prism
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3] + *on[4]) / 4.), 15));
              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4] + *on[5]) / 4.), 16));
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3] + *on[5]) / 4.), 17));
 
              if (type == PRISM20 || type == PRISM21)
                {
                  // Define the triangular face midpoint nodes of the prism
                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 18));
                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4] + *on[5]) / 3.), 19));
 
                  if (type == PRISM21)
                    // Define the interior point of the prism
                    owned_nodes.emplace_back(Node::build(Point((*on[9] + *on[10] + *on[11]) / 3.), 20));
 
                }
            }
        }
 
      else if (type != PRISM6)
        libmesh_error_msg("Unsupported prism element: " << type_str);
 
    } // if Prism
 
  // Elems deriving from Pyramid
  else if (type_str.compare(0, 7, "PYRAMID") == 0)
    {
 
      // The target element is a pyramid with an square base and
      // equilateral triangular sides with volume equal to the volume of the
      // reference element.
 
      // A pyramid with square base sidelength s and equilateral triangular
      // sides has height h = s / sqrt(2).
      // The volume is v = s^2 h / 3 = s^3 / ( 3 sqrt(2)).
      // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
      // non-optimal reference element.
 
      // Side length that preserves the volume of the reference element
      const auto side_length = std::cbrt(3. * sqrt_2 * ref_vol);
      // Pyramid height with the property that all faces are equilateral triangles
      const auto target_height = side_length / sqrt_2;
 
      const auto & s = side_length;
      const auto & h = target_height;
 
      //                                         x        y        z    node_id
      owned_nodes.emplace_back(Node::build(Point(0.,      0.,      0.), 0));
      owned_nodes.emplace_back(Node::build(Point(s,       0.,      0.), 1));
      owned_nodes.emplace_back(Node::build(Point(s,       s,       0.), 2));
      owned_nodes.emplace_back(Node::build(Point(0.,      s,       0.), 3));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * s, h),  4));
 
      if (type == PYRAMID13 || type == PYRAMID14 || type == PYRAMID18)
        {
          const auto & on = owned_nodes;
          // Define the edge midpoint nodes of the pyramid
 
          // Base node to base node midpoint nodes
          owned_nodes.emplace_back(Node::build((*on[0] + *on[1]) / 2., 5));
          owned_nodes.emplace_back(Node::build((*on[1] + *on[2]) / 2., 6));
          owned_nodes.emplace_back(Node::build((*on[2] + *on[3]) / 2., 7));
          owned_nodes.emplace_back(Node::build((*on[3] + *on[0]) / 2., 8));
 
          // Base node to apex node midpoint nodes
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[4]) / 2.), 9));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[4]) / 2.), 10));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[4]) / 2.), 11));
          owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[4]) / 2.), 12));
 
          if (type == PYRAMID14 || type == PYRAMID18)
            {
              // Define the square face midpoint node of the pyramid
              owned_nodes.emplace_back(
                  Node::build(Point((*on[0] + *on[1] + *on[2] + *on[3]) / 4.), 13));
 
              if (type == PYRAMID18)
                {
                  // Define the triangular face nodes
                  owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[4]) / 3.), 14));
                  owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[4]) / 3.), 15));
                  owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3] + *on[4]) / 3.), 16));
                  owned_nodes.emplace_back(Node::build(Point((*on[3] + *on[0] + *on[4]) / 3.), 17));
                }
            }
        }
 
      else if (type != PYRAMID5)
        libmesh_error_msg("Unsupported pyramid element: " << type_str);
 
    } // if Pyramid
 
  // Elems deriving from Tet
  else if (type_str.compare(0, 3, "TET") == 0)
    {
 
      // The ideal target element is a a regular tet with equilateral
      // triangles for all faces, with volume equal to the volume of the
      // reference element.
 
      // The volume of a tet is given by v = b * h / 3, where b is the area of
      // the base face and h is the height of the apex node. The area of an
      // equilateral triangle with side length s is b = sqrt(3) s^2 / 4.
      // For all faces to have side length s, the height of the apex node is
      // h = sqrt(2/3) * s. Then the volume is v = sqrt(2) * s^3 / 12.
      // Solving for s, the side length that will preserve the volume of the
      // reference element is s = (6 * sqrt(2) * v)^(1/3), where v is the volume
      // of the non-optimal reference element (i.e., a right tet).
 
      // Side length that preserves the volume of the reference element
      const auto side_length = std::cbrt(6. * sqrt_2 * ref_vol);
      // tet height with the property that all faces are equilateral triangles
      const auto target_height = sqrt_2 / sqrt_3 * side_length;
 
      const auto & s = side_length;
      const auto & h = target_height;
 
      // For regular tet
      //                                         x        y                z     node_id
      owned_nodes.emplace_back(Node::build(Point(0.,      0.,               0.), 0));
      owned_nodes.emplace_back(Node::build(Point(s,       0.,               0.), 1));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, 0.5 * sqrt_3 * s, 0.), 2));
      owned_nodes.emplace_back(Node::build(Point(0.5 * s, sqrt_3 / 6. * s,  h),  3));
 
      if (type == TET10 || type == TET14)
        {
          const auto & on = owned_nodes;
          // Define the edge midpoint nodes of the tet
 
          // Base node to base node midpoint nodes
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1]) / 2.), 4));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2]) / 2.), 5));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[0]) / 2.), 6));
          // Base node to apex node midpoint nodes
          owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[3]) / 2.), 7));
          owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[3]) / 2.), 8));
          owned_nodes.emplace_back(Node::build(Point((*on[2] + *on[3]) / 2.), 9));
 
          if (type == TET14)
            {
              // Define the face midpoint nodes of the tet
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[2]) / 3.), 10));
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[1] + *on[3]) / 3.), 11));
              owned_nodes.emplace_back(Node::build(Point((*on[1] + *on[2] + *on[3]) / 3.), 12));
              owned_nodes.emplace_back(Node::build(Point((*on[0] + *on[2] + *on[3]) / 3.), 13));
            }
        }
 
      else if (type != TET4)
        libmesh_error_msg("Unsupported tet element: " << type_str);
 
    } // if Tet
 
  // Set the target_elem equal to the reference elem
  else
    for (const auto & node : target_elem->reference_elem()->node_ref_range())
      owned_nodes.emplace_back(Node::build(node, node.id()));
 
  // Set nodes of target element
  for (const auto & node_ptr : owned_nodes)
    target_elem->set_node(node_ptr->id(), node_ptr.get());
 
  libmesh_assert(relative_fuzzy_equals(target_elem->volume(), ref_vol, TOLERANCE));
 
  return std::make_pair(std::move(target_elem), std::move(owned_nodes));
}

References libMesh::Elem::build(), libMesh::Node::build(), libMesh::Utility::enum_to_string(), libMesh::libmesh_assert(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::Real, libMesh::relative_fuzzy_equals(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TOLERANCE, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and libMesh::TRISHELL3.

Referenced by libMesh::VariationalSmootherSystem::prepare_for_smoothing(), and libMesh::Elem::quality().