Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/elem.h"
22 : #include "libmesh/elem_range.h"
23 : #include "libmesh/libmesh_logging.h"
24 : #include "libmesh/mesh_base.h"
25 : #include "libmesh/mesh_communication.h"
26 : #include "libmesh/mesh_serializer.h"
27 : #include "libmesh/mesh_tools.h"
28 : #include "libmesh/node_range.h"
29 : #include "libmesh/parallel.h"
30 : #include "libmesh/parallel_algebra.h"
31 : #include "libmesh/parallel_ghost_sync.h"
32 : #include "libmesh/sphere.h"
33 : #include "libmesh/threads.h"
34 : #include "libmesh/enum_to_string.h"
35 : #include "libmesh/enum_elem_type.h"
36 : #include "libmesh/int_range.h"
37 : #include "libmesh/utility.h"
38 : #include "libmesh/boundary_info.h"
39 :
40 : #ifndef NDEBUG
41 : # include "libmesh/remote_elem.h"
42 : #endif
43 :
44 : // C++ includes
45 : #include <limits>
46 : #include <numeric> // for std::accumulate
47 : #include <set>
48 : #include <unordered_map>
49 : #include <unordered_set>
50 :
51 :
52 :
53 : // ------------------------------------------------------------
54 : // anonymous namespace for helper classes
55 : namespace {
56 :
57 : using namespace libMesh;
58 :
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 2297438 : FindBBox () : _bbox()
108 32814 : {}
109 :
110 40652 : FindBBox (FindBBox & other, Threads::split) :
111 40652 : _bbox(other._bbox)
112 13773 : {}
113 :
114 40636 : void operator()(const ConstNodeRange & range)
115 : {
116 5073100 : for (const auto & node : range)
117 : {
118 457357 : libmesh_assert(node);
119 5032464 : _bbox.union_with(*node);
120 : }
121 40636 : }
122 :
123 2349758 : void operator()(const ConstElemRange & range)
124 : {
125 27081965 : for (const auto & elem : range)
126 : {
127 1274884 : libmesh_assert(elem);
128 24732207 : _bbox.union_with(elem->loose_bounding_box());
129 : }
130 2349758 : }
131 :
132 16688 : Point & min() { return _bbox.min(); }
133 :
134 16688 : 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 13773 : void join (const FindBBox & other)
140 : {
141 40652 : _bbox.union_with(other._bbox);
142 26879 : }
143 : #endif
144 :
145 48940 : libMesh::BoundingBox & bbox ()
146 : {
147 48940 : return _bbox;
148 : }
149 :
150 : private:
151 : BoundingBox _bbox;
152 : };
153 :
154 : #ifdef DEBUG
155 2786430 : void assert_semiverify_dofobj(const Parallel::Communicator & communicator,
156 : const DofObject * d,
157 : unsigned int sysnum = libMesh::invalid_uint)
158 : {
159 2786430 : if (d)
160 : {
161 2710542 : const unsigned int n_sys = d->n_systems();
162 :
163 5421084 : std::vector<unsigned int> n_vars (n_sys, 0);
164 5793166 : for (unsigned int s = 0; s != n_sys; ++s)
165 3082624 : if (sysnum == s ||
166 : sysnum == libMesh::invalid_uint)
167 2710542 : n_vars[s] = d->n_vars(s);
168 :
169 : const unsigned int tot_n_vars =
170 2710542 : std::accumulate(n_vars.begin(), n_vars.end(), 0);
171 :
172 5421084 : std::vector<unsigned int> n_comp (tot_n_vars, 0);
173 5421084 : std::vector<dof_id_type> first_dof (tot_n_vars, 0);
174 :
175 5793166 : for (unsigned int s = 0, i=0; s != n_sys; ++s)
176 : {
177 3082624 : if (sysnum != s &&
178 : sysnum != libMesh::invalid_uint)
179 372082 : continue;
180 :
181 9533250 : for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
182 : {
183 6822708 : n_comp[i] = d->n_comp(s,v);
184 6822708 : first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : DofObject::invalid_id;
185 : }
186 : }
187 :
188 2710542 : libmesh_assert(communicator.semiverify(&n_sys));
189 2710542 : libmesh_assert(communicator.semiverify(&n_vars));
190 2710542 : libmesh_assert(communicator.semiverify(&n_comp));
191 2710542 : libmesh_assert(communicator.semiverify(&first_dof));
192 : }
193 : else
194 : {
195 75888 : const unsigned int * p_ui = nullptr;
196 75888 : const std::vector<unsigned int> * p_vui = nullptr;
197 75888 : const std::vector<dof_id_type> * p_vdid = nullptr;
198 :
199 75888 : libmesh_assert(communicator.semiverify(p_ui));
200 75888 : libmesh_assert(communicator.semiverify(p_vui));
201 75888 : libmesh_assert(communicator.semiverify(p_vui));
202 75888 : libmesh_assert(communicator.semiverify(p_vdid));
203 : }
204 2786430 : }
205 :
206 :
207 :
208 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
209 9197454 : 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 9197454 : if (d)
218 : {
219 8815826 : tempmin = tempmax = d->unique_id();
220 : }
221 : else
222 : {
223 381628 : TIMPI::Attributes<unique_id_type>::set_highest(tempmin);
224 381628 : TIMPI::Attributes<unique_id_type>::set_lowest(tempmax);
225 : }
226 9197454 : comm.min(tempmin);
227 9197454 : comm.max(tempmax);
228 18013280 : bool invalid = d && ((d->unique_id() != tempmin) ||
229 8815826 : (d->unique_id() != tempmax));
230 9197454 : comm.max(invalid);
231 :
232 : // First verify that everything is in sync
233 9197454 : libmesh_assert(!invalid);
234 :
235 : // Then verify that any remote id doesn't duplicate a local one.
236 9197454 : if (!d && tempmin == tempmax)
237 50164 : libmesh_assert(!unique_ids.count(tempmin));
238 9197454 : }
239 : #endif // LIBMESH_ENABLE_UNIQUE_ID
240 : #endif // DEBUG
241 :
242 291668 : 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 16432 : 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 1699882 : for (const auto & elem : node_to_elem_vec)
257 : {
258 : // We only care about active elements...
259 1408214 : if (elem->active())
260 : {
261 : // Which local node number is global_id?
262 1368546 : unsigned local_node_number = elem->local_node(global_id);
263 :
264 : // Make sure it was found
265 39668 : libmesh_assert_not_equal_to(local_node_number, libMesh::invalid_uint);
266 :
267 1408214 : 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 1408214 : if (!n_edges)
271 : {
272 1846 : switch (elem->type())
273 : {
274 710 : case EDGE2:
275 : {
276 20 : switch (local_node_number)
277 : {
278 355 : case 0:
279 : // The other node is a nodal neighbor
280 365 : neighbor_set.insert(elem->node_ptr(1));
281 355 : break;
282 :
283 355 : case 1:
284 : // The other node is a nodal neighbor
285 365 : neighbor_set.insert(elem->node_ptr(0));
286 355 : break;
287 :
288 0 : default:
289 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
290 : }
291 20 : break;
292 : }
293 :
294 852 : case EDGE3:
295 : {
296 24 : switch (local_node_number)
297 : {
298 : // The outside nodes have node 2 as a neighbor
299 710 : case 0:
300 : case 1:
301 730 : neighbor_set.insert(elem->node_ptr(2));
302 710 : 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 24 : 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 1408214 : const auto elem_order = Elem::type_to_default_order_map[elem->type()];
354 :
355 : // Index of the current edge
356 39668 : unsigned current_edge = 0;
357 :
358 1408214 : const unsigned short n_nodes = elem->n_nodes();
359 :
360 5830662 : while (current_edge < n_edges)
361 : {
362 : // Find the edge the node is on
363 124576 : bool found_edge = false;
364 11188890 : for (; current_edge<n_edges; ++current_edge)
365 10298266 : if (elem->is_node_on_edge(local_node_number, current_edge))
366 : {
367 99488 : found_edge = true;
368 99488 : break;
369 : }
370 :
371 : // Did we find one?
372 4422448 : if (found_edge)
373 : {
374 3531824 : const Node * node_to_save = nullptr;
375 :
376 : // Find another node in this element on this edge
377 28524392 : for (unsigned other_node_this_edge = 0; other_node_this_edge != n_nodes; other_node_this_edge++)
378 : {
379 46258488 : const bool both_vertices = elem->is_vertex(local_node_number) &&
380 21265920 : elem->is_vertex(other_node_this_edge);
381 32651016 : if ( elem->is_node_on_edge(other_node_this_edge, current_edge) && // On the current edge
382 25038584 : 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 4272128 : (elem_order == 1 || !both_vertices))
385 : {
386 : // We've found a nodal neighbor! Save a pointer to it..
387 3785578 : node_to_save = elem->node_ptr(other_node_this_edge);
388 :
389 : // Make sure we found something
390 106636 : libmesh_assert(node_to_save != nullptr);
391 :
392 3678942 : neighbor_set.insert(node_to_save);
393 : }
394 : }
395 : }
396 :
397 : // Keep looking for edges, node may be on more than one edge
398 4422448 : 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 283452 : neighbors.assign(neighbor_set.begin(), neighbor_set.end());
407 291668 : }
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 3991 : 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 3991 : if (!mesh.is_serial())
462 : libmesh_deprecated();
463 :
464 3991 : nodes_to_elem_map.resize (mesh.max_node_id());
465 :
466 100916 : for (const auto & elem : mesh.element_ptr_range())
467 192936 : for (auto & node : elem->node_ref_range())
468 : {
469 4104 : libmesh_assert_less (node.id(), nodes_to_elem_map.size());
470 4104 : libmesh_assert_less (elem->id(), mesh.n_elem());
471 :
472 147780 : nodes_to_elem_map[node.id()].push_back(elem->id());
473 3763 : }
474 3991 : }
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 243232 : 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 7220 : nodes_to_elem_map.clear();
503 :
504 52440932 : for (const auto & elem : mesh.element_ptr_range())
505 145263145 : for (auto & node : elem->node_ref_range())
506 116532561 : nodes_to_elem_map[node.id()].push_back(elem->id());
507 243232 : }
508 :
509 :
510 :
511 3840 : 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 110 : nodes_to_elem_map.clear();
515 :
516 565218 : for (const auto & elem : mesh.element_ptr_range())
517 1983679 : for (auto & node : elem->node_ref_range())
518 1691127 : nodes_to_elem_map[node.id()].push_back(elem);
519 3840 : }
520 :
521 :
522 :
523 : std::unordered_set<dof_id_type>
524 4628 : find_boundary_nodes(const MeshBase & mesh)
525 : {
526 136 : 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 959170 : for (const auto & elem : mesh.active_element_ptr_range())
531 2721159 : for (auto s : elem->side_index_range())
532 2294362 : if (elem->neighbor_ptr(s) == nullptr) // on the boundary
533 : {
534 139426 : auto nodes_on_side = elem->nodes_on_side(s);
535 :
536 699937 : for (auto & local_id : nodes_on_side)
537 585299 : boundary_nodes.insert(elem->node_ptr(local_id)->id());
538 4356 : }
539 :
540 4628 : return boundary_nodes;
541 : }
542 :
543 : std::unordered_set<dof_id_type>
544 1930 : find_block_boundary_nodes(const MeshBase & mesh)
545 : {
546 60 : 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 778728 : for (const auto & elem : mesh.active_element_ptr_range())
551 2186902 : for (auto s : elem->side_index_range())
552 1840302 : 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 1810 : }
559 :
560 1930 : return block_boundary_nodes;
561 : }
562 :
563 :
564 :
565 : libMesh::BoundingBox
566 1155381 : create_bounding_box (const MeshBase & mesh)
567 : {
568 : // This function must be run on all processors at once
569 16126 : libmesh_parallel_only(mesh.comm());
570 :
571 16126 : FindBBox find_bbox;
572 :
573 : // Start with any unpartitioned elements we know about locally
574 2294636 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(DofObject::invalid_processor_id),
575 1171507 : mesh.pid_elements_end(DofObject::invalid_processor_id)),
576 : find_bbox);
577 :
578 : // And combine with our local elements
579 1155381 : find_bbox.bbox().union_with(create_local_bounding_box(mesh));
580 :
581 : // Compare the bounding boxes across processors
582 1155381 : mesh.comm().min(find_bbox.min());
583 1155381 : mesh.comm().max(find_bbox.max());
584 :
585 1155381 : return find_bbox.bbox();
586 : }
587 :
588 :
589 :
590 : libMesh::BoundingBox
591 19490 : create_nodal_bounding_box (const MeshBase & mesh)
592 : {
593 : // This function must be run on all processors at once
594 562 : libmesh_parallel_only(mesh.comm());
595 :
596 562 : FindBBox find_bbox;
597 :
598 : // Start with any unpartitioned nodes we know about locally
599 38418 : Threads::parallel_reduce (ConstNodeRange (mesh.pid_nodes_begin(DofObject::invalid_processor_id),
600 20052 : mesh.pid_nodes_end(DofObject::invalid_processor_id)),
601 : find_bbox);
602 :
603 : // Add our local nodes
604 38418 : Threads::parallel_reduce (ConstNodeRange (mesh.local_nodes_begin(),
605 38418 : mesh.local_nodes_end()),
606 : find_bbox);
607 :
608 : // Compare the bounding boxes across processors
609 19490 : mesh.comm().min(find_bbox.min());
610 19490 : mesh.comm().max(find_bbox.max());
611 :
612 19490 : 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 1155381 : create_local_bounding_box (const MeshBase & mesh)
632 : {
633 16126 : FindBBox find_bbox;
634 :
635 2294636 : Threads::parallel_reduce (ConstElemRange (mesh.local_elements_begin(),
636 1171507 : mesh.local_elements_end()),
637 : find_bbox);
638 :
639 1155381 : return find_bbox.bbox();
640 : }
641 :
642 :
643 :
644 : libMesh::BoundingBox
645 0 : create_processor_bounding_box (const MeshBase & mesh,
646 : const processor_id_type pid)
647 : {
648 : // This can only be run in parallel, with consistent arguments.
649 0 : libmesh_parallel_only(mesh.comm());
650 0 : libmesh_assert(mesh.comm().verify(pid));
651 :
652 0 : libmesh_assert_less (pid, mesh.n_processors());
653 :
654 0 : FindBBox find_bbox;
655 :
656 0 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
657 0 : mesh.pid_elements_end(pid)),
658 : find_bbox);
659 :
660 : // Compare the bounding boxes across processors
661 0 : mesh.comm().min(find_bbox.min());
662 0 : mesh.comm().max(find_bbox.max());
663 :
664 0 : return find_bbox.bbox();
665 : }
666 :
667 :
668 :
669 : Sphere
670 0 : processor_bounding_sphere (const MeshBase & mesh,
671 : const processor_id_type pid)
672 : {
673 : libMesh::BoundingBox bbox =
674 0 : create_processor_bounding_box(mesh, pid);
675 :
676 0 : const Real diag = (bbox.second - bbox.first).norm();
677 0 : const Point cent = (bbox.second + bbox.first)/2;
678 :
679 0 : return Sphere (cent, .5*diag);
680 : }
681 :
682 :
683 :
684 : libMesh::BoundingBox
685 0 : create_subdomain_bounding_box (const MeshBase & mesh,
686 : const subdomain_id_type sid)
687 : {
688 : // This can only be run in parallel, with consistent arguments.
689 0 : libmesh_parallel_only(mesh.comm());
690 0 : libmesh_assert(mesh.comm().verify(sid));
691 :
692 0 : FindBBox find_bbox;
693 :
694 : Threads::parallel_reduce
695 0 : (ConstElemRange (mesh.active_local_subdomain_elements_begin(sid),
696 0 : mesh.active_local_subdomain_elements_end(sid)),
697 : find_bbox);
698 :
699 : // Compare the bounding boxes across processors
700 0 : mesh.comm().min(find_bbox.min());
701 0 : mesh.comm().max(find_bbox.max());
702 :
703 0 : return find_bbox.bbox();
704 : }
705 :
706 :
707 :
708 : Sphere
709 0 : subdomain_bounding_sphere (const MeshBase & mesh,
710 : const subdomain_id_type sid)
711 : {
712 : libMesh::BoundingBox bbox =
713 0 : create_subdomain_bounding_box(mesh, sid);
714 :
715 0 : const Real diag = (bbox.second - bbox.first).norm();
716 0 : const Point cent = (bbox.second + bbox.first)/2;
717 :
718 0 : return Sphere (cent, .5*diag);
719 : }
720 :
721 :
722 :
723 0 : void elem_types (const MeshBase & mesh,
724 : std::vector<ElemType> & et)
725 : {
726 : // Loop over the the elements. If the current element type isn't in
727 : // the vector, insert it.
728 0 : for (const auto & elem : mesh.element_ptr_range())
729 0 : if (!std::count(et.begin(), et.end(), elem->type()))
730 0 : et.push_back(elem->type());
731 0 : }
732 :
733 :
734 :
735 0 : dof_id_type n_elem_of_type (const MeshBase & mesh,
736 : const ElemType type)
737 : {
738 0 : return static_cast<dof_id_type>(std::distance(mesh.type_elements_begin(type),
739 0 : mesh.type_elements_end (type)));
740 : }
741 :
742 :
743 :
744 0 : dof_id_type n_active_elem_of_type (const MeshBase & mesh,
745 : const ElemType type)
746 : {
747 0 : return static_cast<dof_id_type>(std::distance(mesh.active_type_elements_begin(type),
748 0 : mesh.active_type_elements_end (type)));
749 : }
750 :
751 0 : dof_id_type n_non_subactive_elem_of_type_at_level(const MeshBase & mesh,
752 : const ElemType type,
753 : const unsigned int level)
754 : {
755 0 : dof_id_type cnt = 0;
756 :
757 : // iterate over the elements of the specified type
758 0 : for (const auto & elem : as_range(mesh.type_elements_begin(type),
759 0 : mesh.type_elements_end(type)))
760 0 : if ((elem->level() == level) && !elem->subactive())
761 0 : cnt++;
762 :
763 0 : return cnt;
764 : }
765 :
766 :
767 1425 : unsigned int n_active_local_levels(const MeshBase & mesh)
768 : {
769 1425 : unsigned int nl = 0;
770 :
771 40866 : for (auto & elem : mesh.active_local_element_ptr_range())
772 41845 : nl = std::max(elem->level() + 1, nl);
773 :
774 1425 : return nl;
775 : }
776 :
777 :
778 :
779 1425 : unsigned int n_active_levels(const MeshBase & mesh)
780 : {
781 42 : libmesh_parallel_only(mesh.comm());
782 :
783 1425 : unsigned int nl = n_active_local_levels(mesh);
784 :
785 2808 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
786 5616 : mesh.unpartitioned_elements_end()))
787 0 : if (elem->active())
788 1341 : nl = std::max(elem->level() + 1, nl);
789 :
790 1425 : mesh.comm().max(nl);
791 1425 : return nl;
792 : }
793 :
794 :
795 :
796 1652430 : unsigned int n_local_levels(const MeshBase & mesh)
797 : {
798 1652430 : unsigned int nl = 0;
799 :
800 4874238 : for (const auto & elem : as_range(mesh.local_elements_begin(),
801 36514058 : mesh.local_elements_end()))
802 33438006 : nl = std::max(elem->level() + 1, nl);
803 :
804 1652430 : return nl;
805 : }
806 :
807 :
808 :
809 1652430 : unsigned int n_levels(const MeshBase & mesh)
810 : {
811 28338 : libmesh_parallel_only(mesh.comm());
812 :
813 1652430 : unsigned int nl = n_local_levels(mesh);
814 :
815 3669858 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
816 20088704 : mesh.unpartitioned_elements_end()))
817 15131430 : nl = std::max(elem->level() + 1, nl);
818 :
819 1652430 : mesh.comm().max(nl);
820 :
821 : // n_levels() is only valid and should only be called in cases where
822 : // the mesh is validly distributed (or serialized). Let's run an
823 : // expensive test in debug mode to make sure this is such a case.
824 : #ifdef DEBUG
825 28338 : const unsigned int paranoid_nl = paranoid_n_levels(mesh);
826 28338 : libmesh_assert_equal_to(nl, paranoid_nl);
827 : #endif
828 1652430 : return nl;
829 : }
830 :
831 :
832 :
833 38446 : unsigned int paranoid_n_levels(const MeshBase & mesh)
834 : {
835 28636 : libmesh_parallel_only(mesh.comm());
836 :
837 38446 : unsigned int nl = 0;
838 4051009 : for (const auto & elem : mesh.element_ptr_range())
839 4012265 : nl = std::max(elem->level() + 1, nl);
840 :
841 38446 : mesh.comm().max(nl);
842 38446 : return nl;
843 : }
844 :
845 :
846 :
847 639 : dof_id_type n_connected_components(const MeshBase & mesh,
848 : Real constraint_tol)
849 : {
850 36 : LOG_SCOPE("n_connected_components()", "MeshTools");
851 :
852 : // Yes, I'm being lazy. This is for mesh analysis before a
853 : // simulation, not anything going in any loops.
854 639 : if (!mesh.is_serial_on_zero())
855 0 : libmesh_not_implemented();
856 :
857 639 : dof_id_type n_components = 0;
858 :
859 657 : if (mesh.processor_id())
860 : {
861 522 : mesh.comm().broadcast(n_components);
862 531 : return n_components;
863 : }
864 :
865 : // All nodes in a set here are connected (at least indirectly) to
866 : // all other nodes in the same set, but have not yet been discovered
867 : // to be connected to nodes in other sets.
868 9 : std::vector<std::unordered_set<const Node *>> components;
869 :
870 : // With a typical mesh with few components and somewhat-contiguous
871 : // ordering, vector performance should be fine. With a mesh with
872 : // many components or completely scrambled ordering, performance
873 : // can be a disaster.
874 3204 : auto find_component = [&components](const Node * n) {
875 267 : std::unordered_set<const Node *> * component = nullptr;
876 :
877 8132 : for (auto & c: components)
878 4928 : if (c.find(n) != c.end())
879 : {
880 105 : libmesh_assert(component == nullptr);
881 105 : component = &c;
882 : }
883 :
884 3204 : return component;
885 108 : };
886 :
887 : auto add_to_component =
888 2010 : [&find_component]
889 2412 : (std::unordered_set<const Node *> & component, const Node * n) {
890 :
891 2412 : auto current_component = find_component(n);
892 : // We may already know we're in the desired component
893 2412 : if (&component == current_component)
894 51 : return;
895 :
896 : // If we're unknown, we should be in the desired component
897 1950 : else if (!current_component)
898 147 : component.insert(n);
899 :
900 : // If we think we're in another component, it should actually be
901 : // part of the desired component
902 : else
903 : {
904 3 : component.merge(*current_component);
905 3 : libmesh_assert(current_component->empty());
906 : }
907 99 : };
908 :
909 9 : auto & constraint_rows = mesh.get_constraint_rows();
910 :
911 1659 : for (const auto & elem : mesh.element_ptr_range())
912 : {
913 792 : const Node * first_node = elem->node_ptr(0);
914 :
915 792 : auto component = find_component(first_node);
916 :
917 : // If we didn't find one, make a new one, reusing an existing
918 : // slot if possible or growing our vector if necessary
919 792 : if (!component)
920 684 : for (auto & c: components)
921 354 : if (c.empty())
922 0 : component = &c;
923 :
924 432 : if (!component)
925 219 : component = &components.emplace_back();
926 :
927 3234 : for (const Node & n : elem->node_ref_range())
928 : {
929 2376 : add_to_component(*component, &n);
930 :
931 2376 : auto it = constraint_rows.find(&n);
932 2376 : if (it == constraint_rows.end())
933 195 : continue;
934 :
935 72 : for (const auto & [pr, val] : it->second)
936 : {
937 : // Ignore too-trivial constraint coefficients if
938 : // we get a non-default-0 constraint_tol
939 39 : if (std::abs(val) < constraint_tol)
940 0 : continue;
941 :
942 36 : const Elem * spline_elem = pr.first;
943 3 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
944 :
945 : const Node * spline_node =
946 36 : spline_elem->node_ptr(pr.second);
947 :
948 36 : add_to_component(*component, spline_node);
949 : }
950 : }
951 90 : }
952 :
953 327 : for (auto & component : components)
954 219 : if (!component.empty())
955 144 : ++n_components;
956 :
957 : // We calculated this on proc 0; now let everyone else know too
958 108 : mesh.comm().broadcast(n_components);
959 :
960 108 : return n_components;
961 90 : }
962 :
963 :
964 :
965 0 : void get_not_subactive_node_ids(const MeshBase & mesh,
966 : std::set<dof_id_type> & not_subactive_node_ids)
967 : {
968 0 : for (const auto & elem : mesh.element_ptr_range())
969 0 : if (!elem->subactive())
970 0 : for (auto & n : elem->node_ref_range())
971 0 : not_subactive_node_ids.insert(n.id());
972 0 : }
973 :
974 :
975 :
976 485036 : dof_id_type n_elem (const MeshBase::const_element_iterator & begin,
977 : const MeshBase::const_element_iterator & end)
978 : {
979 950552 : return cast_int<dof_id_type>(std::distance(begin, end));
980 : }
981 :
982 :
983 :
984 0 : dof_id_type n_nodes (const MeshBase::const_node_iterator & begin,
985 : const MeshBase::const_node_iterator & end)
986 : {
987 0 : return cast_int<dof_id_type>(std::distance(begin, end));
988 : }
989 :
990 :
991 :
992 426 : Real volume (const MeshBase & mesh,
993 : unsigned int dim)
994 : {
995 12 : libmesh_parallel_only(mesh.comm());
996 :
997 426 : if (dim == libMesh::invalid_uint)
998 426 : dim = mesh.mesh_dimension();
999 :
1000 426 : Real vol = 0;
1001 :
1002 : // first my local elements
1003 889 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1004 2807 : mesh.local_elements_end()))
1005 588 : if (elem->dim() == dim)
1006 990 : vol += elem->volume();
1007 :
1008 : // then count any unpartitioned objects, once
1009 438 : if (mesh.processor_id() == 0)
1010 138 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1011 276 : mesh.unpartitioned_elements_end()))
1012 0 : if (elem->dim() == dim)
1013 60 : vol += elem->volume();
1014 :
1015 426 : mesh.comm().sum(vol);
1016 426 : return vol;
1017 : }
1018 :
1019 :
1020 :
1021 1425 : unsigned int n_p_levels (const MeshBase & mesh)
1022 : {
1023 42 : libmesh_parallel_only(mesh.comm());
1024 :
1025 1425 : unsigned int max_p_level = 0;
1026 :
1027 : // first my local elements
1028 7068 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1029 48788 : mesh.local_elements_end()))
1030 48773 : max_p_level = std::max(elem->p_level(), max_p_level);
1031 :
1032 : // then any unpartitioned objects
1033 2808 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1034 5616 : mesh.unpartitioned_elements_end()))
1035 1341 : max_p_level = std::max(elem->p_level(), max_p_level);
1036 :
1037 1425 : mesh.comm().max(max_p_level);
1038 1425 : return max_p_level + 1;
1039 : }
1040 :
1041 :
1042 :
1043 0 : void find_nodal_neighbors(const MeshBase &,
1044 : const Node & node,
1045 : const std::vector<std::vector<const Elem *>> & nodes_to_elem_map,
1046 : std::vector<const Node *> & neighbors)
1047 : {
1048 0 : find_nodal_neighbors_helper(node.id(), nodes_to_elem_map[node.id()],
1049 : neighbors);
1050 0 : }
1051 :
1052 :
1053 :
1054 291668 : void find_nodal_neighbors(const MeshBase &,
1055 : const Node & node,
1056 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1057 : std::vector<const Node *> & neighbors)
1058 : {
1059 : const std::vector<const Elem *> node_to_elem_vec =
1060 299884 : libmesh_map_find(nodes_to_elem_map, node.id());
1061 291668 : find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
1062 291668 : }
1063 :
1064 :
1065 :
1066 0 : void find_hanging_nodes_and_parents(const MeshBase & mesh,
1067 : std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes)
1068 : {
1069 : // Loop through all the elements
1070 0 : for (auto & elem : mesh.active_local_element_ptr_range())
1071 0 : if (elem->type() == QUAD4)
1072 0 : for (auto s : elem->side_index_range())
1073 : {
1074 : // Loop over the sides looking for sides that have hanging nodes
1075 : // This code is inspired by compute_proj_constraints()
1076 0 : const Elem * neigh = elem->neighbor_ptr(s);
1077 :
1078 : // If not a boundary side
1079 0 : if (neigh != nullptr)
1080 : {
1081 : // Is there a coarser element next to this one?
1082 0 : if (neigh->level() < elem->level())
1083 : {
1084 0 : const Elem * ancestor = elem;
1085 0 : while (neigh->level() < ancestor->level())
1086 0 : ancestor = ancestor->parent();
1087 0 : unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
1088 0 : libmesh_assert_less (s_neigh, neigh->n_neighbors());
1089 :
1090 : // Couple of helper uints...
1091 0 : unsigned int local_node1=0;
1092 0 : unsigned int local_node2=0;
1093 :
1094 0 : bool found_in_neighbor = false;
1095 :
1096 : // Find the two vertices that make up this side
1097 0 : while (!elem->is_node_on_side(local_node1++,s)) { }
1098 0 : local_node1--;
1099 :
1100 : // Start looking for the second one with the next node
1101 0 : local_node2=local_node1+1;
1102 :
1103 : // Find the other one
1104 0 : while (!elem->is_node_on_side(local_node2++,s)) { }
1105 0 : local_node2--;
1106 :
1107 : //Pull out their global ids:
1108 0 : dof_id_type node1 = elem->node_id(local_node1);
1109 0 : dof_id_type node2 = elem->node_id(local_node2);
1110 :
1111 : // Now find which node is present in the neighbor
1112 : // FIXME This assumes a level one rule!
1113 : // The _other_ one is the hanging node
1114 :
1115 : // First look for the first one
1116 : // FIXME could be streamlined a bit
1117 0 : for (unsigned int n=0;n<neigh->n_sides();n++)
1118 0 : if (neigh->node_id(n) == node1)
1119 0 : found_in_neighbor=true;
1120 :
1121 0 : dof_id_type hanging_node=0;
1122 :
1123 0 : if (!found_in_neighbor)
1124 0 : hanging_node=node1;
1125 : else // If it wasn't node1 then it must be node2!
1126 0 : hanging_node=node2;
1127 :
1128 : // Reset for reuse
1129 0 : local_node1=0;
1130 :
1131 : // Find the first node that makes up the side in the neighbor (these should be the parent nodes)
1132 0 : while (!neigh->is_node_on_side(local_node1++,s_neigh)) { }
1133 0 : local_node1--;
1134 :
1135 0 : local_node2=local_node1+1;
1136 :
1137 : // Find the second node...
1138 0 : while (!neigh->is_node_on_side(local_node2++,s_neigh)) { }
1139 0 : local_node2--;
1140 :
1141 : // Save them if we haven't already found the parents for this one
1142 0 : if (hanging_nodes[hanging_node].size()<2)
1143 : {
1144 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node1));
1145 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node2));
1146 : }
1147 : }
1148 : }
1149 0 : }
1150 0 : }
1151 :
1152 :
1153 :
1154 36 : void clear_spline_nodes(MeshBase & mesh)
1155 : {
1156 6 : std::vector<Elem *> nodeelem_to_delete;
1157 :
1158 8737 : for (auto & elem : mesh.element_ptr_range())
1159 5073 : if (elem->type() == NODEELEM &&
1160 4140 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
1161 4170 : nodeelem_to_delete.push_back(elem);
1162 :
1163 3 : auto & constraint_rows = mesh.get_constraint_rows();
1164 :
1165 : // All our constraint_rows ought to be for spline constraints we're
1166 : // about to get rid of.
1167 : #ifndef NDEBUG
1168 762 : for (auto & node_row : constraint_rows)
1169 2064 : for (auto pr : node_row.second)
1170 : {
1171 1305 : const Elem * elem = pr.first.first;
1172 1305 : libmesh_assert(elem->type() == NODEELEM);
1173 1305 : libmesh_assert(elem->mapping_type() == RATIONAL_BERNSTEIN_MAP);
1174 : }
1175 : #endif
1176 :
1177 3 : constraint_rows.clear();
1178 :
1179 4176 : for (Elem * elem : nodeelem_to_delete)
1180 : {
1181 690 : Node * node = elem->node_ptr(0);
1182 4140 : mesh.delete_elem(elem);
1183 4140 : mesh.delete_node(node);
1184 : }
1185 36 : }
1186 :
1187 :
1188 :
1189 936 : bool valid_is_prepared (const MeshBase & mesh)
1190 : {
1191 84 : LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools");
1192 :
1193 936 : if (!mesh.is_prepared())
1194 0 : return true;
1195 :
1196 978 : std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
1197 :
1198 : // Try preparing (without allowing repartitioning or renumbering, to
1199 : // avoid false assertion failures)
1200 68 : bool old_allow_renumbering = mesh_clone->allow_renumbering();
1201 42 : mesh_clone->allow_renumbering(false);
1202 : bool old_allow_remote_element_removal =
1203 68 : mesh_clone->allow_remote_element_removal();
1204 68 : bool old_skip_partitioning = mesh_clone->skip_partitioning();
1205 42 : mesh_clone->skip_partitioning(true);
1206 42 : mesh_clone->allow_remote_element_removal(false);
1207 936 : mesh_clone->prepare_for_use();
1208 42 : mesh_clone->allow_renumbering(old_allow_renumbering);
1209 42 : mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal);
1210 42 : mesh_clone->skip_partitioning(old_skip_partitioning);
1211 :
1212 936 : return (mesh == *mesh_clone);
1213 868 : }
1214 :
1215 :
1216 :
1217 : #ifndef NDEBUG
1218 :
1219 122 : void libmesh_assert_equal_n_systems (const MeshBase & mesh)
1220 : {
1221 244 : LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
1222 :
1223 122 : unsigned int n_sys = libMesh::invalid_uint;
1224 :
1225 10842 : for (const auto & elem : mesh.element_ptr_range())
1226 : {
1227 10720 : if (n_sys == libMesh::invalid_uint)
1228 122 : n_sys = elem->n_systems();
1229 : else
1230 10598 : libmesh_assert_equal_to (elem->n_systems(), n_sys);
1231 : }
1232 :
1233 29237 : for (const auto & node : mesh.node_ptr_range())
1234 : {
1235 29115 : if (n_sys == libMesh::invalid_uint)
1236 0 : n_sys = node->n_systems();
1237 : else
1238 29115 : libmesh_assert_equal_to (node->n_systems(), n_sys);
1239 : }
1240 122 : }
1241 :
1242 :
1243 :
1244 : #ifdef LIBMESH_ENABLE_AMR
1245 0 : void libmesh_assert_old_dof_objects (const MeshBase & mesh)
1246 : {
1247 0 : LOG_SCOPE("libmesh_assert_old_dof_objects()", "MeshTools");
1248 :
1249 0 : for (const auto & elem : mesh.element_ptr_range())
1250 : {
1251 0 : if (elem->refinement_flag() == Elem::JUST_REFINED ||
1252 0 : elem->refinement_flag() == Elem::INACTIVE)
1253 0 : continue;
1254 :
1255 0 : if (elem->has_dofs())
1256 0 : libmesh_assert(elem->get_old_dof_object());
1257 :
1258 0 : for (auto & node : elem->node_ref_range())
1259 0 : if (node.has_dofs())
1260 0 : libmesh_assert(node.get_old_dof_object());
1261 : }
1262 0 : }
1263 : #else
1264 : void libmesh_assert_old_dof_objects (const MeshBase &) {}
1265 : #endif // LIBMESH_ENABLE_AMR
1266 :
1267 :
1268 :
1269 0 : void libmesh_assert_valid_node_pointers(const MeshBase & mesh)
1270 : {
1271 0 : LOG_SCOPE("libmesh_assert_valid_node_pointers()", "MeshTools");
1272 :
1273 : // Here we specifically do not want "auto &" because we need to
1274 : // reseat the (temporary) pointer variable in the loop below,
1275 : // without modifying the original.
1276 0 : for (const Elem * elem : mesh.element_ptr_range())
1277 : {
1278 0 : libmesh_assert (elem);
1279 0 : while (elem)
1280 : {
1281 0 : elem->libmesh_assert_valid_node_pointers();
1282 0 : for (auto n : elem->neighbor_ptr_range())
1283 0 : if (n && n != remote_elem)
1284 0 : n->libmesh_assert_valid_node_pointers();
1285 :
1286 0 : libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1287 0 : elem = elem->parent();
1288 : }
1289 : }
1290 0 : }
1291 :
1292 :
1293 :
1294 8698 : void libmesh_assert_valid_remote_elems(const MeshBase & mesh)
1295 : {
1296 17396 : LOG_SCOPE("libmesh_assert_valid_remote_elems()", "MeshTools");
1297 :
1298 645131 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1299 1290262 : mesh.local_elements_end()))
1300 : {
1301 636433 : libmesh_assert (elem);
1302 :
1303 : // We currently don't allow active_local_elements to have
1304 : // remote_elem neighbors
1305 636433 : if (elem->active())
1306 2782245 : for (auto n : elem->neighbor_ptr_range())
1307 2239277 : libmesh_assert_not_equal_to (n, remote_elem);
1308 :
1309 : #ifdef LIBMESH_ENABLE_AMR
1310 636433 : const Elem * parent = elem->parent();
1311 636433 : if (parent)
1312 355748 : libmesh_assert_not_equal_to (parent, remote_elem);
1313 :
1314 : // We can only be strict about active elements' subactive
1315 : // children
1316 636433 : if (elem->active() && elem->has_children())
1317 26661 : for (auto & child : elem->child_ref_range())
1318 21348 : libmesh_assert_not_equal_to (&child, remote_elem);
1319 : #endif
1320 : }
1321 8698 : }
1322 :
1323 :
1324 :
1325 388 : void libmesh_assert_valid_elem_ids(const MeshBase & mesh)
1326 : {
1327 776 : LOG_SCOPE("libmesh_assert_valid_elem_ids()", "MeshTools");
1328 :
1329 388 : processor_id_type lastprocid = 0;
1330 388 : dof_id_type lastelemid = 0;
1331 :
1332 13378 : for (const auto & elem : mesh.active_element_ptr_range())
1333 : {
1334 12990 : libmesh_assert (elem);
1335 12990 : processor_id_type elemprocid = elem->processor_id();
1336 12990 : dof_id_type elemid = elem->id();
1337 :
1338 12990 : libmesh_assert_greater_equal (elemid, lastelemid);
1339 12990 : libmesh_assert_greater_equal (elemprocid, lastprocid);
1340 :
1341 12990 : lastelemid = elemid;
1342 12990 : lastprocid = elemprocid;
1343 : }
1344 388 : }
1345 :
1346 :
1347 :
1348 816 : void libmesh_assert_valid_amr_elem_ids(const MeshBase & mesh)
1349 : {
1350 1632 : LOG_SCOPE("libmesh_assert_valid_amr_elem_ids()", "MeshTools");
1351 :
1352 320464 : for (const auto & elem : mesh.element_ptr_range())
1353 : {
1354 319648 : libmesh_assert (elem);
1355 :
1356 319648 : const Elem * parent = elem->parent();
1357 :
1358 319648 : if (parent)
1359 : {
1360 24548 : libmesh_assert_greater_equal (elem->id(), parent->id());
1361 24548 : libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1362 : }
1363 : }
1364 816 : }
1365 :
1366 :
1367 :
1368 12256 : void libmesh_assert_valid_amr_interior_parents(const MeshBase & mesh)
1369 : {
1370 24512 : LOG_SCOPE("libmesh_assert_valid_amr_interior_parents()", "MeshTools");
1371 :
1372 1395355 : for (const auto & elem : mesh.element_ptr_range())
1373 : {
1374 1383099 : libmesh_assert (elem);
1375 :
1376 : // We can skip to the next element if we're full-dimension
1377 : // and therefore don't have any interior parents
1378 1383099 : if (elem->dim() >= LIBMESH_DIM)
1379 623182 : continue;
1380 :
1381 759917 : const Elem * ip = elem->interior_parent();
1382 :
1383 759917 : const Elem * parent = elem->parent();
1384 :
1385 759917 : if (ip && (ip != remote_elem) && parent)
1386 : {
1387 1000 : libmesh_assert_equal_to (ip->top_parent(),
1388 : elem->top_parent()->interior_parent());
1389 :
1390 1000 : if (ip->level() == elem->level())
1391 1000 : libmesh_assert_equal_to (ip->parent(),
1392 : parent->interior_parent());
1393 : else
1394 : {
1395 0 : libmesh_assert_less (ip->level(), elem->level());
1396 0 : libmesh_assert_equal_to (ip, parent->interior_parent());
1397 : }
1398 : }
1399 : }
1400 12256 : }
1401 :
1402 :
1403 :
1404 0 : void libmesh_assert_contiguous_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1405 : {
1406 0 : LOG_SCOPE("libmesh_assert_contiguous_dof_ids()", "MeshTools");
1407 :
1408 0 : if (mesh.n_processors() == 1)
1409 0 : return;
1410 :
1411 0 : libmesh_parallel_only(mesh.comm());
1412 :
1413 0 : dof_id_type min_dof_id = std::numeric_limits<dof_id_type>::max(),
1414 0 : max_dof_id = std::numeric_limits<dof_id_type>::min();
1415 :
1416 : // Figure out what our local dof id range is
1417 0 : for (const auto * node : mesh.local_node_ptr_range())
1418 : {
1419 0 : for (auto v : make_range(node->n_vars(sysnum)))
1420 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1421 : {
1422 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1423 0 : min_dof_id = std::min (min_dof_id, id);
1424 0 : max_dof_id = std::max (max_dof_id, id);
1425 : }
1426 : }
1427 :
1428 : // Make sure no other processors' ids are inside it
1429 0 : for (const auto * node : mesh.node_ptr_range())
1430 : {
1431 0 : if (node->processor_id() == mesh.processor_id())
1432 0 : continue;
1433 0 : for (auto v : make_range(node->n_vars(sysnum)))
1434 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1435 : {
1436 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1437 0 : libmesh_assert (id < min_dof_id ||
1438 : id > max_dof_id);
1439 : }
1440 : }
1441 : }
1442 :
1443 :
1444 :
1445 : template <>
1446 17694 : void libmesh_assert_topology_consistent_procids<Elem>(const MeshBase & mesh)
1447 : {
1448 35388 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1449 :
1450 : // This parameter is not used when !LIBMESH_ENABLE_AMR
1451 17694 : libmesh_ignore(mesh);
1452 :
1453 : // If we're adaptively refining, check processor ids for consistency
1454 : // between parents and children.
1455 : #ifdef LIBMESH_ENABLE_AMR
1456 :
1457 : // Ancestor elements we won't worry about, but subactive and active
1458 : // elements ought to have parents with consistent processor ids
1459 2642334 : for (const auto & elem : mesh.element_ptr_range())
1460 : {
1461 2624640 : libmesh_assert(elem);
1462 :
1463 2624640 : if (!elem->active() && !elem->subactive())
1464 288060 : continue;
1465 :
1466 2336580 : const Elem * parent = elem->parent();
1467 :
1468 2336580 : if (parent)
1469 : {
1470 1157608 : libmesh_assert(parent->has_children());
1471 1157608 : processor_id_type parent_procid = parent->processor_id();
1472 1157608 : bool matching_child_id = false;
1473 : // If we've got a remote_elem then we don't know whether
1474 : // it's responsible for the parent's processor id; all
1475 : // we can do is assume it is and let its processor fail
1476 : // an assert if there's something wrong.
1477 7222200 : for (auto & child : parent->child_ref_range())
1478 12129184 : if (&child == remote_elem ||
1479 6064592 : child.processor_id() == parent_procid)
1480 5909956 : matching_child_id = true;
1481 1157608 : libmesh_assert(matching_child_id);
1482 : }
1483 : }
1484 : #endif
1485 17694 : }
1486 :
1487 :
1488 :
1489 : template <>
1490 26966 : void libmesh_assert_topology_consistent_procids<Node>(const MeshBase & mesh)
1491 : {
1492 26966 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1493 :
1494 26966 : if (mesh.n_processors() == 1)
1495 0 : return;
1496 :
1497 26966 : libmesh_parallel_only(mesh.comm());
1498 :
1499 : // We want this test to be valid even when called after nodes have
1500 : // been added asynchronously but before they're renumbered.
1501 : //
1502 : // Plus, some code (looking at you, stitch_meshes) modifies
1503 : // DofObject ids without keeping max_elem_id()/max_node_id()
1504 : // consistent, but that's done in a safe way for performance
1505 : // reasons, so we'll play along and just figure out new max ids
1506 : // ourselves.
1507 26966 : dof_id_type parallel_max_node_id = 0;
1508 7279169 : for (const auto & node : mesh.node_ptr_range())
1509 7252203 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
1510 7252203 : node->id()+1);
1511 26966 : mesh.comm().max(parallel_max_node_id);
1512 :
1513 :
1514 53932 : std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1515 :
1516 2223802 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1517 4447604 : mesh.local_elements_end()))
1518 : {
1519 2196836 : libmesh_assert (elem);
1520 :
1521 16558620 : for (auto & node : elem->node_ref_range())
1522 : {
1523 14361784 : dof_id_type nodeid = node.id();
1524 14361784 : node_touched_by_me[nodeid] = true;
1525 : }
1526 : }
1527 53932 : std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1528 26966 : mesh.comm().max(node_touched_by_anyone);
1529 :
1530 3640996 : for (const auto & node : mesh.local_node_ptr_range())
1531 : {
1532 3614030 : libmesh_assert(node);
1533 3614030 : dof_id_type nodeid = node->id();
1534 3614030 : libmesh_assert(!node_touched_by_anyone[nodeid] ||
1535 : node_touched_by_me[nodeid]);
1536 : }
1537 : }
1538 :
1539 :
1540 :
1541 0 : void libmesh_assert_canonical_node_procids (const MeshBase & mesh)
1542 : {
1543 0 : for (const auto & elem : mesh.active_element_ptr_range())
1544 0 : for (auto & node : elem->node_ref_range())
1545 0 : libmesh_assert_equal_to
1546 : (node.processor_id(),
1547 : node.choose_processor_id(node.processor_id(),
1548 : elem->processor_id()));
1549 0 : }
1550 :
1551 :
1552 :
1553 : #ifdef LIBMESH_ENABLE_AMR
1554 1226 : void libmesh_assert_valid_refinement_tree(const MeshBase & mesh)
1555 : {
1556 2452 : LOG_SCOPE("libmesh_assert_valid_refinement_tree()", "MeshTools");
1557 :
1558 59722 : for (const auto & elem : mesh.element_ptr_range())
1559 : {
1560 58496 : libmesh_assert(elem);
1561 58496 : if (elem->has_children())
1562 17748 : for (auto & child : elem->child_ref_range())
1563 14928 : if (&child != remote_elem)
1564 14000 : libmesh_assert_equal_to (child.parent(), elem);
1565 58496 : if (elem->active())
1566 : {
1567 55676 : libmesh_assert(!elem->ancestor());
1568 55676 : libmesh_assert(!elem->subactive());
1569 : }
1570 2820 : else if (elem->ancestor())
1571 : {
1572 2820 : libmesh_assert(!elem->subactive());
1573 : }
1574 : else
1575 0 : libmesh_assert(elem->subactive());
1576 :
1577 58496 : if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1578 0 : libmesh_assert_greater(elem->p_level(), 0);
1579 : }
1580 1226 : }
1581 : #else
1582 : void libmesh_assert_valid_refinement_tree(const MeshBase &)
1583 : {
1584 : }
1585 : #endif // LIBMESH_ENABLE_AMR
1586 :
1587 : #endif // !NDEBUG
1588 :
1589 :
1590 :
1591 : #ifdef DEBUG
1592 :
1593 0 : void libmesh_assert_no_links_to_elem(const MeshBase & mesh,
1594 : const Elem * bad_elem)
1595 : {
1596 0 : for (const auto & elem : mesh.element_ptr_range())
1597 : {
1598 0 : libmesh_assert (elem);
1599 0 : libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1600 0 : for (auto n : elem->neighbor_ptr_range())
1601 0 : libmesh_assert_not_equal_to (n, bad_elem);
1602 :
1603 : #ifdef LIBMESH_ENABLE_AMR
1604 0 : if (elem->has_children())
1605 0 : for (auto & child : elem->child_ref_range())
1606 0 : libmesh_assert_not_equal_to (&child, bad_elem);
1607 : #endif
1608 : }
1609 0 : }
1610 :
1611 :
1612 566 : void libmesh_assert_equal_points (const MeshBase & mesh)
1613 : {
1614 1132 : LOG_SCOPE("libmesh_assert_equal_points()", "MeshTools");
1615 :
1616 566 : dof_id_type pmax_node_id = mesh.max_node_id();
1617 566 : mesh.comm().max(pmax_node_id);
1618 :
1619 915766 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1620 : {
1621 915200 : const Point * p = mesh.query_node_ptr(i);
1622 :
1623 915200 : libmesh_assert(mesh.comm().semiverify(p));
1624 : }
1625 566 : }
1626 :
1627 :
1628 566 : void libmesh_assert_equal_connectivity (const MeshBase & mesh)
1629 : {
1630 1132 : LOG_SCOPE("libmesh_assert_equal_connectivity()", "MeshTools");
1631 :
1632 566 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1633 566 : mesh.comm().max(pmax_elem_id);
1634 :
1635 1142058 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1636 : {
1637 1141492 : const Elem * e = mesh.query_elem_ptr(i);
1638 :
1639 2282984 : std::vector<dof_id_type> nodes;
1640 1141492 : if (e)
1641 5771288 : for (auto n : e->node_index_range())
1642 4629912 : nodes.push_back(e->node_id(n));
1643 :
1644 1141492 : libmesh_assert(mesh.comm().semiverify(e ? &nodes : nullptr));
1645 : }
1646 566 : }
1647 :
1648 :
1649 0 : void libmesh_assert_connected_nodes (const MeshBase & mesh)
1650 : {
1651 0 : LOG_SCOPE("libmesh_assert_connected_nodes()", "MeshTools");
1652 :
1653 0 : std::set<const Node *> used_nodes;
1654 :
1655 0 : for (const auto & elem : mesh.element_ptr_range())
1656 : {
1657 0 : libmesh_assert (elem);
1658 :
1659 0 : for (auto & n : elem->node_ref_range())
1660 0 : used_nodes.insert(&n);
1661 : }
1662 :
1663 0 : for (const auto & node : mesh.node_ptr_range())
1664 : {
1665 0 : libmesh_assert(node);
1666 0 : libmesh_assert(used_nodes.count(node));
1667 : }
1668 0 : }
1669 :
1670 :
1671 :
1672 13244 : void libmesh_assert_valid_constraint_rows (const MeshBase & mesh)
1673 : {
1674 13244 : libmesh_parallel_only(mesh.comm());
1675 :
1676 13244 : const auto & constraint_rows = mesh.get_constraint_rows();
1677 :
1678 13244 : bool have_constraint_rows = !constraint_rows.empty();
1679 13244 : mesh.comm().max(have_constraint_rows);
1680 13244 : if (!have_constraint_rows)
1681 13204 : return;
1682 :
1683 12700 : for (auto & row : constraint_rows)
1684 : {
1685 12660 : const Node * node = row.first;
1686 12660 : libmesh_assert(node == mesh.node_ptr(node->id()));
1687 :
1688 44126 : for (auto & pr : row.second)
1689 : {
1690 31466 : const Elem * spline_elem = pr.first.first;
1691 31466 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1692 : }
1693 : }
1694 :
1695 40 : dof_id_type pmax_node_id = mesh.max_node_id();
1696 40 : mesh.comm().max(pmax_node_id);
1697 :
1698 16990 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1699 : {
1700 16950 : const Node * node = mesh.query_node_ptr(i);
1701 :
1702 16950 : bool have_constraint = constraint_rows.count(node);
1703 :
1704 16950 : const std::size_t my_n_constraints = have_constraint ?
1705 12660 : libmesh_map_find(constraint_rows, node).size() : std::size_t(-1);
1706 16950 : const std::size_t * n_constraints = node ?
1707 16950 : &my_n_constraints : nullptr;
1708 :
1709 16950 : libmesh_assert(mesh.comm().semiverify(n_constraints));
1710 : }
1711 : }
1712 :
1713 :
1714 :
1715 33212 : void libmesh_assert_valid_boundary_ids(const MeshBase & mesh)
1716 : {
1717 33212 : LOG_SCOPE("libmesh_assert_valid_boundary_ids()", "MeshTools");
1718 :
1719 33212 : if (mesh.n_processors() == 1)
1720 688 : return;
1721 :
1722 32524 : libmesh_parallel_only(mesh.comm());
1723 :
1724 32524 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1725 :
1726 32524 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1727 32524 : mesh.comm().max(pmax_elem_id);
1728 :
1729 4421554 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1730 : {
1731 4389030 : const Elem * elem = mesh.query_elem_ptr(i);
1732 4389030 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1733 4389030 : const unsigned int my_n_edges = elem ? elem->n_edges() : 0;
1734 4389030 : const unsigned int my_n_sides = elem ? elem->n_sides() : 0;
1735 : unsigned int
1736 4389030 : n_nodes = my_n_nodes,
1737 4389030 : n_edges = my_n_edges,
1738 4389030 : n_sides = my_n_sides;
1739 :
1740 4389030 : mesh.comm().max(n_nodes);
1741 4389030 : mesh.comm().max(n_edges);
1742 4389030 : mesh.comm().max(n_sides);
1743 :
1744 4389030 : if (elem)
1745 : {
1746 4311263 : libmesh_assert_equal_to(my_n_nodes, n_nodes);
1747 4311263 : libmesh_assert_equal_to(my_n_edges, n_edges);
1748 4311263 : libmesh_assert_equal_to(my_n_sides, n_sides);
1749 : }
1750 :
1751 : // Let's test all IDs on the element with one communication
1752 : // rather than n_nodes + n_edges + n_sides communications, to
1753 : // cut down on latency in dbg modes.
1754 8778060 : std::vector<boundary_id_type> all_bcids;
1755 :
1756 33650182 : for (unsigned int n=0; n != n_nodes; ++n)
1757 : {
1758 58522304 : std::vector<boundary_id_type> bcids;
1759 29261152 : if (elem)
1760 : {
1761 29097358 : boundary_info.boundary_ids(elem->node_ptr(n), bcids);
1762 :
1763 : // Ordering of boundary ids shouldn't matter
1764 29097358 : std::sort(bcids.begin(), bcids.end());
1765 : }
1766 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1767 :
1768 29261152 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1769 58522304 : bcids.end());
1770 : // Separator
1771 29261152 : all_bcids.push_back(BoundaryInfo::invalid_id);
1772 : }
1773 :
1774 27801122 : for (unsigned short e=0; e != n_edges; ++e)
1775 : {
1776 46824184 : std::vector<boundary_id_type> bcids;
1777 :
1778 23412092 : if (elem)
1779 : {
1780 23310828 : boundary_info.edge_boundary_ids(elem, e, bcids);
1781 :
1782 : // Ordering of boundary ids shouldn't matter
1783 23310828 : std::sort(bcids.begin(), bcids.end());
1784 : }
1785 :
1786 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1787 :
1788 23412092 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1789 46824184 : bcids.end());
1790 : // Separator
1791 23412092 : all_bcids.push_back(BoundaryInfo::invalid_id);
1792 :
1793 23412092 : if (elem)
1794 : {
1795 23310828 : boundary_info.raw_edge_boundary_ids(elem, e, bcids);
1796 :
1797 : // Ordering of boundary ids shouldn't matter
1798 23310828 : std::sort(bcids.begin(), bcids.end());
1799 :
1800 23310828 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1801 46621656 : bcids.end());
1802 : // Separator
1803 23310828 : all_bcids.push_back(BoundaryInfo::invalid_id);
1804 : }
1805 :
1806 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1807 : }
1808 :
1809 21801870 : for (unsigned short s=0; s != n_sides; ++s)
1810 : {
1811 34825680 : std::vector<boundary_id_type> bcids;
1812 :
1813 17412840 : if (elem)
1814 : {
1815 17338884 : boundary_info.boundary_ids(elem, s, bcids);
1816 :
1817 : // Ordering of boundary ids shouldn't matter
1818 17338884 : std::sort(bcids.begin(), bcids.end());
1819 :
1820 17338884 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1821 34677768 : bcids.end());
1822 : // Separator
1823 17338884 : all_bcids.push_back(BoundaryInfo::invalid_id);
1824 : }
1825 :
1826 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1827 :
1828 17412840 : if (elem)
1829 : {
1830 17338884 : boundary_info.raw_boundary_ids(elem, s, bcids);
1831 :
1832 : // Ordering of boundary ids shouldn't matter
1833 17338884 : std::sort(bcids.begin(), bcids.end());
1834 :
1835 17338884 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1836 34677768 : bcids.end());
1837 : // Separator
1838 17338884 : all_bcids.push_back(BoundaryInfo::invalid_id);
1839 : }
1840 :
1841 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1842 : }
1843 :
1844 13167090 : for (unsigned short sf=0; sf != 2; ++sf)
1845 : {
1846 17556120 : std::vector<boundary_id_type> bcids;
1847 :
1848 8778060 : if (elem)
1849 : {
1850 8622526 : boundary_info.shellface_boundary_ids(elem, sf, bcids);
1851 :
1852 : // Ordering of boundary ids shouldn't matter
1853 8622526 : std::sort(bcids.begin(), bcids.end());
1854 :
1855 8622526 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1856 17245052 : bcids.end());
1857 : // Separator
1858 8622526 : all_bcids.push_back(BoundaryInfo::invalid_id);
1859 : }
1860 :
1861 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1862 :
1863 8778060 : if (elem)
1864 : {
1865 8622526 : boundary_info.raw_shellface_boundary_ids(elem, sf, bcids);
1866 :
1867 : // Ordering of boundary ids shouldn't matter
1868 8622526 : std::sort(bcids.begin(), bcids.end());
1869 :
1870 8622526 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1871 17245052 : bcids.end());
1872 : // Separator
1873 8622526 : all_bcids.push_back(BoundaryInfo::invalid_id);
1874 : }
1875 :
1876 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1877 : }
1878 :
1879 4389030 : libmesh_assert(mesh.comm().semiverify
1880 : (elem ? &all_bcids : nullptr));
1881 : }
1882 : }
1883 :
1884 :
1885 7652 : void libmesh_assert_valid_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1886 : {
1887 7652 : LOG_SCOPE("libmesh_assert_valid_dof_ids()", "MeshTools");
1888 :
1889 7652 : if (mesh.n_processors() == 1)
1890 0 : return;
1891 :
1892 7652 : libmesh_parallel_only(mesh.comm());
1893 :
1894 7652 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1895 7652 : mesh.comm().max(pmax_elem_id);
1896 :
1897 973224 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1898 965572 : assert_semiverify_dofobj(mesh.comm(),
1899 965572 : mesh.query_elem_ptr(i),
1900 : sysnum);
1901 :
1902 7652 : dof_id_type pmax_node_id = mesh.max_node_id();
1903 7652 : mesh.comm().max(pmax_node_id);
1904 :
1905 1828510 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1906 1820858 : assert_semiverify_dofobj(mesh.comm(),
1907 1820858 : mesh.query_node_ptr(i),
1908 : sysnum);
1909 : }
1910 :
1911 :
1912 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1913 25512 : void libmesh_assert_valid_unique_ids(const MeshBase & mesh)
1914 : {
1915 51024 : LOG_SCOPE("libmesh_assert_valid_unique_ids()", "MeshTools");
1916 :
1917 25512 : libmesh_parallel_only(mesh.comm());
1918 :
1919 : // First collect all the unique_ids we can see and make sure there's
1920 : // no duplicates
1921 51024 : std::unordered_set<unique_id_type> semilocal_unique_ids;
1922 :
1923 2914753 : for (auto const & elem : mesh.active_element_ptr_range())
1924 : {
1925 2889241 : libmesh_assert (!semilocal_unique_ids.count(elem->unique_id()));
1926 2889241 : semilocal_unique_ids.insert(elem->unique_id());
1927 : }
1928 :
1929 5484895 : for (auto const & node : mesh.node_ptr_range())
1930 : {
1931 5459383 : libmesh_assert (!semilocal_unique_ids.count(node->unique_id()));
1932 5459383 : semilocal_unique_ids.insert(node->unique_id());
1933 : }
1934 :
1935 : // Then make sure elements are all in sync and remote elements don't
1936 : // duplicate semilocal
1937 :
1938 25512 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1939 25512 : mesh.comm().max(pmax_elem_id);
1940 :
1941 3457324 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1942 : {
1943 3431812 : const Elem * elem = mesh.query_elem_ptr(i);
1944 3431812 : assert_dofobj_unique_id(mesh.comm(), elem, semilocal_unique_ids);
1945 : }
1946 :
1947 : // Then make sure nodes are all in sync and remote elements don't
1948 : // duplicate semilocal
1949 :
1950 25512 : dof_id_type pmax_node_id = mesh.max_node_id();
1951 25512 : mesh.comm().max(pmax_node_id);
1952 :
1953 5791154 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1954 : {
1955 5765642 : const Node * node = mesh.query_node_ptr(i);
1956 5765642 : assert_dofobj_unique_id(mesh.comm(), node, semilocal_unique_ids);
1957 : }
1958 25512 : }
1959 : #endif
1960 :
1961 0 : void libmesh_assert_consistent_distributed(const MeshBase & mesh)
1962 : {
1963 0 : libmesh_parallel_only(mesh.comm());
1964 :
1965 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
1966 0 : mesh.comm().max(parallel_max_elem_id);
1967 :
1968 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
1969 : {
1970 0 : const Elem * elem = mesh.query_elem_ptr(i);
1971 : processor_id_type pid =
1972 0 : elem ? elem->processor_id() : DofObject::invalid_processor_id;
1973 0 : mesh.comm().min(pid);
1974 0 : libmesh_assert(elem || pid != mesh.processor_id());
1975 : }
1976 :
1977 0 : dof_id_type parallel_max_node_id = mesh.max_node_id();
1978 0 : mesh.comm().max(parallel_max_node_id);
1979 :
1980 0 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
1981 : {
1982 0 : const Node * node = mesh.query_node_ptr(i);
1983 : processor_id_type pid =
1984 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
1985 0 : mesh.comm().min(pid);
1986 0 : libmesh_assert(node || pid != mesh.processor_id());
1987 : }
1988 0 : }
1989 :
1990 :
1991 0 : void libmesh_assert_consistent_distributed_nodes(const MeshBase & mesh)
1992 : {
1993 0 : libmesh_parallel_only(mesh.comm());
1994 0 : auto locator = mesh.sub_point_locator();
1995 :
1996 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
1997 0 : mesh.comm().max(parallel_max_elem_id);
1998 :
1999 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2000 : {
2001 0 : const Elem * elem = mesh.query_elem_ptr(i);
2002 :
2003 0 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2004 0 : unsigned int n_nodes = my_n_nodes;
2005 0 : mesh.comm().max(n_nodes);
2006 :
2007 0 : if (n_nodes)
2008 0 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2009 :
2010 0 : for (unsigned int n=0; n != n_nodes; ++n)
2011 : {
2012 0 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2013 : processor_id_type pid =
2014 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2015 0 : mesh.comm().min(pid);
2016 0 : libmesh_assert(node || pid != mesh.processor_id());
2017 : }
2018 : }
2019 0 : }
2020 :
2021 :
2022 :
2023 : template <>
2024 17694 : void libmesh_assert_parallel_consistent_procids<Elem>(const MeshBase & mesh)
2025 : {
2026 17694 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2027 :
2028 17694 : if (mesh.n_processors() == 1)
2029 0 : return;
2030 :
2031 17694 : libmesh_parallel_only(mesh.comm());
2032 :
2033 : // Some code (looking at you, stitch_meshes) modifies DofObject ids
2034 : // without keeping max_elem_id()/max_node_id() consistent, but
2035 : // that's done in a safe way for performance reasons, so we'll play
2036 : // along and just figure out new max ids ourselves.
2037 17694 : dof_id_type parallel_max_elem_id = 0;
2038 2642334 : for (const auto & elem : mesh.element_ptr_range())
2039 2624640 : parallel_max_elem_id = std::max<dof_id_type>(parallel_max_elem_id,
2040 2624640 : elem->id()+1);
2041 17694 : mesh.comm().max(parallel_max_elem_id);
2042 :
2043 : // Check processor ids for consistency between processors
2044 :
2045 2711048 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2046 : {
2047 2693354 : const Elem * elem = mesh.query_elem_ptr(i);
2048 :
2049 : processor_id_type min_id =
2050 2624640 : elem ? elem->processor_id() :
2051 2693354 : std::numeric_limits<processor_id_type>::max();
2052 2693354 : mesh.comm().min(min_id);
2053 :
2054 : processor_id_type max_id =
2055 2624640 : elem ? elem->processor_id() :
2056 2693354 : std::numeric_limits<processor_id_type>::min();
2057 2693354 : mesh.comm().max(max_id);
2058 :
2059 2693354 : if (elem)
2060 : {
2061 2624640 : libmesh_assert_equal_to (min_id, elem->processor_id());
2062 2624640 : libmesh_assert_equal_to (max_id, elem->processor_id());
2063 : }
2064 :
2065 2693354 : if (min_id == mesh.processor_id())
2066 1272943 : libmesh_assert(elem);
2067 : }
2068 : }
2069 :
2070 :
2071 :
2072 76 : void libmesh_assert_parallel_consistent_new_node_procids(const MeshBase & mesh)
2073 : {
2074 76 : LOG_SCOPE("libmesh_assert_parallel_consistent_new_node_procids()", "MeshTools");
2075 :
2076 76 : if (mesh.n_processors() == 1)
2077 0 : return;
2078 :
2079 76 : libmesh_parallel_only(mesh.comm());
2080 :
2081 : // We want this test to hit every node when called even after nodes
2082 : // have been added asynchronously but before everything has been
2083 : // renumbered.
2084 76 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2085 76 : mesh.comm().max(parallel_max_elem_id);
2086 :
2087 152 : std::vector<bool> elem_touched_by_anyone(parallel_max_elem_id, false);
2088 :
2089 43248 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2090 : {
2091 43172 : const Elem * elem = mesh.query_elem_ptr(i);
2092 :
2093 43172 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2094 43172 : unsigned int n_nodes = my_n_nodes;
2095 43172 : mesh.comm().max(n_nodes);
2096 :
2097 43172 : if (n_nodes)
2098 15840 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2099 :
2100 173404 : for (unsigned int n=0; n != n_nodes; ++n)
2101 : {
2102 130232 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2103 130232 : const processor_id_type pid = node ? node->processor_id() : 0;
2104 130232 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2105 : }
2106 : }
2107 : }
2108 :
2109 : template <>
2110 27188 : void libmesh_assert_parallel_consistent_procids<Node>(const MeshBase & mesh)
2111 : {
2112 27188 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2113 :
2114 27188 : if (mesh.n_processors() == 1)
2115 0 : return;
2116 :
2117 27188 : libmesh_parallel_only(mesh.comm());
2118 :
2119 : // We want this test to be valid even when called even after nodes
2120 : // have been added asynchronously but before they're renumbered
2121 : //
2122 : // Plus, some code (looking at you, stitch_meshes) modifies
2123 : // DofObject ids without keeping max_elem_id()/max_node_id()
2124 : // consistent, but that's done in a safe way for performance
2125 : // reasons, so we'll play along and just figure out new max ids
2126 : // ourselves.
2127 27188 : dof_id_type parallel_max_node_id = 0;
2128 7529771 : for (const auto & node : mesh.node_ptr_range())
2129 7502583 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
2130 7502583 : node->id()+1);
2131 27188 : mesh.comm().max(parallel_max_node_id);
2132 :
2133 54376 : std::vector<bool> node_touched_by_anyone(parallel_max_node_id, false);
2134 :
2135 2269975 : for (const auto & elem : as_range(mesh.local_elements_begin(),
2136 4539950 : mesh.local_elements_end()))
2137 : {
2138 2242787 : libmesh_assert (elem);
2139 :
2140 17008104 : for (auto & node : elem->node_ref_range())
2141 : {
2142 14765317 : dof_id_type nodeid = node.id();
2143 14765317 : node_touched_by_anyone[nodeid] = true;
2144 : }
2145 : }
2146 27188 : mesh.comm().max(node_touched_by_anyone);
2147 :
2148 : // Check processor ids for consistency between processors
2149 : // on any node an element touches
2150 7949022 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2151 : {
2152 7921834 : if (!node_touched_by_anyone[i])
2153 421536 : continue;
2154 :
2155 7500298 : const Node * node = mesh.query_node_ptr(i);
2156 7500298 : const processor_id_type pid = node ? node->processor_id() : 0;
2157 :
2158 7500298 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2159 : }
2160 : }
2161 :
2162 :
2163 :
2164 : #ifdef LIBMESH_ENABLE_AMR
2165 0 : void libmesh_assert_valid_refinement_flags(const MeshBase & mesh)
2166 : {
2167 0 : LOG_SCOPE("libmesh_assert_valid_refinement_flags()", "MeshTools");
2168 :
2169 0 : libmesh_parallel_only(mesh.comm());
2170 0 : if (mesh.n_processors() == 1)
2171 0 : return;
2172 :
2173 0 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2174 0 : mesh.comm().max(pmax_elem_id);
2175 :
2176 0 : std::vector<unsigned char> my_elem_h_state(pmax_elem_id, 255);
2177 0 : std::vector<unsigned char> my_elem_p_state(pmax_elem_id, 255);
2178 :
2179 0 : for (const auto & elem : mesh.element_ptr_range())
2180 : {
2181 0 : libmesh_assert (elem);
2182 0 : dof_id_type elemid = elem->id();
2183 :
2184 0 : my_elem_h_state[elemid] =
2185 0 : static_cast<unsigned char>(elem->refinement_flag());
2186 :
2187 0 : my_elem_p_state[elemid] =
2188 0 : static_cast<unsigned char>(elem->p_refinement_flag());
2189 : }
2190 0 : std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
2191 0 : mesh.comm().min(min_elem_h_state);
2192 :
2193 0 : std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
2194 0 : mesh.comm().min(min_elem_p_state);
2195 :
2196 0 : for (dof_id_type i=0; i!= pmax_elem_id; ++i)
2197 : {
2198 0 : libmesh_assert(my_elem_h_state[i] == 255 ||
2199 : my_elem_h_state[i] == min_elem_h_state[i]);
2200 0 : libmesh_assert(my_elem_p_state[i] == 255 ||
2201 : my_elem_p_state[i] == min_elem_p_state[i]);
2202 : }
2203 : }
2204 : #else
2205 : void libmesh_assert_valid_refinement_flags(const MeshBase &)
2206 : {
2207 : }
2208 : #endif // LIBMESH_ENABLE_AMR
2209 :
2210 :
2211 :
2212 15162 : void libmesh_assert_valid_neighbors(const MeshBase & mesh,
2213 : bool assert_valid_remote_elems)
2214 : {
2215 15162 : LOG_SCOPE("libmesh_assert_valid_neighbors()", "MeshTools");
2216 :
2217 2677812 : for (const auto & elem : mesh.element_ptr_range())
2218 : {
2219 2662650 : libmesh_assert (elem);
2220 2662650 : elem->libmesh_assert_valid_neighbors();
2221 : }
2222 :
2223 15162 : if (mesh.n_processors() == 1)
2224 574 : return;
2225 :
2226 14588 : libmesh_parallel_only(mesh.comm());
2227 :
2228 14588 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2229 14588 : mesh.comm().max(pmax_elem_id);
2230 :
2231 2779384 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
2232 : {
2233 2764796 : const Elem * elem = mesh.query_elem_ptr(i);
2234 :
2235 2764796 : const unsigned int my_n_neigh = elem ? elem->n_neighbors() : 0;
2236 2764796 : unsigned int n_neigh = my_n_neigh;
2237 2764796 : mesh.comm().max(n_neigh);
2238 2764796 : if (elem)
2239 2660502 : libmesh_assert_equal_to (my_n_neigh, n_neigh);
2240 :
2241 13327706 : for (unsigned int n = 0; n != n_neigh; ++n)
2242 : {
2243 10562910 : dof_id_type my_neighbor = DofObject::invalid_id;
2244 10562910 : dof_id_type * p_my_neighbor = nullptr;
2245 :
2246 : // If we have a non-remote_elem neighbor link, then we can
2247 : // verify it.
2248 10562910 : if (elem && elem->neighbor_ptr(n) != remote_elem)
2249 : {
2250 10485156 : p_my_neighbor = &my_neighbor;
2251 10485156 : if (elem->neighbor_ptr(n))
2252 9923316 : my_neighbor = elem->neighbor_ptr(n)->id();
2253 :
2254 : // But wait - if we haven't set remote_elem links yet then
2255 : // some nullptr links on ghost elements might be
2256 : // future-remote_elem links, so we can't verify those.
2257 20977216 : if (!assert_valid_remote_elems &&
2258 10486770 : !elem->neighbor_ptr(n) &&
2259 1614 : elem->processor_id() != mesh.processor_id())
2260 753 : p_my_neighbor = nullptr;
2261 : }
2262 10562910 : libmesh_assert(mesh.comm().semiverify(p_my_neighbor));
2263 : }
2264 : }
2265 : }
2266 : #endif // DEBUG
2267 :
2268 :
2269 :
2270 : // Functors for correct_node_proc_ids
2271 : namespace {
2272 :
2273 : typedef std::unordered_map<dof_id_type, processor_id_type> proc_id_map_type;
2274 :
2275 : struct SyncNodeSet
2276 : {
2277 : typedef unsigned char datum; // bool but without bit twiddling issues
2278 :
2279 158 : SyncNodeSet(std::unordered_set<const Node *> & _set,
2280 8347 : MeshBase & _mesh) :
2281 8347 : node_set(_set), mesh(_mesh) {}
2282 :
2283 : std::unordered_set<const Node *> & node_set;
2284 :
2285 : MeshBase & mesh;
2286 :
2287 : // ------------------------------------------------------------
2288 37825 : void gather_data (const std::vector<dof_id_type> & ids,
2289 : std::vector<datum> & data)
2290 : {
2291 : // Find whether each requested node belongs in the set
2292 37825 : data.resize(ids.size());
2293 :
2294 1830683 : for (auto i : index_range(ids))
2295 : {
2296 1792858 : const dof_id_type id = ids[i];
2297 :
2298 : // We'd better have every node we're asked for
2299 1792858 : Node * node = mesh.node_ptr(id);
2300 :
2301 : // Return if the node is in the set.
2302 2888844 : data[i] = node_set.count(node);
2303 : }
2304 37825 : }
2305 :
2306 : // ------------------------------------------------------------
2307 37825 : bool act_on_data (const std::vector<dof_id_type> & ids,
2308 : const std::vector<datum> in_set)
2309 : {
2310 143 : bool data_changed = false;
2311 :
2312 : // Add nodes we've been informed of to our own set
2313 1830683 : for (auto i : index_range(ids))
2314 : {
2315 1792858 : if (in_set[i])
2316 : {
2317 1162999 : Node * node = mesh.node_ptr(ids[i]);
2318 1162999 : if (!node_set.count(node))
2319 : {
2320 256 : node_set.insert(node);
2321 256 : data_changed = true;
2322 : }
2323 : }
2324 : }
2325 :
2326 37825 : return data_changed;
2327 : }
2328 : };
2329 :
2330 :
2331 8031 : struct NodesNotInSet
2332 : {
2333 158 : NodesNotInSet(const std::unordered_set<const Node *> _set)
2334 8189 : : node_set(_set) {}
2335 :
2336 473220 : bool operator() (const Node * node) const
2337 : {
2338 946440 : if (node_set.count(node))
2339 280716 : return false;
2340 192504 : return true;
2341 : }
2342 :
2343 : const std::unordered_set<const Node *> node_set;
2344 : };
2345 :
2346 :
2347 : struct SyncProcIdsFromMap
2348 : {
2349 : typedef processor_id_type datum;
2350 :
2351 194 : SyncProcIdsFromMap(const proc_id_map_type & _map,
2352 34711 : MeshBase & _mesh) :
2353 34711 : new_proc_ids(_map), mesh(_mesh) {}
2354 :
2355 : const proc_id_map_type & new_proc_ids;
2356 :
2357 : MeshBase & mesh;
2358 :
2359 : // ------------------------------------------------------------
2360 121429 : void gather_data (const std::vector<dof_id_type> & ids,
2361 : std::vector<datum> & data)
2362 : {
2363 : // Find the new processor id of each requested node
2364 121429 : data.resize(ids.size());
2365 :
2366 6983853 : for (auto i : index_range(ids))
2367 : {
2368 6862424 : const dof_id_type id = ids[i];
2369 :
2370 : // Return the node's new processor id if it has one, or its
2371 : // old processor id if not.
2372 6862424 : if (const auto it = new_proc_ids.find(id);
2373 47012 : it != new_proc_ids.end())
2374 6231344 : data[i] = it->second;
2375 : else
2376 : {
2377 : // We'd better find every node we're asked for
2378 631080 : const Node & node = mesh.node_ref(id);
2379 631080 : data[i] = node.processor_id();
2380 : }
2381 : }
2382 121429 : }
2383 :
2384 : // ------------------------------------------------------------
2385 121429 : void act_on_data (const std::vector<dof_id_type> & ids,
2386 : const std::vector<datum> proc_ids)
2387 : {
2388 : // Set the node processor ids we've now been informed of
2389 6983853 : for (auto i : index_range(ids))
2390 : {
2391 6862424 : Node & node = mesh.node_ref(ids[i]);
2392 6862424 : node.processor_id() = proc_ids[i];
2393 : }
2394 121429 : }
2395 : };
2396 : }
2397 :
2398 :
2399 :
2400 34711 : void correct_node_proc_ids (MeshBase & mesh)
2401 : {
2402 388 : LOG_SCOPE("correct_node_proc_ids()","MeshTools");
2403 :
2404 : // This function must be run on all processors at once
2405 194 : libmesh_parallel_only(mesh.comm());
2406 :
2407 : // We require all processors to agree on nodal processor ids before
2408 : // going into this algorithm.
2409 : #ifdef DEBUG
2410 194 : libmesh_assert_parallel_consistent_procids<Node>(mesh);
2411 : #endif
2412 :
2413 : // If we have any unpartitioned elements at this
2414 : // stage there is a problem
2415 194 : libmesh_assert (n_elem(mesh.unpartitioned_elements_begin(),
2416 : mesh.unpartitioned_elements_end()) == 0);
2417 :
2418 : // Fix nodes' processor ids. Coarsening may have left us with nodes
2419 : // which are no longer touched by any elements of the same processor
2420 : // id, and for DofMap to work we need to fix that.
2421 :
2422 : // This is harder now that libMesh no longer requires a distributed
2423 : // mesh to ghost all nodal neighbors: it is possible for two active
2424 : // elements on two different processors to share the same node in
2425 : // such a way that neither processor knows the others' element
2426 : // exists!
2427 :
2428 : // While we're at it, if this mesh is configured to allow
2429 : // repartitioning, we'll repartition *all* the nodes' processor ids
2430 : // using the canonical Node heuristic, to try and improve DoF load
2431 : // balancing. But if the mesh is disallowing repartitioning, we
2432 : // won't touch processor_id on any node where it's valid, regardless
2433 : // of whether or not it's canonical.
2434 194 : bool repartition_all_nodes = !mesh.skip_noncritical_partitioning();
2435 388 : std::unordered_set<const Node *> valid_nodes;
2436 :
2437 : // If we aren't allowed to repartition, then we're going to leave
2438 : // every node we can at its current processor_id, and *only*
2439 : // repartition the nodes whose current processor id is incompatible
2440 : // with DoFMap (because it doesn't touch an active element, e.g. due
2441 : // to coarsening)
2442 34711 : if (!repartition_all_nodes)
2443 : {
2444 1127692 : for (const auto & elem : mesh.active_element_ptr_range())
2445 5812313 : for (const auto & node : elem->node_ref_range())
2446 5188291 : if (elem->processor_id() == node.processor_id())
2447 4818885 : valid_nodes.insert(&node);
2448 :
2449 158 : SyncNodeSet syncv(valid_nodes, mesh);
2450 :
2451 : Parallel::sync_dofobject_data_by_id
2452 16536 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), syncv);
2453 : }
2454 :
2455 : // We build up a set of compatible processor ids for each node
2456 388 : proc_id_map_type new_proc_ids;
2457 :
2458 8783466 : for (auto & elem : mesh.active_element_ptr_range())
2459 : {
2460 4395105 : processor_id_type pid = elem->processor_id();
2461 :
2462 42755747 : for (auto & node : elem->node_ref_range())
2463 : {
2464 38360642 : const dof_id_type id = node.id();
2465 38360642 : if (auto it = new_proc_ids.find(id);
2466 342932 : it == new_proc_ids.end())
2467 152592 : new_proc_ids.emplace(id, pid);
2468 : else
2469 27494590 : it->second = node.choose_processor_id(it->second, pid);
2470 : }
2471 34323 : }
2472 :
2473 : // Sort the new pids to push to each processor
2474 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, processor_id_type>>>
2475 388 : ids_to_push;
2476 :
2477 24904448 : for (const auto & node : mesh.node_ptr_range())
2478 12666690 : if (const auto it = std::as_const(new_proc_ids).find(node->id());
2479 12666690 : it != new_proc_ids.end() && node->processor_id() != DofObject::invalid_processor_id)
2480 10900375 : ids_to_push[node->processor_id()].emplace_back(node->id(), /*pid=*/it->second);
2481 :
2482 : auto action_functor =
2483 156350 : [& mesh, & new_proc_ids]
2484 : (processor_id_type,
2485 11324190 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
2486 : {
2487 11023126 : for (const auto & [id, pid] : data)
2488 : {
2489 10866052 : if (const auto it = new_proc_ids.find(id);
2490 152592 : it == new_proc_ids.end())
2491 0 : new_proc_ids.emplace(id, pid);
2492 : else
2493 : {
2494 10866052 : const Node & node = mesh.node_ref(id);
2495 10866052 : it->second = node.choose_processor_id(it->second, pid);
2496 : }
2497 : }
2498 35241 : };
2499 :
2500 : Parallel::push_parallel_vector_data
2501 34711 : (mesh.comm(), ids_to_push, action_functor);
2502 :
2503 : // Now new_proc_ids is correct for every node we used to own. Let's
2504 : // ask every other processor about the nodes they used to own. But
2505 : // first we'll need to keep track of which nodes we used to own,
2506 : // lest we get them confused with nodes we newly own.
2507 388 : std::unordered_set<Node *> ex_local_nodes;
2508 8853602 : for (auto & node : mesh.local_node_ptr_range())
2509 4527242 : if (const auto it = new_proc_ids.find(node->id());
2510 4527242 : it != new_proc_ids.end() && it->second != mesh.processor_id())
2511 34372 : ex_local_nodes.insert(node);
2512 :
2513 194 : SyncProcIdsFromMap sync(new_proc_ids, mesh);
2514 34711 : if (repartition_all_nodes)
2515 : Parallel::sync_dofobject_data_by_id
2516 52692 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync);
2517 : else
2518 : {
2519 158 : NodesNotInSet nnis(valid_nodes);
2520 :
2521 : Parallel::sync_dofobject_data_by_id
2522 16536 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), nnis, sync);
2523 : }
2524 :
2525 : // And finally let's update the nodes we used to own.
2526 39298 : for (const auto & node : ex_local_nodes)
2527 : {
2528 2103 : if (valid_nodes.count(node))
2529 2083 : continue;
2530 :
2531 2504 : const dof_id_type id = node->id();
2532 10 : const proc_id_map_type::iterator it = new_proc_ids.find(id);
2533 10 : libmesh_assert(it != new_proc_ids.end());
2534 2504 : node->processor_id() = it->second;
2535 : }
2536 :
2537 : // We should still have consistent nodal processor ids coming out of
2538 : // this algorithm, but if we're allowed to repartition the mesh then
2539 : // they should be canonically correct too.
2540 : #ifdef DEBUG
2541 194 : libmesh_assert_valid_procids<Node>(mesh);
2542 : //if (repartition_all_nodes)
2543 : // libmesh_assert_canonical_node_procids(mesh);
2544 : #endif
2545 34711 : }
2546 :
2547 :
2548 :
2549 19490 : void Private::globally_renumber_nodes_and_elements (MeshBase & mesh)
2550 : {
2551 19490 : MeshCommunication().assign_global_indices(mesh);
2552 19490 : }
2553 :
2554 : } // namespace MeshTools
2555 :
2556 : } // namespace libMesh
|
