libMesh
mesh_communication.C
Go to the documentation of this file.
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  SyncNeighbors(MeshBase & _mesh) :
60  mesh(_mesh) {}
61 
62  MeshBase & mesh;
63 
64  // Find the neighbor ids for each requested element
65  void gather_data (const std::vector<dof_id_type> & ids,
66  std::vector<datum> & neighbors) const
67  {
68  neighbors.resize(ids.size());
69 
70  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  const Elem & elem = mesh.elem_ref(ids[i]);
75 
76  // Return the element's neighbors
77  const unsigned int n_neigh = elem.n_neighbors();
78  neighbors[i].resize(n_neigh);
79  for (unsigned int n = 0; n != n_neigh; ++n)
80  {
81  const Elem * neigh = elem.neighbor_ptr(n);
82  if (neigh)
83  {
84  libmesh_assert_not_equal_to(neigh, remote_elem);
85  neighbors[i][n] = neigh->id();
86  }
87  else
88  neighbors[i][n] = DofObject::invalid_id;
89  }
90  }
91  }
92 
93  void act_on_data (const std::vector<dof_id_type> & ids,
94  const std::vector<datum> & neighbors) const
95  {
96  for (auto i : index_range(ids))
97  {
98  Elem & elem = mesh.elem_ref(ids[i]);
99 
100  const datum & new_neigh = neighbors[i];
101 
102  const unsigned int n_neigh = elem.n_neighbors();
103  libmesh_assert_equal_to (n_neigh, new_neigh.size());
104 
105  for (unsigned int n = 0; n != n_neigh; ++n)
106  {
107  const dof_id_type new_neigh_id = new_neigh[n];
108  const Elem * old_neigh = elem.neighbor_ptr(n);
109  if (old_neigh && old_neigh != remote_elem)
110  {
111  libmesh_assert_equal_to(old_neigh->id(), new_neigh_id);
112  }
113  else if (new_neigh_id == DofObject::invalid_id)
114  {
115  libmesh_assert (!old_neigh);
116  }
117  else
118  {
119  Elem * neigh = mesh.query_elem_ptr(new_neigh_id);
120  if (neigh)
121  elem.set_neighbor(n, neigh);
122  else
123  elem.set_neighbor(n, const_cast<RemoteElem *>(remote_elem));
124  }
125  }
126  }
127  }
128 };
129 
130 
131 void
132 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  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  const auto & constraint_rows = mesh->get_constraint_rows();
144 
145  std::unordered_set<const Elem *> constraining_nodes_elems;
146  for (const Elem * elem : connected_elements)
147  {
148  for (const Node & node : elem->node_ref_range())
149  {
150  // Retain all elements containing constraining nodes
151  if (const auto it = constraint_rows.find(&node);
152  it != constraint_rows.end())
153  for (auto & p : it->second)
154  {
155  const Elem * constraining_elem = p.first.first;
156  libmesh_assert(constraining_elem ==
157  mesh->elem_ptr(constraining_elem->id()));
158  if (!connected_elements.count(constraining_elem) &&
159  !new_connected_elements.count(constraining_elem))
160  constraining_nodes_elems.insert(constraining_elem);
161  }
162  }
163  }
164 
165  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  elem_rit = new_connected_elements.rbegin();
188 
189  for (; elem_rit != new_connected_elements.rend(); ++elem_rit)
190  {
191  const Elem * elem = *elem_rit;
192  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  for (const Elem * parent = elem->parent(); parent;
198  parent = parent->parent())
199  if (!connected_elements.count(parent) &&
200  !new_connected_elements.count(parent))
201  newer_connected_elements.insert (parent);
202 
203  auto total_family_insert =
204  [&connected_elements, &new_connected_elements,
205  &newer_connected_elements]
206  (const Elem * e)
207  {
208  if (e->active() && e->has_children())
209  {
210  std::vector<const Elem *> subactive_family;
211  e->total_family_tree(subactive_family);
212  for (const auto & f : subactive_family)
213  {
215  if (!connected_elements.count(f) &&
216  !new_connected_elements.count(f))
217  newer_connected_elements.insert(f);
218  }
219  }
220  };
221 
222  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  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  libmesh_assert(!interior_parent || mesh);
233 
234  if (interior_parent &&
235  interior_parent == mesh->query_elem_ptr(interior_parent->id()) &&
236  !connected_elements.count(interior_parent) &&
237  !new_connected_elements.count(interior_parent))
238  {
239  newer_connected_elements.insert (interior_parent);
240  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  (const Elem * elem)
254  {
255  libmesh_assert(elem);
256  const Elem * parent = elem->parent();
257  if (parent)
258  libmesh_assert(connected_elements.count(parent) ||
259  new_connected_elements.count(parent) ||
260  newer_connected_elements.count(parent));
261  };
262 
263  for (const auto & elem : connected_elements)
264  check_elem(elem);
265  for (const auto & elem : new_connected_elements)
266  check_elem(elem);
267  for (const auto & elem : newer_connected_elements)
268  check_elem(elem);
269 # endif // DEBUG
270 
271 #endif // LIBMESH_ENABLE_AMR
272 }
273 
274 
275 
276 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  for (const auto & elem : new_connected_elements)
282  for (auto & n : elem->node_ref_range())
283  if (!connected_nodes.count(&n) &&
284  !new_connected_nodes.count(&n))
285  newer_connected_nodes.insert(&n);
286 }
287 
288 
289 } // anonymous namespace
290 
291 
292 
293 namespace libMesh
294 {
295 
296 
298  processor_id_type pid,
301  connected_elem_set_type & connected_elements)
302 {
303  for (auto & gf :
306  {
307  GhostingFunctor::map_type elements_to_ghost;
308  libmesh_assert(gf);
309  (*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  for (auto & pr : elements_to_ghost)
314  {
315  const Elem * elem = pr.first;
316  libmesh_assert(elem != remote_elem);
317  libmesh_assert(mesh.elem_ptr(elem->id()) == elem);
318  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  for (; elem_it != elem_end; ++elem_it)
325  connected_elements.insert(*elem_it);
326 }
327 
328 
332  connected_elem_set_type & connected_elements)
333 {
334  // None of these parameters are used when !LIBMESH_ENABLE_AMR.
335  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  for (const auto & elem : as_range(elem_it, elem_end))
342  {
343  if (elem->has_children())
344  for (auto & child : elem->child_ref_range())
345  if (&child != remote_elem)
346  connected_elements.insert(&child);
347  }
348 #endif // LIBMESH_ENABLE_AMR
349 }
350 
351 
352 #ifdef LIBMESH_ENABLE_DEPRECATED
353 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  connect_element_families(connected_elements, connected_elements,
364  connected_elements, mesh);
365 }
366 #endif // LIBMESH_ENABLE_DEPRECATED
367 
368 
369 
370 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  connected_nodes.clear();
376 
377  // Use the newer API
378  connect_nodes(connected_elements, connected_nodes, connected_nodes,
379  connected_nodes);
380 }
381 
382 
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  connected_elem_set_type new_connected_elements;
390  connected_node_set_type new_connected_nodes;
391  new_connected_elements.swap(connected_elements);
392  new_connected_nodes.swap(connected_nodes);
393 
394  while (!new_connected_elements.empty() ||
395  !new_connected_nodes.empty())
396  {
397  auto [newer_connected_elements,
398  newer_connected_nodes] =
400  (mesh, connected_elements, connected_nodes,
401  new_connected_elements, new_connected_nodes);
402 
403  // These have now been examined
404  connected_elements.merge(new_connected_elements);
405  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  libmesh_assert(new_connected_elements.empty());
410  libmesh_assert(new_connected_nodes.empty());
411 
412  // These now need to be examined
413  new_connected_elements.swap(newer_connected_elements);
414  new_connected_nodes.swap(newer_connected_nodes);
415  }
416 }
417 
418 
419 std::pair<connected_elem_set_type, connected_node_set_type>
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  std::pair<connected_elem_set_type, connected_node_set_type> returnval;
427  auto & [newer_connected_elements, newer_connected_nodes] = returnval;
428  connect_element_families(connected_elements, new_connected_elements,
429  newer_connected_elements, &mesh);
430 
431  connect_nodes(new_connected_elements, connected_nodes,
432  new_connected_nodes, newer_connected_nodes);
433 
434  return returnval;
435 }
436 
437 
438 
439 
440 // ------------------------------------------------------------
441 // MeshCommunication class members
443 {
444  // _neighboring_processors.clear();
445 }
446 
447 
448 
449 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
450 // ------------------------------------------------------------
452 {
453  // no MPI == one processor, no redistribution
454  return;
455 }
456 
457 #else
458 // ------------------------------------------------------------
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  libmesh_parallel_only(mesh.comm());
486  libmesh_assert (MeshTools::n_elem(mesh.unpartitioned_elements_begin(),
487  mesh.unpartitioned_elements_end()) == 0);
488 
489  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  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  newly_coarsened_only ?
501  mesh.flagged_elements_begin(Elem::JUST_COARSENED) :
502 #endif
503  mesh.active_elements_begin();
504 
505  const MeshBase::const_element_iterator send_elems_end =
506 #ifdef LIBMESH_ENABLE_AMR
507  newly_coarsened_only ?
508  mesh.flagged_elements_end(Elem::JUST_COARSENED) :
509 #endif
510  mesh.active_elements_end();
511 
512  // See what should get sent where. We don't send to ourselves.
513  for (auto & elem : as_range(send_elems_begin, send_elems_end))
514  if (elem->processor_id() != mesh.processor_id())
515  send_to_pid[elem->processor_id()].push_back(elem);
516 
517  std::map<processor_id_type, std::vector<const Node *>> all_nodes_to_send;
518  std::map<processor_id_type, std::vector<const Elem *>> all_elems_to_send;
519 
520  // We may need to send constraint rows too.
521  auto & constraint_rows = mesh.get_constraint_rows();
522  bool have_constraint_rows = !constraint_rows.empty();
523  mesh.comm().broadcast(have_constraint_rows);
524 
525 #ifdef DEBUG
526  const dof_id_type n_constraint_rows =
527  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  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  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  libmesh_assert(!p_elements.empty());
548 
549  // Be compatible with both deprecated and
550  // corrected MeshBase iterator types
552 
553  v_t * elempp = p_elements.data();
554  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 
564  (elempp, elemend, Predicates::NotNull<v_t *>());
565 
566  const MeshBase::const_element_iterator elem_end =
568  (elemend, elemend, Predicates::NotNull<v_t *>());
569 
570  connected_elem_set_type elements_to_send;
571 
572  // See which to-be-ghosted elements we need to send
573  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  connect_children(mesh, mesh.pid_elements_begin(pid),
579  mesh.pid_elements_end(pid),
580  elements_to_send);
581 
582  // Now see which other elements and nodes they depend on
583  connected_node_set_type connected_nodes;
584  connect_element_dependencies(mesh, elements_to_send,
585  connected_nodes);
586 
587  all_nodes_to_send[pid].assign(connected_nodes.begin(),
588  connected_nodes.end());
589 
590  all_elems_to_send[pid].assign(elements_to_send.begin(),
591  elements_to_send.end());
592 
593  for (auto & [node, row] : constraint_rows)
594  {
595  if (!connected_nodes.count(node))
596  continue;
597 
598  serialized_row_type serialized_row;
599  for (auto [elem_and_node, coef] : row)
600  serialized_row.emplace_back(std::make_pair(elem_and_node.first->id(),
601  elem_and_node.second),
602  coef);
603 
604  all_constraint_rows_to_send[pid].emplace_back
605  (node->id(), std::move(serialized_row));
606  }
607  }
608 
609  // Elem/Node unpack() automatically adds them to the given mesh
610  auto null_node_action = [](processor_id_type, const std::vector<const Node*>&){};
611  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  TIMPI::push_parallel_packed_range(mesh.comm(), all_nodes_to_send, &mesh,
615  null_node_action);
616 
617  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  if (have_constraint_rows)
623  {
624  auto constraint_row_action =
625  [&mesh, &constraint_rows]
626  (processor_id_type /* src_pid */,
627  const std::vector<std::pair<dof_id_type, serialized_row_type>> rows)
628  {
629  for (auto & [node_id, serialized_row] : rows)
630  {
632  for (auto [elem_and_node, coef] : serialized_row)
633  row.emplace_back(std::make_pair(mesh.elem_ptr(elem_and_node.first),
634  elem_and_node.second),
635  coef);
636 
637  constraint_rows[mesh.node_ptr(node_id)] = row;
638  }
639  };
640 
642  all_constraint_rows_to_send,
643  constraint_row_action);
644 
645  }
646 
647  // Check on the redistribution consistency
648 #ifdef DEBUG
650 
652 
653  const dof_id_type new_n_constraint_rows =
654  have_constraint_rows ? mesh.n_constraint_rows() : 0;
655 
656  libmesh_assert_equal_to(n_constraint_rows, new_n_constraint_rows);
657 
659 #endif
660 
661  // If we had a point locator, it's invalid now that there are new
662  // elements it can't locate.
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  mesh.MeshBase::redistribute();
669 }
670 #endif // LIBMESH_HAVE_MPI
671 
672 
673 
674 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
675 // ------------------------------------------------------------
677 {
678  // no MPI == one processor, no need for this method...
679  return;
680 }
681 #else
682 // ------------------------------------------------------------
684 {
685  // Don't need to do anything if there is
686  // only one processor.
687  if (mesh.n_processors() == 1)
688  return;
689 
690  // This function must be run on all processors at once
691  libmesh_parallel_only(mesh.comm());
692 
693  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  mesh.find_neighbors (/* reset_remote_elements = */ true,
724  /* reset_current_list = */ true);
725 
726  // Get a unique message tag to use in communications
728  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  std::vector<processor_id_type> adjacent_processors;
740  for (auto pid : make_range(mesh.n_processors()))
741  if (pid != mesh.processor_id())
742  adjacent_processors.push_back (pid);
743 
744 
745  const processor_id_type n_adjacent_processors =
746  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  std::vector<dof_id_type> my_interface_node_list;
753  std::vector<const Elem *> my_interface_elements;
754  {
755  std::set<dof_id_type> my_interface_node_set;
756 
757  // For avoiding extraneous element side construction
758  ElemSideBuilder side_builder;
759 
760  // since parent nodes are a subset of children nodes, this should be sufficient
761  for (const auto & elem : mesh.active_local_element_ptr_range())
762  {
763  libmesh_assert(elem);
764 
765  if (elem->on_boundary()) // denotes *any* side has a nullptr neighbor
766  {
767  my_interface_elements.push_back(elem); // add the element, but only once, even
768  // if there are multiple nullptr neighbors
769  for (auto s : elem->side_index_range())
770  if (elem->neighbor_ptr(s) == nullptr)
771  {
772  const Elem & side = side_builder(*elem, s);
773 
774  for (auto n : make_range(side.n_vertices()))
775  my_interface_node_set.insert (side.node_id(n));
776  }
777  }
778  }
779 
780  my_interface_node_list.reserve (my_interface_node_set.size());
781  my_interface_node_list.insert (my_interface_node_list.end(),
782  my_interface_node_set.begin(),
783  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  my_interface_node_xfer_buffers (n_adjacent_processors, my_interface_node_list);
795  std::map<processor_id_type, unsigned char> n_comm_steps;
796 
797  std::vector<Parallel::Request> send_requests (3*n_adjacent_processors);
798  unsigned int current_request = 0;
799 
800  for (unsigned int comm_step=0; comm_step<n_adjacent_processors; comm_step++)
801  {
802  n_comm_steps[adjacent_processors[comm_step]]=1;
803  mesh.comm().send (adjacent_processors[comm_step],
804  my_interface_node_xfer_buffers[comm_step],
805  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  adjacent_processors.clear();
818 
819  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  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  status(mesh.comm().probe (Parallel::any_source,
840  element_neighbors_tag));
841  const processor_id_type
842  source_pid_idx = cast_int<processor_id_type>(status.source()),
843  dest_pid_idx = source_pid_idx;
844 
845  //------------------------------------------------------------------
846  // first time - incoming request
847  if (n_comm_steps[source_pid_idx] == 1)
848  {
849  n_comm_steps[source_pid_idx]++;
850 
851  mesh.comm().receive (source_pid_idx,
852  common_interface_node_list,
853  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  (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  common_interface_node_list.begin()),
871  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  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  connected_node_set_type connected_nodes;
897 
898  // Check for quick return?
899  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  mesh.comm().send_packed_range (dest_pid_idx,
908  &mesh,
909  connected_nodes.begin(),
910  connected_nodes.end(),
911  send_requests[current_request++],
912  element_neighbors_tag);
913 
914  mesh.comm().send_packed_range (dest_pid_idx,
915  &mesh,
916  elements_to_send.begin(),
917  elements_to_send.end(),
918  send_requests[current_request++],
919  element_neighbors_tag);
920 
921  continue;
922  }
923  // otherwise, this really *is* an adjacent processor.
924  adjacent_processors.push_back(source_pid_idx);
925 
926  std::vector<const Elem *> family_tree;
927 
928  for (auto & elem : my_interface_elements)
929  {
930  std::size_t n_shared_nodes = 0;
931 
932  for (auto n : make_range(elem->n_vertices()))
933  if (std::binary_search (common_interface_node_list.begin(),
934  common_interface_node_list.end(),
935  elem->node_id(n)))
936  {
937  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  if (n_shared_nodes > 0) break;
943  }
944 
945  if (n_shared_nodes) // share at least one node?
946  {
947  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  if (!elements_to_send.count(elem))
952  {
953 #ifdef LIBMESH_ENABLE_AMR
954  elem->family_tree(family_tree);
955 #else
956  family_tree.clear();
957  family_tree.push_back(elem);
958 #endif
959  for (const auto & f : family_tree)
960  {
961  elem = f;
962  elements_to_send.insert (elem);
963 
964  for (auto & n : elem->node_ref_range())
965  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  libmesh_assert (connected_nodes.empty() || !elements_to_send.empty());
977  libmesh_assert (!connected_nodes.empty() || elements_to_send.empty());
978 
979  // send the nodes off to the destination processor
980  mesh.comm().send_packed_range (dest_pid_idx,
981  &mesh,
982  connected_nodes.begin(),
983  connected_nodes.end(),
984  send_requests[current_request++],
985  element_neighbors_tag);
986 
987  // send the elements off to the destination processor
988  mesh.comm().send_packed_range (dest_pid_idx,
989  &mesh,
990  elements_to_send.begin(),
991  elements_to_send.end(),
992  send_requests[current_request++],
993  element_neighbors_tag);
994  }
995  }
996  //------------------------------------------------------------------
997  // second time - reply of nodes
998  else if (n_comm_steps[source_pid_idx] == 2)
999  {
1000  n_comm_steps[source_pid_idx]++;
1001 
1002  mesh.comm().receive_packed_range (source_pid_idx,
1003  &mesh,
1005  (Node**)nullptr,
1006  element_neighbors_tag);
1007  }
1008  //------------------------------------------------------------------
1009  // third time - reply of elements
1010  else if (n_comm_steps[source_pid_idx] == 3)
1011  {
1012  n_comm_steps[source_pid_idx]++;
1013 
1014  mesh.comm().receive_packed_range (source_pid_idx,
1015  &mesh,
1017  (Elem**)nullptr,
1018  element_neighbors_tag);
1019  }
1020  //------------------------------------------------------------------
1021  // fourth time - shouldn't happen
1022  else
1023  {
1024  libMesh::err << "ERROR: unexpected number of replies: "
1025  << n_comm_steps[source_pid_idx]
1026  << std::endl;
1027  }
1028  } // done catching & processing replies associated with tag ~ 100,000pi
1029 
1030  // allow any pending requests to complete
1031  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.
1036 
1037  // We can now find neighbor information for the interfaces between
1038  // local elements and ghost elements.
1039  mesh.find_neighbors (/* reset_remote_elements = */ true,
1040  /* 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  SyncNeighbors nsync(mesh);
1046 
1048  (mesh.comm(), mesh.elements_begin(), mesh.elements_end(), nsync);
1049 }
1050 #endif // LIBMESH_HAVE_MPI
1051 
1052 
1053 #ifndef LIBMESH_HAVE_MPI // avoid spurious gcc warnings
1054 // ------------------------------------------------------------
1056 {
1057  // no MPI == one processor, no need for this method...
1058  return;
1059 }
1060 #else
1062 {
1063 
1064  // Don't need to do anything if all processors already ghost all non-local
1065  // elements.
1066  if (mesh.is_serial())
1067  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  ghost_map coarsening_elements_to_ghost;
1084 
1085  const processor_id_type proc_id = mesh.processor_id();
1086  // Look for just-coarsened elements
1087  for (auto elem : as_range(mesh.flagged_pid_elements_begin(Elem::COARSEN, proc_id),
1088  mesh.flagged_pid_elements_end(Elem::COARSEN, proc_id)))
1089  {
1090  // If it's flagged for coarsening it had better have a parent
1091  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  const processor_id_type their_proc_id = elem->parent()->processor_id();
1098  if (their_proc_id != proc_id)
1099  coarsening_elements_to_ghost[their_proc_id].push_back(elem);
1100  }
1101 
1102  std::map<processor_id_type, std::vector<const Node *>> all_nodes_to_send;
1103  std::map<processor_id_type, std::vector<const Elem *>> all_elems_to_send;
1104 
1105  const processor_id_type n_proc = mesh.n_processors();
1106 
1107  for (processor_id_type p=0; p != n_proc; ++p)
1108  {
1109  if (p == proc_id)
1110  continue;
1111 
1112  connected_elem_set_type elements_to_send;
1113  std::set<const Node *> nodes_to_send;
1114 
1115  if (const auto it = std::as_const(coarsening_elements_to_ghost).find(p);
1116  it != coarsening_elements_to_ghost.end())
1117  {
1118  const std::vector<Elem *> & elems = it->second;
1119  libmesh_assert(elems.size());
1120 
1121  // Make some fake element iterators defining this vector of
1122  // elements
1123  Elem * const * elempp = const_cast<Elem * const *>(elems.data());
1124  Elem * const * elemend = elempp+elems.size();
1125  const MeshBase::const_element_iterator elem_it =
1127  const MeshBase::const_element_iterator elem_end =
1129 
1130  for (auto & gf : as_range(mesh.ghosting_functors_begin(),
1132  {
1133  GhostingFunctor::map_type elements_to_ghost;
1134  libmesh_assert(gf);
1135  (*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  for (auto & pr : elements_to_ghost)
1140  {
1141  const Elem * elem = pr.first;
1142  libmesh_assert(elem);
1143  while (elem)
1144  {
1145  libmesh_assert(elem != remote_elem);
1146  elements_to_send.insert(elem);
1147  for (auto & n : elem->node_ref_range())
1148  nodes_to_send.insert(&n);
1149  elem = elem->parent();
1150  }
1151  }
1152  }
1153 
1154  all_nodes_to_send[p].assign(nodes_to_send.begin(), nodes_to_send.end());
1155  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  auto null_node_action = [](processor_id_type, const std::vector<const Node*>&){};
1161  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  TIMPI::push_parallel_packed_range(mesh.comm(), all_nodes_to_send, &mesh,
1165  null_node_action);
1166 
1167  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 // ------------------------------------------------------------
1178 {
1179  // no MPI == one processor, no need for this method...
1180  return;
1181 }
1182 #else
1183 // ------------------------------------------------------------
1185 {
1186  // Don't need to do anything if there is
1187  // only one processor.
1188  if (mesh.n_processors() == 1)
1189  return;
1190 
1191  // This function must be run on all processors at once
1192  libmesh_parallel_only(mesh.comm());
1193 
1194  LOG_SCOPE("broadcast()", "MeshCommunication");
1195 
1196  // Explicitly clear the mesh on all but processor 0.
1197  if (mesh.processor_id() != 0)
1198  mesh.clear();
1199 
1200  // We may have set extra data only on processor 0 in a read()
1205 
1206  // We may have set mapping data only on processor 0 in a read()
1207  unsigned char map_type = mesh.default_mapping_type();
1208  unsigned char map_data = mesh.default_mapping_data();
1209  mesh.comm().broadcast(map_type);
1210  mesh.comm().broadcast(map_data);
1212  mesh.set_default_mapping_data(map_data);
1213 
1214  // Broadcast nodes
1216  mesh.nodes_begin(),
1217  mesh.nodes_end(),
1218  &mesh,
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  const unsigned int n_levels = MeshTools::paranoid_n_levels(mesh);
1229 
1230  for (unsigned int l=0; l != n_levels; ++l)
1232  mesh.level_elements_begin(l),
1233  mesh.level_elements_end(l),
1234  &mesh,
1236 
1237  // Make sure mesh_dimension and elem_dimensions are consistent.
1239 
1240  // Make sure mesh id counts are consistent.
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  auto & constraint_rows = mesh.get_constraint_rows();
1248  bool have_constraint_rows = !constraint_rows.empty();
1249  mesh.comm().broadcast(have_constraint_rows);
1250  if (have_constraint_rows)
1251  {
1252  std::map<dof_id_type,
1253  std::vector<std::tuple<dof_id_type, unsigned int, Real>>>
1254  serialized_rows;
1255 
1256  for (auto & row : constraint_rows)
1257  {
1258  const Node * node = row.first;
1259  const dof_id_type rowid = node->id();
1260  libmesh_assert(node == mesh.node_ptr(rowid));
1261 
1262  std::vector<std::tuple<dof_id_type, unsigned int, Real>>
1263  serialized_row;
1264  for (auto & entry : row.second)
1265  serialized_row.push_back
1266  (std::make_tuple(entry.first.first->id(),
1267  entry.first.second, entry.second));
1268 
1269  serialized_rows.emplace(rowid, std::move(serialized_row));
1270  }
1271 
1272  mesh.comm().broadcast(serialized_rows);
1273  if (mesh.processor_id() != 0)
1274  {
1275  constraint_rows.clear();
1276 
1277  for (auto & row : serialized_rows)
1278  {
1279  const dof_id_type rowid = row.first;
1280  const Node * node = mesh.node_ptr(rowid);
1281 
1282  std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>>
1283  deserialized_row;
1284  for (auto & entry : row.second)
1285  deserialized_row.push_back
1286  (std::make_pair(std::make_pair(mesh.elem_ptr(std::get<0>(entry)),
1287  std::get<1>(entry)),
1288  std::get<2>(entry)));
1289 
1290  constraint_rows.emplace(node, deserialized_row);
1291  }
1292  }
1293  }
1294 
1295  // Broadcast all of the named entity information
1299 
1300  // If we had a point locator, it's invalid now that there are new
1301  // elements it can't locate.
1303 
1306 
1307 #ifdef DEBUG
1308  MeshTools::libmesh_assert_valid_procids<Elem>(mesh);
1309  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 // ------------------------------------------------------------
1319 {
1320  // no MPI == one processor, no need for this method...
1321  return;
1322 }
1323 #else
1324 // ------------------------------------------------------------
1325 void MeshCommunication::gather (const processor_id_type root_id, MeshBase & mesh) const
1326 {
1327  // Check for quick return
1328  if (mesh.n_processors() == 1)
1329  return;
1330 
1331  // This function must be run on all processors at once
1332  libmesh_parallel_only(mesh.comm());
1333 
1334  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  approx_total_buffer_size / mesh.comm().size();
1340 
1341  (root_id == DofObject::invalid_processor_id) ?
1342 
1344  mesh.nodes_begin(),
1345  mesh.nodes_end(),
1347  approx_each_buffer_size) :
1348 
1349  mesh.comm().gather_packed_range (root_id,
1350  &mesh,
1351  mesh.nodes_begin(),
1352  mesh.nodes_end(),
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  const unsigned int n_levels = MeshTools::n_levels(mesh);
1359 
1360  for (unsigned int l=0; l != n_levels; ++l)
1361  (root_id == DofObject::invalid_processor_id) ?
1362 
1364  mesh.level_elements_begin(l),
1365  mesh.level_elements_end(l),
1367  approx_each_buffer_size) :
1368 
1369  mesh.comm().gather_packed_range (root_id,
1370  &mesh,
1371  mesh.level_elements_begin(l),
1372  mesh.level_elements_end(l),
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.
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  auto & constraint_rows = mesh.get_constraint_rows();
1385  bool have_constraint_rows = !constraint_rows.empty();
1386  mesh.comm().max(have_constraint_rows);
1387  if (have_constraint_rows)
1388  {
1389  std::map<dof_id_type,
1390  std::vector<std::tuple<dof_id_type, unsigned int, Real>>>
1391  serialized_rows;
1392 
1393  for (auto & row : constraint_rows)
1394  {
1395  const Node * node = row.first;
1396  const dof_id_type rowid = node->id();
1397  libmesh_assert(node == mesh.node_ptr(rowid));
1398 
1399  std::vector<std::tuple<dof_id_type, unsigned int, Real>>
1400  serialized_row;
1401  for (auto & entry : row.second)
1402  serialized_row.push_back
1403  (std::make_tuple(entry.first.first->id(),
1404  entry.first.second, entry.second));
1405 
1406  serialized_rows.emplace(rowid, std::move(serialized_row));
1407  }
1408 
1409  if (root_id == DofObject::invalid_processor_id)
1410  mesh.comm().set_union(serialized_rows);
1411  else
1412  mesh.comm().set_union(serialized_rows, root_id);
1413 
1414  if (root_id == DofObject::invalid_processor_id ||
1415  root_id == mesh.processor_id())
1416  {
1417  for (auto & row : serialized_rows)
1418  {
1419  const dof_id_type rowid = row.first;
1420  const Node * node = mesh.node_ptr(rowid);
1421 
1422  std::vector<std::pair<std::pair<const Elem *, unsigned int>, Real>>
1423  deserialized_row;
1424  for (auto & entry : row.second)
1425  deserialized_row.push_back
1426  (std::make_pair(std::make_pair(mesh.elem_ptr(std::get<0>(entry)),
1427  std::get<1>(entry)),
1428  std::get<2>(entry)));
1429 
1430  constraint_rows.emplace(node, deserialized_row);
1431  }
1432  }
1433 #ifdef DEBUG
1435 #endif
1436  }
1437 
1438 
1439  // If we are doing an allgather(), perform sanity check on the result.
1440  if (root_id == DofObject::invalid_processor_id)
1441  {
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  if (mesh.allow_find_neighbors())
1451  root_id == mesh.processor_id());
1452 
1453  // All done, but let's make sure it's done correctly
1454 
1455 #ifdef DEBUG
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  SyncIds(MeshBase & _mesh, renumber_obj _renumberer) :
1473  mesh(_mesh),
1474  renumber(_renumberer) {}
1475 
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  void gather_data (const std::vector<dof_id_type> & ids,
1483  std::vector<datum> & ids_out) const
1484  {
1485  ids_out = ids;
1486  }
1487 
1488  void act_on_data (const std::vector<dof_id_type> & old_ids,
1489  const std::vector<datum> & new_ids) const
1490  {
1491  for (auto i : index_range(old_ids))
1492  if (old_ids[i] != new_ids[i])
1493  (mesh.*renumber)(old_ids[i], new_ids[i]);
1494  }
1495 };
1496 
1497 
1498 struct SyncNodeIds
1499 {
1500  typedef dof_id_type datum;
1501 
1502  SyncNodeIds(MeshBase & _mesh) :
1503  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;
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  void gather_data (const std::vector<dof_id_type> & ids,
1525  std::vector<datum> & ids_out) const
1526  {
1527  ids_out.clear();
1528  ids_out.resize(ids.size(), DofObject::invalid_id);
1529 
1530  for (auto i : index_range(ids))
1531  {
1532  const dof_id_type id = ids[i];
1533  const Node * node = mesh.query_node_ptr(id);
1534  if (node && (node->processor_id() == mesh.processor_id() ||
1535  definitive_ids.count(node)))
1536  ids_out[i] = id;
1537  }
1538  }
1539 
1540  bool act_on_data (const std::vector<dof_id_type> & old_ids,
1541  const std::vector<datum> & new_ids)
1542  {
1543  bool data_changed = false;
1544  for (auto i : index_range(old_ids))
1545  {
1546  const dof_id_type new_id = new_ids[i];
1547 
1548  const dof_id_type old_id = old_ids[i];
1549 
1550  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  if (!node)
1556  {
1557  // But let's check anyway in debug mode
1558 #ifdef DEBUG
1560  (definitive_renumbering.count(old_id));
1561  libmesh_assert_equal_to
1562  (new_id, definitive_renumbering[old_id]);
1563 #endif
1564  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  if (new_id == DofObject::invalid_id)
1570  {
1571  // But we might have gotten a definitive id from a
1572  // different request
1573  if (!definitive_ids.count(mesh.node_ptr(old_id)))
1574  data_changed = true;
1575  }
1576  else
1577  {
1578  if (node->processor_id() != mesh.processor_id())
1579  definitive_ids.insert(node);
1580  if (old_id != new_id)
1581  {
1582 #ifdef DEBUG
1584  (!definitive_renumbering.count(old_id));
1585  definitive_renumbering[old_id] = new_id;
1586 #endif
1587  mesh.renumber_node(old_id, new_id);
1588  data_changed = true;
1589  }
1590  }
1591  }
1592  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  SyncPLevels(MeshBase & _mesh) :
1603  mesh(_mesh) {}
1604 
1605  MeshBase & mesh;
1606 
1607  // Find the p_level of each requested Elem
1608  void gather_data (const std::vector<dof_id_type> & ids,
1609  std::vector<datum> & ids_out) const
1610  {
1611  ids_out.reserve(ids.size());
1612 
1613  for (const auto & id : ids)
1614  {
1615  Elem & elem = mesh.elem_ref(id);
1616  ids_out.push_back
1617  (std::make_pair(cast_int<unsigned char>(elem.p_level()),
1618  static_cast<unsigned char>(elem.p_refinement_flag())));
1619  }
1620  }
1621 
1622  void act_on_data (const std::vector<dof_id_type> & old_ids,
1623  const std::vector<datum> & new_p_levels) const
1624  {
1625  for (auto i : index_range(old_ids))
1626  {
1627  Elem & elem = mesh.elem_ref(old_ids[i]);
1628  // Make sure these are consistent
1630  (new_p_levels[i].first,
1631  static_cast<Elem::RefinementState>(new_p_levels[i].second));
1632  // Make sure parents' levels are consistent
1633  elem.set_p_level(new_p_levels[i].first);
1634  }
1635  }
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  SyncUniqueIds(MeshBase &_mesh, query_obj _querier) :
1648  mesh(_mesh),
1649  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  void gather_data (const std::vector<dof_id_type> & ids,
1657  std::vector<datum> & ids_out) const
1658  {
1659  ids_out.reserve(ids.size());
1660 
1661  for (const auto & id : ids)
1662  {
1663  DofObjSubclass * d = (mesh.*query)(id);
1664  libmesh_assert(d);
1665  ids_out.push_back(d->unique_id());
1666  }
1667  }
1668 
1669  void act_on_data (const std::vector<dof_id_type> & ids,
1670  const std::vector<datum> & unique_ids) const
1671  {
1672  for (auto i : index_range(ids))
1673  {
1674  DofObjSubclass * d = (mesh.*query)(ids[i]);
1675  libmesh_assert(d);
1676  d->set_unique_id(unique_ids[i]);
1677  }
1678  }
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  SyncBCIds(MeshBase &_mesh) :
1688  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  void gather_data (const std::vector<dof_id_type> & ids,
1695  std::vector<datum> & ids_out) const
1696  {
1697  ids_out.reserve(ids.size());
1698 
1699  const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1700 
1701  for (const auto & id : ids)
1702  {
1703  Node * n = mesh.query_node_ptr(id);
1704  libmesh_assert(n);
1705  std::vector<boundary_id_type> bcids;
1706  boundary_info.boundary_ids(n, bcids);
1707  ids_out.push_back(std::move(bcids));
1708  }
1709  }
1710 
1711  void act_on_data (const std::vector<dof_id_type> & ids,
1712  const std::vector<datum> & bcids) const
1713  {
1714  BoundaryInfo & boundary_info = mesh.get_boundary_info();
1715 
1716  for (auto i : index_range(ids))
1717  {
1718  Node * n = mesh.query_node_ptr(ids[i]);
1719  libmesh_assert(n);
1720  boundary_info.add_node(n, bcids[i]);
1721  }
1722  }
1723 };
1724 
1725 }
1726 
1727 
1728 
1729 // ------------------------------------------------------------
1731 {
1732  // This function must be run on all processors at once
1733  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  LOG_SCOPE ("make_node_ids_parallel_consistent()", "MeshCommunication");
1744 
1745  SyncNodeIds syncids(mesh);
1747  (mesh, mesh.elements_begin(), mesh.elements_end(),
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
1754 #endif
1755 }
1756 
1757 
1758 
1760 {
1761  // Avoid unused variable warnings if unique ids aren't enabled.
1763 
1764  // This function must be run on all processors at once
1765  libmesh_parallel_only(mesh.comm());
1766 
1767 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1768  LOG_SCOPE ("make_node_unique_ids_parallel_consistent()", "MeshCommunication");
1769 
1770  SyncUniqueIds<Node> syncuniqueids(mesh, &MeshBase::query_node_ptr);
1772  mesh.nodes_begin(),
1773  mesh.nodes_end(),
1774  syncuniqueids);
1775 
1776 #endif
1777 }
1778 
1779 
1781 {
1782  // Avoid unused variable warnings if unique ids aren't enabled.
1784 
1785  // This function must be run on all processors at once
1786  libmesh_parallel_only(mesh.comm());
1787 
1788  LOG_SCOPE ("make_node_bcids_parallel_consistent()", "MeshCommunication");
1789 
1790  SyncBCIds<Node> syncbcids(mesh);
1792  mesh.nodes_begin(),
1793  mesh.nodes_end(),
1794  syncbcids);
1795 }
1796 
1797 
1798 
1799 
1800 
1801 // ------------------------------------------------------------
1803 {
1804  // This function must be run on all processors at once
1805  libmesh_parallel_only(mesh.comm());
1806 
1807  LOG_SCOPE ("make_elems_parallel_consistent()", "MeshCommunication");
1808 
1809  SyncIds syncids(mesh, &MeshBase::renumber_elem);
1811  (mesh, mesh.active_elements_begin(),
1812  mesh.active_elements_end(), syncids);
1813 
1814 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1815  SyncUniqueIds<Elem> syncuniqueids(mesh, &MeshBase::query_elem_ptr);
1817  (mesh.comm(), mesh.active_elements_begin(),
1818  mesh.active_elements_end(), syncuniqueids);
1819 #endif
1820 }
1821 
1822 
1823 
1824 // ------------------------------------------------------------
1825 #ifdef LIBMESH_ENABLE_AMR
1827 {
1828  // This function must be run on all processors at once
1829  libmesh_parallel_only(mesh.comm());
1830 
1831  LOG_SCOPE ("make_p_levels_parallel_consistent()", "MeshCommunication");
1832 
1833  SyncPLevels syncplevels(mesh);
1835  (mesh.comm(), mesh.elements_begin(), mesh.elements_end(),
1836  syncplevels);
1837 }
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  SyncProcIds(MeshBase & _mesh) : mesh(_mesh) {}
1850 
1851  MeshBase & mesh;
1852 
1853  // ------------------------------------------------------------
1854  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  data.resize(ids.size());
1859 
1860  for (auto i : index_range(ids))
1861  {
1862  // Look for this point in the mesh
1863  if (ids[i] != DofObject::invalid_id)
1864  {
1865  Node & node = mesh.node_ref(ids[i]);
1866 
1867  // Return the node's correct processor id,
1868  data[i] = node.processor_id();
1869  }
1870  else
1872  }
1873  }
1874 
1875  // ------------------------------------------------------------
1876  bool act_on_data (const std::vector<dof_id_type> & ids,
1877  const std::vector<datum> proc_ids)
1878  {
1879  bool data_changed = false;
1880 
1881  // Set the ghost node processor ids we've now been informed of
1882  for (auto i : index_range(ids))
1883  {
1884  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  if (node.processor_id() > proc_ids[i])
1894  {
1895  data_changed = true;
1896  node.processor_id() = proc_ids[i];
1897  }
1898  }
1899 
1900  return data_changed;
1901  }
1902 };
1903 
1904 
1905 struct ElemNodesMaybeNew
1906 {
1907  ElemNodesMaybeNew() {}
1908 
1909  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  if (elem->refinement_flag() == Elem::JUST_REFINED)
1915  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  for (auto neigh : elem->neighbor_ptr_range())
1922  if (neigh == remote_elem)
1923  return true;
1924  return false;
1925  }
1926 };
1927 
1928 
1929 struct NodeWasNew
1930 {
1931  NodeWasNew(const MeshBase & mesh)
1932  {
1933  for (const auto & node : mesh.node_ptr_range())
1935  was_new.insert(node);
1936  }
1937 
1938  bool operator() (const Elem * elem, unsigned int local_node_num) const
1939  {
1940  if (was_new.count(elem->node_ptr(local_node_num)))
1941  return true;
1942  return false;
1943  }
1944 
1945  std::unordered_set<const Node *> was_new;
1946 };
1947 
1948 }
1949 
1950 
1951 
1952 // ------------------------------------------------------------
1954 {
1955  LOG_SCOPE ("make_node_proc_ids_parallel_consistent()", "MeshCommunication");
1956 
1957  // This function must be run on all processors at once
1958  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  SyncProcIds sync(mesh);
1975  (mesh, mesh.elements_begin(), mesh.elements_end(),
1977 }
1978 
1979 
1980 
1981 // ------------------------------------------------------------
1983 {
1984  LOG_SCOPE ("make_new_node_proc_ids_parallel_consistent()", "MeshCommunication");
1985 
1986  // This function must be run on all processors at once
1987  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
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  [](const Elem * elem, unsigned int local_node_num)
2018  { return elem->node_ref(local_node_num).processor_id() ==
2020 
2021  SyncProcIds sync(mesh);
2022 
2024  (mesh, mesh.not_local_elements_begin(),
2025  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
2032 #endif
2033 
2034  NodeWasNew node_was_new(mesh);
2035 
2036  // Set the lowest processor id we can on truly new nodes
2037  for (auto & elem : mesh.element_ptr_range())
2038  for (auto & node : elem->node_ref_range())
2039  if (node_was_new.was_new.count(&node))
2040  {
2041  processor_id_type & pid = node.processor_id();
2042  pid = std::min(pid, elem->processor_id());
2043  }
2044 
2045  // Then finally see if other processors have a lower option
2047  (mesh, mesh.elements_begin(), mesh.elements_end(),
2048  ElemNodesMaybeNew(), node_was_new, sync);
2049 
2050  // We should have consistent processor ids when we're done.
2051 #ifdef DEBUG
2054 #endif
2055 }
2056 
2057 
2058 
2059 // ------------------------------------------------------------
2061 {
2062  // This function must be run on all processors at once
2063  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 
2091 
2092  // Second, sync up dofobject ids.
2094 
2095  // Third, sync up dofobject unique_ids if applicable.
2097 
2098  // Finally, correct the processor ids to make DofMap happy
2100 }
2101 
2102 
2103 
2104 // ------------------------------------------------------------
2106 {
2107  // This function must be run on all processors at once
2108  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 
2129 
2130  // Second, sync up dofobject ids.
2132 
2133  // Third, sync up dofobject unique_ids if applicable.
2135 
2136  // Fourth, sync up any nodal boundary conditions
2138 
2139  // Finally, correct the processor ids to make DofMap happy
2141 }
2142 
2143 
2144 
2145 // ------------------------------------------------------------
2146 void
2148  const std::set<Elem *> & extra_ghost_elem_ids) const
2149 {
2150  // The mesh should know it's about to be parallelized
2152 
2153  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
2160  const dof_id_type par_max_node_id = mesh.parallel_max_node_id();
2161  const dof_id_type par_max_elem_id = mesh.parallel_max_elem_id();
2162  libmesh_assert_equal_to (par_max_node_id, mesh.max_node_id());
2163  libmesh_assert_equal_to (par_max_elem_id, mesh.max_elem_id());
2164  const dof_id_type n_constraint_rows = mesh.n_constraint_rows();
2165 #endif
2166 
2167  connected_elem_set_type elements_to_keep;
2168 
2169  // Don't delete elements that we were explicitly told not to
2170  for (const auto & elem : extra_ghost_elem_ids)
2171  {
2172  std::vector<const Elem *> active_family;
2173 #ifdef LIBMESH_ENABLE_AMR
2174  if (!elem->subactive())
2175  elem->active_family_tree(active_family);
2176  else
2177 #endif
2178  active_family.push_back(elem);
2179 
2180  for (const auto & f : active_family)
2181  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.
2187  (mesh, mesh.processor_id(),
2188  mesh.active_pid_elements_begin(mesh.processor_id()),
2189  mesh.active_pid_elements_end(mesh.processor_id()),
2190  elements_to_keep);
2193  mesh.active_pid_elements_begin(DofObject::invalid_processor_id),
2194  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  connect_children(mesh, mesh.pid_elements_begin(mesh.processor_id()),
2200  mesh.pid_elements_end(mesh.processor_id()),
2201  elements_to_keep);
2203  mesh.pid_elements_begin(DofObject::invalid_processor_id),
2204  mesh.pid_elements_end(DofObject::invalid_processor_id),
2205  elements_to_keep);
2206 
2207  // And see which elements and nodes they depend on
2208  connected_node_set_type connected_nodes;
2209  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  unsigned int n_levels = MeshTools::n_levels(mesh);
2215 
2216  for (int l = n_levels - 1; l >= 0; --l)
2217  for (auto & elem : as_range(mesh.level_elements_begin(l),
2218  mesh.level_elements_end(l)))
2219  {
2220  libmesh_assert (elem);
2221  // Make sure we don't leave any invalid pointers
2222  const bool keep_me = elements_to_keep.count(elem);
2223 
2224  if (!keep_me)
2225  elem->make_links_to_me_remote();
2226 
2227  // delete_elem doesn't currently invalidate element
2228  // iterators... that had better not change
2229  if (!keep_me)
2230  mesh.delete_elem(elem);
2231  }
2232 
2233  // Delete all the nodes we have no reason to save
2234  for (auto & node : mesh.node_ptr_range())
2235  {
2236  libmesh_assert(node);
2237  if (!connected_nodes.count(node))
2238  {
2239  libmesh_assert_not_equal_to(node->processor_id(),
2240  mesh.processor_id());
2241  mesh.delete_node(node);
2242  }
2243  }
2244 
2245  // If we had a point locator, it's invalid now that some of the
2246  // elements it pointed to have been deleted.
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.
2253  gf->delete_remote_elements();
2254 
2255 #ifdef DEBUG
2256  const dof_id_type n_new_constraint_rows = mesh.n_constraint_rows();
2257  libmesh_assert_equal_to(n_constraint_rows, n_new_constraint_rows);
2258 
2261 #endif
2262 }
2263 
2264 } // namespace libMesh
OStreamProxy err
void set_p_level(const unsigned int p)
Sets the value of the p-refinement level for the element.
void gather_packed_range(const unsigned int root_id, Context *context, Iter range_begin, const Iter range_end, OutputIter out, std::size_t approx_buffer_size=1000000) const
RefinementState refinement_flag() const
Definition: elem.h:3210
void libmesh_assert_topology_consistent_procids< Node >(const MeshBase &mesh)
Definition: mesh_tools.C:1483
const Elem * parent() const
Definition: elem.h:3030
void broadcast_packed_range(const Context *context1, Iter range_begin, const Iter range_end, OutputContext *context2, OutputIter out, const unsigned int root_id=0, std::size_t approx_buffer_size=1000000) const
void make_node_unique_ids_parallel_consistent(MeshBase &)
Assuming all unique_ids on local nodes are globally unique, and assuming all processor ids are parall...
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
Definition: mesh_base.h:1703
A Node is like a Point, but with more information.
Definition: node.h:52
virtual void renumber_node(dof_id_type old_id, dof_id_type new_id)=0
Changes the id of node old_id, both by changing node(old_id)->id() and by moving node(old_id) in the ...
std::vector< std::string > _elem_integer_names
The array of names for integer data associated with each element in the mesh.
Definition: mesh_base.h:2033
const Elem * interior_parent() const
Definition: elem.C:1186
void reconnect_nodes(connected_elem_set_type &connected_elements, connected_node_set_type &connected_nodes)
void correct_node_proc_ids(MeshBase &)
Changes the processor ids on each node so be the same as the id of the lowest element touching that n...
Definition: mesh_tools.C:2393
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2710
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:969
The definition of the const_element_iterator struct.
Definition: mesh_base.h:2216
const Elem * top_parent() const
Definition: elem.h:3056
std::vector< std::pair< std::pair< const Elem *, unsigned int >, Real > > constraint_rows_mapped_type
Definition: mesh_base.h:1697
MessageTag get_unique_tag(int tagvalue=MessageTag::invalid_tag) const
void make_elems_parallel_consistent(MeshBase &)
Copy ids of ghost elements from their local processors.
RefinementState p_refinement_flag() const
Definition: elem.h:3226
void libmesh_assert_valid_refinement_tree(const MeshBase &mesh)
A function for verifying that elements on this processor have valid descendants and consistent active...
Definition: mesh_tools.C:1547
std::unordered_set< const Node * > was_new
void clear()
Clears all data structures and resets to a pristine state.
void family_tree(std::vector< const Elem *> &family, bool reset=true) const
Fills the vector family with the children of this element, recursively.
Definition: elem.C:2101
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void libmesh_assert_valid_constraint_rows(const MeshBase &mesh)
A function for verifying that all mesh constraint rows express relations between nodes and elements t...
Definition: mesh_tools.C:1665
void active_family_tree(std::vector< const Elem *> &active_family, bool reset=true) const
Same as the family_tree() member, but only adds the active children.
Definition: elem.C:2133
void libmesh_assert_parallel_consistent_new_node_procids(const MeshBase &mesh)
A function for verifying that processor assignment is parallel consistent (every processor agrees on ...
Definition: mesh_tools.C:2065
const Parallel::Communicator & comm() const
unsigned int p_level() const
Definition: elem.h:3108
void allgather_packed_range(Context *context, Iter range_begin, const Iter range_end, OutputIter out, std::size_t approx_buffer_size=1000000) const
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
void make_node_bcids_parallel_consistent(MeshBase &)
Assuming all processor ids are parallel consistent, this function makes all ghost boundary ids on nod...
void make_links_to_me_remote()
Resets this element&#39;s neighbors&#39; appropriate neighbor pointers and its parent&#39;s and children&#39;s approp...
Definition: elem.C:1533
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
std::vector< std::string > _node_integer_names
The array of names for integer data associated with each node in the mesh.
Definition: mesh_base.h:2045
void send_packed_range(const unsigned int dest_processor_id, const Context *context, Iter range_begin, const Iter range_end, const MessageTag &tag=no_tag, std::size_t approx_buffer_size=1000000) const
dof_id_type n_constraint_rows() const
Definition: mesh_base.C:2044
Used to iterate over non-nullptr entries in a container.
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:2346
unsigned int paranoid_n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:826
ElemMappingType default_mapping_type() const
Returns the default master space to physical space mapping basis functions to be used on newly added ...
Definition: mesh_base.h:812
void libmesh_assert_equal_n_systems(const MeshBase &mesh)
The following functions, only available in builds with NDEBUG undefined, are for asserting internal c...
Definition: mesh_tools.C:1212
MPI_Status status
processor_id_type size() const
void push_parallel_vector_data(const Communicator &comm, MapToVectors &&data, const ActionFunctor &act_on_data)
uint8_t processor_id_type
std::map< boundary_id_type, std::string > & set_sideset_name_map()
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
void libmesh_ignore(const Args &...)
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
virtual void find_neighbors(const bool reset_remote_elements=false, const bool reset_current_list=true)=0
Locate element face (edge in 2D) neighbors.
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
void make_new_node_proc_ids_parallel_consistent(MeshBase &)
Assuming all processor ids on nodes touching local elements are parallel consistent, this function makes processor ids on new nodes on other processors parallel consistent as well.
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
dof_id_type id() const
Definition: dof_object.h:828
void query_ghosting_functors(const MeshBase &mesh, processor_id_type pid, MeshBase::const_element_iterator elem_it, MeshBase::const_element_iterator elem_end, connected_elem_set_type &connected_elements)
static const processor_id_type invalid_processor_id
An invalid processor_id to distinguish DoFs that have not been assigned to a processor.
Definition: dof_object.h:493
void push_parallel_packed_range(const Communicator &comm, MapToContainers &&data, Context *context, const ActionFunctor &act_on_data)
void receive_packed_range(const unsigned int dest_processor_id, Context *context, OutputIter out, const T *output_type, const MessageTag &tag=any_tag) const
unsigned char default_mapping_data() const
Returns any default data value used by the master space to physical space mapping.
Definition: mesh_base.h:830
virtual void update_parallel_id_counts()=0
Updates parallel caches so that methods like n_elem() accurately reflect changes on other processors...
void sync_node_data_by_element_id(MeshBase &mesh, const MeshBase::const_element_iterator &range_begin, const MeshBase::const_element_iterator &range_end, const ElemCheckFunctor &elem_check, const NodeCheckFunctor &node_check, SyncFunctor &sync)
Synchronize data about a range of ghost nodes uniquely identified by an element id and local node id...
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
status probe(const unsigned int src_processor_id, const MessageTag &tag=any_tag) const
void cache_elem_data()
Definition: mesh_base.C:1733
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:802
virtual const Node * query_node_ptr(const dof_id_type i) const =0
void make_p_levels_parallel_consistent(MeshBase &)
Copy p levels of ghost elements from their local processors.
void set_default_mapping_type(const ElemMappingType type)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
Definition: mesh_base.h:821
virtual dof_id_type max_elem_id() const =0
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1671
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
void family_tree(T elem, std::vector< T > &family, bool reset=true)
Definition: elem_internal.h:41
libmesh_assert(ctx)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
void libmesh_assert_valid_boundary_ids(const MeshBase &mesh)
A function for verifying that boundary condition ids match across processors.
Definition: mesh_tools.C:1708
virtual void delete_node(Node *n)=0
Removes the Node n from the mesh.
Helper for building element sides that minimizes the construction of new elements.
void allow_find_neighbors(bool allow)
If false is passed then this mesh will no longer work to find element neighbors when being prepared f...
Definition: mesh_base.h:1203
void hack_p_level_and_refinement_flag(const unsigned int p, RefinementState pflag)
Sets the value of the p-refinement level for the element without altering the p-level of its ancestor...
Definition: elem.h:3275
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
void set_neighbor(const unsigned int i, Elem *n)
Assigns n as the neighbor.
Definition: elem.h:2618
The DistributedMesh class is derived from the MeshBase class, and is intended to provide identical fu...
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:920
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
umap_type definitive_renumbering
void gather_neighboring_elements(DistributedMesh &) const
bool on_boundary() const
Definition: elem.h:2909
void libmesh_assert_parallel_consistent_procids< Node >(const MeshBase &mesh)
Definition: mesh_tools.C:2103
A do-nothing class for templated methods that expect output iterator arguments.
void connect_children(const MeshBase &mesh, MeshBase::const_element_iterator elem_it, MeshBase::const_element_iterator elem_end, connected_elem_set_type &connected_elements)
std::set< const Node * > connected_node_set_type
SimpleRange< NodeRefIter > node_ref_range()
Returns a range with all nodes of an element, usable in range-based for loops.
Definition: elem.h:2665
virtual const Elem * elem_ptr(const dof_id_type i) const =0
query_obj query
ElemMappingType
Enumeration of possible element master->physical mapping types.
void send_coarse_ghosts(MeshBase &) const
Examine a just-coarsened mesh, and for any newly-coarsened elements, send the associated ghosted elem...
void make_node_ids_parallel_consistent(MeshBase &)
Assuming all ids on local nodes are globally unique, and assuming all processor ids are parallel cons...
std::vector< dof_id_type > _node_integer_default_values
The array of default initialization values for integer data associated with each node in the mesh...
Definition: mesh_base.h:2051
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
std::set< GhostingFunctor * >::const_iterator ghosting_functors_begin() const
Beginning of range of ghosting functors.
Definition: mesh_base.h:1291
std::vector< dof_id_type > _elem_integer_default_values
The array of default initialization values for integer data associated with each element in the mesh...
Definition: mesh_base.h:2039
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void gather(const processor_id_type root_id, MeshBase &) const
This method takes an input DistributedMesh which may be distributed among all the processors...
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
timpi_pure bool verify(const T &r) const
virtual void renumber_elem(dof_id_type old_id, dof_id_type new_id)=0
Changes the id of element old_id, both by changing elem(old_id)->id() and by moving elem(old_id) in t...
void max(const T &r, T &o, Request &req) const
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
std::map< boundary_id_type, std::string > & set_nodeset_name_map()
bool sync_node_data_by_element_id_once(MeshBase &mesh, const MeshBase::const_element_iterator &range_begin, const MeshBase::const_element_iterator &range_end, const ElemCheckFunctor &elem_check, const NodeCheckFunctor &node_check, SyncFunctor &sync)
Synchronize data about a range of ghost nodes uniquely identified by an element id and local node id...
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
unsigned int n_neighbors() const
Definition: elem.h:714
void broadcast(MeshBase &) const
Finds all the processors that may contain elements that neighbor my elements.
void make_node_proc_ids_parallel_consistent(MeshBase &)
Assuming all processor ids on nodes touching local elements are parallel consistent, this function makes all other processor ids parallel consistent as well.
bool subactive() const
Definition: elem.h:2959
virtual const Elem & elem_ref(const dof_id_type i) const
Definition: mesh_base.h:639
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
void connect_element_dependencies(const MeshBase &mesh, connected_elem_set_type &connected_elements, connected_node_set_type &connected_nodes)
std::set< const Elem *, CompareElemIdsByLevel > connected_elem_set_type
void set_default_mapping_data(const unsigned char data)
Set the default master space to physical space mapping basis functions to be used on newly added elem...
Definition: mesh_base.h:839
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
uset_type definitive_ids
SimpleRange< NeighborPtrIter > neighbor_ptr_range()
Returns a range with all neighbors of an element, usable in range-based for loops.
Definition: elem.h:3503
renumber_obj renumber
virtual dof_id_type max_node_id() const =0
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
void make_new_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on new nodes from their local processors.
const DofMap &dof_map LIBMESH_COMMA unsigned int std::string & set_subdomain_name_map()
Definition: mesh_base.h:1692
processor_id_type processor_id() const
Definition: dof_object.h:905
void delete_remote_elements(DistributedMesh &, const std::set< Elem *> &) const
This method takes an input DistributedMesh which may be distributed among all the processors...
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
uint8_t unique_id_type
Definition: id_types.h:86
void connect_families(connected_elem_set_type &connected_elements, const MeshBase *mesh=nullptr)
bool has_children() const
Definition: elem.h:2979
void redistribute(DistributedMesh &mesh, bool newly_coarsened_only=false) const
This method takes a parallel distributed mesh and redistributes the elements.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
virtual dof_id_type n_nodes() const =0
void sync_element_data_by_parent_id(MeshBase &mesh, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost elements uniquely identified by their parent id and which child t...
std::set< GhostingFunctor * >::const_iterator ghosting_functors_end() const
End of range of ghosting functors.
Definition: mesh_base.h:1297
uint8_t dof_id_type
Definition: id_types.h:67
void set_union(T &data, const unsigned int root_id) const
const RemoteElem * remote_elem
Definition: remote_elem.C:57