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