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 2300864 : FindBBox () : _bbox()
108 32882 : {}
109 :
110 40754 : FindBBox (FindBBox & other, Threads::split) :
111 40754 : _bbox(other._bbox)
112 13807 : {}
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 2353354 : void operator()(const ConstElemRange & range)
124 : {
125 27229941 : for (const auto & elem : range)
126 : {
127 1280474 : libmesh_assert(elem);
128 24876587 : _bbox.union_with(elem->loose_bounding_box());
129 : }
130 2353354 : }
131 :
132 16722 : Point & min() { return _bbox.min(); }
133 :
134 16722 : 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 13807 : void join (const FindBBox & other)
140 : {
141 40754 : _bbox.union_with(other._bbox);
142 26947 : }
143 : #endif
144 :
145 49042 : libMesh::BoundingBox & bbox ()
146 : {
147 49042 : return _bbox;
148 : }
149 :
150 : private:
151 : BoundingBox _bbox;
152 : };
153 :
154 : #ifdef DEBUG
155 2792022 : void assert_semiverify_dofobj(const Parallel::Communicator & communicator,
156 : const DofObject * d,
157 : unsigned int sysnum = libMesh::invalid_uint)
158 : {
159 2792022 : if (d)
160 : {
161 2716142 : const unsigned int n_sys = d->n_systems();
162 :
163 5432284 : std::vector<unsigned int> n_vars (n_sys, 0);
164 5804366 : for (unsigned int s = 0; s != n_sys; ++s)
165 3088224 : if (sysnum == s ||
166 : sysnum == libMesh::invalid_uint)
167 2716142 : n_vars[s] = d->n_vars(s);
168 :
169 : const unsigned int tot_n_vars =
170 2716142 : std::accumulate(n_vars.begin(), n_vars.end(), 0);
171 :
172 5432284 : std::vector<unsigned int> n_comp (tot_n_vars, 0);
173 5432284 : std::vector<dof_id_type> first_dof (tot_n_vars, 0);
174 :
175 5804366 : for (unsigned int s = 0, i=0; s != n_sys; ++s)
176 : {
177 3088224 : if (sysnum != s &&
178 : sysnum != libMesh::invalid_uint)
179 372082 : continue;
180 :
181 9555794 : for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
182 : {
183 6839652 : n_comp[i] = d->n_comp(s,v);
184 6839652 : first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : DofObject::invalid_id;
185 : }
186 : }
187 :
188 2716142 : libmesh_assert(communicator.semiverify(&n_sys));
189 2716142 : libmesh_assert(communicator.semiverify(&n_vars));
190 2716142 : libmesh_assert(communicator.semiverify(&n_comp));
191 2716142 : libmesh_assert(communicator.semiverify(&first_dof));
192 : }
193 : else
194 : {
195 75880 : const unsigned int * p_ui = nullptr;
196 75880 : const std::vector<unsigned int> * p_vui = nullptr;
197 75880 : const std::vector<dof_id_type> * p_vdid = nullptr;
198 :
199 75880 : libmesh_assert(communicator.semiverify(p_ui));
200 75880 : libmesh_assert(communicator.semiverify(p_vui));
201 75880 : libmesh_assert(communicator.semiverify(p_vui));
202 75880 : libmesh_assert(communicator.semiverify(p_vdid));
203 : }
204 2792022 : }
205 :
206 :
207 :
208 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
209 9223630 : 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 9223630 : if (d)
218 : {
219 8842018 : tempmin = tempmax = d->unique_id();
220 : }
221 : else
222 : {
223 381612 : TIMPI::Attributes<unique_id_type>::set_highest(tempmin);
224 381612 : TIMPI::Attributes<unique_id_type>::set_lowest(tempmax);
225 : }
226 9223630 : comm.min(tempmin);
227 9223630 : comm.max(tempmax);
228 18065648 : bool invalid = d && ((d->unique_id() != tempmin) ||
229 8842018 : (d->unique_id() != tempmax));
230 9223630 : comm.max(invalid);
231 :
232 : // First verify that everything is in sync
233 9223630 : libmesh_assert(!invalid);
234 :
235 : // Then verify that any remote id doesn't duplicate a local one.
236 9223630 : if (!d && tempmin == tempmax)
237 50164 : libmesh_assert(!unique_ids.count(tempmin));
238 9223630 : }
239 : #endif // LIBMESH_ENABLE_UNIQUE_ID
240 : #endif // DEBUG
241 :
242 363023 : 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 20452 : 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 2115445 : for (const auto & elem : node_to_elem_vec)
257 : {
258 : // We only care about active elements...
259 1752422 : if (elem->active())
260 : {
261 : // Which local node number is global_id?
262 1703058 : unsigned local_node_number = elem->local_node(global_id);
263 :
264 : // Make sure it was found
265 49364 : libmesh_assert_not_equal_to(local_node_number, libMesh::invalid_uint);
266 :
267 1752422 : 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 1752422 : 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 1752422 : const auto elem_order = Elem::type_to_default_order_map[elem->type()];
354 :
355 : // Index of the current edge
356 49364 : unsigned current_edge = 0;
357 :
358 1752422 : const unsigned short n_nodes = elem->n_nodes();
359 :
360 7153818 : while (current_edge < n_edges)
361 : {
362 : // Find the edge the node is on
363 152152 : bool found_edge = false;
364 13500366 : for (; current_edge<n_edges; ++current_edge)
365 12363514 : if (elem->is_node_on_edge(local_node_number, current_edge))
366 : {
367 120128 : found_edge = true;
368 120128 : break;
369 : }
370 :
371 : // Did we find one?
372 5401396 : if (found_edge)
373 : {
374 4264544 : const Node * node_to_save = nullptr;
375 :
376 : // Find another node in this element on this edge
377 34123736 : for (unsigned other_node_this_edge = 0; other_node_this_edge != n_nodes; other_node_this_edge++)
378 : {
379 54565488 : const bool both_vertices = elem->is_vertex(local_node_number) &&
380 24706296 : elem->is_vertex(other_node_this_edge);
381 39206760 : if ( elem->is_node_on_edge(other_node_this_edge, current_edge) && // On the current edge
382 29920568 : 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 5264528 : (elem_order == 1 || !both_vertices))
385 : {
386 : // We've found a nodal neighbor! Save a pointer to it..
387 4650358 : node_to_save = elem->node_ptr(other_node_this_edge);
388 :
389 : // Make sure we found something
390 130996 : libmesh_assert(node_to_save != nullptr);
391 :
392 4519362 : neighbor_set.insert(node_to_save);
393 : }
394 : }
395 : }
396 :
397 : // Keep looking for edges, node may be on more than one edge
398 5401396 : 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 352797 : neighbors.assign(neighbor_set.begin(), neighbor_set.end());
407 363023 : }
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 52441408 : for (const auto & elem : mesh.element_ptr_range())
505 145265427 : for (auto & node : elem->node_ref_range())
506 116534553 : nodes_to_elem_map[node.id()].push_back(elem->id());
507 243232 : }
508 :
509 :
510 :
511 5828 : 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 166 : nodes_to_elem_map.clear();
515 :
516 1169714 : for (const auto & elem : mesh.element_ptr_range())
517 3955391 : for (auto & node : elem->node_ref_range())
518 3347019 : nodes_to_elem_map[node.id()].push_back(elem);
519 5828 : }
520 :
521 :
522 :
523 : std::unordered_set<dof_id_type>
524 4912 : find_boundary_nodes(const MeshBase & mesh)
525 : {
526 144 : 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 1204818 : for (const auto & elem : mesh.active_element_ptr_range())
531 3355191 : for (auto s : elem->side_index_range())
532 2812954 : if (elem->neighbor_ptr(s) == nullptr) // on the boundary
533 : {
534 177970 : auto nodes_on_side = elem->nodes_on_side(s);
535 :
536 883969 : for (auto & local_id : nodes_on_side)
537 735971 : boundary_nodes.insert(elem->node_ptr(local_id)->id());
538 4624 : }
539 :
540 4912 : 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 1157128 : create_bounding_box (const MeshBase & mesh)
567 : {
568 : // This function must be run on all processors at once
569 16160 : libmesh_parallel_only(mesh.comm());
570 :
571 16160 : FindBBox find_bbox;
572 :
573 : // Start with any unpartitioned elements we know about locally
574 2298096 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(DofObject::invalid_processor_id),
575 1173288 : mesh.pid_elements_end(DofObject::invalid_processor_id)),
576 : find_bbox);
577 :
578 : // And combine with our local elements
579 1157128 : find_bbox.bbox().union_with(create_local_bounding_box(mesh));
580 :
581 : // Compare the bounding boxes across processors
582 1157128 : mesh.comm().min(find_bbox.min());
583 1157128 : mesh.comm().max(find_bbox.max());
584 :
585 1157128 : 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 1157128 : create_local_bounding_box (const MeshBase & mesh)
632 : {
633 16160 : FindBBox find_bbox;
634 :
635 2298096 : Threads::parallel_reduce (ConstElemRange (mesh.local_elements_begin(),
636 1173288 : mesh.local_elements_end()),
637 : find_bbox);
638 :
639 1157128 : 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 1654239 : unsigned int n_local_levels(const MeshBase & mesh)
797 : {
798 1654239 : unsigned int nl = 0;
799 :
800 4880010 : for (const auto & elem : as_range(mesh.local_elements_begin(),
801 36545016 : mesh.local_elements_end()))
802 33463638 : nl = std::max(elem->level() + 1, nl);
803 :
804 1654239 : return nl;
805 : }
806 :
807 :
808 :
809 1654239 : unsigned int n_levels(const MeshBase & mesh)
810 : {
811 28378 : libmesh_parallel_only(mesh.comm());
812 :
813 1654239 : unsigned int nl = n_local_levels(mesh);
814 :
815 3677020 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
816 20223092 : mesh.unpartitioned_elements_end()))
817 15260391 : nl = std::max(elem->level() + 1, nl);
818 :
819 1654239 : 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 28378 : const unsigned int paranoid_nl = paranoid_n_levels(mesh);
826 28378 : libmesh_assert_equal_to(nl, paranoid_nl);
827 : #endif
828 1654239 : return nl;
829 : }
830 :
831 :
832 :
833 38486 : unsigned int paranoid_n_levels(const MeshBase & mesh)
834 : {
835 28676 : libmesh_parallel_only(mesh.comm());
836 :
837 38486 : unsigned int nl = 0;
838 4059021 : for (const auto & elem : mesh.element_ptr_range())
839 4020237 : nl = std::max(elem->level() + 1, nl);
840 :
841 38486 : mesh.comm().max(nl);
842 38486 : 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 485942 : dof_id_type n_elem (const MeshBase::const_element_iterator & begin,
977 : const MeshBase::const_element_iterator & end)
978 : {
979 952324 : 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 363023 : 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 373249 : libmesh_map_find(nodes_to_elem_map, node.id());
1061 363023 : find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
1062 363023 : }
1063 :
1064 186588 : void find_nodal_or_face_neighbors(
1065 : const MeshBase & mesh,
1066 : const Node & node,
1067 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1068 : std::vector<const Node *> & neighbors)
1069 : {
1070 : // Find all the nodal neighbors... that is the nodes directly connected
1071 : // to this node through one edge.
1072 186588 : find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
1073 :
1074 : // If no neighbors are found, use all nodes on the containing side as
1075 : // neighbors.
1076 186588 : if (!neighbors.size())
1077 : {
1078 : // Grab the element containing node
1079 20377 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
1080 : // Find the element side containing node
1081 59285 : for (const auto &side : elem->side_index_range())
1082 : {
1083 59285 : const auto &nodes_on_side = elem->nodes_on_side(side);
1084 : const auto it =
1085 69047 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
1086 201975 : return elem->node_id(local_node_id) == node.id();
1087 3340 : });
1088 :
1089 60955 : if (it != nodes_on_side.end())
1090 : {
1091 184174 : for (const auto &local_node_id : nodes_on_side)
1092 : // No need to add node itself as a neighbor
1093 168411 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
1094 4614 : *node_ptr != node)
1095 143420 : neighbors.push_back(node_ptr);
1096 574 : break;
1097 : }
1098 : }
1099 : }
1100 5256 : libmesh_assert(neighbors.size());
1101 186588 : }
1102 :
1103 :
1104 :
1105 0 : void find_hanging_nodes_and_parents(const MeshBase & mesh,
1106 : std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes)
1107 : {
1108 : // Loop through all the elements
1109 0 : for (auto & elem : mesh.active_local_element_ptr_range())
1110 0 : if (elem->type() == QUAD4)
1111 0 : for (auto s : elem->side_index_range())
1112 : {
1113 : // Loop over the sides looking for sides that have hanging nodes
1114 : // This code is inspired by compute_proj_constraints()
1115 0 : const Elem * neigh = elem->neighbor_ptr(s);
1116 :
1117 : // If not a boundary side
1118 0 : if (neigh != nullptr)
1119 : {
1120 : // Is there a coarser element next to this one?
1121 0 : if (neigh->level() < elem->level())
1122 : {
1123 0 : const Elem * ancestor = elem;
1124 0 : while (neigh->level() < ancestor->level())
1125 0 : ancestor = ancestor->parent();
1126 0 : unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
1127 0 : libmesh_assert_less (s_neigh, neigh->n_neighbors());
1128 :
1129 : // Couple of helper uints...
1130 0 : unsigned int local_node1=0;
1131 0 : unsigned int local_node2=0;
1132 :
1133 0 : bool found_in_neighbor = false;
1134 :
1135 : // Find the two vertices that make up this side
1136 0 : while (!elem->is_node_on_side(local_node1++,s)) { }
1137 0 : local_node1--;
1138 :
1139 : // Start looking for the second one with the next node
1140 0 : local_node2=local_node1+1;
1141 :
1142 : // Find the other one
1143 0 : while (!elem->is_node_on_side(local_node2++,s)) { }
1144 0 : local_node2--;
1145 :
1146 : //Pull out their global ids:
1147 0 : dof_id_type node1 = elem->node_id(local_node1);
1148 0 : dof_id_type node2 = elem->node_id(local_node2);
1149 :
1150 : // Now find which node is present in the neighbor
1151 : // FIXME This assumes a level one rule!
1152 : // The _other_ one is the hanging node
1153 :
1154 : // First look for the first one
1155 : // FIXME could be streamlined a bit
1156 0 : for (unsigned int n=0;n<neigh->n_sides();n++)
1157 0 : if (neigh->node_id(n) == node1)
1158 0 : found_in_neighbor=true;
1159 :
1160 0 : dof_id_type hanging_node=0;
1161 :
1162 0 : if (!found_in_neighbor)
1163 0 : hanging_node=node1;
1164 : else // If it wasn't node1 then it must be node2!
1165 0 : hanging_node=node2;
1166 :
1167 : // Reset for reuse
1168 0 : local_node1=0;
1169 :
1170 : // Find the first node that makes up the side in the neighbor (these should be the parent nodes)
1171 0 : while (!neigh->is_node_on_side(local_node1++,s_neigh)) { }
1172 0 : local_node1--;
1173 :
1174 0 : local_node2=local_node1+1;
1175 :
1176 : // Find the second node...
1177 0 : while (!neigh->is_node_on_side(local_node2++,s_neigh)) { }
1178 0 : local_node2--;
1179 :
1180 : // Save them if we haven't already found the parents for this one
1181 0 : if (hanging_nodes[hanging_node].size()<2)
1182 : {
1183 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node1));
1184 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node2));
1185 : }
1186 : }
1187 : }
1188 0 : }
1189 0 : }
1190 :
1191 :
1192 :
1193 36 : void clear_spline_nodes(MeshBase & mesh)
1194 : {
1195 6 : std::vector<Elem *> nodeelem_to_delete;
1196 :
1197 8737 : for (auto & elem : mesh.element_ptr_range())
1198 5073 : if (elem->type() == NODEELEM &&
1199 4140 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
1200 4170 : nodeelem_to_delete.push_back(elem);
1201 :
1202 3 : auto & constraint_rows = mesh.get_constraint_rows();
1203 :
1204 : // All our constraint_rows ought to be for spline constraints we're
1205 : // about to get rid of.
1206 : #ifndef NDEBUG
1207 762 : for (auto & node_row : constraint_rows)
1208 2064 : for (auto pr : node_row.second)
1209 : {
1210 1305 : const Elem * elem = pr.first.first;
1211 1305 : libmesh_assert(elem->type() == NODEELEM);
1212 1305 : libmesh_assert(elem->mapping_type() == RATIONAL_BERNSTEIN_MAP);
1213 : }
1214 : #endif
1215 :
1216 3 : constraint_rows.clear();
1217 :
1218 4176 : for (Elem * elem : nodeelem_to_delete)
1219 : {
1220 690 : Node * node = elem->node_ptr(0);
1221 4140 : mesh.delete_elem(elem);
1222 4140 : mesh.delete_node(node);
1223 : }
1224 36 : }
1225 :
1226 :
1227 :
1228 936 : bool valid_is_prepared (const MeshBase & mesh)
1229 : {
1230 84 : LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools");
1231 :
1232 936 : if (!mesh.is_prepared())
1233 0 : return true;
1234 :
1235 978 : std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
1236 :
1237 : // Try preparing (without allowing repartitioning or renumbering, to
1238 : // avoid false assertion failures)
1239 68 : bool old_allow_renumbering = mesh_clone->allow_renumbering();
1240 42 : mesh_clone->allow_renumbering(false);
1241 : bool old_allow_remote_element_removal =
1242 68 : mesh_clone->allow_remote_element_removal();
1243 68 : bool old_skip_partitioning = mesh_clone->skip_partitioning();
1244 42 : mesh_clone->skip_partitioning(true);
1245 42 : mesh_clone->allow_remote_element_removal(false);
1246 936 : mesh_clone->prepare_for_use();
1247 42 : mesh_clone->allow_renumbering(old_allow_renumbering);
1248 42 : mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal);
1249 42 : mesh_clone->skip_partitioning(old_skip_partitioning);
1250 :
1251 936 : return (mesh == *mesh_clone);
1252 868 : }
1253 :
1254 :
1255 :
1256 : #ifndef NDEBUG
1257 :
1258 122 : void libmesh_assert_equal_n_systems (const MeshBase & mesh)
1259 : {
1260 244 : LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
1261 :
1262 122 : unsigned int n_sys = libMesh::invalid_uint;
1263 :
1264 10842 : for (const auto & elem : mesh.element_ptr_range())
1265 : {
1266 10720 : if (n_sys == libMesh::invalid_uint)
1267 122 : n_sys = elem->n_systems();
1268 : else
1269 10598 : libmesh_assert_equal_to (elem->n_systems(), n_sys);
1270 : }
1271 :
1272 29237 : for (const auto & node : mesh.node_ptr_range())
1273 : {
1274 29115 : if (n_sys == libMesh::invalid_uint)
1275 0 : n_sys = node->n_systems();
1276 : else
1277 29115 : libmesh_assert_equal_to (node->n_systems(), n_sys);
1278 : }
1279 122 : }
1280 :
1281 :
1282 :
1283 : #ifdef LIBMESH_ENABLE_AMR
1284 0 : void libmesh_assert_old_dof_objects (const MeshBase & mesh)
1285 : {
1286 0 : LOG_SCOPE("libmesh_assert_old_dof_objects()", "MeshTools");
1287 :
1288 0 : for (const auto & elem : mesh.element_ptr_range())
1289 : {
1290 0 : if (elem->refinement_flag() == Elem::JUST_REFINED ||
1291 0 : elem->refinement_flag() == Elem::INACTIVE)
1292 0 : continue;
1293 :
1294 0 : if (elem->has_dofs())
1295 0 : libmesh_assert(elem->get_old_dof_object());
1296 :
1297 0 : for (auto & node : elem->node_ref_range())
1298 0 : if (node.has_dofs())
1299 0 : libmesh_assert(node.get_old_dof_object());
1300 : }
1301 0 : }
1302 : #else
1303 : void libmesh_assert_old_dof_objects (const MeshBase &) {}
1304 : #endif // LIBMESH_ENABLE_AMR
1305 :
1306 :
1307 :
1308 0 : void libmesh_assert_valid_node_pointers(const MeshBase & mesh)
1309 : {
1310 0 : LOG_SCOPE("libmesh_assert_valid_node_pointers()", "MeshTools");
1311 :
1312 : // Here we specifically do not want "auto &" because we need to
1313 : // reseat the (temporary) pointer variable in the loop below,
1314 : // without modifying the original.
1315 0 : for (const Elem * elem : mesh.element_ptr_range())
1316 : {
1317 0 : libmesh_assert (elem);
1318 0 : while (elem)
1319 : {
1320 0 : elem->libmesh_assert_valid_node_pointers();
1321 0 : for (auto n : elem->neighbor_ptr_range())
1322 0 : if (n && n != remote_elem)
1323 0 : n->libmesh_assert_valid_node_pointers();
1324 :
1325 0 : libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1326 0 : elem = elem->parent();
1327 : }
1328 : }
1329 0 : }
1330 :
1331 :
1332 :
1333 8718 : void libmesh_assert_valid_remote_elems(const MeshBase & mesh)
1334 : {
1335 17436 : LOG_SCOPE("libmesh_assert_valid_remote_elems()", "MeshTools");
1336 :
1337 647157 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1338 1294314 : mesh.local_elements_end()))
1339 : {
1340 638439 : libmesh_assert (elem);
1341 :
1342 : // We currently don't allow active_local_elements to have
1343 : // remote_elem neighbors
1344 638439 : if (elem->active())
1345 2792382 : for (auto n : elem->neighbor_ptr_range())
1346 2247395 : libmesh_assert_not_equal_to (n, remote_elem);
1347 :
1348 : #ifdef LIBMESH_ENABLE_AMR
1349 638439 : const Elem * parent = elem->parent();
1350 638439 : if (parent)
1351 355722 : libmesh_assert_not_equal_to (parent, remote_elem);
1352 :
1353 : // We can only be strict about active elements' subactive
1354 : // children
1355 638439 : if (elem->active() && elem->has_children())
1356 26661 : for (auto & child : elem->child_ref_range())
1357 21348 : libmesh_assert_not_equal_to (&child, remote_elem);
1358 : #endif
1359 : }
1360 8718 : }
1361 :
1362 :
1363 :
1364 388 : void libmesh_assert_valid_elem_ids(const MeshBase & mesh)
1365 : {
1366 776 : LOG_SCOPE("libmesh_assert_valid_elem_ids()", "MeshTools");
1367 :
1368 388 : processor_id_type lastprocid = 0;
1369 388 : dof_id_type lastelemid = 0;
1370 :
1371 13378 : for (const auto & elem : mesh.active_element_ptr_range())
1372 : {
1373 12990 : libmesh_assert (elem);
1374 12990 : processor_id_type elemprocid = elem->processor_id();
1375 12990 : dof_id_type elemid = elem->id();
1376 :
1377 12990 : libmesh_assert_greater_equal (elemid, lastelemid);
1378 12990 : libmesh_assert_greater_equal (elemprocid, lastprocid);
1379 :
1380 12990 : lastelemid = elemid;
1381 12990 : lastprocid = elemprocid;
1382 : }
1383 388 : }
1384 :
1385 :
1386 :
1387 824 : void libmesh_assert_valid_amr_elem_ids(const MeshBase & mesh)
1388 : {
1389 1648 : LOG_SCOPE("libmesh_assert_valid_amr_elem_ids()", "MeshTools");
1390 :
1391 324000 : for (const auto & elem : mesh.element_ptr_range())
1392 : {
1393 323176 : libmesh_assert (elem);
1394 :
1395 323176 : const Elem * parent = elem->parent();
1396 :
1397 323176 : if (parent)
1398 : {
1399 24524 : libmesh_assert_greater_equal (elem->id(), parent->id());
1400 24524 : libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1401 : }
1402 : }
1403 824 : }
1404 :
1405 :
1406 :
1407 12276 : void libmesh_assert_valid_amr_interior_parents(const MeshBase & mesh)
1408 : {
1409 24552 : LOG_SCOPE("libmesh_assert_valid_amr_interior_parents()", "MeshTools");
1410 :
1411 1399387 : for (const auto & elem : mesh.element_ptr_range())
1412 : {
1413 1387111 : libmesh_assert (elem);
1414 :
1415 : // We can skip to the next element if we're full-dimension
1416 : // and therefore don't have any interior parents
1417 1387111 : if (elem->dim() >= LIBMESH_DIM)
1418 627230 : continue;
1419 :
1420 759881 : const Elem * ip = elem->interior_parent();
1421 :
1422 759881 : const Elem * parent = elem->parent();
1423 :
1424 759881 : if (ip && (ip != remote_elem) && parent)
1425 : {
1426 1000 : libmesh_assert_equal_to (ip->top_parent(),
1427 : elem->top_parent()->interior_parent());
1428 :
1429 1000 : if (ip->level() == elem->level())
1430 1000 : libmesh_assert_equal_to (ip->parent(),
1431 : parent->interior_parent());
1432 : else
1433 : {
1434 0 : libmesh_assert_less (ip->level(), elem->level());
1435 0 : libmesh_assert_equal_to (ip, parent->interior_parent());
1436 : }
1437 : }
1438 : }
1439 12276 : }
1440 :
1441 :
1442 :
1443 0 : void libmesh_assert_contiguous_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1444 : {
1445 0 : LOG_SCOPE("libmesh_assert_contiguous_dof_ids()", "MeshTools");
1446 :
1447 0 : if (mesh.n_processors() == 1)
1448 0 : return;
1449 :
1450 0 : libmesh_parallel_only(mesh.comm());
1451 :
1452 0 : dof_id_type min_dof_id = std::numeric_limits<dof_id_type>::max(),
1453 0 : max_dof_id = std::numeric_limits<dof_id_type>::min();
1454 :
1455 : // Figure out what our local dof id range is
1456 0 : for (const auto * node : mesh.local_node_ptr_range())
1457 : {
1458 0 : for (auto v : make_range(node->n_vars(sysnum)))
1459 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1460 : {
1461 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1462 0 : min_dof_id = std::min (min_dof_id, id);
1463 0 : max_dof_id = std::max (max_dof_id, id);
1464 : }
1465 : }
1466 :
1467 : // Make sure no other processors' ids are inside it
1468 0 : for (const auto * node : mesh.node_ptr_range())
1469 : {
1470 0 : if (node->processor_id() == mesh.processor_id())
1471 0 : continue;
1472 0 : for (auto v : make_range(node->n_vars(sysnum)))
1473 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1474 : {
1475 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1476 0 : libmesh_assert (id < min_dof_id ||
1477 : id > max_dof_id);
1478 : }
1479 : }
1480 : }
1481 :
1482 :
1483 :
1484 : template <>
1485 17734 : void libmesh_assert_topology_consistent_procids<Elem>(const MeshBase & mesh)
1486 : {
1487 35468 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1488 :
1489 : // This parameter is not used when !LIBMESH_ENABLE_AMR
1490 17734 : libmesh_ignore(mesh);
1491 :
1492 : // If we're adaptively refining, check processor ids for consistency
1493 : // between parents and children.
1494 : #ifdef LIBMESH_ENABLE_AMR
1495 :
1496 : // Ancestor elements we won't worry about, but subactive and active
1497 : // elements ought to have parents with consistent processor ids
1498 2650398 : for (const auto & elem : mesh.element_ptr_range())
1499 : {
1500 2632664 : libmesh_assert(elem);
1501 :
1502 2632664 : if (!elem->active() && !elem->subactive())
1503 288008 : continue;
1504 :
1505 2344656 : const Elem * parent = elem->parent();
1506 :
1507 2344656 : if (parent)
1508 : {
1509 1157556 : libmesh_assert(parent->has_children());
1510 1157556 : processor_id_type parent_procid = parent->processor_id();
1511 1157556 : bool matching_child_id = false;
1512 : // If we've got a remote_elem then we don't know whether
1513 : // it's responsible for the parent's processor id; all
1514 : // we can do is assume it is and let its processor fail
1515 : // an assert if there's something wrong.
1516 7222044 : for (auto & child : parent->child_ref_range())
1517 12128976 : if (&child == remote_elem ||
1518 6064488 : child.processor_id() == parent_procid)
1519 5909680 : matching_child_id = true;
1520 1157556 : libmesh_assert(matching_child_id);
1521 : }
1522 : }
1523 : #endif
1524 17734 : }
1525 :
1526 :
1527 :
1528 : template <>
1529 27010 : void libmesh_assert_topology_consistent_procids<Node>(const MeshBase & mesh)
1530 : {
1531 27010 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1532 :
1533 27010 : if (mesh.n_processors() == 1)
1534 0 : return;
1535 :
1536 27010 : libmesh_parallel_only(mesh.comm());
1537 :
1538 : // We want this test to be valid even when called after nodes have
1539 : // been added asynchronously but before they're renumbered.
1540 : //
1541 : // Plus, some code (looking at you, stitch_meshes) modifies
1542 : // DofObject ids without keeping max_elem_id()/max_node_id()
1543 : // consistent, but that's done in a safe way for performance
1544 : // reasons, so we'll play along and just figure out new max ids
1545 : // ourselves.
1546 27010 : dof_id_type parallel_max_node_id = 0;
1547 7288949 : for (const auto & node : mesh.node_ptr_range())
1548 7261939 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
1549 7261939 : node->id()+1);
1550 27010 : mesh.comm().max(parallel_max_node_id);
1551 :
1552 :
1553 54020 : std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1554 :
1555 2231090 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1556 4462180 : mesh.local_elements_end()))
1557 : {
1558 2204080 : libmesh_assert (elem);
1559 :
1560 16603196 : for (auto & node : elem->node_ref_range())
1561 : {
1562 14399116 : dof_id_type nodeid = node.id();
1563 14399116 : node_touched_by_me[nodeid] = true;
1564 : }
1565 : }
1566 54020 : std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1567 27010 : mesh.comm().max(node_touched_by_anyone);
1568 :
1569 3645908 : for (const auto & node : mesh.local_node_ptr_range())
1570 : {
1571 3618898 : libmesh_assert(node);
1572 3618898 : dof_id_type nodeid = node->id();
1573 3618898 : libmesh_assert(!node_touched_by_anyone[nodeid] ||
1574 : node_touched_by_me[nodeid]);
1575 : }
1576 : }
1577 :
1578 :
1579 :
1580 0 : void libmesh_assert_canonical_node_procids (const MeshBase & mesh)
1581 : {
1582 0 : for (const auto & elem : mesh.active_element_ptr_range())
1583 0 : for (auto & node : elem->node_ref_range())
1584 0 : libmesh_assert_equal_to
1585 : (node.processor_id(),
1586 : node.choose_processor_id(node.processor_id(),
1587 : elem->processor_id()));
1588 0 : }
1589 :
1590 :
1591 :
1592 : #ifdef LIBMESH_ENABLE_AMR
1593 1226 : void libmesh_assert_valid_refinement_tree(const MeshBase & mesh)
1594 : {
1595 2452 : LOG_SCOPE("libmesh_assert_valid_refinement_tree()", "MeshTools");
1596 :
1597 59722 : for (const auto & elem : mesh.element_ptr_range())
1598 : {
1599 58496 : libmesh_assert(elem);
1600 58496 : if (elem->has_children())
1601 17748 : for (auto & child : elem->child_ref_range())
1602 14928 : if (&child != remote_elem)
1603 14000 : libmesh_assert_equal_to (child.parent(), elem);
1604 58496 : if (elem->active())
1605 : {
1606 55676 : libmesh_assert(!elem->ancestor());
1607 55676 : libmesh_assert(!elem->subactive());
1608 : }
1609 2820 : else if (elem->ancestor())
1610 : {
1611 2820 : libmesh_assert(!elem->subactive());
1612 : }
1613 : else
1614 0 : libmesh_assert(elem->subactive());
1615 :
1616 58496 : if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1617 0 : libmesh_assert_greater(elem->p_level(), 0);
1618 : }
1619 1226 : }
1620 : #else
1621 : void libmesh_assert_valid_refinement_tree(const MeshBase &)
1622 : {
1623 : }
1624 : #endif // LIBMESH_ENABLE_AMR
1625 :
1626 : #endif // !NDEBUG
1627 :
1628 :
1629 :
1630 : #ifdef DEBUG
1631 :
1632 0 : void libmesh_assert_no_links_to_elem(const MeshBase & mesh,
1633 : const Elem * bad_elem)
1634 : {
1635 0 : for (const auto & elem : mesh.element_ptr_range())
1636 : {
1637 0 : libmesh_assert (elem);
1638 0 : libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1639 0 : for (auto n : elem->neighbor_ptr_range())
1640 0 : libmesh_assert_not_equal_to (n, bad_elem);
1641 :
1642 : #ifdef LIBMESH_ENABLE_AMR
1643 0 : if (elem->has_children())
1644 0 : for (auto & child : elem->child_ref_range())
1645 0 : libmesh_assert_not_equal_to (&child, bad_elem);
1646 : #endif
1647 : }
1648 0 : }
1649 :
1650 :
1651 566 : void libmesh_assert_equal_points (const MeshBase & mesh)
1652 : {
1653 1132 : LOG_SCOPE("libmesh_assert_equal_points()", "MeshTools");
1654 :
1655 566 : dof_id_type pmax_node_id = mesh.max_node_id();
1656 566 : mesh.comm().max(pmax_node_id);
1657 :
1658 915766 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1659 : {
1660 915200 : const Point * p = mesh.query_node_ptr(i);
1661 :
1662 915200 : libmesh_assert(mesh.comm().semiverify(p));
1663 : }
1664 566 : }
1665 :
1666 :
1667 566 : void libmesh_assert_equal_connectivity (const MeshBase & mesh)
1668 : {
1669 1132 : LOG_SCOPE("libmesh_assert_equal_connectivity()", "MeshTools");
1670 :
1671 566 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1672 566 : mesh.comm().max(pmax_elem_id);
1673 :
1674 1142058 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1675 : {
1676 1141492 : const Elem * e = mesh.query_elem_ptr(i);
1677 :
1678 2282984 : std::vector<dof_id_type> nodes;
1679 1141492 : if (e)
1680 5771288 : for (auto n : e->node_index_range())
1681 4629912 : nodes.push_back(e->node_id(n));
1682 :
1683 1141492 : libmesh_assert(mesh.comm().semiverify(e ? &nodes : nullptr));
1684 : }
1685 566 : }
1686 :
1687 :
1688 0 : void libmesh_assert_connected_nodes (const MeshBase & mesh)
1689 : {
1690 0 : LOG_SCOPE("libmesh_assert_connected_nodes()", "MeshTools");
1691 :
1692 0 : std::set<const Node *> used_nodes;
1693 :
1694 0 : for (const auto & elem : mesh.element_ptr_range())
1695 : {
1696 0 : libmesh_assert (elem);
1697 :
1698 0 : for (auto & n : elem->node_ref_range())
1699 0 : used_nodes.insert(&n);
1700 : }
1701 :
1702 0 : for (const auto & node : mesh.node_ptr_range())
1703 : {
1704 0 : libmesh_assert(node);
1705 0 : libmesh_assert(used_nodes.count(node));
1706 : }
1707 0 : }
1708 :
1709 :
1710 :
1711 13272 : void libmesh_assert_valid_constraint_rows (const MeshBase & mesh)
1712 : {
1713 13272 : libmesh_parallel_only(mesh.comm());
1714 :
1715 13272 : const auto & constraint_rows = mesh.get_constraint_rows();
1716 :
1717 13272 : bool have_constraint_rows = !constraint_rows.empty();
1718 13272 : mesh.comm().max(have_constraint_rows);
1719 13272 : if (!have_constraint_rows)
1720 13232 : return;
1721 :
1722 12700 : for (auto & row : constraint_rows)
1723 : {
1724 12660 : const Node * node = row.first;
1725 12660 : libmesh_assert(node == mesh.node_ptr(node->id()));
1726 :
1727 44126 : for (auto & pr : row.second)
1728 : {
1729 31466 : const Elem * spline_elem = pr.first.first;
1730 31466 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1731 : }
1732 : }
1733 :
1734 40 : dof_id_type pmax_node_id = mesh.max_node_id();
1735 40 : mesh.comm().max(pmax_node_id);
1736 :
1737 16990 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1738 : {
1739 16950 : const Node * node = mesh.query_node_ptr(i);
1740 :
1741 16950 : bool have_constraint = constraint_rows.count(node);
1742 :
1743 16950 : const std::size_t my_n_constraints = have_constraint ?
1744 12660 : libmesh_map_find(constraint_rows, node).size() : std::size_t(-1);
1745 16950 : const std::size_t * n_constraints = node ?
1746 16950 : &my_n_constraints : nullptr;
1747 :
1748 16950 : libmesh_assert(mesh.comm().semiverify(n_constraints));
1749 : }
1750 : }
1751 :
1752 :
1753 :
1754 33276 : void libmesh_assert_valid_boundary_ids(const MeshBase & mesh)
1755 : {
1756 33276 : LOG_SCOPE("libmesh_assert_valid_boundary_ids()", "MeshTools");
1757 :
1758 33276 : if (mesh.n_processors() == 1)
1759 688 : return;
1760 :
1761 32588 : libmesh_parallel_only(mesh.comm());
1762 :
1763 32588 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1764 :
1765 32588 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1766 32588 : mesh.comm().max(pmax_elem_id);
1767 :
1768 4440198 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1769 : {
1770 4407610 : const Elem * elem = mesh.query_elem_ptr(i);
1771 4407610 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1772 4407610 : const unsigned int my_n_edges = elem ? elem->n_edges() : 0;
1773 4407610 : const unsigned int my_n_sides = elem ? elem->n_sides() : 0;
1774 : unsigned int
1775 4407610 : n_nodes = my_n_nodes,
1776 4407610 : n_edges = my_n_edges,
1777 4407610 : n_sides = my_n_sides;
1778 :
1779 4407610 : mesh.comm().max(n_nodes);
1780 4407610 : mesh.comm().max(n_edges);
1781 4407610 : mesh.comm().max(n_sides);
1782 :
1783 4407610 : if (elem)
1784 : {
1785 4329843 : libmesh_assert_equal_to(my_n_nodes, n_nodes);
1786 4329843 : libmesh_assert_equal_to(my_n_edges, n_edges);
1787 4329843 : libmesh_assert_equal_to(my_n_sides, n_sides);
1788 : }
1789 :
1790 : // Let's test all IDs on the element with one communication
1791 : // rather than n_nodes + n_edges + n_sides communications, to
1792 : // cut down on latency in dbg modes.
1793 8815220 : std::vector<boundary_id_type> all_bcids;
1794 :
1795 33766334 : for (unsigned int n=0; n != n_nodes; ++n)
1796 : {
1797 58717448 : std::vector<boundary_id_type> bcids;
1798 29358724 : if (elem)
1799 : {
1800 29194930 : boundary_info.boundary_ids(elem->node_ptr(n), bcids);
1801 :
1802 : // Ordering of boundary ids shouldn't matter
1803 29194930 : std::sort(bcids.begin(), bcids.end());
1804 : }
1805 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1806 :
1807 29358724 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1808 58717448 : bcids.end());
1809 : // Separator
1810 29358724 : all_bcids.push_back(BoundaryInfo::invalid_id);
1811 : }
1812 :
1813 27932534 : for (unsigned short e=0; e != n_edges; ++e)
1814 : {
1815 47049848 : std::vector<boundary_id_type> bcids;
1816 :
1817 23524924 : if (elem)
1818 : {
1819 23423660 : boundary_info.edge_boundary_ids(elem, e, bcids);
1820 :
1821 : // Ordering of boundary ids shouldn't matter
1822 23423660 : std::sort(bcids.begin(), bcids.end());
1823 : }
1824 :
1825 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1826 :
1827 23524924 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1828 47049848 : bcids.end());
1829 : // Separator
1830 23524924 : all_bcids.push_back(BoundaryInfo::invalid_id);
1831 :
1832 23524924 : if (elem)
1833 : {
1834 23423660 : boundary_info.raw_edge_boundary_ids(elem, e, bcids);
1835 :
1836 : // Ordering of boundary ids shouldn't matter
1837 23423660 : std::sort(bcids.begin(), bcids.end());
1838 :
1839 23423660 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1840 46847320 : bcids.end());
1841 : // Separator
1842 23423660 : all_bcids.push_back(BoundaryInfo::invalid_id);
1843 : }
1844 :
1845 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1846 : }
1847 :
1848 21895242 : for (unsigned short s=0; s != n_sides; ++s)
1849 : {
1850 34975264 : std::vector<boundary_id_type> bcids;
1851 :
1852 17487632 : if (elem)
1853 : {
1854 17413676 : boundary_info.boundary_ids(elem, s, bcids);
1855 :
1856 : // Ordering of boundary ids shouldn't matter
1857 17413676 : std::sort(bcids.begin(), bcids.end());
1858 :
1859 17413676 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1860 34827352 : bcids.end());
1861 : // Separator
1862 17413676 : all_bcids.push_back(BoundaryInfo::invalid_id);
1863 : }
1864 :
1865 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1866 :
1867 17487632 : if (elem)
1868 : {
1869 17413676 : boundary_info.raw_boundary_ids(elem, s, bcids);
1870 :
1871 : // Ordering of boundary ids shouldn't matter
1872 17413676 : std::sort(bcids.begin(), bcids.end());
1873 :
1874 17413676 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1875 34827352 : bcids.end());
1876 : // Separator
1877 17413676 : all_bcids.push_back(BoundaryInfo::invalid_id);
1878 : }
1879 :
1880 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1881 : }
1882 :
1883 13222830 : for (unsigned short sf=0; sf != 2; ++sf)
1884 : {
1885 17630440 : std::vector<boundary_id_type> bcids;
1886 :
1887 8815220 : if (elem)
1888 : {
1889 8659686 : boundary_info.shellface_boundary_ids(elem, sf, bcids);
1890 :
1891 : // Ordering of boundary ids shouldn't matter
1892 8659686 : std::sort(bcids.begin(), bcids.end());
1893 :
1894 8659686 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1895 17319372 : bcids.end());
1896 : // Separator
1897 8659686 : all_bcids.push_back(BoundaryInfo::invalid_id);
1898 : }
1899 :
1900 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1901 :
1902 8815220 : if (elem)
1903 : {
1904 8659686 : boundary_info.raw_shellface_boundary_ids(elem, sf, bcids);
1905 :
1906 : // Ordering of boundary ids shouldn't matter
1907 8659686 : std::sort(bcids.begin(), bcids.end());
1908 :
1909 8659686 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1910 17319372 : bcids.end());
1911 : // Separator
1912 8659686 : all_bcids.push_back(BoundaryInfo::invalid_id);
1913 : }
1914 :
1915 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1916 : }
1917 :
1918 4407610 : libmesh_assert(mesh.comm().semiverify
1919 : (elem ? &all_bcids : nullptr));
1920 : }
1921 : }
1922 :
1923 :
1924 7660 : void libmesh_assert_valid_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1925 : {
1926 7660 : LOG_SCOPE("libmesh_assert_valid_dof_ids()", "MeshTools");
1927 :
1928 7660 : if (mesh.n_processors() == 1)
1929 0 : return;
1930 :
1931 7660 : libmesh_parallel_only(mesh.comm());
1932 :
1933 7660 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1934 7660 : mesh.comm().max(pmax_elem_id);
1935 :
1936 976732 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1937 969072 : assert_semiverify_dofobj(mesh.comm(),
1938 969072 : mesh.query_elem_ptr(i),
1939 : sysnum);
1940 :
1941 7660 : dof_id_type pmax_node_id = mesh.max_node_id();
1942 7660 : mesh.comm().max(pmax_node_id);
1943 :
1944 1830610 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1945 1822950 : assert_semiverify_dofobj(mesh.comm(),
1946 1822950 : mesh.query_node_ptr(i),
1947 : sysnum);
1948 : }
1949 :
1950 :
1951 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1952 25568 : void libmesh_assert_valid_unique_ids(const MeshBase & mesh)
1953 : {
1954 51136 : LOG_SCOPE("libmesh_assert_valid_unique_ids()", "MeshTools");
1955 :
1956 25568 : libmesh_parallel_only(mesh.comm());
1957 :
1958 : // First collect all the unique_ids we can see and make sure there's
1959 : // no duplicates
1960 51136 : std::unordered_set<unique_id_type> semilocal_unique_ids;
1961 :
1962 2929965 : for (auto const & elem : mesh.active_element_ptr_range())
1963 : {
1964 2904397 : libmesh_assert (!semilocal_unique_ids.count(elem->unique_id()));
1965 2904397 : semilocal_unique_ids.insert(elem->unique_id());
1966 : }
1967 :
1968 5496063 : for (auto const & node : mesh.node_ptr_range())
1969 : {
1970 5470495 : libmesh_assert (!semilocal_unique_ids.count(node->unique_id()));
1971 5470495 : semilocal_unique_ids.insert(node->unique_id());
1972 : }
1973 :
1974 : // Then make sure elements are all in sync and remote elements don't
1975 : // duplicate semilocal
1976 :
1977 25568 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1978 25568 : mesh.comm().max(pmax_elem_id);
1979 :
1980 3472460 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1981 : {
1982 3446892 : const Elem * elem = mesh.query_elem_ptr(i);
1983 3446892 : assert_dofobj_unique_id(mesh.comm(), elem, semilocal_unique_ids);
1984 : }
1985 :
1986 : // Then make sure nodes are all in sync and remote elements don't
1987 : // duplicate semilocal
1988 :
1989 25568 : dof_id_type pmax_node_id = mesh.max_node_id();
1990 25568 : mesh.comm().max(pmax_node_id);
1991 :
1992 5802306 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1993 : {
1994 5776738 : const Node * node = mesh.query_node_ptr(i);
1995 5776738 : assert_dofobj_unique_id(mesh.comm(), node, semilocal_unique_ids);
1996 : }
1997 25568 : }
1998 : #endif
1999 :
2000 0 : void libmesh_assert_consistent_distributed(const MeshBase & mesh)
2001 : {
2002 0 : libmesh_parallel_only(mesh.comm());
2003 :
2004 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2005 0 : mesh.comm().max(parallel_max_elem_id);
2006 :
2007 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2008 : {
2009 0 : const Elem * elem = mesh.query_elem_ptr(i);
2010 : processor_id_type pid =
2011 0 : elem ? elem->processor_id() : DofObject::invalid_processor_id;
2012 0 : mesh.comm().min(pid);
2013 0 : libmesh_assert(elem || pid != mesh.processor_id());
2014 : }
2015 :
2016 0 : dof_id_type parallel_max_node_id = mesh.max_node_id();
2017 0 : mesh.comm().max(parallel_max_node_id);
2018 :
2019 0 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2020 : {
2021 0 : const Node * node = mesh.query_node_ptr(i);
2022 : processor_id_type pid =
2023 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2024 0 : mesh.comm().min(pid);
2025 0 : libmesh_assert(node || pid != mesh.processor_id());
2026 : }
2027 0 : }
2028 :
2029 :
2030 0 : void libmesh_assert_consistent_distributed_nodes(const MeshBase & mesh)
2031 : {
2032 0 : libmesh_parallel_only(mesh.comm());
2033 0 : auto locator = mesh.sub_point_locator();
2034 :
2035 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2036 0 : mesh.comm().max(parallel_max_elem_id);
2037 :
2038 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2039 : {
2040 0 : const Elem * elem = mesh.query_elem_ptr(i);
2041 :
2042 0 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2043 0 : unsigned int n_nodes = my_n_nodes;
2044 0 : mesh.comm().max(n_nodes);
2045 :
2046 0 : if (n_nodes)
2047 0 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2048 :
2049 0 : for (unsigned int n=0; n != n_nodes; ++n)
2050 : {
2051 0 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2052 : processor_id_type pid =
2053 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2054 0 : mesh.comm().min(pid);
2055 0 : libmesh_assert(node || pid != mesh.processor_id());
2056 : }
2057 : }
2058 0 : }
2059 :
2060 :
2061 :
2062 : template <>
2063 17734 : void libmesh_assert_parallel_consistent_procids<Elem>(const MeshBase & mesh)
2064 : {
2065 17734 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2066 :
2067 17734 : if (mesh.n_processors() == 1)
2068 0 : return;
2069 :
2070 17734 : libmesh_parallel_only(mesh.comm());
2071 :
2072 : // Some code (looking at you, stitch_meshes) modifies DofObject ids
2073 : // without keeping max_elem_id()/max_node_id() consistent, but
2074 : // that's done in a safe way for performance reasons, so we'll play
2075 : // along and just figure out new max ids ourselves.
2076 17734 : dof_id_type parallel_max_elem_id = 0;
2077 2650398 : for (const auto & elem : mesh.element_ptr_range())
2078 2632664 : parallel_max_elem_id = std::max<dof_id_type>(parallel_max_elem_id,
2079 2632664 : elem->id()+1);
2080 17734 : mesh.comm().max(parallel_max_elem_id);
2081 :
2082 : // Check processor ids for consistency between processors
2083 :
2084 2719112 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2085 : {
2086 2701378 : const Elem * elem = mesh.query_elem_ptr(i);
2087 :
2088 : processor_id_type min_id =
2089 2632664 : elem ? elem->processor_id() :
2090 2701378 : std::numeric_limits<processor_id_type>::max();
2091 2701378 : mesh.comm().min(min_id);
2092 :
2093 : processor_id_type max_id =
2094 2632664 : elem ? elem->processor_id() :
2095 2701378 : std::numeric_limits<processor_id_type>::min();
2096 2701378 : mesh.comm().max(max_id);
2097 :
2098 2701378 : if (elem)
2099 : {
2100 2632664 : libmesh_assert_equal_to (min_id, elem->processor_id());
2101 2632664 : libmesh_assert_equal_to (max_id, elem->processor_id());
2102 : }
2103 :
2104 2701378 : if (min_id == mesh.processor_id())
2105 1276955 : libmesh_assert(elem);
2106 : }
2107 : }
2108 :
2109 :
2110 :
2111 76 : void libmesh_assert_parallel_consistent_new_node_procids(const MeshBase & mesh)
2112 : {
2113 76 : LOG_SCOPE("libmesh_assert_parallel_consistent_new_node_procids()", "MeshTools");
2114 :
2115 76 : if (mesh.n_processors() == 1)
2116 0 : return;
2117 :
2118 76 : libmesh_parallel_only(mesh.comm());
2119 :
2120 : // We want this test to hit every node when called even after nodes
2121 : // have been added asynchronously but before everything has been
2122 : // renumbered.
2123 76 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2124 76 : mesh.comm().max(parallel_max_elem_id);
2125 :
2126 152 : std::vector<bool> elem_touched_by_anyone(parallel_max_elem_id, false);
2127 :
2128 43248 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2129 : {
2130 43172 : const Elem * elem = mesh.query_elem_ptr(i);
2131 :
2132 43172 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2133 43172 : unsigned int n_nodes = my_n_nodes;
2134 43172 : mesh.comm().max(n_nodes);
2135 :
2136 43172 : if (n_nodes)
2137 15840 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2138 :
2139 173404 : for (unsigned int n=0; n != n_nodes; ++n)
2140 : {
2141 130232 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2142 130232 : const processor_id_type pid = node ? node->processor_id() : 0;
2143 130232 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2144 : }
2145 : }
2146 : }
2147 :
2148 : template <>
2149 27232 : void libmesh_assert_parallel_consistent_procids<Node>(const MeshBase & mesh)
2150 : {
2151 27232 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2152 :
2153 27232 : if (mesh.n_processors() == 1)
2154 0 : return;
2155 :
2156 27232 : libmesh_parallel_only(mesh.comm());
2157 :
2158 : // We want this test to be valid even when called even after nodes
2159 : // have been added asynchronously but before they're renumbered
2160 : //
2161 : // Plus, some code (looking at you, stitch_meshes) modifies
2162 : // DofObject ids without keeping max_elem_id()/max_node_id()
2163 : // consistent, but that's done in a safe way for performance
2164 : // reasons, so we'll play along and just figure out new max ids
2165 : // ourselves.
2166 27232 : dof_id_type parallel_max_node_id = 0;
2167 7539551 : for (const auto & node : mesh.node_ptr_range())
2168 7512319 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
2169 7512319 : node->id()+1);
2170 27232 : mesh.comm().max(parallel_max_node_id);
2171 :
2172 54464 : std::vector<bool> node_touched_by_anyone(parallel_max_node_id, false);
2173 :
2174 2277263 : for (const auto & elem : as_range(mesh.local_elements_begin(),
2175 4554526 : mesh.local_elements_end()))
2176 : {
2177 2250031 : libmesh_assert (elem);
2178 :
2179 17052680 : for (auto & node : elem->node_ref_range())
2180 : {
2181 14802649 : dof_id_type nodeid = node.id();
2182 14802649 : node_touched_by_anyone[nodeid] = true;
2183 : }
2184 : }
2185 27232 : mesh.comm().max(node_touched_by_anyone);
2186 :
2187 : // Check processor ids for consistency between processors
2188 : // on any node an element touches
2189 7958722 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2190 : {
2191 7931490 : if (!node_touched_by_anyone[i])
2192 421456 : continue;
2193 :
2194 7510034 : const Node * node = mesh.query_node_ptr(i);
2195 7510034 : const processor_id_type pid = node ? node->processor_id() : 0;
2196 :
2197 7510034 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2198 : }
2199 : }
2200 :
2201 :
2202 :
2203 : #ifdef LIBMESH_ENABLE_AMR
2204 0 : void libmesh_assert_valid_refinement_flags(const MeshBase & mesh)
2205 : {
2206 0 : LOG_SCOPE("libmesh_assert_valid_refinement_flags()", "MeshTools");
2207 :
2208 0 : libmesh_parallel_only(mesh.comm());
2209 0 : if (mesh.n_processors() == 1)
2210 0 : return;
2211 :
2212 0 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2213 0 : mesh.comm().max(pmax_elem_id);
2214 :
2215 0 : std::vector<unsigned char> my_elem_h_state(pmax_elem_id, 255);
2216 0 : std::vector<unsigned char> my_elem_p_state(pmax_elem_id, 255);
2217 :
2218 0 : for (const auto & elem : mesh.element_ptr_range())
2219 : {
2220 0 : libmesh_assert (elem);
2221 0 : dof_id_type elemid = elem->id();
2222 :
2223 0 : my_elem_h_state[elemid] =
2224 0 : static_cast<unsigned char>(elem->refinement_flag());
2225 :
2226 0 : my_elem_p_state[elemid] =
2227 0 : static_cast<unsigned char>(elem->p_refinement_flag());
2228 : }
2229 0 : std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
2230 0 : mesh.comm().min(min_elem_h_state);
2231 :
2232 0 : std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
2233 0 : mesh.comm().min(min_elem_p_state);
2234 :
2235 0 : for (dof_id_type i=0; i!= pmax_elem_id; ++i)
2236 : {
2237 0 : libmesh_assert(my_elem_h_state[i] == 255 ||
2238 : my_elem_h_state[i] == min_elem_h_state[i]);
2239 0 : libmesh_assert(my_elem_p_state[i] == 255 ||
2240 : my_elem_p_state[i] == min_elem_p_state[i]);
2241 : }
2242 : }
2243 : #else
2244 : void libmesh_assert_valid_refinement_flags(const MeshBase &)
2245 : {
2246 : }
2247 : #endif // LIBMESH_ENABLE_AMR
2248 :
2249 :
2250 :
2251 15182 : void libmesh_assert_valid_neighbors(const MeshBase & mesh,
2252 : bool assert_valid_remote_elems)
2253 : {
2254 15182 : LOG_SCOPE("libmesh_assert_valid_neighbors()", "MeshTools");
2255 :
2256 2681776 : for (const auto & elem : mesh.element_ptr_range())
2257 : {
2258 2666594 : libmesh_assert (elem);
2259 2666594 : elem->libmesh_assert_valid_neighbors();
2260 : }
2261 :
2262 15182 : if (mesh.n_processors() == 1)
2263 574 : return;
2264 :
2265 14608 : libmesh_parallel_only(mesh.comm());
2266 :
2267 14608 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2268 14608 : mesh.comm().max(pmax_elem_id);
2269 :
2270 2783348 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
2271 : {
2272 2768740 : const Elem * elem = mesh.query_elem_ptr(i);
2273 :
2274 2768740 : const unsigned int my_n_neigh = elem ? elem->n_neighbors() : 0;
2275 2768740 : unsigned int n_neigh = my_n_neigh;
2276 2768740 : mesh.comm().max(n_neigh);
2277 2768740 : if (elem)
2278 2664446 : libmesh_assert_equal_to (my_n_neigh, n_neigh);
2279 :
2280 13347698 : for (unsigned int n = 0; n != n_neigh; ++n)
2281 : {
2282 10578958 : dof_id_type my_neighbor = DofObject::invalid_id;
2283 10578958 : dof_id_type * p_my_neighbor = nullptr;
2284 :
2285 : // If we have a non-remote_elem neighbor link, then we can
2286 : // verify it.
2287 10578958 : if (elem && elem->neighbor_ptr(n) != remote_elem)
2288 : {
2289 10501204 : p_my_neighbor = &my_neighbor;
2290 10501204 : if (elem->neighbor_ptr(n))
2291 9937880 : my_neighbor = elem->neighbor_ptr(n)->id();
2292 :
2293 : // But wait - if we haven't set remote_elem links yet then
2294 : // some nullptr links on ghost elements might be
2295 : // future-remote_elem links, so we can't verify those.
2296 21009312 : if (!assert_valid_remote_elems &&
2297 10502818 : !elem->neighbor_ptr(n) &&
2298 1614 : elem->processor_id() != mesh.processor_id())
2299 753 : p_my_neighbor = nullptr;
2300 : }
2301 10578958 : libmesh_assert(mesh.comm().semiverify(p_my_neighbor));
2302 : }
2303 : }
2304 : }
2305 : #endif // DEBUG
2306 :
2307 :
2308 :
2309 : // Functors for correct_node_proc_ids
2310 : namespace {
2311 :
2312 : typedef std::unordered_map<dof_id_type, processor_id_type> proc_id_map_type;
2313 :
2314 : struct SyncNodeSet
2315 : {
2316 : typedef unsigned char datum; // bool but without bit twiddling issues
2317 :
2318 158 : SyncNodeSet(std::unordered_set<const Node *> & _set,
2319 8347 : MeshBase & _mesh) :
2320 8347 : node_set(_set), mesh(_mesh) {}
2321 :
2322 : std::unordered_set<const Node *> & node_set;
2323 :
2324 : MeshBase & mesh;
2325 :
2326 : // ------------------------------------------------------------
2327 37825 : void gather_data (const std::vector<dof_id_type> & ids,
2328 : std::vector<datum> & data)
2329 : {
2330 : // Find whether each requested node belongs in the set
2331 37825 : data.resize(ids.size());
2332 :
2333 1830683 : for (auto i : index_range(ids))
2334 : {
2335 1792858 : const dof_id_type id = ids[i];
2336 :
2337 : // We'd better have every node we're asked for
2338 1792858 : Node * node = mesh.node_ptr(id);
2339 :
2340 : // Return if the node is in the set.
2341 2888844 : data[i] = node_set.count(node);
2342 : }
2343 37825 : }
2344 :
2345 : // ------------------------------------------------------------
2346 37825 : bool act_on_data (const std::vector<dof_id_type> & ids,
2347 : const std::vector<datum> in_set)
2348 : {
2349 143 : bool data_changed = false;
2350 :
2351 : // Add nodes we've been informed of to our own set
2352 1830683 : for (auto i : index_range(ids))
2353 : {
2354 1792858 : if (in_set[i])
2355 : {
2356 1162999 : Node * node = mesh.node_ptr(ids[i]);
2357 1162999 : if (!node_set.count(node))
2358 : {
2359 256 : node_set.insert(node);
2360 256 : data_changed = true;
2361 : }
2362 : }
2363 : }
2364 :
2365 37825 : return data_changed;
2366 : }
2367 : };
2368 :
2369 :
2370 8031 : struct NodesNotInSet
2371 : {
2372 158 : NodesNotInSet(const std::unordered_set<const Node *> _set)
2373 8189 : : node_set(_set) {}
2374 :
2375 473220 : bool operator() (const Node * node) const
2376 : {
2377 946440 : if (node_set.count(node))
2378 280716 : return false;
2379 192504 : return true;
2380 : }
2381 :
2382 : const std::unordered_set<const Node *> node_set;
2383 : };
2384 :
2385 :
2386 : struct SyncProcIdsFromMap
2387 : {
2388 : typedef processor_id_type datum;
2389 :
2390 194 : SyncProcIdsFromMap(const proc_id_map_type & _map,
2391 34711 : MeshBase & _mesh) :
2392 34711 : new_proc_ids(_map), mesh(_mesh) {}
2393 :
2394 : const proc_id_map_type & new_proc_ids;
2395 :
2396 : MeshBase & mesh;
2397 :
2398 : // ------------------------------------------------------------
2399 121429 : void gather_data (const std::vector<dof_id_type> & ids,
2400 : std::vector<datum> & data)
2401 : {
2402 : // Find the new processor id of each requested node
2403 121429 : data.resize(ids.size());
2404 :
2405 6983779 : for (auto i : index_range(ids))
2406 : {
2407 6862350 : const dof_id_type id = ids[i];
2408 :
2409 : // Return the node's new processor id if it has one, or its
2410 : // old processor id if not.
2411 6862350 : if (const auto it = new_proc_ids.find(id);
2412 47012 : it != new_proc_ids.end())
2413 6231270 : data[i] = it->second;
2414 : else
2415 : {
2416 : // We'd better find every node we're asked for
2417 631080 : const Node & node = mesh.node_ref(id);
2418 631080 : data[i] = node.processor_id();
2419 : }
2420 : }
2421 121429 : }
2422 :
2423 : // ------------------------------------------------------------
2424 121429 : void act_on_data (const std::vector<dof_id_type> & ids,
2425 : const std::vector<datum> proc_ids)
2426 : {
2427 : // Set the node processor ids we've now been informed of
2428 6983779 : for (auto i : index_range(ids))
2429 : {
2430 6862350 : Node & node = mesh.node_ref(ids[i]);
2431 6862350 : node.processor_id() = proc_ids[i];
2432 : }
2433 121429 : }
2434 : };
2435 : }
2436 :
2437 :
2438 :
2439 34711 : void correct_node_proc_ids (MeshBase & mesh)
2440 : {
2441 388 : LOG_SCOPE("correct_node_proc_ids()","MeshTools");
2442 :
2443 : // This function must be run on all processors at once
2444 194 : libmesh_parallel_only(mesh.comm());
2445 :
2446 : // We require all processors to agree on nodal processor ids before
2447 : // going into this algorithm.
2448 : #ifdef DEBUG
2449 194 : libmesh_assert_parallel_consistent_procids<Node>(mesh);
2450 : #endif
2451 :
2452 : // If we have any unpartitioned elements at this
2453 : // stage there is a problem
2454 194 : libmesh_assert (n_elem(mesh.unpartitioned_elements_begin(),
2455 : mesh.unpartitioned_elements_end()) == 0);
2456 :
2457 : // Fix nodes' processor ids. Coarsening may have left us with nodes
2458 : // which are no longer touched by any elements of the same processor
2459 : // id, and for DofMap to work we need to fix that.
2460 :
2461 : // This is harder now that libMesh no longer requires a distributed
2462 : // mesh to ghost all nodal neighbors: it is possible for two active
2463 : // elements on two different processors to share the same node in
2464 : // such a way that neither processor knows the others' element
2465 : // exists!
2466 :
2467 : // While we're at it, if this mesh is configured to allow
2468 : // repartitioning, we'll repartition *all* the nodes' processor ids
2469 : // using the canonical Node heuristic, to try and improve DoF load
2470 : // balancing. But if the mesh is disallowing repartitioning, we
2471 : // won't touch processor_id on any node where it's valid, regardless
2472 : // of whether or not it's canonical.
2473 194 : bool repartition_all_nodes = !mesh.skip_noncritical_partitioning();
2474 388 : std::unordered_set<const Node *> valid_nodes;
2475 :
2476 : // If we aren't allowed to repartition, then we're going to leave
2477 : // every node we can at its current processor_id, and *only*
2478 : // repartition the nodes whose current processor id is incompatible
2479 : // with DoFMap (because it doesn't touch an active element, e.g. due
2480 : // to coarsening)
2481 34711 : if (!repartition_all_nodes)
2482 : {
2483 1127692 : for (const auto & elem : mesh.active_element_ptr_range())
2484 5812313 : for (const auto & node : elem->node_ref_range())
2485 5188291 : if (elem->processor_id() == node.processor_id())
2486 4818885 : valid_nodes.insert(&node);
2487 :
2488 158 : SyncNodeSet syncv(valid_nodes, mesh);
2489 :
2490 : Parallel::sync_dofobject_data_by_id
2491 16536 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), syncv);
2492 : }
2493 :
2494 : // We build up a set of compatible processor ids for each node
2495 388 : proc_id_map_type new_proc_ids;
2496 :
2497 8783394 : for (auto & elem : mesh.active_element_ptr_range())
2498 : {
2499 4395069 : processor_id_type pid = elem->processor_id();
2500 :
2501 42755495 : for (auto & node : elem->node_ref_range())
2502 : {
2503 38360426 : const dof_id_type id = node.id();
2504 38360426 : if (auto it = new_proc_ids.find(id);
2505 342932 : it == new_proc_ids.end())
2506 152592 : new_proc_ids.emplace(id, pid);
2507 : else
2508 27494484 : it->second = node.choose_processor_id(it->second, pid);
2509 : }
2510 34323 : }
2511 :
2512 : // Sort the new pids to push to each processor
2513 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, processor_id_type>>>
2514 388 : ids_to_push;
2515 :
2516 24904228 : for (const auto & node : mesh.node_ptr_range())
2517 12666580 : if (const auto it = std::as_const(new_proc_ids).find(node->id());
2518 12666580 : it != new_proc_ids.end() && node->processor_id() != DofObject::invalid_processor_id)
2519 10900265 : ids_to_push[node->processor_id()].emplace_back(node->id(), /*pid=*/it->second);
2520 :
2521 : auto action_functor =
2522 156350 : [& mesh, & new_proc_ids]
2523 : (processor_id_type,
2524 11324080 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
2525 : {
2526 11023016 : for (const auto & [id, pid] : data)
2527 : {
2528 10865942 : if (const auto it = new_proc_ids.find(id);
2529 152592 : it == new_proc_ids.end())
2530 0 : new_proc_ids.emplace(id, pid);
2531 : else
2532 : {
2533 10865942 : const Node & node = mesh.node_ref(id);
2534 10865942 : it->second = node.choose_processor_id(it->second, pid);
2535 : }
2536 : }
2537 35241 : };
2538 :
2539 : Parallel::push_parallel_vector_data
2540 34711 : (mesh.comm(), ids_to_push, action_functor);
2541 :
2542 : // Now new_proc_ids is correct for every node we used to own. Let's
2543 : // ask every other processor about the nodes they used to own. But
2544 : // first we'll need to keep track of which nodes we used to own,
2545 : // lest we get them confused with nodes we newly own.
2546 388 : std::unordered_set<Node *> ex_local_nodes;
2547 8853530 : for (auto & node : mesh.local_node_ptr_range())
2548 4527206 : if (const auto it = new_proc_ids.find(node->id());
2549 4527206 : it != new_proc_ids.end() && it->second != mesh.processor_id())
2550 34372 : ex_local_nodes.insert(node);
2551 :
2552 194 : SyncProcIdsFromMap sync(new_proc_ids, mesh);
2553 34711 : if (repartition_all_nodes)
2554 : Parallel::sync_dofobject_data_by_id
2555 52692 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync);
2556 : else
2557 : {
2558 158 : NodesNotInSet nnis(valid_nodes);
2559 :
2560 : Parallel::sync_dofobject_data_by_id
2561 16536 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), nnis, sync);
2562 : }
2563 :
2564 : // And finally let's update the nodes we used to own.
2565 39298 : for (const auto & node : ex_local_nodes)
2566 : {
2567 2103 : if (valid_nodes.count(node))
2568 2083 : continue;
2569 :
2570 2504 : const dof_id_type id = node->id();
2571 10 : const proc_id_map_type::iterator it = new_proc_ids.find(id);
2572 10 : libmesh_assert(it != new_proc_ids.end());
2573 2504 : node->processor_id() = it->second;
2574 : }
2575 :
2576 : // We should still have consistent nodal processor ids coming out of
2577 : // this algorithm, but if we're allowed to repartition the mesh then
2578 : // they should be canonically correct too.
2579 : #ifdef DEBUG
2580 194 : libmesh_assert_valid_procids<Node>(mesh);
2581 : //if (repartition_all_nodes)
2582 : // libmesh_assert_canonical_node_procids(mesh);
2583 : #endif
2584 34711 : }
2585 :
2586 :
2587 :
2588 19490 : void Private::globally_renumber_nodes_and_elements (MeshBase & mesh)
2589 : {
2590 19490 : MeshCommunication().assign_global_indices(mesh);
2591 19490 : }
2592 :
2593 : } // namespace MeshTools
2594 :
2595 : } // namespace libMesh
|
