WIP (#3385)

Author: dschwen

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_block_boundary_node_map(const MeshBase & mesh);
+
 /**
  * \returns A BoundingBox that bounds the mesh.
  */

src/mesh/mesh_smoother_laplace.C

@@ -54,10 +54,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 block_boundary_map = MeshTools::build_block_boundary_node_map(_mesh);
 
   // We can only update the nodes after all new positions were
   // determined. We store the new positions here
@@ -67,11 +64,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 +97,58 @@ 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, &block_boundary_map, &on_boundary](Node * node)
+      {
+        if (_graph[node->id()].size() == 0)
+          return;
+
+        if (on_boundary.count(node->id()))
+        {
+          // project to boundary
+          return;
+        }
+
+        // check if a node is on a given block boundary
+        auto is_node_on_block_boundary = [&block_boundary_map](dof_id_type node_id, const std::pair<subdomain_id_type, subdomain_id_type> & block_boundary)
+        {
+          auto range = block_boundary_map.equal_range(node_id);
+          for (auto i = range.first; i != range.second; ++i)
+            if (i->second == block_boundary)
+              return true;
+          return false;
+        };
+
+        const auto range = block_boundary_map.equal_range(node->id());
+        const auto num_boundaries = std::distance(range.first, range.second);
+
+        // do not touch nodes at the intersection of two or more boundaries
+        if (num_boundaries > 1)
+          return;
+
+        if (num_boundaries == 1)
+        {
+          // project to block boundary
+          const auto & boundary = range.first->second;
+          // 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))
+              continue;
+
+          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,28 @@ MeshTools::find_block_boundary_nodes(const MeshBase & mesh)
 }
 
 
-std::multimap<dof_id_type, std::pair<subdomain_id_type, subdomain_id_type>>
+std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>>
 MeshTools::build_block_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)
+      if (elem->neighbor_ptr(s) && 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);
+            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;
 }