VariationalMeshSmoother: Support for higher order elements

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.

src/systems/variational_smoother_constraint.C

@@ -521,6 +521,51 @@ 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
+      // It doesn't make sense for a node having no edge neighbors to belong to
+      // multiple elements, right? Otherwise, we will figure out those cases as
+      // they occur.
+      libmesh_assert(nodes_to_elem_map.at(node.id()).size() == 1);
+      const auto *elem = nodes_to_elem_map.at(node.id()).front();
+      // Find the element side containing node
+      for (const auto &side : elem->side_index_range())
+        {
+          bool found_node = false;
+          const auto &nodes_on_side = elem->nodes_on_side(side);
+          for (const auto &local_node_id : nodes_on_side)
+            if (elem->node_id(local_node_id) == node.id())
+            {
+              found_node = true;
+              break;
+            }
+          if (found_node)
+            {
+              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 +633,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
@@ -672,9 +719,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

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,135 @@ 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();
-
-      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();
+      femcontext.pre_fe_reinit(*this, elem);
+      femcontext.elem_fe_reinit();
 
-      // 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 (dynamic_cast<const Tri*>(elem))
             {
-              // 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)
+              };
+
+              // Area of the reference element
+              Real ref_area;
+
+              // Target element
+              std::unique_ptr<Tri> target_elem;
+
+              switch (elem->type())
+                {
+                case TRI3:
+                  {
+                    target_elem = std::make_unique<Tri3>();
+
+                    // equilateral triangle side length that preserves area of reference
+                    // element
+                    ref_area = target_elem->reference_elem()->volume();
+
+                    break;
+                  }
+
+                case TRI6:
+                  {
+                    target_elem = std::make_unique<Tri6>();
+
+                    // equilateral triangle side length that preserves area of reference
+                    // element
+                    ref_area = target_elem->reference_elem()->volume();
+
+                    // 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
+                  auto node_ptr = std::make_unique<Node>(
+                      equilateral_points[node_id] * side_length, node_id);
+                  target_elem->set_node(node_id) = node_ptr.get();
+
+                  // Store the pointer so it stays alive
+                  owned_nodes.push_back(std::move(node_ptr));
+                }
 
               // 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);