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 and subroutines
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 2583974 : FindBBox () : _bbox()
108 36822 : {}
109 :
110 870 : FindBBox (FindBBox & other, Threads::split) :
111 870 : _bbox(other._bbox)
112 298 : {}
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 2600856 : void operator()(const ConstElemRange & range)
124 : {
125 28868534 : for (const auto & elem : range)
126 : {
127 1413140 : libmesh_assert(elem);
128 26267678 : _bbox.union_with(elem->loose_bounding_box());
129 : }
130 2600856 : }
131 :
132 18694 : Point & min() { return _bbox.min(); }
133 :
134 18694 : 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 298 : void join (const FindBBox & other)
140 : {
141 870 : _bbox.union_with(other._bbox);
142 580 : }
143 : #endif
144 :
145 54950 : libMesh::BoundingBox & bbox ()
146 : {
147 54950 : 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 9769010 : for (unsigned int v = 0; v != n_vars[s]; ++v, ++i)
182 : {
183 6863568 : n_comp[i] = d->n_comp(s,v);
184 6863568 : 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 16890522 : 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 16890522 : if (d)
218 : {
219 16223082 : tempmin = tempmax = d->unique_id();
220 : }
221 : else
222 : {
223 667440 : TIMPI::Attributes<unique_id_type>::set_highest(tempmin);
224 667440 : TIMPI::Attributes<unique_id_type>::set_lowest(tempmax);
225 : }
226 16890522 : comm.min(tempmin);
227 16890522 : comm.max(tempmax);
228 33113604 : bool invalid = d && ((d->unique_id() != tempmin) ||
229 16223082 : (d->unique_id() != tempmax));
230 16890522 : comm.max(invalid);
231 :
232 : // First verify that everything is in sync
233 16890522 : libmesh_assert(!invalid);
234 :
235 : // Then verify that any remote id doesn't duplicate a local one.
236 16890522 : if (!d && tempmin == tempmax)
237 81256 : libmesh_assert(!unique_ids.count(tempmin));
238 16890522 : }
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 14138 : std::unique_ptr<MeshBase> reprepared_mesh_clone (const MeshBase & mesh)
412 : {
413 14138 : const MeshBase::Preparation prep = mesh.preparation();
414 :
415 : // If the mesh thinks it's prepared in some way, *re*-preparing in
416 : // that way shouldn't change a clone of it, as long as we disallow
417 : // repartitioning or renumbering or remote element removal.
418 14138 : std::unique_ptr<MeshBase> mesh_clone = mesh.clone();
419 :
420 12466 : const bool old_allow_renumbering = mesh_clone->allow_renumbering();
421 : const bool old_allow_remote_element_removal =
422 12466 : mesh_clone->allow_remote_element_removal();
423 12466 : const bool old_skip_partitioning = mesh_clone->skip_partitioning();
424 12416 : mesh_clone->allow_renumbering(false);
425 12416 : mesh_clone->allow_remote_element_removal(false);
426 12416 : mesh_clone->skip_partitioning(true);
427 :
428 : // If the mesh thinks it's already completely prepared, test that
429 14138 : if (prep)
430 14138 : mesh_clone->prepare_for_use();
431 : // If the mesh thinks it's somewhat prepared, test each way it
432 : // thinks so.
433 : else
434 : {
435 0 : if (prep.has_synched_id_counts)
436 0 : mesh_clone->update_parallel_id_counts();
437 :
438 0 : if (prep.has_neighbor_ptrs)
439 0 : mesh_clone->find_neighbors();
440 :
441 0 : if (prep.has_cached_elem_data)
442 0 : mesh_clone->cache_elem_data();
443 :
444 0 : if (mesh.allow_detect_interior_parents() &&
445 0 : prep.has_interior_parent_ptrs)
446 0 : mesh_clone->detect_interior_parents();
447 :
448 0 : if (old_allow_remote_element_removal &&
449 0 : prep.has_removed_remote_elements)
450 0 : mesh_clone->delete_remote_elements();
451 :
452 0 : if (prep.has_removed_orphaned_nodes)
453 0 : mesh_clone->remove_orphaned_nodes();
454 :
455 0 : if (prep.has_boundary_id_sets)
456 0 : mesh_clone->get_boundary_info().regenerate_id_sets();
457 :
458 : // I don't know how we'll tell if this changes anything, but
459 : // we'll do it for completeness
460 0 : if (prep.has_reinit_ghosting_functors)
461 0 : mesh_clone->reinit_ghosting_functors();
462 : }
463 :
464 : // Restore original flag values
465 12416 : mesh_clone->allow_renumbering(old_allow_renumbering);
466 12416 : mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal);
467 12416 : mesh_clone->skip_partitioning(old_skip_partitioning);
468 :
469 26554 : return mesh_clone;
470 0 : }
471 :
472 :
473 : }
474 :
475 :
476 : namespace libMesh
477 : {
478 :
479 : // ------------------------------------------------------------
480 : // MeshTools functions
481 :
482 : namespace MeshTools
483 : {
484 :
485 0 : dof_id_type total_weight(const MeshBase & mesh)
486 : {
487 0 : if (!mesh.is_serial())
488 : {
489 0 : libmesh_parallel_only(mesh.comm());
490 0 : dof_id_type weight = MeshTools::weight (mesh, mesh.processor_id());
491 0 : mesh.comm().sum(weight);
492 : dof_id_type unpartitioned_weight =
493 0 : MeshTools::weight (mesh, DofObject::invalid_processor_id);
494 0 : return weight + unpartitioned_weight;
495 : }
496 :
497 0 : SumElemWeight sew;
498 :
499 0 : Threads::parallel_reduce (ConstElemRange (mesh.elements_begin(),
500 0 : mesh.elements_end()),
501 : sew);
502 0 : return sew.weight();
503 :
504 : }
505 :
506 :
507 :
508 0 : dof_id_type weight(const MeshBase & mesh, const processor_id_type pid)
509 : {
510 0 : SumElemWeight sew;
511 :
512 0 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
513 0 : mesh.pid_elements_end(pid)),
514 : sew);
515 0 : return sew.weight();
516 : }
517 :
518 :
519 :
520 6289 : void build_nodes_to_elem_map (const MeshBase & mesh,
521 : std::vector<std::vector<dof_id_type>> & nodes_to_elem_map)
522 : {
523 : // A vector indexed over all nodes is too inefficient to use for a
524 : // distributed mesh. Use the unordered_map API instead.
525 6289 : if (!mesh.is_serial())
526 : libmesh_deprecated();
527 :
528 6289 : nodes_to_elem_map.resize (mesh.max_node_id());
529 :
530 159220 : for (const auto & elem : mesh.element_ptr_range())
531 308664 : for (auto & node : elem->node_ref_range())
532 : {
533 12504 : libmesh_assert_less (node.id(), nodes_to_elem_map.size());
534 12504 : libmesh_assert_less (elem->id(), mesh.n_elem());
535 :
536 236268 : nodes_to_elem_map[node.id()].push_back(elem->id());
537 5773 : }
538 6289 : }
539 :
540 :
541 :
542 0 : void build_nodes_to_elem_map (const MeshBase & mesh,
543 : std::vector<std::vector<const Elem *>> & nodes_to_elem_map)
544 : {
545 : // A vector indexed over all nodes is too inefficient to use for a
546 : // distributed mesh. Use the unordered_map API instead.
547 0 : if (!mesh.is_serial())
548 : libmesh_deprecated();
549 :
550 0 : nodes_to_elem_map.resize (mesh.max_node_id());
551 :
552 0 : for (const auto & elem : mesh.element_ptr_range())
553 0 : for (auto & node : elem->node_ref_range())
554 : {
555 0 : libmesh_assert_less (node.id(), nodes_to_elem_map.size());
556 :
557 0 : nodes_to_elem_map[node.id()].push_back(elem);
558 0 : }
559 0 : }
560 :
561 :
562 :
563 256003 : void build_nodes_to_elem_map (const MeshBase & mesh,
564 : std::unordered_map<dof_id_type, std::vector<dof_id_type>> & nodes_to_elem_map)
565 : {
566 7584 : nodes_to_elem_map.clear();
567 :
568 62694564 : for (const auto & elem : mesh.element_ptr_range())
569 190598307 : for (auto & node : elem->node_ref_range())
570 156500029 : nodes_to_elem_map[node.id()].push_back(elem->id());
571 256003 : }
572 :
573 :
574 :
575 5823 : void build_nodes_to_elem_map (const MeshBase & mesh,
576 : std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map)
577 : {
578 164 : nodes_to_elem_map.clear();
579 :
580 769248 : for (const auto & elem : mesh.element_ptr_range())
581 2735901 : for (auto & node : elem->node_ref_range())
582 2328997 : nodes_to_elem_map[node.id()].push_back(elem);
583 5823 : }
584 :
585 :
586 :
587 : std::unordered_set<dof_id_type>
588 6119 : find_boundary_nodes(const MeshBase & mesh)
589 : {
590 178 : std::unordered_set<dof_id_type> boundary_nodes;
591 :
592 : // Loop over elements, find those on boundary, and
593 : // mark them as true in on_boundary.
594 1088800 : for (const auto & elem : mesh.active_element_ptr_range())
595 3109744 : for (auto s : elem->side_index_range())
596 2636242 : if (elem->neighbor_ptr(s) == nullptr) // on the boundary
597 : {
598 143276 : auto nodes_on_side = elem->nodes_on_side(s);
599 :
600 742200 : for (auto & local_id : nodes_on_side)
601 634704 : boundary_nodes.insert(elem->node_ptr(local_id)->id());
602 5763 : }
603 :
604 6119 : return boundary_nodes;
605 : }
606 :
607 : std::unordered_set<dof_id_type>
608 3137 : find_block_boundary_nodes(const MeshBase & mesh)
609 : {
610 94 : std::unordered_set<dof_id_type> block_boundary_nodes;
611 :
612 : // Loop over elements, find those on boundary, and
613 : // mark them as true in on_boundary.
614 864512 : for (const auto & elem : mesh.active_element_ptr_range())
615 2469648 : for (auto s : elem->side_index_range())
616 2090882 : if (elem->neighbor_ptr(s) && (elem->neighbor_ptr(s)->subdomain_id() != elem->subdomain_id()))
617 : {
618 0 : auto nodes_on_side = elem->nodes_on_side(s);
619 :
620 0 : for (auto & local_id : nodes_on_side)
621 0 : block_boundary_nodes.insert(elem->node_ptr(local_id)->id());
622 2949 : }
623 :
624 3137 : return block_boundary_nodes;
625 : }
626 :
627 :
628 :
629 : libMesh::BoundingBox
630 1300299 : create_bounding_box (const MeshBase & mesh)
631 : {
632 : // This function must be run on all processors at once
633 18128 : libmesh_parallel_only(mesh.comm());
634 :
635 18128 : FindBBox find_bbox;
636 :
637 : // Start with any unpartitioned elements we know about locally
638 2600598 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(DofObject::invalid_processor_id),
639 1318427 : mesh.pid_elements_end(DofObject::invalid_processor_id)),
640 : find_bbox);
641 :
642 : // And combine with our local elements
643 1300299 : find_bbox.bbox().union_with(create_local_bounding_box(mesh));
644 :
645 : // Compare the bounding boxes across processors
646 1300299 : mesh.comm().min(find_bbox.min());
647 1300299 : mesh.comm().max(find_bbox.max());
648 :
649 1300299 : return find_bbox.bbox();
650 : }
651 :
652 :
653 :
654 : libMesh::BoundingBox
655 19630 : create_nodal_bounding_box (const MeshBase & mesh)
656 : {
657 : // This function must be run on all processors at once
658 566 : libmesh_parallel_only(mesh.comm());
659 :
660 566 : FindBBox find_bbox;
661 :
662 : // Start with any unpartitioned nodes we know about locally
663 39260 : Threads::parallel_reduce (ConstNodeRange (mesh.pid_nodes_begin(DofObject::invalid_processor_id),
664 39260 : mesh.pid_nodes_end(DofObject::invalid_processor_id)),
665 : find_bbox);
666 :
667 : // Add our local nodes
668 39260 : Threads::parallel_reduce (ConstNodeRange (mesh.local_nodes_begin(),
669 38694 : mesh.local_nodes_end()),
670 : find_bbox);
671 :
672 : // Compare the bounding boxes across processors
673 19630 : mesh.comm().min(find_bbox.min());
674 19630 : mesh.comm().max(find_bbox.max());
675 :
676 19630 : return find_bbox.bbox();
677 : }
678 :
679 :
680 :
681 : Sphere
682 0 : bounding_sphere(const MeshBase & mesh)
683 : {
684 0 : libMesh::BoundingBox bbox = create_bounding_box(mesh);
685 :
686 0 : const Real diag = (bbox.second - bbox.first).norm();
687 0 : const Point cent = (bbox.second + bbox.first)/2;
688 :
689 0 : return Sphere (cent, .5*diag);
690 : }
691 :
692 :
693 :
694 : libMesh::BoundingBox
695 1300299 : create_local_bounding_box (const MeshBase & mesh)
696 : {
697 18128 : FindBBox find_bbox;
698 :
699 2600598 : Threads::parallel_reduce (ConstElemRange (mesh.local_elements_begin(),
700 1318427 : mesh.local_elements_end()),
701 : find_bbox);
702 :
703 1300299 : return find_bbox.bbox();
704 : }
705 :
706 :
707 :
708 : libMesh::BoundingBox
709 0 : create_processor_bounding_box (const MeshBase & mesh,
710 : const processor_id_type pid)
711 : {
712 : // This can only be run in parallel, with consistent arguments.
713 0 : libmesh_parallel_only(mesh.comm());
714 0 : libmesh_assert(mesh.comm().verify(pid));
715 :
716 0 : libmesh_assert_less (pid, mesh.n_processors());
717 :
718 0 : FindBBox find_bbox;
719 :
720 0 : Threads::parallel_reduce (ConstElemRange (mesh.pid_elements_begin(pid),
721 0 : mesh.pid_elements_end(pid)),
722 : find_bbox);
723 :
724 : // Compare the bounding boxes across processors
725 0 : mesh.comm().min(find_bbox.min());
726 0 : mesh.comm().max(find_bbox.max());
727 :
728 0 : return find_bbox.bbox();
729 : }
730 :
731 :
732 :
733 : Sphere
734 0 : processor_bounding_sphere (const MeshBase & mesh,
735 : const processor_id_type pid)
736 : {
737 : libMesh::BoundingBox bbox =
738 0 : create_processor_bounding_box(mesh, pid);
739 :
740 0 : const Real diag = (bbox.second - bbox.first).norm();
741 0 : const Point cent = (bbox.second + bbox.first)/2;
742 :
743 0 : return Sphere (cent, .5*diag);
744 : }
745 :
746 :
747 :
748 : libMesh::BoundingBox
749 0 : create_subdomain_bounding_box (const MeshBase & mesh,
750 : const subdomain_id_type sid)
751 : {
752 : // This can only be run in parallel, with consistent arguments.
753 0 : libmesh_parallel_only(mesh.comm());
754 0 : libmesh_assert(mesh.comm().verify(sid));
755 :
756 0 : FindBBox find_bbox;
757 :
758 : Threads::parallel_reduce
759 0 : (ConstElemRange (mesh.active_local_subdomain_elements_begin(sid),
760 0 : mesh.active_local_subdomain_elements_end(sid)),
761 : find_bbox);
762 :
763 : // Compare the bounding boxes across processors
764 0 : mesh.comm().min(find_bbox.min());
765 0 : mesh.comm().max(find_bbox.max());
766 :
767 0 : return find_bbox.bbox();
768 : }
769 :
770 :
771 :
772 : Sphere
773 0 : subdomain_bounding_sphere (const MeshBase & mesh,
774 : const subdomain_id_type sid)
775 : {
776 : libMesh::BoundingBox bbox =
777 0 : create_subdomain_bounding_box(mesh, sid);
778 :
779 0 : const Real diag = (bbox.second - bbox.first).norm();
780 0 : const Point cent = (bbox.second + bbox.first)/2;
781 :
782 0 : return Sphere (cent, .5*diag);
783 : }
784 :
785 :
786 :
787 0 : void elem_types (const MeshBase & mesh,
788 : std::vector<ElemType> & et)
789 : {
790 : // Loop over the the elements. If the current element type isn't in
791 : // the vector, insert it.
792 0 : for (const auto & elem : mesh.element_ptr_range())
793 0 : if (!std::count(et.begin(), et.end(), elem->type()))
794 0 : et.push_back(elem->type());
795 0 : }
796 :
797 :
798 :
799 0 : dof_id_type n_elem_of_type (const MeshBase & mesh,
800 : const ElemType type)
801 : {
802 0 : return static_cast<dof_id_type>(std::distance(mesh.type_elements_begin(type),
803 0 : mesh.type_elements_end (type)));
804 : }
805 :
806 :
807 :
808 0 : dof_id_type n_active_elem_of_type (const MeshBase & mesh,
809 : const ElemType type)
810 : {
811 0 : return static_cast<dof_id_type>(std::distance(mesh.active_type_elements_begin(type),
812 0 : mesh.active_type_elements_end (type)));
813 : }
814 :
815 0 : dof_id_type n_non_subactive_elem_of_type_at_level(const MeshBase & mesh,
816 : const ElemType type,
817 : const unsigned int level)
818 : {
819 0 : dof_id_type cnt = 0;
820 :
821 : // iterate over the elements of the specified type
822 0 : for (const auto & elem : as_range(mesh.type_elements_begin(type),
823 0 : mesh.type_elements_end(type)))
824 0 : if ((elem->level() == level) && !elem->subactive())
825 0 : cnt++;
826 :
827 0 : return cnt;
828 : }
829 :
830 :
831 5471 : unsigned int n_active_local_levels(const MeshBase & mesh)
832 : {
833 : struct LevelCounter {
834 : unsigned int nl;
835 :
836 5471 : LevelCounter () : nl(0) {}
837 :
838 0 : LevelCounter (LevelCounter &, Threads::split) :
839 0 : nl(0) {}
840 :
841 5471 : void operator()(const ConstElemRange & range) {
842 70636 : for (const Elem * elem : range)
843 69004 : nl = std::max(elem->level() + 1, nl);
844 5471 : }
845 :
846 0 : void join(const LevelCounter & other) {
847 0 : nl = std::max(nl, other.nl);
848 0 : }
849 : };
850 :
851 156 : LevelCounter counter;
852 :
853 5471 : Threads::parallel_reduce(mesh.active_local_element_stored_range(), counter);
854 :
855 5471 : return counter.nl;
856 : }
857 :
858 :
859 :
860 5471 : unsigned int n_active_levels(const MeshBase & mesh)
861 : {
862 156 : libmesh_parallel_only(mesh.comm());
863 :
864 5471 : unsigned int nl = n_active_local_levels(mesh);
865 :
866 10786 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
867 21572 : mesh.unpartitioned_elements_end()))
868 0 : if (elem->active())
869 5159 : nl = std::max(elem->level() + 1, nl);
870 :
871 5471 : mesh.comm().max(nl);
872 5471 : return nl;
873 : }
874 :
875 :
876 :
877 1889446 : unsigned int n_local_levels(const MeshBase & mesh)
878 : {
879 1889446 : unsigned int nl = 0;
880 :
881 6081360 : for (const auto & elem : as_range(mesh.local_elements_begin(),
882 41222961 : mesh.local_elements_end()))
883 37727250 : nl = std::max(elem->level() + 1, nl);
884 :
885 1889446 : return nl;
886 : }
887 :
888 :
889 :
890 1889446 : unsigned int n_levels(const MeshBase & mesh)
891 : {
892 43300 : libmesh_parallel_only(mesh.comm());
893 :
894 1889446 : unsigned int nl = n_local_levels(mesh);
895 :
896 4120560 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
897 21533345 : mesh.unpartitioned_elements_end()))
898 15877541 : nl = std::max(elem->level() + 1, nl);
899 :
900 1889446 : mesh.comm().max(nl);
901 :
902 : // n_levels() is only valid and should only be called in cases where
903 : // the mesh is validly distributed (or serialized). Let's run an
904 : // expensive test in debug mode to make sure this is such a case.
905 : #ifdef DEBUG
906 43300 : const unsigned int paranoid_nl = paranoid_n_levels(mesh);
907 43300 : libmesh_assert_equal_to(nl, paranoid_nl);
908 : #endif
909 1889446 : return nl;
910 : }
911 :
912 :
913 :
914 57479 : unsigned int paranoid_n_levels(const MeshBase & mesh)
915 : {
916 43716 : libmesh_parallel_only(mesh.comm());
917 :
918 57479 : unsigned int nl = 0;
919 5565130 : for (const auto & elem : mesh.element_ptr_range())
920 5507235 : nl = std::max(elem->level() + 1, nl);
921 :
922 57479 : mesh.comm().max(nl);
923 57479 : return nl;
924 : }
925 :
926 :
927 :
928 639 : dof_id_type n_connected_components(const MeshBase & mesh,
929 : Real constraint_tol)
930 : {
931 36 : LOG_SCOPE("n_connected_components()", "MeshTools");
932 :
933 : // Yes, I'm being lazy. This is for mesh analysis before a
934 : // simulation, not anything going in any loops.
935 639 : if (!mesh.is_serial_on_zero())
936 0 : libmesh_not_implemented();
937 :
938 639 : dof_id_type n_components = 0;
939 :
940 657 : if (mesh.processor_id())
941 : {
942 522 : mesh.comm().broadcast(n_components);
943 531 : return n_components;
944 : }
945 :
946 : // All nodes in a set here are connected (at least indirectly) to
947 : // all other nodes in the same set, but have not yet been discovered
948 : // to be connected to nodes in other sets.
949 : //
950 : // Using an unordered_set of ids rather than a set of pointers seems
951 : // to be roughly 150% faster?
952 : // typedef const Node * node_entry_type
953 : typedef dof_id_type node_entry_type;
954 9 : std::vector<std::unordered_set<node_entry_type>> components;
955 :
956 : // With a typical mesh with few components and somewhat-contiguous
957 : // ordering, vector performance should be fine. With a mesh with
958 : // many components or completely scrambled ordering, performance
959 : // can be a disaster.
960 2742 : auto find_component = [&components](node_entry_type n) {
961 6158 : for (auto & c: components)
962 4064 : if (c.find(n) != c.end())
963 54 : return &c;
964 :
965 162 : return (std::unordered_set<node_entry_type> *)(nullptr);
966 108 : };
967 :
968 : auto add_to_component =
969 2010 : [&find_component]
970 2202 : (std::unordered_set<node_entry_type> & component, node_entry_type n) {
971 : // We may already be in the desired component
972 2412 : if (component.find(n) != component.end())
973 51 : return;
974 :
975 1950 : auto current_component = find_component(n);
976 :
977 : // Didn't we *just* check this?
978 150 : libmesh_assert (&component != current_component);
979 :
980 : // If we're unknown, we should be in the desired component
981 1950 : if (!current_component)
982 147 : component.insert(n);
983 :
984 : // If we think we're in another component, it should actually be
985 : // part of the desired component
986 : else
987 : {
988 : // Merge the component likely to be smaller into the one
989 : // likely to be larger - this is orders of magnitude faster
990 : // than the other way around!
991 3 : current_component->merge(component);
992 3 : current_component->swap(component);
993 3 : libmesh_assert(current_component->empty());
994 : }
995 108 : };
996 :
997 9 : auto & constraint_rows = mesh.get_constraint_rows();
998 :
999 1659 : for (const auto & elem : mesh.element_ptr_range())
1000 : {
1001 : // const node_entry_type first_node = elem->node_ptr(0);
1002 792 : const node_entry_type first_node = elem->node_id(0);
1003 :
1004 792 : auto component = find_component(first_node);
1005 :
1006 : // If we didn't find one, make a new one, reusing an existing
1007 : // slot if possible or growing our vector if necessary
1008 792 : if (!component)
1009 684 : for (auto & c: components)
1010 354 : if (c.empty())
1011 0 : component = &c;
1012 :
1013 432 : if (!component)
1014 219 : component = &components.emplace_back();
1015 :
1016 3234 : for (const Node & node : elem->node_ref_range())
1017 : {
1018 : // const node_entry_type n = &node;
1019 396 : const node_entry_type n = node.id();
1020 2376 : add_to_component(*component, n);
1021 :
1022 2376 : auto it = constraint_rows.find(&node);
1023 2376 : if (it == constraint_rows.end())
1024 195 : continue;
1025 :
1026 72 : for (const auto & [pr, val] : it->second)
1027 : {
1028 : // Ignore too-trivial constraint coefficients if
1029 : // we get a non-default-0 constraint_tol
1030 39 : if (std::abs(val) < constraint_tol)
1031 0 : continue;
1032 :
1033 36 : const Elem * spline_elem = pr.first;
1034 3 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1035 :
1036 : const Node * spline_node =
1037 36 : spline_elem->node_ptr(pr.second);
1038 :
1039 : // add_to_component(*component, spline_node);
1040 36 : add_to_component(*component, spline_node->id());
1041 : }
1042 : }
1043 90 : }
1044 :
1045 327 : for (auto & component : components)
1046 219 : if (!component.empty())
1047 144 : ++n_components;
1048 :
1049 : // We calculated this on proc 0; now let everyone else know too
1050 108 : mesh.comm().broadcast(n_components);
1051 :
1052 108 : return n_components;
1053 90 : }
1054 :
1055 :
1056 :
1057 0 : void get_not_subactive_node_ids(const MeshBase & mesh,
1058 : std::set<dof_id_type> & not_subactive_node_ids)
1059 : {
1060 0 : for (const auto & elem : mesh.element_ptr_range())
1061 0 : if (!elem->subactive())
1062 0 : for (auto & n : elem->node_ref_range())
1063 0 : not_subactive_node_ids.insert(n.id());
1064 0 : }
1065 :
1066 :
1067 :
1068 547143 : dof_id_type n_elem (const MeshBase::const_element_iterator & begin,
1069 : const MeshBase::const_element_iterator & end)
1070 : {
1071 1072472 : return cast_int<dof_id_type>(std::distance(begin, end));
1072 : }
1073 :
1074 :
1075 :
1076 0 : dof_id_type n_nodes (const MeshBase::const_node_iterator & begin,
1077 : const MeshBase::const_node_iterator & end)
1078 : {
1079 0 : return cast_int<dof_id_type>(std::distance(begin, end));
1080 : }
1081 :
1082 :
1083 :
1084 4050 : Real volume (const MeshBase & mesh,
1085 : unsigned int dim)
1086 : {
1087 114 : libmesh_parallel_only(mesh.comm());
1088 :
1089 4050 : if (dim == libMesh::invalid_uint)
1090 4050 : dim = mesh.mesh_dimension();
1091 :
1092 4050 : Real vol = 0;
1093 :
1094 : // first my local elements
1095 9733 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1096 56215 : mesh.local_elements_end()))
1097 20995 : if (elem->dim() == dim)
1098 24817 : vol += elem->volume();
1099 :
1100 : // then count any unpartitioned objects, once
1101 4164 : if (mesh.processor_id() == 0)
1102 1317 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1103 2634 : mesh.unpartitioned_elements_end()))
1104 0 : if (elem->dim() == dim)
1105 573 : vol += elem->volume();
1106 :
1107 4050 : mesh.comm().sum(vol);
1108 4050 : return vol;
1109 : }
1110 :
1111 :
1112 :
1113 5471 : unsigned int n_p_levels (const MeshBase & mesh)
1114 : {
1115 156 : libmesh_parallel_only(mesh.comm());
1116 :
1117 5471 : unsigned int max_p_level = 0;
1118 :
1119 : // first my local elements
1120 17480 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1121 93378 : mesh.local_elements_end()))
1122 83659 : max_p_level = std::max(elem->p_level(), max_p_level);
1123 :
1124 : // then any unpartitioned objects
1125 10786 : for (const auto & elem : as_range(mesh.unpartitioned_elements_begin(),
1126 21572 : mesh.unpartitioned_elements_end()))
1127 5159 : max_p_level = std::max(elem->p_level(), max_p_level);
1128 :
1129 5471 : mesh.comm().max(max_p_level);
1130 5471 : return max_p_level + 1;
1131 : }
1132 :
1133 :
1134 :
1135 0 : void find_nodal_neighbors(const MeshBase &,
1136 : const Node & node,
1137 : const std::vector<std::vector<const Elem *>> & nodes_to_elem_map,
1138 : std::vector<const Node *> & neighbors)
1139 : {
1140 0 : find_nodal_neighbors_helper(node.id(), nodes_to_elem_map[node.id()],
1141 : neighbors);
1142 0 : }
1143 :
1144 :
1145 :
1146 225994 : void find_nodal_neighbors(const MeshBase &,
1147 : const Node & node,
1148 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1149 : std::vector<const Node *> & neighbors)
1150 : {
1151 : const std::vector<const Elem *> node_to_elem_vec =
1152 234012 : libmesh_map_find(nodes_to_elem_map, node.id());
1153 225994 : find_nodal_neighbors_helper(node.id(), node_to_elem_vec, neighbors);
1154 225994 : }
1155 :
1156 49559 : void find_nodal_or_face_neighbors(
1157 : const MeshBase & mesh,
1158 : const Node & node,
1159 : const std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem_map,
1160 : std::vector<const Node *> & neighbors)
1161 : {
1162 : // Find all the nodal neighbors... that is the nodes directly connected
1163 : // to this node through one edge.
1164 49559 : find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbors);
1165 :
1166 : // If no neighbors are found, use all nodes on the containing side as
1167 : // neighbors.
1168 49559 : if (!neighbors.size())
1169 : {
1170 : // Grab the element containing node
1171 6804 : const auto * elem = libmesh_map_find(nodes_to_elem_map, node.id()).front();
1172 : // Find the element side containing node
1173 17503 : for (const auto &side : elem->side_index_range())
1174 : {
1175 17503 : const auto &nodes_on_side = elem->nodes_on_side(side);
1176 : const auto it =
1177 23074 : std::find_if(nodes_on_side.begin(), nodes_on_side.end(), [&](auto local_node_id) {
1178 64593 : return elem->node_id(local_node_id) == node.id();
1179 1946 : });
1180 :
1181 18476 : if (it != nodes_on_side.end())
1182 : {
1183 58008 : for (const auto &local_node_id : nodes_on_side)
1184 : // No need to add node itself as a neighbor
1185 53931 : if (const auto *node_ptr = elem->node_ptr(local_node_id);
1186 2727 : *node_ptr != node)
1187 44400 : neighbors.push_back(node_ptr);
1188 347 : break;
1189 : }
1190 : }
1191 : }
1192 3048 : libmesh_assert(neighbors.size());
1193 49559 : }
1194 :
1195 :
1196 :
1197 0 : void find_hanging_nodes_and_parents(const MeshBase & mesh,
1198 : std::map<dof_id_type, std::vector<dof_id_type>> & hanging_nodes)
1199 : {
1200 : // Loop through all the elements
1201 0 : for (auto & elem : mesh.active_local_element_ptr_range())
1202 0 : if (elem->type() == QUAD4)
1203 0 : for (auto s : elem->side_index_range())
1204 : {
1205 : // Loop over the sides looking for sides that have hanging nodes
1206 : // This code is inspired by compute_proj_constraints()
1207 0 : const Elem * neigh = elem->neighbor_ptr(s);
1208 :
1209 : // If not a boundary side
1210 0 : if (neigh != nullptr)
1211 : {
1212 : // Is there a coarser element next to this one?
1213 0 : if (neigh->level() < elem->level())
1214 : {
1215 0 : const Elem * ancestor = elem;
1216 0 : while (neigh->level() < ancestor->level())
1217 0 : ancestor = ancestor->parent();
1218 0 : unsigned int s_neigh = neigh->which_neighbor_am_i(ancestor);
1219 0 : libmesh_assert_less (s_neigh, neigh->n_neighbors());
1220 :
1221 : // Couple of helper uints...
1222 0 : unsigned int local_node1=0;
1223 0 : unsigned int local_node2=0;
1224 :
1225 0 : bool found_in_neighbor = false;
1226 :
1227 : // Find the two vertices that make up this side
1228 0 : while (!elem->is_node_on_side(local_node1++,s)) { }
1229 0 : local_node1--;
1230 :
1231 : // Start looking for the second one with the next node
1232 0 : local_node2=local_node1+1;
1233 :
1234 : // Find the other one
1235 0 : while (!elem->is_node_on_side(local_node2++,s)) { }
1236 0 : local_node2--;
1237 :
1238 : //Pull out their global ids:
1239 0 : dof_id_type node1 = elem->node_id(local_node1);
1240 0 : dof_id_type node2 = elem->node_id(local_node2);
1241 :
1242 : // Now find which node is present in the neighbor
1243 : // FIXME This assumes a level one rule!
1244 : // The _other_ one is the hanging node
1245 :
1246 : // First look for the first one
1247 : // FIXME could be streamlined a bit
1248 0 : for (unsigned int n=0;n<neigh->n_sides();n++)
1249 0 : if (neigh->node_id(n) == node1)
1250 0 : found_in_neighbor=true;
1251 :
1252 0 : dof_id_type hanging_node=0;
1253 :
1254 0 : if (!found_in_neighbor)
1255 0 : hanging_node=node1;
1256 : else // If it wasn't node1 then it must be node2!
1257 0 : hanging_node=node2;
1258 :
1259 : // Reset for reuse
1260 0 : local_node1=0;
1261 :
1262 : // Find the first node that makes up the side in the neighbor (these should be the parent nodes)
1263 0 : while (!neigh->is_node_on_side(local_node1++,s_neigh)) { }
1264 0 : local_node1--;
1265 :
1266 0 : local_node2=local_node1+1;
1267 :
1268 : // Find the second node...
1269 0 : while (!neigh->is_node_on_side(local_node2++,s_neigh)) { }
1270 0 : local_node2--;
1271 :
1272 : // Save them if we haven't already found the parents for this one
1273 0 : if (hanging_nodes[hanging_node].size()<2)
1274 : {
1275 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node1));
1276 0 : hanging_nodes[hanging_node].push_back(neigh->node_id(local_node2));
1277 : }
1278 : }
1279 : }
1280 0 : }
1281 0 : }
1282 :
1283 :
1284 :
1285 36 : void clear_spline_nodes(MeshBase & mesh)
1286 : {
1287 6 : std::vector<Elem *> nodeelem_to_delete;
1288 :
1289 8737 : for (auto & elem : mesh.element_ptr_range())
1290 5073 : if (elem->type() == NODEELEM &&
1291 4140 : elem->mapping_type() == RATIONAL_BERNSTEIN_MAP)
1292 4170 : nodeelem_to_delete.push_back(elem);
1293 :
1294 3 : auto & constraint_rows = mesh.get_constraint_rows();
1295 :
1296 : // All our constraint_rows ought to be for spline constraints we're
1297 : // about to get rid of.
1298 : #ifndef NDEBUG
1299 762 : for (auto & node_row : constraint_rows)
1300 2064 : for (auto pr : node_row.second)
1301 : {
1302 1305 : const Elem * elem = pr.first.first;
1303 1305 : libmesh_assert(elem->type() == NODEELEM);
1304 1305 : libmesh_assert(elem->mapping_type() == RATIONAL_BERNSTEIN_MAP);
1305 : }
1306 : #endif
1307 :
1308 3 : constraint_rows.clear();
1309 :
1310 4176 : for (Elem * elem : nodeelem_to_delete)
1311 : {
1312 690 : Node * node = elem->node_ptr(0);
1313 4140 : mesh.delete_elem(elem);
1314 4140 : mesh.delete_node(node);
1315 : }
1316 36 : }
1317 :
1318 :
1319 :
1320 1772 : bool valid_is_prepared (const MeshBase & mesh)
1321 : {
1322 100 : LOG_SCOPE("valid_is_prepared()", "MeshTools");
1323 :
1324 1772 : const MeshBase::Preparation prep = mesh.preparation();
1325 :
1326 : // If the mesh doesn't think *anything* has been prepared, we have
1327 : // nothing to check.
1328 1772 : if (prep == MeshBase::Preparation())
1329 0 : return true;
1330 :
1331 : // If the mesh thinks it's partitioned, check. These are counts,
1332 : // not caches, so it's a real check.
1333 1772 : if (prep.is_partitioned)
1334 3544 : if (mesh.n_unpartitioned_elem() ||
1335 50 : mesh.n_unpartitioned_nodes())
1336 0 : return false;
1337 :
1338 : // If the mesh thinks it's prepared in some way, *re*-preparing in
1339 : // that way shouldn't change a clone of it, as long as we disallow
1340 : // repartitioning or renumbering or remote element removal.
1341 1822 : std::unique_ptr<MeshBase> mesh_clone = reprepared_mesh_clone(mesh);
1342 :
1343 : // Check whether the original and clone compare equal
1344 1772 : return (mesh == *mesh_clone);
1345 1672 : }
1346 :
1347 :
1348 :
1349 : #ifndef NDEBUG
1350 :
1351 :
1352 12368 : void libmesh_assert_valid_is_prepared (const MeshBase & mesh)
1353 : {
1354 12368 : LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools");
1355 :
1356 12368 : const MeshBase::Preparation prep = mesh.preparation();
1357 :
1358 : // If the mesh doesn't think *anything* has been prepared, we have
1359 : // nothing to check.
1360 12368 : if (prep == MeshBase::Preparation())
1361 2 : return;
1362 :
1363 : // If the mesh thinks it's partitioned, check. These are counts,
1364 : // not caches, so it's a real check.
1365 12366 : if (prep.is_partitioned)
1366 12366 : libmesh_assert_msg
1367 : (!mesh.n_unpartitioned_elem() && !mesh.n_unpartitioned_nodes(),
1368 : "Mesh preparation().is_partitioned does not reflect mesh data.");
1369 :
1370 : // If the mesh thinks it's prepared in some way, *re*-preparing in
1371 : // that way shouldn't change a clone of it, as long as we disallow
1372 : // repartitioning or renumbering or remote element removal.
1373 12366 : std::unique_ptr<MeshBase> mesh_clone = reprepared_mesh_clone(mesh);
1374 :
1375 : mesh.assert_equal_to
1376 12366 : (*mesh_clone,
1377 : "Mesh data does not match mesh preparation().\n"
1378 12366 : "Efficiently-prepared mesh != exhaustively-prepared clone.");
1379 : }
1380 :
1381 :
1382 :
1383 :
1384 :
1385 136 : void libmesh_assert_equal_n_systems (const MeshBase & mesh)
1386 : {
1387 272 : LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
1388 :
1389 136 : unsigned int n_sys = libMesh::invalid_uint;
1390 :
1391 12632 : for (const auto & elem : mesh.element_ptr_range())
1392 : {
1393 12496 : if (n_sys == libMesh::invalid_uint)
1394 136 : n_sys = elem->n_systems();
1395 : else
1396 12360 : libmesh_assert_equal_to (elem->n_systems(), n_sys);
1397 : }
1398 :
1399 31283 : for (const auto & node : mesh.node_ptr_range())
1400 : {
1401 31147 : if (n_sys == libMesh::invalid_uint)
1402 0 : n_sys = node->n_systems();
1403 : else
1404 31147 : libmesh_assert_equal_to (node->n_systems(), n_sys);
1405 : }
1406 136 : }
1407 :
1408 :
1409 :
1410 : #ifdef LIBMESH_ENABLE_AMR
1411 0 : void libmesh_assert_old_dof_objects (const MeshBase & mesh)
1412 : {
1413 0 : LOG_SCOPE("libmesh_assert_old_dof_objects()", "MeshTools");
1414 :
1415 0 : for (const auto & elem : mesh.element_ptr_range())
1416 : {
1417 0 : if (elem->refinement_flag() == Elem::JUST_REFINED ||
1418 0 : elem->refinement_flag() == Elem::INACTIVE)
1419 0 : continue;
1420 :
1421 0 : if (elem->has_dofs())
1422 0 : libmesh_assert(elem->get_old_dof_object());
1423 :
1424 0 : for (auto & node : elem->node_ref_range())
1425 0 : if (node.has_dofs())
1426 0 : libmesh_assert(node.get_old_dof_object());
1427 : }
1428 0 : }
1429 : #else
1430 : void libmesh_assert_old_dof_objects (const MeshBase &) {}
1431 : #endif // LIBMESH_ENABLE_AMR
1432 :
1433 :
1434 :
1435 0 : void libmesh_assert_valid_node_pointers(const MeshBase & mesh)
1436 : {
1437 0 : LOG_SCOPE("libmesh_assert_valid_node_pointers()", "MeshTools");
1438 :
1439 : // Here we specifically do not want "auto &" because we need to
1440 : // reseat the (temporary) pointer variable in the loop below,
1441 : // without modifying the original.
1442 0 : for (const Elem * elem : mesh.element_ptr_range())
1443 : {
1444 0 : libmesh_assert (elem);
1445 0 : while (elem)
1446 : {
1447 0 : elem->libmesh_assert_valid_node_pointers();
1448 0 : for (auto n : elem->neighbor_ptr_range())
1449 0 : if (n && n != remote_elem)
1450 0 : n->libmesh_assert_valid_node_pointers();
1451 :
1452 0 : libmesh_assert_not_equal_to (elem->parent(), remote_elem);
1453 0 : elem = elem->parent();
1454 : }
1455 : }
1456 0 : }
1457 :
1458 :
1459 :
1460 9580 : void libmesh_assert_valid_remote_elems(const MeshBase & mesh)
1461 : {
1462 19160 : LOG_SCOPE("libmesh_assert_valid_remote_elems()", "MeshTools");
1463 :
1464 667815 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1465 1335630 : mesh.local_elements_end()))
1466 : {
1467 658235 : libmesh_assert (elem);
1468 :
1469 : // We currently don't allow active_local_elements to have
1470 : // remote_elem neighbors
1471 658235 : if (elem->active())
1472 2851080 : for (auto n : elem->neighbor_ptr_range())
1473 2293813 : libmesh_assert_not_equal_to (n, remote_elem);
1474 :
1475 : #ifdef LIBMESH_ENABLE_AMR
1476 658235 : const Elem * parent = elem->parent();
1477 658235 : if (parent)
1478 377230 : libmesh_assert_not_equal_to (parent, remote_elem);
1479 :
1480 : // We can only be strict about active elements' subactive
1481 : // children
1482 658235 : if (elem->active() && elem->has_children())
1483 30156 : for (auto & child : elem->child_ref_range())
1484 24144 : libmesh_assert_not_equal_to (&child, remote_elem);
1485 : #endif
1486 : }
1487 9580 : }
1488 :
1489 :
1490 :
1491 392 : void libmesh_assert_valid_elem_ids(const MeshBase & mesh)
1492 : {
1493 784 : LOG_SCOPE("libmesh_assert_valid_elem_ids()", "MeshTools");
1494 :
1495 392 : processor_id_type lastprocid = 0;
1496 392 : dof_id_type lastelemid = 0;
1497 :
1498 13631 : for (const auto & elem : mesh.active_element_ptr_range())
1499 : {
1500 13239 : libmesh_assert (elem);
1501 13239 : processor_id_type elemprocid = elem->processor_id();
1502 13239 : dof_id_type elemid = elem->id();
1503 :
1504 13239 : libmesh_assert_greater_equal (elemid, lastelemid);
1505 13239 : libmesh_assert_greater_equal (elemprocid, lastprocid);
1506 :
1507 13239 : lastelemid = elemid;
1508 13239 : lastprocid = elemprocid;
1509 : }
1510 392 : }
1511 :
1512 :
1513 :
1514 13204 : void libmesh_assert_valid_amr_elem_ids(const MeshBase & mesh)
1515 : {
1516 26408 : LOG_SCOPE("libmesh_assert_valid_amr_elem_ids()", "MeshTools");
1517 :
1518 1676408 : for (const auto & elem : mesh.element_ptr_range())
1519 : {
1520 1663204 : libmesh_assert (elem);
1521 :
1522 1663204 : const Elem * parent = elem->parent();
1523 :
1524 1663204 : if (parent)
1525 : {
1526 838354 : libmesh_assert_greater_equal (elem->id(), parent->id());
1527 838354 : libmesh_assert_greater_equal (elem->processor_id(), parent->processor_id());
1528 : }
1529 : }
1530 13204 : }
1531 :
1532 :
1533 :
1534 25796 : void libmesh_assert_valid_amr_interior_parents(const MeshBase & mesh)
1535 : {
1536 51592 : LOG_SCOPE("libmesh_assert_valid_amr_interior_parents()", "MeshTools");
1537 :
1538 2792707 : for (const auto & elem : mesh.element_ptr_range())
1539 : {
1540 2766911 : libmesh_assert (elem);
1541 :
1542 : // We can skip to the next element if we're full-dimension
1543 : // and therefore don't have any interior parents
1544 2766911 : if (elem->dim() >= LIBMESH_DIM)
1545 1202018 : continue;
1546 :
1547 1564893 : const Elem * ip = elem->interior_parent();
1548 :
1549 1564893 : const Elem * parent = elem->parent();
1550 :
1551 1564893 : if (ip && (ip != remote_elem) && parent)
1552 : {
1553 1952 : libmesh_assert_equal_to (ip->top_parent(),
1554 : elem->top_parent()->interior_parent());
1555 :
1556 1952 : if (ip->level() == elem->level())
1557 1952 : libmesh_assert_equal_to (ip->parent(),
1558 : parent->interior_parent());
1559 : else
1560 : {
1561 0 : libmesh_assert_less (ip->level(), elem->level());
1562 0 : libmesh_assert_equal_to (ip, parent->interior_parent());
1563 : }
1564 : }
1565 : }
1566 25796 : }
1567 :
1568 :
1569 :
1570 0 : void libmesh_assert_contiguous_dof_ids(const MeshBase & mesh, unsigned int sysnum)
1571 : {
1572 0 : LOG_SCOPE("libmesh_assert_contiguous_dof_ids()", "MeshTools");
1573 :
1574 0 : if (mesh.n_processors() == 1)
1575 0 : return;
1576 :
1577 0 : libmesh_parallel_only(mesh.comm());
1578 :
1579 0 : dof_id_type min_dof_id = std::numeric_limits<dof_id_type>::max(),
1580 0 : max_dof_id = std::numeric_limits<dof_id_type>::min();
1581 :
1582 : // Figure out what our local dof id range is
1583 0 : for (const auto * node : mesh.local_node_ptr_range())
1584 : {
1585 0 : for (auto v : make_range(node->n_vars(sysnum)))
1586 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1587 : {
1588 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1589 0 : min_dof_id = std::min (min_dof_id, id);
1590 0 : max_dof_id = std::max (max_dof_id, id);
1591 : }
1592 : }
1593 :
1594 : // Make sure no other processors' ids are inside it
1595 0 : for (const auto * node : mesh.node_ptr_range())
1596 : {
1597 0 : if (node->processor_id() == mesh.processor_id())
1598 0 : continue;
1599 0 : for (auto v : make_range(node->n_vars(sysnum)))
1600 0 : for (auto c : make_range(node->n_comp(sysnum, v)))
1601 : {
1602 0 : dof_id_type id = node->dof_number(sysnum, v, c);
1603 0 : libmesh_assert (id < min_dof_id ||
1604 : id > max_dof_id);
1605 : }
1606 : }
1607 : }
1608 :
1609 :
1610 :
1611 : template <>
1612 19576 : void libmesh_assert_topology_consistent_procids<Elem>(const MeshBase & mesh)
1613 : {
1614 39152 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1615 :
1616 : // This parameter is not used when !LIBMESH_ENABLE_AMR
1617 19576 : libmesh_ignore(mesh);
1618 :
1619 : // If we're adaptively refining, check processor ids for consistency
1620 : // between parents and children.
1621 : #ifdef LIBMESH_ENABLE_AMR
1622 :
1623 : // Ancestor elements we won't worry about, but subactive and active
1624 : // elements ought to have parents with consistent processor ids
1625 2733790 : for (const auto & elem : mesh.element_ptr_range())
1626 : {
1627 2714214 : libmesh_assert(elem);
1628 :
1629 2714214 : if (!elem->active() && !elem->subactive())
1630 306752 : continue;
1631 :
1632 2407462 : const Elem * parent = elem->parent();
1633 :
1634 2407462 : if (parent)
1635 : {
1636 1225116 : libmesh_assert(parent->has_children());
1637 1225116 : processor_id_type parent_procid = parent->processor_id();
1638 1225116 : bool matching_child_id = false;
1639 : // If we've got a remote_elem then we don't know whether
1640 : // it's responsible for the parent's processor id; all
1641 : // we can do is assume it is and let its processor fail
1642 : // an assert if there's something wrong.
1643 7551620 : for (auto & child : parent->child_ref_range())
1644 12650248 : if (&child == remote_elem ||
1645 6323744 : child.processor_id() == parent_procid)
1646 6156932 : matching_child_id = true;
1647 1225116 : libmesh_assert(matching_child_id);
1648 : }
1649 : }
1650 : #endif
1651 19576 : }
1652 :
1653 :
1654 :
1655 : template <>
1656 28044 : void libmesh_assert_topology_consistent_procids<Node>(const MeshBase & mesh)
1657 : {
1658 28044 : LOG_SCOPE("libmesh_assert_topology_consistent_procids()", "MeshTools");
1659 :
1660 28044 : if (mesh.n_processors() == 1)
1661 0 : return;
1662 :
1663 28044 : libmesh_parallel_only(mesh.comm());
1664 :
1665 : // We want this test to be valid even when called after nodes have
1666 : // been added asynchronously but before they're renumbered.
1667 : //
1668 : // Plus, some code (looking at you, stitch_meshes) modifies
1669 : // DofObject ids without keeping max_elem_id()/max_node_id()
1670 : // consistent, but that's done in a safe way for performance
1671 : // reasons, so we'll play along and just figure out new max ids
1672 : // ourselves.
1673 28044 : dof_id_type parallel_max_node_id = 0;
1674 7567583 : for (const auto & node : mesh.node_ptr_range())
1675 7539539 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
1676 7539539 : node->id()+1);
1677 28044 : mesh.comm().max(parallel_max_node_id);
1678 :
1679 :
1680 56088 : std::vector<bool> node_touched_by_me(parallel_max_node_id, false);
1681 :
1682 2188471 : for (const auto & elem : as_range(mesh.local_elements_begin(),
1683 4376942 : mesh.local_elements_end()))
1684 : {
1685 2160427 : libmesh_assert (elem);
1686 :
1687 16486694 : for (auto & node : elem->node_ref_range())
1688 : {
1689 14326267 : dof_id_type nodeid = node.id();
1690 14326267 : node_touched_by_me[nodeid] = true;
1691 : }
1692 : }
1693 56088 : std::vector<bool> node_touched_by_anyone(node_touched_by_me);
1694 28044 : mesh.comm().max(node_touched_by_anyone);
1695 :
1696 3785713 : for (const auto & node : mesh.local_node_ptr_range())
1697 : {
1698 3757669 : libmesh_assert(node);
1699 3757669 : dof_id_type nodeid = node->id();
1700 3757669 : libmesh_assert(!node_touched_by_anyone[nodeid] ||
1701 : node_touched_by_me[nodeid]);
1702 : }
1703 : }
1704 :
1705 :
1706 :
1707 0 : void libmesh_assert_canonical_node_procids (const MeshBase & mesh)
1708 : {
1709 0 : for (const auto & elem : mesh.active_element_ptr_range())
1710 0 : for (auto & node : elem->node_ref_range())
1711 0 : libmesh_assert_equal_to
1712 : (node.processor_id(),
1713 : node.choose_processor_id(node.processor_id(),
1714 : elem->processor_id()));
1715 0 : }
1716 :
1717 :
1718 :
1719 : #ifdef LIBMESH_ENABLE_AMR
1720 1282 : void libmesh_assert_valid_refinement_tree(const MeshBase & mesh)
1721 : {
1722 2564 : LOG_SCOPE("libmesh_assert_valid_refinement_tree()", "MeshTools");
1723 :
1724 67616 : for (const auto & elem : mesh.element_ptr_range())
1725 : {
1726 66334 : libmesh_assert(elem);
1727 66334 : if (elem->has_children())
1728 27268 : for (auto & child : elem->child_ref_range())
1729 23216 : if (&child != remote_elem)
1730 20800 : libmesh_assert_equal_to (child.parent(), elem);
1731 66334 : if (elem->active())
1732 : {
1733 62282 : libmesh_assert(!elem->ancestor());
1734 62282 : libmesh_assert(!elem->subactive());
1735 : }
1736 4052 : else if (elem->ancestor())
1737 : {
1738 4052 : libmesh_assert(!elem->subactive());
1739 : }
1740 : else
1741 0 : libmesh_assert(elem->subactive());
1742 :
1743 66334 : if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1744 0 : libmesh_assert_greater(elem->p_level(), 0);
1745 : }
1746 1282 : }
1747 : #else
1748 : void libmesh_assert_valid_refinement_tree(const MeshBase &)
1749 : {
1750 : }
1751 : #endif // LIBMESH_ENABLE_AMR
1752 :
1753 : #endif // !NDEBUG
1754 :
1755 :
1756 :
1757 : #ifdef DEBUG
1758 :
1759 0 : void libmesh_assert_no_links_to_elem(const MeshBase & mesh,
1760 : const Elem * bad_elem)
1761 : {
1762 0 : for (const auto & elem : mesh.element_ptr_range())
1763 : {
1764 0 : libmesh_assert (elem);
1765 0 : libmesh_assert_not_equal_to (elem->parent(), bad_elem);
1766 0 : for (auto n : elem->neighbor_ptr_range())
1767 0 : libmesh_assert_not_equal_to (n, bad_elem);
1768 :
1769 : #ifdef LIBMESH_ENABLE_AMR
1770 0 : if (elem->has_children())
1771 0 : for (auto & child : elem->child_ref_range())
1772 0 : libmesh_assert_not_equal_to (&child, bad_elem);
1773 : #endif
1774 : }
1775 0 : }
1776 :
1777 :
1778 570 : void libmesh_assert_equal_points (const MeshBase & mesh)
1779 : {
1780 1140 : LOG_SCOPE("libmesh_assert_equal_points()", "MeshTools");
1781 :
1782 570 : dof_id_type pmax_node_id = mesh.max_node_id();
1783 570 : mesh.comm().max(pmax_node_id);
1784 :
1785 922414 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1786 : {
1787 921844 : const Point * p = mesh.query_node_ptr(i);
1788 :
1789 921844 : libmesh_assert(mesh.comm().semiverify(p));
1790 : }
1791 570 : }
1792 :
1793 :
1794 570 : void libmesh_assert_equal_connectivity (const MeshBase & mesh)
1795 : {
1796 1140 : LOG_SCOPE("libmesh_assert_equal_connectivity()", "MeshTools");
1797 :
1798 570 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1799 570 : mesh.comm().max(pmax_elem_id);
1800 :
1801 1144326 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1802 : {
1803 1143756 : const Elem * e = mesh.query_elem_ptr(i);
1804 :
1805 2287512 : std::vector<dof_id_type> nodes;
1806 1143756 : if (e)
1807 5792488 : for (auto n : e->node_index_range())
1808 4648848 : nodes.push_back(e->node_id(n));
1809 :
1810 1143756 : libmesh_assert(mesh.comm().semiverify(e ? &nodes : nullptr));
1811 : }
1812 570 : }
1813 :
1814 :
1815 0 : void libmesh_assert_connected_nodes (const MeshBase & mesh)
1816 : {
1817 0 : LOG_SCOPE("libmesh_assert_connected_nodes()", "MeshTools");
1818 :
1819 0 : std::set<const Node *> used_nodes;
1820 :
1821 0 : for (const auto & elem : mesh.element_ptr_range())
1822 : {
1823 0 : libmesh_assert (elem);
1824 :
1825 0 : for (auto & n : elem->node_ref_range())
1826 0 : used_nodes.insert(&n);
1827 : }
1828 :
1829 0 : for (const auto & node : mesh.node_ptr_range())
1830 : {
1831 0 : libmesh_assert(node);
1832 0 : libmesh_assert(used_nodes.count(node));
1833 : }
1834 0 : }
1835 :
1836 :
1837 :
1838 31618 : void libmesh_assert_valid_constraint_rows (const MeshBase & mesh)
1839 : {
1840 31618 : libmesh_parallel_only(mesh.comm());
1841 :
1842 31618 : const auto & constraint_rows = mesh.get_constraint_rows();
1843 :
1844 31618 : bool have_constraint_rows = !constraint_rows.empty();
1845 31618 : mesh.comm().max(have_constraint_rows);
1846 31618 : if (!have_constraint_rows)
1847 31548 : return;
1848 :
1849 24910 : for (auto & row : constraint_rows)
1850 : {
1851 24840 : const Node * node = row.first;
1852 24840 : libmesh_assert(node == mesh.node_ptr(node->id()));
1853 :
1854 84254 : for (auto & pr : row.second)
1855 : {
1856 59414 : const Elem * spline_elem = pr.first.first;
1857 59414 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1858 : }
1859 : }
1860 :
1861 70 : dof_id_type pmax_node_id = mesh.max_node_id();
1862 70 : mesh.comm().max(pmax_node_id);
1863 :
1864 33172 : for (dof_id_type i=0; i != pmax_node_id; ++i)
1865 : {
1866 33102 : const Node * node = mesh.query_node_ptr(i);
1867 :
1868 33102 : bool have_constraint = constraint_rows.count(node);
1869 :
1870 33102 : const std::size_t my_n_constraints = have_constraint ?
1871 24840 : libmesh_map_find(constraint_rows, node).size() : std::size_t(-1);
1872 33102 : const std::size_t * n_constraints = node ?
1873 33102 : &my_n_constraints : nullptr;
1874 :
1875 33102 : libmesh_assert(mesh.comm().semiverify(n_constraints));
1876 : }
1877 : }
1878 :
1879 :
1880 :
1881 70164 : void libmesh_assert_valid_boundary_ids(const MeshBase & mesh)
1882 : {
1883 70164 : LOG_SCOPE("libmesh_assert_valid_boundary_ids()", "MeshTools");
1884 :
1885 70164 : if (mesh.n_processors() == 1)
1886 1828 : return;
1887 :
1888 68336 : libmesh_parallel_only(mesh.comm());
1889 :
1890 68336 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1891 :
1892 68336 : dof_id_type pmax_elem_id = mesh.max_elem_id();
1893 68336 : mesh.comm().max(pmax_elem_id);
1894 :
1895 6820136 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
1896 : {
1897 6751800 : const Elem * elem = mesh.query_elem_ptr(i);
1898 6751800 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
1899 6751800 : const unsigned int my_n_edges = elem ? elem->n_edges() : 0;
1900 6751800 : const unsigned int my_n_sides = elem ? elem->n_sides() : 0;
1901 : unsigned int
1902 6751800 : n_nodes = my_n_nodes,
1903 6751800 : n_edges = my_n_edges,
1904 6751800 : n_sides = my_n_sides;
1905 :
1906 6751800 : mesh.comm().max(n_nodes);
1907 6751800 : mesh.comm().max(n_edges);
1908 6751800 : mesh.comm().max(n_sides);
1909 :
1910 6751800 : if (elem)
1911 : {
1912 6603715 : libmesh_assert_equal_to(my_n_nodes, n_nodes);
1913 6603715 : libmesh_assert_equal_to(my_n_edges, n_edges);
1914 6603715 : libmesh_assert_equal_to(my_n_sides, n_sides);
1915 : }
1916 :
1917 : // Let's test all IDs on the element with one communication
1918 : // rather than n_nodes + n_edges + n_sides communications, to
1919 : // cut down on latency in dbg modes.
1920 13503600 : std::vector<boundary_id_type> all_bcids;
1921 :
1922 54321782 : for (unsigned int n=0; n != n_nodes; ++n)
1923 : {
1924 95139964 : std::vector<boundary_id_type> bcids;
1925 47569982 : if (elem)
1926 : {
1927 47299718 : boundary_info.boundary_ids(elem->node_ptr(n), bcids);
1928 :
1929 : // Ordering of boundary ids shouldn't matter
1930 47299718 : std::sort(bcids.begin(), bcids.end());
1931 : }
1932 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1933 :
1934 47569982 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1935 95139964 : bcids.end());
1936 : // Separator
1937 47569982 : all_bcids.push_back(BoundaryInfo::invalid_id);
1938 : }
1939 :
1940 42034188 : for (unsigned short e=0; e != n_edges; ++e)
1941 : {
1942 70564776 : std::vector<boundary_id_type> bcids;
1943 :
1944 35282388 : if (elem)
1945 : {
1946 35079474 : boundary_info.edge_boundary_ids(elem, e, bcids);
1947 :
1948 : // Ordering of boundary ids shouldn't matter
1949 35079474 : std::sort(bcids.begin(), bcids.end());
1950 : }
1951 :
1952 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1953 :
1954 35282388 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1955 70564776 : bcids.end());
1956 : // Separator
1957 35282388 : all_bcids.push_back(BoundaryInfo::invalid_id);
1958 :
1959 35282388 : if (elem)
1960 : {
1961 35079474 : boundary_info.raw_edge_boundary_ids(elem, e, bcids);
1962 :
1963 : // Ordering of boundary ids shouldn't matter
1964 35079474 : std::sort(bcids.begin(), bcids.end());
1965 :
1966 35079474 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1967 70158948 : bcids.end());
1968 : // Separator
1969 35079474 : all_bcids.push_back(BoundaryInfo::invalid_id);
1970 : }
1971 :
1972 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1973 : }
1974 :
1975 33432560 : for (unsigned short s=0; s != n_sides; ++s)
1976 : {
1977 53361520 : std::vector<boundary_id_type> bcids;
1978 :
1979 26680760 : if (elem)
1980 : {
1981 26545434 : boundary_info.boundary_ids(elem, s, bcids);
1982 :
1983 : // Ordering of boundary ids shouldn't matter
1984 26545434 : std::sort(bcids.begin(), bcids.end());
1985 :
1986 26545434 : all_bcids.insert(all_bcids.end(), bcids.begin(),
1987 53090868 : bcids.end());
1988 : // Separator
1989 26545434 : all_bcids.push_back(BoundaryInfo::invalid_id);
1990 : }
1991 :
1992 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
1993 :
1994 26680760 : if (elem)
1995 : {
1996 26545434 : boundary_info.raw_boundary_ids(elem, s, bcids);
1997 :
1998 : // Ordering of boundary ids shouldn't matter
1999 26545434 : std::sort(bcids.begin(), bcids.end());
2000 :
2001 26545434 : all_bcids.insert(all_bcids.end(), bcids.begin(),
2002 53090868 : bcids.end());
2003 : // Separator
2004 26545434 : all_bcids.push_back(BoundaryInfo::invalid_id);
2005 : }
2006 :
2007 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
2008 : }
2009 :
2010 20255400 : for (unsigned short sf=0; sf != 2; ++sf)
2011 : {
2012 27007200 : std::vector<boundary_id_type> bcids;
2013 :
2014 13503600 : if (elem)
2015 : {
2016 13207430 : boundary_info.shellface_boundary_ids(elem, sf, bcids);
2017 :
2018 : // Ordering of boundary ids shouldn't matter
2019 13207430 : std::sort(bcids.begin(), bcids.end());
2020 :
2021 13207430 : all_bcids.insert(all_bcids.end(), bcids.begin(),
2022 26414860 : bcids.end());
2023 : // Separator
2024 13207430 : all_bcids.push_back(BoundaryInfo::invalid_id);
2025 : }
2026 :
2027 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
2028 :
2029 13503600 : if (elem)
2030 : {
2031 13207430 : boundary_info.raw_shellface_boundary_ids(elem, sf, bcids);
2032 :
2033 : // Ordering of boundary ids shouldn't matter
2034 13207430 : std::sort(bcids.begin(), bcids.end());
2035 :
2036 13207430 : all_bcids.insert(all_bcids.end(), bcids.begin(),
2037 26414860 : bcids.end());
2038 : // Separator
2039 13207430 : all_bcids.push_back(BoundaryInfo::invalid_id);
2040 : }
2041 :
2042 : // libmesh_assert(mesh.comm().semiverify (elem ? &bcids : nullptr));
2043 : }
2044 :
2045 6751800 : libmesh_assert(mesh.comm().semiverify
2046 : (elem ? &all_bcids : nullptr));
2047 : }
2048 : }
2049 :
2050 :
2051 7910 : void libmesh_assert_valid_dof_ids(const MeshBase & mesh, unsigned int sysnum)
2052 : {
2053 7910 : LOG_SCOPE("libmesh_assert_valid_dof_ids()", "MeshTools");
2054 :
2055 7910 : if (mesh.n_processors() == 1)
2056 0 : return;
2057 :
2058 7910 : libmesh_parallel_only(mesh.comm());
2059 :
2060 7910 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2061 7910 : mesh.comm().max(pmax_elem_id);
2062 :
2063 1045320 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
2064 1037410 : assert_semiverify_dofobj(mesh.comm(),
2065 1037410 : mesh.query_elem_ptr(i),
2066 : sysnum);
2067 :
2068 7910 : dof_id_type pmax_node_id = mesh.max_node_id();
2069 7910 : mesh.comm().max(pmax_node_id);
2070 :
2071 1951402 : for (dof_id_type i=0; i != pmax_node_id; ++i)
2072 1943492 : assert_semiverify_dofobj(mesh.comm(),
2073 1943492 : mesh.query_node_ptr(i),
2074 : sysnum);
2075 : }
2076 :
2077 :
2078 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2079 62204 : void libmesh_assert_valid_unique_ids(const MeshBase & mesh)
2080 : {
2081 124408 : LOG_SCOPE("libmesh_assert_valid_unique_ids()", "MeshTools");
2082 :
2083 62204 : libmesh_parallel_only(mesh.comm());
2084 :
2085 : // Storage for semi-local DofObject ids.
2086 124408 : std::unordered_set<unique_id_type> semilocal_unique_ids;
2087 :
2088 62204 : auto gather_elem_ids = [&]()
2089 : {
2090 4699191 : for (auto const & elem : mesh.active_element_ptr_range())
2091 : {
2092 4636987 : auto [it, inserted] = semilocal_unique_ids.insert(elem->unique_id());
2093 4636987 : libmesh_assert(inserted);
2094 4636987 : libmesh_ignore(it);
2095 : }
2096 62204 : };
2097 :
2098 62204 : auto gather_node_ids = [&]()
2099 : {
2100 10694297 : for (auto const & node : mesh.node_ptr_range())
2101 : {
2102 10632093 : auto [it, inserted] = semilocal_unique_ids.insert(node->unique_id());
2103 10632093 : libmesh_assert(inserted);
2104 10632093 : libmesh_ignore(it);
2105 : }
2106 62204 : };
2107 :
2108 62204 : auto verify_elems = [&]()
2109 : {
2110 62204 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2111 62204 : mesh.comm().max(pmax_elem_id);
2112 :
2113 5798880 : for (auto i : make_range(pmax_elem_id))
2114 : {
2115 5736676 : const Elem * elem = mesh.query_elem_ptr(i);
2116 5736676 : assert_dofobj_unique_id(mesh.comm(), elem, semilocal_unique_ids);
2117 : }
2118 62204 : };
2119 :
2120 62204 : auto verify_nodes = [&]()
2121 : {
2122 62204 : dof_id_type pmax_node_id = mesh.max_node_id();
2123 62204 : mesh.comm().max(pmax_node_id);
2124 :
2125 11216050 : for (auto i : make_range(pmax_node_id))
2126 : {
2127 11153846 : const Node * node = mesh.query_node_ptr(i);
2128 11153846 : assert_dofobj_unique_id(mesh.comm(), node, semilocal_unique_ids);
2129 : }
2130 62204 : };
2131 :
2132 62204 : if (!mesh.allow_node_and_elem_unique_id_overlap())
2133 : {
2134 : // First collect all the unique_ids we can see and make sure there's
2135 : // no duplicates
2136 62204 : gather_elem_ids();
2137 62204 : gather_node_ids();
2138 :
2139 : // Then make sure elements/nodes are all in sync and remote
2140 : // elements don't duplicate semilocal
2141 62204 : verify_elems();
2142 62204 : verify_nodes();
2143 : }
2144 : else
2145 : {
2146 : // If the mesh allows Node and Elem unique_ids to overlap, then we only
2147 : // check for validity and uniqueness of an Elem (resp. Node) unique id
2148 : // within the set of Elem (resp. Node) unique_ids.
2149 0 : gather_elem_ids();
2150 0 : verify_elems();
2151 :
2152 : // Clear id list before checking Nodes
2153 0 : semilocal_unique_ids.clear();
2154 :
2155 : // Finally, check Nodes
2156 0 : gather_node_ids();
2157 0 : verify_nodes();
2158 : }
2159 62204 : }
2160 : #endif
2161 :
2162 0 : void libmesh_assert_consistent_distributed(const MeshBase & mesh)
2163 : {
2164 0 : libmesh_parallel_only(mesh.comm());
2165 :
2166 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2167 0 : mesh.comm().max(parallel_max_elem_id);
2168 :
2169 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2170 : {
2171 0 : const Elem * elem = mesh.query_elem_ptr(i);
2172 : processor_id_type pid =
2173 0 : elem ? elem->processor_id() : DofObject::invalid_processor_id;
2174 0 : mesh.comm().min(pid);
2175 0 : libmesh_assert(elem || pid != mesh.processor_id());
2176 : }
2177 :
2178 0 : dof_id_type parallel_max_node_id = mesh.max_node_id();
2179 0 : mesh.comm().max(parallel_max_node_id);
2180 :
2181 0 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2182 : {
2183 0 : const Node * node = mesh.query_node_ptr(i);
2184 : processor_id_type pid =
2185 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2186 0 : mesh.comm().min(pid);
2187 0 : libmesh_assert(node || pid != mesh.processor_id());
2188 : }
2189 0 : }
2190 :
2191 :
2192 0 : void libmesh_assert_consistent_distributed_nodes(const MeshBase & mesh)
2193 : {
2194 0 : libmesh_parallel_only(mesh.comm());
2195 0 : auto locator = mesh.sub_point_locator();
2196 :
2197 0 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2198 0 : mesh.comm().max(parallel_max_elem_id);
2199 :
2200 0 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2201 : {
2202 0 : const Elem * elem = mesh.query_elem_ptr(i);
2203 :
2204 0 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2205 0 : unsigned int n_nodes = my_n_nodes;
2206 0 : mesh.comm().max(n_nodes);
2207 :
2208 0 : if (n_nodes)
2209 0 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2210 :
2211 0 : for (unsigned int n=0; n != n_nodes; ++n)
2212 : {
2213 0 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2214 : processor_id_type pid =
2215 0 : node ? node->processor_id() : DofObject::invalid_processor_id;
2216 0 : mesh.comm().min(pid);
2217 0 : libmesh_assert(node || pid != mesh.processor_id());
2218 : }
2219 : }
2220 0 : }
2221 :
2222 :
2223 :
2224 : template <>
2225 19576 : void libmesh_assert_parallel_consistent_procids<Elem>(const MeshBase & mesh)
2226 : {
2227 19576 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2228 :
2229 19576 : if (mesh.n_processors() == 1)
2230 0 : return;
2231 :
2232 19576 : libmesh_parallel_only(mesh.comm());
2233 :
2234 : // Some code (looking at you, stitch_meshes) modifies DofObject ids
2235 : // without keeping max_elem_id()/max_node_id() consistent, but
2236 : // that's done in a safe way for performance reasons, so we'll play
2237 : // along and just figure out new max ids ourselves.
2238 19576 : dof_id_type parallel_max_elem_id = 0;
2239 2733790 : for (const auto & elem : mesh.element_ptr_range())
2240 2714214 : parallel_max_elem_id = std::max<dof_id_type>(parallel_max_elem_id,
2241 2714214 : elem->id()+1);
2242 19576 : mesh.comm().max(parallel_max_elem_id);
2243 :
2244 : // Check processor ids for consistency between processors
2245 :
2246 2804840 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2247 : {
2248 2785264 : const Elem * elem = mesh.query_elem_ptr(i);
2249 :
2250 : processor_id_type min_id =
2251 2714214 : elem ? elem->processor_id() :
2252 2785264 : std::numeric_limits<processor_id_type>::max();
2253 2785264 : mesh.comm().min(min_id);
2254 :
2255 : processor_id_type max_id =
2256 2714214 : elem ? elem->processor_id() :
2257 2785264 : std::numeric_limits<processor_id_type>::min();
2258 2785264 : mesh.comm().max(max_id);
2259 :
2260 2785264 : if (elem)
2261 : {
2262 2714214 : libmesh_assert_equal_to (min_id, elem->processor_id());
2263 2714214 : libmesh_assert_equal_to (max_id, elem->processor_id());
2264 : }
2265 :
2266 2785264 : if (min_id == mesh.processor_id())
2267 1318415 : libmesh_assert(elem);
2268 : }
2269 : }
2270 :
2271 :
2272 :
2273 76 : void libmesh_assert_parallel_consistent_new_node_procids(const MeshBase & mesh)
2274 : {
2275 76 : LOG_SCOPE("libmesh_assert_parallel_consistent_new_node_procids()", "MeshTools");
2276 :
2277 76 : if (mesh.n_processors() == 1)
2278 0 : return;
2279 :
2280 76 : libmesh_parallel_only(mesh.comm());
2281 :
2282 : // We want this test to hit every node when called even after nodes
2283 : // have been added asynchronously but before everything has been
2284 : // renumbered.
2285 76 : dof_id_type parallel_max_elem_id = mesh.max_elem_id();
2286 76 : mesh.comm().max(parallel_max_elem_id);
2287 :
2288 152 : std::vector<bool> elem_touched_by_anyone(parallel_max_elem_id, false);
2289 :
2290 43200 : for (dof_id_type i=0; i != parallel_max_elem_id; ++i)
2291 : {
2292 43124 : const Elem * elem = mesh.query_elem_ptr(i);
2293 :
2294 43124 : const unsigned int my_n_nodes = elem ? elem->n_nodes() : 0;
2295 43124 : unsigned int n_nodes = my_n_nodes;
2296 43124 : mesh.comm().max(n_nodes);
2297 :
2298 43124 : if (n_nodes)
2299 15840 : libmesh_assert(mesh.comm().semiverify(elem ? &my_n_nodes : nullptr));
2300 :
2301 173356 : for (unsigned int n=0; n != n_nodes; ++n)
2302 : {
2303 130232 : const Node * node = elem ? elem->node_ptr(n) : nullptr;
2304 130232 : const processor_id_type pid = node ? node->processor_id() : 0;
2305 130232 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2306 : }
2307 : }
2308 : }
2309 :
2310 : template <>
2311 41494 : void libmesh_assert_parallel_consistent_procids<Node>(const MeshBase & mesh)
2312 : {
2313 41494 : LOG_SCOPE("libmesh_assert_parallel_consistent_procids()", "MeshTools");
2314 :
2315 41494 : if (mesh.n_processors() == 1)
2316 342 : return;
2317 :
2318 41152 : libmesh_parallel_only(mesh.comm());
2319 :
2320 : // We want this test to be valid even when called even after nodes
2321 : // have been added asynchronously but before they're renumbered
2322 : //
2323 : // Plus, some code (looking at you, stitch_meshes) modifies
2324 : // DofObject ids without keeping max_elem_id()/max_node_id()
2325 : // consistent, but that's done in a safe way for performance
2326 : // reasons, so we'll play along and just figure out new max ids
2327 : // ourselves.
2328 41152 : dof_id_type parallel_max_node_id = 0;
2329 10505165 : for (const auto & node : mesh.node_ptr_range())
2330 10464013 : parallel_max_node_id = std::max<dof_id_type>(parallel_max_node_id,
2331 10464013 : node->id()+1);
2332 41152 : mesh.comm().max(parallel_max_node_id);
2333 :
2334 82304 : std::vector<bool> node_touched_by_anyone(parallel_max_node_id, false);
2335 :
2336 3080421 : for (const auto & elem : as_range(mesh.local_elements_begin(),
2337 6160842 : mesh.local_elements_end()))
2338 : {
2339 3039269 : libmesh_assert (elem);
2340 :
2341 23359846 : for (auto & node : elem->node_ref_range())
2342 : {
2343 20320577 : dof_id_type nodeid = node.id();
2344 20320577 : node_touched_by_anyone[nodeid] = true;
2345 : }
2346 : }
2347 41152 : mesh.comm().max(node_touched_by_anyone);
2348 :
2349 : // Check processor ids for consistency between processors
2350 : // on any node an element touches
2351 11030202 : for (dof_id_type i=0; i != parallel_max_node_id; ++i)
2352 : {
2353 10989050 : if (!node_touched_by_anyone[i])
2354 518400 : continue;
2355 :
2356 10470650 : const Node * node = mesh.query_node_ptr(i);
2357 10470650 : const processor_id_type pid = node ? node->processor_id() : 0;
2358 :
2359 10470650 : libmesh_assert(mesh.comm().semiverify (node ? &pid : nullptr));
2360 : }
2361 : }
2362 :
2363 :
2364 :
2365 : #ifdef LIBMESH_ENABLE_AMR
2366 0 : void libmesh_assert_valid_refinement_flags(const MeshBase & mesh)
2367 : {
2368 0 : LOG_SCOPE("libmesh_assert_valid_refinement_flags()", "MeshTools");
2369 :
2370 0 : libmesh_parallel_only(mesh.comm());
2371 0 : if (mesh.n_processors() == 1)
2372 0 : return;
2373 :
2374 0 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2375 0 : mesh.comm().max(pmax_elem_id);
2376 :
2377 0 : std::vector<unsigned char> my_elem_h_state(pmax_elem_id, 255);
2378 0 : std::vector<unsigned char> my_elem_p_state(pmax_elem_id, 255);
2379 :
2380 0 : for (const auto & elem : mesh.element_ptr_range())
2381 : {
2382 0 : libmesh_assert (elem);
2383 0 : dof_id_type elemid = elem->id();
2384 :
2385 0 : my_elem_h_state[elemid] =
2386 0 : static_cast<unsigned char>(elem->refinement_flag());
2387 :
2388 0 : my_elem_p_state[elemid] =
2389 0 : static_cast<unsigned char>(elem->p_refinement_flag());
2390 : }
2391 0 : std::vector<unsigned char> min_elem_h_state(my_elem_h_state);
2392 0 : mesh.comm().min(min_elem_h_state);
2393 :
2394 0 : std::vector<unsigned char> min_elem_p_state(my_elem_p_state);
2395 0 : mesh.comm().min(min_elem_p_state);
2396 :
2397 0 : for (dof_id_type i=0; i!= pmax_elem_id; ++i)
2398 : {
2399 0 : libmesh_assert(my_elem_h_state[i] == 255 ||
2400 : my_elem_h_state[i] == min_elem_h_state[i]);
2401 0 : libmesh_assert(my_elem_p_state[i] == 255 ||
2402 : my_elem_p_state[i] == min_elem_p_state[i]);
2403 : }
2404 : }
2405 : #else
2406 : void libmesh_assert_valid_refinement_flags(const MeshBase &)
2407 : {
2408 : }
2409 : #endif // LIBMESH_ENABLE_AMR
2410 :
2411 :
2412 :
2413 29044 : void libmesh_assert_valid_neighbors(const MeshBase & mesh,
2414 : bool assert_valid_remote_elems)
2415 : {
2416 29044 : LOG_SCOPE("libmesh_assert_valid_neighbors()", "MeshTools");
2417 :
2418 4200168 : for (const auto & elem : mesh.element_ptr_range())
2419 : {
2420 4171124 : libmesh_assert (elem);
2421 4171124 : elem->libmesh_assert_valid_neighbors();
2422 : }
2423 :
2424 29044 : if (mesh.n_processors() == 1)
2425 1144 : return;
2426 :
2427 27900 : libmesh_parallel_only(mesh.comm());
2428 :
2429 27900 : dof_id_type pmax_elem_id = mesh.max_elem_id();
2430 27900 : mesh.comm().max(pmax_elem_id);
2431 :
2432 4331582 : for (dof_id_type i=0; i != pmax_elem_id; ++i)
2433 : {
2434 4303682 : const Elem * elem = mesh.query_elem_ptr(i);
2435 :
2436 4303682 : const unsigned int my_n_neigh = elem ? elem->n_neighbors() : 0;
2437 4303682 : unsigned int n_neigh = my_n_neigh;
2438 4303682 : mesh.comm().max(n_neigh);
2439 4303682 : if (elem)
2440 4162008 : libmesh_assert_equal_to (my_n_neigh, n_neigh);
2441 :
2442 20899710 : for (unsigned int n = 0; n != n_neigh; ++n)
2443 : {
2444 16596028 : dof_id_type my_neighbor = DofObject::invalid_id;
2445 16596028 : dof_id_type * p_my_neighbor = nullptr;
2446 :
2447 : // If we have a non-remote_elem neighbor link, then we can
2448 : // verify it.
2449 16596028 : if (elem && elem->neighbor_ptr(n) != remote_elem)
2450 : {
2451 16473255 : p_my_neighbor = &my_neighbor;
2452 16473255 : if (elem->neighbor_ptr(n))
2453 15507102 : my_neighbor = elem->neighbor_ptr(n)->id();
2454 :
2455 : // But wait - if we haven't set remote_elem links yet then
2456 : // some nullptr links on ghost elements might be
2457 : // future-remote_elem links, so we can't verify those.
2458 32953430 : if (!assert_valid_remote_elems &&
2459 16474881 : !elem->neighbor_ptr(n) &&
2460 1626 : elem->processor_id() != mesh.processor_id())
2461 759 : p_my_neighbor = nullptr;
2462 : }
2463 16596028 : libmesh_assert(mesh.comm().semiverify(p_my_neighbor));
2464 : }
2465 : }
2466 : }
2467 : #endif // DEBUG
2468 :
2469 :
2470 :
2471 : // Functors for correct_node_proc_ids
2472 : namespace {
2473 :
2474 : typedef std::unordered_map<dof_id_type, processor_id_type> proc_id_map_type;
2475 :
2476 : struct SyncNodeSet
2477 : {
2478 : typedef unsigned char datum; // bool but without bit twiddling issues
2479 :
2480 182 : SyncNodeSet(std::unordered_set<const Node *> & _set,
2481 9199 : MeshBase & _mesh) :
2482 9199 : node_set(_set), mesh(_mesh) {}
2483 :
2484 : std::unordered_set<const Node *> & node_set;
2485 :
2486 : MeshBase & mesh;
2487 :
2488 : // ------------------------------------------------------------
2489 41122 : void gather_data (const std::vector<dof_id_type> & ids,
2490 : std::vector<datum> & data)
2491 : {
2492 : // Find whether each requested node belongs in the set
2493 41122 : data.resize(ids.size());
2494 :
2495 1867830 : for (auto i : index_range(ids))
2496 : {
2497 1826708 : const dof_id_type id = ids[i];
2498 :
2499 : // We'd better have every node we're asked for
2500 1826708 : Node * node = mesh.node_ptr(id);
2501 :
2502 : // Return if the node is in the set.
2503 2954958 : data[i] = node_set.count(node);
2504 : }
2505 41122 : }
2506 :
2507 : // ------------------------------------------------------------
2508 41122 : bool act_on_data (const std::vector<dof_id_type> & ids,
2509 : const std::vector<datum> in_set)
2510 : {
2511 167 : bool data_changed = false;
2512 :
2513 : // Add nodes we've been informed of to our own set
2514 1867830 : for (auto i : index_range(ids))
2515 : {
2516 1826708 : if (in_set[i])
2517 : {
2518 1196849 : Node * node = mesh.node_ptr(ids[i]);
2519 1196849 : if (!node_set.count(node))
2520 : {
2521 256 : node_set.insert(node);
2522 256 : data_changed = true;
2523 : }
2524 : }
2525 : }
2526 :
2527 41122 : return data_changed;
2528 : }
2529 : };
2530 :
2531 :
2532 8835 : struct NodesNotInSet
2533 : {
2534 182 : NodesNotInSet(const std::unordered_set<const Node *> _set)
2535 9017 : : node_set(_set) {}
2536 :
2537 479564 : bool operator() (const Node * node) const
2538 : {
2539 959128 : if (node_set.count(node))
2540 287060 : return false;
2541 192504 : return true;
2542 : }
2543 :
2544 : const std::unordered_set<const Node *> node_set;
2545 : };
2546 :
2547 :
2548 : struct SyncProcIdsFromMap
2549 : {
2550 : typedef processor_id_type datum;
2551 :
2552 218 : SyncProcIdsFromMap(const proc_id_map_type & _map,
2553 37675 : MeshBase & _mesh) :
2554 37675 : new_proc_ids(_map), mesh(_mesh) {}
2555 :
2556 : const proc_id_map_type & new_proc_ids;
2557 :
2558 : MeshBase & mesh;
2559 :
2560 : // ------------------------------------------------------------
2561 131006 : void gather_data (const std::vector<dof_id_type> & ids,
2562 : std::vector<datum> & data)
2563 : {
2564 : // Find the new processor id of each requested node
2565 131006 : data.resize(ids.size());
2566 :
2567 7260567 : for (auto i : index_range(ids))
2568 : {
2569 7129561 : const dof_id_type id = ids[i];
2570 :
2571 : // Return the node's new processor id if it has one, or its
2572 : // old processor id if not.
2573 7129561 : if (const auto it = new_proc_ids.find(id);
2574 47012 : it != new_proc_ids.end())
2575 6498481 : data[i] = it->second;
2576 : else
2577 : {
2578 : // We'd better find every node we're asked for
2579 631080 : const Node & node = mesh.node_ref(id);
2580 631080 : data[i] = node.processor_id();
2581 : }
2582 : }
2583 131006 : }
2584 :
2585 : // ------------------------------------------------------------
2586 131006 : void act_on_data (const std::vector<dof_id_type> & ids,
2587 : const std::vector<datum> proc_ids)
2588 : {
2589 : // Set the node processor ids we've now been informed of
2590 7260567 : for (auto i : index_range(ids))
2591 : {
2592 7129561 : Node & node = mesh.node_ref(ids[i]);
2593 7129561 : node.processor_id() = proc_ids[i];
2594 : }
2595 131006 : }
2596 : };
2597 : }
2598 :
2599 :
2600 :
2601 37675 : void correct_node_proc_ids (MeshBase & mesh)
2602 : {
2603 436 : LOG_SCOPE("correct_node_proc_ids()","MeshTools");
2604 :
2605 : // This function must be run on all processors at once
2606 218 : libmesh_parallel_only(mesh.comm());
2607 :
2608 : // We require all processors to agree on nodal processor ids before
2609 : // going into this algorithm.
2610 : #ifdef DEBUG
2611 218 : libmesh_assert_parallel_consistent_procids<Node>(mesh);
2612 : #endif
2613 :
2614 : // If we have any unpartitioned elements at this
2615 : // stage there is a problem
2616 218 : libmesh_assert (n_elem(mesh.unpartitioned_elements_begin(),
2617 : mesh.unpartitioned_elements_end()) == 0);
2618 :
2619 : // Fix nodes' processor ids. Coarsening may have left us with nodes
2620 : // which are no longer touched by any elements of the same processor
2621 : // id, and for DofMap to work we need to fix that.
2622 :
2623 : // This is harder now that libMesh no longer requires a distributed
2624 : // mesh to ghost all nodal neighbors: it is possible for two active
2625 : // elements on two different processors to share the same node in
2626 : // such a way that neither processor knows the others' element
2627 : // exists!
2628 :
2629 : // While we're at it, if this mesh is configured to allow
2630 : // repartitioning, we'll repartition *all* the nodes' processor ids
2631 : // using the canonical Node heuristic, to try and improve DoF load
2632 : // balancing. But if the mesh is disallowing repartitioning, we
2633 : // won't touch processor_id on any node where it's valid, regardless
2634 : // of whether or not it's canonical.
2635 218 : bool repartition_all_nodes = !mesh.skip_noncritical_partitioning();
2636 436 : std::unordered_set<const Node *> valid_nodes;
2637 :
2638 : // If we aren't allowed to repartition, then we're going to leave
2639 : // every node we can at its current processor_id, and *only*
2640 : // repartition the nodes whose current processor id is incompatible
2641 : // with DoFMap (because it doesn't touch an active element, e.g. due
2642 : // to coarsening)
2643 37675 : if (!repartition_all_nodes)
2644 : {
2645 1158938 : for (const auto & elem : mesh.active_element_ptr_range())
2646 5946053 : for (const auto & node : elem->node_ref_range())
2647 5305120 : if (elem->processor_id() == node.processor_id())
2648 4909802 : valid_nodes.insert(&node);
2649 :
2650 182 : SyncNodeSet syncv(valid_nodes, mesh);
2651 :
2652 : Parallel::sync_dofobject_data_by_id
2653 18216 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), syncv);
2654 : }
2655 :
2656 : // We build up a set of compatible processor ids for each node
2657 436 : proc_id_map_type new_proc_ids;
2658 :
2659 8996070 : for (auto & elem : mesh.active_element_ptr_range())
2660 : {
2661 4499519 : processor_id_type pid = elem->processor_id();
2662 :
2663 44096702 : for (auto & node : elem->node_ref_range())
2664 : {
2665 39597183 : const dof_id_type id = node.id();
2666 39597183 : if (auto it = new_proc_ids.find(id);
2667 350660 : it == new_proc_ids.end())
2668 155764 : new_proc_ids.emplace(id, pid);
2669 : else
2670 28313309 : it->second = node.choose_processor_id(it->second, pid);
2671 : }
2672 37239 : }
2673 :
2674 : // Sort the new pids to push to each processor
2675 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, processor_id_type>>>
2676 436 : ids_to_push;
2677 :
2678 25740532 : for (const auto & node : mesh.node_ptr_range())
2679 13084952 : if (const auto it = std::as_const(new_proc_ids).find(node->id());
2680 13084952 : it != new_proc_ids.end() && node->processor_id() != DofObject::invalid_processor_id)
2681 11321113 : ids_to_push[node->processor_id()].emplace_back(node->id(), /*pid=*/it->second);
2682 :
2683 : auto action_functor =
2684 171277 : [& mesh, & new_proc_ids]
2685 : (processor_id_type,
2686 11751576 : const std::vector<std::pair<dof_id_type, processor_id_type>> & data)
2687 : {
2688 11455971 : for (const auto & [id, pid] : data)
2689 : {
2690 11283874 : if (const auto it = new_proc_ids.find(id);
2691 155764 : it == new_proc_ids.end())
2692 0 : new_proc_ids.emplace(id, pid);
2693 : else
2694 : {
2695 11283874 : const Node & node = mesh.node_ref(id);
2696 11283874 : it->second = node.choose_processor_id(it->second, pid);
2697 : }
2698 : }
2699 38277 : };
2700 :
2701 : Parallel::push_parallel_vector_data
2702 37675 : (mesh.comm(), ids_to_push, action_functor);
2703 :
2704 : // Now new_proc_ids is correct for every node we used to own. Let's
2705 : // ask every other processor about the nodes they used to own. But
2706 : // first we'll need to keep track of which nodes we used to own,
2707 : // lest we get them confused with nodes we newly own.
2708 436 : std::unordered_set<Node *> ex_local_nodes;
2709 9087712 : for (auto & node : mesh.local_node_ptr_range())
2710 4642931 : if (const auto it = new_proc_ids.find(node->id());
2711 4642931 : it != new_proc_ids.end() && it->second != mesh.processor_id())
2712 37288 : ex_local_nodes.insert(node);
2713 :
2714 218 : SyncProcIdsFromMap sync(new_proc_ids, mesh);
2715 37675 : if (repartition_all_nodes)
2716 : Parallel::sync_dofobject_data_by_id
2717 56916 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync);
2718 : else
2719 : {
2720 182 : NodesNotInSet nnis(valid_nodes);
2721 :
2722 : Parallel::sync_dofobject_data_by_id
2723 18216 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), nnis, sync);
2724 : }
2725 :
2726 : // And finally let's update the nodes we used to own.
2727 42262 : for (const auto & node : ex_local_nodes)
2728 : {
2729 2103 : if (valid_nodes.count(node))
2730 2083 : continue;
2731 :
2732 2504 : const dof_id_type id = node->id();
2733 10 : const proc_id_map_type::iterator it = new_proc_ids.find(id);
2734 10 : libmesh_assert(it != new_proc_ids.end());
2735 2504 : node->processor_id() = it->second;
2736 : }
2737 :
2738 : // We should still have consistent nodal processor ids coming out of
2739 : // this algorithm, but if we're allowed to repartition the mesh then
2740 : // they should be canonically correct too.
2741 : #ifdef DEBUG
2742 218 : libmesh_assert_valid_procids<Node>(mesh);
2743 : //if (repartition_all_nodes)
2744 : // libmesh_assert_canonical_node_procids(mesh);
2745 : #endif
2746 37675 : }
2747 :
2748 :
2749 :
2750 19630 : void Private::globally_renumber_nodes_and_elements (MeshBase & mesh)
2751 : {
2752 19630 : MeshCommunication().assign_global_indices(mesh);
2753 19630 : }
2754 :
2755 : } // namespace MeshTools
2756 :
2757 : } // namespace libMesh
|
