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