Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/elem.h"
22 : #include "libmesh/elem_range.h"
23 : #include "libmesh/libmesh_logging.h"
24 : #include "libmesh/mesh_base.h"
25 : #include "libmesh/mesh_communication.h"
26 : #include "libmesh/mesh_serializer.h"
27 : #include "libmesh/mesh_tools.h"
28 : #include "libmesh/node_range.h"
29 : #include "libmesh/parallel.h"
30 : #include "libmesh/parallel_algebra.h"
31 : #include "libmesh/parallel_ghost_sync.h"
32 : #include "libmesh/sphere.h"
33 : #include "libmesh/threads.h"
34 : #include "libmesh/enum_to_string.h"
35 : #include "libmesh/enum_elem_type.h"
36 : #include "libmesh/int_range.h"
37 : #include "libmesh/utility.h"
38 : #include "libmesh/boundary_info.h"
39 :
40 : #ifndef NDEBUG
41 : # include "libmesh/remote_elem.h"
42 : #endif
43 :
44 : // C++ includes
45 : #include <limits>
46 : #include <numeric> // for std::accumulate
47 : #include <set>
48 : #include <unordered_map>
49 : #include <unordered_set>
50 :
51 :
52 :
53 : // ------------------------------------------------------------
54 : // anonymous namespace for helper classes
55 : namespace {
56 :
57 : using namespace libMesh;
58 :
59 : /**
60 : * SumElemWeight(Range) sums the number of nodes per element
61 : * for each element in the provided range. The join() method
62 : * defines how to combine the reduction operation from two
63 : * distinct instances of this class which may be executed on
64 : * separate threads.
65 : */
66 : class SumElemWeight
67 : {
68 : public:
69 0 : SumElemWeight () :
70 0 : _weight(0)
71 0 : {}
72 :
73 0 : SumElemWeight (SumElemWeight &, Threads::split) :
74 0 : _weight(0)
75 0 : {}
76 :
77 0 : void operator()(const ConstElemRange & range)
78 : {
79 0 : for (const auto & elem : range)
80 0 : _weight += elem->n_nodes();
81 0 : }
82 :
83 0 : dof_id_type weight() const
84 0 : { return _weight; }
85 :
86 : // If we don't have threads we never need a join, and icpc yells a
87 : // warning if it sees an anonymous function that's never used
88 : #if LIBMESH_USING_THREADS
89 0 : void join (const SumElemWeight & other)
90 0 : { _weight += other.weight(); }
91 : #endif
92 :
93 : private:
94 : dof_id_type _weight;
95 : };
96 :
97 :
98 : /**
99 : * FindBBox(Range) computes the bounding box for the objects
100 : * in the specified range. This class may be split and subranges
101 : * can be executed on separate threads. The join() method
102 : * defines how the results from two separate threads are combined.
103 : */
104 : class FindBBox
105 : {
106 : public:
107 2562222 : FindBBox () : _bbox()
108 36086 : {}
109 :
110 862 : FindBBox (FindBBox & other, Threads::split) :
111 862 : _bbox(other._bbox)
112 290 : {}
113 :
114 39872 : void operator()(const ConstNodeRange & range)
115 : {
116 5108830 : for (const auto & node : range)
117 : {
118 460679 : libmesh_assert(node);
119 5068958 : _bbox.union_with(*node);
120 : }
121 39872 : }
122 :
123 2578924 : void operator()(const ConstElemRange & range)
124 : {
125 28584708 : for (const auto & elem : range)
126 : {
127 1288196 : libmesh_assert(elem);
128 26005784 : _bbox.union_with(elem->loose_bounding_box());
129 : }
130 2578924 : }
131 :
132 18326 : Point & min() { return _bbox.min(); }
133 :
134 18326 : 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 290 : void join (const FindBBox & other)
140 : {
141 862 : _bbox.union_with(other._bbox);
142 572 : }
143 : #endif
144 :
145 53846 : libMesh::BoundingBox & bbox ()
146 : {
147 53846 : return _bbox;
148 : }
149 :
150 : private:
151 : BoundingBox _bbox;
152 : };
153 :
154 : #ifdef DEBUG
155 2980902 : void assert_semiverify_dofobj(const Parallel::Communicator & communicator,
156 : const DofObject * d,
157 : unsigned int sysnum = libMesh::invalid_uint)
158 : {
159 2980902 : if (d)
160 : {
161 2905442 : const unsigned int n_sys = d->n_systems();
162 :
163 5810884 : std::vector<unsigned int> n_vars (n_sys, 0);
164 6509094 : for (unsigned int s = 0; s != n_sys; ++s)
165 3603652 : if (sysnum == s ||
166 : sysnum == libMesh::invalid_uint)
167 2905442 : n_vars[s] = d->n_vars(s);
168 :
169 : const unsigned int tot_n_vars =
170 2905442 : std::accumulate(n_vars.begin(), n_vars.end(), 0);
171 :
172 5810884 : std::vector<unsigned int> n_comp (tot_n_vars, 0);
173 5810884 : std::vector<dof_id_type> first_dof (tot_n_vars, 0);
174 :
175 6509094 : for (unsigned int s = 0, i=0; s != n_sys; ++s)
176 : {
177 3603652 : if (sysnum != s &&
178 : sysnum != libMesh::invalid_uint)
179 698210 : continue;
180 :
181 9767050 : for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
182 : {
183 6861608 : n_comp[i] = d->n_comp(s,v);
184 6861608 : first_dof[i] = n_comp[i] ? d->dof_number(s,v,0) : DofObject::invalid_id;
185 : }
186 : }
187 :
188 2905442 : libmesh_assert(communicator.semiverify(&n_sys));
189 2905442 : libmesh_assert(communicator.semiverify(&n_vars));
190 2905442 : libmesh_assert(communicator.semiverify(&n_comp));
191 2905442 : libmesh_assert(communicator.semiverify(&first_dof));
192 : }
193 : else
194 : {
195 75460 : const unsigned int * p_ui = nullptr;
196 75460 : const std::vector<unsigned int> * p_vui = nullptr;
197 75460 : const std::vector<dof_id_type> * p_vdid = nullptr;
198 :
199 75460 : libmesh_assert(communicator.semiverify(p_ui));
200 75460 : libmesh_assert(communicator.semiverify(p_vui));
201 75460 : libmesh_assert(communicator.semiverify(p_vui));
202 75460 : libmesh_assert(communicator.semiverify(p_vdid));
203 : }
204 2980902 : }
205 :
206 :
207 :
208 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
209 8378462 : 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 8378462 : if (d)
218 : {
219 7980330 : tempmin = tempmax = d->unique_id();
220 : }
221 : else
222 : {
223 398132 : TIMPI::Attributes<unique_id_type>::set_highest(tempmin);
224 398132 : TIMPI::Attributes<unique_id_type>::set_lowest(tempmax);
225 : }
226 8378462 : comm.min(tempmin);
227 8378462 : comm.max(tempmax);
228 16358792 : bool invalid = d && ((d->unique_id() != tempmin) ||
229 7980330 : (d->unique_id() != tempmax));
230 8378462 : comm.max(invalid);
231 :
232 : // First verify that everything is in sync
233 8378462 : libmesh_assert(!invalid);
234 :
235 : // Then verify that any remote id doesn't duplicate a local one.
236 8378462 : if (!d && tempmin == tempmax)
237 52668 : libmesh_assert(!unique_ids.count(tempmin));
238 8378462 : }
239 : #endif // LIBMESH_ENABLE_UNIQUE_ID
240 : #endif // DEBUG
241 :
242 225994 : 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 16036 : 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 1460302 : for (const auto & elem : node_to_elem_vec)
257 : {
258 : // We only care about active elements...
259 1234308 : if (elem->active())
260 : {
261 : // Which local node number is global_id?
262 1193204 : unsigned local_node_number = elem->local_node(global_id);
263 :
264 : // Make sure it was found
265 41104 : libmesh_assert_not_equal_to(local_node_number, libMesh::invalid_uint);
266 :
267 1234308 : 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 1234308 : if (!n_edges)
271 : {
272 1374 : switch (elem->type())
273 : {
274 592 : case EDGE2:
275 : {
276 18 : switch (local_node_number)
277 : {
278 296 : case 0:
279 : // The other node is a nodal neighbor
280 305 : neighbor_set.insert(elem->node_ptr(1));
281 296 : break;
282 :
283 296 : case 1:
284 : // The other node is a nodal neighbor
285 305 : neighbor_set.insert(elem->node_ptr(0));
286 296 : break;
287 :
288 0 : default:
289 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
290 : }
291 18 : break;
292 : }
293 :
294 498 : case EDGE3:
295 : {
296 18 : switch (local_node_number)
297 : {
298 : // The outside nodes have node 2 as a neighbor
299 356 : case 0:
300 : case 1:
301 370 : neighbor_set.insert(elem->node_ptr(2));
302 356 : break;
303 :
304 : // The middle node has the outer nodes as neighbors
305 142 : case 2:
306 146 : neighbor_set.insert(elem->node_ptr(0));
307 146 : neighbor_set.insert(elem->node_ptr(1));
308 142 : break;
309 :
310 0 : default:
311 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
312 : }
313 18 : break;
314 : }
315 :
316 284 : case EDGE4:
317 : {
318 8 : switch (local_node_number)
319 : {
320 71 : case 0:
321 : // The left-middle node is a nodal neighbor
322 73 : neighbor_set.insert(elem->node_ptr(2));
323 71 : break;
324 :
325 71 : case 1:
326 : // The right-middle node is a nodal neighbor
327 73 : neighbor_set.insert(elem->node_ptr(3));
328 71 : break;
329 :
330 : // The left-middle node
331 71 : case 2:
332 73 : neighbor_set.insert(elem->node_ptr(0));
333 73 : neighbor_set.insert(elem->node_ptr(3));
334 71 : break;
335 :
336 : // The right-middle node
337 71 : case 3:
338 73 : neighbor_set.insert(elem->node_ptr(1));
339 73 : neighbor_set.insert(elem->node_ptr(2));
340 71 : break;
341 :
342 0 : default:
343 0 : libmesh_error_msg("Invalid local node number: " << local_node_number << " found." << std::endl);
344 : }
345 8 : break;
346 : }
347 :
348 0 : default:
349 0 : libmesh_error_msg("Unrecognized ElemType: " << Utility::enum_to_string(elem->type()) << std::endl);
350 : }
351 : }
352 :
353 1234308 : const auto elem_order = Elem::type_to_default_order_map[elem->type()];
354 :
355 : // Index of the current edge
356 41104 : unsigned current_edge = 0;
357 :
358 1234308 : const unsigned short n_nodes = elem->n_nodes();
359 :
360 5024403 : while (current_edge < n_edges)
361 : {
362 : // Find the edge the node is on
363 125794 : bool found_edge = false;
364 9166991 : for (; current_edge<n_edges; ++current_edge)
365 8404404 : if (elem->is_node_on_edge(local_node_number, current_edge))
366 : {
367 99678 : found_edge = true;
368 99678 : break;
369 : }
370 :
371 : // Did we find one?
372 3790095 : if (found_edge)
373 : {
374 3027508 : const Node * node_to_save = nullptr;
375 :
376 : // Find another node in this element on this edge
377 22530524 : for (unsigned other_node_this_edge = 0; other_node_this_edge != n_nodes; other_node_this_edge++)
378 : {
379 35692072 : const bool both_vertices = elem->is_vertex(local_node_number) &&
380 16189056 : elem->is_vertex(other_node_this_edge);
381 26015146 : if ( elem->is_node_on_edge(other_node_this_edge, current_edge) && // On the current edge
382 19551572 : 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 3632856 : (elem_order == 1 || !both_vertices))
385 : {
386 : // We've found a nodal neighbor! Save a pointer to it..
387 3288908 : node_to_save = elem->node_ptr(other_node_this_edge);
388 :
389 : // Make sure we found something
390 108954 : libmesh_assert(node_to_save != nullptr);
391 :
392 3179954 : neighbor_set.insert(node_to_save);
393 : }
394 : }
395 : }
396 :
397 : // Keep looking for edges, node may be on more than one edge
398 3790095 : 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 217976 : neighbors.assign(neighbor_set.begin(), neighbor_set.end());
407 225994 : }
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 6121 : 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 6121 : if (!mesh.is_serial())
462 : libmesh_deprecated();
463 :
464 6121 : nodes_to_elem_map.resize (mesh.max_node_id());
465 :
466 157004 : for (const auto & elem : mesh.element_ptr_range())
467 300472 : for (auto & node : elem->node_ref_range())
468 : {
469 6360 : libmesh_assert_less (node.id(), nodes_to_elem_map.size());
470 6360 : libmesh_assert_less (elem->id(), mesh.n_elem());
471 :
472 230124 : nodes_to_elem_map[node.id()].push_back(elem->id());
473 5773 : }
474 6121 : }
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 256003 : 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 7584 : nodes_to_elem_map.clear();
503 :
504 62694614 : for (const auto & elem : mesh.element_ptr_range())
505 190598482 : for (auto & node : elem->node_ref_range())
506 156500179 : nodes_to_elem_map[node.id()].push_back(elem->id());
507 256003 : }
508 :
509 :
510 :
511 5823 : 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 164 : nodes_to_elem_map.clear();
515 :
516 769248 : for (const auto & elem : mesh.element_ptr_range())
517 2735901 : for (auto & node : elem->node_ref_range())
518 2328997 : nodes_to_elem_map[node.id()].push_back(elem);
519 5823 : }
520 :
521 :
522 :
523 : std::unordered_set<dof_id_type>
524 5693 : find_boundary_nodes(const MeshBase & mesh)
525 : {
526 166 : 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 1029376 : for (const auto & elem : mesh.active_element_ptr_range())
531 2888348 : for (auto s : elem->side_index_range())
532 2435962 : if (elem->neighbor_ptr(s) == nullptr) // on the boundary
533 : {
534 135188 : auto nodes_on_side = elem->nodes_on_side(s);
535 :
536 666360 : for (auto & local_id : nodes_on_side)
537 561912 : boundary_nodes.insert(elem->node_ptr(local_id)->id());
538 5361 : }
539 :
540 5693 : return boundary_nodes;
541 : }
542 :
543 : std::unordered_set<dof_id_type>
544 2711 : find_block_boundary_nodes(const MeshBase & mesh)
545 : {
546 82 : 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 805088 : for (const auto & elem : mesh.active_element_ptr_range())
551 2248252 : for (auto s : elem->side_index_range())
552 1890602 : 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 2547 : }
559 :
560 2711 : return block_boundary_nodes;
561 : }
562 :
563 :
564 :
565 : libMesh::BoundingBox
566 1289337 : create_bounding_box (const MeshBase & mesh)
567 : {
568 : // This function must be run on all processors at once
569 17760 : libmesh_parallel_only(mesh.comm());
570 :
571 17760 : FindBBox find_bbox;
572 :
573 : // Start with any unpartitioned elements we know about locally
574 2578674 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(DofObject::invalid_processor_id),
575 1307097 : mesh.pid_elements_end(DofObject::invalid_processor_id)),
576 : find_bbox);
577 :
578 : // And combine with our local elements
579 1289337 : find_bbox.bbox().union_with(create_local_bounding_box(mesh));
580 :
581 : // Compare the bounding boxes across processors
582 1289337 : mesh.comm().min(find_bbox.min());
583 1289337 : mesh.comm().max(find_bbox.max());
584 :
585 1289337 : return find_bbox.bbox();
586 : }
587 :
588 :
589 :
590 : libMesh::BoundingBox
591 19630 : create_nodal_bounding_box (const MeshBase & mesh)
592 : {
593 : // This function must be run on all processors at once
594 566 : libmesh_parallel_only(mesh.comm());
595 :
596 566 : FindBBox find_bbox;
597 :
598 : // Start with any unpartitioned nodes we know about locally
599 39260 : Threads::parallel_reduce (ConstNodeRange (mesh.pid_nodes_begin(DofObject::invalid_processor_id),
600 39260 : mesh.pid_nodes_end(DofObject::invalid_processor_id)),
601 : find_bbox);
602 :
603 : // Add our local nodes
604 39260 : Threads::parallel_reduce (ConstNodeRange (mesh.local_nodes_begin(),
605 38694 : mesh.local_nodes_end()),
606 : find_bbox);
607 :
608 : // Compare the bounding boxes across processors
609 19630 : mesh.comm().min(find_bbox.min());
610 19630 : mesh.comm().max(find_bbox.max());
611 :
612 19630 : 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 1289337 : create_local_bounding_box (const MeshBase & mesh)
632 : {
633 17760 : FindBBox find_bbox;
634 :
635 2578674 : Threads::parallel_reduce (ConstElemRange (mesh.local_elements_begin(),
636 1307097 : mesh.local_elements_end()),
637 : find_bbox);
638 :
639 1289337 : 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 5471 : unsigned int n_active_local_levels(const MeshBase & mesh)
768 : {
769 : struct LevelCounter {
770 : unsigned int nl;
771 :
772 5471 : LevelCounter () : nl(0) {}
773 :
774 0 : LevelCounter (LevelCounter &, Threads::split) :
775 0 : nl(0) {}
776 :
777 5471 : void operator()(const ConstElemRange & range) {
778 70636 : for (const Elem * elem : range)
779 69004 : nl = std::max(elem->level() + 1, nl);
780 5471 : }
781 :
782 0 : void join(const LevelCounter & other) {
783 0 : nl = std::max(nl, other.nl);
784 0 : }
785 : };
786 :
787 156 : LevelCounter counter;
788 :
789 5471 : Threads::parallel_reduce(mesh.active_local_element_stored_range(), counter);
790 :
791 5471 : return counter.nl;
792 : }
793 :
794 :
795 :
796 5471 : unsigned int n_active_levels(const MeshBase & mesh)
797 : {
798 156 : libmesh_parallel_only(mesh.comm());
799 :
800 5471 : unsigned int nl = n_active_local_levels(mesh);
801 :
802 10786 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
803 21572 : mesh.unpartitioned_elements_end()))
804 0 : if (elem->active())
805 5159 : nl = std::max(elem->level() + 1, nl);
806 :
807 5471 : mesh.comm().max(nl);
808 5471 : return nl;
809 : }
810 :
811 :
812 :
813 1795001 : unsigned int n_local_levels(const MeshBase & mesh)
814 : {
815 1795001 : unsigned int nl = 0;
816 :
817 5212851 : for (const auto & elem : as_range(mesh.local_elements_begin(),
818 38344365 : mesh.local_elements_end()))
819 34961356 : nl = std::max(elem->level() + 1, nl);
820 :
821 1795001 : return nl;
822 : }
823 :
824 :
825 :
826 1795001 : unsigned int n_levels(const MeshBase & mesh)
827 : {
828 30488 : libmesh_parallel_only(mesh.comm());
829 :
830 1795001 : unsigned int nl = n_local_levels(mesh);
831 :
832 3941826 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
833 21171893 : mesh.unpartitioned_elements_end()))
834 15786906 : nl = std::max(elem->level() + 1, nl);
835 :
836 1795001 : mesh.comm().max(nl);
837 :
838 : // n_levels() is only valid and should only be called in cases where
839 : // the mesh is validly distributed (or serialized). Let's run an
840 : // expensive test in debug mode to make sure this is such a case.
841 : #ifdef DEBUG
842 30488 : const unsigned int paranoid_nl = paranoid_n_levels(mesh);
843 30488 : libmesh_assert_equal_to(nl, paranoid_nl);
844 : #endif
845 1795001 : return nl;
846 : }
847 :
848 :
849 :
850 44667 : unsigned int paranoid_n_levels(const MeshBase & mesh)
851 : {
852 30904 : libmesh_parallel_only(mesh.comm());
853 :
854 44667 : unsigned int nl = 0;
855 4178626 : for (const auto & elem : mesh.element_ptr_range())
856 4133543 : nl = std::max(elem->level() + 1, nl);
857 :
858 44667 : mesh.comm().max(nl);
859 44667 : return nl;
860 : }
861 :
862 :
863 :
864 639 : dof_id_type n_connected_components(const MeshBase & mesh,
865 : Real constraint_tol)
866 : {
867 36 : LOG_SCOPE("n_connected_components()", "MeshTools");
868 :
869 : // Yes, I'm being lazy. This is for mesh analysis before a
870 : // simulation, not anything going in any loops.
871 639 : if (!mesh.is_serial_on_zero())
872 0 : libmesh_not_implemented();
873 :
874 639 : dof_id_type n_components = 0;
875 :
876 657 : if (mesh.processor_id())
877 : {
878 522 : mesh.comm().broadcast(n_components);
879 531 : return n_components;
880 : }
881 :
882 : // All nodes in a set here are connected (at least indirectly) to
883 : // all other nodes in the same set, but have not yet been discovered
884 : // to be connected to nodes in other sets.
885 : //
886 : // Using an unordered_set of ids rather than a set of pointers seems
887 : // to be roughly 150% faster?
888 : // typedef const Node * node_entry_type
889 : typedef dof_id_type node_entry_type;
890 9 : std::vector<std::unordered_set<node_entry_type>> components;
891 :
892 : // With a typical mesh with few components and somewhat-contiguous
893 : // ordering, vector performance should be fine. With a mesh with
894 : // many components or completely scrambled ordering, performance
895 : // can be a disaster.
896 2742 : auto find_component = [&components](node_entry_type n) {
897 6158 : for (auto & c: components)
898 4064 : if (c.find(n) != c.end())
899 54 : return &c;
900 :
901 162 : return (std::unordered_set<node_entry_type> *)(nullptr);
902 108 : };
903 :
904 : auto add_to_component =
905 2010 : [&find_component]
906 2202 : (std::unordered_set<node_entry_type> & component, node_entry_type n) {
907 : // We may already be in the desired component
908 2412 : if (component.find(n) != component.end())
909 51 : return;
910 :
911 1950 : auto current_component = find_component(n);
912 :
913 : // Didn't we *just* check this?
914 150 : libmesh_assert (&component != current_component);
915 :
916 : // If we're unknown, we should be in the desired component
917 1950 : if (!current_component)
918 147 : component.insert(n);
919 :
920 : // If we think we're in another component, it should actually be
921 : // part of the desired component
922 : else
923 : {
924 : // Merge the component likely to be smaller into the one
925 : // likely to be larger - this is orders of magnitude faster
926 : // than the other way around!
927 3 : current_component->merge(component);
928 3 : current_component->swap(component);
929 3 : libmesh_assert(current_component->empty());
930 : }
931 108 : };
932 :
933 9 : auto & constraint_rows = mesh.get_constraint_rows();
934 :
935 1659 : for (const auto & elem : mesh.element_ptr_range())
936 : {
937 : // const node_entry_type first_node = elem->node_ptr(0);
938 792 : const node_entry_type first_node = elem->node_id(0);
939 :
940 792 : auto component = find_component(first_node);
941 :
942 : // If we didn't find one, make a new one, reusing an existing
943 : // slot if possible or growing our vector if necessary
944 792 : if (!component)
945 684 : for (auto & c: components)
946 354 : if (c.empty())
947 0 : component = &c;
948 :
949 432 : if (!component)
950 219 : component = &components.emplace_back();
951 :
952 3234 : for (const Node & node : elem->node_ref_range())
953 : {
954 : // const node_entry_type n = &node;
955 396 : const node_entry_type n = node.id();
956 2376 : add_to_component(*component, n);
957 :
958 2376 : auto it = constraint_rows.find(&node);
959 2376 : if (it == constraint_rows.end())
960 195 : continue;
961 :
962 72 : for (const auto & [pr, val] : it->second)
963 : {
964 : // Ignore too-trivial constraint coefficients if
965 : // we get a non-default-0 constraint_tol
966 39 : if (std::abs(val) < constraint_tol)
967 0 : continue;
968 :
969 36 : const Elem * spline_elem = pr.first;
970 3 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
971 :
972 : const Node * spline_node =
973 36 : spline_elem->node_ptr(pr.second);
974 :
975 : // add_to_component(*component, spline_node);
976 36 : add_to_component(*component, spline_node->id());
977 : }
978 : }
979 90 : }
980 :
981 327 : for (auto & component : components)
982 219 : if (!component.empty())
983 144 : ++n_components;
984 :
985 : // We calculated this on proc 0; now let everyone else know too
986 108 : mesh.comm().broadcast(n_components);
987 :
988 108 : return n_components;
989 90 : }
990 :
991 :
992 :
993 0 : void get_not_subactive_node_ids(const MeshBase & mesh,
994 : std::set<dof_id_type> & not_subactive_node_ids)
995 : {
996 0 : for (const auto & elem : mesh.element_ptr_range())
997 0 : if (!elem->subactive())
998 0 : for (auto & n : elem->node_ref_range())
999 0 : not_subactive_node_ids.insert(n.id());
1000 0 : }
1001 :
1002 :
1003 :
1004 542613 : dof_id_type n_elem (const MeshBase::const_element_iterator & begin,
1005 : const MeshBase::const_element_iterator & end)
1006 : {
1007 1063570 : return cast_int<dof_id_type>(std::distance(begin, end));
1008 : }
1009 :
1010 :
1011 :
1012 0 : dof_id_type n_nodes (const MeshBase::const_node_iterator & begin,
1013 : const MeshBase::const_node_iterator & end)
1014 : {
1015 0 : return cast_int<dof_id_type>(std::distance(begin, end));
1016 : }
1017 :
1018 :
1019 :
1020 3198 : Real volume (const MeshBase & mesh,
1021 : unsigned int dim)
1022 : {
1023 90 : libmesh_parallel_only(mesh.comm());
1024 :
1025 3198 : if (dim == libMesh::invalid_uint)
1026 3198 : dim = mesh.mesh_dimension();
1027 :
1028 3198 : Real vol = 0;
1029 :
1030 : // first my local elements
1031 7297 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1032 35467 : mesh.local_elements_end()))
1033 11923 : if (elem->dim() == dim)
1034 14941 : vol += elem->volume();
1035 :
1036 : // then count any unpartitioned objects, once
1037 3288 : if (mesh.processor_id() == 0)
1038 1041 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1039 2082 : mesh.unpartitioned_elements_end()))
1040 0 : if (elem->dim() == dim)
1041 453 : vol += elem->volume();
1042 :
1043 3198 : mesh.comm().sum(vol);
1044 3198 : return vol;
1045 : }
1046 :
1047 :
1048 :
1049 5471 : unsigned int n_p_levels (const MeshBase & mesh)
1050 : {
1051 156 : libmesh_parallel_only(mesh.comm());
1052 :
1053 5471 : unsigned int max_p_level = 0;
1054 :
1055 : // first my local elements
1056 17480 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1057 93378 : mesh.local_elements_end()))
1058 83659 : max_p_level = std::max(elem->p_level(), max_p_level);
1059 :
1060 : // then any unpartitioned objects
1061 10786 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1062 21572 : mesh.unpartitioned_elements_end()))
1063 5159 : max_p_level = std::max(elem->p_level(), max_p_level);
1064 :
1065 5471 : mesh.comm().max(max_p_level);
1066 5471 : return max_p_level + 1;
1067 : }
1068 :
1069 :
1070 :
1071 0 : void find_nodal_neighbors(const MeshBase &,
1072 : const Node & node,
1073 : const std::vector<std::vector<const Elem *>> & nodes_to_elem_map,
1074 : std::vector<const Node *> & neighbors)
1075 : {
1076 0 : find_nodal_neighbors_helper(node.id(), nodes_to_elem_map[node.id()],
1077 : neighbors);
1078 0 : }
1079 :
1080 :
1081 :
1082 225994 : void find_nodal_neighbors(const MeshBase &,
1083 : const Node & node,
1084 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1085 : std::vector<const Node *> & neighbors)
1086 : {
1087 : const std::vector<const Elem *> node_to_elem_vec =
1088 234012 : libmesh_map_find(nodes_to_elem_map, node.id());
1089 225994 : find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
1090 225994 : }
1091 :
1092 49559 : void find_nodal_or_face_neighbors(
1093 : const MeshBase & mesh,
1094 : const Node & node,
1095 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1096 : std::vector<const Node *> & neighbors)
1097 : {
1098 : // Find all the nodal neighbors... that is the nodes directly connected
1099 : // to this node through one edge.
1100 49559 : find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
1101 :
1102 : // If no neighbors are found, use all nodes on the containing side as
1103 : // neighbors.
1104 49559 : if (!neighbors.size())
1105 : {
1106 : // Grab the element containing node
1107 6804 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
1108 : // Find the element side containing node
1109 17503 : for (const auto &side : elem->side_index_range())
1110 : {
1111 17503 : const auto &nodes_on_side = elem->nodes_on_side(side);
1112 : const auto it =
1113 23074 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
1114 64593 : return elem->node_id(local_node_id) == node.id();
1115 1946 : });
1116 :
1117 18476 : if (it != nodes_on_side.end())
1118 : {
1119 58008 : for (const auto &local_node_id : nodes_on_side)
1120 : // No need to add node itself as a neighbor
1121 53931 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
1122 2727 : *node_ptr != node)
1123 44400 : neighbors.push_back(node_ptr);
1124 347 : break;
1125 : }
1126 : }
1127 : }
1128 3048 : libmesh_assert(neighbors.size());
1129 49559 : }
1130 :
1131 :
1132 :
1133 0 : void find_hanging_nodes_and_parents(const MeshBase & mesh,
1134 : std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes)
1135 : {
1136 : // Loop through all the elements
1137 0 : for (auto & elem : mesh.active_local_element_ptr_range())
1138 0 : if (elem->type() == QUAD4)
1139 0 : for (auto s : elem->side_index_range())
1140 : {
1141 : // Loop over the sides looking for sides that have hanging nodes
1142 : // This code is inspired by compute_proj_constraints()
1143 0 : const Elem * neigh = elem->neighbor_ptr(s);
1144 :
1145 : // If not a boundary side
1146 0 : if (neigh != nullptr)
1147 : {
1148 : // Is there a coarser element next to this one?
1149 0 : if (neigh->level() < elem->level())
1150 : {
1151 0 : const Elem * ancestor = elem;
1152 0 : while (neigh->level() < ancestor->level())
1153 0 : ancestor = ancestor->parent();
1154 0 : unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
1155 0 : libmesh_assert_less (s_neigh, neigh->n_neighbors());
1156 :
1157 : // Couple of helper uints...
1158 0 : unsigned int local_node1=0;
1159 0 : unsigned int local_node2=0;
1160 :
1161 0 : bool found_in_neighbor = false;
1162 :
1163 : // Find the two vertices that make up this side
1164 0 : while (!elem->is_node_on_side(local_node1++,s)) { }
1165 0 : local_node1--;
1166 :
1167 : // Start looking for the second one with the next node
1168 0 : local_node2=local_node1+1;
1169 :
1170 : // Find the other one
1171 0 : while (!elem->is_node_on_side(local_node2++,s)) { }
1172 0 : local_node2--;
1173 :
1174 : //Pull out their global ids:
1175 0 : dof_id_type node1 = elem->node_id(local_node1);
1176 0 : dof_id_type node2 = elem->node_id(local_node2);
1177 :
1178 : // Now find which node is present in the neighbor
1179 : // FIXME This assumes a level one rule!
1180 : // The _other_ one is the hanging node
1181 :
1182 : // First look for the first one
1183 : // FIXME could be streamlined a bit
1184 0 : for (unsigned int n=0;n<neigh->n_sides();n++)
1185 0 : if (neigh->node_id(n) == node1)
1186 0 : found_in_neighbor=true;
1187 :
1188 0 : dof_id_type hanging_node=0;
1189 :
1190 0 : if (!found_in_neighbor)
1191 0 : hanging_node=node1;
1192 : else // If it wasn't node1 then it must be node2!
1193 0 : hanging_node=node2;
1194 :
1195 : // Reset for reuse
1196 0 : local_node1=0;
1197 :
1198 : // Find the first node that makes up the side in the neighbor (these should be the parent nodes)
1199 0 : while (!neigh->is_node_on_side(local_node1++,s_neigh)) { }
1200 0 : local_node1--;
1201 :
1202 0 : local_node2=local_node1+1;
1203 :
1204 : // Find the second node...
1205 0 : while (!neigh->is_node_on_side(local_node2++,s_neigh)) { }
1206 0 : local_node2--;
1207 :
1208 : // Save them if we haven't already found the parents for this one
1209 0 : if (hanging_nodes[hanging_node].size()<2)
1210 : {
1211 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node1));
1212 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node2));
1213 : }
1214 : }
1215 : }
1216 0 : }
1217 0 : }
1218 :
1219 :
1220 :
1221 36 : void clear_spline_nodes(MeshBase & mesh)
1222 : {
1223 6 : std::vector<Elem *> nodeelem_to_delete;
1224 :
1225 8737 : for (auto & elem : mesh.element_ptr_range())
1226 5073 : if (elem->type() == NODEELEM &&
1227 4140 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
1228 4170 : nodeelem_to_delete.push_back(elem);
1229 :
1230 3 : auto & constraint_rows = mesh.get_constraint_rows();
1231 :
1232 : // All our constraint_rows ought to be for spline constraints we're
1233 : // about to get rid of.
1234 : #ifndef NDEBUG
1235 762 : for (auto & node_row : constraint_rows)
1236 2064 : for (auto pr : node_row.second)
1237 : {
1238 1305 : const Elem * elem = pr.first.first;
1239 1305 : libmesh_assert(elem->type() == NODEELEM);
1240 1305 : libmesh_assert(elem->mapping_type() == RATIONAL_BERNSTEIN_MAP);
1241 : }
1242 : #endif
1243 :
1244 3 : constraint_rows.clear();
1245 :
1246 4176 : for (Elem * elem : nodeelem_to_delete)
1247 : {
1248 690 : Node * node = elem->node_ptr(0);
1249 4140 : mesh.delete_elem(elem);
1250 4140 : mesh.delete_node(node);
1251 : }
1252 36 : }
1253 :
1254 :
1255 :
1256 1790 : bool valid_is_prepared (const MeshBase & mesh)
1257 : {
1258 136 : LOG_SCOPE("valid_is_prepared()", "MeshTools");
1259 :
1260 1790 : if (!mesh.is_prepared())
1261 2 : return true;
1262 :
1263 1854 : std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
1264 :
1265 : // Try preparing (without allowing repartitioning or renumbering, to
1266 : // avoid false assertion failures)
1267 116 : bool old_allow_renumbering = mesh_clone->allow_renumbering();
1268 66 : mesh_clone->allow_renumbering(false);
1269 :
1270 116 : bool old_skip_partitioning = mesh_clone->skip_partitioning();
1271 66 : mesh_clone->skip_partitioning(true);
1272 :
1273 116 : bool old_allow_remote_element_removal = mesh_clone->allow_remote_element_removal();
1274 66 : mesh_clone->allow_remote_element_removal(false);
1275 :
1276 : // Call prepare_for_use() on the clone
1277 1788 : mesh_clone->prepare_for_use();
1278 :
1279 : // Restore original flag values
1280 66 : mesh_clone->allow_renumbering(old_allow_renumbering);
1281 66 : mesh_clone->skip_partitioning(old_skip_partitioning);
1282 66 : mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal);
1283 :
1284 : // Check whether the original and clone compare equal
1285 1788 : return (mesh == *mesh_clone);
1286 1672 : }
1287 :
1288 :
1289 :
1290 : #ifndef NDEBUG
1291 :
1292 136 : void libmesh_assert_equal_n_systems (const MeshBase & mesh)
1293 : {
1294 272 : LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
1295 :
1296 136 : unsigned int n_sys = libMesh::invalid_uint;
1297 :
1298 12632 : for (const auto & elem : mesh.element_ptr_range())
1299 : {
1300 12496 : if (n_sys == libMesh::invalid_uint)
1301 136 : n_sys = elem->n_systems();
1302 : else
1303 12360 : libmesh_assert_equal_to (elem->n_systems(), n_sys);
1304 : }
1305 :
1306 31283 : for (const auto & node : mesh.node_ptr_range())
1307 : {
1308 31147 : if (n_sys == libMesh::invalid_uint)
1309 0 : n_sys = node->n_systems();
1310 : else
1311 31147 : libmesh_assert_equal_to (node->n_systems(), n_sys);
1312 : }
1313 136 : }
1314 :
1315 :
1316 :
1317 : #ifdef LIBMESH_ENABLE_AMR
1318 0 : void libmesh_assert_old_dof_objects (const MeshBase & mesh)
1319 : {
1320 0 : LOG_SCOPE("libmesh_assert_old_dof_objects()", "MeshTools");
1321 :
1322 0 : for (const auto & elem : mesh.element_ptr_range())
1323 : {
1324 0 : if (elem->refinement_flag() == Elem::JUST_REFINED ||
1325 0 : elem->refinement_flag() == Elem::INACTIVE)
1326 0 : continue;
1327 :
1328 0 : if (elem->has_dofs())
1329 0 : libmesh_assert(elem->get_old_dof_object());
1330 :
1331 0 : for (auto & node : elem->node_ref_range())
1332 0 : if (node.has_dofs())
1333 0 : libmesh_assert(node.get_old_dof_object());
1334 : }
1335 0 : }
1336 : #else
1337 : void libmesh_assert_old_dof_objects (const MeshBase &) {}
1338 : #endif // LIBMESH_ENABLE_AMR
1339 :
1340 :
1341 :
1342 0 : void libmesh_assert_valid_node_pointers(const MeshBase & mesh)
1343 : {
1344 0 : LOG_SCOPE("libmesh_assert_valid_node_pointers()", "MeshTools");
1345 :
1346 : // Here we specifically do not want "auto &" because we need to
1347 : // reseat the (temporary) pointer variable in the loop below,
1348 : // without modifying the original.
1349 0 : for (const Elem * elem : mesh.element_ptr_range())
1350 : {
1351 0 : libmesh_assert (elem);
1352 0 : while (elem)
1353 : {
1354 0 : elem->libmesh_assert_valid_node_pointers();
1355 0 : for (auto n : elem->neighbor_ptr_range())
1356 0 : if (n && n != remote_elem)
1357 0 : n->libmesh_assert_valid_node_pointers();
1358 :
1359 0 : libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1360 0 : elem = elem->parent();
1361 : }
1362 : }
1363 0 : }
1364 :
1365 :
1366 :
1367 9518 : void libmesh_assert_valid_remote_elems(const MeshBase & mesh)
1368 : {
1369 19036 : LOG_SCOPE("libmesh_assert_valid_remote_elems()", "MeshTools");
1370 :
1371 663375 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1372 1326750 : mesh.local_elements_end()))
1373 : {
1374 653857 : libmesh_assert (elem);
1375 :
1376 : // We currently don't allow active_local_elements to have
1377 : // remote_elem neighbors
1378 653857 : if (elem->active())
1379 2821606 : for (auto n : elem->neighbor_ptr_range())
1380 2268563 : libmesh_assert_not_equal_to (n, remote_elem);
1381 :
1382 : #ifdef LIBMESH_ENABLE_AMR
1383 653857 : const Elem * parent = elem->parent();
1384 653857 : if (parent)
1385 375998 : libmesh_assert_not_equal_to (parent, remote_elem);
1386 :
1387 : // We can only be strict about active elements' subactive
1388 : // children
1389 653857 : if (elem->active() && elem->has_children())
1390 30156 : for (auto & child : elem->child_ref_range())
1391 24144 : libmesh_assert_not_equal_to (&child, remote_elem);
1392 : #endif
1393 : }
1394 9518 : }
1395 :
1396 :
1397 :
1398 392 : void libmesh_assert_valid_elem_ids(const MeshBase & mesh)
1399 : {
1400 784 : LOG_SCOPE("libmesh_assert_valid_elem_ids()", "MeshTools");
1401 :
1402 392 : processor_id_type lastprocid = 0;
1403 392 : dof_id_type lastelemid = 0;
1404 :
1405 13631 : for (const auto & elem : mesh.active_element_ptr_range())
1406 : {
1407 13239 : libmesh_assert (elem);
1408 13239 : processor_id_type elemprocid = elem->processor_id();
1409 13239 : dof_id_type elemid = elem->id();
1410 :
1411 13239 : libmesh_assert_greater_equal (elemid, lastelemid);
1412 13239 : libmesh_assert_greater_equal (elemprocid, lastprocid);
1413 :
1414 13239 : lastelemid = elemid;
1415 13239 : lastprocid = elemprocid;
1416 : }
1417 392 : }
1418 :
1419 :
1420 :
1421 854 : void libmesh_assert_valid_amr_elem_ids(const MeshBase & mesh)
1422 : {
1423 1708 : LOG_SCOPE("libmesh_assert_valid_amr_elem_ids()", "MeshTools");
1424 :
1425 324068 : for (const auto & elem : mesh.element_ptr_range())
1426 : {
1427 323214 : libmesh_assert (elem);
1428 :
1429 323214 : const Elem * parent = elem->parent();
1430 :
1431 323214 : if (parent)
1432 : {
1433 24404 : libmesh_assert_greater_equal (elem->id(), parent->id());
1434 24404 : libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1435 : }
1436 : }
1437 854 : }
1438 :
1439 :
1440 :
1441 13198 : void libmesh_assert_valid_amr_interior_parents(const MeshBase & mesh)
1442 : {
1443 26396 : LOG_SCOPE("libmesh_assert_valid_amr_interior_parents()", "MeshTools");
1444 :
1445 1430217 : for (const auto & elem : mesh.element_ptr_range())
1446 : {
1447 1417019 : libmesh_assert (elem);
1448 :
1449 : // We can skip to the next element if we're full-dimension
1450 : // and therefore don't have any interior parents
1451 1417019 : if (elem->dim() >= LIBMESH_DIM)
1452 614732 : continue;
1453 :
1454 802287 : const Elem * ip = elem->interior_parent();
1455 :
1456 802287 : const Elem * parent = elem->parent();
1457 :
1458 802287 : if (ip && (ip != remote_elem) && parent)
1459 : {
1460 1000 : libmesh_assert_equal_to (ip->top_parent(),
1461 : elem->top_parent()->interior_parent());
1462 :
1463 1000 : if (ip->level() == elem->level())
1464 1000 : libmesh_assert_equal_to (ip->parent(),
1465 : parent->interior_parent());
1466 : else
1467 : {
1468 0 : libmesh_assert_less (ip->level(), elem->level());
1469 0 : libmesh_assert_equal_to (ip, parent->interior_parent());
1470 : }
1471 : }
1472 : }
1473 13198 : }
1474 :
1475 :
1476 :
1477 0 : void libmesh_assert_contiguous_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1478 : {
1479 0 : LOG_SCOPE("libmesh_assert_contiguous_dof_ids()", "MeshTools");
1480 :
1481 0 : if (mesh.n_processors() == 1)
1482 0 : return;
1483 :
1484 0 : libmesh_parallel_only(mesh.comm());
1485 :
1486 0 : dof_id_type min_dof_id = std::numeric_limits<dof_id_type>::max(),
1487 0 : max_dof_id = std::numeric_limits<dof_id_type>::min();
1488 :
1489 : // Figure out what our local dof id range is
1490 0 : for (const auto * node : mesh.local_node_ptr_range())
1491 : {
1492 0 : for (auto v : make_range(node->n_vars(sysnum)))
1493 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1494 : {
1495 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1496 0 : min_dof_id = std::min (min_dof_id, id);
1497 0 : max_dof_id = std::max (max_dof_id, id);
1498 : }
1499 : }
1500 :
1501 : // Make sure no other processors' ids are inside it
1502 0 : for (const auto * node : mesh.node_ptr_range())
1503 : {
1504 0 : if (node->processor_id() == mesh.processor_id())
1505 0 : continue;
1506 0 : for (auto v : make_range(node->n_vars(sysnum)))
1507 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1508 : {
1509 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1510 0 : libmesh_assert (id < min_dof_id ||
1511 : id > max_dof_id);
1512 : }
1513 : }
1514 : }
1515 :
1516 :
1517 :
1518 : template <>
1519 19452 : void libmesh_assert_topology_consistent_procids<Elem>(const MeshBase & mesh)
1520 : {
1521 38904 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1522 :
1523 : // This parameter is not used when !LIBMESH_ENABLE_AMR
1524 19452 : libmesh_ignore(mesh);
1525 :
1526 : // If we're adaptively refining, check processor ids for consistency
1527 : // between parents and children.
1528 : #ifdef LIBMESH_ENABLE_AMR
1529 :
1530 : // Ancestor elements we won't worry about, but subactive and active
1531 : // elements ought to have parents with consistent processor ids
1532 2716154 : for (const auto & elem : mesh.element_ptr_range())
1533 : {
1534 2696702 : libmesh_assert(elem);
1535 :
1536 2696702 : if (!elem->active() && !elem->subactive())
1537 306136 : continue;
1538 :
1539 2390566 : const Elem * parent = elem->parent();
1540 :
1541 2390566 : if (parent)
1542 : {
1543 1220636 : libmesh_assert(parent->has_children());
1544 1220636 : processor_id_type parent_procid = parent->processor_id();
1545 1220636 : bool matching_child_id = false;
1546 : // If we've got a remote_elem then we don't know whether
1547 : // it's responsible for the parent's processor id; all
1548 : // we can do is assume it is and let its processor fail
1549 : // an assert if there's something wrong.
1550 7511300 : for (auto & child : parent->child_ref_range())
1551 12578568 : if (&child == remote_elem ||
1552 6287904 : child.processor_id() == parent_procid)
1553 6128772 : matching_child_id = true;
1554 1220636 : libmesh_assert(matching_child_id);
1555 : }
1556 : }
1557 : #endif
1558 19452 : }
1559 :
1560 :
1561 :
1562 : template <>
1563 27948 : void libmesh_assert_topology_consistent_procids<Node>(const MeshBase & mesh)
1564 : {
1565 27948 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1566 :
1567 27948 : if (mesh.n_processors() == 1)
1568 0 : return;
1569 :
1570 27948 : libmesh_parallel_only(mesh.comm());
1571 :
1572 : // We want this test to be valid even when called after nodes have
1573 : // been added asynchronously but before they're renumbered.
1574 : //
1575 : // Plus, some code (looking at you, stitch_meshes) modifies
1576 : // DofObject ids without keeping max_elem_id()/max_node_id()
1577 : // consistent, but that's done in a safe way for performance
1578 : // reasons, so we'll play along and just figure out new max ids
1579 : // ourselves.
1580 27948 : dof_id_type parallel_max_node_id = 0;
1581 7519699 : for (const auto & node : mesh.node_ptr_range())
1582 7491751 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
1583 7491751 : node->id()+1);
1584 27948 : mesh.comm().max(parallel_max_node_id);
1585 :
1586 :
1587 55896 : std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1588 :
1589 2182219 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1590 4364438 : mesh.local_elements_end()))
1591 : {
1592 2154271 : libmesh_assert (elem);
1593 :
1594 16392601 : for (auto & node : elem->node_ref_range())
1595 : {
1596 14238330 : dof_id_type nodeid = node.id();
1597 14238330 : node_touched_by_me[nodeid] = true;
1598 : }
1599 : }
1600 55896 : std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1601 27948 : mesh.comm().max(node_touched_by_anyone);
1602 :
1603 3761723 : for (const auto & node : mesh.local_node_ptr_range())
1604 : {
1605 3733775 : libmesh_assert(node);
1606 3733775 : dof_id_type nodeid = node->id();
1607 3733775 : libmesh_assert(!node_touched_by_anyone[nodeid] ||
1608 : node_touched_by_me[nodeid]);
1609 : }
1610 : }
1611 :
1612 :
1613 :
1614 0 : void libmesh_assert_canonical_node_procids (const MeshBase & mesh)
1615 : {
1616 0 : for (const auto & elem : mesh.active_element_ptr_range())
1617 0 : for (auto & node : elem->node_ref_range())
1618 0 : libmesh_assert_equal_to
1619 : (node.processor_id(),
1620 : node.choose_processor_id(node.processor_id(),
1621 : elem->processor_id()));
1622 0 : }
1623 :
1624 :
1625 :
1626 : #ifdef LIBMESH_ENABLE_AMR
1627 1282 : void libmesh_assert_valid_refinement_tree(const MeshBase & mesh)
1628 : {
1629 2564 : LOG_SCOPE("libmesh_assert_valid_refinement_tree()", "MeshTools");
1630 :
1631 67616 : for (const auto & elem : mesh.element_ptr_range())
1632 : {
1633 66334 : libmesh_assert(elem);
1634 66334 : if (elem->has_children())
1635 27268 : for (auto & child : elem->child_ref_range())
1636 23216 : if (&child != remote_elem)
1637 20800 : libmesh_assert_equal_to (child.parent(), elem);
1638 66334 : if (elem->active())
1639 : {
1640 62282 : libmesh_assert(!elem->ancestor());
1641 62282 : libmesh_assert(!elem->subactive());
1642 : }
1643 4052 : else if (elem->ancestor())
1644 : {
1645 4052 : libmesh_assert(!elem->subactive());
1646 : }
1647 : else
1648 0 : libmesh_assert(elem->subactive());
1649 :
1650 66334 : if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1651 0 : libmesh_assert_greater(elem->p_level(), 0);
1652 : }
1653 1282 : }
1654 : #else
1655 : void libmesh_assert_valid_refinement_tree(const MeshBase &)
1656 : {
1657 : }
1658 : #endif // LIBMESH_ENABLE_AMR
1659 :
1660 : #endif // !NDEBUG
1661 :
1662 :
1663 :
1664 : #ifdef DEBUG
1665 :
1666 0 : void libmesh_assert_no_links_to_elem(const MeshBase & mesh,
1667 : const Elem * bad_elem)
1668 : {
1669 0 : for (const auto & elem : mesh.element_ptr_range())
1670 : {
1671 0 : libmesh_assert (elem);
1672 0 : libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1673 0 : for (auto n : elem->neighbor_ptr_range())
1674 0 : libmesh_assert_not_equal_to (n, bad_elem);
1675 :
1676 : #ifdef LIBMESH_ENABLE_AMR
1677 0 : if (elem->has_children())
1678 0 : for (auto & child : elem->child_ref_range())
1679 0 : libmesh_assert_not_equal_to (&child, bad_elem);
1680 : #endif
1681 : }
1682 0 : }
1683 :
1684 :
1685 570 : void libmesh_assert_equal_points (const MeshBase & mesh)
1686 : {
1687 1140 : LOG_SCOPE("libmesh_assert_equal_points()", "MeshTools");
1688 :
1689 570 : dof_id_type pmax_node_id = mesh.max_node_id();
1690 570 : mesh.comm().max(pmax_node_id);
1691 :
1692 922414 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1693 : {
1694 921844 : const Point * p = mesh.query_node_ptr(i);
1695 :
1696 921844 : libmesh_assert(mesh.comm().semiverify(p));
1697 : }
1698 570 : }
1699 :
1700 :
1701 570 : void libmesh_assert_equal_connectivity (const MeshBase & mesh)
1702 : {
1703 1140 : LOG_SCOPE("libmesh_assert_equal_connectivity()", "MeshTools");
1704 :
1705 570 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1706 570 : mesh.comm().max(pmax_elem_id);
1707 :
1708 1144326 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1709 : {
1710 1143756 : const Elem * e = mesh.query_elem_ptr(i);
1711 :
1712 2287512 : std::vector<dof_id_type> nodes;
1713 1143756 : if (e)
1714 5792488 : for (auto n : e->node_index_range())
1715 4648848 : nodes.push_back(e->node_id(n));
1716 :
1717 1143756 : libmesh_assert(mesh.comm().semiverify(e ? &nodes : nullptr));
1718 : }
1719 570 : }
1720 :
1721 :
1722 0 : void libmesh_assert_connected_nodes (const MeshBase & mesh)
1723 : {
1724 0 : LOG_SCOPE("libmesh_assert_connected_nodes()", "MeshTools");
1725 :
1726 0 : std::set<const Node *> used_nodes;
1727 :
1728 0 : for (const auto & elem : mesh.element_ptr_range())
1729 : {
1730 0 : libmesh_assert (elem);
1731 :
1732 0 : for (auto & n : elem->node_ref_range())
1733 0 : used_nodes.insert(&n);
1734 : }
1735 :
1736 0 : for (const auto & node : mesh.node_ptr_range())
1737 : {
1738 0 : libmesh_assert(node);
1739 0 : libmesh_assert(used_nodes.count(node));
1740 : }
1741 0 : }
1742 :
1743 :
1744 :
1745 13556 : void libmesh_assert_valid_constraint_rows (const MeshBase & mesh)
1746 : {
1747 13556 : libmesh_parallel_only(mesh.comm());
1748 :
1749 13556 : const auto & constraint_rows = mesh.get_constraint_rows();
1750 :
1751 13556 : bool have_constraint_rows = !constraint_rows.empty();
1752 13556 : mesh.comm().max(have_constraint_rows);
1753 13556 : if (!have_constraint_rows)
1754 13516 : return;
1755 :
1756 12700 : for (auto & row : constraint_rows)
1757 : {
1758 12660 : const Node * node = row.first;
1759 12660 : libmesh_assert(node == mesh.node_ptr(node->id()));
1760 :
1761 44126 : for (auto & pr : row.second)
1762 : {
1763 31466 : const Elem * spline_elem = pr.first.first;
1764 31466 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1765 : }
1766 : }
1767 :
1768 40 : dof_id_type pmax_node_id = mesh.max_node_id();
1769 40 : mesh.comm().max(pmax_node_id);
1770 :
1771 16990 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1772 : {
1773 16950 : const Node * node = mesh.query_node_ptr(i);
1774 :
1775 16950 : bool have_constraint = constraint_rows.count(node);
1776 :
1777 16950 : const std::size_t my_n_constraints = have_constraint ?
1778 12660 : libmesh_map_find(constraint_rows, node).size() : std::size_t(-1);
1779 16950 : const std::size_t * n_constraints = node ?
1780 16950 : &my_n_constraints : nullptr;
1781 :
1782 16950 : libmesh_assert(mesh.comm().semiverify(n_constraints));
1783 : }
1784 : }
1785 :
1786 :
1787 :
1788 34040 : void libmesh_assert_valid_boundary_ids(const MeshBase & mesh)
1789 : {
1790 34040 : LOG_SCOPE("libmesh_assert_valid_boundary_ids()", "MeshTools");
1791 :
1792 34040 : if (mesh.n_processors() == 1)
1793 808 : return;
1794 :
1795 33232 : libmesh_parallel_only(mesh.comm());
1796 :
1797 33232 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1798 :
1799 33232 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1800 33232 : mesh.comm().max(pmax_elem_id);
1801 :
1802 3941484 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1803 : {
1804 3908252 : const Elem * elem = mesh.query_elem_ptr(i);
1805 3908252 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1806 3908252 : const unsigned int my_n_edges = elem ? elem->n_edges() : 0;
1807 3908252 : const unsigned int my_n_sides = elem ? elem->n_sides() : 0;
1808 : unsigned int
1809 3908252 : n_nodes = my_n_nodes,
1810 3908252 : n_edges = my_n_edges,
1811 3908252 : n_sides = my_n_sides;
1812 :
1813 3908252 : mesh.comm().max(n_nodes);
1814 3908252 : mesh.comm().max(n_edges);
1815 3908252 : mesh.comm().max(n_sides);
1816 :
1817 3908252 : if (elem)
1818 : {
1819 3828115 : libmesh_assert_equal_to(my_n_nodes, n_nodes);
1820 3828115 : libmesh_assert_equal_to(my_n_edges, n_edges);
1821 3828115 : libmesh_assert_equal_to(my_n_sides, n_sides);
1822 : }
1823 :
1824 : // Let's test all IDs on the element with one communication
1825 : // rather than n_nodes + n_edges + n_sides communications, to
1826 : // cut down on latency in dbg modes.
1827 7816504 : std::vector<boundary_id_type> all_bcids;
1828 :
1829 30837818 : for (unsigned int n=0; n != n_nodes; ++n)
1830 : {
1831 53859132 : std::vector<boundary_id_type> bcids;
1832 26929566 : if (elem)
1833 : {
1834 26754870 : boundary_info.boundary_ids(elem->node_ptr(n), bcids);
1835 :
1836 : // Ordering of boundary ids shouldn't matter
1837 26754870 : std::sort(bcids.begin(), bcids.end());
1838 : }
1839 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1840 :
1841 26929566 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1842 53859132 : bcids.end());
1843 : // Separator
1844 26929566 : all_bcids.push_back(BoundaryInfo::invalid_id);
1845 : }
1846 :
1847 24039536 : for (unsigned short e=0; e != n_edges; ++e)
1848 : {
1849 40262568 : std::vector<boundary_id_type> bcids;
1850 :
1851 20131284 : if (elem)
1852 : {
1853 20014102 : boundary_info.edge_boundary_ids(elem, e, bcids);
1854 :
1855 : // Ordering of boundary ids shouldn't matter
1856 20014102 : std::sort(bcids.begin(), bcids.end());
1857 : }
1858 :
1859 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1860 :
1861 20131284 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1862 40262568 : bcids.end());
1863 : // Separator
1864 20131284 : all_bcids.push_back(BoundaryInfo::invalid_id);
1865 :
1866 20131284 : if (elem)
1867 : {
1868 20014102 : boundary_info.raw_edge_boundary_ids(elem, e, bcids);
1869 :
1870 : // Ordering of boundary ids shouldn't matter
1871 20014102 : std::sort(bcids.begin(), bcids.end());
1872 :
1873 20014102 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1874 40028204 : bcids.end());
1875 : // Separator
1876 20014102 : all_bcids.push_back(BoundaryInfo::invalid_id);
1877 : }
1878 :
1879 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1880 : }
1881 :
1882 19307032 : for (unsigned short s=0; s != n_sides; ++s)
1883 : {
1884 30797560 : std::vector<boundary_id_type> bcids;
1885 :
1886 15398780 : if (elem)
1887 : {
1888 15316122 : boundary_info.boundary_ids(elem, s, bcids);
1889 :
1890 : // Ordering of boundary ids shouldn't matter
1891 15316122 : std::sort(bcids.begin(), bcids.end());
1892 :
1893 15316122 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1894 30632244 : bcids.end());
1895 : // Separator
1896 15316122 : all_bcids.push_back(BoundaryInfo::invalid_id);
1897 : }
1898 :
1899 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1900 :
1901 15398780 : if (elem)
1902 : {
1903 15316122 : boundary_info.raw_boundary_ids(elem, s, bcids);
1904 :
1905 : // Ordering of boundary ids shouldn't matter
1906 15316122 : std::sort(bcids.begin(), bcids.end());
1907 :
1908 15316122 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1909 30632244 : bcids.end());
1910 : // Separator
1911 15316122 : all_bcids.push_back(BoundaryInfo::invalid_id);
1912 : }
1913 :
1914 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1915 : }
1916 :
1917 11724756 : for (unsigned short sf=0; sf != 2; ++sf)
1918 : {
1919 15633008 : std::vector<boundary_id_type> bcids;
1920 :
1921 7816504 : if (elem)
1922 : {
1923 7656230 : boundary_info.shellface_boundary_ids(elem, sf, bcids);
1924 :
1925 : // Ordering of boundary ids shouldn't matter
1926 7656230 : std::sort(bcids.begin(), bcids.end());
1927 :
1928 7656230 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1929 15312460 : bcids.end());
1930 : // Separator
1931 7656230 : all_bcids.push_back(BoundaryInfo::invalid_id);
1932 : }
1933 :
1934 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1935 :
1936 7816504 : if (elem)
1937 : {
1938 7656230 : boundary_info.raw_shellface_boundary_ids(elem, sf, bcids);
1939 :
1940 : // Ordering of boundary ids shouldn't matter
1941 7656230 : std::sort(bcids.begin(), bcids.end());
1942 :
1943 7656230 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1944 15312460 : bcids.end());
1945 : // Separator
1946 7656230 : all_bcids.push_back(BoundaryInfo::invalid_id);
1947 : }
1948 :
1949 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1950 : }
1951 :
1952 3908252 : libmesh_assert(mesh.comm().semiverify
1953 : (elem ? &all_bcids : nullptr));
1954 : }
1955 : }
1956 :
1957 :
1958 7910 : void libmesh_assert_valid_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1959 : {
1960 7910 : LOG_SCOPE("libmesh_assert_valid_dof_ids()", "MeshTools");
1961 :
1962 7910 : if (mesh.n_processors() == 1)
1963 0 : return;
1964 :
1965 7910 : libmesh_parallel_only(mesh.comm());
1966 :
1967 7910 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1968 7910 : mesh.comm().max(pmax_elem_id);
1969 :
1970 1045320 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1971 1037410 : assert_semiverify_dofobj(mesh.comm(),
1972 1037410 : mesh.query_elem_ptr(i),
1973 : sysnum);
1974 :
1975 7910 : dof_id_type pmax_node_id = mesh.max_node_id();
1976 7910 : mesh.comm().max(pmax_node_id);
1977 :
1978 1951402 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1979 1943492 : assert_semiverify_dofobj(mesh.comm(),
1980 1943492 : mesh.query_node_ptr(i),
1981 : sysnum);
1982 : }
1983 :
1984 :
1985 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1986 26080 : void libmesh_assert_valid_unique_ids(const MeshBase & mesh)
1987 : {
1988 52160 : LOG_SCOPE("libmesh_assert_valid_unique_ids()", "MeshTools");
1989 :
1990 26080 : libmesh_parallel_only(mesh.comm());
1991 :
1992 : // Storage for semi-local DofObject ids.
1993 52160 : std::unordered_set<unique_id_type> semilocal_unique_ids;
1994 :
1995 26080 : auto gather_elem_ids = [&]()
1996 : {
1997 2346139 : for (auto const & elem : mesh.active_element_ptr_range())
1998 : {
1999 2320059 : auto [it, inserted] = semilocal_unique_ids.insert(elem->unique_id());
2000 2320059 : libmesh_assert(inserted);
2001 2320059 : libmesh_ignore(it);
2002 : }
2003 26080 : };
2004 :
2005 26080 : auto gather_node_ids = [&]()
2006 : {
2007 5203453 : for (auto const & node : mesh.node_ptr_range())
2008 : {
2009 5177373 : auto [it, inserted] = semilocal_unique_ids.insert(node->unique_id());
2010 5177373 : libmesh_assert(inserted);
2011 5177373 : libmesh_ignore(it);
2012 : }
2013 26080 : };
2014 :
2015 26080 : auto verify_elems = [&]()
2016 : {
2017 26080 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2018 26080 : mesh.comm().max(pmax_elem_id);
2019 :
2020 2906776 : for (auto i : make_range(pmax_elem_id))
2021 : {
2022 2880696 : const Elem * elem = mesh.query_elem_ptr(i);
2023 2880696 : assert_dofobj_unique_id(mesh.comm(), elem, semilocal_unique_ids);
2024 : }
2025 26080 : };
2026 :
2027 26080 : auto verify_nodes = [&]()
2028 : {
2029 26080 : dof_id_type pmax_node_id = mesh.max_node_id();
2030 26080 : mesh.comm().max(pmax_node_id);
2031 :
2032 5523846 : for (auto i : make_range(pmax_node_id))
2033 : {
2034 5497766 : const Node * node = mesh.query_node_ptr(i);
2035 5497766 : assert_dofobj_unique_id(mesh.comm(), node, semilocal_unique_ids);
2036 : }
2037 26080 : };
2038 :
2039 26080 : if (!mesh.allow_node_and_elem_unique_id_overlap())
2040 : {
2041 : // First collect all the unique_ids we can see and make sure there's
2042 : // no duplicates
2043 26080 : gather_elem_ids();
2044 26080 : gather_node_ids();
2045 :
2046 : // Then make sure elements/nodes are all in sync and remote
2047 : // elements don't duplicate semilocal
2048 26080 : verify_elems();
2049 26080 : verify_nodes();
2050 : }
2051 : else
2052 : {
2053 : // If the mesh allows Node and Elem unique_ids to overlap, then we only
2054 : // check for validity and uniqueness of an Elem (resp. Node) unique id
2055 : // within the set of Elem (resp. Node) unique_ids.
2056 0 : gather_elem_ids();
2057 0 : verify_elems();
2058 :
2059 : // Clear id list before checking Nodes
2060 0 : semilocal_unique_ids.clear();
2061 :
2062 : // Finally, check Nodes
2063 0 : gather_node_ids();
2064 0 : verify_nodes();
2065 : }
2066 26080 : }
2067 : #endif
2068 :
2069 0 : void libmesh_assert_consistent_distributed(const MeshBase & mesh)
2070 : {
2071 0 : libmesh_parallel_only(mesh.comm());
2072 :
2073 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2074 0 : mesh.comm().max(parallel_max_elem_id);
2075 :
2076 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2077 : {
2078 0 : const Elem * elem = mesh.query_elem_ptr(i);
2079 : processor_id_type pid =
2080 0 : elem ? elem->processor_id() : DofObject::invalid_processor_id;
2081 0 : mesh.comm().min(pid);
2082 0 : libmesh_assert(elem || pid != mesh.processor_id());
2083 : }
2084 :
2085 0 : dof_id_type parallel_max_node_id = mesh.max_node_id();
2086 0 : mesh.comm().max(parallel_max_node_id);
2087 :
2088 0 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2089 : {
2090 0 : const Node * node = mesh.query_node_ptr(i);
2091 : processor_id_type pid =
2092 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2093 0 : mesh.comm().min(pid);
2094 0 : libmesh_assert(node || pid != mesh.processor_id());
2095 : }
2096 0 : }
2097 :
2098 :
2099 0 : void libmesh_assert_consistent_distributed_nodes(const MeshBase & mesh)
2100 : {
2101 0 : libmesh_parallel_only(mesh.comm());
2102 0 : auto locator = mesh.sub_point_locator();
2103 :
2104 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2105 0 : mesh.comm().max(parallel_max_elem_id);
2106 :
2107 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2108 : {
2109 0 : const Elem * elem = mesh.query_elem_ptr(i);
2110 :
2111 0 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2112 0 : unsigned int n_nodes = my_n_nodes;
2113 0 : mesh.comm().max(n_nodes);
2114 :
2115 0 : if (n_nodes)
2116 0 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2117 :
2118 0 : for (unsigned int n=0; n != n_nodes; ++n)
2119 : {
2120 0 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2121 : processor_id_type pid =
2122 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2123 0 : mesh.comm().min(pid);
2124 0 : libmesh_assert(node || pid != mesh.processor_id());
2125 : }
2126 : }
2127 0 : }
2128 :
2129 :
2130 :
2131 : template <>
2132 19452 : void libmesh_assert_parallel_consistent_procids<Elem>(const MeshBase & mesh)
2133 : {
2134 19452 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2135 :
2136 19452 : if (mesh.n_processors() == 1)
2137 0 : return;
2138 :
2139 19452 : libmesh_parallel_only(mesh.comm());
2140 :
2141 : // Some code (looking at you, stitch_meshes) modifies DofObject ids
2142 : // without keeping max_elem_id()/max_node_id() consistent, but
2143 : // that's done in a safe way for performance reasons, so we'll play
2144 : // along and just figure out new max ids ourselves.
2145 19452 : dof_id_type parallel_max_elem_id = 0;
2146 2716154 : for (const auto & elem : mesh.element_ptr_range())
2147 2696702 : parallel_max_elem_id = std::max<dof_id_type>(parallel_max_elem_id,
2148 2696702 : elem->id()+1);
2149 19452 : mesh.comm().max(parallel_max_elem_id);
2150 :
2151 : // Check processor ids for consistency between processors
2152 :
2153 2787204 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2154 : {
2155 2767752 : const Elem * elem = mesh.query_elem_ptr(i);
2156 :
2157 : processor_id_type min_id =
2158 2696702 : elem ? elem->processor_id() :
2159 2767752 : std::numeric_limits<processor_id_type>::max();
2160 2767752 : mesh.comm().min(min_id);
2161 :
2162 : processor_id_type max_id =
2163 2696702 : elem ? elem->processor_id() :
2164 2767752 : std::numeric_limits<processor_id_type>::min();
2165 2767752 : mesh.comm().max(max_id);
2166 :
2167 2767752 : if (elem)
2168 : {
2169 2696702 : libmesh_assert_equal_to (min_id, elem->processor_id());
2170 2696702 : libmesh_assert_equal_to (max_id, elem->processor_id());
2171 : }
2172 :
2173 2767752 : if (min_id == mesh.processor_id())
2174 1309659 : libmesh_assert(elem);
2175 : }
2176 : }
2177 :
2178 :
2179 :
2180 76 : void libmesh_assert_parallel_consistent_new_node_procids(const MeshBase & mesh)
2181 : {
2182 76 : LOG_SCOPE("libmesh_assert_parallel_consistent_new_node_procids()", "MeshTools");
2183 :
2184 76 : if (mesh.n_processors() == 1)
2185 0 : return;
2186 :
2187 76 : libmesh_parallel_only(mesh.comm());
2188 :
2189 : // We want this test to hit every node when called even after nodes
2190 : // have been added asynchronously but before everything has been
2191 : // renumbered.
2192 76 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2193 76 : mesh.comm().max(parallel_max_elem_id);
2194 :
2195 152 : std::vector<bool> elem_touched_by_anyone(parallel_max_elem_id, false);
2196 :
2197 43200 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2198 : {
2199 43124 : const Elem * elem = mesh.query_elem_ptr(i);
2200 :
2201 43124 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2202 43124 : unsigned int n_nodes = my_n_nodes;
2203 43124 : mesh.comm().max(n_nodes);
2204 :
2205 43124 : if (n_nodes)
2206 15840 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2207 :
2208 173356 : for (unsigned int n=0; n != n_nodes; ++n)
2209 : {
2210 130232 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2211 130232 : const processor_id_type pid = node ? node->processor_id() : 0;
2212 130232 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2213 : }
2214 : }
2215 : }
2216 :
2217 : template <>
2218 29030 : void libmesh_assert_parallel_consistent_procids<Node>(const MeshBase & mesh)
2219 : {
2220 29030 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2221 :
2222 29030 : if (mesh.n_processors() == 1)
2223 0 : return;
2224 :
2225 29030 : libmesh_parallel_only(mesh.comm());
2226 :
2227 : // We want this test to be valid even when called even after nodes
2228 : // have been added asynchronously but before they're renumbered
2229 : //
2230 : // Plus, some code (looking at you, stitch_meshes) modifies
2231 : // DofObject ids without keeping max_elem_id()/max_node_id()
2232 : // consistent, but that's done in a safe way for performance
2233 : // reasons, so we'll play along and just figure out new max ids
2234 : // ourselves.
2235 29030 : dof_id_type parallel_max_node_id = 0;
2236 7956845 : for (const auto & node : mesh.node_ptr_range())
2237 7927815 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
2238 7927815 : node->id()+1);
2239 29030 : mesh.comm().max(parallel_max_node_id);
2240 :
2241 58060 : std::vector<bool> node_touched_by_anyone(parallel_max_node_id, false);
2242 :
2243 2390977 : for (const auto & elem : as_range(mesh.local_elements_begin(),
2244 4781954 : mesh.local_elements_end()))
2245 : {
2246 2361947 : libmesh_assert (elem);
2247 :
2248 17736362 : for (auto & node : elem->node_ref_range())
2249 : {
2250 15374415 : dof_id_type nodeid = node.id();
2251 15374415 : node_touched_by_anyone[nodeid] = true;
2252 : }
2253 : }
2254 29030 : mesh.comm().max(node_touched_by_anyone);
2255 :
2256 : // Check processor ids for consistency between processors
2257 : // on any node an element touches
2258 8381202 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2259 : {
2260 8352172 : if (!node_touched_by_anyone[i])
2261 426512 : continue;
2262 :
2263 7925660 : const Node * node = mesh.query_node_ptr(i);
2264 7925660 : const processor_id_type pid = node ? node->processor_id() : 0;
2265 :
2266 7925660 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2267 : }
2268 : }
2269 :
2270 :
2271 :
2272 : #ifdef LIBMESH_ENABLE_AMR
2273 0 : void libmesh_assert_valid_refinement_flags(const MeshBase & mesh)
2274 : {
2275 0 : LOG_SCOPE("libmesh_assert_valid_refinement_flags()", "MeshTools");
2276 :
2277 0 : libmesh_parallel_only(mesh.comm());
2278 0 : if (mesh.n_processors() == 1)
2279 0 : return;
2280 :
2281 0 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2282 0 : mesh.comm().max(pmax_elem_id);
2283 :
2284 0 : std::vector<unsigned char> my_elem_h_state(pmax_elem_id, 255);
2285 0 : std::vector<unsigned char> my_elem_p_state(pmax_elem_id, 255);
2286 :
2287 0 : for (const auto & elem : mesh.element_ptr_range())
2288 : {
2289 0 : libmesh_assert (elem);
2290 0 : dof_id_type elemid = elem->id();
2291 :
2292 0 : my_elem_h_state[elemid] =
2293 0 : static_cast<unsigned char>(elem->refinement_flag());
2294 :
2295 0 : my_elem_p_state[elemid] =
2296 0 : static_cast<unsigned char>(elem->p_refinement_flag());
2297 : }
2298 0 : std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
2299 0 : mesh.comm().min(min_elem_h_state);
2300 :
2301 0 : std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
2302 0 : mesh.comm().min(min_elem_p_state);
2303 :
2304 0 : for (dof_id_type i=0; i!= pmax_elem_id; ++i)
2305 : {
2306 0 : libmesh_assert(my_elem_h_state[i] == 255 ||
2307 : my_elem_h_state[i] == min_elem_h_state[i]);
2308 0 : libmesh_assert(my_elem_p_state[i] == 255 ||
2309 : my_elem_p_state[i] == min_elem_p_state[i]);
2310 : }
2311 : }
2312 : #else
2313 : void libmesh_assert_valid_refinement_flags(const MeshBase &)
2314 : {
2315 : }
2316 : #endif // LIBMESH_ENABLE_AMR
2317 :
2318 :
2319 :
2320 16446 : void libmesh_assert_valid_neighbors(const MeshBase & mesh,
2321 : bool assert_valid_remote_elems)
2322 : {
2323 16446 : LOG_SCOPE("libmesh_assert_valid_neighbors()", "MeshTools");
2324 :
2325 2837678 : for (const auto & elem : mesh.element_ptr_range())
2326 : {
2327 2821232 : libmesh_assert (elem);
2328 2821232 : elem->libmesh_assert_valid_neighbors();
2329 : }
2330 :
2331 16446 : if (mesh.n_processors() == 1)
2332 634 : return;
2333 :
2334 15812 : libmesh_parallel_only(mesh.comm());
2335 :
2336 15812 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2337 15812 : mesh.comm().max(pmax_elem_id);
2338 :
2339 2941844 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
2340 : {
2341 2926032 : const Elem * elem = mesh.query_elem_ptr(i);
2342 :
2343 2926032 : const unsigned int my_n_neigh = elem ? elem->n_neighbors() : 0;
2344 2926032 : unsigned int n_neigh = my_n_neigh;
2345 2926032 : mesh.comm().max(n_neigh);
2346 2926032 : if (elem)
2347 2818332 : libmesh_assert_equal_to (my_n_neigh, n_neigh);
2348 :
2349 14058598 : for (unsigned int n = 0; n != n_neigh; ++n)
2350 : {
2351 11132566 : dof_id_type my_neighbor = DofObject::invalid_id;
2352 11132566 : dof_id_type * p_my_neighbor = nullptr;
2353 :
2354 : // If we have a non-remote_elem neighbor link, then we can
2355 : // verify it.
2356 11132566 : if (elem && elem->neighbor_ptr(n) != remote_elem)
2357 : {
2358 11039926 : p_my_neighbor = &my_neighbor;
2359 11039926 : if (elem->neighbor_ptr(n))
2360 10449048 : my_neighbor = elem->neighbor_ptr(n)->id();
2361 :
2362 : // But wait - if we haven't set remote_elem links yet then
2363 : // some nullptr links on ghost elements might be
2364 : // future-remote_elem links, so we can't verify those.
2365 22086772 : if (!assert_valid_remote_elems &&
2366 11041552 : !elem->neighbor_ptr(n) &&
2367 1626 : elem->processor_id() != mesh.processor_id())
2368 759 : p_my_neighbor = nullptr;
2369 : }
2370 11132566 : libmesh_assert(mesh.comm().semiverify(p_my_neighbor));
2371 : }
2372 : }
2373 : }
2374 : #endif // DEBUG
2375 :
2376 :
2377 :
2378 : // Functors for correct_node_proc_ids
2379 : namespace {
2380 :
2381 : typedef std::unordered_map<dof_id_type, processor_id_type> proc_id_map_type;
2382 :
2383 : struct SyncNodeSet
2384 : {
2385 : typedef unsigned char datum; // bool but without bit twiddling issues
2386 :
2387 164 : SyncNodeSet(std::unordered_set<const Node *> & _set,
2388 8560 : MeshBase & _mesh) :
2389 8560 : node_set(_set), mesh(_mesh) {}
2390 :
2391 : std::unordered_set<const Node *> & node_set;
2392 :
2393 : MeshBase & mesh;
2394 :
2395 : // ------------------------------------------------------------
2396 38196 : void gather_data (const std::vector<dof_id_type> & ids,
2397 : std::vector<datum> & data)
2398 : {
2399 : // Find whether each requested node belongs in the set
2400 38196 : data.resize(ids.size());
2401 :
2402 1831968 : for (auto i : index_range(ids))
2403 : {
2404 1793772 : const dof_id_type id = ids[i];
2405 :
2406 : // We'd better have every node we're asked for
2407 1793772 : Node * node = mesh.node_ptr(id);
2408 :
2409 : // Return if the node is in the set.
2410 2890616 : data[i] = node_set.count(node);
2411 : }
2412 38196 : }
2413 :
2414 : // ------------------------------------------------------------
2415 38196 : bool act_on_data (const std::vector<dof_id_type> & ids,
2416 : const std::vector<datum> in_set)
2417 : {
2418 149 : bool data_changed = false;
2419 :
2420 : // Add nodes we've been informed of to our own set
2421 1831968 : for (auto i : index_range(ids))
2422 : {
2423 1793772 : if (in_set[i])
2424 : {
2425 1163913 : Node * node = mesh.node_ptr(ids[i]);
2426 1163913 : if (!node_set.count(node))
2427 : {
2428 256 : node_set.insert(node);
2429 256 : data_changed = true;
2430 : }
2431 : }
2432 : }
2433 :
2434 38196 : return data_changed;
2435 : }
2436 : };
2437 :
2438 :
2439 8232 : struct NodesNotInSet
2440 : {
2441 164 : NodesNotInSet(const std::unordered_set<const Node *> _set)
2442 8396 : : node_set(_set) {}
2443 :
2444 473444 : bool operator() (const Node * node) const
2445 : {
2446 946888 : if (node_set.count(node))
2447 280940 : return false;
2448 192504 : return true;
2449 : }
2450 :
2451 : const std::unordered_set<const Node *> node_set;
2452 : };
2453 :
2454 :
2455 : struct SyncProcIdsFromMap
2456 : {
2457 : typedef processor_id_type datum;
2458 :
2459 200 : SyncProcIdsFromMap(const proc_id_map_type & _map,
2460 36268 : MeshBase & _mesh) :
2461 36268 : new_proc_ids(_map), mesh(_mesh) {}
2462 :
2463 : const proc_id_map_type & new_proc_ids;
2464 :
2465 : MeshBase & mesh;
2466 :
2467 : // ------------------------------------------------------------
2468 126094 : void gather_data (const std::vector<dof_id_type> & ids,
2469 : std::vector<datum> & data)
2470 : {
2471 : // Find the new processor id of each requested node
2472 126094 : data.resize(ids.size());
2473 :
2474 7024650 : for (auto i : index_range(ids))
2475 : {
2476 6898556 : const dof_id_type id = ids[i];
2477 :
2478 : // Return the node's new processor id if it has one, or its
2479 : // old processor id if not.
2480 6898556 : if (const auto it = new_proc_ids.find(id);
2481 47012 : it != new_proc_ids.end())
2482 6267476 : data[i] = it->second;
2483 : else
2484 : {
2485 : // We'd better find every node we're asked for
2486 631080 : const Node & node = mesh.node_ref(id);
2487 631080 : data[i] = node.processor_id();
2488 : }
2489 : }
2490 126094 : }
2491 :
2492 : // ------------------------------------------------------------
2493 126094 : void act_on_data (const std::vector<dof_id_type> & ids,
2494 : const std::vector<datum> proc_ids)
2495 : {
2496 : // Set the node processor ids we've now been informed of
2497 7024650 : for (auto i : index_range(ids))
2498 : {
2499 6898556 : Node & node = mesh.node_ref(ids[i]);
2500 6898556 : node.processor_id() = proc_ids[i];
2501 : }
2502 126094 : }
2503 : };
2504 : }
2505 :
2506 :
2507 :
2508 36268 : void correct_node_proc_ids (MeshBase & mesh)
2509 : {
2510 400 : LOG_SCOPE("correct_node_proc_ids()","MeshTools");
2511 :
2512 : // This function must be run on all processors at once
2513 200 : libmesh_parallel_only(mesh.comm());
2514 :
2515 : // We require all processors to agree on nodal processor ids before
2516 : // going into this algorithm.
2517 : #ifdef DEBUG
2518 200 : libmesh_assert_parallel_consistent_procids<Node>(mesh);
2519 : #endif
2520 :
2521 : // If we have any unpartitioned elements at this
2522 : // stage there is a problem
2523 200 : libmesh_assert (n_elem(mesh.unpartitioned_elements_begin(),
2524 : mesh.unpartitioned_elements_end()) == 0);
2525 :
2526 : // Fix nodes' processor ids. Coarsening may have left us with nodes
2527 : // which are no longer touched by any elements of the same processor
2528 : // id, and for DofMap to work we need to fix that.
2529 :
2530 : // This is harder now that libMesh no longer requires a distributed
2531 : // mesh to ghost all nodal neighbors: it is possible for two active
2532 : // elements on two different processors to share the same node in
2533 : // such a way that neither processor knows the others' element
2534 : // exists!
2535 :
2536 : // While we're at it, if this mesh is configured to allow
2537 : // repartitioning, we'll repartition *all* the nodes' processor ids
2538 : // using the canonical Node heuristic, to try and improve DoF load
2539 : // balancing. But if the mesh is disallowing repartitioning, we
2540 : // won't touch processor_id on any node where it's valid, regardless
2541 : // of whether or not it's canonical.
2542 200 : bool repartition_all_nodes = !mesh.skip_noncritical_partitioning();
2543 400 : std::unordered_set<const Node *> valid_nodes;
2544 :
2545 : // If we aren't allowed to repartition, then we're going to leave
2546 : // every node we can at its current processor_id, and *only*
2547 : // repartition the nodes whose current processor id is incompatible
2548 : // with DoFMap (because it doesn't touch an active element, e.g. due
2549 : // to coarsening)
2550 36268 : if (!repartition_all_nodes)
2551 : {
2552 1129514 : for (const auto & elem : mesh.active_element_ptr_range())
2553 5815397 : for (const auto & node : elem->node_ref_range())
2554 5190562 : if (elem->processor_id() == node.processor_id())
2555 4820981 : valid_nodes.insert(&node);
2556 :
2557 164 : SyncNodeSet syncv(valid_nodes, mesh);
2558 :
2559 : Parallel::sync_dofobject_data_by_id
2560 16956 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), syncv);
2561 : }
2562 :
2563 : // We build up a set of compatible processor ids for each node
2564 400 : proc_id_map_type new_proc_ids;
2565 :
2566 8828594 : for (auto & elem : mesh.active_element_ptr_range())
2567 : {
2568 4416171 : processor_id_type pid = elem->processor_id();
2569 :
2570 42894840 : for (auto & node : elem->node_ref_range())
2571 : {
2572 38478669 : const dof_id_type id = node.id();
2573 38478669 : if (auto it = new_proc_ids.find(id);
2574 343100 : it == new_proc_ids.end())
2575 152704 : new_proc_ids.emplace(id, pid);
2576 : else
2577 27559949 : it->second = node.choose_processor_id(it->second, pid);
2578 : }
2579 35868 : }
2580 :
2581 : // Sort the new pids to push to each processor
2582 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, processor_id_type>>>
2583 400 : ids_to_push;
2584 :
2585 25013328 : for (const auto & node : mesh.node_ptr_range())
2586 12719688 : if (const auto it = std::as_const(new_proc_ids).find(node->id());
2587 12719688 : it != new_proc_ids.end() && node->processor_id() != DofObject::invalid_processor_id)
2588 10954588 : ids_to_push[node->processor_id()].emplace_back(node->id(), /*pid=*/it->second);
2589 :
2590 : auto action_functor =
2591 162336 : [& mesh, & new_proc_ids]
2592 : (processor_id_type,
2593 11377206 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
2594 : {
2595 11081804 : for (const auto & [id, pid] : data)
2596 : {
2597 10918720 : if (const auto it = new_proc_ids.find(id);
2598 152704 : it == new_proc_ids.end())
2599 0 : new_proc_ids.emplace(id, pid);
2600 : else
2601 : {
2602 10918720 : const Node & node = mesh.node_ref(id);
2603 10918720 : it->second = node.choose_processor_id(it->second, pid);
2604 : }
2605 : }
2606 36816 : };
2607 :
2608 : Parallel::push_parallel_vector_data
2609 36268 : (mesh.comm(), ids_to_push, action_functor);
2610 :
2611 : // Now new_proc_ids is correct for every node we used to own. Let's
2612 : // ask every other processor about the nodes they used to own. But
2613 : // first we'll need to keep track of which nodes we used to own,
2614 : // lest we get them confused with nodes we newly own.
2615 400 : std::unordered_set<Node *> ex_local_nodes;
2616 8888390 : for (auto & node : mesh.local_node_ptr_range())
2617 4543138 : if (const auto it = new_proc_ids.find(node->id());
2618 4543138 : it != new_proc_ids.end() && it->second != mesh.processor_id())
2619 35917 : ex_local_nodes.insert(node);
2620 :
2621 200 : SyncProcIdsFromMap sync(new_proc_ids, mesh);
2622 36268 : if (repartition_all_nodes)
2623 : Parallel::sync_dofobject_data_by_id
2624 55380 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync);
2625 : else
2626 : {
2627 164 : NodesNotInSet nnis(valid_nodes);
2628 :
2629 : Parallel::sync_dofobject_data_by_id
2630 16956 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), nnis, sync);
2631 : }
2632 :
2633 : // And finally let's update the nodes we used to own.
2634 40855 : for (const auto & node : ex_local_nodes)
2635 : {
2636 2103 : if (valid_nodes.count(node))
2637 2083 : continue;
2638 :
2639 2504 : const dof_id_type id = node->id();
2640 10 : const proc_id_map_type::iterator it = new_proc_ids.find(id);
2641 10 : libmesh_assert(it != new_proc_ids.end());
2642 2504 : node->processor_id() = it->second;
2643 : }
2644 :
2645 : // We should still have consistent nodal processor ids coming out of
2646 : // this algorithm, but if we're allowed to repartition the mesh then
2647 : // they should be canonically correct too.
2648 : #ifdef DEBUG
2649 200 : libmesh_assert_valid_procids<Node>(mesh);
2650 : //if (repartition_all_nodes)
2651 : // libmesh_assert_canonical_node_procids(mesh);
2652 : #endif
2653 36268 : }
2654 :
2655 :
2656 :
2657 19630 : void Private::globally_renumber_nodes_and_elements (MeshBase & mesh)
2658 : {
2659 19630 : MeshCommunication().assign_global_indices(mesh);
2660 19630 : }
2661 :
2662 : } // namespace MeshTools
2663 :
2664 : } // namespace libMesh
|
