Merge pull request #4218 from pbehne/variational_smoother_tri6

roystgnr · web-flow · commit 231f9e666cc7 · 2025-07-31T08:25:15.000-05:00
VariationalMeshSmoother: Support for higher order elements
diff --git a/include/systems/variational_smoother_constraint.h b/include/systems/variational_smoother_constraint.h
@@ -57,7 +57,7 @@ class PointConstraint
    * @param point The point defining the constraint.
    * @param tol The tolerance to use for numerical comparisons.
    */
-  PointConstraint(const Point &point, const Real &tol = TOLERANCE);
+  PointConstraint(const Point &point, const Real &tol = TOLERANCE * TOLERANCE);
 
   /**
    * Comparison operator for ordering PointConstraint objects.
@@ -132,7 +132,7 @@ class LineConstraint
    * @param tol The tolerance to use for numerical comparisons.
    */
   LineConstraint(const Point &point, const Point &direction,
-                 const Real &tol = TOLERANCE);
+                 const Real &tol = TOLERANCE * TOLERANCE);
 
   /**
    * Comparison operator for ordering LineConstraint objects.
@@ -232,7 +232,7 @@ class PlaneConstraint
    * @param tol The tolerance to use for numerical comparisons.
    */
   PlaneConstraint(const Point &point, const Point &normal,
-                  const Real &tol = TOLERANCE);
+                  const Real &tol = TOLERANCE * TOLERANCE);
 
   /**
    * Comparison operator for ordering PlaneConstraint objects.
@@ -421,6 +421,19 @@ 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.
diff --git a/src/systems/variational_smoother_constraint.C b/src/systems/variational_smoother_constraint.C
@@ -521,6 +521,45 @@ 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
@@ -588,9 +627,11 @@ VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
 {
 
   // Find all the nodal neighbors... that is the nodes directly connected
-  // to this node through one edge
+  // to this node through one edge, or if none exists, use other nodes on the
+  // containing face
   std::vector<const Node *> neighbors;
-  MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
+  VariationalSmootherConstraint::find_nodal_or_face_neighbors(
+      mesh, node, nodes_to_elem_map, neighbors);
 
   // Each constituent set corresponds to neighbors sharing a face on the
   // subdomain boundary
@@ -600,8 +641,8 @@ VariationalSmootherConstraint::get_neighbors_for_subdomain_constraint(
     {
       // Determine whether the neighbor is on the subdomain boundary
       // First, find the common elements that both node and neigh belong to
-      const auto & elems_containing_node = nodes_to_elem_map.at(node.id());
-      const auto & elems_containing_neigh = nodes_to_elem_map.at(neigh->id());
+      const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
+      const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
       const Elem * common_elem = nullptr;
       for (const auto * neigh_elem : elems_containing_neigh)
         {
@@ -672,9 +713,11 @@ VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
 {
 
   // Find all the nodal neighbors... that is the nodes directly connected
-  // to this node through one edge
+  // to this node through one edge, or if none exists, use other nodes on the
+  // containing face
   std::vector<const Node *> neighbors;
-  MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
+  VariationalSmootherConstraint::find_nodal_or_face_neighbors(
+      mesh, node, nodes_to_elem_map, neighbors);
 
   // Each constituent set corresponds to neighbors sharing a face on the
   // boundary
@@ -689,8 +732,8 @@ VariationalSmootherConstraint::get_neighbors_for_boundary_constraint(
 
       // Determine whether nodes share a common boundary id
       // First, find the common element that both node and neigh belong to
-      const auto & elems_containing_node = nodes_to_elem_map.at(node.id());
-      const auto & elems_containing_neigh = nodes_to_elem_map.at(neigh->id());
+      const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
+      const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
       const Elem * common_elem = nullptr;
       for (const auto * neigh_elem : elems_containing_neigh)
         {
diff --git a/src/systems/variational_smoother_system.C b/src/systems/variational_smoother_system.C
@@ -19,6 +19,7 @@
 
 #include "libmesh/elem.h"
 #include "libmesh/face_tri3.h"
+#include "libmesh/face_tri6.h"
 #include "libmesh/fe_base.h"
 #include "libmesh/fe_interface.h"
 #include "libmesh/fem_context.h"
@@ -113,90 +114,121 @@ void VariationalSmootherSystem::prepare_for_smoothing()
   std::map<ElemType, std::vector<Real>> target_elem_inverse_jacobian_dets;
 
   for (const auto * elem : mesh.active_local_element_ptr_range())
-  {
-    femcontext.pre_fe_reinit(*this, elem);
-    femcontext.elem_fe_reinit();
-
-    // Add target element info, if applicable
-    if (_target_inverse_jacobians.find(elem->type()) == _target_inverse_jacobians.end())
     {
-      // Create FEMap to compute target_element mapping information
-      FEMap fe_map_target;
-
-      // pre-request mapping derivatives
-      const auto & dxyzdxi = fe_map_target.get_dxyzdxi();
-      const auto & dxyzdeta = fe_map_target.get_dxyzdeta();
-      // const auto & dxyzdzeta = fe_map_target.get_dxyzdzeta();
+      femcontext.pre_fe_reinit(*this, elem);
+      femcontext.elem_fe_reinit();
 
-      const auto & qrule_points = femcontext.get_element_qrule().get_points();
-      const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
-      const auto nq_points = femcontext.get_element_qrule().n_points();
-
-      // If the target element is the reference element, Jacobian matrix is
-      // identity, det of inverse is 1. This will only be overwritten if a
-      // different target elemen is explicitly specified.
-      target_elem_inverse_jacobian_dets[elem->type()] =
-          std::vector<Real>(nq_points, 1.0);
-
-      switch (elem->type())
+      // Add target element info, if applicable
+      if (_target_inverse_jacobians.find(elem->type()) == _target_inverse_jacobians.end())
         {
-          case TRI3:
+          // Create FEMap to compute target_element mapping information
+          FEMap fe_map_target;
+
+          // pre-request mapping derivatives
+          const auto & dxyzdxi = fe_map_target.get_dxyzdxi();
+          const auto & dxyzdeta = fe_map_target.get_dxyzdeta();
+          // const auto & dxyzdzeta = fe_map_target.get_dxyzdzeta();
+
+          const auto & qrule_points = femcontext.get_element_qrule().get_points();
+          const auto & qrule_weights = femcontext.get_element_qrule().get_weights();
+          const auto nq_points = femcontext.get_element_qrule().n_points();
+
+          // If the target element is the reference element, Jacobian matrix is
+          // identity, det of inverse is 1. This will only be overwritten if a
+          // different target elemen is explicitly specified.
+          target_elem_inverse_jacobian_dets[elem->type()] =
+              std::vector<Real>(nq_points, 1.0);
+
+          // Elems deriving from Tri
+          if (elem->type() == TRI3 || elem->type() == TRI6)
             {
-              // Build target element: an equilateral triangle
-              Tri3 target_elem;
 
-              // equilateral triangle side length that preserves area of reference
-              // element
+              // The target element will be an equilateral triangle with area equal to
+              // the area of the reference element.
+
+              // First, let's define the nodal locations of the vertices
               const Real sqrt_3 = std::sqrt(3.);
-              const auto ref_volume = target_elem.reference_elem()->volume();
-              const auto s = std::sqrt(4. / sqrt_3 * ref_volume);
+
+              std::vector<Point> equilateral_points{
+                Point(0.,  0.),
+                Point(1.,  0.),
+                Point(0.5, 0.5 * sqrt_3)
+              };
+
+              // Target element
+              const auto target_elem = Elem::build(elem->type());
+
+              // Area of the reference element
+              const auto ref_area = target_elem->reference_elem()->volume();
+
+              switch (elem->type())
+                {
+                case TRI3:
+                  {
+                    // Nothing to do here, vertices already defined in equilateral_points
+                    break;
+                  }
+
+                case TRI6:
+                  {
+                    // Define the midpoint nodes of the equilateral triangle
+                    equilateral_points.emplace_back(0.50, 0.          );
+                    equilateral_points.emplace_back(0.75, 0.25 * sqrt_3);
+                    equilateral_points.emplace_back(0.25, 0.25 * sqrt_3);
+
+                    break;
+                  }
+
+                default:
+                  libmesh_error_msg("Unsupported triangular element!");
+                  break;
+                }
+
+              // Equilateral triangle side length preserving area of the reference element
+              const auto side_length = std::sqrt(4. / sqrt_3 * ref_area);
+
+              std::vector<std::unique_ptr<Node>> owned_nodes;
 
               // Set nodes of target element to form an equilateral triangle
-              Node node_0 = Node(0., 0.);
-              Node node_1 = Node(s, 0.);
-              Node node_2 = Node(0.5 * s, 0.5 * s * sqrt_3);
-              target_elem.set_node(0, &node_0);
-              target_elem.set_node(1, &node_1);
-              target_elem.set_node(2, &node_2);
+              for (const dof_id_type node_id : index_range(equilateral_points))
+                {
+                  // Scale the nodal positions to conserve area, store the pointer to keep it alive
+                  owned_nodes.emplace_back(
+                      Node::build(equilateral_points[node_id] * side_length, node_id));
+                  target_elem->set_node(node_id, owned_nodes.back().get());
+                }
 
               // build map
-              fe_map_target.init_reference_to_physical_map(dim, qrule_points,
-                                                           &target_elem);
-              fe_map_target.compute_map(dim, qrule_weights, &target_elem,
-                                        /*d2phi=*/false);
+              fe_map_target.init_reference_to_physical_map(dim, qrule_points, target_elem.get());
+              fe_map_target.compute_map(dim, qrule_weights, target_elem.get(), /*d2phi=*/false);
 
               // Yes, triangle-to-triangle mappings have constant Jacobians, but we
               // will keep things general for now
