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 : // Local Includes
21 : #include "libmesh/boundary_info.h"
22 : #include "libmesh/distributed_mesh.h"
23 : #include "libmesh/elem.h"
24 : #include "libmesh/ghosting_functor.h"
25 : #include "libmesh/libmesh_config.h"
26 : #include "libmesh/libmesh_common.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_communication.h"
30 : #include "libmesh/null_output_iterator.h"
31 : #include "libmesh/mesh_tools.h"
32 : #include "libmesh/parallel.h"
33 : #include "libmesh/parallel_elem.h"
34 : #include "libmesh/parallel_node.h"
35 : #include "libmesh/parallel_ghost_sync.h"
36 : #include "libmesh/utility.h"
37 : #include "libmesh/remote_elem.h"
38 : #include "libmesh/int_range.h"
39 : #include "libmesh/elem_side_builder.h"
40 :
41 : // C++ Includes
42 : #include <numeric>
43 : #include <set>
44 : #include <unordered_set>
45 : #include <unordered_map>
46 :
47 :
48 :
49 : //-----------------------------------------------
50 : // anonymous namespace for implementation details
51 : namespace {
52 :
53 : using namespace libMesh;
54 :
55 : struct SyncNeighbors
56 : {
57 : typedef std::vector<dof_id_type> datum;
58 :
59 36 : SyncNeighbors(MeshBase & _mesh) :
60 36 : mesh(_mesh) {}
61 :
62 : MeshBase & mesh;
63 :
64 : // Find the neighbor ids for each requested element
65 20 : void gather_data (const std::vector<dof_id_type> & ids,
66 : std::vector<datum> & neighbors) const
67 : {
68 20 : neighbors.resize(ids.size());
69 :
70 92 : for (auto i : index_range(ids))
71 : {
72 : // Look for this element in the mesh
73 : // We'd better find every element we're asked for
74 72 : const Elem & elem = mesh.elem_ref(ids[i]);
75 :
76 : // Return the element's neighbors
77 38 : const unsigned int n_neigh = elem.n_neighbors();
78 72 : neighbors[i].resize(n_neigh);
79 356 : for (unsigned int n = 0; n != n_neigh; ++n)
80 : {
81 150 : const Elem * neigh = elem.neighbor_ptr(n);
82 284 : if (neigh)
83 : {
84 98 : libmesh_assert_not_equal_to(neigh, remote_elem);
85 188 : neighbors[i][n] = neigh->id();
86 : }
87 : else
88 96 : neighbors[i][n] = DofObject::invalid_id;
89 : }
90 : }
91 20 : }
92 :
93 20 : void act_on_data (const std::vector<dof_id_type> & ids,
94 : const std::vector<datum> & neighbors) const
95 : {
96 92 : for (auto i : index_range(ids))
97 : {
98 72 : Elem & elem = mesh.elem_ref(ids[i]);
99 :
100 38 : const datum & new_neigh = neighbors[i];
101 :
102 38 : const unsigned int n_neigh = elem.n_neighbors();
103 38 : libmesh_assert_equal_to (n_neigh, new_neigh.size());
104 :
105 356 : for (unsigned int n = 0; n != n_neigh; ++n)
106 : {
107 284 : const dof_id_type new_neigh_id = new_neigh[n];
108 150 : const Elem * old_neigh = elem.neighbor_ptr(n);
109 284 : if (old_neigh && old_neigh != remote_elem)
110 : {
111 98 : libmesh_assert_equal_to(old_neigh->id(), new_neigh_id);
112 : }
113 104 : else if (new_neigh_id == DofObject::invalid_id)
114 : {
115 52 : libmesh_assert (!old_neigh);
116 : }
117 : else
118 : {
119 8 : Elem * neigh = mesh.query_elem_ptr(new_neigh_id);
120 8 : if (neigh)
121 0 : elem.set_neighbor(n, neigh);
122 : else
123 8 : elem.set_neighbor(n, const_cast<RemoteElem *>(remote_elem));
124 : }
125 : }
126 : }
127 20 : }
128 : };
129 :
130 :
131 : void
132 3810 : connect_element_families(const connected_elem_set_type & connected_elements,
133 : const connected_elem_set_type & new_connected_elements,
134 : connected_elem_set_type & newer_connected_elements,
135 : const MeshBase * mesh)
136 : {
137 : // mesh was an optional parameter for API backwards compatibility
138 3810 : if (mesh && !mesh->get_constraint_rows().empty())
139 : {
140 : // We start with the constraint connections, not ancestors,
141 : // because we don't need constraining nodes of elements'
142 : // ancestors' constrained nodes.
143 0 : const auto & constraint_rows = mesh->get_constraint_rows();
144 :
145 0 : std::unordered_set<const Elem *> constraining_nodes_elems;
146 0 : for (const Elem * elem : connected_elements)
147 : {
148 0 : for (const Node & node : elem->node_ref_range())
149 : {
150 : // Retain all elements containing constraining nodes
151 0 : if (const auto it = constraint_rows.find(&node);
152 0 : it != constraint_rows.end())
153 0 : for (auto & p : it->second)
154 : {
155 0 : const Elem * constraining_elem = p.first.first;
156 0 : libmesh_assert(constraining_elem ==
157 : mesh->elem_ptr(constraining_elem->id()));
158 0 : if (!connected_elements.count(constraining_elem) &&
159 0 : !new_connected_elements.count(constraining_elem))
160 0 : constraining_nodes_elems.insert(constraining_elem);
161 : }
162 : }
163 : }
164 :
165 0 : newer_connected_elements.insert(constraining_nodes_elems.begin(),
166 : constraining_nodes_elems.end());
167 : }
168 :
169 : #ifdef LIBMESH_ENABLE_AMR
170 :
171 : // Because our set is sorted by ascending level, we can traverse it
172 : // in reverse order, adding parents as we go, and end up with all
173 : // ancestors added. This is safe for std::set where insert doesn't
174 : // invalidate iterators.
175 : //
176 : // This only works because we do *not* cache
177 : // connected_elements.rend(), whose value can change when we insert
178 : // elements which are sorted before the original rend.
179 : //
180 : // We're also going to get subactive descendents here, when any
181 : // exist. We're iterating in the wrong direction to do that
182 : // non-recursively, so we'll cop out and rely on total_family_tree.
183 : // Iterating backwards does mean that we won't be querying the newly
184 : // inserted subactive elements redundantly.
185 :
186 : connected_elem_set_type::reverse_iterator
187 1130 : elem_rit = new_connected_elements.rbegin();
188 :
189 91088 : for (; elem_rit != new_connected_elements.rend(); ++elem_rit)
190 : {
191 87278 : const Elem * elem = *elem_rit;
192 25434 : libmesh_assert(elem);
193 :
194 : // We let ghosting functors worry about only active elements,
195 : // but the remote processor needs all its semilocal elements'
196 : // ancestors and active semilocal elements' descendants too.
197 158876 : for (const Elem * parent = elem->parent(); parent;
198 59848 : parent = parent->parent())
199 27368 : if (!connected_elements.count(parent) &&
200 13684 : !new_connected_elements.count(parent))
201 22284 : newer_connected_elements.insert (parent);
202 :
203 : auto total_family_insert =
204 36410 : [&connected_elements, &new_connected_elements,
205 : &newer_connected_elements]
206 50868 : (const Elem * e)
207 : {
208 25434 : if (e->active() && e->has_children())
209 : {
210 0 : std::vector<const Elem *> subactive_family;
211 0 : e->total_family_tree(subactive_family);
212 0 : for (const auto & f : subactive_family)
213 : {
214 0 : libmesh_assert(f != remote_elem);
215 0 : if (!connected_elements.count(f) &&
216 0 : !new_connected_elements.count(f))
217 0 : newer_connected_elements.insert(f);
218 : }
219 : }
220 87278 : };
221 :
222 87278 : total_family_insert(elem);
223 :
224 : // We also need any interior parents on this mesh, which will
225 : // then need their own ancestors and descendants.
226 87278 : const Elem * interior_parent = elem->interior_parent();
227 :
228 : // Don't try to grab interior parents from other meshes, e.g. if
229 : // this was a BoundaryMesh associated with a separate Mesh.
230 :
231 : // We can't test this if someone's using the pre-mesh-ptr API
232 25434 : libmesh_assert(!interior_parent || mesh);
233 :
234 25810 : if (interior_parent &&
235 376 : interior_parent == mesh->query_elem_ptr(interior_parent->id()) &&
236 87372 : !connected_elements.count(interior_parent) &&
237 0 : !new_connected_elements.count(interior_parent))
238 : {
239 0 : newer_connected_elements.insert (interior_parent);
240 0 : total_family_insert(interior_parent);
241 : }
242 : }
243 :
244 : # ifdef DEBUG
245 : // Let's be paranoid and make sure that all our ancestors
246 : // really did get inserted. I screwed this up the first time
247 : // by caching rend, and I can easily imagine screwing it up in
248 : // the future by changing containers.
249 : auto check_elem =
250 : [&connected_elements,
251 : &new_connected_elements,
252 : &newer_connected_elements]
253 92236 : (const Elem * elem)
254 : {
255 55118 : libmesh_assert(elem);
256 55118 : const Elem * parent = elem->parent();
257 55118 : if (parent)
258 20106 : libmesh_assert(connected_elements.count(parent) ||
259 : new_connected_elements.count(parent) ||
260 : newer_connected_elements.count(parent));
261 55118 : };
262 :
263 29870 : for (const auto & elem : connected_elements)
264 28740 : check_elem(elem);
265 26564 : for (const auto & elem : new_connected_elements)
266 25434 : check_elem(elem);
267 2074 : for (const auto & elem : newer_connected_elements)
268 944 : check_elem(elem);
269 : # endif // DEBUG
270 :
271 : #endif // LIBMESH_ENABLE_AMR
272 3810 : }
273 :
274 :
275 :
276 3810 : void connect_nodes (const connected_elem_set_type & new_connected_elements,
277 : const connected_node_set_type & connected_nodes,
278 : const connected_node_set_type & new_connected_nodes,
279 : connected_node_set_type & newer_connected_nodes)
280 : {
281 91088 : for (const auto & elem : new_connected_elements)
282 976896 : for (auto & n : elem->node_ref_range())
283 1116075 : if (!connected_nodes.count(&n) &&
284 1126385 : !new_connected_nodes.count(&n))
285 849106 : newer_connected_nodes.insert(&n);
286 3810 : }
287 :
288 :
289 : } // anonymous namespace
290 :
291 :
292 :
293 : namespace libMesh
294 : {
295 :
296 :
297 3235 : void query_ghosting_functors(const MeshBase & mesh,
298 : processor_id_type pid,
299 : MeshBase::const_element_iterator elem_it,
300 : MeshBase::const_element_iterator elem_end,
301 : connected_elem_set_type & connected_elements)
302 : {
303 2424 : for (auto & gf :
304 949 : as_range(mesh.ghosting_functors_begin(),
305 8968 : mesh.ghosting_functors_end()))
306 : {
307 2174 : GhostingFunctor::map_type elements_to_ghost;
308 1087 : libmesh_assert(gf);
309 3697 : (*gf)(elem_it, elem_end, pid, elements_to_ghost);
310 :
311 : // We can ignore the CouplingMatrix in ->second, but we
312 : // need to ghost all the elements in ->first.
313 30021 : for (auto & pr : elements_to_ghost)
314 : {
315 26324 : const Elem * elem = pr.first;
316 8773 : libmesh_assert(elem != remote_elem);
317 8773 : libmesh_assert(mesh.elem_ptr(elem->id()) == elem);
318 17551 : connected_elements.insert(elem);
319 : }
320 : }
321 :
322 : // The GhostingFunctors won't be telling us about the elements from
323 : // pid; we need to add those ourselves.
324 83939 : for (; elem_it != elem_end; ++elem_it)
325 55328 : connected_elements.insert(*elem_it);
326 3235 : }
327 :
328 :
329 3235 : void connect_children(const MeshBase & mesh,
330 : MeshBase::const_element_iterator elem_it,
331 : MeshBase::const_element_iterator elem_end,
332 : connected_elem_set_type & connected_elements)
333 : {
334 : // None of these parameters are used when !LIBMESH_ENABLE_AMR.
335 949 : libmesh_ignore(mesh, elem_it, elem_end, connected_elements);
336 :
337 : #ifdef LIBMESH_ENABLE_AMR
338 : // Our XdrIO output needs inactive local elements to not have any
339 : // remote_elem children. Let's make sure that doesn't happen.
340 : //
341 91301 : for (const auto & elem : as_range(elem_it, elem_end))
342 : {
343 58892 : if (elem->has_children())
344 24076 : for (auto & child : elem->child_ref_range())
345 20512 : if (&child != remote_elem)
346 19396 : connected_elements.insert(&child);
347 1337 : }
348 : #endif // LIBMESH_ENABLE_AMR
349 3235 : }
350 :
351 :
352 0 : void reconnect_nodes (connected_elem_set_type & connected_elements,
353 : connected_node_set_type & connected_nodes)
354 : {
355 : // We're done using the nodes list for element decisions; now
356 : // let's reuse it for nodes of the elements we've decided on.
357 0 : connected_nodes.clear();
358 :
359 : // Use the newer API
360 0 : connect_nodes(connected_elements, connected_nodes, connected_nodes,
361 : connected_nodes);
362 0 : }
363 :
364 :
365 1899 : void connect_element_dependencies(const MeshBase & mesh,
366 : connected_elem_set_type & connected_elements,
367 : connected_node_set_type & connected_nodes)
368 : {
369 : // We haven't examined any of these inputs for dependencies yet, so
370 : // let's mark them all as to be examined now.
371 1126 : connected_elem_set_type new_connected_elements;
372 1126 : connected_node_set_type new_connected_nodes;
373 563 : new_connected_elements.swap(connected_elements);
374 563 : new_connected_nodes.swap(connected_nodes);
375 :
376 6742 : while (!new_connected_elements.empty() ||
377 1033 : !new_connected_nodes.empty())
378 : {
379 : auto [newer_connected_elements,
380 1130 : newer_connected_nodes] =
381 : connect_element_dependencies
382 : (mesh, connected_elements, connected_nodes,
383 6070 : new_connected_elements, new_connected_nodes);
384 :
385 : // These have now been examined
386 1130 : connected_elements.merge(new_connected_elements);
387 1130 : connected_nodes.merge(new_connected_nodes);
388 :
389 : // merge() doesn't guarantee empty() unless there are no
390 : // duplicates, which there shouldn't be
391 1130 : libmesh_assert(new_connected_elements.empty());
392 1130 : libmesh_assert(new_connected_nodes.empty());
393 :
394 : // These now need to be examined
395 1130 : new_connected_elements.swap(newer_connected_elements);
396 1130 : new_connected_nodes.swap(newer_connected_nodes);
397 1550 : }
398 1899 : }
399 :
400 :
401 : std::pair<connected_elem_set_type, connected_node_set_type>
402 3810 : connect_element_dependencies(const MeshBase & mesh,
403 : const connected_elem_set_type & connected_elements,
404 : const connected_node_set_type & connected_nodes,
405 : const connected_elem_set_type & new_connected_elements,
406 : const connected_node_set_type & new_connected_nodes)
407 : {
408 1130 : std::pair<connected_elem_set_type, connected_node_set_type> returnval;
409 1130 : auto & [newer_connected_elements, newer_connected_nodes] = returnval;
410 3810 : connect_element_families(connected_elements, new_connected_elements,
411 : newer_connected_elements, &mesh);
412 :
413 3810 : connect_nodes(new_connected_elements, connected_nodes,
414 : new_connected_nodes, newer_connected_nodes);
415 :
416 3810 : return returnval;
417 0 : }
418 :
419 :
420 :
421 :
422 : // ------------------------------------------------------------
423 : // MeshCommunication class members
424 0 : void MeshCommunication::clear ()
425 : {
426 : // _neighboring_processors.clear();
427 0 : }
428 :
429 :
430 :
431 : #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
432 : // ------------------------------------------------------------
433 : void MeshCommunication::redistribute (DistributedMesh &, bool) const
434 : {
435 : // no MPI == one processor, no redistribution
436 : return;
437 : }
438 :
439 : #else
440 : // ------------------------------------------------------------
441 414 : void MeshCommunication::redistribute (DistributedMesh & mesh,
442 : bool newly_coarsened_only) const
443 : {
444 : // This method will be called after a new partitioning has been
445 : // assigned to the elements. This partitioning was defined in
446 : // terms of the active elements, and "trickled down" to the
447 : // parents and nodes as to be consistent.
448 : //
449 : // The point is that the entire concept of local elements is
450 : // kinda shaky in this method. Elements which were previously
451 : // local may now be assigned to other processors, so we need to
452 : // send those off. Similarly, we need to accept elements from
453 : // other processors.
454 :
455 : // This method is also useful in the more limited case of
456 : // post-coarsening redistribution: if elements are only ghosting
457 : // neighbors of their active elements, but adaptive coarsening
458 : // causes an inactive element to become active, then we may need a
459 : // copy of that inactive element's neighbors.
460 :
461 : // The approach is as follows:
462 : // (1) send all relevant elements we have stored to their proper homes
463 : // (2) receive elements from all processors, watching for duplicates
464 : // (3) deleting all nonlocal elements elements
465 : // (4) obtaining required ghost elements from neighboring processors
466 138 : libmesh_parallel_only(mesh.comm());
467 138 : libmesh_assert (!mesh.is_serial());
468 138 : libmesh_assert (MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
469 : mesh.unpartitioned_elements_end()) == 0);
470 :
471 414 : LOG_SCOPE("redistribute()", "MeshCommunication");
472 :
473 : // Keep track of whether we're redistributing so that we can
474 : // relax requirements such as we must have the periodic neighbor
475 : // of an owned element
476 414 : mesh._is_redistributing = true;
477 :
478 : // Be compatible with both deprecated and corrected MeshBase iterator types
479 : typedef std::remove_const<MeshBase::const_element_iterator::value_type>::type nc_v_t;
480 :
481 : // We're going to sort elements-to-send by pid in one pass, to avoid
482 : // sending predicated iterators through the whole mesh N_p times
483 276 : std::unordered_map<processor_id_type, std::vector<nc_v_t>> send_to_pid;
484 :
485 : const MeshBase::const_element_iterator send_elems_begin =
486 : #ifdef LIBMESH_ENABLE_AMR
487 414 : newly_coarsened_only ?
488 : mesh.flagged_elements_begin(Elem::JUST_COARSENED) :
489 : #endif
490 276 : mesh.active_elements_begin();
491 :
492 : const MeshBase::const_element_iterator send_elems_end =
493 : #ifdef LIBMESH_ENABLE_AMR
494 414 : newly_coarsened_only ?
495 : mesh.flagged_elements_end(Elem::JUST_COARSENED) :
496 : #endif
497 276 : mesh.active_elements_end();
498 :
499 : // See what should get sent where. We don't send to ourselves.
500 46912 : for (auto & elem : as_range(send_elems_begin, send_elems_end))
501 46223 : if (elem->processor_id() != mesh.processor_id())
502 11627 : send_to_pid[elem->processor_id()].push_back(elem);
503 :
504 276 : std::map<processor_id_type, std::vector<const Node *>> all_nodes_to_send;
505 276 : std::map<processor_id_type, std::vector<const Elem *>> all_elems_to_send;
506 :
507 : // We may need to send constraint rows too.
508 : //
509 : // We know we're not changing n_constraint_rows(), so we're using
510 : // friend class status here rather than using the setter that
511 : // invalidates caches.
512 414 : auto & constraint_rows = mesh._constraint_rows;
513 414 : bool have_constraint_rows = !constraint_rows.empty();
514 414 : mesh.comm().broadcast(have_constraint_rows);
515 :
516 : #ifdef DEBUG
517 : const dof_id_type n_constraint_rows =
518 138 : have_constraint_rows ? mesh.n_constraint_rows() : 0;
519 : #endif
520 :
521 : typedef std::vector<std::pair<std::pair<dof_id_type, unsigned int>, Real>>
522 : serialized_row_type;
523 : std::unordered_map<processor_id_type,
524 : std::vector<std::pair<dof_id_type, serialized_row_type>>>
525 276 : all_constraint_rows_to_send;
526 :
527 : // If we don't have any just-coarsened elements to send to a
528 : // pid, then there won't be any nodes or any elements pulled
529 : // in by ghosting either, and we're done with this pid.
530 819 : for (const auto & [pid, p_elements] : send_to_pid)
531 : {
532 : // Build up a list of nodes and elements to send to processor pid.
533 : // We will certainly send all the elements assigned to this processor,
534 : // but we will also ship off any elements which are required
535 : // to be ghosted and any nodes which are used by any of the
536 : // above.
537 :
538 135 : libmesh_assert(!p_elements.empty());
539 :
540 : // Be compatible with both deprecated and
541 : // corrected MeshBase iterator types
542 : typedef MeshBase::const_element_iterator::value_type v_t;
543 :
544 405 : v_t * elempp = p_elements.data();
545 405 : v_t * elemend = elempp + p_elements.size();
546 :
547 : #ifndef LIBMESH_ENABLE_AMR
548 : // This parameter is not used when !LIBMESH_ENABLE_AMR.
549 : libmesh_ignore(newly_coarsened_only);
550 : libmesh_assert(!newly_coarsened_only);
551 : #endif
552 :
553 : MeshBase::const_element_iterator elem_it =
554 : MeshBase::const_element_iterator
555 540 : (elempp, elemend, Predicates::NotNull<v_t *>());
556 :
557 : const MeshBase::const_element_iterator elem_end =
558 : MeshBase::const_element_iterator
559 810 : (elemend, elemend, Predicates::NotNull<v_t *>());
560 :
561 270 : connected_elem_set_type elements_to_send;
562 :
563 : // See which to-be-ghosted elements we need to send
564 675 : query_ghosting_functors (mesh, pid, elem_it, elem_end,
565 : elements_to_send);
566 :
567 : // The inactive elements we need to send should have their
568 : // immediate children present.
569 945 : connect_children(mesh, mesh.pid_elements_begin(pid),
570 810 : mesh.pid_elements_end(pid),
571 : elements_to_send);
572 :
573 : // Now see which other elements and nodes they depend on
574 270 : connected_node_set_type connected_nodes;
575 405 : connect_element_dependencies(mesh, elements_to_send,
576 : connected_nodes);
577 :
578 405 : all_nodes_to_send[pid].assign(connected_nodes.begin(),
579 : connected_nodes.end());
580 :
581 405 : all_elems_to_send[pid].assign(elements_to_send.begin(),
582 : elements_to_send.end());
583 :
584 405 : for (auto & [node, row] : constraint_rows)
585 : {
586 0 : if (!connected_nodes.count(node))
587 0 : continue;
588 :
589 0 : serialized_row_type serialized_row;
590 0 : for (auto [elem_and_node, coef] : row)
591 0 : serialized_row.emplace_back(std::make_pair(elem_and_node.first->id(),
592 0 : elem_and_node.second),
593 0 : coef);
594 :
595 0 : all_constraint_rows_to_send[pid].emplace_back
596 0 : (node->id(), std::move(serialized_row));
597 : }
598 : }
599 :
600 : // Elem/Node unpack() automatically adds them to the given mesh
601 135 : auto null_node_action = [](processor_id_type, const std::vector<const Node*>&){};
602 135 : auto null_elem_action = [](processor_id_type, const std::vector<const Elem*>&){};
603 :
604 : // Communicate nodes first since elements will need to attach to them
605 414 : TIMPI::push_parallel_packed_range(mesh.comm(), all_nodes_to_send, &mesh,
606 : null_node_action);
607 :
608 414 : TIMPI::push_parallel_packed_range(mesh.comm(), all_elems_to_send, &mesh,
609 : null_elem_action);
610 :
611 : // Every element and ghost we need has now arrived
612 414 : mesh._is_redistributing = false;
613 :
614 : // At this point we have all the nodes and elems we need, so we can
615 : // communicate any constraint rows that our targets will need.
616 414 : if (have_constraint_rows)
617 : {
618 : auto constraint_row_action =
619 0 : [&mesh, &constraint_rows]
620 : (processor_id_type /* src_pid */,
621 0 : const std::vector<std::pair<dof_id_type, serialized_row_type>> rows)
622 : {
623 0 : for (auto & [node_id, serialized_row] : rows)
624 : {
625 0 : MeshBase::constraint_rows_mapped_type row;
626 0 : for (auto [elem_and_node, coef] : serialized_row)
627 0 : row.emplace_back(std::make_pair(mesh.elem_ptr(elem_and_node.first),
628 0 : elem_and_node.second),
629 0 : coef);
630 :
631 0 : constraint_rows[mesh.node_ptr(node_id)] = row;
632 : }
633 0 : };
634 :
635 0 : TIMPI::push_parallel_vector_data(mesh.comm(),
636 : all_constraint_rows_to_send,
637 : constraint_row_action);
638 :
639 : }
640 :
641 : // Check on the redistribution consistency
642 : #ifdef DEBUG
643 138 : MeshTools::libmesh_assert_equal_n_systems(mesh);
644 :
645 138 : MeshTools::libmesh_assert_valid_refinement_tree(mesh);
646 :
647 : const dof_id_type new_n_constraint_rows =
648 138 : have_constraint_rows ? mesh.n_constraint_rows() : 0;
649 :
650 138 : libmesh_assert_equal_to(n_constraint_rows, new_n_constraint_rows);
651 :
652 138 : MeshTools::libmesh_assert_valid_constraint_rows(mesh);
653 : #endif
654 :
655 : // If we had a point locator, it's invalid now that there are new
656 : // elements it can't locate.
657 414 : mesh.clear_point_locator();
658 :
659 : // Let the mesh handle any other post-redistribute() tasks, like
660 : // notifying GhostingFunctors. Be sure we're just calling the base
661 : // class method so we don't recurse back into ourselves here.
662 414 : mesh.MeshBase::redistribute();
663 414 : }
664 : #endif // LIBMESH_HAVE_MPI
665 :
666 :
667 :
668 : #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
669 : // ------------------------------------------------------------
670 : void MeshCommunication::gather_neighboring_elements (DistributedMesh &) const
671 : {
672 : // no MPI == one processor, no need for this method...
673 : return;
674 : }
675 : #else
676 : // ------------------------------------------------------------
677 41 : void MeshCommunication::gather_neighboring_elements (DistributedMesh & mesh) const
678 : {
679 : // Don't need to do anything if there is
680 : // only one processor.
681 53 : if (mesh.n_processors() == 1)
682 5 : return;
683 :
684 : // This function must be run on all processors at once
685 12 : libmesh_parallel_only(mesh.comm());
686 :
687 36 : LOG_SCOPE("gather_neighboring_elements()", "MeshCommunication");
688 :
689 : //------------------------------------------------------------------
690 : // The purpose of this function is to provide neighbor data structure
691 : // consistency for a parallel, distributed mesh. In libMesh we require
692 : // that each local element have access to a full set of valid face
693 : // neighbors. In some cases this requires us to store "ghost elements" -
694 : // elements that belong to other processors but we store to provide
695 : // data structure consistency. Also, it is assumed that any element
696 : // with a nullptr neighbor resides on a physical domain boundary. So,
697 : // even our "ghost elements" must have non-nullptr neighbors. To handle
698 : // this the concept of "RemoteElem" is used - a special construct which
699 : // is used to denote that an element has a face neighbor, but we do
700 : // not actually store detailed information about that neighbor. This
701 : // is required to prevent data structure explosion.
702 : //
703 : // So when this method is called we should have only local elements.
704 : // These local elements will then find neighbors among the local
705 : // element set. After this is completed, any element with a nullptr
706 : // neighbor has either (i) a face on the physical boundary of the mesh,
707 : // or (ii) a neighboring element which lives on a remote processor.
708 : // To handle case (ii), we communicate the global node indices connected
709 : // to all such faces to our neighboring processors. They then send us
710 : // all their elements with a nullptr neighbor that are connected to any
711 : // of the nodes in our list.
712 : //------------------------------------------------------------------
713 :
714 : // Let's begin with finding consistent neighbor data information
715 : // for all the elements we currently have. We'll use a clean
716 : // slate here - clear any existing information, including RemoteElem's.
717 60 : mesh.find_neighbors (/* reset_remote_elements = */ true,
718 24 : /* reset_current_list = */ true);
719 :
720 : // Get a unique message tag to use in communications
721 : Parallel::MessageTag
722 60 : element_neighbors_tag = mesh.comm().get_unique_tag();
723 :
724 : // Now any element with a nullptr neighbor either
725 : // (i) lives on the physical domain boundary, or
726 : // (ii) lives on an inter-processor boundary.
727 : // We will now gather all the elements from adjacent processors
728 : // which are of the same state, which should address all the type (ii)
729 : // elements.
730 :
731 : // A list of all the processors which *may* contain neighboring elements.
732 : // (for development simplicity, just make this the identity map)
733 24 : std::vector<processor_id_type> adjacent_processors;
734 108 : for (auto pid : make_range(mesh.n_processors()))
735 96 : if (pid != mesh.processor_id())
736 36 : adjacent_processors.push_back (pid);
737 :
738 :
739 : const processor_id_type n_adjacent_processors =
740 24 : cast_int<processor_id_type>(adjacent_processors.size());
741 :
742 : //-------------------------------------------------------------------------
743 : // Let's build a list of all nodes which live on nullptr-neighbor sides.
744 : // For simplicity, we will use a set to build the list, then transfer
745 : // it to a vector for communication.
746 24 : std::vector<dof_id_type> my_interface_node_list;
747 24 : std::vector<const Elem *> my_interface_elements;
748 : {
749 24 : std::set<dof_id_type> my_interface_node_set;
750 :
751 : // For avoiding extraneous element side construction
752 12 : ElemSideBuilder side_builder;
753 :
754 : // since parent nodes are a subset of children nodes, this should be sufficient
755 216 : for (const auto & elem : mesh.active_local_element_ptr_range())
756 : {
757 39 : libmesh_assert(elem);
758 :
759 156 : if (elem->on_boundary()) // denotes *any* side has a nullptr neighbor
760 : {
761 117 : my_interface_elements.push_back(elem); // add the element, but only once, even
762 : // if there are multiple nullptr neighbors
763 618 : for (auto s : elem->side_index_range())
764 616 : if (elem->neighbor_ptr(s) == nullptr)
765 : {
766 286 : const Elem & side = side_builder(*elem, s);
767 :
768 858 : for (auto n : make_range(side.n_vertices()))
769 768 : my_interface_node_set.insert (side.node_id(n));
770 : }
771 : }
772 12 : }
773 :
774 36 : my_interface_node_list.reserve (my_interface_node_set.size());
775 24 : my_interface_node_list.insert (my_interface_node_list.end(),
776 : my_interface_node_set.begin(),
777 36 : my_interface_node_set.end());
778 : }
779 :
780 : // we will now send my_interface_node_list to all of the adjacent processors.
781 : // note that for the time being we will copy the list to a unique buffer for
782 : // each processor so that we can use a nonblocking send and not access the
783 : // buffer again until the send completes. it is my understanding that the
784 : // MPI 2.1 standard seeks to remove this restriction as unnecessary, so in
785 : // the future we should change this to send the same buffer to each of the
786 : // adjacent processors. - BSK 11/17/2008
787 : std::vector<std::vector<dof_id_type>>
788 60 : my_interface_node_xfer_buffers (n_adjacent_processors, my_interface_node_list);
789 24 : std::map<processor_id_type, unsigned char> n_comm_steps;
790 :
791 60 : std::vector<Parallel::Request> send_requests (3*n_adjacent_processors);
792 12 : unsigned int current_request = 0;
793 :
794 72 : for (unsigned int comm_step=0; comm_step<n_adjacent_processors; comm_step++)
795 : {
796 48 : n_comm_steps[adjacent_processors[comm_step]]=1;
797 60 : mesh.comm().send (adjacent_processors[comm_step],
798 24 : my_interface_node_xfer_buffers[comm_step],
799 36 : send_requests[current_request++],
800 : element_neighbors_tag);
801 : }
802 :
803 : //-------------------------------------------------------------------------
804 : // processor pairings are symmetric - I expect to receive an interface node
805 : // list from each processor in adjacent_processors as well!
806 : // now we will catch an incoming node list for each of our adjacent processors.
807 : //
808 : // we are done with the adjacent_processors list - note that it is in general
809 : // a superset of the processors we truly share elements with. so let's
810 : // clear the superset list, and we will fill it with the true list.
811 12 : adjacent_processors.clear();
812 :
813 12 : std::vector<dof_id_type> common_interface_node_list;
814 :
815 : // we expect two classes of messages -
816 : // (1) incoming interface node lists, to which we will reply with our elements
817 : // touching nodes in the list, and
818 : // (2) replies from the requests we sent off previously.
819 : // (2.a) - nodes
820 : // (2.b) - elements
821 : // so we expect 3 communications from each adjacent processor.
822 : // by structuring the communication in this way we hopefully impose no
823 : // order on the handling of the arriving messages. in particular, we
824 : // should be able to handle the case where we receive a request and
825 : // all replies from processor A before even receiving a request from
826 : // processor B.
827 :
828 144 : for (unsigned int comm_step=0; comm_step<3*n_adjacent_processors; comm_step++)
829 : {
830 : //------------------------------------------------------------------
831 : // catch incoming node list
832 : Parallel::Status
833 108 : status(mesh.comm().probe (Parallel::any_source,
834 72 : element_neighbors_tag));
835 : const processor_id_type
836 108 : source_pid_idx = cast_int<processor_id_type>(status.source()),
837 36 : dest_pid_idx = source_pid_idx;
838 :
839 : //------------------------------------------------------------------
840 : // first time - incoming request
841 108 : if (n_comm_steps[source_pid_idx] == 1)
842 : {
843 36 : n_comm_steps[source_pid_idx]++;
844 :
845 36 : mesh.comm().receive (source_pid_idx,
846 : common_interface_node_list,
847 24 : element_neighbors_tag);
848 :
849 : // const std::size_t
850 : // their_interface_node_list_size = common_interface_node_list.size();
851 :
852 : // we now have the interface node list from processor source_pid_idx.
853 : // now we can find all of our elements which touch any of these nodes
854 : // and send copies back to this processor. however, we can make our
855 : // search more efficient by first excluding all the nodes in
856 : // their list which are not also contained in
857 : // my_interface_node_list. we can do this in place as a set
858 : // intersection.
859 : common_interface_node_list.erase
860 12 : (std::set_intersection (my_interface_node_list.begin(),
861 : my_interface_node_list.end(),
862 : common_interface_node_list.begin(),
863 : common_interface_node_list.end(),
864 12 : common_interface_node_list.begin()),
865 24 : common_interface_node_list.end());
866 :
867 : // if (false)
868 : // libMesh::out << "[" << mesh.processor_id() << "] "
869 : // << "my_interface_node_list.size()=" << my_interface_node_list.size()
870 : // << ", [" << source_pid_idx << "] "
871 : // << "their_interface_node_list.size()=" << their_interface_node_list_size
872 : // << ", common_interface_node_list.size()=" << common_interface_node_list.size()
873 : // << std::endl;
874 :
875 : // Now we need to see which of our elements touch the nodes in the list.
876 : // We will certainly send all the active elements which intersect source_pid_idx,
877 : // but we will also ship off the other elements in the same family tree
878 : // as the active ones for data structure consistency.
879 : //
880 : // FIXME - shipping full family trees is unnecessary and inefficient.
881 : //
882 : // We also ship any nodes connected to these elements. Note
883 : // some of these nodes and elements may be replicated from
884 : // other processors, but that is OK.
885 12 : connected_elem_set_type elements_to_send;
886 :
887 : // Technically we're not doing reconnect_nodes on this, but
888 : // we might as well use the same data structure since the
889 : // performance questions ought to be similar
890 12 : connected_node_set_type connected_nodes;
891 :
892 : // Check for quick return?
893 38 : if (common_interface_node_list.empty())
894 : {
895 : // let's try to be smart here - if we have no nodes in common,
896 : // we cannot share elements. so post the messages expected
897 : // from us here and go on about our business.
898 : // note that even though these are nonblocking sends
899 : // they should complete essentially instantly, because
900 : // in all cases the send buffers are empty
901 8 : mesh.comm().send_packed_range (dest_pid_idx,
902 : &mesh,
903 : connected_nodes.begin(),
904 : connected_nodes.end(),
905 6 : send_requests[current_request++],
906 : element_neighbors_tag);
907 :
908 8 : mesh.comm().send_packed_range (dest_pid_idx,
909 : &mesh,
910 : elements_to_send.begin(),
911 : elements_to_send.end(),
912 6 : send_requests[current_request++],
913 : element_neighbors_tag);
914 :
915 2 : continue;
916 : }
917 : // otherwise, this really *is* an adjacent processor.
918 30 : adjacent_processors.push_back(source_pid_idx);
919 :
920 20 : std::vector<const Elem *> family_tree;
921 :
922 144 : for (auto & elem : my_interface_elements)
923 : {
924 38 : std::size_t n_shared_nodes = 0;
925 :
926 190 : for (auto n : make_range(elem->n_vertices()))
927 186 : if (std::binary_search (common_interface_node_list.begin(),
928 : common_interface_node_list.end(),
929 302 : elem->node_id(n)))
930 : {
931 38 : n_shared_nodes++;
932 :
933 : // TBD - how many nodes do we need to share
934 : // before we care? certainly 2, but 1? not
935 : // sure, so let's play it safe...
936 38 : if (n_shared_nodes > 0) break;
937 : }
938 :
939 114 : if (n_shared_nodes) // share at least one node?
940 : {
941 182 : elem = elem->top_parent();
942 :
943 : // avoid a lot of duplicated effort -- if we already have elem
944 : // in the set its entire family tree is already in the set.
945 38 : if (!elements_to_send.count(elem))
946 : {
947 : #ifdef LIBMESH_ENABLE_AMR
948 110 : elem->family_tree(family_tree);
949 : #else
950 : family_tree.clear();
951 : family_tree.push_back(elem);
952 : #endif
953 220 : for (const auto & f : family_tree)
954 : {
955 110 : elem = f;
956 72 : elements_to_send.insert (elem);
957 :
958 600 : for (auto & n : elem->node_ref_range())
959 452 : connected_nodes.insert (&n);
960 : }
961 : }
962 : }
963 : }
964 :
965 : // The elements_to_send and connected_nodes sets now contain all
966 : // the elements and nodes we need to send to this processor.
967 : // All that remains is to pack up the objects (along with
968 : // any boundary conditions) and send the messages off.
969 : {
970 10 : libmesh_assert (connected_nodes.empty() || !elements_to_send.empty());
971 10 : libmesh_assert (!connected_nodes.empty() || elements_to_send.empty());
972 :
973 : // send the nodes off to the destination processor
974 40 : mesh.comm().send_packed_range (dest_pid_idx,
975 : &mesh,
976 : connected_nodes.begin(),
977 : connected_nodes.end(),
978 30 : send_requests[current_request++],
979 : element_neighbors_tag);
980 :
981 : // send the elements off to the destination processor
982 40 : mesh.comm().send_packed_range (dest_pid_idx,
983 : &mesh,
984 : elements_to_send.begin(),
985 : elements_to_send.end(),
986 30 : send_requests[current_request++],
987 : element_neighbors_tag);
988 : }
989 : }
990 : //------------------------------------------------------------------
991 : // second time - reply of nodes
992 72 : else if (n_comm_steps[source_pid_idx] == 2)
993 : {
994 36 : n_comm_steps[source_pid_idx]++;
995 :
996 36 : mesh.comm().receive_packed_range (source_pid_idx,
997 : &mesh,
998 : null_output_iterator<Node>(),
999 : (Node**)nullptr,
1000 : element_neighbors_tag);
1001 : }
1002 : //------------------------------------------------------------------
1003 : // third time - reply of elements
1004 36 : else if (n_comm_steps[source_pid_idx] == 3)
1005 : {
1006 36 : n_comm_steps[source_pid_idx]++;
1007 :
1008 36 : mesh.comm().receive_packed_range (source_pid_idx,
1009 : &mesh,
1010 : null_output_iterator<Elem>(),
1011 : (Elem**)nullptr,
1012 : element_neighbors_tag);
1013 : }
1014 : //------------------------------------------------------------------
1015 : // fourth time - shouldn't happen
1016 : else
1017 : {
1018 0 : libMesh::err << "ERROR: unexpected number of replies: "
1019 0 : << n_comm_steps[source_pid_idx]
1020 0 : << std::endl;
1021 : }
1022 : } // done catching & processing replies associated with tag ~ 100,000pi
1023 :
1024 : // allow any pending requests to complete
1025 36 : Parallel::wait (send_requests);
1026 :
1027 : // If we had a point locator, it's invalid now that there are new
1028 : // elements it can't locate.
1029 36 : mesh.clear_point_locator();
1030 :
1031 : // We can now find neighbor information for the interfaces between
1032 : // local elements and ghost elements.
1033 48 : mesh.find_neighbors (/* reset_remote_elements = */ true,
1034 24 : /* reset_current_list = */ false);
1035 :
1036 : // Ghost elements may not have correct remote_elem neighbor links,
1037 : // and we may not be able to locally infer correct neighbor links to
1038 : // remote elements. So we synchronize ghost element neighbor links.
1039 12 : SyncNeighbors nsync(mesh);
1040 :
1041 : Parallel::sync_dofobject_data_by_id
1042 60 : (mesh.comm(), mesh.elements_begin(), mesh.elements_end(), nsync);
1043 36 : }
1044 : #endif // LIBMESH_HAVE_MPI
1045 :
1046 :
1047 : #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
1048 : // ------------------------------------------------------------
1049 : void MeshCommunication::send_coarse_ghosts(MeshBase &) const
1050 : {
1051 : // no MPI == one processor, no need for this method...
1052 : return;
1053 : }
1054 : #else
1055 885 : void MeshCommunication::send_coarse_ghosts(MeshBase & mesh) const
1056 : {
1057 :
1058 : // Don't need to do anything if all processors already ghost all non-local
1059 : // elements.
1060 885 : if (mesh.is_serial())
1061 873 : return;
1062 :
1063 : // This algorithm uses the MeshBase::flagged_pid_elements_begin/end iterators
1064 : // which are only available when AMR is enabled.
1065 : #ifndef LIBMESH_ENABLE_AMR
1066 : libmesh_error_msg("Calling MeshCommunication::send_coarse_ghosts() requires AMR to be enabled. "
1067 : "Please configure libmesh with --enable-amr.");
1068 : #else
1069 : // When we coarsen elements on a DistributedMesh, we make their
1070 : // parents active. This may increase the ghosting requirements on
1071 : // the processor which owns the newly-activated parent element. To
1072 : // ensure ghosting requirements are satisfied, processors which
1073 : // coarsen an element will send all the associated ghosted elements
1074 : // to all processors which own any of the coarsened-away-element's
1075 : // siblings.
1076 : typedef std::unordered_map<processor_id_type, std::vector<Elem *>> ghost_map;
1077 8 : ghost_map coarsening_elements_to_ghost;
1078 :
1079 8 : const processor_id_type proc_id = mesh.processor_id();
1080 : // Look for just-coarsened elements
1081 740 : for (auto elem : as_range(mesh.flagged_pid_elements_begin(Elem::COARSEN, proc_id),
1082 3640 : mesh.flagged_pid_elements_end(Elem::COARSEN, proc_id)))
1083 : {
1084 : // If it's flagged for coarsening it had better have a parent
1085 720 : libmesh_assert(elem->parent());
1086 :
1087 : // On a distributed mesh:
1088 : // If we don't own this element's parent but we do own it, then
1089 : // there is a chance that we are aware of ghost elements which
1090 : // the parent's owner needs us to send them.
1091 2160 : const processor_id_type their_proc_id = elem->parent()->processor_id();
1092 2160 : if (their_proc_id != proc_id)
1093 0 : coarsening_elements_to_ghost[their_proc_id].push_back(elem);
1094 4 : }
1095 :
1096 8 : std::map<processor_id_type, std::vector<const Node *>> all_nodes_to_send;
1097 8 : std::map<processor_id_type, std::vector<const Elem *>> all_elems_to_send;
1098 :
1099 8 : const processor_id_type n_proc = mesh.n_processors();
1100 :
1101 36 : for (processor_id_type p=0; p != n_proc; ++p)
1102 : {
1103 24 : if (p == proc_id)
1104 12 : continue;
1105 :
1106 8 : connected_elem_set_type elements_to_send;
1107 8 : std::set<const Node *> nodes_to_send;
1108 :
1109 12 : if (const auto it = std::as_const(coarsening_elements_to_ghost).find(p);
1110 4 : it != coarsening_elements_to_ghost.end())
1111 : {
1112 0 : const std::vector<Elem *> & elems = it->second;
1113 0 : libmesh_assert(elems.size());
1114 :
1115 : // Make some fake element iterators defining this vector of
1116 : // elements
1117 0 : Elem * const * elempp = const_cast<Elem * const *>(elems.data());
1118 0 : Elem * const * elemend = elempp+elems.size();
1119 : const MeshBase::const_element_iterator elem_it =
1120 0 : MeshBase::const_element_iterator(elempp, elemend, Predicates::NotNull<Elem * const *>());
1121 : const MeshBase::const_element_iterator elem_end =
1122 0 : MeshBase::const_element_iterator(elemend, elemend, Predicates::NotNull<Elem * const *>());
1123 :
1124 0 : for (auto & gf : as_range(mesh.ghosting_functors_begin(),
1125 0 : mesh.ghosting_functors_end()))
1126 : {
1127 0 : GhostingFunctor::map_type elements_to_ghost;
1128 0 : libmesh_assert(gf);
1129 0 : (*gf)(elem_it, elem_end, p, elements_to_ghost);
1130 :
1131 : // We can ignore the CouplingMatrix in ->second, but we
1132 : // need to ghost all the elements in ->first.
1133 0 : for (auto & pr : elements_to_ghost)
1134 : {
1135 0 : const Elem * elem = pr.first;
1136 0 : libmesh_assert(elem);
1137 0 : while (elem)
1138 : {
1139 0 : libmesh_assert(elem != remote_elem);
1140 0 : elements_to_send.insert(elem);
1141 0 : for (auto & n : elem->node_ref_range())
1142 0 : nodes_to_send.insert(&n);
1143 0 : elem = elem->parent();
1144 : }
1145 : }
1146 : }
1147 :
1148 0 : all_nodes_to_send[p].assign(nodes_to_send.begin(), nodes_to_send.end());
1149 0 : all_elems_to_send[p].assign(elements_to_send.begin(), elements_to_send.end());
1150 : }
1151 : }
1152 :
1153 : // Elem/Node unpack() automatically adds them to the given mesh
1154 0 : auto null_node_action = [](processor_id_type, const std::vector<const Node*>&){};
1155 0 : auto null_elem_action = [](processor_id_type, const std::vector<const Elem*>&){};
1156 :
1157 : // Communicate nodes first since elements will need to attach to them
1158 12 : TIMPI::push_parallel_packed_range(mesh.comm(), all_nodes_to_send, &mesh,
1159 : null_node_action);
1160 :
1161 12 : TIMPI::push_parallel_packed_range(mesh.comm(), all_elems_to_send, &mesh,
1162 : null_elem_action);
1163 :
1164 : #endif // LIBMESH_ENABLE_AMR
1165 : }
1166 :
1167 : #endif // LIBMESH_HAVE_MPI
1168 :
1169 : #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
1170 : // ------------------------------------------------------------
1171 : void MeshCommunication::broadcast (MeshBase &) const
1172 : {
1173 : // no MPI == one processor, no need for this method...
1174 : return;
1175 : }
1176 : #else
1177 : // ------------------------------------------------------------
1178 1576 : void MeshCommunication::broadcast (MeshBase & mesh) const
1179 : {
1180 : // Don't need to do anything if there is
1181 : // only one processor.
1182 2024 : if (mesh.n_processors() == 1)
1183 228 : return;
1184 :
1185 : // This function must be run on all processors at once
1186 444 : libmesh_parallel_only(mesh.comm());
1187 :
1188 1332 : LOG_SCOPE("broadcast()", "MeshCommunication");
1189 :
1190 : // Explicitly clear the mesh on all but processor 0.
1191 1792 : if (mesh.processor_id() != 0)
1192 674 : mesh.clear();
1193 :
1194 : // We may have set extra data only on processor 0 in a read()
1195 1348 : mesh.comm().broadcast(mesh._node_integer_names);
1196 1348 : mesh.comm().broadcast(mesh._node_integer_default_values);
1197 1348 : mesh.comm().broadcast(mesh._elem_integer_names);
1198 1348 : mesh.comm().broadcast(mesh._elem_integer_default_values);
1199 :
1200 : // We may have set mapping data only on processor 0 in a read()
1201 1348 : unsigned char map_type = mesh.default_mapping_type();
1202 1348 : unsigned char map_data = mesh.default_mapping_data();
1203 1348 : mesh.comm().broadcast(map_type);
1204 1348 : mesh.comm().broadcast(map_data);
1205 1348 : mesh.set_default_mapping_type(ElemMappingType(map_type));
1206 1348 : mesh.set_default_mapping_data(map_data);
1207 :
1208 : // Broadcast nodes
1209 1792 : mesh.comm().broadcast_packed_range(&mesh,
1210 2696 : mesh.nodes_begin(),
1211 1792 : mesh.nodes_end(),
1212 : &mesh,
1213 : null_output_iterator<Node>());
1214 :
1215 : // Broadcast elements from coarsest to finest, so that child
1216 : // elements will see their parents already in place.
1217 : //
1218 : // When restarting from a checkpoint, we may have elements which are
1219 : // assigned to a processor but which have not yet been sent to that
1220 : // processor, so we need to use a paranoid n_levels() count and not
1221 : // the usual fast algorithm.
1222 1348 : const unsigned int n_levels = MeshTools::paranoid_n_levels(mesh);
1223 :
1224 2696 : for (unsigned int l=0; l != n_levels; ++l)
1225 1792 : mesh.comm().broadcast_packed_range(&mesh,
1226 2696 : mesh.level_elements_begin(l),
1227 2696 : mesh.level_elements_end(l),
1228 : &mesh,
1229 : null_output_iterator<Elem>());
1230 :
1231 : // Make sure mesh_dimension and elem_dimensions are consistent.
1232 1348 : mesh.cache_elem_data();
1233 :
1234 : // Make sure mesh id counts are consistent.
1235 1348 : mesh.update_parallel_id_counts();
1236 :
1237 : // We may have constraint rows on IsoGeometric Analysis meshes. We
1238 : // don't want to send these along with constrained nodes (like we
1239 : // send boundary info for those nodes) because the associated rows'
1240 : // elements may not exist at that point.
1241 : //
1242 : // We know we're not changing n_constraint_rows(), so we're using
1243 : // friend class status here rather than using the setter that
1244 : // invalidates caches.
1245 904 : auto & constraint_rows = mesh._constraint_rows;
1246 :
1247 1348 : bool have_constraint_rows = !constraint_rows.empty();
1248 1348 : mesh.comm().broadcast(have_constraint_rows);
1249 1348 : if (have_constraint_rows)
1250 : {
1251 : std::map<dof_id_type,
1252 : std::vector<std::tuple<dof_id_type, unsigned int, Real>>>
1253 40 : serialized_rows;
1254 :
1255 18494 : for (auto & row : constraint_rows)
1256 : {
1257 18416 : const Node * node = row.first;
1258 18416 : const dof_id_type rowid = node->id();
1259 5850 : libmesh_assert(node == mesh.node_ptr(rowid));
1260 :
1261 : std::vector<std::tuple<dof_id_type, unsigned int, Real>>
1262 11700 : serialized_row;
1263 56918 : for (auto & entry : row.second)
1264 : serialized_row.push_back
1265 38502 : (std::make_tuple(entry.first.first->id(),
1266 12215 : entry.first.second, entry.second));
1267 :
1268 12566 : serialized_rows.emplace(rowid, std::move(serialized_row));
1269 : }
1270 :
1271 78 : mesh.comm().broadcast(serialized_rows);
1272 98 : if (mesh.processor_id() != 0)
1273 : {
1274 10 : constraint_rows.clear();
1275 :
1276 18455 : for (auto & row : serialized_rows)
1277 : {
1278 18416 : const dof_id_type rowid = row.first;
1279 18416 : const Node * node = mesh.node_ptr(rowid);
1280 :
1281 : std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>>
1282 11700 : deserialized_row;
1283 56918 : for (auto & entry : row.second)
1284 : deserialized_row.push_back
1285 77004 : (std::make_pair(std::make_pair(mesh.elem_ptr(std::get<0>(entry)),
1286 36645 : std::get<1>(entry)),
1287 36645 : std::get<2>(entry)));
1288 :
1289 12566 : constraint_rows.emplace(node, deserialized_row);
1290 : }
1291 : }
1292 : }
1293 :
1294 : // Broadcast all of the named entity information
1295 1348 : mesh.comm().broadcast(mesh.set_subdomain_name_map());
1296 1348 : mesh.comm().broadcast(mesh.get_boundary_info().set_sideset_name_map());
1297 1348 : mesh.comm().broadcast(mesh.get_boundary_info().set_nodeset_name_map());
1298 :
1299 : // If we had a point locator, it's invalid now that there are new
1300 : // elements it can't locate.
1301 1348 : mesh.clear_point_locator();
1302 :
1303 444 : libmesh_assert (mesh.comm().verify(mesh.n_elem()));
1304 444 : libmesh_assert (mesh.comm().verify(mesh.max_elem_id()));
1305 444 : libmesh_assert (mesh.comm().verify(mesh.n_nodes()));
1306 444 : libmesh_assert (mesh.comm().verify(mesh.max_node_id()));
1307 :
1308 : #ifdef DEBUG
1309 444 : MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
1310 444 : MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1311 : #endif
1312 : }
1313 : #endif // LIBMESH_HAVE_MPI
1314 :
1315 :
1316 :
1317 : #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
1318 : // ------------------------------------------------------------
1319 : void MeshCommunication::gather (const processor_id_type, MeshBase &) const
1320 : {
1321 : // no MPI == one processor, no need for this method...
1322 : return;
1323 : }
1324 : #else
1325 : // ------------------------------------------------------------
1326 253 : void MeshCommunication::gather (const processor_id_type root_id, MeshBase & mesh) const
1327 : {
1328 : // Check for quick return
1329 329 : if (mesh.n_processors() == 1)
1330 25 : return;
1331 :
1332 : // This function must be run on all processors at once
1333 76 : libmesh_parallel_only(mesh.comm());
1334 :
1335 228 : LOG_SCOPE("(all)gather()", "MeshCommunication");
1336 :
1337 : // Ensure we don't build too big a buffer at once
1338 : static const std::size_t approx_total_buffer_size = 1e8;
1339 : const std::size_t approx_each_buffer_size =
1340 228 : approx_total_buffer_size / mesh.comm().size();
1341 :
1342 228 : (root_id == DofObject::invalid_processor_id) ?
1343 :
1344 248 : mesh.comm().allgather_packed_range (&mesh,
1345 248 : mesh.nodes_begin(),
1346 414 : mesh.nodes_end(),
1347 : null_output_iterator<Node>(),
1348 : approx_each_buffer_size) :
1349 :
1350 56 : mesh.comm().gather_packed_range (root_id,
1351 : &mesh,
1352 270 : mesh.nodes_begin(),
1353 298 : mesh.nodes_end(),
1354 : null_output_iterator<Node>(),
1355 : approx_each_buffer_size);
1356 :
1357 : // Gather elements from coarsest to finest, so that child
1358 : // elements will see their parents already in place.
1359 228 : const unsigned int n_levels = MeshTools::n_levels(mesh);
1360 :
1361 504 : for (unsigned int l=0; l != n_levels; ++l)
1362 276 : (root_id == DofObject::invalid_processor_id) ?
1363 :
1364 296 : mesh.comm().allgather_packed_range (&mesh,
1365 296 : mesh.level_elements_begin(l),
1366 646 : mesh.level_elements_end(l),
1367 : null_output_iterator<Elem>(),
1368 : approx_each_buffer_size) :
1369 :
1370 72 : mesh.comm().gather_packed_range (root_id,
1371 : &mesh,
1372 330 : mesh.level_elements_begin(l),
1373 366 : mesh.level_elements_end(l),
1374 : null_output_iterator<Elem>(),
1375 : approx_each_buffer_size);
1376 :
1377 : // If we had a point locator, it's invalid now that there are new
1378 : // elements it can't locate.
1379 228 : mesh.clear_point_locator();
1380 :
1381 : // We may have constraint rows on IsoGeometric Analysis meshes. We
1382 : // don't want to send these along with constrained nodes (like we
1383 : // send boundary info for those nodes) because the associated rows'
1384 : // elements may not exist at that point.
1385 : //
1386 : // We know we're not changing n_constraint_rows(), so we're using
1387 : // friend class status here rather than using the setter that
1388 : // invalidates caches.
1389 152 : auto & constraint_rows = mesh._constraint_rows;
1390 228 : bool have_constraint_rows = !constraint_rows.empty();
1391 228 : mesh.comm().max(have_constraint_rows);
1392 228 : if (have_constraint_rows)
1393 : {
1394 : std::map<dof_id_type,
1395 : std::vector<std::tuple<dof_id_type, unsigned int, Real>>>
1396 0 : serialized_rows;
1397 :
1398 0 : for (auto & row : constraint_rows)
1399 : {
1400 0 : const Node * node = row.first;
1401 0 : const dof_id_type rowid = node->id();
1402 0 : libmesh_assert(node == mesh.node_ptr(rowid));
1403 :
1404 : std::vector<std::tuple<dof_id_type, unsigned int, Real>>
1405 0 : serialized_row;
1406 0 : for (auto & entry : row.second)
1407 : serialized_row.push_back
1408 0 : (std::make_tuple(entry.first.first->id(),
1409 0 : entry.first.second, entry.second));
1410 :
1411 0 : serialized_rows.emplace(rowid, std::move(serialized_row));
1412 : }
1413 :
1414 0 : if (root_id == DofObject::invalid_processor_id)
1415 0 : mesh.comm().set_union(serialized_rows);
1416 : else
1417 0 : mesh.comm().set_union(serialized_rows, root_id);
1418 :
1419 0 : if (root_id == DofObject::invalid_processor_id ||
1420 0 : root_id == mesh.processor_id())
1421 : {
1422 0 : for (auto & row : serialized_rows)
1423 : {
1424 0 : const dof_id_type rowid = row.first;
1425 0 : const Node * node = mesh.node_ptr(rowid);
1426 :
1427 : std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>>
1428 0 : deserialized_row;
1429 0 : for (auto & entry : row.second)
1430 : deserialized_row.push_back
1431 0 : (std::make_pair(std::make_pair(mesh.elem_ptr(std::get<0>(entry)),
1432 0 : std::get<1>(entry)),
1433 0 : std::get<2>(entry)));
1434 :
1435 0 : constraint_rows.emplace(node, deserialized_row);
1436 : }
1437 : }
1438 : #ifdef DEBUG
1439 0 : MeshTools::libmesh_assert_valid_constraint_rows(mesh);
1440 : #endif
1441 : }
1442 :
1443 :
1444 : // If we are doing an allgather(), perform sanity check on the result.
1445 76 : if (root_id == DofObject::invalid_processor_id)
1446 : {
1447 62 : libmesh_assert (mesh.comm().verify(mesh.n_elem()));
1448 62 : libmesh_assert (mesh.comm().verify(mesh.n_nodes()));
1449 : }
1450 :
1451 : // Inform new elements of their neighbors,
1452 : // while resetting all remote_elem links on
1453 : // the ranks which did the gather.
1454 228 : if (mesh.allow_find_neighbors())
1455 304 : mesh.find_neighbors(root_id == DofObject::invalid_processor_id ||
1456 152 : root_id == mesh.processor_id());
1457 :
1458 : // All done, but let's make sure it's done correctly
1459 :
1460 : #ifdef DEBUG
1461 76 : MeshTools::libmesh_assert_valid_boundary_ids(mesh);
1462 : #endif
1463 : }
1464 : #endif // LIBMESH_HAVE_MPI
1465 :
1466 :
1467 :
1468 : // Functor for make_elems_parallel_consistent and
1469 : // make_node_ids_parallel_consistent
1470 : namespace {
1471 :
1472 : struct SyncIds
1473 : {
1474 : typedef dof_id_type datum;
1475 : typedef void (MeshBase::*renumber_obj)(dof_id_type, dof_id_type);
1476 :
1477 129 : SyncIds(MeshBase & _mesh, renumber_obj _renumberer) :
1478 53 : mesh(_mesh),
1479 129 : renumber(_renumberer) {}
1480 :
1481 : MeshBase & mesh;
1482 : renumber_obj renumber;
1483 : // renumber_obj & renumber;
1484 :
1485 : // Find the id of each requested DofObject -
1486 : // Parallel::sync_* already did the work for us
1487 34 : void gather_data (const std::vector<dof_id_type> & ids,
1488 : std::vector<datum> & ids_out) const
1489 : {
1490 102 : ids_out = ids;
1491 68 : }
1492 :
1493 102 : void act_on_data (const std::vector<dof_id_type> & old_ids,
1494 : const std::vector<datum> & new_ids) const
1495 : {
1496 3798 : for (auto i : index_range(old_ids))
1497 6160 : if (old_ids[i] != new_ids[i])
1498 3480 : (mesh.*renumber)(old_ids[i], new_ids[i]);
1499 102 : }
1500 : };
1501 :
1502 :
1503 66 : struct SyncNodeIds
1504 : {
1505 : typedef dof_id_type datum;
1506 :
1507 162 : SyncNodeIds(MeshBase & _mesh) :
1508 162 : mesh(_mesh) {}
1509 :
1510 : MeshBase & mesh;
1511 :
1512 : // We only know a Node id() is definitive if we own the Node or if
1513 : // we're told it's definitive. We keep track of the latter cases by
1514 : // putting definitively id'd ghost nodes into this set.
1515 : typedef std::unordered_set<const Node *> uset_type;
1516 : uset_type definitive_ids;
1517 :
1518 : // We should never be told two different definitive ids for the same
1519 : // node, but let's check on that in debug mode.
1520 : #ifdef DEBUG
1521 : typedef std::unordered_map<dof_id_type, dof_id_type> umap_type;
1522 : umap_type definitive_renumbering;
1523 : #endif
1524 :
1525 : // Find the id of each requested DofObject -
1526 : // Parallel::sync_* already tried to do the work for us, but we can
1527 : // only say the result is definitive if we own the DofObject or if
1528 : // we were given the definitive result from another processor.
1529 396 : void gather_data (const std::vector<dof_id_type> & ids,
1530 : std::vector<datum> & ids_out) const
1531 : {
1532 132 : ids_out.clear();
1533 528 : ids_out.resize(ids.size(), DofObject::invalid_id);
1534 :
1535 102168 : for (auto i : index_range(ids))
1536 : {
1537 101772 : const dof_id_type id = ids[i];
1538 101772 : const Node * node = mesh.query_node_ptr(id);
1539 101772 : if (node && (node->processor_id() == mesh.processor_id() ||
1540 33924 : definitive_ids.count(node)))
1541 135696 : ids_out[i] = id;
1542 : }
1543 396 : }
1544 :
1545 396 : bool act_on_data (const std::vector<dof_id_type> & old_ids,
1546 : const std::vector<datum> & new_ids)
1547 : {
1548 132 : bool data_changed = false;
1549 102168 : for (auto i : index_range(old_ids))
1550 : {
1551 101772 : const dof_id_type new_id = new_ids[i];
1552 :
1553 101772 : const dof_id_type old_id = old_ids[i];
1554 :
1555 101772 : Node * node = mesh.query_node_ptr(old_id);
1556 :
1557 : // If we can't find the node we were asking about, another
1558 : // processor must have already given us the definitive id
1559 : // for it
1560 101772 : if (!node)
1561 : {
1562 : // But let's check anyway in debug mode
1563 : #ifdef DEBUG
1564 5626 : libmesh_assert
1565 : (definitive_renumbering.count(old_id));
1566 5626 : libmesh_assert_equal_to
1567 : (new_id, definitive_renumbering[old_id]);
1568 : #endif
1569 11252 : continue;
1570 : }
1571 :
1572 : // If we asked for an id but there's no definitive id ready
1573 : // for us yet, then we can't quit trying to sync yet.
1574 84894 : if (new_id == DofObject::invalid_id)
1575 : {
1576 : // But we might have gotten a definitive id from a
1577 : // different request
1578 0 : if (!definitive_ids.count(mesh.node_ptr(old_id)))
1579 0 : data_changed = true;
1580 : }
1581 : else
1582 : {
1583 113192 : if (node->processor_id() != mesh.processor_id())
1584 46812 : definitive_ids.insert(node);
1585 84894 : if (old_id != new_id)
1586 : {
1587 : #ifdef DEBUG
1588 3234 : libmesh_assert
1589 : (!definitive_renumbering.count(old_id));
1590 3234 : definitive_renumbering[old_id] = new_id;
1591 : #endif
1592 9702 : mesh.renumber_node(old_id, new_id);
1593 3234 : data_changed = true;
1594 : }
1595 : }
1596 : }
1597 396 : return data_changed;
1598 : }
1599 : };
1600 :
1601 :
1602 : #ifdef LIBMESH_ENABLE_AMR
1603 : struct SyncPLevels
1604 : {
1605 : typedef std::pair<unsigned char,unsigned char> datum;
1606 :
1607 0 : SyncPLevels(MeshBase & _mesh) :
1608 0 : mesh(_mesh) {}
1609 :
1610 : MeshBase & mesh;
1611 :
1612 : // Find the p_level of each requested Elem
1613 0 : void gather_data (const std::vector<dof_id_type> & ids,
1614 : std::vector<datum> & ids_out) const
1615 : {
1616 0 : ids_out.reserve(ids.size());
1617 :
1618 0 : for (const auto & id : ids)
1619 : {
1620 0 : Elem & elem = mesh.elem_ref(id);
1621 : ids_out.push_back
1622 0 : (std::make_pair(cast_int<unsigned char>(elem.p_level()),
1623 0 : static_cast<unsigned char>(elem.p_refinement_flag())));
1624 : }
1625 0 : }
1626 :
1627 0 : void act_on_data (const std::vector<dof_id_type> & old_ids,
1628 : const std::vector<datum> & new_p_levels) const
1629 : {
1630 0 : for (auto i : index_range(old_ids))
1631 : {
1632 0 : Elem & elem = mesh.elem_ref(old_ids[i]);
1633 : // Make sure these are consistent
1634 : elem.hack_p_level_and_refinement_flag
1635 0 : (new_p_levels[i].first,
1636 0 : static_cast<Elem::RefinementState>(new_p_levels[i].second));
1637 : // Make sure parents' levels are consistent
1638 0 : elem.set_p_level(new_p_levels[i].first);
1639 : }
1640 0 : }
1641 : };
1642 : #endif // LIBMESH_ENABLE_AMR
1643 :
1644 :
1645 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1646 : template <typename DofObjSubclass>
1647 : struct SyncUniqueIds
1648 : {
1649 : typedef unique_id_type datum;
1650 : typedef DofObjSubclass* (MeshBase::*query_obj)(const dof_id_type);
1651 :
1652 422 : SyncUniqueIds(MeshBase &_mesh, query_obj _querier) :
1653 174 : mesh(_mesh),
1654 422 : query(_querier) {}
1655 :
1656 : MeshBase & mesh;
1657 : query_obj query;
1658 :
1659 : // Find the id of each requested DofObject -
1660 : // Parallel::sync_* already did the work for us
1661 288 : void gather_data (const std::vector<dof_id_type> & ids,
1662 : std::vector<datum> & ids_out) const
1663 : {
1664 384 : ids_out.reserve(ids.size());
1665 :
1666 20420 : for (const auto & id : ids)
1667 : {
1668 20132 : DofObjSubclass * d = (mesh.*query)(id);
1669 6712 : libmesh_assert(d);
1670 20132 : ids_out.push_back(d->unique_id());
1671 : }
1672 288 : }
1673 :
1674 288 : void act_on_data (const std::vector<dof_id_type> & ids,
1675 : const std::vector<datum> & unique_ids) const
1676 : {
1677 20420 : for (auto i : index_range(ids))
1678 : {
1679 20132 : DofObjSubclass * d = (mesh.*query)(ids[i]);
1680 6712 : libmesh_assert(d);
1681 26844 : d->set_unique_id(unique_ids[i]);
1682 : }
1683 288 : }
1684 : };
1685 : #endif // LIBMESH_ENABLE_UNIQUE_ID
1686 :
1687 : template <typename DofObjSubclass>
1688 : struct SyncBCIds
1689 : {
1690 : typedef std::vector<boundary_id_type> datum;
1691 :
1692 129 : SyncBCIds(MeshBase &_mesh) :
1693 129 : mesh(_mesh) {}
1694 :
1695 : MeshBase & mesh;
1696 :
1697 : // Find the id of each requested DofObject -
1698 : // Parallel::sync_* already did the work for us
1699 68 : void gather_data (const std::vector<dof_id_type> & ids,
1700 : std::vector<datum> & ids_out) const
1701 : {
1702 68 : ids_out.reserve(ids.size());
1703 :
1704 68 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1705 :
1706 6030 : for (const auto & id : ids)
1707 : {
1708 5962 : Node * n = mesh.query_node_ptr(id);
1709 2981 : libmesh_assert(n);
1710 5962 : std::vector<boundary_id_type> bcids;
1711 5962 : boundary_info.boundary_ids(n, bcids);
1712 2981 : ids_out.push_back(std::move(bcids));
1713 : }
1714 68 : }
1715 :
1716 68 : void act_on_data (const std::vector<dof_id_type> & ids,
1717 : const std::vector<datum> & bcids) const
1718 : {
1719 68 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1720 :
1721 6030 : for (auto i : index_range(ids))
1722 : {
1723 5962 : Node * n = mesh.query_node_ptr(ids[i]);
1724 2981 : libmesh_assert(n);
1725 5962 : boundary_info.add_node(n, bcids[i]);
1726 : }
1727 68 : }
1728 : };
1729 :
1730 : }
1731 :
1732 :
1733 :
1734 : // ------------------------------------------------------------
1735 162 : void MeshCommunication::make_node_ids_parallel_consistent (MeshBase & mesh)
1736 : {
1737 : // This function must be run on all processors at once
1738 48 : libmesh_parallel_only(mesh.comm());
1739 :
1740 : // We need to agree on which processor owns every node, but we can't
1741 : // easily assert that here because we don't currently agree on which
1742 : // id every node has, and some of our temporary ids on unrelated
1743 : // nodes will "overlap".
1744 : //#ifdef DEBUG
1745 : // MeshTools::libmesh_assert_parallel_consistent_procids<Node> (mesh);
1746 : //#endif // DEBUG
1747 :
1748 144 : LOG_SCOPE ("make_node_ids_parallel_consistent()", "MeshCommunication");
1749 :
1750 96 : SyncNodeIds syncids(mesh);
1751 : Parallel::sync_node_data_by_element_id
1752 390 : (mesh, mesh.elements_begin(), mesh.elements_end(),
1753 114 : Parallel::SyncEverything(), Parallel::SyncEverything(), syncids);
1754 :
1755 : // At this point, with both ids and processor ids synced, we can
1756 : // finally check for topological consistency of node processor ids.
1757 : #ifdef DEBUG
1758 48 : MeshTools::libmesh_assert_topology_consistent_procids<Node> (mesh);
1759 : #endif
1760 162 : }
1761 :
1762 :
1763 :
1764 293 : void MeshCommunication::make_node_unique_ids_parallel_consistent (MeshBase & mesh)
1765 : {
1766 : // Avoid unused variable warnings if unique ids aren't enabled.
1767 86 : libmesh_ignore(mesh);
1768 :
1769 : // This function must be run on all processors at once
1770 86 : libmesh_parallel_only(mesh.comm());
1771 :
1772 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1773 172 : LOG_SCOPE ("make_node_unique_ids_parallel_consistent()", "MeshCommunication");
1774 :
1775 86 : SyncUniqueIds<Node> syncuniqueids(mesh, &MeshBase::query_node_ptr);
1776 172 : Parallel::sync_dofobject_data_by_id(mesh.comm(),
1777 586 : mesh.nodes_begin(),
1778 465 : mesh.nodes_end(),
1779 : syncuniqueids);
1780 :
1781 : #endif
1782 293 : }
1783 :
1784 :
1785 129 : void MeshCommunication::make_node_bcids_parallel_consistent (MeshBase & mesh)
1786 : {
1787 : // Avoid unused variable warnings if unique ids aren't enabled.
1788 38 : libmesh_ignore(mesh);
1789 :
1790 : // This function must be run on all processors at once
1791 38 : libmesh_parallel_only(mesh.comm());
1792 :
1793 76 : LOG_SCOPE ("make_node_bcids_parallel_consistent()", "MeshCommunication");
1794 :
1795 38 : SyncBCIds<Node> syncbcids(mesh);
1796 76 : Parallel::sync_dofobject_data_by_id(mesh.comm(),
1797 258 : mesh.nodes_begin(),
1798 205 : mesh.nodes_end(),
1799 : syncbcids);
1800 129 : }
1801 :
1802 :
1803 :
1804 :
1805 :
1806 : // ------------------------------------------------------------
1807 129 : void MeshCommunication::make_elems_parallel_consistent(MeshBase & mesh)
1808 : {
1809 : // This function must be run on all processors at once
1810 38 : libmesh_parallel_only(mesh.comm());
1811 :
1812 76 : LOG_SCOPE ("make_elems_parallel_consistent()", "MeshCommunication");
1813 :
1814 38 : SyncIds syncids(mesh, &MeshBase::renumber_elem);
1815 : Parallel::sync_element_data_by_parent_id
1816 220 : (mesh, mesh.active_elements_begin(),
1817 167 : mesh.active_elements_end(), syncids);
1818 :
1819 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1820 38 : SyncUniqueIds<Elem> syncuniqueids(mesh, &MeshBase::query_elem_ptr);
1821 : Parallel::sync_dofobject_data_by_id
1822 220 : (mesh.comm(), mesh.active_elements_begin(),
1823 205 : mesh.active_elements_end(), syncuniqueids);
1824 : #endif
1825 129 : }
1826 :
1827 :
1828 :
1829 : // ------------------------------------------------------------
1830 : #ifdef LIBMESH_ENABLE_AMR
1831 0 : void MeshCommunication::make_p_levels_parallel_consistent(MeshBase & mesh)
1832 : {
1833 : // This function must be run on all processors at once
1834 0 : libmesh_parallel_only(mesh.comm());
1835 :
1836 0 : LOG_SCOPE ("make_p_levels_parallel_consistent()", "MeshCommunication");
1837 :
1838 0 : SyncPLevels syncplevels(mesh);
1839 : Parallel::sync_dofobject_data_by_id
1840 0 : (mesh.comm(), mesh.elements_begin(), mesh.elements_end(),
1841 : syncplevels);
1842 0 : }
1843 : #endif // LIBMESH_ENABLE_AMR
1844 :
1845 :
1846 :
1847 : // Functors for make_node_proc_ids_parallel_consistent
1848 : namespace {
1849 :
1850 : struct SyncProcIds
1851 : {
1852 : typedef processor_id_type datum;
1853 :
1854 162 : SyncProcIds(MeshBase & _mesh) : mesh(_mesh) {}
1855 :
1856 : MeshBase & mesh;
1857 :
1858 : // ------------------------------------------------------------
1859 172 : void gather_data (const std::vector<dof_id_type> & ids,
1860 : std::vector<datum> & data)
1861 : {
1862 : // Find the processor id of each requested node
1863 172 : data.resize(ids.size());
1864 :
1865 40120 : for (auto i : index_range(ids))
1866 : {
1867 : // Look for this point in the mesh
1868 39948 : if (ids[i] != DofObject::invalid_id)
1869 : {
1870 39948 : Node & node = mesh.node_ref(ids[i]);
1871 :
1872 : // Return the node's correct processor id,
1873 39948 : data[i] = node.processor_id();
1874 : }
1875 : else
1876 0 : data[i] = DofObject::invalid_processor_id;
1877 : }
1878 172 : }
1879 :
1880 : // ------------------------------------------------------------
1881 172 : bool act_on_data (const std::vector<dof_id_type> & ids,
1882 : const std::vector<datum> proc_ids)
1883 : {
1884 86 : bool data_changed = false;
1885 :
1886 : // Set the ghost node processor ids we've now been informed of
1887 40120 : for (auto i : index_range(ids))
1888 : {
1889 39948 : Node & node = mesh.node_ref(ids[i]);
1890 :
1891 : // We may not have ids synched when this synchronization is done, so we
1892 : // *can't* use id to load-balance processor id properly; we have to use
1893 : // the old heuristic of choosing the smallest valid processor id.
1894 : //
1895 : // If someone tells us our node processor id is too low, then
1896 : // they're wrong. If they tell us our node processor id is
1897 : // too high, then we're wrong.
1898 39948 : if (node.processor_id() > proc_ids[i])
1899 : {
1900 0 : data_changed = true;
1901 0 : node.processor_id() = proc_ids[i];
1902 : }
1903 : }
1904 :
1905 172 : return data_changed;
1906 : }
1907 : };
1908 :
1909 :
1910 : struct ElemNodesMaybeNew
1911 : {
1912 38 : ElemNodesMaybeNew() {}
1913 :
1914 25096 : bool operator() (const Elem * elem) const
1915 : {
1916 : // If this element was just refined then it may have new nodes we
1917 : // need to work on
1918 : #ifdef LIBMESH_ENABLE_AMR
1919 25096 : if (elem->refinement_flag() == Elem::JUST_REFINED)
1920 7896 : return true;
1921 : #endif
1922 :
1923 : // If this element has remote_elem neighbors then there may have
1924 : // been refinement of those neighbors that affect its nodes'
1925 : // processor_id()
1926 33108 : for (auto neigh : elem->neighbor_ptr_range())
1927 27516 : if (neigh == remote_elem)
1928 1160 : return true;
1929 5592 : return false;
1930 : }
1931 : };
1932 :
1933 :
1934 53 : struct NodeWasNew
1935 : {
1936 129 : NodeWasNew(const MeshBase & mesh)
1937 129 : {
1938 52678 : for (const auto & node : mesh.node_ptr_range())
1939 38575 : if (node->processor_id() == DofObject::invalid_processor_id)
1940 10531 : was_new.insert(node);
1941 129 : }
1942 :
1943 46948 : bool operator() (const Elem * elem, unsigned int local_node_num) const
1944 : {
1945 46948 : if (was_new.count(elem->node_ptr(local_node_num)))
1946 27744 : return true;
1947 9602 : return false;
1948 : }
1949 :
1950 : std::unordered_set<const Node *> was_new;
1951 : };
1952 :
1953 : }
1954 :
1955 :
1956 :
1957 : // ------------------------------------------------------------
1958 33 : void MeshCommunication::make_node_proc_ids_parallel_consistent(MeshBase & mesh)
1959 : {
1960 10 : LOG_SCOPE ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");
1961 :
1962 : // This function must be run on all processors at once
1963 10 : libmesh_parallel_only(mesh.comm());
1964 :
1965 : // When this function is called, each section of a parallelized mesh
1966 : // should be in the following state:
1967 : //
1968 : // All nodes should have the exact same physical location on every
1969 : // processor where they exist.
1970 : //
1971 : // Local nodes should have unique authoritative ids,
1972 : // and processor ids consistent with all processors which own
1973 : // an element touching them.
1974 : //
1975 : // Ghost nodes touching local elements should have processor ids
1976 : // consistent with all processors which own an element touching
1977 : // them.
1978 10 : SyncProcIds sync(mesh);
1979 : Parallel::sync_node_data_by_element_id
1980 79 : (mesh, mesh.elements_begin(), mesh.elements_end(),
1981 33 : Parallel::SyncEverything(), Parallel::SyncEverything(), sync);
1982 33 : }
1983 :
1984 :
1985 :
1986 : // ------------------------------------------------------------
1987 129 : void MeshCommunication::make_new_node_proc_ids_parallel_consistent(MeshBase & mesh)
1988 : {
1989 114 : LOG_SCOPE ("make_new_node_proc_ids_parallel_consistent()", "MeshCommunication");
1990 :
1991 : // This function must be run on all processors at once
1992 38 : libmesh_parallel_only(mesh.comm());
1993 :
1994 : // When this function is called, each section of a parallelized mesh
1995 : // should be in the following state:
1996 : //
1997 : // Local nodes should have unique authoritative ids,
1998 : // and new nodes should be unpartitioned.
1999 : //
2000 : // New ghost nodes touching local elements should be unpartitioned.
2001 :
2002 : // We may not have consistent processor ids for new nodes (because a
2003 : // node may be old and partitioned on one processor but new and
2004 : // unpartitioned on another) when we start
2005 : #ifdef DEBUG
2006 38 : MeshTools::libmesh_assert_parallel_consistent_procids<Node>(mesh);
2007 : // MeshTools::libmesh_assert_parallel_consistent_new_node_procids(mesh);
2008 : #endif
2009 :
2010 : // We have two kinds of new nodes. *NEW* nodes are unpartitioned on
2011 : // all processors: we need to use a id-independent (i.e. dumb)
2012 : // heuristic to partition them. But "new" nodes are newly created
2013 : // on some processors (when ghost elements are refined) yet
2014 : // correspond to existing nodes on other processors: we need to use
2015 : // the existing processor id for them.
2016 : //
2017 : // A node which is "new" on one processor will be associated with at
2018 : // least one ghost element, and we can just query that ghost
2019 : // element's owner to find out the correct processor id.
2020 :
2021 : auto node_unpartitioned =
2022 24034 : [](const Elem * elem, unsigned int local_node_num)
2023 24034 : { return elem->node_ref(local_node_num).processor_id() ==
2024 24034 : DofObject::invalid_processor_id; };
2025 :
2026 38 : SyncProcIds sync(mesh);
2027 :
2028 : sync_node_data_by_element_id_once
2029 311 : (mesh, mesh.not_local_elements_begin(),
2030 167 : mesh.not_local_elements_end(), Parallel::SyncEverything(),
2031 : node_unpartitioned, sync);
2032 :
2033 : // Nodes should now be unpartitioned iff they are truly new; those
2034 : // are the *only* nodes we will touch.
2035 : #ifdef DEBUG
2036 38 : MeshTools::libmesh_assert_parallel_consistent_new_node_procids(mesh);
2037 : #endif
2038 :
2039 167 : NodeWasNew node_was_new(mesh);
2040 :
2041 : // Set the lowest processor id we can on truly new nodes
2042 25316 : for (auto & elem : mesh.element_ptr_range())
2043 167835 : for (auto & node : elem->node_ref_range())
2044 144251 : if (node_was_new.was_new.count(&node))
2045 : {
2046 18292 : processor_id_type & pid = node.processor_id();
2047 74469 : pid = std::min(pid, elem->processor_id());
2048 53 : }
2049 :
2050 : // Then finally see if other processors have a lower option
2051 : Parallel::sync_node_data_by_element_id
2052 311 : (mesh, mesh.elements_begin(), mesh.elements_end(),
2053 91 : ElemNodesMaybeNew(), node_was_new, sync);
2054 :
2055 : // We should have consistent processor ids when we're done.
2056 : #ifdef DEBUG
2057 38 : MeshTools::libmesh_assert_parallel_consistent_procids<Node>(mesh);
2058 38 : MeshTools::libmesh_assert_parallel_consistent_new_node_procids(mesh);
2059 : #endif
2060 129 : }
2061 :
2062 :
2063 :
2064 : // ------------------------------------------------------------
2065 33 : void MeshCommunication::make_nodes_parallel_consistent (MeshBase & mesh)
2066 : {
2067 : // This function must be run on all processors at once
2068 10 : libmesh_parallel_only(mesh.comm());
2069 :
2070 : // When this function is called, each section of a parallelized mesh
2071 : // should be in the following state:
2072 : //
2073 : // Element ids and locations should be consistent on every processor
2074 : // where they exist.
2075 : //
2076 : // All nodes should have the exact same physical location on every
2077 : // processor where they exist.
2078 : //
2079 : // Local nodes should have unique authoritative ids,
2080 : // and processor ids consistent with all processors which own
2081 : // an element touching them.
2082 : //
2083 : // Ghost nodes touching local elements should have processor ids
2084 : // consistent with all processors which own an element touching
2085 : // them.
2086 : //
2087 : // Ghost nodes should have ids which are either already correct
2088 : // or which are in the "unpartitioned" id space.
2089 :
2090 : // First, let's sync up processor ids. Some of these processor ids
2091 : // may be "wrong" from coarsening, but they're right in the sense
2092 : // that they'll tell us who has the authoritative dofobject ids for
2093 : // each node.
2094 :
2095 33 : this->make_node_proc_ids_parallel_consistent(mesh);
2096 :
2097 : // Second, sync up dofobject ids.
2098 33 : this->make_node_ids_parallel_consistent(mesh);
2099 :
2100 : // Third, sync up dofobject unique_ids if applicable.
2101 33 : this->make_node_unique_ids_parallel_consistent(mesh);
2102 :
2103 : // Finally, correct the processor ids to make DofMap happy
2104 33 : MeshTools::correct_node_proc_ids(mesh);
2105 33 : }
2106 :
2107 :
2108 :
2109 : // ------------------------------------------------------------
2110 129 : void MeshCommunication::make_new_nodes_parallel_consistent (MeshBase & mesh)
2111 : {
2112 : // This function must be run on all processors at once
2113 38 : libmesh_parallel_only(mesh.comm());
2114 :
2115 : // When this function is called, each section of a parallelized mesh
2116 : // should be in the following state:
2117 : //
2118 : // All nodes should have the exact same physical location on every
2119 : // processor where they exist.
2120 : //
2121 : // Local nodes should have unique authoritative ids,
2122 : // and new nodes should be unpartitioned.
2123 : //
2124 : // New ghost nodes touching local elements should be unpartitioned.
2125 : //
2126 : // New ghost nodes should have ids which are either already correct
2127 : // or which are in the "unpartitioned" id space.
2128 : //
2129 : // Non-new nodes should have correct ids and processor ids already.
2130 :
2131 : // First, let's sync up new nodes' processor ids.
2132 :
2133 129 : this->make_new_node_proc_ids_parallel_consistent(mesh);
2134 :
2135 : // Second, sync up dofobject ids.
2136 129 : this->make_node_ids_parallel_consistent(mesh);
2137 :
2138 : // Third, sync up dofobject unique_ids if applicable.
2139 129 : this->make_node_unique_ids_parallel_consistent(mesh);
2140 :
2141 : // Fourth, sync up any nodal boundary conditions
2142 129 : this->make_node_bcids_parallel_consistent(mesh);
2143 :
2144 : // Finally, correct the processor ids to make DofMap happy
2145 129 : MeshTools::correct_node_proc_ids(mesh);
2146 129 : }
2147 :
2148 :
2149 :
2150 : // ------------------------------------------------------------
2151 : void
2152 1396 : MeshCommunication::delete_remote_elements (DistributedMesh & mesh,
2153 : const std::set<Elem *> & extra_ghost_elem_ids) const
2154 : {
2155 : // The mesh should know it's about to be parallelized
2156 402 : libmesh_assert (!mesh.is_serial());
2157 :
2158 1206 : LOG_SCOPE("delete_remote_elements()", "MeshCommunication");
2159 :
2160 : #ifdef DEBUG
2161 : // We expect maximum ids to be in sync so we can use them to size
2162 : // vectors
2163 402 : libmesh_assert(mesh.comm().verify(mesh.max_node_id()));
2164 402 : libmesh_assert(mesh.comm().verify(mesh.max_elem_id()));
2165 402 : const dof_id_type par_max_node_id = mesh.parallel_max_node_id();
2166 402 : const dof_id_type par_max_elem_id = mesh.parallel_max_elem_id();
2167 402 : libmesh_assert_equal_to (par_max_node_id, mesh.max_node_id());
2168 402 : libmesh_assert_equal_to (par_max_elem_id, mesh.max_elem_id());
2169 :
2170 : // n_constraint_rows might not be cached, but we're in dbg mode,
2171 : // let's just (re)count them
2172 402 : const dof_id_type n_constraint_rows = mesh.count_constraint_rows();
2173 : #endif
2174 :
2175 804 : connected_elem_set_type elements_to_keep;
2176 :
2177 : // Don't delete elements that we were explicitly told not to
2178 1396 : for (const auto & elem : extra_ghost_elem_ids)
2179 : {
2180 0 : std::vector<const Elem *> active_family;
2181 : #ifdef LIBMESH_ENABLE_AMR
2182 0 : if (!elem->subactive())
2183 0 : elem->active_family_tree(active_family);
2184 : else
2185 : #endif
2186 0 : active_family.push_back(elem);
2187 :
2188 0 : for (const auto & f : active_family)
2189 0 : elements_to_keep.insert(f);
2190 : }
2191 :
2192 : // See which elements we still need to keep ghosted, given that
2193 : // we're keeping local and unpartitioned elements.
2194 : query_ghosting_functors
2195 3004 : (mesh, mesh.processor_id(),
2196 3194 : mesh.active_pid_elements_begin(mesh.processor_id()),
2197 2200 : mesh.active_pid_elements_end(mesh.processor_id()),
2198 : elements_to_keep);
2199 : query_ghosting_functors
2200 2602 : (mesh, DofObject::invalid_processor_id,
2201 2792 : mesh.active_pid_elements_begin(DofObject::invalid_processor_id),
2202 2390 : mesh.active_pid_elements_end(DofObject::invalid_processor_id),
2203 : elements_to_keep);
2204 :
2205 : // The inactive elements we need to send should have their
2206 : // immediate children present.
2207 3786 : connect_children(mesh, mesh.pid_elements_begin(mesh.processor_id()),
2208 2200 : mesh.pid_elements_end(mesh.processor_id()),
2209 : elements_to_keep);
2210 2602 : connect_children(mesh,
2211 2792 : mesh.pid_elements_begin(DofObject::invalid_processor_id),
2212 2792 : mesh.pid_elements_end(DofObject::invalid_processor_id),
2213 : elements_to_keep);
2214 :
2215 : // And see which elements and nodes they depend on
2216 804 : connected_node_set_type connected_nodes;
2217 1396 : connect_element_dependencies(mesh, elements_to_keep, connected_nodes);
2218 :
2219 : // Delete all the elements we have no reason to save,
2220 : // starting with the most refined so that the mesh
2221 : // is valid at all intermediate steps
2222 1396 : unsigned int n_levels = MeshTools::n_levels(mesh);
2223 :
2224 3119 : for (int l = n_levels - 1; l >= 0; --l)
2225 22282 : for (auto & elem : as_range(mesh.level_elements_begin(l),
2226 124192 : mesh.level_elements_end(l)))
2227 : {
2228 19332 : libmesh_assert (elem);
2229 : // Make sure we don't leave any invalid pointers
2230 56175 : const bool keep_me = elements_to_keep.count(elem);
2231 :
2232 19332 : if (!keep_me)
2233 5474 : elem->make_links_to_me_remote();
2234 :
2235 : // delete_elem doesn't currently invalidate element
2236 : // iterators... that had better not change
2237 38664 : if (!keep_me)
2238 5474 : mesh.delete_elem(elem);
2239 731 : }
2240 :
2241 : // Delete all the nodes we have no reason to save
2242 265166 : for (auto & node : mesh.node_ptr_range())
2243 : {
2244 51998 : libmesh_assert(node);
2245 152816 : if (!connected_nodes.count(node))
2246 : {
2247 3178 : libmesh_assert_not_equal_to(node->processor_id(),
2248 : mesh.processor_id());
2249 9545 : mesh.delete_node(node);
2250 : }
2251 592 : }
2252 :
2253 : // If we had a point locator, it's invalid now that some of the
2254 : // elements it pointed to have been deleted.
2255 1396 : mesh.clear_point_locator();
2256 :
2257 : // We now have all remote elements and nodes deleted; our ghosting
2258 : // functors should be ready to delete any now-redundant cached data
2259 : // they use too.
2260 3014 : for (auto & gf : as_range(mesh.ghosting_functors_begin(), mesh.ghosting_functors_end()))
2261 1618 : gf->delete_remote_elements();
2262 :
2263 : #ifdef DEBUG
2264 402 : const dof_id_type n_new_constraint_rows = mesh.count_constraint_rows();
2265 402 : libmesh_assert_equal_to(n_constraint_rows, n_new_constraint_rows);
2266 :
2267 402 : MeshTools::libmesh_assert_valid_refinement_tree(mesh);
2268 402 : MeshTools::libmesh_assert_valid_constraint_rows(mesh);
2269 : #endif
2270 1396 : }
2271 :
2272 : } // namespace libMesh
|
