Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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/threads.h"
32 : #include "libmesh/elem_range.h"
33 : #include "libmesh/int_range.h"
34 : #include "libmesh/node_range.h"
35 : #include "libmesh/elem_side_builder.h"
36 : #include "libmesh/libmesh_logging.h"
37 :
38 : namespace libMesh
39 : {
40 : // LaplaceMeshSmoother member functions
41 2143 : LaplaceMeshSmoother::LaplaceMeshSmoother(UnstructuredMesh &mesh,
42 2143 : const unsigned int n_iterations)
43 2209 : : MeshSmoother(mesh), _initialized(false), _n_iterations(n_iterations) {}
44 :
45 3137 : void LaplaceMeshSmoother::smooth()
46 : {
47 282 : LOG_SCOPE("smooth()", "LaplaceMeshSmoother");
48 :
49 3137 : if (!_initialized)
50 2143 : this->init();
51 :
52 : // Don't smooth the nodes on the boundary...
53 : // this would change the mesh geometry which
54 : // is probably not something we want!
55 3231 : auto on_boundary = MeshTools::find_boundary_nodes(_mesh);
56 :
57 : // Also: don't smooth block boundary nodes
58 3231 : auto on_block_boundary = MeshTools::find_block_boundary_nodes(_mesh);
59 :
60 : // Merge them
61 3043 : on_boundary.insert(on_block_boundary.begin(), on_block_boundary.end());
62 :
63 : // We can only update the nodes after all new positions were
64 : // determined. We store the new positions here
65 94 : std::vector<Point> new_positions;
66 :
67 16077 : for (unsigned int n=0; n<_n_iterations; n++)
68 : {
69 12940 : new_positions.resize(_mesh.max_node_id());
70 :
71 1097256 : auto calculate_new_position = [this, &on_boundary, &new_positions](const Node * node) {
72 : // leave the boundary intact
73 : // Only relocate the nodes which are vertices of an element
74 : // All other entries of _graph (the secondary nodes) are empty
75 623602 : if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
76 : {
77 11208 : Point avg_position(0.,0.,0.);
78 :
79 915226 : for (const auto & connected_id : _graph[node->id()])
80 : {
81 : // Will these nodal positions always be available
82 : // or will they refer to remote nodes? This will
83 : // fail an assertion in the latter case, which
84 : // shouldn't occur if DistributedMesh is working
85 : // correctly.
86 799072 : const Point & connected_node = _mesh.point(connected_id);
87 :
88 76024 : avg_position.add( connected_node );
89 : } // end for (j)
90 :
91 : // Compute the average, store in the new_positions vector
92 160986 : new_positions[node->id()] = avg_position / static_cast<Real>(_graph[node->id()].size());
93 : } // end if
94 119092 : };
95 :
96 : // calculate new node positions (local and unpartitioned nodes only)
97 : Threads::parallel_for_each
98 12940 : (_mesh.local_node_stored_range(), calculate_new_position);
99 :
100 : Threads::parallel_for_each
101 12940 : (_mesh.unpartitioned_node_stored_range(), calculate_new_position);
102 :
103 : // now update the node positions (local and unpartitioned nodes only)
104 : Threads::parallel_for_each
105 12940 : (_mesh.local_node_stored_range(),
106 682046 : [this, &on_boundary, &new_positions](Node * node)
107 : {
108 619778 : if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
109 20256 : *node = new_positions[node->id()];
110 549726 : });
111 :
112 : Threads::parallel_for_each
113 12940 : (_mesh.unpartitioned_node_stored_range(),
114 20340 : [this, &on_boundary, &new_positions](Node * node)
115 : {
116 3824 : if (!on_boundary.count(node->id()) && (_graph[node->id()].size() > 0))
117 2160 : *node = new_positions[node->id()];
118 3824 : });
119 :
120 : // Now the nodes which are ghosts on this processor may have been moved on
121 : // the processors which own them. So we need to synchronize with our neighbors
122 : // and get the most up-to-date positions for the ghosts.
123 12940 : SyncNodalPositions sync_object(_mesh);
124 : Parallel::sync_dofobject_data_by_id
125 25508 : (_mesh.comm(), _mesh.nodes_begin(), _mesh.nodes_end(), sync_object);
126 :
127 : } // end for _n_iterations
128 :
129 : // finally adjust the second order nodes (those located between vertices)
130 : // these nodes will be located between their adjacent nodes
131 : // do this element-wise
132 : Threads::parallel_for_each
133 3137 : (_mesh.element_stored_range(),
134 2880770 : [this, &on_boundary](Elem * elem)
135 : {
136 : // get the second order nodes (son)
137 : // their element indices start at n_vertices and go to n_nodes
138 449900 : const unsigned int son_begin = elem->n_vertices();
139 449900 : const unsigned int son_end = elem->n_nodes();
140 :
141 : // loop over all second order nodes (son)
142 1975258 : for (unsigned int son=son_begin; son<son_end; son++)
143 : {
144 : // Don't smooth second-order nodes which are on the boundary
145 1613804 : if (on_boundary.count(elem->node_id(son)))
146 5304 : continue;
147 :
148 : const unsigned int n_adjacent_vertices =
149 1436618 : elem->n_second_order_adjacent_vertices(son);
150 :
151 : // calculate the new position which is the average of the
152 : // position of the adjacent vertices
153 83142 : Point avg_position(0,0,0);
154 5681930 : for (unsigned int v=0; v<n_adjacent_vertices; v++)
155 : avg_position +=
156 4245312 : _mesh.point( elem->node_id( elem->second_order_adjacent_vertex(son,v) ) );
157 :
158 1519760 : _mesh.node_ref(elem->node_id(son)) = avg_position / n_adjacent_vertices;
159 : }
160 449900 : });
161 3137 : }
162 :
163 : #ifdef LIBMESH_ENABLE_DEPRECATED
164 0 : void LaplaceMeshSmoother::smooth(unsigned int n_iterations)
165 : {
166 : libmesh_deprecated();
167 0 : _n_iterations = n_iterations;
168 0 : this->smooth();
169 0 : }
170 : #endif
171 :
172 :
173 2143 : void LaplaceMeshSmoother::init()
174 : {
175 : // For avoiding extraneous element side construction
176 132 : ElemSideBuilder side_builder;
177 :
178 2143 : switch (_mesh.mesh_dimension())
179 : {
180 :
181 : // TODO:[BSK] Fix this to work for refined meshes... I think
182 : // the implementation was done quickly for Damien, who did not have
183 : // refined grids. Fix it here and in the original Mesh member.
184 :
185 1286 : case 2: // Stolen directly from build_L_graph in mesh_base.C
186 : {
187 : // Initialize space in the graph. It is indexed by node id.
188 : // Each node may be connected to an arbitrary number of other
189 : // nodes via edges.
190 1286 : _graph.resize(_mesh.max_node_id());
191 :
192 : auto elem_to_graph =
193 40753 : [this, &side_builder](const Elem & elem) {
194 59052 : for (auto s : elem.side_index_range())
195 : {
196 : // Only operate on sides which are on the
197 : // boundary or for which the current element's
198 : // id is greater than its neighbor's.
199 : // Sides get only built once.
200 56668 : if ((elem.neighbor_ptr(s) == nullptr) ||
201 10660 : (elem.id() > elem.neighbor_ptr(s)->id()))
202 : {
203 24212 : const Elem & side = side_builder(elem, s);
204 30374 : _graph[side.node_id(0)].push_back(side.node_id(1));
205 30374 : _graph[side.node_id(1)].push_back(side.node_id(0));
206 : }
207 : }
208 14706 : };
209 :
210 24418 : for (auto & elem : _mesh.active_local_element_ptr_range())
211 13162 : elem_to_graph(*elem);
212 :
213 1326 : if (!_mesh.processor_id())
214 1924 : for (auto & elem : _mesh.active_unpartitioned_element_ptr_range())
215 1684 : elem_to_graph(*elem);
216 :
217 1286 : _initialized = true;
218 40 : break;
219 : } // case 2
220 :
221 26 : case 3: // Stolen blatantly from build_L_graph in mesh_base.C
222 : {
223 : // Extra builder for the face elements
224 52 : ElemSideBuilder face_builder;
225 :
226 : // Initialize space in the graph.
227 857 : _graph.resize(_mesh.max_node_id());
228 :
229 : auto elem_to_graph =
230 1021646 : [this, &side_builder, &face_builder](const Elem & elem) {
231 565714 : for (auto f : elem.side_index_range()) // Loop over faces
232 547652 : if ((elem.neighbor_ptr(f) == nullptr) ||
233 79920 : (elem.id() > elem.neighbor_ptr(f)->id()))
234 : {
235 239572 : const Elem & face = face_builder(elem, f);
236 :
237 1067584 : for (auto s : face.side_index_range()) // Loop over face's edges
238 : {
239 828012 : const Elem & side = side_builder(face, s);
240 :
241 : // At this point, we just insert the node numbers
242 : // again. At the end we'll call sort and unique
243 : // to make sure there are no duplicates
244 974892 : _graph[side.node_id(0)].push_back(side.node_id(1));
245 974892 : _graph[side.node_id(1)].push_back(side.node_id(0));
246 : }
247 : }
248 99929 : };
249 :
250 182608 : for (auto & elem : _mesh.active_local_element_ptr_range())
251 99903 : elem_to_graph(*elem);
252 :
253 883 : if (!_mesh.processor_id())
254 281 : for (auto & elem : _mesh.active_unpartitioned_element_ptr_range())
255 121 : elem_to_graph(*elem);
256 :
257 857 : _initialized = true;
258 26 : break;
259 : } // case 3
260 :
261 0 : default:
262 0 : libmesh_error_msg("At this time it is not possible to smooth a dimension " << _mesh.mesh_dimension() << "mesh. Aborting...");
263 : }
264 :
265 : // Done building graph from local and/or unpartitioned elements.
266 : // Let's now allgather the graph so that it is available on all
267 : // processors for the actual smoothing operation.
268 2143 : this->allgather_graph();
269 :
270 : // In 3D, it's possible for > 2 processor partitions to meet
271 : // at a single edge, while in 2D only 2 processor partitions
272 : // share an edge. Therefore the allgather'd graph in 3D may
273 : // now have duplicate entries and we need to remove them so
274 : // they don't foul up the averaging algorithm employed by the
275 : // Laplace smoother.
276 2893915 : for (auto & id_vec : _graph)
277 : {
278 : // The std::unique algorithm removes duplicate *consecutive* elements from a range,
279 : // so it only makes sense to call it on a sorted range...
280 2891772 : std::sort(id_vec.begin(), id_vec.end());
281 2891772 : id_vec.erase(std::unique(id_vec.begin(), id_vec.end()), id_vec.end());
282 : }
283 :
284 2143 : } // init()
285 :
286 :
287 :
288 :
289 0 : void LaplaceMeshSmoother::print_graph(std::ostream & out_stream) const
290 : {
291 0 : for (auto i : index_range(_graph))
292 : {
293 0 : out_stream << i << ": ";
294 0 : std::copy(_graph[i].begin(),
295 0 : _graph[i].end(),
296 0 : std::ostream_iterator<unsigned>(out_stream, " "));
297 0 : out_stream << std::endl;
298 : }
299 0 : }
300 :
301 :
302 :
303 2143 : void LaplaceMeshSmoother::allgather_graph()
304 : {
305 : // The graph data structure is not well-suited for parallel communication,
306 : // so copy the graph into a single vector defined by:
307 : // NA A_0 A_1 ... A_{NA} | NB B_0 B_1 ... B_{NB} | NC C_0 C_1 ... C_{NC}
308 : // where:
309 : // * NA is the number of graph connections for node A
310 : // * A_0, A_1, etc. are the IDs connected to node A
311 132 : std::vector<dof_id_type> flat_graph;
312 :
313 : // Reserve at least enough space for each node to have zero entries
314 2209 : flat_graph.reserve(_graph.size());
315 :
316 2893915 : for (const auto & id_vec : _graph)
317 : {
318 : // First push back the number of entries for this node
319 2964180 : flat_graph.push_back (cast_int<dof_id_type>(id_vec.size()));
320 :
321 : // Then push back all the IDs
322 4596220 : for (const auto & dof : id_vec)
323 1704448 : flat_graph.push_back(dof);
324 : }
325 :
326 : // // A copy of the flat graph (for printing only, delete me later)
327 : // std::vector<unsigned> copy_of_flat_graph(flat_graph);
328 :
329 : // Use the allgather routine to combine all the flat graphs on all processors
330 2143 : _mesh.comm().allgather(flat_graph);
331 :
332 : // Now reconstruct _graph from the allgathered flat_graph.
333 :
334 : // // (Delete me later, the copy is just for printing purposes.)
335 : // std::vector<std::vector<unsigned >> copy_of_graph(_graph);
336 :
337 : // Make sure the old graph is cleared out
338 2077 : _graph.clear();
339 2143 : const auto max_node_id = _mesh.max_node_id();
340 2143 : _graph.resize(max_node_id);
341 :
342 : // Our current position in the allgather'd flat_graph
343 66 : std::size_t cursor=0;
344 :
345 : // There are max_node_id * n_processors entries to read in total
346 2143 : const auto n_procs = _mesh.n_processors();
347 22958 : for (processor_id_type p = 0; p != n_procs; ++p)
348 29475671 : for (dof_id_type node_ctr : make_range(max_node_id))
349 : {
350 : // Read the number of entries for this node, move cursor
351 29454856 : std::size_t n_entries = flat_graph[cursor++];
352 :
353 : // Reserve space for that many more entries, then push back
354 29744488 : _graph[node_ctr].reserve(_graph[node_ctr].size() + n_entries);
355 :
356 : // Read all graph connections for this node, move the cursor each time
357 : // Note: there might be zero entries but that's fine
358 39371366 : for (std::size_t i=0; i<n_entries; ++i)
359 10528678 : _graph[node_ctr].push_back(flat_graph[cursor++]);
360 : }
361 :
362 : // // Print local graph to uniquely named file (debugging)
363 : // {
364 : // // Generate unique filename for this processor
365 : // std::ostringstream oss;
366 : // oss << "graph_filename_" << _mesh.processor_id() << ".txt";
367 : // std::ofstream graph_stream(oss.str().c_str());
368 : //
369 : // // Print the local non-flat graph
370 : // std::swap(_graph, copy_of_graph);
371 : // print_graph(graph_stream);
372 : //
373 : // // Print the (local) flat graph for verification
374 : // for (const auto & dof : copy_of_flat_graph)
375 : // graph_stream << dof << " ";
376 : // graph_stream << "\n";
377 : //
378 : // // Print the allgather'd grap for verification
379 : // for (const auto & dof : flat_graph)
380 : // graph_stream << dof << " ";
381 : // graph_stream << "\n";
382 : //
383 : // // Print the global non-flat graph
384 : // std::swap(_graph, copy_of_graph);
385 : // print_graph(graph_stream);
386 : // }
387 2143 : } // allgather_graph()
388 :
389 : } // namespace libMesh
|
