Skip to content

Commit b55cd94

Browse files
dschwendschwen
committed
WIP (#3385)
1 parent fd2b4f7 commit b55cd94

File tree

4 files changed

+156
-28
lines changed

4 files changed

+156
-28
lines changed

‎include/mesh/mesh_smoother_laplace.h‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ class LaplaceMeshSmoother : public MeshSmoother
4949
* in the protected data section of the class.
5050
*/
5151
explicit
52-
LaplaceMeshSmoother(UnstructuredMesh & mesh);
52+
LaplaceMeshSmoother(UnstructuredMesh & mesh, bool smooth_boundary = false);
5353

5454
/**
5555
* Destructor.
@@ -98,6 +98,11 @@ class LaplaceMeshSmoother : public MeshSmoother
9898
*/
9999
bool _initialized;
100100

101+
/**
102+
* Allow boundary nodes to move tangentially
103+
*/
104+
bool _smooth_boundary;
105+
101106
/**
102107
* Data structure for holding the L-graph
103108
*/

‎include/mesh/mesh_tools.h‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -129,6 +129,13 @@ std::unordered_set<dof_id_type> find_boundary_nodes(const MeshBase & mesh);
129129
*/
130130
std::unordered_set<dof_id_type> find_block_boundary_nodes(const MeshBase & mesh);
131131

132+
/**
133+
* Returns a std::multimap for all block boundary nodes, listing all subdomain id pairs
134+
* for each block boundary the node is on
135+
*/
136+
std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>>
137+
build_subdomain_boundary_node_map(const MeshBase & mesh);
138+
132139
/**
133140
* \returns A BoundingBox that bounds the mesh.
134141
*/

‎src/mesh/mesh_smoother_laplace.C‎

Lines changed: 124 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -34,9 +34,10 @@
3434
namespace libMesh
3535
{
3636
// LaplaceMeshSmoother member functions
37-
LaplaceMeshSmoother::LaplaceMeshSmoother(UnstructuredMesh & mesh)
37+
LaplaceMeshSmoother::LaplaceMeshSmoother(UnstructuredMesh & mesh, bool smooth_boundary)
3838
: MeshSmoother(mesh),
39-
_initialized(false)
39+
_initialized(false),
40+
_smooth_boundary(smooth_boundary)
4041
{
4142
}
4243

@@ -54,10 +55,7 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
5455
auto on_boundary = MeshTools::find_boundary_nodes(_mesh);
5556

5657
// Also: don't smooth block boundary nodes
57-
auto on_block_boundary = MeshTools::find_block_boundary_nodes(_mesh);
58-
59-
// Merge them
60-
on_boundary.insert(on_block_boundary.begin(), on_block_boundary.end());
58+
auto subdomain_boundary_map = MeshTools::build_subdomain_boundary_node_map(_mesh);
6159

6260
// We can only update the nodes after all new positions were
6361
// determined. We store the new positions here
@@ -67,11 +65,11 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
6765
{
6866
new_positions.resize(_mesh.max_node_id());
6967

70-
auto calculate_new_position = [this, &on_boundary, &new_positions](const Node * node) {
68+
auto calculate_new_position = [this, &new_positions](const Node * node) {
7169
// leave the boundary intact
7270
// Only relocate the nodes which are vertices of an element
7371
// All other entries of _graph (the secondary nodes) are empty
74-
if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
72+
if (_graph[node->id()].size() > 0)
7573
{
7674
Point avg_position(0.,0.,0.);
7775

@@ -100,16 +98,130 @@ void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
10098
_mesh.pid_nodes_end(DofObject::invalid_processor_id)))
10199
calculate_new_position(node);
102100

101+
// update node position
102+
auto update_node_position = [this, &new_positions, &subdomain_boundary_map, &on_boundary](Node * node)
103+
{
104+
if (_graph[node->id()].size() == 0)
105+
return;
106+
107+
// check if a node is on a given block boundary
108+
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)
109+
{
110+
const auto it = subdomain_boundary_map.find(node_id);
111+
if (it == subdomain_boundary_map.end())
112+
return false;
113+
return it->second.count(block_boundary) > 0;
114+
};
115+
116+
// check if a series of vectors is collinear/coplanar
117+
RealVectorValue base;
118+
std::size_t n_edges = 0;
119+
auto has_zero_curvature = [this, &node, &base, &n_edges](dof_id_type connected_id) {
120+
const auto connected_node = _mesh.point(connected_id);
121+
RealVectorValue vec = connected_node - *node;
122+
const auto vec_norm_sq = vec.norm_sq();
123+
// if any connected edge has length zero we give up and don't move the node
124+
if (vec_norm_sq < libMesh::TOLERANCE * libMesh::TOLERANCE)
125+
return false;
126+
127+
// normalize
128+
vec /= std::sqrt(vec_norm_sq);
129+
130+
// count edges
131+
n_edges++;
132+
133+
// store first two edges for later projection
134+
if (n_edges == 1)
135+
{
136+
base = vec;
137+
return true;
138+
}
139+
140+
// 2D - collinear - check cross product of first edge with all other edges
141+
if (_mesh.mesh_dimension() == 2)
142+
return (base.cross(vec).norm_sq() < libMesh::TOLERANCE * libMesh::TOLERANCE);
143+
144+
// 3D
145+
if (n_edges == 2)
146+
{
147+
const auto cross = base.cross(vec);
148+
const auto cross_norm_sq = cross.norm_sq();
149+
if (cross_norm_sq < libMesh::TOLERANCE * libMesh::TOLERANCE)
150+
n_edges--;
151+
else
152+
base = cross / std::sqrt(cross_norm_sq);
153+
return true;
154+
}
155+
156+
return (base * vec < libMesh::TOLERANCE);
157+
};
158+
159+
if (on_boundary.count(node->id()))
160+
{
161+
if (!_smooth_boundary)
162+
return;
163+
164+
// project to boundary
165+
for (const auto & connected_id : _graph[node->id()])
166+
if (on_boundary.count(connected_id) && !has_zero_curvature(connected_id))
167+
return;
168+
// we have not failed the curvature check, but do we have enough edges?
169+
if (n_edges < _mesh.mesh_dimension())
170+
return;
171+
172+
// now project
173+
base /= base.norm();
174+
auto delta = new_positions[node->id()] - *node;
175+
if (_mesh.mesh_dimension() == 2)
176+
delta = base * (delta * base);
177+
else
178+
delta -= base * (delta * base);
179+
*node += delta;
180+
return;
181+
}
182+
183+
const auto it = subdomain_boundary_map.find(node->id());
184+
if (it != subdomain_boundary_map.end())
185+
{
186+
if (!_smooth_boundary)
187+
return;
188+
189+
const auto num_boundaries = it->second.size();
190+
191+
// do not touch nodes at the intersection of two or more boundaries
192+
if (num_boundaries > 1)
193+
return;
194+
195+
if (num_boundaries == 1)
196+
{
197+
// project to block boundary
198+
const auto & boundary = *(it->second.begin());
199+
// iterate over neighboring nodes that share the same boundary and check if all edges are coplanar/collinear
200+
for (const auto & connected_id : _graph[node->id()])
201+
if (is_node_on_block_boundary(connected_id, boundary) && !has_zero_curvature(connected_id))
202+
return;
203+
// we have not failed the curvature check, but do we have enough edges?
204+
if (n_edges < _mesh.mesh_dimension())
205+
return;
206+
207+
// now project
208+
// if (_mesh.mesh_dimension() == 2)
209+
210+
return;
211+
}
212+
}
213+
214+
// otherwise just move the node freely
215+
*node = new_positions[node->id()];
216+
};
103217

104218
// now update the node positions (local and unpartitioned nodes only)
105219
for (auto & node : _mesh.local_node_ptr_range())
106-
if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
107-
*node = new_positions[node->id()];
220+
update_node_position(node);
108221

109222
for (auto & node : as_range(_mesh.pid_nodes_begin(DofObject::invalid_processor_id),
110223
_mesh.pid_nodes_end(DofObject::invalid_processor_id)))
111-
if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
112-
*node = new_positions[node->id()];
224+
update_node_position(node);
113225

114226
// Now the nodes which are ghosts on this processor may have been moved on
115227
// the processors which own them. So we need to synchronize with our neighbors

‎src/mesh/mesh_tools.C‎

Lines changed: 19 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -537,7 +537,7 @@ MeshTools::find_block_boundary_nodes(const MeshBase & mesh)
537537
// mark them as true in on_boundary.
538538
for (const auto & elem : mesh.active_element_ptr_range())
539539
for (auto s : elem->side_index_range())
540-
if (elem->neighbor_ptr(s) && elem->neighbor_ptr(s)->subdomain_id() > elem->subdomain_id())
540+
if (elem->neighbor_ptr(s) && elem->neighbor_ptr(s)->subdomain_id() != elem->subdomain_id())
541541
{
542542
auto nodes_on_side = elem->nodes_on_side(s);
543543

@@ -549,27 +549,31 @@ MeshTools::find_block_boundary_nodes(const MeshBase & mesh)
549549
}
550550

551551

552-
std::multimap<dof_id_type, std::pair<subdomain_id_type, subdomain_id_type>>
553-
MeshTools::build_block_boundary_node_map(const MeshBase & mesh)
552+
std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>>
553+
MeshTools::build_subdomain_boundary_node_map(const MeshBase & mesh)
554554
{
555-
std::multimap<dof_id_type, std::pair<subdomain_id_type, subdomain_id_type>> block_boundary_node_map;
555+
std::map<dof_id_type, std::set<std::pair<subdomain_id_type, subdomain_id_type>>> block_boundary_node_map;
556556

557-
// Loop over elements, find those on boundary, and
558-
// mark them as true in on_boundary.
557+
// Loop over elements, find those on a block boundary, and
558+
// add each boundary the node is on as a subdomain id pair
559559
for (const auto & elem : mesh.active_element_ptr_range())
560560
{
561561
const auto id1 = elem->subdomain_id();
562562
for (auto s : elem->side_index_range())
563-
{
564-
const auto id2 = elem->neighbor_ptr(s)->subdomain_id();
565-
if (elem->neighbor_ptr(s) && id2 > id1)
566-
{
567-
auto nodes_on_side = elem->nodes_on_side(s);
563+
{
564+
if (elem->neighbor_ptr(s))
565+
{
566+
const auto id2 = elem->neighbor_ptr(s)->subdomain_id();
567+
if (id2 != id1)
568+
{
569+
auto nodes_on_side = elem->nodes_on_side(s);
568570

569-
for (auto & local_id : nodes_on_side)
570-
block_boundary_node_map.emplace(elem->node_ptr(local_id)->id(), id1, id2);
571-
}
572-
}
571+
for (auto & local_id : nodes_on_side)
572+
block_boundary_node_map[elem->node_ptr(local_id)->id()].insert(std::make_pair(std::min(id1, id2), std::max(id1, id2)));
573+
}
574+
}
575+
}
576+
}
573577

574578
return block_boundary_node_map;
575579
}

0 commit comments

Comments
 (0)