Move zero curvature nodes tangentially (libMesh#3385)

Author: dschwen

include/mesh/mesh_smoother_laplace.h

@@ -49,7 +49,7 @@ class LaplaceMeshSmoother : public MeshSmoother
    * in the protected data section of the class.
    */
   explicit
-  LaplaceMeshSmoother(UnstructuredMesh & mesh);
+  LaplaceMeshSmoother(UnstructuredMesh & mesh, bool smooth_boundary = false);
 
   /**
    * Destructor.
@@ -98,6 +98,11 @@ class LaplaceMeshSmoother : public MeshSmoother
    */
   bool _initialized;
 
+  /**
+   * Allow boundary nodes to move tangentially
+   */
+  bool _smooth_boundary;
+
   /**
    * Data structure for holding the L-graph
    */

include/mesh/mesh_tools.h

@@ -129,6 +129,13 @@ std::unordered_set<dof_id_type> find_boundary_nodes(const MeshBase & mesh);
  */
 std::unordered_set<dof_id_type> find_block_boundary_nodes(const MeshBase & mesh);
 
+/**
+ * Returns a std::multimap for all block boundary nodes, listing all subdomain id pairs
+ * for each block boundary the node is on
+ */
+std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>>
+build_subdomain_boundary_node_map(const MeshBase & mesh);
+
 /**
  * \returns A BoundingBox that bounds the mesh.
  */

src/mesh/mesh_smoother_laplace.C

@@ -34,9 +34,10 @@
 namespace libMesh
 {
 // LaplaceMeshSmoother member functions
-LaplaceMeshSmoother::LaplaceMeshSmoother(UnstructuredMesh & mesh)
+LaplaceMeshSmoother::LaplaceMeshSmoother(UnstructuredMesh & mesh, bool smooth_boundary)
   : MeshSmoother(mesh),
-    _initialized(false)
+    _initialized(false),
+    _smooth_boundary(smooth_boundary)
 {
 }
 
@@ -54,10 +55,7 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
   auto on_boundary = MeshTools::find_boundary_nodes(_mesh);
 
   // Also: don't smooth block boundary nodes
-  auto on_block_boundary = MeshTools::find_block_boundary_nodes(_mesh);
-
-  // Merge them
-  on_boundary.insert(on_block_boundary.begin(), on_block_boundary.end());
+  auto subdomain_boundary_map = MeshTools::build_subdomain_boundary_node_map(_mesh);
 
   // We can only update the nodes after all new positions were
   // determined. We store the new positions here
@@ -67,11 +65,11 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
     {
       new_positions.resize(_mesh.max_node_id());
 
-      auto calculate_new_position = [this, &on_boundary, &new_positions](const Node * node) {
+      auto calculate_new_position = [this, &new_positions](const Node * node) {
         // leave the boundary intact
         // Only relocate the nodes which are vertices of an element
         // All other entries of _graph (the secondary nodes) are empty
-        if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
+        if (_graph[node->id()].size() > 0)
           {
             Point avg_position(0.,0.,0.);
 
@@ -100,16 +98,134 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
                                   _mesh.pid_nodes_end(DofObject::invalid_processor_id)))
         calculate_new_position(node);
 
+      // update node position
+      auto update_node_position = [this, &new_positions, &subdomain_boundary_map, &on_boundary](Node * node)
+      {
+        if (_graph[node->id()].size() == 0)
+          return;
+
+        // check if a node is on a given block boundary
+        auto is_node_on_block_boundary = [&subdomain_boundary_map](dof_id_type node_id, const std::pair<subdomain_id_type, subdomain_id_type> & block_boundary)
+        {
+          const auto it = subdomain_boundary_map.find(node_id);
+          if (it == subdomain_boundary_map.end())
+            return false;
+          return it->second.count(block_boundary) > 0;
+        };
+
+        // check if a series of vectors is collinear/coplanar
+        RealVectorValue base;
+        std::size_t n_edges = 0;
+        auto has_zero_curvature = [this, &node, &base, &n_edges](dof_id_type connected_id) {
+          const auto connected_node = _mesh.point(connected_id);
+          RealVectorValue vec = connected_node - *node;
+          const auto vec_norm_sq = vec.norm_sq();
+          // if any connected edge has length zero we give up and don't move the node
+          if (vec_norm_sq < libMesh::TOLERANCE * libMesh::TOLERANCE)
+            return false;
+
+          // normalize
+          vec /= std::sqrt(vec_norm_sq);
+
+          // count edges
+          n_edges++;
+
+          // store first two edges for later projection
+          if (n_edges == 1)
+            {
+              base = vec;
+              return true;
+            }
+
+          // 2D - collinear - check cross product of first edge with all other edges
+          if (_mesh.mesh_dimension() == 2)
+            return (base.cross(vec).norm_sq() < libMesh::TOLERANCE * libMesh::TOLERANCE);
+
+          // 3D
+          if (n_edges == 2)
+            {
+              const auto cross = base.cross(vec);
+              const auto cross_norm_sq = cross.norm_sq();
+              if (cross_norm_sq < libMesh::TOLERANCE * libMesh::TOLERANCE)
+                n_edges--;
+              else
+                base = cross / std::sqrt(cross_norm_sq);
+              return true;
+            }
+
+          return (base * vec < libMesh::TOLERANCE);
+        };
+
+        const auto it = subdomain_boundary_map.find(node->id());
+        if (it != subdomain_boundary_map.end())
+        {
+          if (!_smooth_boundary)
+            return;
+
+          const auto num_boundaries = it->second.size();
+
+          // do not touch nodes at the intersection of two or more boundaries
+          if (num_boundaries > 1 || on_boundary.count(node->id()) > 0)
+            return;
+
+          if (num_boundaries == 1)
+            {
+              // project to block boundary
+              const auto & boundary = *(it->second.begin());
+              // iterate over neighboring nodes that share the same boundary and check if all edges are coplanar/collinear
+              for (const auto & connected_id : _graph[node->id()])
+                if (is_node_on_block_boundary(connected_id, boundary) && !has_zero_curvature(connected_id))
+                  return;
+              // we have not failed the curvature check, but do we have enough edges?
+              if (n_edges < _mesh.mesh_dimension())
+                return;
+
+              // now project
+              auto delta = new_positions[node->id()] - *node;
+              if (_mesh.mesh_dimension() == 2)
+                delta = base * (delta * base);
+              else
+                delta -= base * (delta * base);
+              *node += delta;
+              return;
+            }
+        }
+
+        if (on_boundary.count(node->id()))
+        {
+          if (!_smooth_boundary)
+            return;
+
+          // project to boundary
+          for (const auto & connected_id : _graph[node->id()])
+            if (on_boundary.count(connected_id) && !has_zero_curvature(connected_id))
+              return;
+          // we have not failed the curvature check, but do we have enough edges?
+          if (n_edges < _mesh.mesh_dimension())
+            return;
+
+          // now project
+          auto delta = new_positions[node->id()] - *node;
+          if (_mesh.mesh_dimension() == 2)
+            delta = base * (delta * base);
+          else
+            delta -= base * (delta * base);
+          *node += delta;
+          return;
+        }
+
+
+        // otherwise just move the node freely
+        *node = new_positions[node->id()];
+      };
 
       // now update the node positions (local and unpartitioned nodes only)
       for (auto & node : _mesh.local_node_ptr_range())
-        if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
-          *node = new_positions[node->id()];
+        update_node_position(node);
 
       for (auto & node : as_range(_mesh.pid_nodes_begin(DofObject::invalid_processor_id),
                                   _mesh.pid_nodes_end(DofObject::invalid_processor_id)))
-        if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
-          *node = new_positions[node->id()];
+        update_node_position(node);
 
       // Now the nodes which are ghosts on this processor may have been moved on
       // the processors which own them.  So we need to synchronize with our neighbors

src/mesh/mesh_tools.C

@@ -537,7 +537,7 @@ MeshTools::find_block_boundary_nodes(const MeshBase & mesh)
   // mark them as true in on_boundary.
   for (const auto & elem : mesh.active_element_ptr_range())
     for (auto s : elem->side_index_range())
-      if (elem->neighbor_ptr(s) && elem->neighbor_ptr(s)->subdomain_id() > elem->subdomain_id())
+      if (elem->neighbor_ptr(s) && elem->neighbor_ptr(s)->subdomain_id() != elem->subdomain_id())
         {
           auto nodes_on_side = elem->nodes_on_side(s);
 
@@ -549,27 +549,31 @@ MeshTools::find_block_boundary_nodes(const MeshBase & mesh)
 }
 
 
-std::multimap<dof_id_type, std::pair<subdomain_id_type, subdomain_id_type>>
-MeshTools::build_block_boundary_node_map(const MeshBase & mesh)
+std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>>
+MeshTools::build_subdomain_boundary_node_map(const MeshBase & mesh)
 {
-  std::multimap<dof_id_type, std::pair<subdomain_id_type, subdomain_id_type>> block_boundary_node_map;
+  std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>> block_boundary_node_map;
 
-  // Loop over elements, find those on boundary, and
-  // mark them as true in on_boundary.
+  // Loop over elements, find those on a block boundary, and
+  // add each boundary the node is on as a subdomain id pair
   for (const auto & elem : mesh.active_element_ptr_range())
   {
     const auto id1 = elem->subdomain_id();
     for (auto s : elem->side_index_range())
-    {
-      const auto id2 = elem->neighbor_ptr(s)->subdomain_id();
-      if (elem->neighbor_ptr(s) && id2 > id1)
-        {
-          auto nodes_on_side = elem->nodes_on_side(s);
+      {
+        if (elem->neighbor_ptr(s))
+          {
+            const auto id2 = elem->neighbor_ptr(s)->subdomain_id();
+            if (id2 != id1)
+              {
+                auto nodes_on_side = elem->nodes_on_side(s);
 
-          for (auto & local_id : nodes_on_side)
-            block_boundary_node_map.emplace(elem->node_ptr(local_id)->id(), id1, id2);
-        }
-    }
+                for (auto & local_id : nodes_on_side)
+                  block_boundary_node_map[elem->node_ptr(local_id)->id()].insert(std::make_pair(std::min(id1, id2), std::max(id1, id2)));
+              }
+          }
+      }
+  }
 
   return block_boundary_node_map;
 }