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