Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/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 :
59 : /**
60 : * SumElemWeight(Range) sums the number of nodes per element
61 : * for each element in the provided range. The join() method
62 : * defines how to combine the reduction operation from two
63 : * distinct instances of this class which may be executed on
64 : * separate threads.
65 : */
66 : class SumElemWeight
67 : {
68 : public:
69 0 : SumElemWeight () :
70 0 : _weight(0)
71 0 : {}
72 :
73 0 : SumElemWeight (SumElemWeight &, Threads::split) :
74 0 : _weight(0)
75 0 : {}
76 :
77 0 : void operator()(const ConstElemRange & range)
78 : {
79 0 : for (const auto & elem : range)
80 0 : _weight += elem->n_nodes();
81 0 : }
82 :
83 0 : dof_id_type weight() const
84 0 : { 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 0 : void join (const SumElemWeight & other)
90 0 : { _weight += other.weight(); }
91 : #endif
92 :
93 : private:
94 : dof_id_type _weight;
95 : };
96 :
97 :
98 : /**
99 : * FindBBox(Range) computes the bounding box for the objects
100 : * in the specified range. This class may be split and subranges
101 : * can be executed on separate threads. The join() method
102 : * defines how the results from two separate threads are combined.
103 : */
104 : class FindBBox
105 : {
106 : public:
107 3188843 : FindBBox () : _bbox()
108 28870 : {}
109 :
110 746 : FindBBox (FindBBox & other, Threads::split) :
111 746 : _bbox(other._bbox)
112 250 : {}
113 :
114 41718 : void operator()(const ConstNodeRange & range)
115 : {
116 5123960 : for (const auto & node : range)
117 : {
118 461782 : libmesh_assert(node);
119 5082242 : _bbox.union_with(*node);
120 : }
121 41718 : }
122 :
123 3197290 : void operator()(const ConstElemRange & range)
124 : {
125 26415734 : for (const auto & elem : range)
126 : {
127 792029 : libmesh_assert(elem);
128 23218444 : _bbox.union_with(elem->loose_bounding_box());
129 : }
130 3197290 : }
131 :
132 14731 : Point & min() { return _bbox.min(); }
133 :
134 14731 : 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 250 : void join (const FindBBox & other)
140 : {
141 746 : _bbox.union_with(other._bbox);
142 496 : }
143 : #endif
144 :
145 43009 : libMesh::BoundingBox & bbox ()
146 : {
147 43009 : return _bbox;
148 : }
149 :
150 : private:
151 : BoundingBox _bbox;
152 : };
153 :
154 : #ifdef DEBUG
155 3027608 : void assert_semiverify_dofobj(const Parallel::Communicator & communicator,
156 : const DofObject * d,
157 : unsigned int sysnum = libMesh::invalid_uint)
158 : {
159 3027608 : if (d)
160 : {
161 2981720 : const unsigned int n_sys = d->n_systems();
162 :
163 5963440 : std::vector<unsigned int> n_vars (n_sys, 0);
164 6660558 : for (unsigned int s = 0; s != n_sys; ++s)
165 3678838 : if (sysnum == s ||
166 : sysnum == libMesh::invalid_uint)
167 2981720 : n_vars[s] = d->n_vars(s);
168 :
169 : const unsigned int tot_n_vars =
170 2981720 : std::accumulate(n_vars.begin(), n_vars.end(), 0);
171 :
172 5963440 : std::vector<unsigned int> n_comp (tot_n_vars, 0);
173 5963440 : std::vector<dof_id_type> first_dof (tot_n_vars, 0);
174 :
175 6660558 : for (unsigned int s = 0, i=0; s != n_sys; ++s)
176 : {
177 3678838 : if (sysnum != s &&
178 : sysnum != libMesh::invalid_uint)
179 697118 : continue;
180 :
181 9922636 : for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
182 : {
183 6940916 : n_comp[i] = d->n_comp(s,v);
184 6940916 : first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : DofObject::invalid_id;
185 : }
186 : }
187 :
188 2981720 : libmesh_assert(communicator.semiverify(&n_sys));
189 2981720 : libmesh_assert(communicator.semiverify(&n_vars));
190 2981720 : libmesh_assert(communicator.semiverify(&n_comp));
191 2981720 : libmesh_assert(communicator.semiverify(&first_dof));
192 : }
193 : else // ranks w/o a DofObject d must still perform matching calls to the semiverify() collective
194 : {
195 45888 : const unsigned int * p_ui = nullptr;
196 45888 : const std::vector<unsigned int> * p_vui = nullptr;
197 45888 : const std::vector<dof_id_type> * p_vdid = nullptr;
198 :
199 45888 : libmesh_assert(communicator.semiverify(p_ui));
200 45888 : libmesh_assert(communicator.semiverify(p_vui));
201 45888 : libmesh_assert(communicator.semiverify(p_vui));
202 45888 : libmesh_assert(communicator.semiverify(p_vdid));
203 : }
204 3027608 : }
205 :
206 :
207 :
208 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
209 9471096 : 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 9471096 : if (d)
218 : {
219 9141624 : tempmin = tempmax = d->unique_id();
220 : }
221 : else
222 : {
223 329472 : TIMPI::Attributes<unique_id_type>::set_highest(tempmin);
224 329472 : TIMPI::Attributes<unique_id_type>::set_lowest(tempmax);
225 : }
226 9471096 : comm.min(tempmin);
227 9471096 : comm.max(tempmax);
228 18612720 : bool invalid = d && ((d->unique_id() != tempmin) ||
229 9141624 : (d->unique_id() != tempmax));
230 9471096 : comm.max(invalid);
231 :
232 : // First verify that everything is in sync
233 9471096 : libmesh_assert(!invalid);
234 :
235 : // Then verify that any remote id doesn't duplicate a local one.
236 9471096 : if (!d && tempmin == tempmax)
237 52832 : libmesh_assert(!unique_ids.count(tempmin));
238 9471096 : }
239 : #endif // LIBMESH_ENABLE_UNIQUE_ID
240 : #endif // DEBUG
241 :
242 227170 : 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 16232 : 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 1464478 : for (const auto & elem : node_to_elem_vec)
257 : {
258 : // We only care about active elements...
259 1237308 : if (elem->active())
260 : {
261 : // Which local node number is global_id?
262 1195954 : unsigned local_node_number = elem->local_node(global_id);
263 :
264 : // Make sure it was found
265 41354 : libmesh_assert_not_equal_to(local_node_number, libMesh::invalid_uint);
266 :
267 1237308 : 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 1237308 : if (!n_edges)
271 : {
272 1374 : switch (elem->type())
273 : {
274 592 : case EDGE2:
275 : {
276 18 : switch (local_node_number)
277 : {
278 296 : case 0:
279 : // The other node is a nodal neighbor
280 305 : neighbor_set.insert(elem->node_ptr(1));
281 296 : break;
282 :
283 296 : case 1:
284 : // The other node is a nodal neighbor
285 305 : neighbor_set.insert(elem->node_ptr(0));
286 296 : break;
287 :
288 0 : default:
289 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
290 : }
291 18 : break;
292 : }
293 :
294 498 : case EDGE3:
295 : {
296 18 : switch (local_node_number)
297 : {
298 : // The outside nodes have node 2 as a neighbor
299 356 : case 0:
300 : case 1:
301 370 : neighbor_set.insert(elem->node_ptr(2));
302 356 : break;
303 :
304 : // The middle node has the outer nodes as neighbors
305 142 : case 2:
306 146 : neighbor_set.insert(elem->node_ptr(0));
307 146 : neighbor_set.insert(elem->node_ptr(1));
308 142 : break;
309 :
310 0 : default:
311 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
312 : }
313 18 : break;
314 : }
315 :
316 284 : case EDGE4:
317 : {
318 8 : switch (local_node_number)
319 : {
320 71 : case 0:
321 : // The left-middle node is a nodal neighbor
322 73 : neighbor_set.insert(elem->node_ptr(2));
323 71 : break;
324 :
325 71 : case 1:
326 : // The right-middle node is a nodal neighbor
327 73 : neighbor_set.insert(elem->node_ptr(3));
328 71 : break;
329 :
330 : // The left-middle node
331 71 : case 2:
332 73 : neighbor_set.insert(elem->node_ptr(0));
333 73 : neighbor_set.insert(elem->node_ptr(3));
334 71 : break;
335 :
336 : // The right-middle node
337 71 : case 3:
338 73 : neighbor_set.insert(elem->node_ptr(1));
339 73 : neighbor_set.insert(elem->node_ptr(2));
340 71 : break;
341 :
342 0 : default:
343 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
344 : }
345 8 : break;
346 : }
347 :
348 0 : default:
349 0 : libmesh_error_msg("Unrecognized ElemType: " << Utility::enum_to_string(elem->type()) << std::endl);
350 : }
351 : }
352 :
353 1237308 : const auto elem_order = Elem::type_to_default_order_map[elem->type()];
354 :
355 : // Index of the current edge
356 41354 : unsigned current_edge = 0;
357 :
358 1237308 : const unsigned short n_nodes = elem->n_nodes();
359 :
360 5034507 : while (current_edge < n_edges)
361 : {
362 : // Find the edge the node is on
363 126386 : bool found_edge = false;
364 9178151 : for (; current_edge<n_edges; ++current_edge)
365 8414460 : if (elem->is_node_on_edge(local_node_number, current_edge))
366 : {
367 100178 : found_edge = true;
368 100178 : break;
369 : }
370 :
371 : // Did we find one?
372 3797199 : if (found_edge)
373 : {
374 3033508 : const Node * node_to_save = nullptr;
375 :
376 : // Find another node in this element on this edge
377 22556636 : for (unsigned other_node_this_edge = 0; other_node_this_edge != n_nodes; other_node_this_edge++)
378 : {
379 35732296 : const bool both_vertices = elem->is_vertex(local_node_number) &&
380 16209168 : elem->is_vertex(other_node_this_edge);
381 26046258 : if ( elem->is_node_on_edge(other_node_this_edge, current_edge) && // On the current edge
382 19571684 : elem->node_id(other_node_this_edge) != global_id && // But not the original node
383 : // vertex nodes on the same edge of higher order elements are not nodal neighbors
384 3638356 : (elem_order == 1 || !both_vertices))
385 : {
386 : // We've found a nodal neighbor! Save a pointer to it..
387 3294908 : node_to_save = elem->node_ptr(other_node_this_edge);
388 :
389 : // Make sure we found something
390 109454 : libmesh_assert(node_to_save != nullptr);
391 :
392 3185454 : neighbor_set.insert(node_to_save);
393 : }
394 : }
395 : }
396 :
397 : // Keep looking for edges, node may be on more than one edge
398 3797199 : current_edge++;
399 : }
400 : } // if (elem->active())
401 : } // for
402 :
403 : // Assign the entries from the set to the vector. Note: this
404 : // replaces any existing contents in neighbors and modifies its size
405 : // accordingly.
406 219054 : neighbors.assign(neighbor_set.begin(), neighbor_set.end());
407 227170 : }
408 :
409 : }
410 :
411 :
412 : namespace libMesh
413 : {
414 :
415 : // ------------------------------------------------------------
416 : // MeshTools functions
417 :
418 : namespace MeshTools
419 : {
420 :
421 0 : dof_id_type total_weight(const MeshBase & mesh)
422 : {
423 0 : if (!mesh.is_serial())
424 : {
425 0 : libmesh_parallel_only(mesh.comm());
426 0 : dof_id_type weight = MeshTools::weight (mesh, mesh.processor_id());
427 0 : mesh.comm().sum(weight);
428 : dof_id_type unpartitioned_weight =
429 0 : MeshTools::weight (mesh, DofObject::invalid_processor_id);
430 0 : return weight + unpartitioned_weight;
431 : }
432 :
433 0 : SumElemWeight sew;
434 :
435 0 : Threads::parallel_reduce (ConstElemRange (mesh.elements_begin(),
436 0 : mesh.elements_end()),
437 : sew);
438 0 : return sew.weight();
439 :
440 : }
441 :
442 :
443 :
444 0 : dof_id_type weight(const MeshBase & mesh, const processor_id_type pid)
445 : {
446 0 : SumElemWeight sew;
447 :
448 0 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
449 0 : mesh.pid_elements_end(pid)),
450 : sew);
451 0 : return sew.weight();
452 : }
453 :
454 :
455 :
456 12143 : void build_nodes_to_elem_map (const MeshBase & mesh,
457 : std::vector<std::vector<dof_id_type>> & nodes_to_elem_map)
458 : {
459 : // A vector indexed over all nodes is too inefficient to use for a
460 : // distributed mesh. Use the unordered_map API instead.
461 12143 : if (!mesh.is_serial())
462 : libmesh_deprecated();
463 :
464 12143 : nodes_to_elem_map.resize (mesh.max_node_id());
465 :
466 315619 : for (const auto & elem : mesh.element_ptr_range())
467 604940 : for (auto & node : elem->node_ref_range())
468 : {
469 12948 : libmesh_assert_less (node.id(), nodes_to_elem_map.size());
470 12948 : libmesh_assert_less (elem->id(), mesh.n_elem());
471 :
472 463416 : nodes_to_elem_map[node.id()].push_back(elem->id());
473 11449 : }
474 12143 : }
475 :
476 :
477 :
478 0 : void build_nodes_to_elem_map (const MeshBase & mesh,
479 : std::vector<std::vector<const Elem *>> & nodes_to_elem_map)
480 : {
481 : // A vector indexed over all nodes is too inefficient to use for a
482 : // distributed mesh. Use the unordered_map API instead.
483 0 : if (!mesh.is_serial())
484 : libmesh_deprecated();
485 :
486 0 : nodes_to_elem_map.resize (mesh.max_node_id());
487 :
488 0 : for (const auto & elem : mesh.element_ptr_range())
489 0 : for (auto & node : elem->node_ref_range())
490 : {
491 0 : libmesh_assert_less (node.id(), nodes_to_elem_map.size());
492 :
493 0 : nodes_to_elem_map[node.id()].push_back(elem);
494 0 : }
495 0 : }
496 :
497 :
498 :
499 369958 : void build_nodes_to_elem_map (const MeshBase & mesh,
500 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & nodes_to_elem_map)
501 : {
502 10794 : nodes_to_elem_map.clear();
503 :
504 62221890 : for (const auto & elem : mesh.element_ptr_range())
505 188562945 : for (auto & node : elem->node_ref_range())
506 154901065 : nodes_to_elem_map[node.id()].push_back(elem->id());
507 369958 : }
508 :
509 :
510 :
511 5965 : void build_nodes_to_elem_map (const MeshBase & mesh,
512 : std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
513 : {
514 168 : nodes_to_elem_map.clear();
515 :
516 776980 : for (const auto & elem : mesh.element_ptr_range())
517 2752339 : for (auto & node : elem->node_ref_range())
518 2341627 : nodes_to_elem_map[node.id()].push_back(elem);
519 5965 : }
520 :
521 :
522 :
523 : std::unordered_set<dof_id_type>
524 6261 : find_boundary_nodes(const MeshBase & mesh)
525 : {
526 182 : std::unordered_set<dof_id_type> boundary_nodes;
527 :
528 : // Loop over elements, find those on boundary, and
529 : // mark them as true in on_boundary.
530 1096532 : for (const auto & elem : mesh.active_element_ptr_range())
531 3126182 : for (auto s : elem->side_index_range())
532 2649090 : if (elem->neighbor_ptr(s) == nullptr) // on the boundary
533 : {
534 148240 : auto nodes_on_side = elem->nodes_on_side(s);
535 :
536 756684 : for (auto & local_id : nodes_on_side)
537 644632 : boundary_nodes.insert(elem->node_ptr(local_id)->id());
538 5897 : }
539 :
540 6261 : return boundary_nodes;
541 : }
542 :
543 : std::unordered_set<dof_id_type>
544 3137 : find_block_boundary_nodes(const MeshBase & mesh)
545 : {
546 94 : std::unordered_set<dof_id_type> block_boundary_nodes;
547 :
548 : // Loop over elements, find those on boundary, and
549 : // mark them as true in on_boundary.
550 864512 : for (const auto & elem : mesh.active_element_ptr_range())
551 2469648 : for (auto s : elem->side_index_range())
552 2090882 : if (elem->neighbor_ptr(s) && (elem->neighbor_ptr(s)->subdomain_id() != elem->subdomain_id()))
553 : {
554 0 : auto nodes_on_side = elem->nodes_on_side(s);
555 :
556 0 : for (auto & local_id : nodes_on_side)
557 0 : block_boundary_nodes.insert(elem->node_ptr(local_id)->id());
558 2949 : }
559 :
560 3137 : return block_boundary_nodes;
561 : }
562 :
563 :
564 :
565 : libMesh::BoundingBox
566 1598578 : create_bounding_box (const MeshBase & mesh)
567 : {
568 : // This function must be run on all processors at once
569 14139 : libmesh_parallel_only(mesh.comm());
570 :
571 14139 : FindBBox find_bbox;
572 :
573 : // Start with any unpartitioned elements we know about locally
574 3197156 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(DofObject::invalid_processor_id),
575 1612717 : mesh.pid_elements_end(DofObject::invalid_processor_id)),
576 : find_bbox);
577 :
578 : // And combine with our local elements
579 1598578 : find_bbox.bbox().union_with(create_local_bounding_box(mesh));
580 :
581 : // Compare the bounding boxes across processors
582 1598578 : mesh.comm().min(find_bbox.min());
583 1598578 : mesh.comm().max(find_bbox.max());
584 :
585 1598578 : return find_bbox.bbox();
586 : }
587 :
588 :
589 :
590 : libMesh::BoundingBox
591 20553 : create_nodal_bounding_box (const MeshBase & mesh)
592 : {
593 : // This function must be run on all processors at once
594 592 : libmesh_parallel_only(mesh.comm());
595 :
596 592 : FindBBox find_bbox;
597 :
598 : // Start with any unpartitioned nodes we know about locally
599 41106 : Threads::parallel_reduce (ConstNodeRange (mesh.pid_nodes_begin(DofObject::invalid_processor_id),
600 41106 : mesh.pid_nodes_end(DofObject::invalid_processor_id)),
601 : find_bbox);
602 :
603 : // Add our local nodes
604 41106 : Threads::parallel_reduce (ConstNodeRange (mesh.local_nodes_begin(),
605 40514 : mesh.local_nodes_end()),
606 : find_bbox);
607 :
608 : // Compare the bounding boxes across processors
609 20553 : mesh.comm().min(find_bbox.min());
610 20553 : mesh.comm().max(find_bbox.max());
611 :
612 20553 : return find_bbox.bbox();
613 : }
614 :
615 :
616 :
617 : Sphere
618 0 : bounding_sphere(const MeshBase & mesh)
619 : {
620 0 : libMesh::BoundingBox bbox = create_bounding_box(mesh);
621 :
622 0 : const Real diag = (bbox.second - bbox.first).norm();
623 0 : const Point cent = (bbox.second + bbox.first)/2;
624 :
625 0 : return Sphere (cent, .5*diag);
626 : }
627 :
628 :
629 :
630 : libMesh::BoundingBox
631 1598578 : create_local_bounding_box (const MeshBase & mesh)
632 : {
633 14139 : FindBBox find_bbox;
634 :
635 1598578 : Threads::parallel_reduce (mesh.active_local_element_stored_range(),
636 : find_bbox);
637 :
638 1598578 : return find_bbox.bbox();
639 : }
640 :
641 :
642 :
643 : libMesh::BoundingBox
644 0 : create_processor_bounding_box (const MeshBase & mesh,
645 : const processor_id_type pid)
646 : {
647 : // This can only be run in parallel, with consistent arguments.
648 0 : libmesh_parallel_only(mesh.comm());
649 0 : libmesh_assert(mesh.comm().verify(pid));
650 :
651 0 : libmesh_assert_less (pid, mesh.n_processors());
652 :
653 0 : FindBBox find_bbox;
654 :
655 0 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
656 0 : mesh.pid_elements_end(pid)),
657 : find_bbox);
658 :
659 : // Compare the bounding boxes across processors
660 0 : mesh.comm().min(find_bbox.min());
661 0 : mesh.comm().max(find_bbox.max());
662 :
663 0 : return find_bbox.bbox();
664 : }
665 :
666 :
667 :
668 : Sphere
669 0 : processor_bounding_sphere (const MeshBase & mesh,
670 : const processor_id_type pid)
671 : {
672 : libMesh::BoundingBox bbox =
673 0 : create_processor_bounding_box(mesh, pid);
674 :
675 0 : const Real diag = (bbox.second - bbox.first).norm();
676 0 : const Point cent = (bbox.second + bbox.first)/2;
677 :
678 0 : return Sphere (cent, .5*diag);
679 : }
680 :
681 :
682 :
683 : libMesh::BoundingBox
684 0 : create_subdomain_bounding_box (const MeshBase & mesh,
685 : const subdomain_id_type sid)
686 : {
687 : // This can only be run in parallel, with consistent arguments.
688 0 : libmesh_parallel_only(mesh.comm());
689 0 : libmesh_assert(mesh.comm().verify(sid));
690 :
691 0 : FindBBox find_bbox;
692 :
693 : Threads::parallel_reduce
694 0 : (ConstElemRange (mesh.active_local_subdomain_elements_begin(sid),
695 0 : mesh.active_local_subdomain_elements_end(sid)),
696 : find_bbox);
697 :
698 : // Compare the bounding boxes across processors
699 0 : mesh.comm().min(find_bbox.min());
700 0 : mesh.comm().max(find_bbox.max());
701 :
702 0 : return find_bbox.bbox();
703 : }
704 :
705 :
706 :
707 : Sphere
708 0 : subdomain_bounding_sphere (const MeshBase & mesh,
709 : const subdomain_id_type sid)
710 : {
711 : libMesh::BoundingBox bbox =
712 0 : create_subdomain_bounding_box(mesh, sid);
713 :
714 0 : const Real diag = (bbox.second - bbox.first).norm();
715 0 : const Point cent = (bbox.second + bbox.first)/2;
716 :
717 0 : return Sphere (cent, .5*diag);
718 : }
719 :
720 :
721 :
722 0 : void elem_types (const MeshBase & mesh,
723 : std::vector<ElemType> & et)
724 : {
725 : // Loop over the the elements. If the current element type isn't in
726 : // the vector, insert it.
727 0 : for (const auto & elem : mesh.element_ptr_range())
728 0 : if (!std::count(et.begin(), et.end(), elem->type()))
729 0 : et.push_back(elem->type());
730 0 : }
731 :
732 :
733 :
734 0 : dof_id_type n_elem_of_type (const MeshBase & mesh,
735 : const ElemType type)
736 : {
737 0 : return static_cast<dof_id_type>(std::distance(mesh.type_elements_begin(type),
738 0 : mesh.type_elements_end (type)));
739 : }
740 :
741 :
742 :
743 0 : dof_id_type n_active_elem_of_type (const MeshBase & mesh,
744 : const ElemType type)
745 : {
746 0 : return static_cast<dof_id_type>(std::distance(mesh.active_type_elements_begin(type),
747 0 : mesh.active_type_elements_end (type)));
748 : }
749 :
750 0 : dof_id_type n_non_subactive_elem_of_type_at_level(const MeshBase & mesh,
751 : const ElemType type,
752 : const unsigned int level)
753 : {
754 0 : dof_id_type cnt = 0;
755 :
756 : // iterate over the elements of the specified type
757 0 : for (const auto & elem : as_range(mesh.type_elements_begin(type),
758 0 : mesh.type_elements_end(type)))
759 0 : if ((elem->level() == level) && !elem->subactive())
760 0 : cnt++;
761 :
762 0 : return cnt;
763 : }
764 :
765 :
766 5471 : unsigned int n_active_local_levels(const MeshBase & mesh)
767 : {
768 : struct LevelCounter {
769 : unsigned int nl;
770 :
771 5471 : LevelCounter () : nl(0) {}
772 :
773 0 : LevelCounter (LevelCounter &, Threads::split) :
774 0 : nl(0) {}
775 :
776 5471 : void operator()(const ConstElemRange & range) {
777 70642 : for (const Elem * elem : range)
778 69039 : nl = std::max(elem->level() + 1, nl);
779 5471 : }
780 :
781 0 : void join(const LevelCounter & other) {
782 0 : nl = std::max(nl, other.nl);
783 0 : }
784 : };
785 :
786 156 : LevelCounter counter;
787 :
788 5471 : Threads::parallel_reduce(mesh.active_local_element_stored_range(), counter);
789 :
790 5471 : return counter.nl;
791 : }
792 :
793 :
794 :
795 5471 : unsigned int n_active_levels(const MeshBase & mesh)
796 : {
797 156 : libmesh_parallel_only(mesh.comm());
798 :
799 5471 : unsigned int nl = n_active_local_levels(mesh);
800 :
801 10786 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
802 21572 : mesh.unpartitioned_elements_end()))
803 0 : if (elem->active())
804 5159 : nl = std::max(elem->level() + 1, nl);
805 :
806 5471 : mesh.comm().max(nl);
807 5471 : return nl;
808 : }
809 :
810 :
811 :
812 2441849 : unsigned int n_local_levels(const MeshBase & mesh)
813 : {
814 2441849 : unsigned int nl = 0;
815 :
816 6545837 : for (const auto & elem : as_range(mesh.local_elements_begin(),
817 43415282 : mesh.local_elements_end()))
818 38240297 : nl = std::max(elem->level() + 1, nl);
819 :
820 2441849 : return nl;
821 : }
822 :
823 :
824 :
825 2441849 : unsigned int n_levels(const MeshBase & mesh)
826 : {
827 40673 : libmesh_parallel_only(mesh.comm());
828 :
829 2441849 : unsigned int nl = n_local_levels(mesh);
830 :
831 5255865 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
832 24802460 : mesh.unpartitioned_elements_end()))
833 17476937 : nl = std::max(elem->level() + 1, nl);
834 :
835 2441849 : mesh.comm().max(nl);
836 :
837 : // n_levels() is only valid and should only be called in cases where
838 : // the mesh is validly distributed (or serialized). Let's run an
839 : // expensive test in debug mode to make sure this is such a case.
840 : #ifdef DEBUG
841 40673 : const unsigned int paranoid_nl = paranoid_n_levels(mesh);
842 40673 : libmesh_assert_equal_to(nl, paranoid_nl);
843 : #endif
844 2441849 : return nl;
845 : }
846 :
847 :
848 :
849 55692 : unsigned int paranoid_n_levels(const MeshBase & mesh)
850 : {
851 41117 : libmesh_parallel_only(mesh.comm());
852 :
853 55692 : unsigned int nl = 0;
854 4316161 : for (const auto & elem : mesh.element_ptr_range())
855 4260025 : nl = std::max(elem->level() + 1, nl);
856 :
857 55692 : mesh.comm().max(nl);
858 55692 : return nl;
859 : }
860 :
861 :
862 :
863 639 : dof_id_type n_connected_components(const MeshBase & mesh,
864 : Real constraint_tol)
865 : {
866 54 : LOG_SCOPE("n_connected_components()", "MeshTools");
867 :
868 : // Yes, I'm being lazy. This is for mesh analysis before a
869 : // simulation, not anything going in any loops.
870 639 : if (!mesh.is_serial_on_zero())
871 0 : libmesh_not_implemented();
872 :
873 639 : dof_id_type n_components = 0;
874 :
875 657 : if (mesh.processor_id())
876 : {
877 522 : mesh.comm().broadcast(n_components);
878 531 : return n_components;
879 : }
880 :
881 : // All nodes in a set here are connected (at least indirectly) to
882 : // all other nodes in the same set, but have not yet been discovered
883 : // to be connected to nodes in other sets.
884 : //
885 : // Using an unordered_set of ids rather than a set of pointers seems
886 : // to be roughly 150% faster?
887 : // typedef const Node * node_entry_type
888 : typedef dof_id_type node_entry_type;
889 9 : std::vector<std::unordered_set<node_entry_type>> components;
890 :
891 : // With a typical mesh with few components and somewhat-contiguous
892 : // ordering, vector performance should be fine. With a mesh with
893 : // many components or completely scrambled ordering, performance
894 : // can be a disaster.
895 2742 : auto find_component = [&components](node_entry_type n) {
896 6158 : for (auto & c: components)
897 4064 : if (c.find(n) != c.end())
898 54 : return &c;
899 :
900 162 : return (std::unordered_set<node_entry_type> *)(nullptr);
901 108 : };
902 :
903 : auto add_to_component =
904 2010 : [&find_component]
905 2202 : (std::unordered_set<node_entry_type> & component, node_entry_type n) {
906 : // We may already be in the desired component
907 2412 : if (component.find(n) != component.end())
908 51 : return;
909 :
910 1950 : auto current_component = find_component(n);
911 :
912 : // Didn't we *just* check this?
913 150 : libmesh_assert (&component != current_component);
914 :
915 : // If we're unknown, we should be in the desired component
916 1950 : if (!current_component)
917 147 : component.insert(n);
918 :
919 : // If we think we're in another component, it should actually be
920 : // part of the desired component
921 : else
922 : {
923 : // Merge the component likely to be smaller into the one
924 : // likely to be larger - this is orders of magnitude faster
925 : // than the other way around!
926 3 : current_component->merge(component);
927 3 : current_component->swap(component);
928 3 : libmesh_assert(current_component->empty());
929 : }
930 108 : };
931 :
932 9 : const auto & constraint_rows = mesh.get_constraint_rows();
933 :
934 1659 : for (const auto & elem : mesh.element_ptr_range())
935 : {
936 : // const node_entry_type first_node = elem->node_ptr(0);
937 792 : const node_entry_type first_node = elem->node_id(0);
938 :
939 792 : auto component = find_component(first_node);
940 :
941 : // If we didn't find one, make a new one, reusing an existing
942 : // slot if possible or growing our vector if necessary
943 792 : if (!component)
944 684 : for (auto & c: components)
945 354 : if (c.empty())
946 0 : component = &c;
947 :
948 432 : if (!component)
949 219 : component = &components.emplace_back();
950 :
951 3234 : for (const Node & node : elem->node_ref_range())
952 : {
953 : // const node_entry_type n = &node;
954 396 : const node_entry_type n = node.id();
955 2376 : add_to_component(*component, n);
956 :
957 2376 : auto it = constraint_rows.find(&node);
958 2376 : if (it == constraint_rows.end())
959 195 : continue;
960 :
961 72 : for (const auto & [pr, val] : it->second)
962 : {
963 : // Ignore too-trivial constraint coefficients if
964 : // we get a non-default-0 constraint_tol
965 39 : if (std::abs(val) < constraint_tol)
966 0 : continue;
967 :
968 36 : const Elem * spline_elem = pr.first;
969 3 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
970 :
971 : const Node * spline_node =
972 36 : spline_elem->node_ptr(pr.second);
973 :
974 : // add_to_component(*component, spline_node);
975 36 : add_to_component(*component, spline_node->id());
976 : }
977 : }
978 90 : }
979 :
980 327 : for (auto & component : components)
981 219 : if (!component.empty())
982 144 : ++n_components;
983 :
984 : // We calculated this on proc 0; now let everyone else know too
985 108 : mesh.comm().broadcast(n_components);
986 :
987 108 : return n_components;
988 90 : }
989 :
990 :
991 :
992 0 : void get_not_subactive_node_ids(const MeshBase & mesh,
993 : std::set<dof_id_type> & not_subactive_node_ids)
994 : {
995 0 : for (const auto & elem : mesh.element_ptr_range())
996 0 : if (!elem->subactive())
997 0 : for (auto & n : elem->node_ref_range())
998 0 : not_subactive_node_ids.insert(n.id());
999 0 : }
1000 :
1001 :
1002 :
1003 686411 : dof_id_type n_elem (const MeshBase::const_element_iterator & begin,
1004 : const MeshBase::const_element_iterator & end)
1005 : {
1006 1346584 : return cast_int<dof_id_type>(std::distance(begin, end));
1007 : }
1008 :
1009 :
1010 :
1011 0 : dof_id_type n_nodes (const MeshBase::const_node_iterator & begin,
1012 : const MeshBase::const_node_iterator & end)
1013 : {
1014 0 : return cast_int<dof_id_type>(std::distance(begin, end));
1015 : }
1016 :
1017 :
1018 :
1019 5268 : Real volume (const MeshBase & mesh,
1020 : unsigned int dim)
1021 : {
1022 148 : libmesh_parallel_only(mesh.comm());
1023 :
1024 5268 : if (dim == libMesh::invalid_uint)
1025 5268 : dim = mesh.mesh_dimension();
1026 :
1027 5268 : Real vol = 0;
1028 :
1029 : // first my local elements
1030 12243 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1031 63639 : mesh.local_elements_end()))
1032 22359 : if (elem->dim() == dim)
1033 27331 : vol += elem->volume();
1034 :
1035 : // then count any unpartitioned objects, once
1036 5416 : if (mesh.processor_id() == 0)
1037 1730 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1038 3460 : mesh.unpartitioned_elements_end()))
1039 0 : if (elem->dim() == dim)
1040 754 : vol += elem->volume();
1041 :
1042 5268 : mesh.comm().sum(vol);
1043 5268 : return vol;
1044 : }
1045 :
1046 :
1047 :
1048 5471 : unsigned int n_p_levels (const MeshBase & mesh)
1049 : {
1050 156 : libmesh_parallel_only(mesh.comm());
1051 :
1052 5471 : unsigned int max_p_level = 0;
1053 :
1054 : // first my local elements
1055 17480 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1056 93386 : mesh.local_elements_end()))
1057 83667 : max_p_level = std::max(elem->p_level(), max_p_level);
1058 :
1059 : // then any unpartitioned objects
1060 10786 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1061 21572 : mesh.unpartitioned_elements_end()))
1062 5159 : max_p_level = std::max(elem->p_level(), max_p_level);
1063 :
1064 5471 : mesh.comm().max(max_p_level);
1065 5471 : return max_p_level + 1;
1066 : }
1067 :
1068 :
1069 :
1070 0 : void find_nodal_neighbors(const MeshBase &,
1071 : const Node & node,
1072 : const std::vector<std::vector<const Elem *>> & nodes_to_elem_map,
1073 : std::vector<const Node *> & neighbors)
1074 : {
1075 0 : find_nodal_neighbors_helper(node.id(), nodes_to_elem_map[node.id()],
1076 : neighbors);
1077 0 : }
1078 :
1079 :
1080 :
1081 227170 : void find_nodal_neighbors(const MeshBase &,
1082 : const Node & node,
1083 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1084 : std::vector<const Node *> & neighbors)
1085 : {
1086 : const std::vector<const Elem *> node_to_elem_vec =
1087 235286 : libmesh_map_find(nodes_to_elem_map, node.id());
1088 227170 : find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
1089 227170 : }
1090 :
1091 50735 : void find_nodal_or_face_neighbors(
1092 : const MeshBase & mesh,
1093 : const Node & node,
1094 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1095 : std::vector<const Node *> & neighbors)
1096 : {
1097 : // Find all the nodal neighbors... that is the nodes directly connected
1098 : // to this node through one edge.
1099 50735 : find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
1100 :
1101 : // If no neighbors are found, use all nodes on the containing side as
1102 : // neighbors.
1103 50735 : if (!neighbors.size())
1104 : {
1105 : // Grab the element containing node
1106 6804 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
1107 : // Find the element side containing node
1108 17503 : for (const auto &side : elem->side_index_range())
1109 : {
1110 17503 : const auto &nodes_on_side = elem->nodes_on_side(side);
1111 : const auto it =
1112 23074 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
1113 64593 : return elem->node_id(local_node_id) == node.id();
1114 1946 : });
1115 :
1116 18476 : if (it != nodes_on_side.end())
1117 : {
1118 58008 : for (const auto &local_node_id : nodes_on_side)
1119 : // No need to add node itself as a neighbor
1120 53931 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
1121 2727 : *node_ptr != node)
1122 44400 : neighbors.push_back(node_ptr);
1123 347 : break;
1124 : }
1125 : }
1126 : }
1127 3146 : libmesh_assert(neighbors.size());
1128 50735 : }
1129 :
1130 :
1131 :
1132 0 : void find_hanging_nodes_and_parents(const MeshBase & mesh,
1133 : std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes)
1134 : {
1135 : // Loop through all the elements
1136 0 : for (auto & elem : mesh.active_local_element_ptr_range())
1137 0 : if (elem->type() == QUAD4)
1138 0 : for (auto s : elem->side_index_range())
1139 : {
1140 : // Loop over the sides looking for sides that have hanging nodes
1141 : // This code is inspired by compute_proj_constraints()
1142 0 : const Elem * neigh = elem->neighbor_ptr(s);
1143 :
1144 : // If not a boundary side
1145 0 : if (neigh != nullptr)
1146 : {
1147 : // Is there a coarser element next to this one?
1148 0 : if (neigh->level() < elem->level())
1149 : {
1150 0 : const Elem * ancestor = elem;
1151 0 : while (neigh->level() < ancestor->level())
1152 0 : ancestor = ancestor->parent();
1153 0 : unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
1154 0 : libmesh_assert_less (s_neigh, neigh->n_neighbors());
1155 :
1156 : // Couple of helper uints...
1157 0 : unsigned int local_node1=0;
1158 0 : unsigned int local_node2=0;
1159 :
1160 0 : bool found_in_neighbor = false;
1161 :
1162 : // Find the two vertices that make up this side
1163 0 : while (!elem->is_node_on_side(local_node1++,s)) { }
1164 0 : local_node1--;
1165 :
1166 : // Start looking for the second one with the next node
1167 0 : local_node2=local_node1+1;
1168 :
1169 : // Find the other one
1170 0 : while (!elem->is_node_on_side(local_node2++,s)) { }
1171 0 : local_node2--;
1172 :
1173 : //Pull out their global ids:
1174 0 : dof_id_type node1 = elem->node_id(local_node1);
1175 0 : dof_id_type node2 = elem->node_id(local_node2);
1176 :
1177 : // Now find which node is present in the neighbor
1178 : // FIXME This assumes a level one rule!
1179 : // The _other_ one is the hanging node
1180 :
1181 : // First look for the first one
1182 : // FIXME could be streamlined a bit
1183 0 : for (unsigned int n=0;n<neigh->n_sides();n++)
1184 0 : if (neigh->node_id(n) == node1)
1185 0 : found_in_neighbor=true;
1186 :
1187 0 : dof_id_type hanging_node=0;
1188 :
1189 0 : if (!found_in_neighbor)
1190 0 : hanging_node=node1;
1191 : else // If it wasn't node1 then it must be node2!
1192 0 : hanging_node=node2;
1193 :
1194 : // Reset for reuse
1195 0 : local_node1=0;
1196 :
1197 : // Find the first node that makes up the side in the neighbor (these should be the parent nodes)
1198 0 : while (!neigh->is_node_on_side(local_node1++,s_neigh)) { }
1199 0 : local_node1--;
1200 :
1201 0 : local_node2=local_node1+1;
1202 :
1203 : // Find the second node...
1204 0 : while (!neigh->is_node_on_side(local_node2++,s_neigh)) { }
1205 0 : local_node2--;
1206 :
1207 : // Save them if we haven't already found the parents for this one
1208 0 : if (hanging_nodes[hanging_node].size()<2)
1209 : {
1210 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node1));
1211 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node2));
1212 : }
1213 : }
1214 : }
1215 0 : }
1216 0 : }
1217 :
1218 :
1219 :
1220 36 : void clear_spline_nodes(MeshBase & mesh)
1221 : {
1222 6 : std::vector<Elem *> nodeelem_to_delete;
1223 :
1224 8737 : for (auto & elem : mesh.element_ptr_range())
1225 5073 : if (elem->type() == NODEELEM &&
1226 4140 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
1227 4170 : nodeelem_to_delete.push_back(elem);
1228 :
1229 3 : auto & constraint_rows = mesh.set_constraint_rows();
1230 :
1231 : // All our constraint_rows ought to be for spline constraints we're
1232 : // about to get rid of.
1233 : #ifndef NDEBUG
1234 762 : for (auto & node_row : constraint_rows)
1235 2064 : for (auto pr : node_row.second)
1236 : {
1237 1305 : const Elem * elem = pr.first.first;
1238 1305 : libmesh_assert(elem->type() == NODEELEM);
1239 1305 : libmesh_assert(elem->mapping_type() == RATIONAL_BERNSTEIN_MAP);
1240 : }
1241 : #endif
1242 :
1243 3 : constraint_rows.clear();
1244 :
1245 4176 : for (Elem * elem : nodeelem_to_delete)
1246 : {
1247 690 : Node * node = elem->node_ptr(0);
1248 4140 : mesh.delete_elem(elem);
1249 4140 : mesh.delete_node(node);
1250 : }
1251 36 : }
1252 :
1253 :
1254 :
1255 1798 : bool valid_is_prepared (const MeshBase & mesh)
1256 : {
1257 202 : LOG_SCOPE("valid_is_prepared()", "MeshTools");
1258 :
1259 1798 : if (!mesh.is_prepared())
1260 2 : return true;
1261 :
1262 1870 : std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
1263 :
1264 : // Try preparing (without allowing repartitioning or renumbering, to
1265 : // avoid false assertion failures)
1266 124 : bool old_allow_renumbering = mesh_clone->allow_renumbering();
1267 74 : mesh_clone->allow_renumbering(false);
1268 :
1269 124 : bool old_skip_partitioning = mesh_clone->skip_partitioning();
1270 74 : mesh_clone->skip_partitioning(true);
1271 :
1272 124 : bool old_allow_remote_element_removal = mesh_clone->allow_remote_element_removal();
1273 74 : mesh_clone->allow_remote_element_removal(false);
1274 :
1275 : // Call prepare_for_use() on the clone
1276 1796 : mesh_clone->prepare_for_use();
1277 :
1278 : // Restore original flag values
1279 74 : mesh_clone->allow_renumbering(old_allow_renumbering);
1280 74 : mesh_clone->skip_partitioning(old_skip_partitioning);
1281 74 : mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal);
1282 :
1283 : // Check whether the original and clone compare equal
1284 1796 : return (mesh == *mesh_clone);
1285 1672 : }
1286 :
1287 :
1288 :
1289 : #ifndef NDEBUG
1290 :
1291 136 : void libmesh_assert_equal_n_systems (const MeshBase & mesh)
1292 : {
1293 272 : LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
1294 :
1295 136 : unsigned int n_sys = libMesh::invalid_uint;
1296 :
1297 12632 : for (const auto & elem : mesh.element_ptr_range())
1298 : {
1299 12496 : if (n_sys == libMesh::invalid_uint)
1300 136 : n_sys = elem->n_systems();
1301 : else
1302 12360 : libmesh_assert_equal_to (elem->n_systems(), n_sys);
1303 : }
1304 :
1305 31283 : for (const auto & node : mesh.node_ptr_range())
1306 : {
1307 31147 : if (n_sys == libMesh::invalid_uint)
1308 0 : n_sys = node->n_systems();
1309 : else
1310 31147 : libmesh_assert_equal_to (node->n_systems(), n_sys);
1311 : }
1312 136 : }
1313 :
1314 :
1315 :
1316 : #ifdef LIBMESH_ENABLE_AMR
1317 0 : void libmesh_assert_old_dof_objects (const MeshBase & mesh)
1318 : {
1319 0 : LOG_SCOPE("libmesh_assert_old_dof_objects()", "MeshTools");
1320 :
1321 0 : for (const auto & elem : mesh.element_ptr_range())
1322 : {
1323 0 : if (elem->refinement_flag() == Elem::JUST_REFINED ||
1324 0 : elem->refinement_flag() == Elem::INACTIVE)
1325 0 : continue;
1326 :
1327 0 : if (elem->has_dofs())
1328 0 : libmesh_assert(elem->get_old_dof_object());
1329 :
1330 0 : for (auto & node : elem->node_ref_range())
1331 0 : if (node.has_dofs())
1332 0 : libmesh_assert(node.get_old_dof_object());
1333 : }
1334 0 : }
1335 : #else
1336 : void libmesh_assert_old_dof_objects (const MeshBase &) {}
1337 : #endif // LIBMESH_ENABLE_AMR
1338 :
1339 :
1340 :
1341 0 : void libmesh_assert_valid_node_pointers(const MeshBase & mesh)
1342 : {
1343 0 : LOG_SCOPE("libmesh_assert_valid_node_pointers()", "MeshTools");
1344 :
1345 : // Here we specifically do not want "auto &" because we need to
1346 : // reseat the (temporary) pointer variable in the loop below,
1347 : // without modifying the original.
1348 0 : for (const Elem * elem : mesh.element_ptr_range())
1349 : {
1350 0 : libmesh_assert (elem);
1351 0 : while (elem)
1352 : {
1353 0 : elem->libmesh_assert_valid_node_pointers();
1354 0 : for (auto n : elem->neighbor_ptr_range())
1355 0 : if (n && n != remote_elem)
1356 0 : n->libmesh_assert_valid_node_pointers();
1357 :
1358 0 : libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1359 0 : elem = elem->parent();
1360 : }
1361 : }
1362 0 : }
1363 :
1364 :
1365 :
1366 11787 : void libmesh_assert_valid_remote_elems(const MeshBase & mesh)
1367 : {
1368 23574 : LOG_SCOPE("libmesh_assert_valid_remote_elems()", "MeshTools");
1369 :
1370 689781 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1371 1379562 : mesh.local_elements_end()))
1372 : {
1373 677994 : libmesh_assert (elem);
1374 :
1375 : // We currently don't allow active_local_elements to have
1376 : // remote_elem neighbors
1377 677994 : if (elem->active())
1378 2953671 : for (auto n : elem->neighbor_ptr_range())
1379 2376655 : libmesh_assert_not_equal_to (n, remote_elem);
1380 :
1381 : #ifdef LIBMESH_ENABLE_AMR
1382 677994 : const Elem * parent = elem->parent();
1383 677994 : if (parent)
1384 377198 : libmesh_assert_not_equal_to (parent, remote_elem);
1385 :
1386 : // We can only be strict about active elements' subactive
1387 : // children
1388 677994 : if (elem->active() && elem->has_children())
1389 30186 : for (auto & child : elem->child_ref_range())
1390 24168 : libmesh_assert_not_equal_to (&child, remote_elem);
1391 : #endif
1392 : }
1393 11787 : }
1394 :
1395 :
1396 :
1397 424 : void libmesh_assert_valid_elem_ids(const MeshBase & mesh)
1398 : {
1399 848 : LOG_SCOPE("libmesh_assert_valid_elem_ids()", "MeshTools");
1400 :
1401 424 : processor_id_type lastprocid = 0;
1402 424 : dof_id_type lastelemid = 0;
1403 :
1404 13881 : for (const auto & elem : mesh.active_element_ptr_range())
1405 : {
1406 13457 : libmesh_assert (elem);
1407 13457 : processor_id_type elemprocid = elem->processor_id();
1408 13457 : dof_id_type elemid = elem->id();
1409 :
1410 13457 : libmesh_assert_greater_equal (elemid, lastelemid);
1411 13457 : libmesh_assert_greater_equal (elemprocid, lastprocid);
1412 :
1413 13457 : lastelemid = elemid;
1414 13457 : lastprocid = elemprocid;
1415 : }
1416 424 : }
1417 :
1418 :
1419 :
1420 866 : void libmesh_assert_valid_amr_elem_ids(const MeshBase & mesh)
1421 : {
1422 1732 : LOG_SCOPE("libmesh_assert_valid_amr_elem_ids()", "MeshTools");
1423 :
1424 307434 : for (const auto & elem : mesh.element_ptr_range())
1425 : {
1426 306568 : libmesh_assert (elem);
1427 :
1428 306568 : const Elem * parent = elem->parent();
1429 :
1430 306568 : if (parent)
1431 : {
1432 24404 : libmesh_assert_greater_equal (elem->id(), parent->id());
1433 24404 : libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1434 : }
1435 : }
1436 866 : }
1437 :
1438 :
1439 :
1440 17830 : void libmesh_assert_valid_amr_interior_parents(const MeshBase & mesh)
1441 : {
1442 35660 : LOG_SCOPE("libmesh_assert_valid_amr_interior_parents()", "MeshTools");
1443 :
1444 1486673 : for (const auto & elem : mesh.element_ptr_range())
1445 : {
1446 1468843 : libmesh_assert (elem);
1447 :
1448 : // We can skip to the next element if we're full-dimension
1449 : // and therefore don't have any interior parents
1450 1468843 : if (elem->dim() >= LIBMESH_DIM)
1451 657582 : continue;
1452 :
1453 811261 : const Elem * ip = elem->interior_parent();
1454 :
1455 811261 : const Elem * parent = elem->parent();
1456 :
1457 811261 : if (ip && (ip != remote_elem) && parent)
1458 : {
1459 1000 : libmesh_assert_equal_to (ip->top_parent(),
1460 : elem->top_parent()->interior_parent());
1461 :
1462 1000 : if (ip->level() == elem->level())
1463 1000 : libmesh_assert_equal_to (ip->parent(),
1464 : parent->interior_parent());
1465 : else
1466 : {
1467 0 : libmesh_assert_less (ip->level(), elem->level());
1468 0 : libmesh_assert_equal_to (ip, parent->interior_parent());
1469 : }
1470 : }
1471 : }
1472 17830 : }
1473 :
1474 :
1475 :
1476 0 : void libmesh_assert_contiguous_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1477 : {
1478 0 : LOG_SCOPE("libmesh_assert_contiguous_dof_ids()", "MeshTools");
1479 :
1480 0 : if (mesh.n_processors() == 1)
1481 0 : return;
1482 :
1483 0 : libmesh_parallel_only(mesh.comm());
1484 :
1485 0 : dof_id_type min_dof_id = std::numeric_limits<dof_id_type>::max(),
1486 0 : max_dof_id = std::numeric_limits<dof_id_type>::min();
1487 :
1488 : // Figure out what our local dof id range is
1489 0 : for (const auto * node : mesh.local_node_ptr_range())
1490 : {
1491 0 : for (auto v : make_range(node->n_vars(sysnum)))
1492 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1493 : {
1494 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1495 0 : min_dof_id = std::min (min_dof_id, id);
1496 0 : max_dof_id = std::max (max_dof_id, id);
1497 : }
1498 : }
1499 :
1500 : // Make sure no other processors' ids are inside it
1501 0 : for (const auto * node : mesh.node_ptr_range())
1502 : {
1503 0 : if (node->processor_id() == mesh.processor_id())
1504 0 : continue;
1505 0 : for (auto v : make_range(node->n_vars(sysnum)))
1506 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1507 : {
1508 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1509 0 : libmesh_assert (id < min_dof_id ||
1510 : id > max_dof_id);
1511 : }
1512 : }
1513 : }
1514 :
1515 :
1516 :
1517 : template <>
1518 24018 : void libmesh_assert_topology_consistent_procids<Elem>(const MeshBase & mesh)
1519 : {
1520 48036 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1521 :
1522 : // This parameter is not used when !LIBMESH_ENABLE_AMR
1523 24018 : libmesh_ignore(mesh);
1524 :
1525 : // If we're adaptively refining, check processor ids for consistency
1526 : // between parents and children.
1527 : #ifdef LIBMESH_ENABLE_AMR
1528 :
1529 : // Ancestor elements we won't worry about, but subactive and active
1530 : // elements ought to have parents with consistent processor ids
1531 2817350 : for (const auto & elem : mesh.element_ptr_range())
1532 : {
1533 2793332 : libmesh_assert(elem);
1534 :
1535 2793332 : if (!elem->active() && !elem->subactive())
1536 306696 : continue;
1537 :
1538 2486636 : const Elem * parent = elem->parent();
1539 :
1540 2486636 : if (parent)
1541 : {
1542 1225076 : libmesh_assert(parent->has_children());
1543 1225076 : processor_id_type parent_procid = parent->processor_id();
1544 1225076 : bool matching_child_id = false;
1545 : // If we've got a remote_elem then we don't know whether
1546 : // it's responsible for the parent's processor id; all
1547 : // we can do is assume it is and let its processor fail
1548 : // an assert if there's something wrong.
1549 7551420 : for (auto & child : parent->child_ref_range())
1550 12649928 : if (&child == remote_elem ||
1551 6323584 : child.processor_id() == parent_procid)
1552 6154620 : matching_child_id = true;
1553 1225076 : libmesh_assert(matching_child_id);
1554 : }
1555 : }
1556 : #endif
1557 24018 : }
1558 :
1559 :
1560 :
1561 : template <>
1562 36821 : void libmesh_assert_topology_consistent_procids<Node>(const MeshBase & mesh)
1563 : {
1564 36821 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1565 :
1566 36821 : if (mesh.n_processors() == 1)
1567 1 : return;
1568 :
1569 36820 : libmesh_parallel_only(mesh.comm());
1570 :
1571 : // We want this test to be valid even when called after nodes have
1572 : // been added asynchronously but before they're renumbered.
1573 : //
1574 : // Plus, some code (looking at you, stitch_meshes) modifies
1575 : // DofObject ids without keeping max_elem_id()/max_node_id()
1576 : // consistent, but that's done in a safe way for performance
1577 : // reasons, so we'll play along and just figure out new max ids
1578 : // ourselves.
1579 36820 : dof_id_type parallel_max_node_id = 0;
1580 7850961 : for (const auto & node : mesh.node_ptr_range())
1581 7814141 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
1582 7814141 : node->id()+1);
1583 36820 : mesh.comm().max(parallel_max_node_id);
1584 :
1585 :
1586 73640 : std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1587 :
1588 2229608 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1589 4459216 : mesh.local_elements_end()))
1590 : {
1591 2192788 : libmesh_assert (elem);
1592 :
1593 16866120 : for (auto & node : elem->node_ref_range())
1594 : {
1595 14673332 : dof_id_type nodeid = node.id();
1596 14673332 : node_touched_by_me[nodeid] = true;
1597 : }
1598 : }
1599 73640 : std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1600 36820 : mesh.comm().max(node_touched_by_anyone);
1601 :
1602 3931677 : for (const auto & node : mesh.local_node_ptr_range())
1603 : {
1604 3894857 : libmesh_assert(node);
1605 3894857 : dof_id_type nodeid = node->id();
1606 3894857 : libmesh_assert(!node_touched_by_anyone[nodeid] ||
1607 : node_touched_by_me[nodeid]);
1608 : }
1609 : }
1610 :
1611 :
1612 :
1613 0 : void libmesh_assert_canonical_node_procids (const MeshBase & mesh)
1614 : {
1615 0 : for (const auto & elem : mesh.active_element_ptr_range())
1616 0 : for (auto & node : elem->node_ref_range())
1617 0 : libmesh_assert_equal_to
1618 : (node.processor_id(),
1619 : node.choose_processor_id(node.processor_id(),
1620 : elem->processor_id()));
1621 0 : }
1622 :
1623 :
1624 :
1625 : #ifdef LIBMESH_ENABLE_AMR
1626 1336 : void libmesh_assert_valid_refinement_tree(const MeshBase & mesh)
1627 : {
1628 2672 : LOG_SCOPE("libmesh_assert_valid_refinement_tree()", "MeshTools");
1629 :
1630 68018 : for (const auto & elem : mesh.element_ptr_range())
1631 : {
1632 66682 : libmesh_assert(elem);
1633 66682 : if (elem->has_children())
1634 27268 : for (auto & child : elem->child_ref_range())
1635 23216 : if (&child != remote_elem)
1636 20800 : libmesh_assert_equal_to (child.parent(), elem);
1637 66682 : if (elem->active())
1638 : {
1639 62630 : libmesh_assert(!elem->ancestor());
1640 62630 : libmesh_assert(!elem->subactive());
1641 : }
1642 4052 : else if (elem->ancestor())
1643 : {
1644 4052 : libmesh_assert(!elem->subactive());
1645 : }
1646 : else
1647 0 : libmesh_assert(elem->subactive());
1648 :
1649 66682 : if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1650 0 : libmesh_assert_greater(elem->p_level(), 0);
1651 : }
1652 1336 : }
1653 : #else
1654 : void libmesh_assert_valid_refinement_tree(const MeshBase &)
1655 : {
1656 : }
1657 : #endif // LIBMESH_ENABLE_AMR
1658 :
1659 : #endif // !NDEBUG
1660 :
1661 :
1662 :
1663 : #ifdef DEBUG
1664 :
1665 0 : void libmesh_assert_no_links_to_elem(const MeshBase & mesh,
1666 : const Elem * bad_elem)
1667 : {
1668 0 : for (const auto & elem : mesh.element_ptr_range())
1669 : {
1670 0 : libmesh_assert (elem);
1671 0 : libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1672 0 : for (auto n : elem->neighbor_ptr_range())
1673 0 : libmesh_assert_not_equal_to (n, bad_elem);
1674 :
1675 : #ifdef LIBMESH_ENABLE_AMR
1676 0 : if (elem->has_children())
1677 0 : for (auto & child : elem->child_ref_range())
1678 0 : libmesh_assert_not_equal_to (&child, bad_elem);
1679 : #endif
1680 : }
1681 0 : }
1682 :
1683 :
1684 570 : void libmesh_assert_equal_points (const MeshBase & mesh)
1685 : {
1686 1140 : LOG_SCOPE("libmesh_assert_equal_points()", "MeshTools");
1687 :
1688 570 : dof_id_type pmax_node_id = mesh.max_node_id();
1689 570 : mesh.comm().max(pmax_node_id);
1690 :
1691 922414 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1692 : {
1693 921844 : const Point * p = mesh.query_node_ptr(i);
1694 :
1695 921844 : libmesh_assert(mesh.comm().semiverify(p));
1696 : }
1697 570 : }
1698 :
1699 :
1700 570 : void libmesh_assert_equal_connectivity (const MeshBase & mesh)
1701 : {
1702 1140 : LOG_SCOPE("libmesh_assert_equal_connectivity()", "MeshTools");
1703 :
1704 570 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1705 570 : mesh.comm().max(pmax_elem_id);
1706 :
1707 1144326 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1708 : {
1709 1143756 : const Elem * e = mesh.query_elem_ptr(i);
1710 :
1711 2287512 : std::vector<dof_id_type> nodes;
1712 1143756 : if (e)
1713 5792488 : for (auto n : e->node_index_range())
1714 4648848 : nodes.push_back(e->node_id(n));
1715 :
1716 1143756 : libmesh_assert(mesh.comm().semiverify(e ? &nodes : nullptr));
1717 : }
1718 570 : }
1719 :
1720 :
1721 0 : void libmesh_assert_connected_nodes (const MeshBase & mesh)
1722 : {
1723 0 : LOG_SCOPE("libmesh_assert_connected_nodes()", "MeshTools");
1724 :
1725 0 : std::set<const Node *> used_nodes;
1726 :
1727 0 : for (const auto & elem : mesh.element_ptr_range())
1728 : {
1729 0 : libmesh_assert (elem);
1730 :
1731 0 : for (auto & n : elem->node_ref_range())
1732 0 : used_nodes.insert(&n);
1733 : }
1734 :
1735 0 : for (const auto & node : mesh.node_ptr_range())
1736 : {
1737 0 : libmesh_assert(node);
1738 0 : libmesh_assert(used_nodes.count(node));
1739 : }
1740 0 : }
1741 :
1742 :
1743 :
1744 26922 : void libmesh_assert_valid_constraint_rows (const MeshBase & mesh)
1745 : {
1746 26922 : libmesh_parallel_only(mesh.comm());
1747 :
1748 26922 : const auto & constraint_rows = mesh.get_constraint_rows();
1749 :
1750 26922 : bool have_constraint_rows = !constraint_rows.empty();
1751 26922 : mesh.comm().max(have_constraint_rows);
1752 26922 : if (!have_constraint_rows)
1753 26882 : return;
1754 :
1755 12700 : for (auto & row : constraint_rows)
1756 : {
1757 12660 : const Node * node = row.first;
1758 12660 : libmesh_assert(node == mesh.node_ptr(node->id()));
1759 :
1760 44126 : for (auto & pr : row.second)
1761 : {
1762 31466 : const Elem * spline_elem = pr.first.first;
1763 31466 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1764 : }
1765 : }
1766 :
1767 40 : dof_id_type pmax_node_id = mesh.max_node_id();
1768 40 : mesh.comm().max(pmax_node_id);
1769 :
1770 16990 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1771 : {
1772 16950 : const Node * node = mesh.query_node_ptr(i);
1773 :
1774 16950 : bool have_constraint = constraint_rows.count(node);
1775 :
1776 16950 : const std::size_t my_n_constraints = have_constraint ?
1777 12660 : libmesh_map_find(constraint_rows, node).size() : std::size_t(-1);
1778 16950 : const std::size_t * n_constraints = node ?
1779 16950 : &my_n_constraints : nullptr;
1780 :
1781 16950 : libmesh_assert(mesh.comm().semiverify(n_constraints));
1782 : }
1783 : }
1784 :
1785 :
1786 :
1787 64006 : void libmesh_assert_valid_boundary_ids(const MeshBase & mesh)
1788 : {
1789 64006 : LOG_SCOPE("libmesh_assert_valid_boundary_ids()", "MeshTools");
1790 :
1791 64006 : if (mesh.n_processors() == 1)
1792 1164 : return;
1793 :
1794 62842 : libmesh_parallel_only(mesh.comm());
1795 :
1796 62842 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1797 :
1798 62842 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1799 62842 : mesh.comm().max(pmax_elem_id);
1800 :
1801 4191636 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1802 : {
1803 4128794 : const Elem * elem = mesh.query_elem_ptr(i);
1804 4128794 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1805 4128794 : const unsigned int my_n_edges = elem ? elem->n_edges() : 0;
1806 4128794 : const unsigned int my_n_sides = elem ? elem->n_sides() : 0;
1807 : unsigned int
1808 4128794 : n_nodes = my_n_nodes,
1809 4128794 : n_edges = my_n_edges,
1810 4128794 : n_sides = my_n_sides;
1811 :
1812 4128794 : mesh.comm().max(n_nodes);
1813 4128794 : mesh.comm().max(n_edges);
1814 4128794 : mesh.comm().max(n_sides);
1815 :
1816 4128794 : if (elem)
1817 : {
1818 4052609 : libmesh_assert_equal_to(my_n_nodes, n_nodes);
1819 4052609 : libmesh_assert_equal_to(my_n_edges, n_edges);
1820 4052609 : libmesh_assert_equal_to(my_n_sides, n_sides);
1821 : }
1822 :
1823 : // Let's test all IDs on the element with one communication
1824 : // rather than n_nodes + n_edges + n_sides communications, to
1825 : // cut down on latency in dbg modes.
1826 8257588 : std::vector<boundary_id_type> all_bcids;
1827 :
1828 33598754 : for (unsigned int n=0; n != n_nodes; ++n)
1829 : {
1830 58939920 : std::vector<boundary_id_type> bcids;
1831 29469960 : if (elem)
1832 : {
1833 29295094 : boundary_info.boundary_ids(elem->node_ptr(n), bcids);
1834 :
1835 : // Ordering of boundary ids shouldn't matter
1836 29295094 : std::sort(bcids.begin(), bcids.end());
1837 : }
1838 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1839 :
1840 29469960 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1841 58939920 : bcids.end());
1842 : // Separator
1843 29469960 : all_bcids.push_back(BoundaryInfo::invalid_id);
1844 : }
1845 :
1846 25685780 : for (unsigned short e=0; e != n_edges; ++e)
1847 : {
1848 43113972 : std::vector<boundary_id_type> bcids;
1849 :
1850 21556986 : if (elem)
1851 : {
1852 21439690 : boundary_info.edge_boundary_ids(elem, e, bcids);
1853 :
1854 : // Ordering of boundary ids shouldn't matter
1855 21439690 : std::sort(bcids.begin(), bcids.end());
1856 : }
1857 :
1858 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1859 :
1860 21556986 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1861 43113972 : bcids.end());
1862 : // Separator
1863 21556986 : all_bcids.push_back(BoundaryInfo::invalid_id);
1864 :
1865 21556986 : if (elem)
1866 : {
1867 21439690 : boundary_info.raw_edge_boundary_ids(elem, e, bcids);
1868 :
1869 : // Ordering of boundary ids shouldn't matter
1870 21439690 : std::sort(bcids.begin(), bcids.end());
1871 :
1872 21439690 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1873 42879380 : bcids.end());
1874 : // Separator
1875 21439690 : all_bcids.push_back(BoundaryInfo::invalid_id);
1876 : }
1877 :
1878 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1879 : }
1880 :
1881 20471752 : for (unsigned short s=0; s != n_sides; ++s)
1882 : {
1883 32685916 : std::vector<boundary_id_type> bcids;
1884 :
1885 16342958 : if (elem)
1886 : {
1887 16260150 : boundary_info.boundary_ids(elem, s, bcids);
1888 :
1889 : // Ordering of boundary ids shouldn't matter
1890 16260150 : std::sort(bcids.begin(), bcids.end());
1891 :
1892 16260150 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1893 32520300 : bcids.end());
1894 : // Separator
1895 16260150 : all_bcids.push_back(BoundaryInfo::invalid_id);
1896 : }
1897 :
1898 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1899 :
1900 16342958 : if (elem)
1901 : {
1902 16260150 : boundary_info.raw_boundary_ids(elem, s, bcids);
1903 :
1904 : // Ordering of boundary ids shouldn't matter
1905 16260150 : std::sort(bcids.begin(), bcids.end());
1906 :
1907 16260150 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1908 32520300 : bcids.end());
1909 : // Separator
1910 16260150 : all_bcids.push_back(BoundaryInfo::invalid_id);
1911 : }
1912 :
1913 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1914 : }
1915 :
1916 12386382 : for (unsigned short sf=0; sf != 2; ++sf)
1917 : {
1918 16515176 : std::vector<boundary_id_type> bcids;
1919 :
1920 8257588 : if (elem)
1921 : {
1922 8105218 : boundary_info.shellface_boundary_ids(elem, sf, bcids);
1923 :
1924 : // Ordering of boundary ids shouldn't matter
1925 8105218 : std::sort(bcids.begin(), bcids.end());
1926 :
1927 8105218 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1928 16210436 : bcids.end());
1929 : // Separator
1930 8105218 : all_bcids.push_back(BoundaryInfo::invalid_id);
1931 : }
1932 :
1933 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1934 :
1935 8257588 : if (elem)
1936 : {
1937 8105218 : boundary_info.raw_shellface_boundary_ids(elem, sf, bcids);
1938 :
1939 : // Ordering of boundary ids shouldn't matter
1940 8105218 : std::sort(bcids.begin(), bcids.end());
1941 :
1942 8105218 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1943 16210436 : bcids.end());
1944 : // Separator
1945 8105218 : all_bcids.push_back(BoundaryInfo::invalid_id);
1946 : }
1947 :
1948 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1949 : }
1950 :
1951 4128794 : libmesh_assert(mesh.comm().semiverify
1952 : (elem ? &all_bcids : nullptr));
1953 : }
1954 : }
1955 :
1956 :
1957 11176 : void libmesh_assert_valid_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1958 : {
1959 11176 : LOG_SCOPE("libmesh_assert_valid_dof_ids()", "MeshTools");
1960 :
1961 11176 : if (mesh.n_processors() == 1)
1962 0 : return;
1963 :
1964 11176 : libmesh_parallel_only(mesh.comm());
1965 :
1966 11176 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1967 11176 : mesh.comm().max(pmax_elem_id);
1968 :
1969 1061176 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1970 1050000 : assert_semiverify_dofobj(mesh.comm(),
1971 1050000 : mesh.query_elem_ptr(i),
1972 : sysnum);
1973 :
1974 11176 : dof_id_type pmax_node_id = mesh.max_node_id();
1975 11176 : mesh.comm().max(pmax_node_id);
1976 :
1977 1988784 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1978 1977608 : assert_semiverify_dofobj(mesh.comm(),
1979 1977608 : mesh.query_node_ptr(i),
1980 : sysnum);
1981 : }
1982 :
1983 :
1984 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1985 52776 : void libmesh_assert_valid_unique_ids(const MeshBase & mesh)
1986 : {
1987 105552 : LOG_SCOPE("libmesh_assert_valid_unique_ids()", "MeshTools");
1988 :
1989 52776 : libmesh_parallel_only(mesh.comm());
1990 :
1991 : // Storage for semi-local DofObject ids.
1992 105552 : std::unordered_set<unique_id_type> semilocal_unique_ids;
1993 :
1994 52776 : auto gather_elem_ids = [&]()
1995 : {
1996 2587985 : for (auto const & elem : mesh.active_element_ptr_range())
1997 : {
1998 2535209 : auto [it, inserted] = semilocal_unique_ids.insert(elem->unique_id());
1999 2535209 : libmesh_assert(inserted);
2000 2535209 : libmesh_ignore(it);
2001 : }
2002 52776 : };
2003 :
2004 52776 : auto gather_node_ids = [&]()
2005 : {
2006 6175137 : for (auto const & node : mesh.node_ptr_range())
2007 : {
2008 6122361 : auto [it, inserted] = semilocal_unique_ids.insert(node->unique_id());
2009 6122361 : libmesh_assert(inserted);
2010 6122361 : libmesh_ignore(it);
2011 : }
2012 52776 : };
2013 :
2014 52776 : auto verify_elems = [&]()
2015 : {
2016 52776 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2017 52776 : mesh.comm().max(pmax_elem_id);
2018 :
2019 3145816 : for (auto i : make_range(pmax_elem_id))
2020 : {
2021 3093040 : const Elem * elem = mesh.query_elem_ptr(i);
2022 3093040 : assert_dofobj_unique_id(mesh.comm(), elem, semilocal_unique_ids);
2023 : }
2024 52776 : };
2025 :
2026 52776 : auto verify_nodes = [&]()
2027 : {
2028 52776 : dof_id_type pmax_node_id = mesh.max_node_id();
2029 52776 : mesh.comm().max(pmax_node_id);
2030 :
2031 6430832 : for (auto i : make_range(pmax_node_id))
2032 : {
2033 6378056 : const Node * node = mesh.query_node_ptr(i);
2034 6378056 : assert_dofobj_unique_id(mesh.comm(), node, semilocal_unique_ids);
2035 : }
2036 52776 : };
2037 :
2038 52776 : if (!mesh.allow_node_and_elem_unique_id_overlap())
2039 : {
2040 : // First collect all the unique_ids we can see and make sure there's
2041 : // no duplicates
2042 52776 : gather_elem_ids();
2043 52776 : gather_node_ids();
2044 :
2045 : // Then make sure elements/nodes are all in sync and remote
2046 : // elements don't duplicate semilocal
2047 52776 : verify_elems();
2048 52776 : verify_nodes();
2049 : }
2050 : else
2051 : {
2052 : // If the mesh allows Node and Elem unique_ids to overlap, then we only
2053 : // check for validity and uniqueness of an Elem (resp. Node) unique id
2054 : // within the set of Elem (resp. Node) unique_ids.
2055 0 : gather_elem_ids();
2056 0 : verify_elems();
2057 :
2058 : // Clear id list before checking Nodes
2059 0 : semilocal_unique_ids.clear();
2060 :
2061 : // Finally, check Nodes
2062 0 : gather_node_ids();
2063 0 : verify_nodes();
2064 : }
2065 52776 : }
2066 : #endif
2067 :
2068 0 : void libmesh_assert_consistent_distributed(const MeshBase & mesh)
2069 : {
2070 0 : libmesh_parallel_only(mesh.comm());
2071 :
2072 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2073 0 : mesh.comm().max(parallel_max_elem_id);
2074 :
2075 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2076 : {
2077 0 : const Elem * elem = mesh.query_elem_ptr(i);
2078 : processor_id_type pid =
2079 0 : elem ? elem->processor_id() : DofObject::invalid_processor_id;
2080 0 : mesh.comm().min(pid);
2081 0 : libmesh_assert(elem || pid != mesh.processor_id());
2082 : }
2083 :
2084 0 : dof_id_type parallel_max_node_id = mesh.max_node_id();
2085 0 : mesh.comm().max(parallel_max_node_id);
2086 :
2087 0 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2088 : {
2089 0 : const Node * node = mesh.query_node_ptr(i);
2090 : processor_id_type pid =
2091 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2092 0 : mesh.comm().min(pid);
2093 0 : libmesh_assert(node || pid != mesh.processor_id());
2094 : }
2095 0 : }
2096 :
2097 :
2098 0 : void libmesh_assert_consistent_distributed_nodes(const MeshBase & mesh)
2099 : {
2100 0 : libmesh_parallel_only(mesh.comm());
2101 0 : auto locator = mesh.sub_point_locator();
2102 :
2103 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2104 0 : mesh.comm().max(parallel_max_elem_id);
2105 :
2106 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2107 : {
2108 0 : const Elem * elem = mesh.query_elem_ptr(i);
2109 :
2110 0 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2111 0 : unsigned int n_nodes = my_n_nodes;
2112 0 : mesh.comm().max(n_nodes);
2113 :
2114 0 : if (n_nodes)
2115 0 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2116 :
2117 0 : for (unsigned int n=0; n != n_nodes; ++n)
2118 : {
2119 0 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2120 : processor_id_type pid =
2121 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2122 0 : mesh.comm().min(pid);
2123 0 : libmesh_assert(node || pid != mesh.processor_id());
2124 : }
2125 : }
2126 0 : }
2127 :
2128 :
2129 :
2130 : template <>
2131 24018 : void libmesh_assert_parallel_consistent_procids<Elem>(const MeshBase & mesh)
2132 : {
2133 24018 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2134 :
2135 24018 : if (mesh.n_processors() == 1)
2136 2 : return;
2137 :
2138 24016 : libmesh_parallel_only(mesh.comm());
2139 :
2140 : // Some code (looking at you, stitch_meshes) modifies DofObject ids
2141 : // without keeping max_elem_id()/max_node_id() consistent, but
2142 : // that's done in a safe way for performance reasons, so we'll play
2143 : // along and just figure out new max ids ourselves.
2144 24016 : dof_id_type parallel_max_elem_id = 0;
2145 2817340 : for (const auto & elem : mesh.element_ptr_range())
2146 2793324 : parallel_max_elem_id = std::max<dof_id_type>(parallel_max_elem_id,
2147 2793324 : elem->id()+1);
2148 24016 : mesh.comm().max(parallel_max_elem_id);
2149 :
2150 : // Check processor ids for consistency between processors
2151 :
2152 2884390 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2153 : {
2154 2860374 : const Elem * elem = mesh.query_elem_ptr(i);
2155 :
2156 : processor_id_type min_id =
2157 2793324 : elem ? elem->processor_id() :
2158 2860374 : std::numeric_limits<processor_id_type>::max();
2159 2860374 : mesh.comm().min(min_id);
2160 :
2161 : processor_id_type max_id =
2162 2793324 : elem ? elem->processor_id() :
2163 2860374 : std::numeric_limits<processor_id_type>::min();
2164 2860374 : mesh.comm().max(max_id);
2165 :
2166 2860374 : if (elem)
2167 : {
2168 2793324 : libmesh_assert_equal_to (min_id, elem->processor_id());
2169 2793324 : libmesh_assert_equal_to (max_id, elem->processor_id());
2170 : }
2171 :
2172 2860374 : if (min_id == mesh.processor_id())
2173 1357933 : libmesh_assert(elem);
2174 : }
2175 : }
2176 :
2177 :
2178 :
2179 76 : void libmesh_assert_parallel_consistent_new_node_procids(const MeshBase & mesh)
2180 : {
2181 76 : LOG_SCOPE("libmesh_assert_parallel_consistent_new_node_procids()", "MeshTools");
2182 :
2183 76 : if (mesh.n_processors() == 1)
2184 0 : return;
2185 :
2186 76 : libmesh_parallel_only(mesh.comm());
2187 :
2188 : // We want this test to hit every node when called even after nodes
2189 : // have been added asynchronously but before everything has been
2190 : // renumbered.
2191 76 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2192 76 : mesh.comm().max(parallel_max_elem_id);
2193 :
2194 152 : std::vector<bool> elem_touched_by_anyone(parallel_max_elem_id, false);
2195 :
2196 42408 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2197 : {
2198 42332 : const Elem * elem = mesh.query_elem_ptr(i);
2199 :
2200 42332 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2201 42332 : unsigned int n_nodes = my_n_nodes;
2202 42332 : mesh.comm().max(n_nodes);
2203 :
2204 42332 : if (n_nodes)
2205 15840 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2206 :
2207 172564 : for (unsigned int n=0; n != n_nodes; ++n)
2208 : {
2209 130232 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2210 130232 : const processor_id_type pid = node ? node->processor_id() : 0;
2211 130232 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2212 : }
2213 : }
2214 : }
2215 :
2216 : template <>
2217 37933 : void libmesh_assert_parallel_consistent_procids<Node>(const MeshBase & mesh)
2218 : {
2219 37933 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2220 :
2221 37933 : if (mesh.n_processors() == 1)
2222 1 : return;
2223 :
2224 37932 : libmesh_parallel_only(mesh.comm());
2225 :
2226 : // We want this test to be valid even when called even after nodes
2227 : // have been added asynchronously but before they're renumbered
2228 : //
2229 : // Plus, some code (looking at you, stitch_meshes) modifies
2230 : // DofObject ids without keeping max_elem_id()/max_node_id()
2231 : // consistent, but that's done in a safe way for performance
2232 : // reasons, so we'll play along and just figure out new max ids
2233 : // ourselves.
2234 37932 : dof_id_type parallel_max_node_id = 0;
2235 8287727 : for (const auto & node : mesh.node_ptr_range())
2236 8249795 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
2237 8249795 : node->id()+1);
2238 37932 : mesh.comm().max(parallel_max_node_id);
2239 :
2240 75864 : std::vector<bool> node_touched_by_anyone(parallel_max_node_id, false);
2241 :
2242 2430581 : for (const auto & elem : as_range(mesh.local_elements_begin(),
2243 4861162 : mesh.local_elements_end()))
2244 : {
2245 2392649 : libmesh_assert (elem);
2246 :
2247 18172622 : for (auto & node : elem->node_ref_range())
2248 : {
2249 15779973 : dof_id_type nodeid = node.id();
2250 15779973 : node_touched_by_anyone[nodeid] = true;
2251 : }
2252 : }
2253 37932 : mesh.comm().max(node_touched_by_anyone);
2254 :
2255 : // Check processor ids for consistency between processors
2256 : // on any node an element touches
2257 8706252 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2258 : {
2259 8668320 : if (!node_touched_by_anyone[i])
2260 420946 : continue;
2261 :
2262 8247374 : const Node * node = mesh.query_node_ptr(i);
2263 8247374 : const processor_id_type pid = node ? node->processor_id() : 0;
2264 :
2265 8247374 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2266 : }
2267 : }
2268 :
2269 :
2270 :
2271 : #ifdef LIBMESH_ENABLE_AMR
2272 0 : void libmesh_assert_valid_refinement_flags(const MeshBase & mesh)
2273 : {
2274 0 : LOG_SCOPE("libmesh_assert_valid_refinement_flags()", "MeshTools");
2275 :
2276 0 : libmesh_parallel_only(mesh.comm());
2277 0 : if (mesh.n_processors() == 1)
2278 0 : return;
2279 :
2280 0 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2281 0 : mesh.comm().max(pmax_elem_id);
2282 :
2283 0 : std::vector<unsigned char> my_elem_h_state(pmax_elem_id, 255);
2284 0 : std::vector<unsigned char> my_elem_p_state(pmax_elem_id, 255);
2285 :
2286 0 : for (const auto & elem : mesh.element_ptr_range())
2287 : {
2288 0 : libmesh_assert (elem);
2289 0 : dof_id_type elemid = elem->id();
2290 :
2291 0 : my_elem_h_state[elemid] =
2292 0 : static_cast<unsigned char>(elem->refinement_flag());
2293 :
2294 0 : my_elem_p_state[elemid] =
2295 0 : static_cast<unsigned char>(elem->p_refinement_flag());
2296 : }
2297 0 : std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
2298 0 : mesh.comm().min(min_elem_h_state);
2299 :
2300 0 : std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
2301 0 : mesh.comm().min(min_elem_p_state);
2302 :
2303 0 : for (dof_id_type i=0; i!= pmax_elem_id; ++i)
2304 : {
2305 0 : libmesh_assert(my_elem_h_state[i] == 255 ||
2306 : my_elem_h_state[i] == min_elem_h_state[i]);
2307 0 : libmesh_assert(my_elem_p_state[i] == 255 ||
2308 : my_elem_p_state[i] == min_elem_p_state[i]);
2309 : }
2310 : }
2311 : #else
2312 : void libmesh_assert_valid_refinement_flags(const MeshBase &)
2313 : {
2314 : }
2315 : #endif // LIBMESH_ENABLE_AMR
2316 :
2317 :
2318 :
2319 21120 : void libmesh_assert_valid_neighbors(const MeshBase & mesh,
2320 : bool assert_valid_remote_elems)
2321 : {
2322 21120 : LOG_SCOPE("libmesh_assert_valid_neighbors()", "MeshTools");
2323 :
2324 2893986 : for (const auto & elem : mesh.element_ptr_range())
2325 : {
2326 2872866 : libmesh_assert (elem);
2327 2872866 : elem->libmesh_assert_valid_neighbors();
2328 : }
2329 :
2330 21120 : if (mesh.n_processors() == 1)
2331 812 : return;
2332 :
2333 20308 : libmesh_parallel_only(mesh.comm());
2334 :
2335 20308 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2336 20308 : mesh.comm().max(pmax_elem_id);
2337 :
2338 2989962 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
2339 : {
2340 2969654 : const Elem * elem = mesh.query_elem_ptr(i);
2341 :
2342 2969654 : const unsigned int my_n_neigh = elem ? elem->n_neighbors() : 0;
2343 2969654 : unsigned int n_neigh = my_n_neigh;
2344 2969654 : mesh.comm().max(n_neigh);
2345 2969654 : if (elem)
2346 2867762 : libmesh_assert_equal_to (my_n_neigh, n_neigh);
2347 :
2348 14320524 : for (unsigned int n = 0; n != n_neigh; ++n)
2349 : {
2350 11350870 : dof_id_type my_neighbor = DofObject::invalid_id;
2351 11350870 : dof_id_type * p_my_neighbor = nullptr;
2352 :
2353 : // If we have a non-remote_elem neighbor link, then we can
2354 : // verify it.
2355 11350870 : if (elem && elem->neighbor_ptr(n) != remote_elem)
2356 : {
2357 11258076 : p_my_neighbor = &my_neighbor;
2358 11258076 : if (elem->neighbor_ptr(n))
2359 10621260 : my_neighbor = elem->neighbor_ptr(n)->id();
2360 :
2361 : // But wait - if we haven't set remote_elem links yet then
2362 : // some nullptr links on ghost elements might be
2363 : // future-remote_elem links, so we can't verify those.
2364 22523096 : if (!assert_valid_remote_elems &&
2365 11259726 : !elem->neighbor_ptr(n) &&
2366 1650 : elem->processor_id() != mesh.processor_id())
2367 773 : p_my_neighbor = nullptr;
2368 : }
2369 11350870 : libmesh_assert(mesh.comm().semiverify(p_my_neighbor));
2370 : }
2371 : }
2372 : }
2373 : #endif // DEBUG
2374 :
2375 :
2376 :
2377 : // Functors for correct_node_proc_ids
2378 : namespace {
2379 :
2380 : typedef std::unordered_map<dof_id_type, processor_id_type> proc_id_map_type;
2381 :
2382 : struct SyncNodeSet
2383 : {
2384 : typedef unsigned char datum; // bool but without bit twiddling issues
2385 :
2386 182 : SyncNodeSet(std::unordered_set<const Node *> & _set,
2387 9199 : MeshBase & _mesh) :
2388 9199 : node_set(_set), mesh(_mesh) {}
2389 :
2390 : std::unordered_set<const Node *> & node_set;
2391 :
2392 : MeshBase & mesh;
2393 :
2394 : // ------------------------------------------------------------
2395 41122 : void gather_data (const std::vector<dof_id_type> & ids,
2396 : std::vector<datum> & data)
2397 : {
2398 : // Find whether each requested node belongs in the set
2399 41122 : data.resize(ids.size());
2400 :
2401 1867830 : for (auto i : index_range(ids))
2402 : {
2403 1826708 : const dof_id_type id = ids[i];
2404 :
2405 : // We'd better have every node we're asked for
2406 1826708 : Node * node = mesh.node_ptr(id);
2407 :
2408 : // Return if the node is in the set.
2409 2954958 : data[i] = node_set.count(node);
2410 : }
2411 41122 : }
2412 :
2413 : // ------------------------------------------------------------
2414 41122 : bool act_on_data (const std::vector<dof_id_type> & ids,
2415 : const std::vector<datum> in_set)
2416 : {
2417 167 : bool data_changed = false;
2418 :
2419 : // Add nodes we've been informed of to our own set
2420 1867830 : for (auto i : index_range(ids))
2421 : {
2422 1826708 : if (in_set[i])
2423 : {
2424 1196849 : Node * node = mesh.node_ptr(ids[i]);
2425 1196849 : if (!node_set.count(node))
2426 : {
2427 256 : node_set.insert(node);
2428 256 : data_changed = true;
2429 : }
2430 : }
2431 : }
2432 :
2433 41122 : return data_changed;
2434 : }
2435 : };
2436 :
2437 :
2438 8835 : struct NodesNotInSet
2439 : {
2440 182 : NodesNotInSet(const std::unordered_set<const Node *> _set)
2441 9017 : : node_set(_set) {}
2442 :
2443 479564 : bool operator() (const Node * node) const
2444 : {
2445 959128 : if (node_set.count(node))
2446 287060 : return false;
2447 192504 : return true;
2448 : }
2449 :
2450 : const std::unordered_set<const Node *> node_set;
2451 : };
2452 :
2453 :
2454 : struct SyncProcIdsFromMap
2455 : {
2456 : typedef processor_id_type datum;
2457 :
2458 218 : SyncProcIdsFromMap(const proc_id_map_type & _map,
2459 38187 : MeshBase & _mesh) :
2460 38187 : new_proc_ids(_map), mesh(_mesh) {}
2461 :
2462 : const proc_id_map_type & new_proc_ids;
2463 :
2464 : MeshBase & mesh;
2465 :
2466 : // ------------------------------------------------------------
2467 132325 : void gather_data (const std::vector<dof_id_type> & ids,
2468 : std::vector<datum> & data)
2469 : {
2470 : // Find the new processor id of each requested node
2471 132325 : data.resize(ids.size());
2472 :
2473 7264017 : for (auto i : index_range(ids))
2474 : {
2475 7131692 : const dof_id_type id = ids[i];
2476 :
2477 : // Return the node's new processor id if it has one, or its
2478 : // old processor id if not.
2479 7131692 : if (const auto it = new_proc_ids.find(id);
2480 47012 : it != new_proc_ids.end())
2481 6500176 : data[i] = it->second;
2482 : else
2483 : {
2484 : // We'd better find every node we're asked for
2485 631516 : const Node & node = mesh.node_ref(id);
2486 631516 : data[i] = node.processor_id();
2487 : }
2488 : }
2489 132325 : }
2490 :
2491 : // ------------------------------------------------------------
2492 132325 : void act_on_data (const std::vector<dof_id_type> & ids,
2493 : const std::vector<datum> proc_ids)
2494 : {
2495 : // Set the node processor ids we've now been informed of
2496 7264017 : for (auto i : index_range(ids))
2497 : {
2498 7131692 : Node & node = mesh.node_ref(ids[i]);
2499 7131692 : node.processor_id() = proc_ids[i];
2500 : }
2501 132325 : }
2502 : };
2503 : }
2504 :
2505 :
2506 :
2507 38187 : void correct_node_proc_ids (MeshBase & mesh)
2508 : {
2509 654 : LOG_SCOPE("correct_node_proc_ids()","MeshTools");
2510 :
2511 : // This function must be run on all processors at once
2512 218 : libmesh_parallel_only(mesh.comm());
2513 :
2514 : // We require all processors to agree on nodal processor ids before
2515 : // going into this algorithm.
2516 : #ifdef DEBUG
2517 218 : libmesh_assert_parallel_consistent_procids<Node>(mesh);
2518 : #endif
2519 :
2520 : // If we have any unpartitioned elements at this
2521 : // stage there is a problem
2522 218 : libmesh_assert (n_elem(mesh.unpartitioned_elements_begin(),
2523 : mesh.unpartitioned_elements_end()) == 0);
2524 :
2525 : // Fix nodes' processor ids. Coarsening may have left us with nodes
2526 : // which are no longer touched by any elements of the same processor
2527 : // id, and for DofMap to work we need to fix that.
2528 :
2529 : // This is harder now that libMesh no longer requires a distributed
2530 : // mesh to ghost all nodal neighbors: it is possible for two active
2531 : // elements on two different processors to share the same node in
2532 : // such a way that neither processor knows the others' element
2533 : // exists!
2534 :
2535 : // While we're at it, if this mesh is configured to allow
2536 : // repartitioning, we'll repartition *all* the nodes' processor ids
2537 : // using the canonical Node heuristic, to try and improve DoF load
2538 : // balancing. But if the mesh is disallowing repartitioning, we
2539 : // won't touch processor_id on any node where it's valid, regardless
2540 : // of whether or not it's canonical.
2541 218 : bool repartition_all_nodes = !mesh.skip_noncritical_partitioning();
2542 436 : std::unordered_set<const Node *> valid_nodes;
2543 :
2544 : // If we aren't allowed to repartition, then we're going to leave
2545 : // every node we can at its current processor_id, and *only*
2546 : // repartition the nodes whose current processor id is incompatible
2547 : // with DoFMap (because it doesn't touch an active element, e.g. due
2548 : // to coarsening)
2549 38187 : if (!repartition_all_nodes)
2550 : {
2551 1158938 : for (const auto & elem : mesh.active_element_ptr_range())
2552 5946053 : for (const auto & node : elem->node_ref_range())
2553 5305120 : if (elem->processor_id() == node.processor_id())
2554 4909802 : valid_nodes.insert(&node);
2555 :
2556 182 : SyncNodeSet syncv(valid_nodes, mesh);
2557 :
2558 : Parallel::sync_dofobject_data_by_id
2559 18216 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), syncv);
2560 : }
2561 :
2562 : // We build up a set of compatible processor ids for each node
2563 436 : proc_id_map_type new_proc_ids;
2564 :
2565 8999106 : for (auto & elem : mesh.active_element_ptr_range())
2566 : {
2567 4500525 : processor_id_type pid = elem->processor_id();
2568 :
2569 44099773 : for (auto & node : elem->node_ref_range())
2570 : {
2571 39599248 : const dof_id_type id = node.id();
2572 39599248 : if (auto it = new_proc_ids.find(id);
2573 350660 : it == new_proc_ids.end())
2574 155764 : new_proc_ids.emplace(id, pid);
2575 : else
2576 28311685 : it->second = node.choose_processor_id(it->second, pid);
2577 : }
2578 37751 : }
2579 :
2580 : // Sort the new pids to push to each processor
2581 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, processor_id_type>>>
2582 436 : ids_to_push;
2583 :
2584 25751072 : for (const auto & node : mesh.node_ptr_range())
2585 13089710 : if (const auto it = std::as_const(new_proc_ids).find(node->id());
2586 13089710 : it != new_proc_ids.end() && node->processor_id() != DofObject::invalid_processor_id)
2587 11325314 : ids_to_push[node->processor_id()].emplace_back(node->id(), /*pid=*/it->second);
2588 :
2589 : auto action_functor =
2590 172786 : [& mesh, & new_proc_ids]
2591 : (processor_id_type,
2592 11755265 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
2593 : {
2594 11461169 : for (const auto & [id, pid] : data)
2595 : {
2596 11287563 : if (const auto it = new_proc_ids.find(id);
2597 155764 : it == new_proc_ids.end())
2598 0 : new_proc_ids.emplace(id, pid);
2599 : else
2600 : {
2601 11287563 : const Node & node = mesh.node_ref(id);
2602 11287563 : it->second = node.choose_processor_id(it->second, pid);
2603 : }
2604 : }
2605 38789 : };
2606 :
2607 : Parallel::push_parallel_vector_data
2608 38187 : (mesh.comm(), ids_to_push, action_functor);
2609 :
2610 : // Now new_proc_ids is correct for every node we used to own. Let's
2611 : // ask every other processor about the nodes they used to own. But
2612 : // first we'll need to keep track of which nodes we used to own,
2613 : // lest we get them confused with nodes we newly own.
2614 436 : std::unordered_set<Node *> ex_local_nodes;
2615 9093990 : for (auto & node : mesh.local_node_ptr_range())
2616 4645558 : if (const auto it = new_proc_ids.find(node->id());
2617 4645558 : it != new_proc_ids.end() && it->second != mesh.processor_id())
2618 37800 : ex_local_nodes.insert(node);
2619 :
2620 218 : SyncProcIdsFromMap sync(new_proc_ids, mesh);
2621 38187 : if (repartition_all_nodes)
2622 : Parallel::sync_dofobject_data_by_id
2623 57940 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync);
2624 : else
2625 : {
2626 182 : NodesNotInSet nnis(valid_nodes);
2627 :
2628 : Parallel::sync_dofobject_data_by_id
2629 18216 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), nnis, sync);
2630 : }
2631 :
2632 : // And finally let's update the nodes we used to own.
2633 42818 : for (const auto & node : ex_local_nodes)
2634 : {
2635 2103 : if (valid_nodes.count(node))
2636 2083 : continue;
2637 :
2638 2548 : const dof_id_type id = node->id();
2639 10 : const proc_id_map_type::iterator it = new_proc_ids.find(id);
2640 10 : libmesh_assert(it != new_proc_ids.end());
2641 2548 : node->processor_id() = it->second;
2642 : }
2643 :
2644 : // We should still have consistent nodal processor ids coming out of
2645 : // this algorithm, but if we're allowed to repartition the mesh then
2646 : // they should be canonically correct too.
2647 : #ifdef DEBUG
2648 218 : libmesh_assert_valid_procids<Node>(mesh);
2649 : //if (repartition_all_nodes)
2650 : // libmesh_assert_canonical_node_procids(mesh);
2651 : #endif
2652 38187 : }
2653 :
2654 :
2655 :
2656 19630 : void Private::globally_renumber_nodes_and_elements (MeshBase & mesh)
2657 : {
2658 19630 : MeshCommunication().assign_global_indices(mesh);
2659 19630 : }
2660 :
2661 : } // namespace MeshTools
2662 :
2663 : } // namespace libMesh
|
