libMesh
mesh_tools.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/elem.h"
22 #include "libmesh/elem_range.h"
23 #include "libmesh/libmesh_logging.h"
24 #include "libmesh/mesh_base.h"
25 #include "libmesh/mesh_communication.h"
26 #include "libmesh/mesh_serializer.h"
27 #include "libmesh/mesh_tools.h"
28 #include "libmesh/node_range.h"
29 #include "libmesh/parallel.h"
30 #include "libmesh/parallel_algebra.h"
31 #include "libmesh/parallel_ghost_sync.h"
32 #include "libmesh/sphere.h"
33 #include "libmesh/threads.h"
34 #include "libmesh/enum_to_string.h"
35 #include "libmesh/enum_elem_type.h"
36 #include "libmesh/int_range.h"
37 #include "libmesh/utility.h"
38 #include "libmesh/boundary_info.h"
39 
40 #ifndef NDEBUG
41 # include "libmesh/remote_elem.h"
42 #endif
43 
44 // C++ includes
45 #include <limits>
46 #include <numeric> // for std::accumulate
47 #include <set>
48 #include <unordered_map>
49 #include <unordered_set>
50 
51 
52 
53 // ------------------------------------------------------------
54 // anonymous namespace for helper classes
55 namespace {
56 
57 using namespace libMesh;
58 
66 class SumElemWeight
67 {
68 public:
69  SumElemWeight () :
70  _weight(0)
71  {}
72 
73  SumElemWeight (SumElemWeight &, Threads::split) :
74  _weight(0)
75  {}
76 
77  void operator()(const ConstElemRange & range)
78  {
79  for (const auto & elem : range)
80  _weight += elem->n_nodes();
81  }
82 
83  dof_id_type weight() const
84  { return _weight; }
85 
86  // If we don't have threads we never need a join, and icpc yells a
87  // warning if it sees an anonymous function that's never used
88 #if LIBMESH_USING_THREADS
89  void join (const SumElemWeight & other)
90  { _weight += other.weight(); }
91 #endif
92 
93 private:
94  dof_id_type _weight;
95 };
96 
97 
104 class FindBBox
105 {
106 public:
107  FindBBox () : _bbox()
108  {}
109 
110  FindBBox (FindBBox & other, Threads::split) :
111  _bbox(other._bbox)
112  {}
113 
114  void operator()(const ConstNodeRange & range)
115  {
116  for (const auto & node : range)
117  {
118  libmesh_assert(node);
119  _bbox.union_with(*node);
120  }
121  }
122 
123  void operator()(const ConstElemRange & range)
124  {
125  for (const auto & elem : range)
126  {
127  libmesh_assert(elem);
128  _bbox.union_with(elem->loose_bounding_box());
129  }
130  }
131 
132  Point & min() { return _bbox.min(); }
133 
134  Point & max() { return _bbox.max(); }
135 
136  // If we don't have threads we never need a join, and icpc yells a
137  // warning if it sees an anonymous function that's never used
138 #if LIBMESH_USING_THREADS
139  void join (const FindBBox & other)
140  {
141  _bbox.union_with(other._bbox);
142  }
143 #endif
144 
145  libMesh::BoundingBox & bbox ()
146  {
147  return _bbox;
148  }
149 
150 private:
151  BoundingBox _bbox;
152 };
153 
154 #ifdef DEBUG
155 void assert_semiverify_dofobj(const Parallel::Communicator & communicator,
156  const DofObject * d,
157  unsigned int sysnum = libMesh::invalid_uint)
158 {
159  if (d)
160  {
161  const unsigned int n_sys = d->n_systems();
162 
163  std::vector<unsigned int> n_vars (n_sys, 0);
164  for (unsigned int s = 0; s != n_sys; ++s)
165  if (sysnum == s ||
166  sysnum == libMesh::invalid_uint)
167  n_vars[s] = d->n_vars(s);
168 
169  const unsigned int tot_n_vars =
170  std::accumulate(n_vars.begin(), n_vars.end(), 0);
171 
172  std::vector<unsigned int> n_comp (tot_n_vars, 0);
173  std::vector<dof_id_type> first_dof (tot_n_vars, 0);
174 
175  for (unsigned int s = 0, i=0; s != n_sys; ++s)
176  {
177  if (sysnum != s &&
178  sysnum != libMesh::invalid_uint)
179  continue;
180 
181  for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
182  {
183  n_comp[i] = d->n_comp(s,v);
184  first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : DofObject::invalid_id;
185  }
186  }
187 
188  libmesh_assert(communicator.semiverify(&n_sys));
189  libmesh_assert(communicator.semiverify(&n_vars));
190  libmesh_assert(communicator.semiverify(&n_comp));
191  libmesh_assert(communicator.semiverify(&first_dof));
192  }
193  else
194  {
195  const unsigned int * p_ui = nullptr;
196  const std::vector<unsigned int> * p_vui = nullptr;
197  const std::vector<dof_id_type> * p_vdid = nullptr;
198 
199  libmesh_assert(communicator.semiverify(p_ui));
200  libmesh_assert(communicator.semiverify(p_vui));
201  libmesh_assert(communicator.semiverify(p_vui));
202  libmesh_assert(communicator.semiverify(p_vdid));
203  }
204 }
205 
206 
207 
208 #ifdef LIBMESH_ENABLE_UNIQUE_ID
209 void assert_dofobj_unique_id(const Parallel::Communicator & comm,
210  const DofObject * d,
211  const std::unordered_set<unique_id_type> & unique_ids)
212 {
213  // Duplicating some semiverify code here so we can reuse
214  // tempmin,tempmax afterward
215 
216  unique_id_type tempmin, tempmax;
217  if (d)
218  {
219  tempmin = tempmax = d->unique_id();
220  }
221  else
222  {
225  }
226  comm.min(tempmin);
227  comm.max(tempmax);
228  bool invalid = d && ((d->unique_id() != tempmin) ||
229  (d->unique_id() != tempmax));
230  comm.max(invalid);
231 
232  // First verify that everything is in sync
233  libmesh_assert(!invalid);
234 
235  // Then verify that any remote id doesn't duplicate a local one.
236  if (!d && tempmin == tempmax)
237  libmesh_assert(!unique_ids.count(tempmin));
238 }
239 #endif // LIBMESH_ENABLE_UNIQUE_ID
240 #endif // DEBUG
241 
242 void find_nodal_neighbors_helper(const dof_id_type global_id,
243  const std::vector<const Elem *> & node_to_elem_vec,
244  std::vector<const Node *> & neighbors)
245 {
246  // We'll construct a std::set<const Node *> for more efficient
247  // searching while finding the nodal neighbors, and return it to the
248  // user in a std::vector.
249  std::set<const Node *> neighbor_set;
250 
251  // Look through the elements that contain this node
252  // find the local node id... then find the side that
253  // node lives on in the element
254  // next, look for the _other_ node on that side
255  // That other node is a "nodal_neighbor"... save it
256  for (const auto & elem : node_to_elem_vec)
257  {
258  // We only care about active elements...
259  if (elem->active())
260  {
261  // Which local node number is global_id?
262  unsigned local_node_number = elem->local_node(global_id);
263 
264  // Make sure it was found
265  libmesh_assert_not_equal_to(local_node_number, libMesh::invalid_uint);
266 
267  const unsigned short n_edges = elem->n_edges();
268 
269  // If this element has no edges, the edge-based algorithm below doesn't make sense.
270  if (!n_edges)
271  {
272  switch (elem->type())
273  {
274  case EDGE2:
275  {
276  switch (local_node_number)
277  {
278  case 0:
279  // The other node is a nodal neighbor
280  neighbor_set.insert(elem->node_ptr(1));
281  break;
282 
283  case 1:
284  // The other node is a nodal neighbor
285  neighbor_set.insert(elem->node_ptr(0));
286  break;
287 
288  default:
289  libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
290  }
291  break;
292  }
293 
294  case EDGE3:
295  {
296  switch (local_node_number)
297  {
298  // The outside nodes have node 2 as a neighbor
299  case 0:
300  case 1:
301  neighbor_set.insert(elem->node_ptr(2));
302  break;
303 
304  // The middle node has the outer nodes as neighbors
305  case 2:
306  neighbor_set.insert(elem->node_ptr(0));
307  neighbor_set.insert(elem->node_ptr(1));
308  break;
309 
310  default:
311  libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
312  }
313  break;
314  }
315 
316  case EDGE4:
317  {
318  switch (local_node_number)
319  {
320  case 0:
321  // The left-middle node is a nodal neighbor
322  neighbor_set.insert(elem->node_ptr(2));
323  break;
324 
325  case 1:
326  // The right-middle node is a nodal neighbor
327  neighbor_set.insert(elem->node_ptr(3));
328  break;
329 
330  // The left-middle node
331  case 2:
332  neighbor_set.insert(elem->node_ptr(0));
333  neighbor_set.insert(elem->node_ptr(3));
334  break;
335 
336  // The right-middle node
337  case 3:
338  neighbor_set.insert(elem->node_ptr(1));
339  neighbor_set.insert(elem->node_ptr(2));
340  break;
341 
342  default:
343  libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
344  }
345  break;
346  }
347 
348  default:
349  libmesh_error_msg("Unrecognized ElemType: " << Utility::enum_to_string(elem->type()) << std::endl);
350  }
351  }
352 
353  // Index of the current edge
354  unsigned current_edge = 0;
355 
356  const unsigned short n_nodes = elem->n_nodes();
357 
358  while (current_edge < n_edges)
359  {
360  // Find the edge the node is on
361  bool found_edge = false;
362  for (; current_edge<n_edges; ++current_edge)
363  if (elem->is_node_on_edge(local_node_number, current_edge))
364  {
365  found_edge = true;
366  break;
367  }
368 
369  // Did we find one?
370  if (found_edge)
371  {
372  const Node * node_to_save = nullptr;
373 
374  // Find another node in this element on this edge
375  for (unsigned other_node_this_edge = 0; other_node_this_edge != n_nodes; other_node_this_edge++)
376  if ( (elem->is_node_on_edge(other_node_this_edge, current_edge)) && // On the current edge
377  (elem->node_id(other_node_this_edge) != global_id)) // But not the original node
378  {
379  // We've found a nodal neighbor! Save a pointer to it..
380  node_to_save = elem->node_ptr(other_node_this_edge);
381  break;
382  }
383 
384  // Make sure we found something
385  libmesh_assert(node_to_save != nullptr);
386 
387  neighbor_set.insert(node_to_save);
388  }
389 
390  // Keep looking for edges, node may be on more than one edge
391  current_edge++;
392  }
393  } // if (elem->active())
394  } // for
395 
396  // Assign the entries from the set to the vector. Note: this
397  // replaces any existing contents in neighbors and modifies its size
398  // accordingly.
399  neighbors.assign(neighbor_set.begin(), neighbor_set.end());
400 }
401 
402 }
403 
404 
405 namespace libMesh
406 {
407 
408 // ------------------------------------------------------------
409 // MeshTools functions
410 
411 namespace MeshTools
412 {
413 
415 {
416  if (!mesh.is_serial())
417  {
418  libmesh_parallel_only(mesh.comm());
420  mesh.comm().sum(weight);
421  dof_id_type unpartitioned_weight =
423  return weight + unpartitioned_weight;
424  }
425 
426  SumElemWeight sew;
427 
428  Threads::parallel_reduce (ConstElemRange (mesh.elements_begin(),
429  mesh.elements_end()),
430  sew);
431  return sew.weight();
432 
433 }
434 
435 
436 
438 {
439  SumElemWeight sew;
440 
441  Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
442  mesh.pid_elements_end(pid)),
443  sew);
444  return sew.weight();
445 }
446 
447 
448 
450  std::vector<std::vector<dof_id_type>> & nodes_to_elem_map)
451 {
452  // A vector indexed over all nodes is too inefficient to use for a
453  // distributed mesh. Use the unordered_map API instead.
454  if (!mesh.is_serial())
455  libmesh_deprecated();
456 
457  nodes_to_elem_map.resize (mesh.max_node_id());
458 
459  for (const auto & elem : mesh.element_ptr_range())
460  for (auto & node : elem->node_ref_range())
461  {
462  libmesh_assert_less (node.id(), nodes_to_elem_map.size());
463  libmesh_assert_less (elem->id(), mesh.n_elem());
464 
465  nodes_to_elem_map[node.id()].push_back(elem->id());
466  }
467 }
468 
469 
470 
472  std::vector<std::vector<const Elem *>> & nodes_to_elem_map)
473 {
474  // A vector indexed over all nodes is too inefficient to use for a
475  // distributed mesh. Use the unordered_map API instead.
476  if (!mesh.is_serial())
477  libmesh_deprecated();
478 
479  nodes_to_elem_map.resize (mesh.max_node_id());
480 
481  for (const auto & elem : mesh.element_ptr_range())
482  for (auto & node : elem->node_ref_range())
483  {
484  libmesh_assert_less (node.id(), nodes_to_elem_map.size());
485 
486  nodes_to_elem_map[node.id()].push_back(elem);
487  }
488 }
489 
490 
491 
493  std::unordered_map<dof_id_type, std::vector<dof_id_type>> & nodes_to_elem_map)
494 {
495  nodes_to_elem_map.clear();
496 
497  for (const auto & elem : mesh.element_ptr_range())
498  for (auto & node : elem->node_ref_range())
499  nodes_to_elem_map[node.id()].push_back(elem->id());
500 }
501 
502 
503 
505  std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
506 {
507  nodes_to_elem_map.clear();
508 
509  for (const auto & elem : mesh.element_ptr_range())
510  for (auto & node : elem->node_ref_range())
511  nodes_to_elem_map[node.id()].push_back(elem);
512 }
513 
514 
515 
516 std::unordered_set<dof_id_type>
518 {
519  std::unordered_set<dof_id_type> boundary_nodes;
520 
521  // Loop over elements, find those on boundary, and
522  // mark them as true in on_boundary.
523  for (const auto & elem : mesh.active_element_ptr_range())
524  for (auto s : elem->side_index_range())
525  if (elem->neighbor_ptr(s) == nullptr) // on the boundary
526  {
527  auto nodes_on_side = elem->nodes_on_side(s);
528 
529  for (auto & local_id : nodes_on_side)
530  boundary_nodes.insert(elem->node_ptr(local_id)->id());
531  }
532 
533  return boundary_nodes;
534 }
535 
536 std::unordered_set<dof_id_type>
538 {
539  std::unordered_set<dof_id_type> block_boundary_nodes;
540 
541  // Loop over elements, find those on boundary, and
542  // mark them as true in on_boundary.
543  for (const auto & elem : mesh.active_element_ptr_range())
544  for (auto s : elem->side_index_range())
545  if (elem->neighbor_ptr(s) && (elem->neighbor_ptr(s)->subdomain_id() != elem->subdomain_id()))
546  {
547  auto nodes_on_side = elem->nodes_on_side(s);
548 
549  for (auto & local_id : nodes_on_side)
550  block_boundary_nodes.insert(elem->node_ptr(local_id)->id());
551  }
552 
553  return block_boundary_nodes;
554 }
555 
556 
557 
560 {
561  // This function must be run on all processors at once
562  libmesh_parallel_only(mesh.comm());
563 
564  FindBBox find_bbox;
565 
566  // Start with any unpartitioned elements we know about locally
568  mesh.pid_elements_end(DofObject::invalid_processor_id)),
569  find_bbox);
570 
571  // And combine with our local elements
572  find_bbox.bbox().union_with(create_local_bounding_box(mesh));
573 
574  // Compare the bounding boxes across processors
575  mesh.comm().min(find_bbox.min());
576  mesh.comm().max(find_bbox.max());
577 
578  return find_bbox.bbox();
579 }
580 
581 
582 
585 {
586  // This function must be run on all processors at once
587  libmesh_parallel_only(mesh.comm());
588 
589  FindBBox find_bbox;
590 
591  // Start with any unpartitioned nodes we know about locally
593  mesh.pid_nodes_end(DofObject::invalid_processor_id)),
594  find_bbox);
595 
596  // Add our local nodes
597  Threads::parallel_reduce (ConstNodeRange (mesh.local_nodes_begin(),
598  mesh.local_nodes_end()),
599  find_bbox);
600 
601  // Compare the bounding boxes across processors
602  mesh.comm().min(find_bbox.min());
603  mesh.comm().max(find_bbox.max());
604 
605  return find_bbox.bbox();
606 }
607 
608 
609 
610 Sphere
612 {
614 
615  const Real diag = (bbox.second - bbox.first).norm();
616  const Point cent = (bbox.second + bbox.first)/2;
617 
618  return Sphere (cent, .5*diag);
619 }
620 
621 
622 
625 {
626  FindBBox find_bbox;
627 
628  Threads::parallel_reduce (ConstElemRange (mesh.local_elements_begin(),
629  mesh.local_elements_end()),
630  find_bbox);
631 
632  return find_bbox.bbox();
633 }
634 
635 
636 
639  const processor_id_type pid)
640 {
641  // This can only be run in parallel, with consistent arguments.
642  libmesh_parallel_only(mesh.comm());
643  libmesh_assert(mesh.comm().verify(pid));
644 
645  libmesh_assert_less (pid, mesh.n_processors());
646 
647  FindBBox find_bbox;
648 
649  Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
650  mesh.pid_elements_end(pid)),
651  find_bbox);
652 
653  // Compare the bounding boxes across processors
654  mesh.comm().min(find_bbox.min());
655  mesh.comm().max(find_bbox.max());
656 
657  return find_bbox.bbox();
658 }
659 
660 
661 
662 Sphere
664  const processor_id_type pid)
665 {
666  libMesh::BoundingBox bbox =
668 
669  const Real diag = (bbox.second - bbox.first).norm();
670  const Point cent = (bbox.second + bbox.first)/2;
671 
672  return Sphere (cent, .5*diag);
673 }
674 
675 
676 
679  const subdomain_id_type sid)
680 {
681  // This can only be run in parallel, with consistent arguments.
682  libmesh_parallel_only(mesh.comm());
683  libmesh_assert(mesh.comm().verify(sid));
684 
685  FindBBox find_bbox;
686 
688  (ConstElemRange (mesh.active_local_subdomain_elements_begin(sid),
689  mesh.active_local_subdomain_elements_end(sid)),
690  find_bbox);
691 
692  // Compare the bounding boxes across processors
693  mesh.comm().min(find_bbox.min());
694  mesh.comm().max(find_bbox.max());
695 
696  return find_bbox.bbox();
697 }
698 
699 
700 
701 Sphere
703  const subdomain_id_type sid)
704 {
705  libMesh::BoundingBox bbox =
707 
708  const Real diag = (bbox.second - bbox.first).norm();
709  const Point cent = (bbox.second + bbox.first)/2;
710 
711  return Sphere (cent, .5*diag);
712 }
713 
714 
715 
716 void elem_types (const MeshBase & mesh,
717  std::vector<ElemType> & et)
718 {
719  // Loop over the the elements. If the current element type isn't in
720  // the vector, insert it.
721  for (const auto & elem : mesh.element_ptr_range())
722  if (!std::count(et.begin(), et.end(), elem->type()))
723  et.push_back(elem->type());
724 }
725 
726 
727 
729  const ElemType type)
730 {
731  return static_cast<dof_id_type>(std::distance(mesh.type_elements_begin(type),
732  mesh.type_elements_end (type)));
733 }
734 
735 
736 
738  const ElemType type)
739 {
740  return static_cast<dof_id_type>(std::distance(mesh.active_type_elements_begin(type),
741  mesh.active_type_elements_end (type)));
742 }
743 
745  const ElemType type,
746  const unsigned int level)
747 {
748  dof_id_type cnt = 0;
749 
750  // iterate over the elements of the specified type
751  for (const auto & elem : as_range(mesh.type_elements_begin(type),
752  mesh.type_elements_end(type)))
753  if ((elem->level() == level) && !elem->subactive())
754  cnt++;
755 
756  return cnt;
757 }
758 
759 
760 unsigned int n_active_local_levels(const MeshBase & mesh)
761 {
762  unsigned int nl = 0;
763 
764  for (auto & elem : mesh.active_local_element_ptr_range())
765  nl = std::max(elem->level() + 1, nl);
766 
767  return nl;
768 }
769 
770 
771 
772 unsigned int n_active_levels(const MeshBase & mesh)
773 {
774  libmesh_parallel_only(mesh.comm());
775 
776  unsigned int nl = n_active_local_levels(mesh);
777 
778  for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
779  mesh.unpartitioned_elements_end()))
780  if (elem->active())
781  nl = std::max(elem->level() + 1, nl);
782 
783  mesh.comm().max(nl);
784  return nl;
785 }
786 
787 
788 
789 unsigned int n_local_levels(const MeshBase & mesh)
790 {
791  unsigned int nl = 0;
792 
793  for (const auto & elem : as_range(mesh.local_elements_begin(),
794  mesh.local_elements_end()))
795  nl = std::max(elem->level() + 1, nl);
796 
797  return nl;
798 }
799 
800 
801 
802 unsigned int n_levels(const MeshBase & mesh)
803 {
804  libmesh_parallel_only(mesh.comm());
805 
806  unsigned int nl = n_local_levels(mesh);
807 
808  for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
809  mesh.unpartitioned_elements_end()))
810  nl = std::max(elem->level() + 1, nl);
811 
812  mesh.comm().max(nl);
813 
814  // n_levels() is only valid and should only be called in cases where
815  // the mesh is validly distributed (or serialized). Let's run an
816  // expensive test in debug mode to make sure this is such a case.
817 #ifdef DEBUG
818  const unsigned int paranoid_nl = paranoid_n_levels(mesh);
819  libmesh_assert_equal_to(nl, paranoid_nl);
820 #endif
821  return nl;
822 }
823 
824 
825 
826 unsigned int paranoid_n_levels(const MeshBase & mesh)
827 {
828  libmesh_parallel_only(mesh.comm());
829 
830  unsigned int nl = 0;
831  for (const auto & elem : mesh.element_ptr_range())
832  nl = std::max(elem->level() + 1, nl);
833 
834  mesh.comm().max(nl);
835  return nl;
836 }
837 
838 
839 
841  Real constraint_tol)
842 {
843  LOG_SCOPE("n_connected_components()", "MeshTools");
844 
845  // Yes, I'm being lazy. This is for mesh analysis before a
846  // simulation, not anything going in any loops.
847  if (!mesh.is_serial_on_zero())
848  libmesh_not_implemented();
849 
850  dof_id_type n_components = 0;
851 
852  if (mesh.processor_id())
853  {
854  mesh.comm().broadcast(n_components);
855  return n_components;
856  }
857 
858  // All nodes in a set here are connected (at least indirectly) to
859  // all other nodes in the same set, but have not yet been discovered
860  // to be connected to nodes in other sets.
861  std::vector<std::unordered_set<const Node *>> components;
862 
863  // With a typical mesh with few components and somewhat-contiguous
864  // ordering, vector performance should be fine. With a mesh with
865  // many components or completely scrambled ordering, performance
866  // can be a disaster.
867  auto find_component = [&components](const Node * n) {
868  std::unordered_set<const Node *> * component = nullptr;
869 
870  for (auto & c: components)
871  if (c.find(n) != c.end())
872  {
873  libmesh_assert(component == nullptr);
874  component = &c;
875  }
876 
877  return component;
878  };
879 
880  auto add_to_component =
881  [&find_component]
882  (std::unordered_set<const Node *> & component, const Node * n) {
883 
884  auto current_component = find_component(n);
885  // We may already know we're in the desired component
886  if (&component == current_component)
887  return;
888 
889  // If we're unknown, we should be in the desired component
890  else if (!current_component)
891  component.insert(n);
892 
893  // If we think we're in another component, it should actually be
894  // part of the desired component
895  else
896  {
897  component.merge(*current_component);
898  libmesh_assert(current_component->empty());
899  }
900  };
901 
902  auto & constraint_rows = mesh.get_constraint_rows();
903 
904  for (const auto & elem : mesh.element_ptr_range())
905  {
906  const Node * first_node = elem->node_ptr(0);
907 
908  auto component = find_component(first_node);
909 
910  // If we didn't find one, make a new one, reusing an existing
911  // slot if possible or growing our vector if necessary
912  if (!component)
913  for (auto & c: components)
914  if (c.empty())
915  component = &c;
916 
917  if (!component)
918  component = &components.emplace_back();
919 
920  for (const Node & n : elem->node_ref_range())
921  {
922  add_to_component(*component, &n);
923 
924  auto it = constraint_rows.find(&n);
925  if (it == constraint_rows.end())
926  continue;
927 
928  for (const auto & [pr, val] : it->second)
929  {
930  // Ignore too-trivial constraint coefficients if
931  // we get a non-default-0 constraint_tol
932  if (std::abs(val) < constraint_tol)
933  continue;
934 
935  const Elem * spline_elem = pr.first;
936  libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
937 
938  const Node * spline_node =
939  spline_elem->node_ptr(pr.second);
940 
941  add_to_component(*component, spline_node);
942  }
943  }
944  }
945 
946  for (auto & component : components)
947  if (!component.empty())
948  ++n_components;
949 
950  // We calculated this on proc 0; now let everyone else know too
951  mesh.comm().broadcast(n_components);
952 
953  return n_components;
954 }
955 
956 
957 
959  std::set<dof_id_type> & not_subactive_node_ids)
960 {
961  for (const auto & elem : mesh.element_ptr_range())
962  if (!elem->subactive())
963  for (auto & n : elem->node_ref_range())
964  not_subactive_node_ids.insert(n.id());
965 }
966 
967 
968 
971 {
972  return cast_int<dof_id_type>(std::distance(begin, end));
973 }
974 
975 
976 
978  const MeshBase::const_node_iterator & end)
979 {
980  return cast_int<dof_id_type>(std::distance(begin, end));
981 }
982 
983 
984 
986  unsigned int dim)
987 {
988  libmesh_parallel_only(mesh.comm());
989 
990  if (dim == libMesh::invalid_uint)
991  dim = mesh.mesh_dimension();
992 
993  Real vol = 0;
994 
995  // first my local elements
996  for (const auto & elem : as_range(mesh.local_elements_begin(),
997  mesh.local_elements_end()))
998  if (elem->dim() == dim)
999  vol += elem->volume();
1000 
1001  // then count any unpartitioned objects, once
1002  if (mesh.processor_id() == 0)
1003  for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1004  mesh.unpartitioned_elements_end()))
1005  if (elem->dim() == dim)
1006  vol += elem->volume();
1007 
1008  mesh.comm().sum(vol);
1009  return vol;
1010 }
1011 
1012 
1013 
1014 unsigned int n_p_levels (const MeshBase & mesh)
1015 {
1016  libmesh_parallel_only(mesh.comm());
1017 
1018  unsigned int max_p_level = 0;
1019 
1020  // first my local elements
1021  for (const auto & elem : as_range(mesh.local_elements_begin(),
1022  mesh.local_elements_end()))
1023  max_p_level = std::max(elem->p_level(), max_p_level);
1024 
1025  // then any unpartitioned objects
1026  for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1027  mesh.unpartitioned_elements_end()))
1028  max_p_level = std::max(elem->p_level(), max_p_level);
1029 
1030  mesh.comm().max(max_p_level);
1031  return max_p_level + 1;
1032 }
1033 
1034 
1035 
1037  const Node & node,
1038  const std::vector<std::vector<const Elem *>> & nodes_to_elem_map,
1039  std::vector<const Node *> & neighbors)
1040 {
1041  find_nodal_neighbors_helper(node.id(), nodes_to_elem_map[node.id()],
1042  neighbors);
1043 }
1044 
1045 
1046 
1048  const Node & node,
1049  const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1050  std::vector<const Node *> & neighbors)
1051 {
1052  const std::vector<const Elem *> node_to_elem_vec =
1053  libmesh_map_find(nodes_to_elem_map, node.id());
1054  find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
1055 }
1056 
1057 
1058 
1060  std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes)
1061 {
1062  // Loop through all the elements
1063  for (auto & elem : mesh.active_local_element_ptr_range())
1064  if (elem->type() == QUAD4)
1065  for (auto s : elem->side_index_range())
1066  {
1067  // Loop over the sides looking for sides that have hanging nodes
1068  // This code is inspired by compute_proj_constraints()
1069  const Elem * neigh = elem->neighbor_ptr(s);
1070 
1071  // If not a boundary side
1072  if (neigh != nullptr)
1073  {
1074  // Is there a coarser element next to this one?
1075  if (neigh->level() < elem->level())
1076  {
1077  const Elem * ancestor = elem;
1078  while (neigh->level() < ancestor->level())
1079  ancestor = ancestor->parent();
1080  unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
1081  libmesh_assert_less (s_neigh, neigh->n_neighbors());
1082 
1083  // Couple of helper uints...
1084  unsigned int local_node1=0;
1085  unsigned int local_node2=0;
1086 
1087  bool found_in_neighbor = false;
1088 
1089  // Find the two vertices that make up this side
1090  while (!elem->is_node_on_side(local_node1++,s)) { }
1091  local_node1--;
1092 
1093  // Start looking for the second one with the next node
1094  local_node2=local_node1+1;
1095 
1096  // Find the other one
1097  while (!elem->is_node_on_side(local_node2++,s)) { }
1098  local_node2--;
1099 
1100  //Pull out their global ids:
1101  dof_id_type node1 = elem->node_id(local_node1);
1102  dof_id_type node2 = elem->node_id(local_node2);
1103 
1104  // Now find which node is present in the neighbor
1105  // FIXME This assumes a level one rule!
1106  // The _other_ one is the hanging node
1107 
1108  // First look for the first one
1109  // FIXME could be streamlined a bit
1110  for (unsigned int n=0;n<neigh->n_sides();n++)
1111  if (neigh->node_id(n) == node1)
1112  found_in_neighbor=true;
1113 
1114  dof_id_type hanging_node=0;
1115 
1116  if (!found_in_neighbor)
1117  hanging_node=node1;
1118  else // If it wasn't node1 then it must be node2!
1119  hanging_node=node2;
1120 
1121  // Reset for reuse
1122  local_node1=0;
1123 
1124  // Find the first node that makes up the side in the neighbor (these should be the parent nodes)
1125  while (!neigh->is_node_on_side(local_node1++,s_neigh)) { }
1126  local_node1--;
1127 
1128  local_node2=local_node1+1;
1129 
1130  // Find the second node...
1131  while (!neigh->is_node_on_side(local_node2++,s_neigh)) { }
1132  local_node2--;
1133 
1134  // Save them if we haven't already found the parents for this one
1135  if (hanging_nodes[hanging_node].size()<2)
1136  {
1137  hanging_nodes[hanging_node].push_back(neigh->node_id(local_node1));
1138  hanging_nodes[hanging_node].push_back(neigh->node_id(local_node2));
1139  }
1140  }
1141  }
1142  }
1143 }
1144 
1145 
1146 
1148 {
1149  std::vector<Elem *> nodeelem_to_delete;
1150 
1151  for (auto & elem : mesh.element_ptr_range())
1152  if (elem->type() == NODEELEM &&
1153  elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
1154  nodeelem_to_delete.push_back(elem);
1155 
1156  auto & constraint_rows = mesh.get_constraint_rows();
1157 
1158  // All our constraint_rows ought to be for spline constraints we're
1159  // about to get rid of.
1160 #ifndef NDEBUG
1161  for (auto & node_row : constraint_rows)
1162  for (auto pr : node_row.second)
1163  {
1164  const Elem * elem = pr.first.first;
1165  libmesh_assert(elem->type() == NODEELEM);
1167  }
1168 #endif
1169 
1170  constraint_rows.clear();
1171 
1172  for (Elem * elem : nodeelem_to_delete)
1173  {
1174  Node * node = elem->node_ptr(0);
1175  mesh.delete_elem(elem);
1176  mesh.delete_node(node);
1177  }
1178 }
1179 
1180 
1181 
1183 {
1184  LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools");
1185 
1186  if (!mesh.is_prepared())
1187  return true;
1188 
1189  std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
1190 
1191  // Try preparing (without allowing repartitioning or renumbering, to
1192  // avoid false assertion failures)
1193  bool old_allow_renumbering = mesh_clone->allow_renumbering();
1194  mesh_clone->allow_renumbering(false);
1195  bool old_allow_remote_element_removal =
1196  mesh_clone->allow_remote_element_removal();
1197  bool old_skip_partitioning = mesh_clone->skip_partitioning();
1198  mesh_clone->skip_partitioning(true);
1199  mesh_clone->allow_remote_element_removal(false);
1200  mesh_clone->prepare_for_use();
1201  mesh_clone->allow_renumbering(old_allow_renumbering);
1202  mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal);
1203  mesh_clone->skip_partitioning(old_skip_partitioning);
1204 
1205  return (mesh == *mesh_clone);
1206 }
1207 
1208 
1209 
1210 #ifndef NDEBUG
1211 
1213 {
1214  LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
1215 
1216  unsigned int n_sys = libMesh::invalid_uint;
1217 
1218  for (const auto & elem : mesh.element_ptr_range())
1219  {
1220  if (n_sys == libMesh::invalid_uint)
1221  n_sys = elem->n_systems();
1222  else
1223  libmesh_assert_equal_to (elem->n_systems(), n_sys);
1224  }
1225 
1226  for (const auto & node : mesh.node_ptr_range())
1227  {
1228  if (n_sys == libMesh::invalid_uint)
1229  n_sys = node->n_systems();
1230  else
1231  libmesh_assert_equal_to (node->n_systems(), n_sys);
1232  }
1233 }
1234 
1235 
1236 
1237 #ifdef LIBMESH_ENABLE_AMR
1239 {
1240  LOG_SCOPE("libmesh_assert_old_dof_objects()", "MeshTools");
1241 
1242  for (const auto & elem : mesh.element_ptr_range())
1243  {
1244  if (elem->refinement_flag() == Elem::JUST_REFINED ||
1245  elem->refinement_flag() == Elem::INACTIVE)
1246  continue;
1247 
1248  if (elem->has_dofs())
1249  libmesh_assert(elem->get_old_dof_object());
1250 
1251  for (auto & node : elem->node_ref_range())
1252  if (node.has_dofs())
1253  libmesh_assert(node.get_old_dof_object());
1254  }
1255 }
1256 #else
1257 void libmesh_assert_old_dof_objects (const MeshBase &) {}
1258 #endif // LIBMESH_ENABLE_AMR
1259 
1260 
1261 
1263 {
1264  LOG_SCOPE("libmesh_assert_valid_node_pointers()", "MeshTools");
1265 
1266  // Here we specifically do not want "auto &" because we need to
1267  // reseat the (temporary) pointer variable in the loop below,
1268  // without modifying the original.
1269  for (const Elem * elem : mesh.element_ptr_range())
1270  {
1271  libmesh_assert (elem);
1272  while (elem)
1273  {
1274  elem->libmesh_assert_valid_node_pointers();
1275  for (auto n : elem->neighbor_ptr_range())
1276  if (n && n != remote_elem)
1277  n->libmesh_assert_valid_node_pointers();
1278 
1279  libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1280  elem = elem->parent();
1281  }
1282  }
1283 }
1284 
1285 
1286 
1288 {
1289  LOG_SCOPE("libmesh_assert_valid_remote_elems()", "MeshTools");
1290 
1291  for (const auto & elem : as_range(mesh.local_elements_begin(),
1292  mesh.local_elements_end()))
1293  {
1294  libmesh_assert (elem);
1295 
1296  // We currently don't allow active_local_elements to have
1297  // remote_elem neighbors
1298  if (elem->active())
1299  for (auto n : elem->neighbor_ptr_range())
1300  libmesh_assert_not_equal_to (n, remote_elem);
1301 
1302 #ifdef LIBMESH_ENABLE_AMR
1303  const Elem * parent = elem->parent();
1304  if (parent)
1305  libmesh_assert_not_equal_to (parent, remote_elem);
1306 
1307  // We can only be strict about active elements' subactive
1308  // children
1309  if (elem->active() && elem->has_children())
1310  for (auto & child : elem->child_ref_range())
1311  libmesh_assert_not_equal_to (&child, remote_elem);
1312 #endif
1313  }
1314 }
1315 
1316 
1317 
1319 {
1320  LOG_SCOPE("libmesh_assert_valid_elem_ids()", "MeshTools");
1321 
1322  processor_id_type lastprocid = 0;
1323  dof_id_type lastelemid = 0;
1324 
1325  for (const auto & elem : mesh.active_element_ptr_range())
1326  {
1327  libmesh_assert (elem);
1328  processor_id_type elemprocid = elem->processor_id();
1329  dof_id_type elemid = elem->id();
1330 
1331  libmesh_assert_greater_equal (elemid, lastelemid);
1332  libmesh_assert_greater_equal (elemprocid, lastprocid);
1333 
1334  lastelemid = elemid;
1335  lastprocid = elemprocid;
1336  }
1337 }
1338 
1339 
1340 
1342 {
1343  LOG_SCOPE("libmesh_assert_valid_amr_elem_ids()", "MeshTools");
1344 
1345  for (const auto & elem : mesh.element_ptr_range())
1346  {
1347  libmesh_assert (elem);
1348 
1349  const Elem * parent = elem->parent();
1350 
1351  if (parent)
1352  {
1353  libmesh_assert_greater_equal (elem->id(), parent->id());
1354  libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1355  }
1356  }
1357 }
1358 
1359 
1360 
1362 {
1363  LOG_SCOPE("libmesh_assert_valid_amr_interior_parents()", "MeshTools");
1364 
1365  for (const auto & elem : mesh.element_ptr_range())
1366  {
1367  libmesh_assert (elem);
1368 
1369  // We can skip to the next element if we're full-dimension
1370  // and therefore don't have any interior parents
1371  if (elem->dim() >= LIBMESH_DIM)
1372  continue;
1373 
1374  const Elem * ip = elem->interior_parent();
1375 
1376  const Elem * parent = elem->parent();
1377 
1378  if (ip && (ip != remote_elem) && parent)
1379  {
1380  libmesh_assert_equal_to (ip->top_parent(),
1381  elem->top_parent()->interior_parent());
1382 
1383  if (ip->level() == elem->level())
1384  libmesh_assert_equal_to (ip->parent(),
1385  parent->interior_parent());
1386  else
1387  {
1388  libmesh_assert_less (ip->level(), elem->level());
1389  libmesh_assert_equal_to (ip, parent->interior_parent());
1390  }
1391  }
1392  }
1393 }
1394 
1395 
1396 
1397 void libmesh_assert_contiguous_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1398 {
1399  LOG_SCOPE("libmesh_assert_contiguous_dof_ids()", "MeshTools");
1400 
1401  if (mesh.n_processors() == 1)
1402  return;
1403 
1404  libmesh_parallel_only(mesh.comm());
1405 
1406  dof_id_type min_dof_id = std::numeric_limits<dof_id_type>::max(),
1407  max_dof_id = std::numeric_limits<dof_id_type>::min();
1408 
1409  // Figure out what our local dof id range is
1410  for (const auto * node : mesh.local_node_ptr_range())
1411  {
1412  for (auto v : make_range(node->n_vars(sysnum)))
1413  for (auto c : make_range(node->n_comp(sysnum, v)))
1414  {
1415  dof_id_type id = node->dof_number(sysnum, v, c);
1416  min_dof_id = std::min (min_dof_id, id);
1417  max_dof_id = std::max (max_dof_id, id);
1418  }
1419  }
1420 
1421  // Make sure no other processors' ids are inside it
1422  for (const auto * node : mesh.node_ptr_range())
1423  {
1424  if (node->processor_id() == mesh.processor_id())
1425  continue;
1426  for (auto v : make_range(node->n_vars(sysnum)))
1427  for (auto c : make_range(node->n_comp(sysnum, v)))
1428  {
1429  dof_id_type id = node->dof_number(sysnum, v, c);
1430  libmesh_assert (id < min_dof_id ||
1431  id > max_dof_id);
1432  }
1433  }
1434 }
1435 
1436 
1437 
1438 template <>
1440 {
1441  LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1442 
1443  // This parameter is not used when !LIBMESH_ENABLE_AMR
1445 
1446  // If we're adaptively refining, check processor ids for consistency
1447  // between parents and children.
1448 #ifdef LIBMESH_ENABLE_AMR
1449 
1450  // Ancestor elements we won't worry about, but subactive and active
1451  // elements ought to have parents with consistent processor ids
1452  for (const auto & elem : mesh.element_ptr_range())
1453  {
1454  libmesh_assert(elem);
1455 
1456  if (!elem->active() && !elem->subactive())
1457  continue;
1458 
1459  const Elem * parent = elem->parent();
1460 
1461  if (parent)
1462  {
1463  libmesh_assert(parent->has_children());
1464  processor_id_type parent_procid = parent->processor_id();
1465  bool matching_child_id = false;
1466  // If we've got a remote_elem then we don't know whether
1467  // it's responsible for the parent's processor id; all
1468  // we can do is assume it is and let its processor fail
1469  // an assert if there's something wrong.
1470  for (auto & child : parent->child_ref_range())
1471  if (&child == remote_elem ||
1472  child.processor_id() == parent_procid)
1473  matching_child_id = true;
1474  libmesh_assert(matching_child_id);
1475  }
1476  }
1477 #endif
1478 }
1479 
1480 
1481 
1482 template <>
1484 {
1485  LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1486 
1487  if (mesh.n_processors() == 1)
1488  return;
1489 
1490  libmesh_parallel_only(mesh.comm());
1491 
1492  // We want this test to be valid even when called after nodes have
1493  // been added asynchronously but before they're renumbered.
1494  //
1495  // Plus, some code (looking at you, stitch_meshes) modifies
1496  // DofObject ids without keeping max_elem_id()/max_node_id()
1497  // consistent, but that's done in a safe way for performance
1498  // reasons, so we'll play along and just figure out new max ids
1499  // ourselves.
1500  dof_id_type parallel_max_node_id = 0;
1501  for (const auto & node : mesh.node_ptr_range())
1502  parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
1503  node->id()+1);
1504  mesh.comm().max(parallel_max_node_id);
1505 
1506 
1507  std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1508 
1509  for (const auto & elem : as_range(mesh.local_elements_begin(),
1510  mesh.local_elements_end()))
1511  {
1512  libmesh_assert (elem);
1513 
1514  for (auto & node : elem->node_ref_range())
1515  {
1516  dof_id_type nodeid = node.id();
1517  node_touched_by_me[nodeid] = true;
1518  }
1519  }
1520  std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1521  mesh.comm().max(node_touched_by_anyone);
1522 
1523  for (const auto & node : mesh.local_node_ptr_range())
1524  {
1525  libmesh_assert(node);
1526  dof_id_type nodeid = node->id();
1527  libmesh_assert(!node_touched_by_anyone[nodeid] ||
1528  node_touched_by_me[nodeid]);
1529  }
1530 }
1531 
1532 
1533 
1535 {
1536  for (const auto & elem : mesh.active_element_ptr_range())
1537  for (auto & node : elem->node_ref_range())
1538  libmesh_assert_equal_to
1539  (node.processor_id(),
1540  node.choose_processor_id(node.processor_id(),
1541  elem->processor_id()));
1542 }
1543 
1544 
1545 
1546 #ifdef LIBMESH_ENABLE_AMR
1548 {
1549  LOG_SCOPE("libmesh_assert_valid_refinement_tree()", "MeshTools");
1550 
1551  for (const auto & elem : mesh.element_ptr_range())
1552  {
1553  libmesh_assert(elem);
1554  if (elem->has_children())
1555  for (auto & child : elem->child_ref_range())
1556  if (&child != remote_elem)
1557  libmesh_assert_equal_to (child.parent(), elem);
1558  if (elem->active())
1559  {
1560  libmesh_assert(!elem->ancestor());
1561  libmesh_assert(!elem->subactive());
1562  }
1563  else if (elem->ancestor())
1564  {
1565  libmesh_assert(!elem->subactive());
1566  }
1567  else
1568  libmesh_assert(elem->subactive());
1569 
1570  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1571  libmesh_assert_greater(elem->p_level(), 0);
1572  }
1573 }
1574 #else
1576 {
1577 }
1578 #endif // LIBMESH_ENABLE_AMR
1579 
1580 #endif // !NDEBUG
1581 
1582 
1583 
1584 #ifdef DEBUG
1585 
1587  const Elem * bad_elem)
1588 {
1589  for (const auto & elem : mesh.element_ptr_range())
1590  {
1591  libmesh_assert (elem);
1592  libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1593  for (auto n : elem->neighbor_ptr_range())
1594  libmesh_assert_not_equal_to (n, bad_elem);
1595 
1596 #ifdef LIBMESH_ENABLE_AMR
1597  if (elem->has_children())
1598  for (auto & child : elem->child_ref_range())
1599  libmesh_assert_not_equal_to (&child, bad_elem);
1600 #endif
1601  }
1602 }
1603 
1604 
1606 {
1607  LOG_SCOPE("libmesh_assert_equal_points()", "MeshTools");
1608 
1609  dof_id_type pmax_node_id = mesh.max_node_id();
1610  mesh.comm().max(pmax_node_id);
1611 
1612  for (dof_id_type i=0; i != pmax_node_id; ++i)
1613  {
1614  const Point * p = mesh.query_node_ptr(i);
1615 
1617  }
1618 }
1619 
1620 
1622 {
1623  LOG_SCOPE("libmesh_assert_equal_connectivity()", "MeshTools");
1624 
1625  dof_id_type pmax_elem_id = mesh.max_elem_id();
1626  mesh.comm().max(pmax_elem_id);
1627 
1628  for (dof_id_type i=0; i != pmax_elem_id; ++i)
1629  {
1630  const Elem * e = mesh.query_elem_ptr(i);
1631 
1632  std::vector<dof_id_type> nodes;
1633  if (e)
1634  for (auto n : e->node_index_range())
1635  nodes.push_back(e->node_id(n));
1636 
1637  libmesh_assert(mesh.comm().semiverify(e ? &nodes : nullptr));
1638  }
1639 }
1640 
1641 
1643 {
1644  LOG_SCOPE("libmesh_assert_connected_nodes()", "MeshTools");
1645 
1646  std::set<const Node *> used_nodes;
1647 
1648  for (const auto & elem : mesh.element_ptr_range())
1649  {
1650  libmesh_assert (elem);
1651 
1652  for (auto & n : elem->node_ref_range())
1653  used_nodes.insert(&n);
1654  }
1655 
1656  for (const auto & node : mesh.node_ptr_range())
1657  {
1658  libmesh_assert(node);
1659  libmesh_assert(used_nodes.count(node));
1660  }
1661 }
1662 
1663 
1664 
1666 {
1667  libmesh_parallel_only(mesh.comm());
1668 
1669  const auto & constraint_rows = mesh.get_constraint_rows();
1670 
1671  bool have_constraint_rows = !constraint_rows.empty();
1672  mesh.comm().max(have_constraint_rows);
1673  if (!have_constraint_rows)
1674  return;
1675 
1676  for (auto & row : constraint_rows)
1677  {
1678  const Node * node = row.first;
1679  libmesh_assert(node == mesh.node_ptr(node->id()));
1680 
1681  for (auto & pr : row.second)
1682  {
1683  const Elem * spline_elem = pr.first.first;
1684  libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1685  }
1686  }
1687 
1688  dof_id_type pmax_node_id = mesh.max_node_id();
1689  mesh.comm().max(pmax_node_id);
1690 
1691  for (dof_id_type i=0; i != pmax_node_id; ++i)
1692  {
1693  const Node * node = mesh.query_node_ptr(i);
1694 
1695  bool have_constraint = constraint_rows.count(node);
1696 
1697  const std::size_t my_n_constraints = have_constraint ?
1698  libmesh_map_find(constraint_rows, node).size() : std::size_t(-1);
1699  const std::size_t * n_constraints = node ?
1700  &my_n_constraints : nullptr;
1701 
1702  libmesh_assert(mesh.comm().semiverify(n_constraints));
1703  }
1704 }
1705 
1706 
1707 
1709 {
1710  LOG_SCOPE("libmesh_assert_valid_boundary_ids()", "MeshTools");
1711 
1712  if (mesh.n_processors() == 1)
1713  return;
1714 
1715  libmesh_parallel_only(mesh.comm());
1716 
1717  const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1718 
1719  dof_id_type pmax_elem_id = mesh.max_elem_id();
1720  mesh.comm().max(pmax_elem_id);
1721 
1722  for (dof_id_type i=0; i != pmax_elem_id; ++i)
1723  {
1724  const Elem * elem = mesh.query_elem_ptr(i);
1725  const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1726  const unsigned int my_n_edges = elem ? elem->n_edges() : 0;
1727  const unsigned int my_n_sides = elem ? elem->n_sides() : 0;
1728  unsigned int
1729  n_nodes = my_n_nodes,
1730  n_edges = my_n_edges,
1731  n_sides = my_n_sides;
1732 
1733  mesh.comm().max(n_nodes);
1734  mesh.comm().max(n_edges);
1735  mesh.comm().max(n_sides);
1736 
1737  if (elem)
1738  {
1739  libmesh_assert_equal_to(my_n_nodes, n_nodes);
1740  libmesh_assert_equal_to(my_n_edges, n_edges);
1741  libmesh_assert_equal_to(my_n_sides, n_sides);
1742  }
1743 
1744  // Let's test all IDs on the element with one communication
1745  // rather than n_nodes + n_edges + n_sides communications, to
1746  // cut down on latency in dbg modes.
1747  std::vector<boundary_id_type> all_bcids;
1748 
1749  for (unsigned int n=0; n != n_nodes; ++n)
1750  {
1751  std::vector<boundary_id_type> bcids;
1752  if (elem)
1753  {
1754  boundary_info.boundary_ids(elem->node_ptr(n), bcids);
1755 
1756  // Ordering of boundary ids shouldn't matter
1757  std::sort(bcids.begin(), bcids.end());
1758  }
1759  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1760 
1761  all_bcids.insert(all_bcids.end(), bcids.begin(),
1762  bcids.end());
1763  // Separator
1764  all_bcids.push_back(BoundaryInfo::invalid_id);
1765  }
1766 
1767  for (unsigned short e=0; e != n_edges; ++e)
1768  {
1769  std::vector<boundary_id_type> bcids;
1770 
1771  if (elem)
1772  {
1773  boundary_info.edge_boundary_ids(elem, e, bcids);
1774 
1775  // Ordering of boundary ids shouldn't matter
1776  std::sort(bcids.begin(), bcids.end());
1777  }
1778 
1779  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1780 
1781  all_bcids.insert(all_bcids.end(), bcids.begin(),
1782  bcids.end());
1783  // Separator
1784  all_bcids.push_back(BoundaryInfo::invalid_id);
1785 
1786  if (elem)
1787  {
1788  boundary_info.raw_edge_boundary_ids(elem, e, bcids);
1789 
1790  // Ordering of boundary ids shouldn't matter
1791  std::sort(bcids.begin(), bcids.end());
1792 
1793  all_bcids.insert(all_bcids.end(), bcids.begin(),
1794  bcids.end());
1795  // Separator
1796  all_bcids.push_back(BoundaryInfo::invalid_id);
1797  }
1798 
1799  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1800  }
1801 
1802  for (unsigned short s=0; s != n_sides; ++s)
1803  {
1804  std::vector<boundary_id_type> bcids;
1805 
1806  if (elem)
1807  {
1808  boundary_info.boundary_ids(elem, s, bcids);
1809 
1810  // Ordering of boundary ids shouldn't matter
1811  std::sort(bcids.begin(), bcids.end());
1812 
1813  all_bcids.insert(all_bcids.end(), bcids.begin(),
1814  bcids.end());
1815  // Separator
1816  all_bcids.push_back(BoundaryInfo::invalid_id);
1817  }
1818 
1819  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1820 
1821  if (elem)
1822  {
1823  boundary_info.raw_boundary_ids(elem, s, bcids);
1824 
1825  // Ordering of boundary ids shouldn't matter
1826  std::sort(bcids.begin(), bcids.end());
1827 
1828  all_bcids.insert(all_bcids.end(), bcids.begin(),
1829  bcids.end());
1830  // Separator
1831  all_bcids.push_back(BoundaryInfo::invalid_id);
1832  }
1833 
1834  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1835  }
1836 
1837  for (unsigned short sf=0; sf != 2; ++sf)
1838  {
1839  std::vector<boundary_id_type> bcids;
1840 
1841  if (elem)
1842  {
1843  boundary_info.shellface_boundary_ids(elem, sf, bcids);
1844 
1845  // Ordering of boundary ids shouldn't matter
1846  std::sort(bcids.begin(), bcids.end());
1847 
1848  all_bcids.insert(all_bcids.end(), bcids.begin(),
1849  bcids.end());
1850  // Separator
1851  all_bcids.push_back(BoundaryInfo::invalid_id);
1852  }
1853 
1854  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1855 
1856  if (elem)
1857  {
1858  boundary_info.raw_shellface_boundary_ids(elem, sf, bcids);
1859 
1860  // Ordering of boundary ids shouldn't matter
1861  std::sort(bcids.begin(), bcids.end());
1862 
1863  all_bcids.insert(all_bcids.end(), bcids.begin(),
1864  bcids.end());
1865  // Separator
1866  all_bcids.push_back(BoundaryInfo::invalid_id);
1867  }
1868 
1869  // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1870  }
1871 
1873  (elem ? &all_bcids : nullptr));
1874  }
1875 }
1876 
1877 
1878 void libmesh_assert_valid_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1879 {
1880  LOG_SCOPE("libmesh_assert_valid_dof_ids()", "MeshTools");
1881 
1882  if (mesh.n_processors() == 1)
1883  return;
1884 
1885  libmesh_parallel_only(mesh.comm());
1886 
1887  dof_id_type pmax_elem_id = mesh.max_elem_id();
1888  mesh.comm().max(pmax_elem_id);
1889 
1890  for (dof_id_type i=0; i != pmax_elem_id; ++i)
1891  assert_semiverify_dofobj(mesh.comm(),
1892  mesh.query_elem_ptr(i),
1893  sysnum);
1894 
1895  dof_id_type pmax_node_id = mesh.max_node_id();
1896  mesh.comm().max(pmax_node_id);
1897 
1898  for (dof_id_type i=0; i != pmax_node_id; ++i)
1899  assert_semiverify_dofobj(mesh.comm(),
1900  mesh.query_node_ptr(i),
1901  sysnum);
1902 }
1903 
1904 
1905 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1907 {
1908  LOG_SCOPE("libmesh_assert_valid_unique_ids()", "MeshTools");
1909 
1910  libmesh_parallel_only(mesh.comm());
1911 
1912  // First collect all the unique_ids we can see and make sure there's
1913  // no duplicates
1914  std::unordered_set<unique_id_type> semilocal_unique_ids;
1915 
1916  for (auto const & elem : mesh.active_element_ptr_range())
1917  {
1918  libmesh_assert (!semilocal_unique_ids.count(elem->unique_id()));
1919  semilocal_unique_ids.insert(elem->unique_id());
1920  }
1921 
1922  for (auto const & node : mesh.node_ptr_range())
1923  {
1924  libmesh_assert (!semilocal_unique_ids.count(node->unique_id()));
1925  semilocal_unique_ids.insert(node->unique_id());
1926  }
1927 
1928  // Then make sure elements are all in sync and remote elements don't
1929  // duplicate semilocal
1930 
1931  dof_id_type pmax_elem_id = mesh.max_elem_id();
1932  mesh.comm().max(pmax_elem_id);
1933 
1934  for (dof_id_type i=0; i != pmax_elem_id; ++i)
1935  {
1936  const Elem * elem = mesh.query_elem_ptr(i);
1937  assert_dofobj_unique_id(mesh.comm(), elem, semilocal_unique_ids);
1938  }
1939 
1940  // Then make sure nodes are all in sync and remote elements don't
1941  // duplicate semilocal
1942 
1943  dof_id_type pmax_node_id = mesh.max_node_id();
1944  mesh.comm().max(pmax_node_id);
1945 
1946  for (dof_id_type i=0; i != pmax_node_id; ++i)
1947  {
1948  const Node * node = mesh.query_node_ptr(i);
1949  assert_dofobj_unique_id(mesh.comm(), node, semilocal_unique_ids);
1950  }
1951 }
1952 #endif
1953 
1955 {
1956  libmesh_parallel_only(mesh.comm());
1957 
1958  dof_id_type parallel_max_elem_id = mesh.max_elem_id();
1959  mesh.comm().max(parallel_max_elem_id);
1960 
1961  for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
1962  {
1963  const Elem * elem = mesh.query_elem_ptr(i);
1964  processor_id_type pid =
1966  mesh.comm().min(pid);
1967  libmesh_assert(elem || pid != mesh.processor_id());
1968  }
1969 
1970  dof_id_type parallel_max_node_id = mesh.max_node_id();
1971  mesh.comm().max(parallel_max_node_id);
1972 
1973  for (dof_id_type i=0; i != parallel_max_node_id; ++i)
1974  {
1975  const Node * node = mesh.query_node_ptr(i);
1976  processor_id_type pid =
1978  mesh.comm().min(pid);
1979  libmesh_assert(node || pid != mesh.processor_id());
1980  }
1981 }
1982 
1983 
1985 {
1986  libmesh_parallel_only(mesh.comm());
1987  auto locator = mesh.sub_point_locator();
1988 
1989  dof_id_type parallel_max_elem_id = mesh.max_elem_id();
1990  mesh.comm().max(parallel_max_elem_id);
1991 
1992  for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
1993  {
1994  const Elem * elem = mesh.query_elem_ptr(i);
1995 
1996  const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1997  unsigned int n_nodes = my_n_nodes;
1998  mesh.comm().max(n_nodes);
1999 
2000  if (n_nodes)
2001  libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2002 
2003  for (unsigned int n=0; n != n_nodes; ++n)
2004  {
2005  const Node * node = elem ? elem->node_ptr(n) : nullptr;
2006  processor_id_type pid =
2008  mesh.comm().min(pid);
2009  libmesh_assert(node || pid != mesh.processor_id());
2010  }
2011  }
2012 }
2013 
2014 
2015 
2016 template <>
2018 {
2019  LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2020 
2021  if (mesh.n_processors() == 1)
2022  return;
2023 
2024  libmesh_parallel_only(mesh.comm());
2025 
2026  // Some code (looking at you, stitch_meshes) modifies DofObject ids
2027  // without keeping max_elem_id()/max_node_id() consistent, but
2028  // that's done in a safe way for performance reasons, so we'll play
2029  // along and just figure out new max ids ourselves.
2030  dof_id_type parallel_max_elem_id = 0;
2031  for (const auto & elem : mesh.element_ptr_range())
2032  parallel_max_elem_id = std::max<dof_id_type>(parallel_max_elem_id,
2033  elem->id()+1);
2034  mesh.comm().max(parallel_max_elem_id);
2035 
2036  // Check processor ids for consistency between processors
2037 
2038  for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2039  {
2040  const Elem * elem = mesh.query_elem_ptr(i);
2041 
2042  processor_id_type min_id =
2043  elem ? elem->processor_id() :
2044  std::numeric_limits<processor_id_type>::max();
2045  mesh.comm().min(min_id);
2046 
2047  processor_id_type max_id =
2048  elem ? elem->processor_id() :
2049  std::numeric_limits<processor_id_type>::min();
2050  mesh.comm().max(max_id);
2051 
2052  if (elem)
2053  {
2054  libmesh_assert_equal_to (min_id, elem->processor_id());
2055  libmesh_assert_equal_to (max_id, elem->processor_id());
2056  }
2057 
2058  if (min_id == mesh.processor_id())
2059  libmesh_assert(elem);
2060  }
2061 }
2062 
2063 
2064 
2066 {
2067  LOG_SCOPE("libmesh_assert_parallel_consistent_new_node_procids()", "MeshTools");
2068 
2069  if (mesh.n_processors() == 1)
2070  return;
2071 
2072  libmesh_parallel_only(mesh.comm());
2073 
2074  // We want this test to hit every node when called even after nodes
2075  // have been added asynchronously but before everything has been
2076  // renumbered.
2077  dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2078  mesh.comm().max(parallel_max_elem_id);
2079 
2080  std::vector<bool> elem_touched_by_anyone(parallel_max_elem_id, false);
2081 
2082  for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2083  {
2084  const Elem * elem = mesh.query_elem_ptr(i);
2085 
2086  const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2087  unsigned int n_nodes = my_n_nodes;
2088  mesh.comm().max(n_nodes);
2089 
2090  if (n_nodes)
2091  libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2092 
2093  for (unsigned int n=0; n != n_nodes; ++n)
2094  {
2095  const Node * node = elem ? elem->node_ptr(n) : nullptr;
2096  const processor_id_type pid = node ? node->processor_id() : 0;
2097  libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2098  }
2099  }
2100 }
2101 
2102 template <>
2104 {
2105  LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2106 
2107  if (mesh.n_processors() == 1)
2108  return;
2109 
2110  libmesh_parallel_only(mesh.comm());
2111 
2112  // We want this test to be valid even when called even after nodes
2113  // have been added asynchronously but before they're renumbered
2114  //
2115  // Plus, some code (looking at you, stitch_meshes) modifies
2116  // DofObject ids without keeping max_elem_id()/max_node_id()
2117  // consistent, but that's done in a safe way for performance
2118  // reasons, so we'll play along and just figure out new max ids
2119  // ourselves.
2120  dof_id_type parallel_max_node_id = 0;
2121  for (const auto & node : mesh.node_ptr_range())
2122  parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
2123  node->id()+1);
2124  mesh.comm().max(parallel_max_node_id);
2125 
2126  std::vector<bool> node_touched_by_anyone(parallel_max_node_id, false);
2127 
2128  for (const auto & elem : as_range(mesh.local_elements_begin(),
2129  mesh.local_elements_end()))
2130  {
2131  libmesh_assert (elem);
2132 
2133  for (auto & node : elem->node_ref_range())
2134  {
2135  dof_id_type nodeid = node.id();
2136  node_touched_by_anyone[nodeid] = true;
2137  }
2138  }
2139  mesh.comm().max(node_touched_by_anyone);
2140 
2141  // Check processor ids for consistency between processors
2142  // on any node an element touches
2143  for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2144  {
2145  if (!node_touched_by_anyone[i])
2146  continue;
2147 
2148  const Node * node = mesh.query_node_ptr(i);
2149  const processor_id_type pid = node ? node->processor_id() : 0;
2150 
2151  libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2152  }
2153 }
2154 
2155 
2156 
2157 #ifdef LIBMESH_ENABLE_AMR
2159 {
2160  LOG_SCOPE("libmesh_assert_valid_refinement_flags()", "MeshTools");
2161 
2162  libmesh_parallel_only(mesh.comm());
2163  if (mesh.n_processors() == 1)
2164  return;
2165 
2166  dof_id_type pmax_elem_id = mesh.max_elem_id();
2167  mesh.comm().max(pmax_elem_id);
2168 
2169  std::vector<unsigned char> my_elem_h_state(pmax_elem_id, 255);
2170  std::vector<unsigned char> my_elem_p_state(pmax_elem_id, 255);
2171 
2172  for (const auto & elem : mesh.element_ptr_range())
2173  {
2174  libmesh_assert (elem);
2175  dof_id_type elemid = elem->id();
2176 
2177  my_elem_h_state[elemid] =
2178  static_cast<unsigned char>(elem->refinement_flag());
2179 
2180  my_elem_p_state[elemid] =
2181  static_cast<unsigned char>(elem->p_refinement_flag());
2182  }
2183  std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
2184  mesh.comm().min(min_elem_h_state);
2185 
2186  std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
2187  mesh.comm().min(min_elem_p_state);
2188 
2189  for (dof_id_type i=0; i!= pmax_elem_id; ++i)
2190  {
2191  libmesh_assert(my_elem_h_state[i] == 255 ||
2192  my_elem_h_state[i] == min_elem_h_state[i]);
2193  libmesh_assert(my_elem_p_state[i] == 255 ||
2194  my_elem_p_state[i] == min_elem_p_state[i]);
2195  }
2196 }
2197 #else
2199 {
2200 }
2201 #endif // LIBMESH_ENABLE_AMR
2202 
2203 
2204 
2206  bool assert_valid_remote_elems)
2207 {
2208  LOG_SCOPE("libmesh_assert_valid_neighbors()", "MeshTools");
2209 
2210  for (const auto & elem : mesh.element_ptr_range())
2211  {
2212  libmesh_assert (elem);
2213  elem->libmesh_assert_valid_neighbors();
2214  }
2215 
2216  if (mesh.n_processors() == 1)
2217  return;
2218 
2219  libmesh_parallel_only(mesh.comm());
2220 
2221  dof_id_type pmax_elem_id = mesh.max_elem_id();
2222  mesh.comm().max(pmax_elem_id);
2223 
2224  for (dof_id_type i=0; i != pmax_elem_id; ++i)
2225  {
2226  const Elem * elem = mesh.query_elem_ptr(i);
2227 
2228  const unsigned int my_n_neigh = elem ? elem->n_neighbors() : 0;
2229  unsigned int n_neigh = my_n_neigh;
2230  mesh.comm().max(n_neigh);
2231  if (elem)
2232  libmesh_assert_equal_to (my_n_neigh, n_neigh);
2233 
2234  for (unsigned int n = 0; n != n_neigh; ++n)
2235  {
2236  dof_id_type my_neighbor = DofObject::invalid_id;
2237  dof_id_type * p_my_neighbor = nullptr;
2238 
2239  // If we have a non-remote_elem neighbor link, then we can
2240  // verify it.
2241  if (elem && elem->neighbor_ptr(n) != remote_elem)
2242  {
2243  p_my_neighbor = &my_neighbor;
2244  if (elem->neighbor_ptr(n))
2245  my_neighbor = elem->neighbor_ptr(n)->id();
2246 
2247  // But wait - if we haven't set remote_elem links yet then
2248  // some nullptr links on ghost elements might be
2249  // future-remote_elem links, so we can't verify those.
2250  if (!assert_valid_remote_elems &&
2251  !elem->neighbor_ptr(n) &&
2252  elem->processor_id() != mesh.processor_id())
2253  p_my_neighbor = nullptr;
2254  }
2255  libmesh_assert(mesh.comm().semiverify(p_my_neighbor));
2256  }
2257  }
2258 }
2259 #endif // DEBUG
2260 
2261 
2262 
2263 // Functors for correct_node_proc_ids
2264 namespace {
2265 
2266 typedef std::unordered_map<dof_id_type, processor_id_type> proc_id_map_type;
2267 
2268 struct SyncNodeSet
2269 {
2270  typedef unsigned char datum; // bool but without bit twiddling issues
2271 
2272  SyncNodeSet(std::unordered_set<const Node *> & _set,
2273  MeshBase & _mesh) :
2274  node_set(_set), mesh(_mesh) {}
2275 
2276  std::unordered_set<const Node *> & node_set;
2277 
2279 
2280  // ------------------------------------------------------------
2281  void gather_data (const std::vector<dof_id_type> & ids,
2282  std::vector<datum> & data)
2283  {
2284  // Find whether each requested node belongs in the set
2285  data.resize(ids.size());
2286 
2287  for (auto i : index_range(ids))
2288  {
2289  const dof_id_type id = ids[i];
2290 
2291  // We'd better have every node we're asked for
2292  Node * node = mesh.node_ptr(id);
2293 
2294  // Return if the node is in the set.
2295  data[i] = node_set.count(node);
2296  }
2297  }
2298 
2299  // ------------------------------------------------------------
2300  bool act_on_data (const std::vector<dof_id_type> & ids,
2301  const std::vector<datum> in_set)
2302  {
2303  bool data_changed = false;
2304 
2305  // Add nodes we've been informed of to our own set
2306  for (auto i : index_range(ids))
2307  {
2308  if (in_set[i])
2309  {
2310  Node * node = mesh.node_ptr(ids[i]);
2311  if (!node_set.count(node))
2312  {
2313  node_set.insert(node);
2314  data_changed = true;
2315  }
2316  }
2317  }
2318 
2319  return data_changed;
2320  }
2321 };
2322 
2323 
2324 struct NodesNotInSet
2325 {
2326  NodesNotInSet(const std::unordered_set<const Node *> _set)
2327  : node_set(_set) {}
2328 
2329  bool operator() (const Node * node) const
2330  {
2331  if (node_set.count(node))
2332  return false;
2333  return true;
2334  }
2335 
2336  const std::unordered_set<const Node *> node_set;
2337 };
2338 
2339 
2340 struct SyncProcIdsFromMap
2341 {
2342  typedef processor_id_type datum;
2343 
2344  SyncProcIdsFromMap(const proc_id_map_type & _map,
2345  MeshBase & _mesh) :
2346  new_proc_ids(_map), mesh(_mesh) {}
2347 
2348  const proc_id_map_type & new_proc_ids;
2349 
2350  MeshBase & mesh;
2351 
2352  // ------------------------------------------------------------
2353  void gather_data (const std::vector<dof_id_type> & ids,
2354  std::vector<datum> & data)
2355  {
2356  // Find the new processor id of each requested node
2357  data.resize(ids.size());
2358 
2359  for (auto i : index_range(ids))
2360  {
2361  const dof_id_type id = ids[i];
2362 
2363  // Return the node's new processor id if it has one, or its
2364  // old processor id if not.
2365  if (const auto it = new_proc_ids.find(id);
2366  it != new_proc_ids.end())
2367  data[i] = it->second;
2368  else
2369  {
2370  // We'd better find every node we're asked for
2371  const Node & node = mesh.node_ref(id);
2372  data[i] = node.processor_id();
2373  }
2374  }
2375  }
2376 
2377  // ------------------------------------------------------------
2378  void act_on_data (const std::vector<dof_id_type> & ids,
2379  const std::vector<datum> proc_ids)
2380  {
2381  // Set the node processor ids we've now been informed of
2382  for (auto i : index_range(ids))
2383  {
2384  Node & node = mesh.node_ref(ids[i]);
2385  node.processor_id() = proc_ids[i];
2386  }
2387  }
2388 };
2389 }
2390 
2391 
2392 
2394 {
2395  LOG_SCOPE("correct_node_proc_ids()","MeshTools");
2396 
2397  // This function must be run on all processors at once
2398  libmesh_parallel_only(mesh.comm());
2399 
2400  // We require all processors to agree on nodal processor ids before
2401  // going into this algorithm.
2402 #ifdef DEBUG
2404 #endif
2405 
2406  // If we have any unpartitioned elements at this
2407  // stage there is a problem
2408  libmesh_assert (n_elem(mesh.unpartitioned_elements_begin(),
2409  mesh.unpartitioned_elements_end()) == 0);
2410 
2411  // Fix nodes' processor ids. Coarsening may have left us with nodes
2412  // which are no longer touched by any elements of the same processor
2413  // id, and for DofMap to work we need to fix that.
2414 
2415  // This is harder now that libMesh no longer requires a distributed
2416  // mesh to ghost all nodal neighbors: it is possible for two active
2417  // elements on two different processors to share the same node in
2418  // such a way that neither processor knows the others' element
2419  // exists!
2420 
2421  // While we're at it, if this mesh is configured to allow
2422  // repartitioning, we'll repartition *all* the nodes' processor ids
2423  // using the canonical Node heuristic, to try and improve DoF load
2424  // balancing. But if the mesh is disallowing repartitioning, we
2425  // won't touch processor_id on any node where it's valid, regardless
2426  // of whether or not it's canonical.
2427  bool repartition_all_nodes = !mesh.skip_noncritical_partitioning();
2428  std::unordered_set<const Node *> valid_nodes;
2429 
2430  // If we aren't allowed to repartition, then we're going to leave
2431  // every node we can at its current processor_id, and *only*
2432  // repartition the nodes whose current processor id is incompatible
2433  // with DoFMap (because it doesn't touch an active element, e.g. due
2434  // to coarsening)
2435  if (!repartition_all_nodes)
2436  {
2437  for (const auto & elem : mesh.active_element_ptr_range())
2438  for (const auto & node : elem->node_ref_range())
2439  if (elem->processor_id() == node.processor_id())
2440  valid_nodes.insert(&node);
2441 
2442  SyncNodeSet syncv(valid_nodes, mesh);
2443 
2445  (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), syncv);
2446  }
2447 
2448  // We build up a set of compatible processor ids for each node
2449  proc_id_map_type new_proc_ids;
2450 
2451  for (auto & elem : mesh.active_element_ptr_range())
2452  {
2453  processor_id_type pid = elem->processor_id();
2454 
2455  for (auto & node : elem->node_ref_range())
2456  {
2457  const dof_id_type id = node.id();
2458  if (auto it = new_proc_ids.find(id);
2459  it == new_proc_ids.end())
2460  new_proc_ids.emplace(id, pid);
2461  else
2462  it->second = node.choose_processor_id(it->second, pid);
2463  }
2464  }
2465 
2466  // Sort the new pids to push to each processor
2467  std::map<processor_id_type, std::vector<std::pair<dof_id_type, processor_id_type>>>
2468  ids_to_push;
2469 
2470  for (const auto & node : mesh.node_ptr_range())
2471  if (const auto it = std::as_const(new_proc_ids).find(node->id());
2472  it != new_proc_ids.end() && node->processor_id() != DofObject::invalid_processor_id)
2473  ids_to_push[node->processor_id()].emplace_back(node->id(), /*pid=*/it->second);
2474 
2475  auto action_functor =
2476  [& mesh, & new_proc_ids]
2478  const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
2479  {
2480  for (const auto & [id, pid] : data)
2481  {
2482  if (const auto it = new_proc_ids.find(id);
2483  it == new_proc_ids.end())
2484  new_proc_ids.emplace(id, pid);
2485  else
2486  {
2487  const Node & node = mesh.node_ref(id);
2488  it->second = node.choose_processor_id(it->second, pid);
2489  }
2490  }
2491  };
2492 
2493  Parallel::push_parallel_vector_data
2494  (mesh.comm(), ids_to_push, action_functor);
2495 
2496  // Now new_proc_ids is correct for every node we used to own. Let's
2497  // ask every other processor about the nodes they used to own. But
2498  // first we'll need to keep track of which nodes we used to own,
2499  // lest we get them confused with nodes we newly own.
2500  std::unordered_set<Node *> ex_local_nodes;
2501  for (auto & node : mesh.local_node_ptr_range())
2502  if (const auto it = new_proc_ids.find(node->id());
2503  it != new_proc_ids.end() && it->second != mesh.processor_id())
2504  ex_local_nodes.insert(node);
2505 
2506  SyncProcIdsFromMap sync(new_proc_ids, mesh);
2507  if (repartition_all_nodes)
2509  (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync);
2510  else
2511  {
2512  NodesNotInSet nnis(valid_nodes);
2513 
2515  (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), nnis, sync);
2516  }
2517 
2518  // And finally let's update the nodes we used to own.
2519  for (const auto & node : ex_local_nodes)
2520  {
2521  if (valid_nodes.count(node))
2522  continue;
2523 
2524  const dof_id_type id = node->id();
2525  const proc_id_map_type::iterator it = new_proc_ids.find(id);
2526  libmesh_assert(it != new_proc_ids.end());
2527  node->processor_id() = it->second;
2528  }
2529 
2530  // We should still have consistent nodal processor ids coming out of
2531  // this algorithm, but if we're allowed to repartition the mesh then
2532  // they should be canonically correct too.
2533 #ifdef DEBUG
2534  libmesh_assert_valid_procids<Node>(mesh);
2535  //if (repartition_all_nodes)
2536  // libmesh_assert_canonical_node_procids(mesh);
2537 #endif
2538 }
2539 
2540 
2541 
2543 {
2545 }
2546 
2547 } // namespace MeshTools
2548 
2549 } // namespace libMesh
unsigned int n_active_local_levels(const MeshBase &mesh)
Definition: mesh_tools.C:760
void libmesh_assert_equal_points(const MeshBase &mesh)
A function for testing that node locations match across processors.
Definition: mesh_tools.C:1605
ElemType
Defines an enum for geometric element types.
void libmesh_assert_topology_consistent_procids< Node >(const MeshBase &mesh)
Definition: mesh_tools.C:1483
void elem_types(const MeshBase &mesh, std::vector< ElemType > &et)
Fills in a vector of all element types in the mesh.
Definition: mesh_tools.C:716
dof_id_type n_nodes(const MeshBase::const_node_iterator &begin, const MeshBase::const_node_iterator &end)
Count up the number of nodes of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:977
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
const Elem * parent() const
Definition: elem.h:3030
bool is_prepared() const
Definition: mesh_base.h:198
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
void libmesh_assert_valid_amr_interior_parents(const MeshBase &mesh)
A function for verifying that any interior_parent pointers on elements are consistent with AMR (paren...
Definition: mesh_tools.C:1361
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:1638
This class defines a sphere.
Definition: sphere.h:72
void libmesh_assert_valid_remote_elems(const MeshBase &mesh)
A function for verifying that active local elements&#39; neighbors are never remote elements.
Definition: mesh_tools.C:1287
void libmesh_assert_consistent_distributed_nodes(const MeshBase &mesh)
A function for verifying that distribution of nodes is parallel consistent (every processor can see e...
Definition: mesh_tools.C:1984
const Elem * interior_parent() const
Definition: elem.C:1186
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
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
libMesh::BoundingBox create_bounding_box(const MeshBase &mesh)
Definition: mesh_tools.C:559
The definition of the const_element_iterator struct.
Definition: mesh_base.h:2216
bool valid_is_prepared(const MeshBase &mesh)
A function for testing whether a mesh&#39;s cached is_prepared() setting is not a false positive...
Definition: mesh_tools.C:1182
unsigned int dim
const Elem * top_parent() const
Definition: elem.h:3056
Sphere processor_bounding_sphere(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:663
Dummy "splitting object" used to distinguish splitting constructors from copy constructors.
Definition: threads_none.h:63
void clear_spline_nodes(MeshBase &)
Remove spline node (for IsoGeometric Analysis meshes) elements and nodes and constraints from the mes...
Definition: mesh_tools.C:1147
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
void sum(T &r) const
std::unordered_set< dof_id_type > find_block_boundary_nodes(const MeshBase &mesh)
Returns a std::set containing Node IDs for all of the block boundary nodes.
Definition: mesh_tools.C:537
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
void assign_global_indices(MeshBase &) const
This method assigns globally unique, partition-agnostic indices to the nodes and elements in the mesh...
void find_nodal_neighbors(const MeshBase &mesh, const Node &n, const std::vector< std::vector< const Elem *>> &nodes_to_elem_map, std::vector< const Node *> &neighbors)
Given a mesh and a node in the mesh, the vector will be filled with every node directly attached to t...
Definition: mesh_tools.C:1036
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
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
unsigned int n_local_levels(const MeshBase &mesh)
Definition: mesh_tools.C:789
void skip_noncritical_partitioning(bool skip)
If true is passed in then the elements on this mesh will no longer be (re)partitioned, and the nodes on this mesh will only be repartitioned if they are found "orphaned" via coarsening or other removal of the last element responsible for their node/element processor id consistency.
Definition: mesh_base.h:1236
unique_id_type unique_id() const
Definition: dof_object.h:844
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
void build_nodes_to_elem_map(const MeshBase &mesh, std::vector< std::vector< dof_id_type >> &nodes_to_elem_map)
After calling this function the input vector nodes_to_elem_map will contain the node to element conne...
Definition: mesh_tools.C:449
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
The libMesh namespace provides an interface to certain functionality in the library.
MeshBase & mesh
Definition: mesh_tools.C:2278
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
virtual std::unique_ptr< MeshBase > clone() const =0
Virtual "copy constructor".
Real distance(const Point &p)
dof_id_type n_elem_of_type(const MeshBase &mesh, const ElemType type)
Definition: mesh_tools.C:728
processor_id_type choose_processor_id(processor_id_type pid1, processor_id_type pid2) const
Return which of pid1 and pid2 would be preferred by the current load-balancing heuristic applied to t...
Definition: node.C:78
libMesh::BoundingBox create_local_bounding_box(const MeshBase &mesh)
Definition: mesh_tools.C:624
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
void libmesh_assert_consistent_distributed(const MeshBase &mesh)
A function for verifying that distribution of dof objects is parallel consistent (every processor can...
Definition: mesh_tools.C:1954
unsigned int paranoid_n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:826
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:517
virtual bool is_serial_on_zero() const
Definition: mesh_base.h:218
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
void libmesh_assert_valid_dof_ids(const MeshBase &mesh, unsigned int sysnum=libMesh::invalid_uint)
A function for verifying that degree of freedom indexing matches across processors.
Definition: mesh_tools.C:1878
uint8_t processor_id_type
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
void libmesh_ignore(const Args &...)
const dof_id_type n_nodes
Definition: tecplot_io.C:67
void libmesh_assert_valid_node_pointers(const MeshBase &mesh)
A function for walking across the mesh to try and ferret out invalidated or misassigned pointers...
Definition: mesh_tools.C:1262
static const boundary_id_type invalid_id
Number used for internal use.
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
ElemMappingType mapping_type() const
Definition: elem.h:3120
This is the MeshCommunication class.
void libmesh_assert_canonical_node_procids(const MeshBase &mesh)
A function for verifying that processor assignment of nodes matches the heuristic specified in Node::...
Definition: mesh_tools.C:1534
dof_id_type id() const
Definition: dof_object.h:828
void min(const T &r, T &o, Request &req) const
unsigned int n_vars
virtual unsigned int n_nodes() const =0
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 libmesh_assert_no_links_to_elem(const MeshBase &mesh, const Elem *bad_elem)
The following functions, only available in builds with DEBUG defined, typically have surprisingly slo...
Definition: mesh_tools.C:1586
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:437
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:967
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition: elem.h:2919
unsigned int n_systems() const
Definition: dof_object.h:937
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
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
virtual dof_id_type max_elem_id() const =0
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
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.
void libmesh_assert_valid_amr_elem_ids(const MeshBase &mesh)
A function for verifying that ids of elements are correctly sorted for AMR (parents have lower ids th...
Definition: mesh_tools.C:1341
void libmesh_assert_parallel_consistent_procids< Elem >(const MeshBase &mesh)
Definition: mesh_tools.C:2017
Sphere subdomain_bounding_sphere(const MeshBase &mesh, const subdomain_id_type sid)
Definition: mesh_tools.C:702
const proc_id_map_type & new_proc_ids
Definition: mesh_tools.C:2348
virtual unsigned int n_edges() const =0
timpi_pure bool semiverify(const T *r) const
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
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 globally_renumber_nodes_and_elements(MeshBase &)
There is no reason for a user to ever call this function.
Definition: mesh_tools.C:2542
Real volume(const MeshBase &mesh, unsigned int dim=libMesh::invalid_uint)
Find the total volume of a mesh (interpreting that as area for dim = 2, or total arc length for dim =...
Definition: mesh_tools.C:985
void find_hanging_nodes_and_parents(const MeshBase &mesh, std::map< dof_id_type, std::vector< dof_id_type >> &hanging_nodes)
Given a mesh hanging_nodes will be filled with an associative array keyed off the global id of all th...
Definition: mesh_tools.C:1059
dof_id_type total_weight(const MeshBase &mesh)
Definition: mesh_tools.C:414
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
void libmesh_assert_parallel_consistent_procids< Node >(const MeshBase &mesh)
Definition: mesh_tools.C:2103
dof_id_type n_connected_components(const MeshBase &mesh, Real constraint_tol=0)
Definition: mesh_tools.C:840
libMesh::BoundingBox create_subdomain_bounding_box(const MeshBase &mesh, const subdomain_id_type sid)
Definition: mesh_tools.C:678
std::string enum_to_string(const T e)
Sphere bounding_sphere(const MeshBase &mesh)
Definition: mesh_tools.C:611
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual unsigned int n_sides() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
void libmesh_assert_equal_connectivity(const MeshBase &mesh)
A function for testing that element nodal connectivities match across processors. ...
Definition: mesh_tools.C:1621
std::unordered_set< const Node * > & node_set
Definition: mesh_tools.C:2276
unsigned int level() const
Definition: elem.h:3074
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
timpi_pure bool verify(const T &r) const
void libmesh_assert_connected_nodes(const MeshBase &mesh)
A function for verifying that all nodes are connected to at least one element.
Definition: mesh_tools.C:1642
void max(const T &r, T &o, Request &req) const
DIE A HORRIBLE DEATH HERE typedef MPI_Comm communicator
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
unsigned int n_neighbors() const
Definition: elem.h:714
void get_not_subactive_node_ids(const MeshBase &mesh, std::set< dof_id_type > &not_subactive_node_ids)
Builds a set of node IDs for nodes which belong to non-subactive elements.
Definition: mesh_tools.C:958
The definition of the const_node_iterator struct.
Definition: mesh_base.h:2267
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
The DofObject defines an abstract base class for objects that have degrees of freedom associated with...
Definition: dof_object.h:54
dof_id_type n_active_elem_of_type(const MeshBase &mesh, const ElemType type)
Definition: mesh_tools.C:737
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
void parallel_reduce(const Range &range, Body &body)
Execute the provided reduction operation in parallel on the specified range.
Definition: threads_none.h:101
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2683
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
void libmesh_assert_valid_neighbors(const MeshBase &mesh, bool assert_valid_remote_elems=true)
A function for verifying that neighbor connectivity is correct (each element is a neighbor of or desc...
Definition: mesh_tools.C:2205
libMesh::BoundingBox create_nodal_bounding_box(const MeshBase &mesh)
Definition: mesh_tools.C:584
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
void libmesh_assert_valid_unique_ids(const MeshBase &mesh)
A function for verifying that unique ids match across processors.
Definition: mesh_tools.C:1906
processor_id_type processor_id() const
libMesh::BoundingBox create_processor_bounding_box(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:638
unsigned int n_p_levels(const MeshBase &mesh)
Definition: mesh_tools.C:1014
processor_id_type processor_id() const
Definition: dof_object.h:905
void libmesh_assert_old_dof_objects(const MeshBase &mesh)
A function for testing that all non-recently-created DofObjects within a mesh have old_dof_object dat...
Definition: mesh_tools.C:1238
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
uint8_t unique_id_type
Definition: id_types.h:86
bool has_children() const
Definition: elem.h:2979
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
void libmesh_assert_contiguous_dof_ids(const MeshBase &mesh, unsigned int sysnum)
A function for verifying that degree of freedom indexes are contiguous on each processor, as is required by libMesh numeric classes.
Definition: mesh_tools.C:1397
dof_id_type n_non_subactive_elem_of_type_at_level(const MeshBase &mesh, const ElemType type, const unsigned int level)
Definition: mesh_tools.C:744
void libmesh_assert_topology_consistent_procids< Elem >(const MeshBase &mesh)
Definition: mesh_tools.C:1439
uint8_t dof_id_type
Definition: id_types.h:67
unsigned int n_active_levels(const MeshBase &mesh)
Definition: mesh_tools.C:772
void libmesh_assert_valid_elem_ids(const MeshBase &mesh)
A function for verifying that ids and processor assignment of elements are correctly sorted (monotone...
Definition: mesh_tools.C:1318
void libmesh_assert_valid_refinement_flags(const MeshBase &mesh)
A function for verifying that refinement flags on elements are consistent between processors...
Definition: mesh_tools.C:2158
const RemoteElem * remote_elem
Definition: remote_elem.C:57