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,117 @@ 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;
+
+ // 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;
+ }
+
+ const auto it = block_boundary_map.find(node->id());
+ if (it != block_boundary_map.end())
+ {
+ const auto num_boundaries = it->second.size();
+
+ // 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;
+
+ // now project
+ // if (_mesh.mesh_dimension() == 2)
+
+ 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;
 }