-              _target_inverse_jacobians[target_elem.type()] =
+              _target_inverse_jacobians[target_elem->type()] =
                   std::vector<RealTensor>(nq_points);
               for (const auto qp : make_range(nq_points))
-              {
-                const RealTensor H_inv =
-                    RealTensor(dxyzdxi[qp](0), dxyzdeta[qp](0), 0, dxyzdxi[qp](1),
-                               dxyzdeta[qp](1), 0, 0, 0, 1)
-                        .inverse();
-
-                _target_inverse_jacobians[target_elem.type()][qp] = H_inv;
-                target_elem_inverse_jacobian_dets[target_elem.type()][qp] =
-                    H_inv.det();
-              }
-
-              break;
-            }
-
-          default:
-            break;
-        }
-    }
-
-    Real elem_integrated_det_J(0.);
-    for (const auto qp : index_range(JxW))
-      elem_integrated_det_J +=
-          JxW[qp] * target_elem_inverse_jacobian_dets[elem->type()][qp];
-    const auto ref_elem_vol = elem->reference_elem()->volume();
-    elem_averaged_det_J_sum += elem_integrated_det_J / ref_elem_vol;
-
-  } // for elem
+                {
+                  const RealTensor H_inv =
+                      RealTensor(dxyzdxi[qp](0), dxyzdeta[qp](0), 0, dxyzdxi[qp](1),
+                                 dxyzdeta[qp](1), 0, 0, 0, 1)
+                          .inverse();
+
+                  _target_inverse_jacobians[target_elem->type()][qp] = H_inv;
+                  target_elem_inverse_jacobian_dets[target_elem->type()][qp] =
+                      H_inv.det();
+                }
+            }// if Tri
+        }// if find == end()
+
+      // Reference volume computation
+      Real elem_integrated_det_J = 0.;
+      for (const auto qp : index_range(JxW))
+        elem_integrated_det_J +=
+            JxW[qp] * target_elem_inverse_jacobian_dets[elem->type()][qp];
+      const auto ref_elem_vol = elem->reference_elem()->volume();
+      elem_averaged_det_J_sum += elem_integrated_det_J / ref_elem_vol;
+
+    } // for elem
 
   // Get contributions from elements on other processors
   mesh.comm().sum(elem_averaged_det_J_sum);
@@ -395,7 +427,7 @@ bool VariationalSmootherSystem::element_time_derivative (bool request_jacobian,
       const Real chi_prime = 0.5 * (1. + det / sqrt_term);
       const Real chi_prime_sq = chi_prime * chi_prime;
       // d2chi(x) / dx2
-      const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / std::pow(sqrt_term, 3));
+      const Real chi_2prime = 0.5 * (1. / sqrt_term - det_sq / Utility::pow<3>(sqrt_term));
 
       // Distortion metric (beta)
       //const Real beta = trace_powers[0] / chi;
diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C
