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