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