Merge pull request #4287 from pbehne/variational_smoother_tets

Variational smoother: tet support

Author: roystgnr

include/mesh/mesh_tools.h

@@ -347,6 +347,19 @@ void find_nodal_neighbors(const MeshBase & mesh,
                           const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
                           std::vector<const Node *> & neighbors);
 
+/**
+ * Given a mesh and a node in the mesh, the vector will be filled with
+ * every node directly attached to the given one. IF NO nodal neighbors are found,
+ * all nodes on the containing side are treated as neighbors. This is useful
+ * when the node does not lie on an edge, such as the central face node in
+ * HEX27, TET14, etc. elements.
+ */
+void find_nodal_or_face_neighbors(
+    const MeshBase & mesh,
+    const Node & node,
+    const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
+    std::vector<const Node *> & neighbors);
+
 /**
  * Given a mesh hanging_nodes will be filled with an associative array keyed off the
  * global id of all the hanging nodes in the mesh.  It will hold an array of the

include/systems/variational_smoother_constraint.h

@@ -433,19 +433,6 @@ class VariationalSmootherConstraint : public System::Constraint
    */
   void constrain_node_to_line(const Node & node, const Point & line_vec);
 
-  /**
-   * Given a mesh and a node in the mesh, the vector will be filled with
-   * every node directly attached to the given one. IF NO neighbors are found,
-   * all nodes on the containing side are treated as neighbors. This is useful
-   * when the node does not lie on an edge, such as the central face node in
-   * HEX27 elements.
-   */
-  static void find_nodal_or_face_neighbors(
-      const MeshBase & mesh,
-      const Node & node,
-      const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
-      std::vector<const Node *> & neighbors);
-
   /**
    * Determines whether two neighboring nodes share a common boundary id.
    * @param boundary_node The first of the two prospective nodes.

src/mesh/mesh_tools.C

@@ -1061,6 +1061,45 @@ void find_nodal_neighbors(const MeshBase &,
   find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
 }
 
+void find_nodal_or_face_neighbors(
+    const MeshBase & mesh,
+    const Node & node,
+    const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
+    std::vector<const Node *> & neighbors)
+{
+  // Find all the nodal neighbors... that is the nodes directly connected
+  // to this node through one edge.
+  find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
+
+  // If no neighbors are found, use all nodes on the containing side as
+  // neighbors.
+  if (!neighbors.size())
+    {
+      // Grab the element containing node
+      const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
+      // Find the element side containing node
+      for (const auto &side : elem->side_index_range())
+        {
+          const auto &nodes_on_side = elem->nodes_on_side(side);
+          const auto it =
+              std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
+                return elem->node_id(local_node_id) == node.id();
+              });
+
+          if (it != nodes_on_side.end())
+            {
+              for (const auto &local_node_id : nodes_on_side)
+                // No need to add node itself as a neighbor
+                if (const auto *node_ptr = elem->node_ptr(local_node_id);
+                    *node_ptr != node)
+                  neighbors.push_back(node_ptr);
+              break;
+            }
+        }
+    }
+  libmesh_assert(neighbors.size());
+}
+
 
 
 void find_hanging_nodes_and_parents(const MeshBase & mesh,

src/systems/variational_smoother_constraint.C

@@ -562,45 +562,6 @@ void VariationalSmootherConstraint::constrain_node_to_line(const Node & node,
     }
 }
 
-void VariationalSmootherConstraint::find_nodal_or_face_neighbors(
-    const MeshBase & mesh,
-    const Node & node,
-    const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
-    std::vector<const Node *> & neighbors)
-{
-  // Find all the nodal neighbors... that is the nodes directly connected
-  // to this node through one edge.
-  MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
-
-  // If no neighbors are found, use all nodes on the containing side as
-  // neighbors.
-  if (!neighbors.size())
-    {
-      // Grab the element containing node
-      const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
-      // Find the element side containing node
-      for (const auto &side : elem->side_index_range())
-        {
-          const auto &nodes_on_side = elem->nodes_on_side(side);
-          const auto it =
-              std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
-                return elem->node_id(local_node_id) == node.id();
-              });
-
-          if (it != nodes_on_side.end())
-            {
-              for (const auto &local_node_id : nodes_on_side)
-                // No need to add node itself as a neighbor
-                if (const auto *node_ptr = elem->node_ptr(local_node_id);
-                    *node_ptr != node)
-                  neighbors.push_back(node_ptr);
-              break;
-            }
-        }
-    }
-  libmesh_assert(neighbors.size());
-}
-
 // Utility function to determine whether two nodes share a boundary ID.
 // The motivation for this is that a sliding boundary node on a triangular
 // element can have a neighbor boundary node in the same element that is not
@@ -671,7 +632,7 @@ VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
   // to this node through one edge, or if none exists, use other nodes on the
   // containing face
   std::vector<const Node *> neighbors;
-  VariationalSmootherConstraint::find_nodal_or_face_neighbors(
+  MeshTools::find_nodal_or_face_neighbors(
       mesh, node, nodes_to_elem_map, neighbors);
 
   // Each constituent set corresponds to neighbors sharing a face on the
@@ -757,7 +718,7 @@ VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
   // to this node through one edge, or if none exists, use other nodes on the
   // containing face
   std::vector<const Node *> neighbors;
-  VariationalSmootherConstraint::find_nodal_or_face_neighbors(
+  MeshTools::find_nodal_or_face_neighbors(
       mesh, node, nodes_to_elem_map, neighbors);
 
   // Each constituent set corresponds to neighbors sharing a face on the

src/systems/variational_smoother_system.C

@@ -923,6 +923,7 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type)
   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;
 
@@ -1055,7 +1056,6 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type)
       // Solving for s: s = (3 sqrt(2) v)^(1/3), where v is the volume of the
       // non-optimal reference element.
 
-      const Real sqrt_2 = std::sqrt(Real(2));
       // Side length that preserves the volume of the reference element
       const auto side_length = std::pow(3. * sqrt_2 * ref_vol, 1. / 3.);
       // Pyramid height with the property that all faces are equilateral triangles
@@ -1110,6 +1110,67 @@ VariationalSmootherSystem::get_target_elem(const ElemType & type)
 
     } // 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())