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