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