Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // C++ includes
21 : #include <algorithm> // for std::copy, std::sort
22 :
23 : // Local includes
24 : #include "libmesh/mesh_smoother_laplace.h"
25 : #include "libmesh/mesh_tools.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/unstructured_mesh.h"
28 : #include "libmesh/parallel.h"
29 : #include "libmesh/parallel_ghost_sync.h" // sync_dofobject_data_by_id()
30 : #include "libmesh/parallel_algebra.h" // StandardType<Point>
31 : #include "libmesh/int_range.h"
32 : #include "libmesh/elem_side_builder.h"
33 :
34 : namespace libMesh
35 : {
36 : // LaplaceMeshSmoother member functions
37 936 : LaplaceMeshSmoother::LaplaceMeshSmoother(UnstructuredMesh & mesh)
38 : : MeshSmoother(mesh),
39 936 : _initialized(false)
40 : {
41 936 : }
42 :
43 :
44 :
45 :
46 1930 : void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
47 : {
48 1930 : if (!_initialized)
49 936 : this->init();
50 :
51 : // Don't smooth the nodes on the boundary...
52 : // this would change the mesh geometry which
53 : // is probably not something we want!
54 1990 : auto on_boundary = MeshTools::find_boundary_nodes(_mesh);
55 :
56 : // Also: don't smooth block boundary nodes
57 1990 : auto on_block_boundary = MeshTools::find_block_boundary_nodes(_mesh);
58 :
59 : // Merge them
60 1870 : on_boundary.insert(on_block_boundary.begin(), on_block_boundary.end());
61 :
62 : // We can only update the nodes after all new positions were
63 : // determined. We store the new positions here
64 120 : std::vector<Point> new_positions;
65 :
66 4504 : for (unsigned int n=0; n<n_iterations; n++)
67 : {
68 2574 : new_positions.resize(_mesh.max_node_id());
69 :
70 650282 : auto calculate_new_position = [this, &on_boundary, &new_positions](const Node * node) {
71 : // leave the boundary intact
72 : // Only relocate the nodes which are vertices of an element
73 : // All other entries of _graph (the secondary nodes) are empty
74 300810 : if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
75 : {
76 7674 : Point avg_position(0.,0.,0.);
77 :
78 647818 : for (const auto & connected_id : _graph[node->id()])
79 : {
80 : // Will these nodal positions always be available
81 : // or will they refer to remote nodes? This will
82 : // fail an assertion in the latter case, which
83 : // shouldn't occur if DistributedMesh is working
84 : // correctly.
85 574072 : const Point & connected_node = _mesh.point(connected_id);
86 :
87 57276 : avg_position.add( connected_node );
88 : } // end for (j)
89 :
90 : // Compute the average, store in the new_positions vector
91 104442 : new_positions[node->id()] = avg_position / static_cast<Real>(_graph[node->id()].size());
92 : } // end if
93 275412 : };
94 :
95 : // calculate new node positions (local and unpartitioned nodes only)
96 486956 : for (auto & node : _mesh.local_node_ptr_range())
97 271508 : calculate_new_position(node);
98 :
99 6980 : for (auto & node : as_range(_mesh.pid_nodes_begin(DofObject::invalid_processor_id),
100 15872 : _mesh.pid_nodes_end(DofObject::invalid_processor_id)))
101 6238 : calculate_new_position(node);
102 :
103 :
104 : // now update the node positions (local and unpartitioned nodes only)
105 486956 : for (auto & node : _mesh.local_node_ptr_range())
106 296986 : if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
107 15602 : *node = new_positions[node->id()];
108 :
109 6980 : for (auto & node : as_range(_mesh.pid_nodes_begin(DofObject::invalid_processor_id),
110 15872 : _mesh.pid_nodes_end(DofObject::invalid_processor_id)))
111 3824 : if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
112 4574 : *node = new_positions[node->id()];
113 :
114 : // Now the nodes which are ghosts on this processor may have been moved on
115 : // the processors which own them. So we need to synchronize with our neighbors
116 : // and get the most up-to-date positions for the ghosts.
117 2574 : SyncNodalPositions sync_object(_mesh);
118 : Parallel::sync_dofobject_data_by_id
119 5068 : (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
120 :
121 : } // end for n_iterations
122 :
123 : // finally adjust the second order nodes (those located between vertices)
124 : // these nodes will be located between their adjacent nodes
125 : // do this element-wise
126 778728 : for (auto & elem : _mesh.active_element_ptr_range())
127 : {
128 : // get the second order nodes (son)
129 : // their element indices start at n_vertices and go to n_nodes
130 405264 : const unsigned int son_begin = elem->n_vertices();
131 405264 : const unsigned int son_end = elem->n_nodes();
132 :
133 : // loop over all second order nodes (son)
134 1281228 : for (unsigned int son=son_begin; son<son_end; son++)
135 : {
136 : // Don't smooth second-order nodes which are on the boundary
137 921964 : if (!on_boundary.count(elem->node_id(son)))
138 : {
139 : const unsigned int n_adjacent_vertices =
140 827895 : elem->n_second_order_adjacent_vertices(son);
141 :
142 : // calculate the new position which is the average of the
143 : // position of the adjacent vertices
144 43408 : Point avg_position(0,0,0);
145 3292179 : for (unsigned int v=0; v<n_adjacent_vertices; v++)
146 : avg_position +=
147 2464284 : _mesh.point( elem->node_id( elem->second_order_adjacent_vertex(son,v) ) );
148 :
149 871303 : _mesh.node_ref(elem->node_id(son)) = avg_position / n_adjacent_vertices;
150 : }
151 : }
152 1810 : }
153 1930 : }
154 :
155 :
156 :
157 :
158 :
159 936 : void LaplaceMeshSmoother::init()
160 : {
161 : // For avoiding extraneous element side construction
162 64 : ElemSideBuilder side_builder;
163 :
164 936 : switch (_mesh.mesh_dimension())
165 : {
166 :
167 : // TODO:[BSK] Fix this to work for refined meshes... I think
168 : // the implementation was done quickly for Damien, who did not have
169 : // refined grids. Fix it here and in the original Mesh member.
170 :
171 505 : case 2: // Stolen directly from build_L_graph in mesh_base.C
172 : {
173 : // Initialize space in the graph. It is indexed by node id.
174 : // Each node may be connected to an arbitrary number of other
175 : // nodes via edges.
176 505 : _graph.resize(_mesh.max_node_id());
177 :
178 : auto elem_to_graph =
179 24584 : [this, &side_builder](const Elem & elem) {
180 34872 : for (auto s : elem.side_index_range())
181 : {
182 : // Only operate on sides which are on the
183 : // boundary or for which the current element's
184 : // id is greater than its neighbor's.
185 : // Sides get only built once.
186 34784 : if ((elem.neighbor_ptr(s) == nullptr) ||
187 7712 : (elem.id() > elem.neighbor_ptr(s)->id()))
188 : {
189 14216 : const Elem & side = side_builder(elem, s);
190 18712 : _graph[side.node_id(0)].push_back(side.node_id(1));
191 18712 : _graph[side.node_id(1)].push_back(side.node_id(0));
192 : }
193 : }
194 8607 : };
195 :
196 13088 : for (auto & elem : _mesh.active_local_element_ptr_range())
197 7085 : elem_to_graph(*elem);
198 :
199 523 : if (!_mesh.processor_id())
200 1671 : for (auto & elem : _mesh.active_unpartitioned_element_ptr_range())
201 1574 : elem_to_graph(*elem);
202 :
203 505 : _initialized = true;
204 18 : break;
205 : } // case 2
206 :
207 14 : case 3: // Stolen blatantly from build_L_graph in mesh_base.C
208 : {
209 : // Extra builder for the face elements
210 28 : ElemSideBuilder face_builder;
211 :
212 : // Initialize space in the graph.
213 431 : _graph.resize(_mesh.max_node_id());
214 :
215 : auto elem_to_graph =
216 840164 : [this, &side_builder, &face_builder](const Elem & elem) {
217 479836 : for (auto f : elem.side_index_range()) // Loop over faces
218 462032 : if ((elem.neighbor_ptr(f) == nullptr) ||
219 68160 : (elem.id() > elem.neighbor_ptr(f)->id()))
220 : {
221 201256 : const Elem & face = face_builder(elem, f);
222 :
223 876016 : for (auto s : face.side_index_range()) // Loop over face's edges
224 : {
225 674760 : const Elem & side = side_builder(face, s);
226 :
227 : // At this point, we just insert the node numbers
228 : // again. At the end we'll call sort and unique
229 : // to make sure there are no duplicates
230 796104 : _graph[side.node_id(0)].push_back(side.node_id(1));
231 796104 : _graph[side.node_id(1)].push_back(side.node_id(0));
232 : }
233 : }
234 87245 : };
235 :
236 159272 : for (auto & elem : _mesh.active_local_element_ptr_range())
237 87231 : elem_to_graph(*elem);
238 :
239 445 : if (!_mesh.processor_id())
240 143 : for (auto & elem : _mesh.active_unpartitioned_element_ptr_range())
241 61 : elem_to_graph(*elem);
242 :
243 431 : _initialized = true;
244 14 : break;
245 : } // case 3
246 :
247 0 : default:
248 0 : libmesh_error_msg("At this time it is not possible to smooth a dimension " << _mesh.mesh_dimension() << "mesh. Aborting...");
249 : }
250 :
251 : // Done building graph from local and/or unpartitioned elements.
252 : // Let's now allgather the graph so that it is available on all
253 : // processors for the actual smoothing operation.
254 936 : this->allgather_graph();
255 :
256 : // In 3D, it's possible for > 2 processor partitions to meet
257 : // at a single edge, while in 2D only 2 processor partitions
258 : // share an edge. Therefore the allgather'd graph in 3D may
259 : // now have duplicate entries and we need to remove them so
260 : // they don't foul up the averaging algorithm employed by the
261 : // Laplace smoother.
262 2446095 : for (auto & id_vec : _graph)
263 : {
264 : // The std::unique algorithm removes duplicate *consecutive* elements from a range,
265 : // so it only makes sense to call it on a sorted range...
266 2445159 : std::sort(id_vec.begin(), id_vec.end());
267 2445159 : id_vec.erase(std::unique(id_vec.begin(), id_vec.end()), id_vec.end());
268 : }
269 :
270 936 : } // init()
271 :
272 :
273 :
274 :
275 0 : void LaplaceMeshSmoother::print_graph(std::ostream & out_stream) const
276 : {
277 0 : for (auto i : index_range(_graph))
278 : {
279 0 : out_stream << i << ": ";
280 0 : std::copy(_graph[i].begin(),
281 0 : _graph[i].end(),
282 0 : std::ostream_iterator<unsigned>(out_stream, " "));
283 0 : out_stream << std::endl;
284 : }
285 0 : }
286 :
287 :
288 :
289 936 : void LaplaceMeshSmoother::allgather_graph()
290 : {
291 : // The graph data structure is not well-suited for parallel communication,
292 : // so copy the graph into a single vector defined by:
293 : // NA A_0 A_1 ... A_{NA} | NB B_0 B_1 ... B_{NB} | NC C_0 C_1 ... C_{NC}
294 : // where:
295 : // * NA is the number of graph connections for node A
296 : // * A_0, A_1, etc. are the IDs connected to node A
297 64 : std::vector<dof_id_type> flat_graph;
298 :
299 : // Reserve at least enough space for each node to have zero entries
300 968 : flat_graph.reserve(_graph.size());
301 :
302 2446095 : for (const auto & id_vec : _graph)
303 : {
304 : // First push back the number of entries for this node
305 2498133 : flat_graph.push_back (cast_int<dof_id_type>(id_vec.size()));
306 :
307 : // Then push back all the IDs
308 3823111 : for (const auto & dof : id_vec)
309 1377952 : flat_graph.push_back(dof);
310 : }
311 :
312 : // // A copy of the flat graph (for printing only, delete me later)
313 : // std::vector<unsigned> copy_of_flat_graph(flat_graph);
314 :
315 : // Use the allgather routine to combine all the flat graphs on all processors
316 936 : _mesh.comm().allgather(flat_graph);
317 :
318 : // Now reconstruct _graph from the allgathered flat_graph.
319 :
320 : // // (Delete me later, the copy is just for printing purposes.)
321 : // std::vector<std::vector<unsigned >> copy_of_graph(_graph);
322 :
323 : // Make sure the old graph is cleared out
324 904 : _graph.clear();
325 936 : const auto max_node_id = _mesh.max_node_id();
326 936 : _graph.resize(max_node_id);
327 :
328 : // Our current position in the allgather'd flat_graph
329 32 : std::size_t cursor=0;
330 :
331 : // There are max_node_id * n_processors entries to read in total
332 936 : const auto n_procs = _mesh.n_processors();
333 9970 : for (processor_id_type p = 0; p != n_procs; ++p)
334 25706843 : for (dof_id_type node_ctr : make_range(max_node_id))
335 : {
336 : // Read the number of entries for this node, move cursor
337 25697809 : std::size_t n_entries = flat_graph[cursor++];
338 :
339 : // Reserve space for that many more entries, then push back
340 25909705 : _graph[node_ctr].reserve(_graph[node_ctr].size() + n_entries);
341 :
342 : // Read all graph connections for this node, move the cursor each time
343 : // Note: there might be zero entries but that's fine
344 33681897 : for (std::size_t i=0; i<n_entries; ++i)
345 8487448 : _graph[node_ctr].push_back(flat_graph[cursor++]);
346 : }
347 :
348 : // // Print local graph to uniquely named file (debugging)
349 : // {
350 : // // Generate unique filename for this processor
351 : // std::ostringstream oss;
352 : // oss << "graph_filename_" << _mesh.processor_id() << ".txt";
353 : // std::ofstream graph_stream(oss.str().c_str());
354 : //
355 : // // Print the local non-flat graph
356 : // std::swap(_graph, copy_of_graph);
357 : // print_graph(graph_stream);
358 : //
359 : // // Print the (local) flat graph for verification
360 : // for (const auto & dof : copy_of_flat_graph)
361 : // graph_stream << dof << " ";
362 : // graph_stream << "\n";
363 : //
364 : // // Print the allgather'd grap for verification
365 : // for (const auto & dof : flat_graph)
366 : // graph_stream << dof << " ";
367 : // graph_stream << "\n";
368 : //
369 : // // Print the global non-flat graph
370 : // std::swap(_graph, copy_of_graph);
371 : // print_graph(graph_stream);
372 : // }
373 936 : } // allgather_graph()
374 :
375 : } // namespace libMesh
|
