WIP2

Author: dschwen

src/mesh/mesh_smoother_laplace.C

@@ -103,39 +103,98 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
         if (_graph[node->id()].size() == 0)
           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)
+        {
+          const auto it = block_boundary_map.find(node_id);
+          if (it == block_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);
+          const RealVectorValue vec = connected_node - *node;
+          // if any connected edge has length zero we give up and don't move the node
+          if (vec.norm_sq() < libMesh::TOLERANCE)
+            return false;
+
+          // 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);
+
+          // 3D
+          if (n_edges == 2) {
+            const auto cross = base.cross(vec);
+            if (cross.norm_sq() < libMesh::TOLERANCE)
+              n_edges--;
+            else
+              base = cross;
+            return true;
+          }
+
+          return (base * vec < libMesh::TOLERANCE);
+        };
+
         if (on_boundary.count(node->id()))
         {
           // 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
+          // base /= base.norm();
+          // auto delta = new_positions[node->id()] - *node;
+          // if (_mesh.mesh_dimension() == 2)
+          //   delta = base * (delta * base);
+          // else
+          //   delta -= base * (delta * base);
+          // *node += delta;
           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)
+        const auto it = block_boundary_map.find(node->id());
+        if (it != block_boundary_map.end())
         {
-          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 num_boundaries = it->second.size();
 
-        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;
 
-        // 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 = *(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;
 
-        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;
+            // now project
+            // if (_mesh.mesh_dimension() == 2)
 
-          return;
+            return;
+          }
         }
 
         // otherwise just move the node freely