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 1254569 : 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 1254569 : 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 674446 : for (const Elem * elem : connected_elements)
147 : {
148 3107538 : for (const Node & node : elem->node_ref_range())
149 : {
150 : // Retain all elements containing constraining nodes
151 2445264 : if (const auto it = constraint_rows.find(&node);
152 0 : it != constraint_rows.end())
153 26143658 : for (auto & p : it->second)
154 : {
155 23923962 : 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 12172 : 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 1070 : elem_rit = new_connected_elements.rbegin();
188 :
189 30480050 : for (; elem_rit != new_connected_elements.rend(); ++elem_rit)
190 : {
191 29225481 : const Elem * elem = *elem_rit;
192 25220 : 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 63305664 : for (const Elem * parent = elem->parent(); parent;
198 34068647 : parent = parent->parent())
199 27368 : if (!connected_elements.count(parent) &&
200 13684 : !new_connected_elements.count(parent))
201 18382446 : newer_connected_elements.insert (parent);
202 :
203 : auto total_family_insert =
204 29175041 : [&connected_elements, &new_connected_elements,
205 : &newer_connected_elements]
206 380288 : (const Elem * e)
207 : {
208 25220 : if (e->active() && e->has_children())
209 : {
210 0 : std::vector<const Elem *> subactive_family;
211 37444 : e->total_family_tree(subactive_family);
212 202368 : 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 29784 : newer_connected_elements.insert(f);
218 : }
219 : }
220 29225481 : };
221 :
222 29225481 : 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 29225481 : 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 25220 : libmesh_assert(!interior_parent || mesh);
233 :
234 38304 : if (interior_parent &&
235 13084 : interior_parent == mesh->query_elem_ptr(interior_parent->id()) &&
236 29225575 : !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 91808 : (const Elem * elem)
254 : {
255 54690 : libmesh_assert(elem);
256 54690 : const Elem * parent = elem->parent();
257 54690 : if (parent)
258 20106 : libmesh_assert(connected_elements.count(parent) ||
259 : new_connected_elements.count(parent) ||
260 : newer_connected_elements.count(parent));
261 54690 : };
262 :
263 29596 : for (const auto & elem : connected_elements)
264 28526 : check_elem(elem);
265 26290 : for (const auto & elem : new_connected_elements)
266 25220 : check_elem(elem);
267 2014 : for (const auto & elem : newer_connected_elements)
268 944 : check_elem(elem);
269 : # endif // DEBUG
270 :
271 : #endif // LIBMESH_ENABLE_AMR
272 1254569 : }
273 :
274 :
275 :
276 1254569 : 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 30480050 : for (const auto & elem : new_connected_elements)
282 313040014 : for (auto & n : elem->node_ref_range())
283 284040274 : if (!connected_nodes.count(&n) &&
284 294291697 : !new_connected_nodes.count(&n))
285 273523435 : newer_connected_nodes.insert(&n);
286 1254569 : }
287 :
288 :
289 : } // anonymous namespace
290 :
291 :
292 :
293 : namespace libMesh
294 : {
295 :
296 :
297 1158328 : 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 1157559 : for (auto & gf :
304 905 : as_range(mesh.ghosting_functors_begin(),
305 2730167 : mesh.ghosting_functors_end()))
306 : {
307 2082 : GhostingFunctor::map_type elements_to_ghost;
308 1041 : libmesh_assert(gf);
309 1569893 : (*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 18606609 : for (auto & pr : elements_to_ghost)
314 : {
315 17036716 : const Elem * elem = pr.first;
316 8651 : libmesh_assert(elem != remote_elem);
317 8651 : libmesh_assert(mesh.elem_ptr(elem->id()) == elem);
318 17028065 : 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 20399534 : for (; elem_it != elem_end; ++elem_it)
325 9635447 : connected_elements.insert(*elem_it);
326 1158328 : }
327 :
328 :
329 1158328 : 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 905 : 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 24991479 : for (const auto & elem : as_range(elem_it, elem_end))
342 : {
343 11353734 : if (elem->has_children())
344 9539345 : for (auto & child : elem->child_ref_range())
345 7894502 : if (&child != remote_elem)
346 7637226 : connected_elements.insert(&child);
347 1156518 : }
348 : #endif // LIBMESH_ENABLE_AMR
349 1158328 : }
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 758572 : 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 1066 : connected_elem_set_type new_connected_elements;
372 1066 : connected_node_set_type new_connected_nodes;
373 533 : new_connected_elements.swap(connected_elements);
374 533 : new_connected_nodes.swap(connected_nodes);
375 :
376 2014114 : while (!new_connected_elements.empty() ||
377 973 : !new_connected_nodes.empty())
378 : {
379 : auto [newer_connected_elements,
380 1070 : newer_connected_nodes] =
381 : connect_element_dependencies
382 : (mesh, connected_elements, connected_nodes,
383 1256709 : new_connected_elements, new_connected_nodes);
384 :
385 : // These have now been examined
386 1070 : connected_elements.merge(new_connected_elements);
387 1070 : 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 1070 : libmesh_assert(new_connected_elements.empty());
392 1070 : libmesh_assert(new_connected_nodes.empty());
393 :
394 : // These now need to be examined
395 1070 : new_connected_elements.swap(newer_connected_elements);
396 1070 : new_connected_nodes.swap(newer_connected_nodes);
397 1252429 : }
398 758572 : }
399 :
400 :
401 : std::pair<connected_elem_set_type, connected_node_set_type>
402 1254569 : 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 1070 : std::pair<connected_elem_set_type, connected_node_set_type> returnval;
409 1070 : auto & [newer_connected_elements, newer_connected_nodes] = returnval;
410 1254569 : connect_element_families(connected_elements, new_connected_elements,
411 : newer_connected_elements, &mesh);
412 :
413 1254569 : connect_nodes(new_connected_elements, connected_nodes,
414 : new_connected_nodes, newer_connected_nodes);
415 :
416 1254569 : 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 66243 : 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 66243 : 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 66243 : 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 10032400 : for (auto & elem : as_range(send_elems_begin, send_elems_end))
496 4973033 : if (elem->processor_id() != mesh.processor_id())
497 3629010 : 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 66243 : bool have_constraint_rows = !constraint_rows.empty();
505 66243 : 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 424519 : 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 358276 : v_t * elempp = p_elements.data();
536 358276 : 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 358409 : (elempp, elemend, Predicates::NotNull<v_t *>());
547 :
548 : const MeshBase::const_element_iterator elem_end =
549 : MeshBase::const_element_iterator
550 716552 : (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 716419 : 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 1074562 : connect_children(mesh, mesh.pid_elements_begin(pid),
561 716552 : 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 358276 : connect_element_dependencies(mesh, elements_to_send,
567 : connected_nodes);
568 :
569 358276 : all_nodes_to_send[pid].assign(connected_nodes.begin(),
570 : connected_nodes.end());
571 :
572 358276 : all_elems_to_send[pid].assign(elements_to_send.begin(),
573 : elements_to_send.end());
574 :
575 448534 : for (auto & [node, row] : constraint_rows)
576 : {
577 29562 : if (!connected_nodes.count(node))
578 29562 : continue;
579 :
580 0 : serialized_row_type serialized_row;
581 648454 : for (auto [elem_and_node, coef] : row)
582 587758 : 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 60696 : (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 66243 : TIMPI::push_parallel_packed_range(mesh.comm(), all_nodes_to_send, &mesh,
597 : null_node_action);
598 :
599 66243 : 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 66243 : if (have_constraint_rows)
605 : {
606 : auto constraint_row_action =
607 1740 : [&mesh, &constraint_rows]
608 : (processor_id_type /* src_pid */,
609 647866 : const std::vector<std::pair<dof_id_type, serialized_row_type>> rows)
610 : {
611 62240 : for (auto & [node_id, serialized_row] : rows)
612 : {
613 0 : MeshBase::constraint_rows_mapped_type row;
614 647866 : for (auto [elem_and_node, coef] : serialized_row)
615 587366 : row.emplace_back(std::make_pair(mesh.elem_ptr(elem_and_node.first),
616 0 : elem_and_node.second),
617 0 : coef);
618 :
619 60500 : constraint_rows[mesh.node_ptr(node_id)] = row;
620 : }
621 1959 : };
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 66243 : 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 66243 : mesh.MeshBase::redistribute();
651 66243 : }
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 7015 : 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 7015 : if (mesh.is_serial())
1049 5851 : 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 3044 : for (auto elem : as_range(mesh.flagged_pid_elements_begin(Elem::COARSEN, proc_id),
1070 212399 : 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 103920 : const processor_id_type their_proc_id = elem->parent()->processor_id();
1080 103920 : if (their_proc_id != proc_id)
1081 631 : coarsening_elements_to_ghost[their_proc_id].push_back(elem);
1082 1156 : }
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 13428 : for (processor_id_type p=0; p != n_proc; ++p)
1090 : {
1091 12264 : if (p == proc_id)
1092 1164 : continue;
1093 :
1094 8 : connected_elem_set_type elements_to_send;
1095 8 : std::set<const Node *> nodes_to_send;
1096 :
1097 11100 : 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 201 : Elem * const * elempp = const_cast<Elem * const *>(elems.data());
1106 201 : Elem * const * elemend = elempp+elems.size();
1107 : const MeshBase::const_element_iterator elem_it =
1108 201 : MeshBase::const_element_iterator(elempp, elemend, Predicates::NotNull<Elem * const *>());
1109 : const MeshBase::const_element_iterator elem_end =
1110 402 : MeshBase::const_element_iterator(elemend, elemend, Predicates::NotNull<Elem * const *>());
1111 :
1112 201 : for (auto & gf : as_range(mesh.ghosting_functors_begin(),
1113 828 : mesh.ghosting_functors_end()))
1114 : {
1115 0 : GhostingFunctor::map_type elements_to_ghost;
1116 0 : libmesh_assert(gf);
1117 627 : (*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 4941 : for (auto & pr : elements_to_ghost)
1122 : {
1123 4314 : const Elem * elem = pr.first;
1124 0 : libmesh_assert(elem);
1125 15059 : while (elem)
1126 : {
1127 0 : libmesh_assert(elem != remote_elem);
1128 10745 : elements_to_send.insert(elem);
1129 104948 : for (auto & n : elem->node_ref_range())
1130 94203 : nodes_to_send.insert(&n);
1131 10745 : elem = elem->parent();
1132 : }
1133 : }
1134 : }
1135 :
1136 201 : all_nodes_to_send[p].assign(nodes_to_send.begin(), nodes_to_send.end());
1137 201 : 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 1164 : TIMPI::push_parallel_packed_range(mesh.comm(), all_nodes_to_send, &mesh,
1147 : null_node_action);
1148 :
1149 1164 : 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 14582 : void MeshCommunication::broadcast (MeshBase & mesh) const
1167 : {
1168 : // Don't need to do anything if there is
1169 : // only one processor.
1170 14998 : if (mesh.n_processors() == 1)
1171 403 : return;
1172 :
1173 : // This function must be run on all processors at once
1174 416 : libmesh_parallel_only(mesh.comm());
1175 :
1176 832 : LOG_SCOPE("broadcast()", "MeshCommunication");
1177 :
1178 : // Explicitly clear the mesh on all but processor 0.
1179 14595 : if (mesh.processor_id() != 0)
1180 12112 : mesh.clear();
1181 :
1182 : // We may have set extra data only on processor 0 in a read()
1183 14179 : mesh.comm().broadcast(mesh._node_integer_names);
1184 14179 : mesh.comm().broadcast(mesh._node_integer_default_values);
1185 14179 : mesh.comm().broadcast(mesh._elem_integer_names);
1186 14179 : mesh.comm().broadcast(mesh._elem_integer_default_values);
1187 :
1188 : // We may have set mapping data only on processor 0 in a read()
1189 14179 : unsigned char map_type = mesh.default_mapping_type();
1190 14179 : unsigned char map_data = mesh.default_mapping_data();
1191 14179 : mesh.comm().broadcast(map_type);
1192 14179 : mesh.comm().broadcast(map_data);
1193 14179 : mesh.set_default_mapping_type(ElemMappingType(map_type));
1194 14179 : mesh.set_default_mapping_data(map_data);
1195 :
1196 : // Broadcast nodes
1197 14595 : mesh.comm().broadcast_packed_range(&mesh,
1198 28358 : mesh.nodes_begin(),
1199 14595 : 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 14179 : const unsigned int n_levels = MeshTools::paranoid_n_levels(mesh);
1211 :
1212 28358 : for (unsigned int l=0; l != n_levels; ++l)
1213 14595 : mesh.comm().broadcast_packed_range(&mesh,
1214 28358 : mesh.level_elements_begin(l),
1215 28358 : mesh.level_elements_end(l),
1216 : &mesh,
1217 : null_output_iterator<Elem>());
1218 :
1219 : // Make sure mesh_dimension and elem_dimensions are consistent.
1220 14179 : mesh.cache_elem_data();
1221 :
1222 : // Make sure mesh id counts are consistent.
1223 14179 : 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 416 : auto & constraint_rows = mesh.get_constraint_rows();
1230 14179 : bool have_constraint_rows = !constraint_rows.empty();
1231 14179 : mesh.comm().broadcast(have_constraint_rows);
1232 14179 : 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 14179 : mesh.comm().broadcast(mesh.set_subdomain_name_map());
1279 14179 : mesh.comm().broadcast(mesh.get_boundary_info().set_sideset_name_map());
1280 14179 : 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 14179 : mesh.clear_point_locator();
1285 :
1286 416 : libmesh_assert (mesh.comm().verify(mesh.n_elem()));
1287 416 : libmesh_assert (mesh.comm().verify(mesh.n_nodes()));
1288 :
1289 : #ifdef DEBUG
1290 416 : MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
1291 416 : MeshTools::libmesh_assert_valid_procids<Node>(mesh);
1292 : #endif
1293 : }
1294 : #endif // LIBMESH_HAVE_MPI
1295 :
1296 :
1297 :
1298 : #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
1299 : // ------------------------------------------------------------
1300 : void MeshCommunication::gather (const processor_id_type, MeshBase &) const
1301 : {
1302 : // no MPI == one processor, no need for this method...
1303 : return;
1304 : }
1305 : #else
1306 : // ------------------------------------------------------------
1307 110460 : void MeshCommunication::gather (const processor_id_type root_id, MeshBase & mesh) const
1308 : {
1309 : // Check for quick return
1310 110532 : if (mesh.n_processors() == 1)
1311 1093 : return;
1312 :
1313 : // This function must be run on all processors at once
1314 72 : libmesh_parallel_only(mesh.comm());
1315 :
1316 144 : LOG_SCOPE("(all)gather()", "MeshCommunication");
1317 :
1318 : // Ensure we don't build too big a buffer at once
1319 : static const std::size_t approx_total_buffer_size = 1e8;
1320 : const std::size_t approx_each_buffer_size =
1321 109367 : approx_total_buffer_size / mesh.comm().size();
1322 :
1323 109367 : (root_id == DofObject::invalid_processor_id) ?
1324 :
1325 95791 : mesh.comm().allgather_packed_range (&mesh,
1326 95791 : mesh.nodes_begin(),
1327 205096 : mesh.nodes_end(),
1328 : null_output_iterator<Node>(),
1329 : approx_each_buffer_size) :
1330 :
1331 13648 : mesh.comm().gather_packed_range (root_id,
1332 : &mesh,
1333 123005 : mesh.nodes_begin(),
1334 136633 : mesh.nodes_end(),
1335 : null_output_iterator<Node>(),
1336 : approx_each_buffer_size);
1337 :
1338 : // Gather elements from coarsest to finest, so that child
1339 : // elements will see their parents already in place.
1340 109367 : const unsigned int n_levels = MeshTools::n_levels(mesh);
1341 :
1342 298073 : for (unsigned int l=0; l != n_levels; ++l)
1343 188706 : (root_id == DofObject::invalid_processor_id) ?
1344 :
1345 168452 : mesh.comm().allgather_packed_range (&mesh,
1346 168452 : mesh.level_elements_begin(l),
1347 525388 : mesh.level_elements_end(l),
1348 : null_output_iterator<Elem>(),
1349 : approx_each_buffer_size) :
1350 :
1351 20342 : mesh.comm().gather_packed_range (root_id,
1352 : &mesh,
1353 209034 : mesh.level_elements_begin(l),
1354 229348 : mesh.level_elements_end(l),
1355 : null_output_iterator<Elem>(),
1356 : approx_each_buffer_size);
1357 :
1358 : // If we had a point locator, it's invalid now that there are new
1359 : // elements it can't locate.
1360 109367 : mesh.clear_point_locator();
1361 :
1362 : // We may have constraint rows on IsoGeometric Analysis meshes. We
1363 : // don't want to send these along with constrained nodes (like we
1364 : // send boundary info for those nodes) because the associated rows'
1365 : // elements may not exist at that point.
1366 72 : auto & constraint_rows = mesh.get_constraint_rows();
1367 109367 : bool have_constraint_rows = !constraint_rows.empty();
1368 109367 : mesh.comm().max(have_constraint_rows);
1369 109367 : if (have_constraint_rows)
1370 : {
1371 : std::map<dof_id_type,
1372 : std::vector<std::tuple<dof_id_type, unsigned int, Real>>>
1373 0 : serialized_rows;
1374 :
1375 361 : for (auto & row : constraint_rows)
1376 : {
1377 238 : const Node * node = row.first;
1378 238 : const dof_id_type rowid = node->id();
1379 0 : libmesh_assert(node == mesh.node_ptr(rowid));
1380 :
1381 : std::vector<std::tuple<dof_id_type, unsigned int, Real>>
1382 0 : serialized_row;
1383 714 : for (auto & entry : row.second)
1384 : serialized_row.push_back
1385 476 : (std::make_tuple(entry.first.first->id(),
1386 0 : entry.first.second, entry.second));
1387 :
1388 238 : serialized_rows.emplace(rowid, std::move(serialized_row));
1389 : }
1390 :
1391 123 : if (root_id == DofObject::invalid_processor_id)
1392 123 : mesh.comm().set_union(serialized_rows);
1393 : else
1394 0 : mesh.comm().set_union(serialized_rows, root_id);
1395 :
1396 123 : if (root_id == DofObject::invalid_processor_id ||
1397 0 : root_id == mesh.processor_id())
1398 : {
1399 918 : for (auto & row : serialized_rows)
1400 : {
1401 795 : const dof_id_type rowid = row.first;
1402 795 : const Node * node = mesh.node_ptr(rowid);
1403 :
1404 : std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>>
1405 0 : deserialized_row;
1406 2385 : for (auto & entry : row.second)
1407 : deserialized_row.push_back
1408 3180 : (std::make_pair(std::make_pair(mesh.elem_ptr(std::get<0>(entry)),
1409 0 : std::get<1>(entry)),
1410 0 : std::get<2>(entry)));
1411 :
1412 795 : constraint_rows.emplace(node, deserialized_row);
1413 : }
1414 : }
1415 : #ifdef DEBUG
1416 0 : MeshTools::libmesh_assert_valid_constraint_rows(mesh);
1417 : #endif
1418 : }
1419 :
1420 :
1421 : // If we are doing an allgather(), perform sanity check on the result.
1422 72 : if (root_id == DofObject::invalid_processor_id)
1423 : {
1424 62 : libmesh_assert (mesh.comm().verify(mesh.n_elem()));
1425 62 : libmesh_assert (mesh.comm().verify(mesh.n_nodes()));
1426 : }
1427 :
1428 : // Inform new elements of their neighbors,
1429 : // while resetting all remote_elem links on
1430 : // the ranks which did the gather.
1431 109367 : if (mesh.allow_find_neighbors())
1432 196506 : mesh.find_neighbors(root_id == DofObject::invalid_processor_id ||
1433 144 : root_id == mesh.processor_id());
1434 :
1435 : // All done, but let's make sure it's done correctly
1436 :
1437 : #ifdef DEBUG
1438 72 : MeshTools::libmesh_assert_valid_boundary_ids(mesh);
1439 : #endif
1440 : }
1441 : #endif // LIBMESH_HAVE_MPI
1442 :
1443 :
1444 :
1445 : // Functor for make_elems_parallel_consistent and
1446 : // make_node_ids_parallel_consistent
1447 : namespace {
1448 :
1449 : struct SyncIds
1450 : {
1451 : typedef dof_id_type datum;
1452 : typedef void (MeshBase::*renumber_obj)(dof_id_type, dof_id_type);
1453 :
1454 22145 : SyncIds(MeshBase & _mesh, renumber_obj _renumberer) :
1455 22069 : mesh(_mesh),
1456 22145 : renumber(_renumberer) {}
1457 :
1458 : MeshBase & mesh;
1459 : renumber_obj renumber;
1460 : // renumber_obj & renumber;
1461 :
1462 : // Find the id of each requested DofObject -
1463 : // Parallel::sync_* already did the work for us
1464 34 : void gather_data (const std::vector<dof_id_type> & ids,
1465 : std::vector<datum> & ids_out) const
1466 : {
1467 65902 : ids_out = ids;
1468 65868 : }
1469 :
1470 65902 : void act_on_data (const std::vector<dof_id_type> & old_ids,
1471 : const std::vector<datum> & new_ids) const
1472 : {
1473 2716038 : for (auto i : index_range(old_ids))
1474 2652600 : if (old_ids[i] != new_ids[i])
1475 2604934 : (mesh.*renumber)(old_ids[i], new_ids[i]);
1476 65902 : }
1477 : };
1478 :
1479 :
1480 30210 : struct SyncNodeIds
1481 : {
1482 : typedef dof_id_type datum;
1483 :
1484 30258 : SyncNodeIds(MeshBase & _mesh) :
1485 30258 : mesh(_mesh) {}
1486 :
1487 : MeshBase & mesh;
1488 :
1489 : // We only know a Node id() is definitive if we own the Node or if
1490 : // we're told it's definitive. We keep track of the latter cases by
1491 : // putting definitively id'd ghost nodes into this set.
1492 : typedef std::unordered_set<const Node *> uset_type;
1493 : uset_type definitive_ids;
1494 :
1495 : // We should never be told two different definitive ids for the same
1496 : // node, but let's check on that in debug mode.
1497 : #ifdef DEBUG
1498 : typedef std::unordered_map<dof_id_type, dof_id_type> umap_type;
1499 : umap_type definitive_renumbering;
1500 : #endif
1501 :
1502 : // Find the id of each requested DofObject -
1503 : // Parallel::sync_* already tried to do the work for us, but we can
1504 : // only say the result is definitive if we own the DofObject or if
1505 : // we were given the definitive result from another processor.
1506 286238 : void gather_data (const std::vector<dof_id_type> & ids,
1507 : std::vector<datum> & ids_out) const
1508 : {
1509 132 : ids_out.clear();
1510 286370 : ids_out.resize(ids.size(), DofObject::invalid_id);
1511 :
1512 59729949 : for (auto i : index_range(ids))
1513 : {
1514 59443711 : const dof_id_type id = ids[i];
1515 59443711 : const Node * node = mesh.query_node_ptr(id);
1516 59443711 : if (node && (node->processor_id() == mesh.processor_id() ||
1517 33924 : definitive_ids.count(node)))
1518 59477351 : ids_out[i] = id;
1519 : }
1520 286238 : }
1521 :
1522 286238 : bool act_on_data (const std::vector<dof_id_type> & old_ids,
1523 : const std::vector<datum> & new_ids)
1524 : {
1525 132 : bool data_changed = false;
1526 59729949 : for (auto i : index_range(old_ids))
1527 : {
1528 59443711 : const dof_id_type new_id = new_ids[i];
1529 :
1530 59443711 : const dof_id_type old_id = old_ids[i];
1531 :
1532 59443711 : Node * node = mesh.query_node_ptr(old_id);
1533 :
1534 : // If we can't find the node we were asking about, another
1535 : // processor must have already given us the definitive id
1536 : // for it
1537 59443711 : if (!node)
1538 : {
1539 : // But let's check anyway in debug mode
1540 : #ifdef DEBUG
1541 5626 : libmesh_assert
1542 : (definitive_renumbering.count(old_id));
1543 5626 : libmesh_assert_equal_to
1544 : (new_id, definitive_renumbering[old_id]);
1545 : #endif
1546 11346984 : continue;
1547 : }
1548 :
1549 : // If we asked for an id but there's no definitive id ready
1550 : // for us yet, then we can't quit trying to sync yet.
1551 48091101 : if (new_id == DofObject::invalid_id)
1552 : {
1553 : // But we might have gotten a definitive id from a
1554 : // different request
1555 284 : if (!definitive_ids.count(mesh.node_ptr(old_id)))
1556 0 : data_changed = true;
1557 : }
1558 : else
1559 : {
1560 48119115 : if (node->processor_id() != mesh.processor_id())
1561 46812 : definitive_ids.insert(node);
1562 48090817 : if (old_id != new_id)
1563 : {
1564 : #ifdef DEBUG
1565 3234 : libmesh_assert
1566 : (!definitive_renumbering.count(old_id));
1567 3234 : definitive_renumbering[old_id] = new_id;
1568 : #endif
1569 4745707 : mesh.renumber_node(old_id, new_id);
1570 3234 : data_changed = true;
1571 : }
1572 : }
1573 : }
1574 286238 : return data_changed;
1575 : }
1576 : };
1577 :
1578 :
1579 : #ifdef LIBMESH_ENABLE_AMR
1580 : struct SyncPLevels
1581 : {
1582 : typedef std::pair<unsigned char,unsigned char> datum;
1583 :
1584 2560 : SyncPLevels(MeshBase & _mesh) :
1585 2560 : mesh(_mesh) {}
1586 :
1587 : MeshBase & mesh;
1588 :
1589 : // Find the p_level of each requested Elem
1590 3346 : void gather_data (const std::vector<dof_id_type> & ids,
1591 : std::vector<datum> & ids_out) const
1592 : {
1593 3346 : ids_out.reserve(ids.size());
1594 :
1595 18571 : for (const auto & id : ids)
1596 : {
1597 15225 : Elem & elem = mesh.elem_ref(id);
1598 : ids_out.push_back
1599 15225 : (std::make_pair(cast_int<unsigned char>(elem.p_level()),
1600 0 : static_cast<unsigned char>(elem.p_refinement_flag())));
1601 : }
1602 3346 : }
1603 :
1604 3346 : void act_on_data (const std::vector<dof_id_type> & old_ids,
1605 : const std::vector<datum> & new_p_levels) const
1606 : {
1607 18571 : for (auto i : index_range(old_ids))
1608 : {
1609 15225 : Elem & elem = mesh.elem_ref(old_ids[i]);
1610 : // Make sure these are consistent
1611 : elem.hack_p_level_and_refinement_flag
1612 15225 : (new_p_levels[i].first,
1613 15225 : static_cast<Elem::RefinementState>(new_p_levels[i].second));
1614 : // Make sure parents' levels are consistent
1615 15225 : elem.set_p_level(new_p_levels[i].first);
1616 : }
1617 3346 : }
1618 : };
1619 : #endif // LIBMESH_ENABLE_AMR
1620 :
1621 :
1622 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1623 : template <typename DofObjSubclass>
1624 : struct SyncUniqueIds
1625 : {
1626 : typedef unique_id_type datum;
1627 : typedef DofObjSubclass* (MeshBase::*query_obj)(const dof_id_type);
1628 :
1629 55286 : SyncUniqueIds(MeshBase &_mesh, query_obj _querier) :
1630 55038 : mesh(_mesh),
1631 55286 : query(_querier) {}
1632 :
1633 : MeshBase & mesh;
1634 : query_obj query;
1635 :
1636 : // Find the id of each requested DofObject -
1637 : // Parallel::sync_* already did the work for us
1638 203964 : void gather_data (const std::vector<dof_id_type> & ids,
1639 : std::vector<datum> & ids_out) const
1640 : {
1641 204060 : ids_out.reserve(ids.size());
1642 :
1643 9687481 : for (const auto & id : ids)
1644 : {
1645 9483517 : DofObjSubclass * d = (mesh.*query)(id);
1646 6712 : libmesh_assert(d);
1647 9483517 : ids_out.push_back(d->unique_id());
1648 : }
1649 203964 : }
1650 :
1651 203964 : void act_on_data (const std::vector<dof_id_type> & ids,
1652 : const std::vector<datum> & unique_ids) const
1653 : {
1654 9687481 : for (auto i : index_range(ids))
1655 : {
1656 9483517 : DofObjSubclass * d = (mesh.*query)(ids[i]);
1657 6712 : libmesh_assert(d);
1658 9490229 : d->set_unique_id(unique_ids[i]);
1659 : }
1660 203964 : }
1661 : };
1662 : #endif // LIBMESH_ENABLE_UNIQUE_ID
1663 :
1664 : template <typename DofObjSubclass>
1665 : struct SyncBCIds
1666 : {
1667 : typedef std::vector<boundary_id_type> datum;
1668 :
1669 22145 : SyncBCIds(MeshBase &_mesh) :
1670 22145 : mesh(_mesh) {}
1671 :
1672 : MeshBase & mesh;
1673 :
1674 : // Find the id of each requested DofObject -
1675 : // Parallel::sync_* already did the work for us
1676 90569 : void gather_data (const std::vector<dof_id_type> & ids,
1677 : std::vector<datum> & ids_out) const
1678 : {
1679 90569 : ids_out.reserve(ids.size());
1680 :
1681 90569 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1682 :
1683 6122064 : for (const auto & id : ids)
1684 : {
1685 6031495 : Node * n = mesh.query_node_ptr(id);
1686 2981 : libmesh_assert(n);
1687 5962 : std::vector<boundary_id_type> bcids;
1688 6031495 : boundary_info.boundary_ids(n, bcids);
1689 2981 : ids_out.push_back(std::move(bcids));
1690 : }
1691 90569 : }
1692 :
1693 90569 : void act_on_data (const std::vector<dof_id_type> & ids,
1694 : const std::vector<datum> & bcids) const
1695 : {
1696 90569 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
1697 :
1698 6122064 : for (auto i : index_range(ids))
1699 : {
1700 6031495 : Node * n = mesh.query_node_ptr(ids[i]);
1701 2981 : libmesh_assert(n);
1702 6031495 : boundary_info.add_node(n, bcids[i]);
1703 : }
1704 90569 : }
1705 : };
1706 :
1707 : }
1708 :
1709 :
1710 :
1711 : // ------------------------------------------------------------
1712 30258 : void MeshCommunication::make_node_ids_parallel_consistent (MeshBase & mesh)
1713 : {
1714 : // This function must be run on all processors at once
1715 48 : libmesh_parallel_only(mesh.comm());
1716 :
1717 : // We need to agree on which processor owns every node, but we can't
1718 : // easily assert that here because we don't currently agree on which
1719 : // id every node has, and some of our temporary ids on unrelated
1720 : // nodes will "overlap".
1721 : //#ifdef DEBUG
1722 : // MeshTools::libmesh_assert_parallel_consistent_procids<Node> (mesh);
1723 : //#endif // DEBUG
1724 :
1725 96 : LOG_SCOPE ("make_node_ids_parallel_consistent()", "MeshCommunication");
1726 :
1727 96 : SyncNodeIds syncids(mesh);
1728 : Parallel::sync_node_data_by_element_id
1729 90678 : (mesh, mesh.elements_begin(), mesh.elements_end(),
1730 30210 : Parallel::SyncEverything(), Parallel::SyncEverything(), syncids);
1731 :
1732 : // At this point, with both ids and processor ids synced, we can
1733 : // finally check for topological consistency of node processor ids.
1734 : #ifdef DEBUG
1735 48 : MeshTools::libmesh_assert_topology_consistent_procids<Node> (mesh);
1736 : #endif
1737 30258 : }
1738 :
1739 :
1740 :
1741 33141 : void MeshCommunication::make_node_unique_ids_parallel_consistent (MeshBase & mesh)
1742 : {
1743 : // Avoid unused variable warnings if unique ids aren't enabled.
1744 86 : libmesh_ignore(mesh);
1745 :
1746 : // This function must be run on all processors at once
1747 86 : libmesh_parallel_only(mesh.comm());
1748 :
1749 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1750 86 : LOG_SCOPE ("make_node_unique_ids_parallel_consistent()", "MeshCommunication");
1751 :
1752 86 : SyncUniqueIds<Node> syncuniqueids(mesh, &MeshBase::query_node_ptr);
1753 172 : Parallel::sync_dofobject_data_by_id(mesh.comm(),
1754 66282 : mesh.nodes_begin(),
1755 33313 : mesh.nodes_end(),
1756 : syncuniqueids);
1757 :
1758 : #endif
1759 33141 : }
1760 :
1761 :
1762 22145 : void MeshCommunication::make_node_bcids_parallel_consistent (MeshBase & mesh)
1763 : {
1764 : // Avoid unused variable warnings if unique ids aren't enabled.
1765 38 : libmesh_ignore(mesh);
1766 :
1767 : // This function must be run on all processors at once
1768 38 : libmesh_parallel_only(mesh.comm());
1769 :
1770 38 : LOG_SCOPE ("make_node_bcids_parallel_consistent()", "MeshCommunication");
1771 :
1772 38 : SyncBCIds<Node> syncbcids(mesh);
1773 76 : Parallel::sync_dofobject_data_by_id(mesh.comm(),
1774 44290 : mesh.nodes_begin(),
1775 22221 : mesh.nodes_end(),
1776 : syncbcids);
1777 22145 : }
1778 :
1779 :
1780 :
1781 :
1782 :
1783 : // ------------------------------------------------------------
1784 22145 : void MeshCommunication::make_elems_parallel_consistent(MeshBase & mesh)
1785 : {
1786 : // This function must be run on all processors at once
1787 38 : libmesh_parallel_only(mesh.comm());
1788 :
1789 38 : LOG_SCOPE ("make_elems_parallel_consistent()", "MeshCommunication");
1790 :
1791 38 : SyncIds syncids(mesh, &MeshBase::renumber_elem);
1792 : Parallel::sync_element_data_by_parent_id
1793 44252 : (mesh, mesh.active_elements_begin(),
1794 22183 : mesh.active_elements_end(), syncids);
1795 :
1796 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1797 38 : SyncUniqueIds<Elem> syncuniqueids(mesh, &MeshBase::query_elem_ptr);
1798 : Parallel::sync_dofobject_data_by_id
1799 44252 : (mesh.comm(), mesh.active_elements_begin(),
1800 22221 : mesh.active_elements_end(), syncuniqueids);
1801 : #endif
1802 22145 : }
1803 :
1804 :
1805 :
1806 : // ------------------------------------------------------------
1807 : #ifdef LIBMESH_ENABLE_AMR
1808 2560 : void MeshCommunication::make_p_levels_parallel_consistent(MeshBase & mesh)
1809 : {
1810 : // This function must be run on all processors at once
1811 0 : libmesh_parallel_only(mesh.comm());
1812 :
1813 0 : LOG_SCOPE ("make_p_levels_parallel_consistent()", "MeshCommunication");
1814 :
1815 0 : SyncPLevels syncplevels(mesh);
1816 : Parallel::sync_dofobject_data_by_id
1817 5120 : (mesh.comm(), mesh.elements_begin(), mesh.elements_end(),
1818 : syncplevels);
1819 2560 : }
1820 : #endif // LIBMESH_ENABLE_AMR
1821 :
1822 :
1823 :
1824 : // Functors for make_node_proc_ids_parallel_consistent
1825 : namespace {
1826 :
1827 : struct SyncProcIds
1828 : {
1829 : typedef processor_id_type datum;
1830 :
1831 30258 : SyncProcIds(MeshBase & _mesh) : mesh(_mesh) {}
1832 :
1833 : MeshBase & mesh;
1834 :
1835 : // ------------------------------------------------------------
1836 236654 : void gather_data (const std::vector<dof_id_type> & ids,
1837 : std::vector<datum> & data)
1838 : {
1839 : // Find the processor id of each requested node
1840 236654 : data.resize(ids.size());
1841 :
1842 49509857 : for (auto i : index_range(ids))
1843 : {
1844 : // Look for this point in the mesh
1845 49273203 : if (ids[i] != DofObject::invalid_id)
1846 : {
1847 49273203 : Node & node = mesh.node_ref(ids[i]);
1848 :
1849 : // Return the node's correct processor id,
1850 49273203 : data[i] = node.processor_id();
1851 : }
1852 : else
1853 0 : data[i] = DofObject::invalid_processor_id;
1854 : }
1855 236654 : }
1856 :
1857 : // ------------------------------------------------------------
1858 236654 : bool act_on_data (const std::vector<dof_id_type> & ids,
1859 : const std::vector<datum> proc_ids)
1860 : {
1861 86 : bool data_changed = false;
1862 :
1863 : // Set the ghost node processor ids we've now been informed of
1864 49509857 : for (auto i : index_range(ids))
1865 : {
1866 49273203 : Node & node = mesh.node_ref(ids[i]);
1867 :
1868 : // We may not have ids synched when this synchronization is done, so we
1869 : // *can't* use id to load-balance processor id properly; we have to use
1870 : // the old heuristic of choosing the smallest valid processor id.
1871 : //
1872 : // If someone tells us our node processor id is too low, then
1873 : // they're wrong. If they tell us our node processor id is
1874 : // too high, then we're wrong.
1875 49273203 : if (node.processor_id() > proc_ids[i])
1876 : {
1877 0 : data_changed = true;
1878 80126 : node.processor_id() = proc_ids[i];
1879 : }
1880 : }
1881 :
1882 236654 : return data_changed;
1883 : }
1884 : };
1885 :
1886 :
1887 : struct ElemNodesMaybeNew
1888 : {
1889 38 : ElemNodesMaybeNew() {}
1890 :
1891 16740018 : bool operator() (const Elem * elem) const
1892 : {
1893 : // If this element was just refined then it may have new nodes we
1894 : // need to work on
1895 : #ifdef LIBMESH_ENABLE_AMR
1896 16740018 : if (elem->refinement_flag() == Elem::JUST_REFINED)
1897 7896 : return true;
1898 : #endif
1899 :
1900 : // If this element has remote_elem neighbors then there may have
1901 : // been refinement of those neighbors that affect its nodes'
1902 : // processor_id()
1903 13926636 : for (auto neigh : elem->neighbor_ptr_range())
1904 11589570 : if (neigh == remote_elem)
1905 880054 : return true;
1906 2337066 : return false;
1907 : }
1908 : };
1909 :
1910 :
1911 22107 : struct NodeWasNew
1912 : {
1913 22145 : NodeWasNew(const MeshBase & mesh)
1914 22145 : {
1915 17142230 : for (const auto & node : mesh.node_ptr_range())
1916 8561335 : if (node->processor_id() == DofObject::invalid_processor_id)
1917 32547 : was_new.insert(node);
1918 22145 : }
1919 :
1920 93800260 : bool operator() (const Elem * elem, unsigned int local_node_num) const
1921 : {
1922 93800260 : if (was_new.count(elem->node_ptr(local_node_num)))
1923 60794062 : return true;
1924 9602 : return false;
1925 : }
1926 :
1927 : std::unordered_set<const Node *> was_new;
1928 : };
1929 :
1930 : }
1931 :
1932 :
1933 :
1934 : // ------------------------------------------------------------
1935 8113 : void MeshCommunication::make_node_proc_ids_parallel_consistent(MeshBase & mesh)
1936 : {
1937 10 : LOG_SCOPE ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");
1938 :
1939 : // This function must be run on all processors at once
1940 10 : libmesh_parallel_only(mesh.comm());
1941 :
1942 : // When this function is called, each section of a parallelized mesh
1943 : // should be in the following state:
1944 : //
1945 : // All nodes should have the exact same physical location on every
1946 : // processor where they exist.
1947 : //
1948 : // Local nodes should have unique authoritative ids,
1949 : // and processor ids consistent with all processors which own
1950 : // an element touching them.
1951 : //
1952 : // Ghost nodes touching local elements should have processor ids
1953 : // consistent with all processors which own an element touching
1954 : // them.
1955 10 : SyncProcIds sync(mesh);
1956 : Parallel::sync_node_data_by_element_id
1957 24319 : (mesh, mesh.elements_begin(), mesh.elements_end(),
1958 8113 : Parallel::SyncEverything(), Parallel::SyncEverything(), sync);
1959 8113 : }
1960 :
1961 :
1962 :
1963 : // ------------------------------------------------------------
1964 22145 : void MeshCommunication::make_new_node_proc_ids_parallel_consistent(MeshBase & mesh)
1965 : {
1966 76 : LOG_SCOPE ("make_new_node_proc_ids_parallel_consistent()", "MeshCommunication");
1967 :
1968 : // This function must be run on all processors at once
1969 38 : libmesh_parallel_only(mesh.comm());
1970 :
1971 : // When this function is called, each section of a parallelized mesh
1972 : // should be in the following state:
1973 : //
1974 : // Local nodes should have unique authoritative ids,
1975 : // and new nodes should be unpartitioned.
1976 : //
1977 : // New ghost nodes touching local elements should be unpartitioned.
1978 :
1979 : // We may not have consistent processor ids for new nodes (because a
1980 : // node may be old and partitioned on one processor but new and
1981 : // unpartitioned on another) when we start
1982 : #ifdef DEBUG
1983 38 : MeshTools::libmesh_assert_parallel_consistent_procids<Node>(mesh);
1984 : // MeshTools::libmesh_assert_parallel_consistent_new_node_procids(mesh);
1985 : #endif
1986 :
1987 : // We have two kinds of new nodes. *NEW* nodes are unpartitioned on
1988 : // all processors: we need to use a id-independent (i.e. dumb)
1989 : // heuristic to partition them. But "new" nodes are newly created
1990 : // on some processors (when ghost elements are refined) yet
1991 : // correspond to existing nodes on other processors: we need to use
1992 : // the existing processor id for them.
1993 : //
1994 : // A node which is "new" on one processor will be associated with at
1995 : // least one ghost element, and we can just query that ghost
1996 : // element's owner to find out the correct processor id.
1997 :
1998 : auto node_unpartitioned =
1999 24034 : [](const Elem * elem, unsigned int local_node_num)
2000 24034 : { return elem->node_ref(local_node_num).processor_id() ==
2001 24034 : DofObject::invalid_processor_id; };
2002 :
2003 38 : SyncProcIds sync(mesh);
2004 :
2005 : sync_node_data_by_element_id_once
2006 66359 : (mesh, mesh.not_local_elements_begin(),
2007 22183 : mesh.not_local_elements_end(), Parallel::SyncEverything(),
2008 : node_unpartitioned, sync);
2009 :
2010 : // Nodes should now be unpartitioned iff they are truly new; those
2011 : // are the *only* nodes we will touch.
2012 : #ifdef DEBUG
2013 38 : MeshTools::libmesh_assert_parallel_consistent_new_node_procids(mesh);
2014 : #endif
2015 :
2016 22183 : NodeWasNew node_was_new(mesh);
2017 :
2018 : // Set the lowest processor id we can on truly new nodes
2019 9419710 : for (auto & elem : mesh.element_ptr_range())
2020 43950696 : for (auto & node : elem->node_ref_range())
2021 39251931 : if (node_was_new.was_new.count(&node))
2022 : {
2023 18292 : processor_id_type & pid = node.processor_id();
2024 29155819 : pid = std::min(pid, elem->processor_id());
2025 22069 : }
2026 :
2027 : // Then finally see if other processors have a lower option
2028 : Parallel::sync_node_data_by_element_id
2029 66359 : (mesh, mesh.elements_begin(), mesh.elements_end(),
2030 22107 : ElemNodesMaybeNew(), node_was_new, sync);
2031 :
2032 : // We should have consistent processor ids when we're done.
2033 : #ifdef DEBUG
2034 38 : MeshTools::libmesh_assert_parallel_consistent_procids<Node>(mesh);
2035 38 : MeshTools::libmesh_assert_parallel_consistent_new_node_procids(mesh);
2036 : #endif
2037 22145 : }
2038 :
2039 :
2040 :
2041 : // ------------------------------------------------------------
2042 8033 : void MeshCommunication::make_nodes_parallel_consistent (MeshBase & mesh)
2043 : {
2044 : // This function must be run on all processors at once
2045 10 : libmesh_parallel_only(mesh.comm());
2046 :
2047 : // When this function is called, each section of a parallelized mesh
2048 : // should be in the following state:
2049 : //
2050 : // Element ids and locations should be consistent on every processor
2051 : // where they exist.
2052 : //
2053 : // All nodes should have the exact same physical location on every
2054 : // processor where they exist.
2055 : //
2056 : // Local nodes should have unique authoritative ids,
2057 : // and processor ids consistent with all processors which own
2058 : // an element touching them.
2059 : //
2060 : // Ghost nodes touching local elements should have processor ids
2061 : // consistent with all processors which own an element touching
2062 : // them.
2063 : //
2064 : // Ghost nodes should have ids which are either already correct
2065 : // or which are in the "unpartitioned" id space.
2066 :
2067 : // First, let's sync up processor ids. Some of these processor ids
2068 : // may be "wrong" from coarsening, but they're right in the sense
2069 : // that they'll tell us who has the authoritative dofobject ids for
2070 : // each node.
2071 :
2072 8033 : this->make_node_proc_ids_parallel_consistent(mesh);
2073 :
2074 : // Second, sync up dofobject ids.
2075 8033 : this->make_node_ids_parallel_consistent(mesh);
2076 :
2077 : // Third, sync up dofobject unique_ids if applicable.
2078 8033 : this->make_node_unique_ids_parallel_consistent(mesh);
2079 :
2080 : // Finally, correct the processor ids to make DofMap happy
2081 8033 : MeshTools::correct_node_proc_ids(mesh);
2082 8033 : }
2083 :
2084 :
2085 :
2086 : // ------------------------------------------------------------
2087 22145 : void MeshCommunication::make_new_nodes_parallel_consistent (MeshBase & mesh)
2088 : {
2089 : // This function must be run on all processors at once
2090 38 : libmesh_parallel_only(mesh.comm());
2091 :
2092 : // When this function is called, each section of a parallelized mesh
2093 : // should be in the following state:
2094 : //
2095 : // All nodes should have the exact same physical location on every
2096 : // processor where they exist.
2097 : //
2098 : // Local nodes should have unique authoritative ids,
2099 : // and new nodes should be unpartitioned.
2100 : //
2101 : // New ghost nodes touching local elements should be unpartitioned.
2102 : //
2103 : // New ghost nodes should have ids which are either already correct
2104 : // or which are in the "unpartitioned" id space.
2105 : //
2106 : // Non-new nodes should have correct ids and processor ids already.
2107 :
2108 : // First, let's sync up new nodes' processor ids.
2109 :
2110 22145 : this->make_new_node_proc_ids_parallel_consistent(mesh);
2111 :
2112 : // Second, sync up dofobject ids.
2113 22145 : this->make_node_ids_parallel_consistent(mesh);
2114 :
2115 : // Third, sync up dofobject unique_ids if applicable.
2116 22145 : this->make_node_unique_ids_parallel_consistent(mesh);
2117 :
2118 : // Fourth, sync up any nodal boundary conditions
2119 22145 : this->make_node_bcids_parallel_consistent(mesh);
2120 :
2121 : // Finally, correct the processor ids to make DofMap happy
2122 22145 : MeshTools::correct_node_proc_ids(mesh);
2123 22145 : }
2124 :
2125 :
2126 :
2127 : // ------------------------------------------------------------
2128 : void
2129 399980 : MeshCommunication::delete_remote_elements (DistributedMesh & mesh,
2130 : const std::set<Elem *> & extra_ghost_elem_ids) const
2131 : {
2132 : // The mesh should know it's about to be parallelized
2133 382 : libmesh_assert (!mesh.is_serial());
2134 :
2135 764 : LOG_SCOPE("delete_remote_elements()", "MeshCommunication");
2136 :
2137 : #ifdef DEBUG
2138 : // We expect maximum ids to be in sync so we can use them to size
2139 : // vectors
2140 382 : libmesh_assert(mesh.comm().verify(mesh.max_node_id()));
2141 382 : libmesh_assert(mesh.comm().verify(mesh.max_elem_id()));
2142 382 : const dof_id_type par_max_node_id = mesh.parallel_max_node_id();
2143 382 : const dof_id_type par_max_elem_id = mesh.parallel_max_elem_id();
2144 382 : libmesh_assert_equal_to (par_max_node_id, mesh.max_node_id());
2145 382 : libmesh_assert_equal_to (par_max_elem_id, mesh.max_elem_id());
2146 382 : const dof_id_type n_constraint_rows = mesh.n_constraint_rows();
2147 : #endif
2148 :
2149 764 : connected_elem_set_type elements_to_keep;
2150 :
2151 : // Don't delete elements that we were explicitly told not to
2152 399980 : for (const auto & elem : extra_ghost_elem_ids)
2153 : {
2154 0 : std::vector<const Elem *> active_family;
2155 : #ifdef LIBMESH_ENABLE_AMR
2156 0 : if (!elem->subactive())
2157 0 : elem->active_family_tree(active_family);
2158 : else
2159 : #endif
2160 0 : active_family.push_back(elem);
2161 :
2162 0 : for (const auto & f : active_family)
2163 0 : elements_to_keep.insert(f);
2164 : }
2165 :
2166 : // See which elements we still need to keep ghosted, given that
2167 : // we're keeping local and unpartitioned elements.
2168 : query_ghosting_functors
2169 401508 : (mesh, mesh.processor_id(),
2170 800342 : mesh.active_pid_elements_begin(mesh.processor_id()),
2171 400744 : mesh.active_pid_elements_end(mesh.processor_id()),
2172 : elements_to_keep);
2173 : query_ghosting_functors
2174 401126 : (mesh, DofObject::invalid_processor_id,
2175 799960 : mesh.active_pid_elements_begin(DofObject::invalid_processor_id),
2176 799578 : mesh.active_pid_elements_end(DofObject::invalid_processor_id),
2177 : elements_to_keep);
2178 :
2179 : // The inactive elements we need to send should have their
2180 : // immediate children present.
2181 1199558 : connect_children(mesh, mesh.pid_elements_begin(mesh.processor_id()),
2182 400744 : mesh.pid_elements_end(mesh.processor_id()),
2183 : elements_to_keep);
2184 401126 : connect_children(mesh,
2185 799960 : mesh.pid_elements_begin(DofObject::invalid_processor_id),
2186 799960 : mesh.pid_elements_end(DofObject::invalid_processor_id),
2187 : elements_to_keep);
2188 :
2189 : // And see which elements and nodes they depend on
2190 764 : connected_node_set_type connected_nodes;
2191 399980 : connect_element_dependencies(mesh, elements_to_keep, connected_nodes);
2192 :
2193 : // Delete all the elements we have no reason to save,
2194 : // starting with the most refined so that the mesh
2195 : // is valid at all intermediate steps
2196 399980 : unsigned int n_levels = MeshTools::n_levels(mesh);
2197 :
2198 913877 : for (int l = n_levels - 1; l >= 0; --l)
2199 1046434 : for (auto & elem : as_range(mesh.level_elements_begin(l),
2200 100475114 : mesh.level_elements_end(l)))
2201 : {
2202 19116 : libmesh_assert (elem);
2203 : // Make sure we don't leave any invalid pointers
2204 55593 : const bool keep_me = elements_to_keep.count(elem);
2205 :
2206 19116 : if (!keep_me)
2207 35844140 : elem->make_links_to_me_remote();
2208 :
2209 : // delete_elem doesn't currently invalidate element
2210 : // iterators... that had better not change
2211 38232 : if (!keep_me)
2212 35844140 : mesh.delete_elem(elem);
2213 512945 : }
2214 :
2215 : // Delete all the nodes we have no reason to save
2216 165485436 : for (auto & node : mesh.node_ptr_range())
2217 : {
2218 51556 : libmesh_assert(node);
2219 151615 : if (!connected_nodes.count(node))
2220 : {
2221 3053 : libmesh_assert_not_equal_to(node->processor_id(),
2222 : mesh.processor_id());
2223 49936854 : mesh.delete_node(node);
2224 : }
2225 399216 : }
2226 :
2227 : // If we had a point locator, it's invalid now that some of the
2228 : // elements it pointed to have been deleted.
2229 399980 : mesh.clear_point_locator();
2230 :
2231 : // We now have all remote elements and nodes deleted; our ghosting
2232 : // functors should be ready to delete any now-redundant cached data
2233 : // they use too.
2234 932984 : for (auto & gf : as_range(mesh.ghosting_functors_begin(), mesh.ghosting_functors_end()))
2235 533004 : gf->delete_remote_elements();
2236 :
2237 : #ifdef DEBUG
2238 382 : const dof_id_type n_new_constraint_rows = mesh.n_constraint_rows();
2239 382 : libmesh_assert_equal_to(n_constraint_rows, n_new_constraint_rows);
2240 :
2241 382 : MeshTools::libmesh_assert_valid_refinement_tree(mesh);
2242 382 : MeshTools::libmesh_assert_valid_constraint_rows(mesh);
2243 : #endif
2244 399980 : }
2245 :
2246 : } // namespace libMesh
|
