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/libmesh_config.h"
22 : #include "libmesh/libmesh_logging.h"
23 : #include "libmesh/boundary_info.h"
24 : #include "libmesh/distributed_mesh.h"
25 : #include "libmesh/elem.h"
26 : #include "libmesh/mesh_communication.h"
27 : #include "libmesh/mesh_serializer.h"
28 : #include "libmesh/parallel.h"
29 : #include "libmesh/partitioner.h"
30 : #include "libmesh/remote_elem.h"
31 : #include "libmesh/unstructured_mesh.h"
32 : #include "libmesh/elem_side_builder.h"
33 : #include "libmesh/utility.h"
34 :
35 : // TIMPI includes
36 : #include "timpi/parallel_sync.h"
37 :
38 : // C++ includes
39 : #include <iterator> // std::distance
40 : #include <algorithm> // std::max_element
41 :
42 : namespace
43 : {
44 :
45 : // Templated helper function for removing a subset of keys from a
46 : // multimap that further satisfy a given predicate on the
47 : // corresponding values.
48 : template <class Key, class T, class Pred>
49 1232403 : void erase_if(std::multimap<Key,T> & map, Key k, Pred pred)
50 : {
51 60344 : auto rng = map.equal_range(k);
52 60344 : auto it = rng.first;
53 1502306 : while (it != rng.second)
54 : {
55 252067 : if (pred(it->second))
56 65051 : it = map.erase(it);
57 : else
58 8807 : ++it;
59 : }
60 1232403 : }
61 :
62 : // Similar to the helper function above but doesn't take a key,
63 : // instead it applies the predicate to every value in the map.
64 : template <class Key, class T, class Pred>
65 2840 : void erase_if(std::multimap<Key,T> & map, Pred pred)
66 : {
67 80 : auto it = map.begin();
68 8335 : while (it != map.end())
69 : {
70 5495 : if (pred(it->second))
71 1947 : it = map.erase(it);
72 : else
73 196 : ++it;
74 : }
75 2840 : }
76 :
77 : // Helper func for renumber_id
78 : template <typename Map, typename T>
79 6358 : void renumber_name(Map & m, T old_id, T new_id)
80 : {
81 6358 : if (const auto it = std::as_const(m).find(old_id);
82 180 : it != m.end())
83 : {
84 4370 : m[new_id] = it->second;
85 4246 : m.erase(it);
86 : }
87 6358 : }
88 :
89 : }
90 :
91 : namespace libMesh
92 : {
93 :
94 :
95 :
96 : //------------------------------------------------------
97 : // BoundaryInfo static member initializations
98 : const boundary_id_type BoundaryInfo::invalid_id = -123;
99 :
100 :
101 :
102 : //------------------------------------------------------
103 : // BoundaryInfo functions
104 442693 : BoundaryInfo::BoundaryInfo(MeshBase & m) :
105 : ParallelObject(m.comm()),
106 416753 : _mesh (&m),
107 468609 : _children_on_boundary(false)
108 : {
109 442693 : }
110 :
111 24684 : BoundaryInfo & BoundaryInfo::operator=(const BoundaryInfo & other_boundary_info)
112 : {
113 : // Overwrite any preexisting boundary info
114 24684 : this->clear();
115 :
116 : /**
117 : * We're going to attempt to pull _new_ pointers out of the mesh
118 : * assigned to this boundary info.
119 : *
120 : * This will only work if the mesh assigned to this BoundaryInfo is
121 : * the same mesh object as other_boundary_info _or_ was constructed
122 : * in exactly the same way (or constructed as a copy, or a refined
123 : * copy without renumbering, etc.).
124 : */
125 :
126 : // Copy node boundary info
127 485927 : for (const auto & [node, bid] : other_boundary_info._boundary_node_id)
128 461243 : _boundary_node_id.emplace(_mesh->node_ptr(node->id()), bid);
129 :
130 : // Copy edge boundary info
131 24798 : for (const auto & [elem, id_pair] : other_boundary_info._boundary_edge_id)
132 114 : _boundary_edge_id.emplace(_mesh->elem_ptr(elem->id()), id_pair);
133 :
134 : // Copy shellface boundary info
135 24728 : for (const auto & [elem, id_pair] : other_boundary_info._boundary_shellface_id)
136 44 : _boundary_shellface_id.emplace(_mesh->elem_ptr(elem->id()), id_pair);
137 :
138 : // Copy side boundary info
139 542427 : for (const auto & [elem, id_pair] : other_boundary_info._boundary_side_id)
140 517743 : _boundary_side_id.emplace(_mesh->elem_ptr(elem->id()), id_pair);
141 :
142 24684 : _children_on_boundary = other_boundary_info._children_on_boundary;
143 :
144 820 : _boundary_ids = other_boundary_info._boundary_ids;
145 820 : _global_boundary_ids = other_boundary_info._global_boundary_ids;
146 820 : _side_boundary_ids = other_boundary_info._side_boundary_ids;
147 820 : _node_boundary_ids = other_boundary_info._node_boundary_ids;
148 820 : _edge_boundary_ids = other_boundary_info._edge_boundary_ids;
149 820 : _shellface_boundary_ids = other_boundary_info._shellface_boundary_ids;
150 :
151 820 : _ss_id_to_name = other_boundary_info._ss_id_to_name;
152 820 : _ns_id_to_name = other_boundary_info._ns_id_to_name;
153 820 : _es_id_to_name = other_boundary_info._es_id_to_name;
154 :
155 24684 : return *this;
156 : }
157 :
158 :
159 15238 : bool BoundaryInfo::operator==(const BoundaryInfo & other_boundary_info) const
160 : {
161 769392 : for (const auto & [other_node, bid] : other_boundary_info._boundary_node_id)
162 : {
163 754154 : const Node * node = this->_mesh->query_node_ptr(other_node->id());
164 754154 : if (!node)
165 0 : return false;
166 754154 : if (!this->has_boundary_id(node, bid))
167 0 : return false;
168 : }
169 767758 : for (const auto & [node, bid] : this->_boundary_node_id)
170 : {
171 : const Node * other_node =
172 752696 : other_boundary_info._mesh->query_node_ptr(node->id());
173 752696 : if (!other_node)
174 8 : return false;
175 752555 : if (!other_boundary_info.has_boundary_id(other_node, bid))
176 2 : return false;
177 : }
178 :
179 188 : auto compare_edges = [&](const Elem * elem,
180 : const Elem * other_elem,
181 60 : unsigned short int edge)
182 : {
183 248 : if (!elem)
184 0 : return false;
185 248 : if (!other_elem)
186 0 : return false;
187 :
188 80 : std::vector<boundary_id_type> our_edges, other_edges;
189 248 : this->edge_boundary_ids(elem, edge, our_edges);
190 248 : other_boundary_info.edge_boundary_ids(other_elem, edge, other_edges);
191 288 : if (our_edges.size() != other_edges.size())
192 0 : return false;
193 :
194 248 : std::sort(our_edges.begin(), our_edges.end());
195 248 : std::sort(other_edges.begin(), other_edges.end());
196 496 : for (auto i : index_range(our_edges))
197 268 : if (our_edges[i] != other_edges[i])
198 0 : return false;
199 60 : return true;
200 15062 : };
201 :
202 15186 : for (const auto & [other_elem, edge_id_pair] : other_boundary_info._boundary_edge_id)
203 : {
204 124 : const Elem * elem = this->_mesh->query_elem_ptr(other_elem->id());
205 124 : if (!compare_edges(elem, other_elem, edge_id_pair.first))
206 0 : return false;
207 : }
208 :
209 15186 : for (const auto & [elem, edge_id_pair] : this->_boundary_edge_id)
210 : {
211 124 : const Elem * other_elem = other_boundary_info._mesh->query_elem_ptr(elem->id());
212 124 : if (!compare_edges(elem, other_elem, edge_id_pair.first))
213 0 : return false;
214 : }
215 :
216 583588 : auto compare_sides = [&](const Elem * elem,
217 : const Elem * other_elem,
218 67872 : unsigned short int side)
219 : {
220 651460 : if (!elem)
221 0 : return false;
222 651460 : if (!other_elem)
223 0 : return false;
224 :
225 99680 : std::vector<boundary_id_type> our_sides, other_sides;
226 651460 : this->boundary_ids(elem, side, our_sides);
227 651460 : other_boundary_info.boundary_ids(other_elem, side, other_sides);
228 687524 : if (our_sides.size() != other_sides.size())
229 0 : return false;
230 :
231 651460 : std::sort(our_sides.begin(), our_sides.end());
232 651460 : std::sort(other_sides.begin(), other_sides.end());
233 1302920 : for (auto i : index_range(our_sides))
234 669492 : if (our_sides[i] != other_sides[i])
235 0 : return false;
236 67872 : return true;
237 15062 : };
238 :
239 340792 : for (const auto & [other_elem, side_id_pair] : other_boundary_info._boundary_side_id)
240 : {
241 325730 : const Elem * elem = this->_mesh->query_elem_ptr(other_elem->id());
242 325730 : if (!compare_sides(elem, other_elem, side_id_pair.first))
243 0 : return false;
244 : }
245 :
246 340792 : for (const auto & [elem, side_id_pair] : this->_boundary_side_id)
247 : {
248 325730 : const Elem * other_elem = other_boundary_info._mesh->query_elem_ptr(elem->id());
249 325730 : if (!compare_sides(elem, other_elem, side_id_pair.first))
250 0 : return false;
251 : }
252 :
253 72 : auto compare_shellfaces = [&](const Elem * elem,
254 : const Elem * other_elem,
255 24 : unsigned short int shellface)
256 : {
257 96 : if (!elem)
258 0 : return false;
259 96 : if (!other_elem)
260 0 : return false;
261 :
262 32 : std::vector<boundary_id_type> our_shellfaces, other_shellfaces;
263 96 : this->shellface_boundary_ids(elem, shellface, our_shellfaces);
264 96 : other_boundary_info.shellface_boundary_ids(other_elem, shellface, other_shellfaces);
265 112 : if (our_shellfaces.size() != other_shellfaces.size())
266 0 : return false;
267 :
268 96 : std::sort(our_shellfaces.begin(), our_shellfaces.end());
269 96 : std::sort(other_shellfaces.begin(), other_shellfaces.end());
270 192 : for (auto i : index_range(our_shellfaces))
271 104 : if (our_shellfaces[i] != other_shellfaces[i])
272 0 : return false;
273 24 : return true;
274 15062 : };
275 :
276 15110 : for (const auto & [other_elem, shellface_id_pair] : other_boundary_info._boundary_shellface_id)
277 : {
278 48 : const Elem * elem = this->_mesh->query_elem_ptr(other_elem->id());
279 48 : if (!compare_shellfaces(elem, other_elem, shellface_id_pair.first))
280 0 : return false;
281 : }
282 :
283 15110 : for (const auto & [elem, shellface_id_pair] : this->_boundary_shellface_id)
284 : {
285 48 : const Elem * other_elem = other_boundary_info._mesh->query_elem_ptr(elem->id());
286 48 : if (!compare_shellfaces(elem, other_elem, shellface_id_pair.first))
287 0 : return false;
288 : }
289 :
290 15062 : if (_children_on_boundary != other_boundary_info._children_on_boundary)
291 0 : return false;
292 :
293 90084 : auto compare_sets = [](const auto & set1, const auto & set2)
294 : {
295 90084 : if (set1.size() != set2.size())
296 0 : return false;
297 348728 : for (boundary_id_type bid : set1)
298 12502 : if (!set2.count(bid))
299 0 : return false;
300 :
301 87600 : return true;
302 : };
303 :
304 17206 : if (!compare_sets(_boundary_ids,
305 29052 : other_boundary_info._boundary_ids) ||
306 15062 : !compare_sets(_global_boundary_ids,
307 15464 : other_boundary_info._global_boundary_ids) ||
308 14990 : !compare_sets(_edge_boundary_ids,
309 15392 : other_boundary_info._edge_boundary_ids) ||
310 14990 : !compare_sets(_node_boundary_ids,
311 15392 : other_boundary_info._node_boundary_ids) ||
312 14990 : !compare_sets(_shellface_boundary_ids,
313 31526 : other_boundary_info._shellface_boundary_ids) ||
314 14990 : !compare_sets(_side_boundary_ids,
315 14990 : other_boundary_info._side_boundary_ids))
316 0 : return false;
317 :
318 44898 : auto compare_maps = [](const auto & map1, const auto & map2)
319 : {
320 44898 : if (map1.size() != map2.size())
321 0 : return false;
322 176096 : for (const auto & pair : map1)
323 258824 : if (!map2.count(pair.first) ||
324 131234 : map2.at(pair.first) != pair.second)
325 0 : return false;
326 :
327 3216 : return true;
328 : };
329 :
330 17134 : if (!compare_maps(_ss_id_to_name,
331 28872 : other_boundary_info._ss_id_to_name) ||
332 14954 : !compare_maps(_ns_id_to_name,
333 29944 : other_boundary_info._ns_id_to_name) ||
334 14954 : !compare_maps(_es_id_to_name,
335 14954 : other_boundary_info._es_id_to_name))
336 36 : return false;
337 :
338 1072 : return true;
339 : }
340 :
341 :
342 833573 : BoundaryInfo::~BoundaryInfo() = default;
343 :
344 :
345 :
346 1293450 : void BoundaryInfo::clear()
347 : {
348 38520 : _boundary_node_id.clear();
349 38520 : _boundary_side_id.clear();
350 38520 : _boundary_edge_id.clear();
351 38520 : _boundary_shellface_id.clear();
352 38520 : _boundary_ids.clear();
353 38520 : _side_boundary_ids.clear();
354 38520 : _node_boundary_ids.clear();
355 38520 : _edge_boundary_ids.clear();
356 38520 : _shellface_boundary_ids.clear();
357 38520 : _ss_id_to_name.clear();
358 38520 : _ns_id_to_name.clear();
359 38520 : _es_id_to_name.clear();
360 1293450 : }
361 :
362 :
363 :
364 603028 : void BoundaryInfo::regenerate_id_sets()
365 : {
366 32714 : const auto old_ss_id_to_name = _ss_id_to_name;
367 32714 : const auto old_ns_id_to_name = _ns_id_to_name;
368 32714 : const auto old_es_id_to_name = _es_id_to_name;
369 :
370 : // Clear the old caches
371 16357 : _boundary_ids.clear();
372 16357 : _side_boundary_ids.clear();
373 16357 : _node_boundary_ids.clear();
374 16357 : _edge_boundary_ids.clear();
375 16357 : _shellface_boundary_ids.clear();
376 32690 : _ss_id_to_name.clear();
377 32690 : _ns_id_to_name.clear();
378 32690 : _es_id_to_name.clear();
379 :
380 : // Loop over id maps to regenerate each set.
381 13162941 : for (const auto & pr : _boundary_node_id)
382 : {
383 12559913 : const boundary_id_type id = pr.second;
384 11971780 : _boundary_ids.insert(id);
385 11971780 : _node_boundary_ids.insert(id);
386 12559913 : if (const auto it = old_ns_id_to_name.find(id);
387 588181 : it != old_ns_id_to_name.end())
388 11055910 : _ns_id_to_name.emplace(id, it->second);
389 : }
390 :
391 605385 : for (const auto & pr : _boundary_edge_id)
392 : {
393 2357 : const boundary_id_type id = pr.second.second;
394 2211 : _boundary_ids.insert(id);
395 2211 : _edge_boundary_ids.insert(id);
396 2357 : if (const auto it = old_es_id_to_name.find(id);
397 146 : it != old_es_id_to_name.end())
398 416 : _es_id_to_name.emplace(id, it->second);
399 : }
400 :
401 6296841 : for (const auto & pr : _boundary_side_id)
402 : {
403 5693813 : const boundary_id_type id = pr.second.second;
404 5416383 : _boundary_ids.insert(id);
405 5416383 : _side_boundary_ids.insert(id);
406 5693813 : if (const auto it = old_ss_id_to_name.find(id);
407 278210 : it != old_ss_id_to_name.end())
408 4306648 : _ss_id_to_name.emplace(id, it->second);
409 : }
410 :
411 729076 : for (const auto & pr : _boundary_shellface_id)
412 : {
413 126048 : const boundary_id_type id = pr.second.second;
414 112732 : _boundary_ids.insert(id);
415 112732 : _shellface_boundary_ids.insert(id);
416 : }
417 :
418 : // Handle global data
419 16357 : libmesh_assert(_mesh);
420 603028 : if (!_mesh->is_serial())
421 : {
422 494786 : _communicator.set_union(_ss_id_to_name);
423 494786 : _communicator.set_union(_ns_id_to_name);
424 494786 : _communicator.set_union(_es_id_to_name);
425 : }
426 :
427 603028 : this->synchronize_global_id_set();
428 603028 : }
429 :
430 :
431 :
432 603099 : void BoundaryInfo::synchronize_global_id_set()
433 : {
434 : // Handle global data
435 16359 : _global_boundary_ids = _boundary_ids;
436 16359 : libmesh_assert(_mesh);
437 603099 : if (!_mesh->is_serial())
438 494850 : _communicator.set_union(_global_boundary_ids);
439 :
440 603099 : _mesh->_preparation.has_boundary_id_sets = true;
441 603099 : }
442 :
443 :
444 :
445 639 : void BoundaryInfo::sync (UnstructuredMesh & boundary_mesh)
446 : {
447 36 : std::set<boundary_id_type> request_boundary_ids(_boundary_ids);
448 621 : request_boundary_ids.insert(invalid_id);
449 639 : if (!_mesh->is_serial())
450 576 : this->comm().set_union(request_boundary_ids);
451 :
452 639 : this->sync(request_boundary_ids,
453 : boundary_mesh);
454 639 : }
455 :
456 :
457 :
458 2130 : void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_ids,
459 : UnstructuredMesh & boundary_mesh)
460 : {
461 : // Call the 3 argument version of this function with a dummy value for the third set.
462 120 : std::set<subdomain_id_type> subdomains_relative_to;
463 2070 : subdomains_relative_to.insert(Elem::invalid_subdomain_id);
464 :
465 2130 : this->sync(requested_boundary_ids,
466 : boundary_mesh,
467 : subdomains_relative_to);
468 2130 : }
469 :
470 :
471 :
472 2556 : void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_ids,
473 : UnstructuredMesh & boundary_mesh,
474 : const std::set<subdomain_id_type> & subdomains_relative_to)
475 : {
476 216 : LOG_SCOPE("sync()", "BoundaryInfo");
477 :
478 2556 : boundary_mesh.clear();
479 :
480 : /**
481 : * Deleting 0 elements seems weird, but it's better encapsulating
482 : * than exposing a set_is_serial(false) capability that might be
483 : * easily misused.
484 : */
485 2556 : if (!_mesh->is_serial())
486 2304 : boundary_mesh.delete_remote_elements();
487 :
488 : /**
489 : * If the boundary_mesh is still serial, that means we *can't*
490 : * parallelize it, so to make sure we can construct it in full on
491 : * every processor we'll serialize the interior mesh.
492 : *
493 : * We'll use a temporary MeshSerializer here, but as soon as we
494 : * unserialize we'll be turning a bunch of interior_parent() links
495 : * into dangling pointers, and it won't be easy to tell which. So
496 : * we'll keep around a distributed copy for that case, and query it
497 : * to fix up interior_parent() links as necessary.
498 : *
499 : * We'll also need to make sure to unserialize the mesh *before* we
500 : * prepare the boundary mesh for use, or the prepare_for_use() call
501 : * on a refined boundary mesh will happily notice that it can find
502 : * and restore some refined elements' interior_parent pointers, not
503 : * knowing that those interior parents are about to go remote.
504 : */
505 2484 : std::unique_ptr<MeshBase> mesh_copy;
506 2556 : if (boundary_mesh.is_serial() && !_mesh->is_serial())
507 896 : mesh_copy = _mesh->clone();
508 :
509 : auto serializer = std::make_unique<MeshSerializer>
510 2628 : (const_cast<MeshBase &>(*_mesh), boundary_mesh.is_serial());
511 :
512 : /**
513 : * Re-create the boundary mesh.
514 : */
515 :
516 2556 : boundary_mesh.set_n_partitions() = _mesh->n_partitions();
517 :
518 144 : std::map<dof_id_type, dof_id_type> node_id_map;
519 :
520 2556 : this->_find_id_maps(requested_boundary_ids, 0, &node_id_map, 0, nullptr, subdomains_relative_to);
521 :
522 : // Let's add all the boundary nodes we found to the boundary mesh
523 909642 : for (const auto & node : _mesh->node_ptr_range())
524 : {
525 477729 : dof_id_type node_id = node->id();
526 25428 : if (node_id_map.count(node_id))
527 : {
528 62649 : boundary_mesh.add_point(*node, node_id_map[node_id], node->processor_id());
529 :
530 : // Copy over all the node's boundary IDs to boundary_mesh
531 6788 : std::vector<boundary_id_type> node_boundary_ids;
532 62649 : this->boundary_ids(node, node_boundary_ids);
533 80661 : for (const auto & node_bid : node_boundary_ids)
534 18012 : boundary_mesh.get_boundary_info().add_node(node_id_map[node_id], node_bid);
535 : }
536 2412 : }
537 :
538 : // Add the elements. When syncing a boundary mesh, we also store the
539 : // parent side ids in addition to the interior_parent pointers,
540 : // since this information is frequently needed on boundary meshes.
541 2556 : this->add_elements(requested_boundary_ids,
542 : boundary_mesh,
543 : subdomains_relative_to,
544 : /*store_parent_side_ids=*/true);
545 :
546 : // The new elements are currently using the interior mesh's nodes;
547 : // we want them to use the boundary mesh's nodes instead.
548 :
549 : // This side's Node pointers still point to the nodes of the original mesh.
550 : // We need to re-point them to the boundary mesh's nodes! Since we copied *ALL* of
551 : // the original mesh's nodes over, we should be guaranteed to have the same ordering.
552 46026 : for (auto & new_elem : boundary_mesh.element_ptr_range())
553 : {
554 136142 : for (auto nn : new_elem->node_index_range())
555 : {
556 : // Get the correct node pointer, based on the id()
557 : Node * new_node =
558 119963 : boundary_mesh.node_ptr(node_id_map[new_elem->node_id(nn)]);
559 :
560 : // sanity check: be sure that the new Node exists and its
561 : // global id really matches
562 6934 : libmesh_assert (new_node);
563 6934 : libmesh_assert_equal_to (new_node->id(),
564 : node_id_map[new_elem->node_id(nn)]);
565 :
566 : // Assign the new node pointer
567 113029 : new_elem->set_node(nn, new_node);
568 : }
569 2412 : }
570 :
571 : // The new elements might have interior parent pointers aimed at
572 : // _mesh elements which are about to go remote, and we don't to
573 : // leave those pointers dangling. Fix them up if needed.
574 2556 : if (mesh_copy.get())
575 : {
576 7552 : for (auto & new_elem : boundary_mesh.element_ptr_range())
577 : {
578 : const dof_id_type interior_parent_id =
579 3328 : new_elem->interior_parent()->id();
580 :
581 3328 : if (!mesh_copy->query_elem_ptr(interior_parent_id))
582 : new_elem->set_interior_parent
583 2318 : (const_cast<RemoteElem *>(remote_elem));
584 448 : }
585 : }
586 :
587 : // Don't repartition this mesh; we want it to stay in sync with the
588 : // interior partitioning.
589 2556 : boundary_mesh.partitioner().reset(nullptr);
590 :
591 : // Deserialize the interior mesh before the boundary mesh
592 : // prepare_for_use() can come to erroneous conclusions about which
593 : // of its elements are semilocal
594 72 : serializer.reset();
595 :
596 : // Make boundary_mesh nodes and elements contiguous
597 2556 : boundary_mesh.prepare_for_use();
598 :
599 : // and finally distribute element partitioning to the nodes
600 2556 : Partitioner::set_node_processor_ids(boundary_mesh);
601 2556 : }
602 :
603 :
604 0 : void BoundaryInfo::get_side_and_node_maps (UnstructuredMesh & boundary_mesh,
605 : std::map<dof_id_type, dof_id_type> & node_id_map,
606 : std::map<dof_id_type, unsigned char> & side_id_map,
607 : Real tolerance)
608 : {
609 0 : LOG_SCOPE("get_side_and_node_maps()", "BoundaryInfo");
610 :
611 0 : node_id_map.clear();
612 0 : side_id_map.clear();
613 :
614 : // For building element sides without extraneous allocation
615 0 : ElemSideBuilder side_builder;
616 : // Pull objects out of the loop to reduce heap operations
617 : const Elem * interior_parent_side;
618 :
619 0 : for (const auto & boundary_elem : boundary_mesh.active_element_ptr_range())
620 : {
621 0 : const Elem * interior_parent = boundary_elem->interior_parent();
622 :
623 : // Find out which side of interior_parent boundary_elem corresponds to.
624 : // Use distance between average vertex location as a way to check.
625 0 : unsigned char interior_parent_side_index = 0;
626 0 : bool found_matching_sides = false;
627 0 : for (auto side : interior_parent->side_index_range())
628 : {
629 0 : interior_parent_side = &side_builder(*interior_parent, side);
630 0 : Real va_distance = (boundary_elem->vertex_average() - interior_parent_side->vertex_average()).norm();
631 :
632 0 : if (va_distance < (tolerance * boundary_elem->hmin()))
633 : {
634 0 : interior_parent_side_index = cast_int<unsigned char>(side);
635 0 : found_matching_sides = true;
636 0 : break;
637 : }
638 : }
639 :
640 0 : libmesh_error_msg_if(!found_matching_sides, "No matching side found within the specified tolerance");
641 :
642 0 : side_id_map[boundary_elem->id()] = interior_parent_side_index;
643 :
644 0 : for (auto local_node_index : boundary_elem->node_index_range())
645 : {
646 0 : dof_id_type boundary_node_id = boundary_elem->node_id(local_node_index);
647 0 : dof_id_type interior_node_id = interior_parent_side->node_id(local_node_index);
648 :
649 0 : node_id_map[interior_node_id] = boundary_node_id;
650 : }
651 0 : }
652 0 : }
653 :
654 :
655 :
656 781 : void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_boundary_ids,
657 : UnstructuredMesh & boundary_mesh,
658 : bool store_parent_side_ids,
659 : const std::vector<subdomain_id_type> & new_subdomain_ids)
660 : {
661 : // Call the 3 argument version of this function with a dummy value for the third arg.
662 44 : std::set<subdomain_id_type> subdomains_relative_to;
663 759 : subdomains_relative_to.insert(Elem::invalid_subdomain_id);
664 :
665 781 : this->add_elements(requested_boundary_ids,
666 : boundary_mesh,
667 : subdomains_relative_to,
668 : store_parent_side_ids,
669 : new_subdomain_ids);
670 781 : }
671 :
672 :
673 :
674 3337 : void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_boundary_ids,
675 : UnstructuredMesh & boundary_mesh,
676 : const std::set<subdomain_id_type> & subdomains_relative_to,
677 : bool store_parent_side_ids,
678 : const std::vector<subdomain_id_type> & new_subdomain_ids)
679 : {
680 94 : libmesh_assert(new_subdomain_ids.empty() ||
681 : new_subdomain_ids.size() == 1 ||
682 : new_subdomain_ids.size() == requested_boundary_ids.size());
683 282 : LOG_SCOPE("add_elements()", "BoundaryInfo");
684 :
685 : // We're not prepared to mix serial and distributed meshes in this
686 : // method, so make sure their statuses match from the start.
687 : //
688 : // Specifically test *is_serial* here - we can handle a mix of
689 : // ReplicatedMesh and serialized DistributedMesh.
690 94 : libmesh_assert_equal_to(_mesh->is_serial(),
691 : boundary_mesh.is_serial());
692 :
693 : // If the boundary mesh already has interior pointers pointing at
694 : // elements in a third mesh then we're in trouble
695 94 : libmesh_assert(&boundary_mesh.interior_mesh() == &boundary_mesh ||
696 : &boundary_mesh.interior_mesh() == _mesh);
697 :
698 : // And now we're going to add interior pointers to elements from
699 : // this mesh
700 3337 : boundary_mesh.set_interior_mesh(*_mesh);
701 :
702 188 : std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> side_id_map;
703 3337 : this->_find_id_maps(requested_boundary_ids,
704 : 0,
705 : nullptr,
706 3337 : boundary_mesh.max_elem_id(),
707 : &side_id_map,
708 : subdomains_relative_to);
709 :
710 : // We have to add sides *outside* any element loop, because if
711 : // boundary_mesh and _mesh are the same then those additions can
712 : // invalidate our element iterators. So we just use the element
713 : // loop to make a list of sides to add.
714 : typedef std::vector<std::tuple<dof_id_type, unsigned char, boundary_id_type>>
715 : side_container;
716 188 : side_container sides_to_add;
717 :
718 168002 : for (const auto & elem : _mesh->element_ptr_range())
719 : {
720 : // If the subdomains_relative_to container has the
721 : // invalid_subdomain_id, we fall back on the "old" behavior of
722 : // adding sides regardless of this Elem's subdomain. Otherwise,
723 : // if the subdomains_relative_to container doesn't contain the
724 : // current Elem's subdomain_id(), we won't add any sides from
725 : // it.
726 10476 : if (!subdomains_relative_to.count(Elem::invalid_subdomain_id) &&
727 7232 : !subdomains_relative_to.count(elem->subdomain_id()))
728 5212 : continue;
729 :
730 461907 : for (auto s : elem->side_index_range())
731 : {
732 22816 : bool add_this_side = false;
733 376920 : boundary_id_type triggering_bcid = invalid_id;
734 :
735 : // Find all the boundary side ids for this Elem side.
736 45632 : std::vector<boundary_id_type> bcids;
737 376920 : this->boundary_ids(elem, s, bcids);
738 :
739 404530 : for (const boundary_id_type bcid : bcids)
740 : {
741 : // if the user wants this id, we want this side
742 2868 : if (requested_boundary_ids.count(bcid))
743 : {
744 1512 : add_this_side = true;
745 24512 : triggering_bcid = bcid;
746 24512 : break;
747 : }
748 : }
749 :
750 : // We may still want to add this side if the user called
751 : // sync() with no requested_boundary_ids. This corresponds
752 : // to the "old" style of calling sync() in which the entire
753 : // boundary was copied to the BoundaryMesh, and handles the
754 : // case where elements on the geometric boundary are not in
755 : // any sidesets.
756 226732 : if (requested_boundary_ids.count(invalid_id) &&
757 203916 : elem->neighbor_ptr(s) == nullptr)
758 560 : add_this_side = true;
759 :
760 368717 : if (add_this_side)
761 24512 : sides_to_add.emplace_back(elem->id(), s, triggering_bcid);
762 : }
763 3149 : }
764 :
765 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
766 3337 : unique_id_type old_max_unique_id = boundary_mesh.parallel_max_unique_id();
767 : #endif
768 :
769 : // Add an "extra" integer for storing the side index of the parent
770 : // Elem which each boundary Elem corresponds to. We do this once
771 : // before any Elems have been added.
772 3337 : unsigned int parent_side_index_tag = store_parent_side_ids ?
773 5821 : boundary_mesh.add_elem_integer("parent_side_index") : libMesh::invalid_uint;
774 :
775 : // When new elements are assigned subdomain ids in one-to-one
776 : // correspondence with the requested boundary ids, we need to map a
777 : // boundary id to its index in the (sorted) set. Copying the set into a
778 : // contiguous vector lets us do that index lookup in O(1) (after an
779 : // O(log N) Utility::binary_find()) rather than the O(N) std::distance()
780 : // between std::set iterators.
781 188 : std::vector<boundary_id_type> requested_boundary_ids_vec;
782 3431 : if (new_subdomain_ids.size() > 1)
783 69 : requested_boundary_ids_vec.assign(requested_boundary_ids.begin(),
784 : requested_boundary_ids.end());
785 :
786 27849 : for (const auto & [elem_id, s, triggering_bcid] : sides_to_add)
787 : {
788 24512 : Elem * elem = _mesh->elem_ptr(elem_id);
789 :
790 26024 : std::unique_ptr<Elem> side = elem->build_side_ptr(s);
791 :
792 26024 : side->processor_id() = elem->processor_id();
793 :
794 1512 : const std::pair<dof_id_type, unsigned char> side_pair(elem_id, s);
795 :
796 1512 : libmesh_assert(side_id_map.count(side_pair));
797 :
798 24512 : const dof_id_type new_side_id = side_id_map[side_pair];
799 :
800 1512 : side->set_id(new_side_id);
801 :
802 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
803 24512 : side->set_unique_id(old_max_unique_id + new_side_id);
804 : #endif
805 :
806 : // Add the side
807 26024 : Elem * new_elem = boundary_mesh.add_elem(std::move(side));
808 :
809 : // If requested, new_elem gets an "extra" integer equal to the
810 : // side id "s" of the interior_parent it corresponds to.
811 24512 : if (store_parent_side_ids)
812 21803 : new_elem->set_extra_integer(parent_side_index_tag, s);
813 :
814 : // Assign subdomain id if requested.
815 26024 : if (new_subdomain_ids.size() == 1)
816 370 : new_elem->subdomain_id() = new_subdomain_ids[0];
817 24142 : else if (new_subdomain_ids.size() > 1)
818 : {
819 161 : auto it = Utility::binary_find(requested_boundary_ids_vec.begin(),
820 : requested_boundary_ids_vec.end(),
821 24 : triggering_bcid);
822 12 : libmesh_assert(it != requested_boundary_ids_vec.end());
823 185 : new_elem->subdomain_id() =
824 24 : new_subdomain_ids[std::distance(requested_boundary_ids_vec.begin(), it)];
825 : }
826 :
827 : #ifdef LIBMESH_ENABLE_AMR
828 3024 : new_elem->set_refinement_flag(elem->refinement_flag());
829 3024 : new_elem->set_p_refinement_flag(elem->p_refinement_flag());
830 :
831 : // Set parent links
832 26024 : if (elem->parent())
833 : {
834 672 : const std::pair<dof_id_type, unsigned char> parent_side_pair(elem->parent()->id(), s);
835 :
836 336 : libmesh_assert(side_id_map.count(parent_side_pair));
837 :
838 4432 : Elem * side_parent = boundary_mesh.elem_ptr(side_id_map[parent_side_pair]);
839 :
840 336 : libmesh_assert(side_parent);
841 :
842 672 : new_elem->set_parent(side_parent);
843 :
844 : // Figuring out which child we are of our parent
845 : // is a trick. Due to libMesh child numbering
846 : // conventions, if we are an element on a vertex,
847 : // then we share that vertex with our parent, with
848 : // the same local index.
849 336 : bool found_child = false;
850 13296 : for (auto v : make_range(new_elem->n_vertices()))
851 10208 : if (new_elem->node_ptr(v) == side_parent->node_ptr(v))
852 : {
853 4432 : side_parent->add_child(new_elem, v);
854 336 : found_child = true;
855 : }
856 :
857 : // If we don't share any vertex with our parent,
858 : // then we're the fourth child (index 3) of a
859 : // triangle.
860 4432 : if (!found_child)
861 : {
862 0 : libmesh_assert_equal_to (new_elem->n_vertices(), 3);
863 0 : side_parent->add_child(new_elem, 3);
864 : }
865 : }
866 :
867 : // Set remote_elem child links if necessary. Rather than
868 : // worrying about which interior child corresponds to which side
869 : // child we'll just set all null links to be remote and we'll
870 : // rely on our detection of actual semilocal children to
871 : // overwrite the links that shouldn't be remote.
872 1512 : if (elem->has_children())
873 16380 : for (auto c : make_range(elem->n_children()))
874 16486 : if (elem->child_ptr(c) == remote_elem &&
875 3382 : elem->is_child_on_side(c, s))
876 : {
877 6360 : for (auto sc : make_range(new_elem->n_children()))
878 2880 : if (!new_elem->raw_child_ptr(sc))
879 : new_elem->add_child
880 2720 : (const_cast<RemoteElem*>(remote_elem), sc);
881 : }
882 : #endif
883 :
884 24512 : new_elem->set_interior_parent (elem);
885 :
886 : // On non-local elements on DistributedMesh we might have
887 : // RemoteElem neighbor links to construct
888 24512 : if (!_mesh->is_serial() &&
889 15892 : (elem->processor_id() != this->processor_id()))
890 : {
891 10260 : const unsigned short n_nodes = elem->n_nodes();
892 :
893 10260 : const unsigned short bdy_n_sides = new_elem->n_sides();
894 10260 : const unsigned short bdy_n_nodes = new_elem->n_nodes();
895 :
896 : // Check every interior side for a RemoteElem
897 59696 : for (auto interior_side : elem->side_index_range())
898 : {
899 : // Might this interior side have a RemoteElem that
900 : // needs a corresponding Remote on a boundary side?
901 49436 : if (elem->neighbor_ptr(interior_side) != remote_elem)
902 41047 : continue;
903 :
904 : // Which boundary side?
905 9283 : for (unsigned short boundary_side = 0;
906 17672 : boundary_side != bdy_n_sides; ++boundary_side)
907 : {
908 : // Look for matching node points. This is safe in
909 : // *this* context.
910 0 : bool found_all_nodes = true;
911 53569 : for (unsigned short boundary_node = 0;
912 70473 : boundary_node != bdy_n_nodes; ++boundary_node)
913 : {
914 62852 : if (!new_elem->is_node_on_side(boundary_node,
915 0 : boundary_side))
916 36608 : continue;
917 :
918 0 : bool found_this_node = false;
919 286233 : for (unsigned short interior_node = 0;
920 312477 : interior_node != n_nodes; ++interior_node)
921 : {
922 303194 : if (!elem->is_node_on_side(interior_node,
923 0 : interior_side))
924 183756 : continue;
925 :
926 0 : if (new_elem->point(boundary_node) ==
927 0 : elem->point(interior_node))
928 : {
929 0 : found_this_node = true;
930 0 : break;
931 : }
932 : }
933 26244 : if (!found_this_node)
934 : {
935 0 : found_all_nodes = false;
936 0 : break;
937 : }
938 : }
939 :
940 16904 : if (found_all_nodes)
941 : {
942 : new_elem->set_neighbor
943 7621 : (boundary_side,
944 : const_cast<RemoteElem *>(remote_elem));
945 0 : break;
946 : }
947 : }
948 : }
949 : }
950 21488 : }
951 :
952 : // We haven't been bothering to keep unique ids consistent on ghost
953 : // elements or nodes, unless we're doing everything the same on
954 : // every processor.
955 3337 : if (!boundary_mesh.is_replicated())
956 2609 : MeshCommunication().make_node_unique_ids_parallel_consistent(boundary_mesh);
957 :
958 : // Make sure we didn't add ids inconsistently
959 : #ifdef DEBUG
960 : # ifdef LIBMESH_HAVE_RTTI
961 94 : DistributedMesh * parmesh = dynamic_cast<DistributedMesh *>(&boundary_mesh);
962 94 : if (parmesh)
963 14 : parmesh->libmesh_assert_valid_parallel_ids();
964 : # endif
965 : #endif
966 :
967 : // global containers are not synced
968 94 : boundary_mesh.unset_has_boundary_id_sets();
969 :
970 : // New elements were added; mark the affected preparation flags invalid.
971 94 : boundary_mesh.unset_has_neighbor_ptrs();
972 94 : boundary_mesh.unset_has_cached_elem_data();
973 94 : boundary_mesh.unset_has_reinit_ghosting_functors();
974 3337 : }
975 :
976 :
977 :
978 284570 : void BoundaryInfo::add_node(const dof_id_type node_id,
979 : const boundary_id_type id)
980 : {
981 284570 : const Node * node_ptr = _mesh->query_node_ptr(node_id);
982 :
983 : // The user could easily ask for an invalid node id, so let's throw
984 : // an easy-to-understand error message when this happens.
985 284570 : libmesh_error_msg_if(!node_ptr,
986 : "BoundaryInfo::add_node(): Could not retrieve pointer for node "
987 : << node_id
988 : << ", no boundary id was added.");
989 :
990 284570 : this->add_node (node_ptr, id);
991 284570 : }
992 :
993 :
994 :
995 54064523 : void BoundaryInfo::add_node(const Node * node,
996 : const boundary_id_type id)
997 : {
998 54064523 : libmesh_error_msg_if(id == invalid_id,
999 : "ERROR: You may not set a boundary ID of "
1000 : << invalid_id
1001 : << "\n That is reserved for internal use.");
1002 :
1003 : // Don't add the same ID twice
1004 80907578 : for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
1005 59335412 : if (pr.second == id)
1006 8766 : return;
1007 :
1008 558303 : _boundary_node_id.emplace(node, id);
1009 21013863 : _boundary_ids.insert(id);
1010 21013863 : _node_boundary_ids.insert(id); // Also add this ID to the set of node boundary IDs
1011 : }
1012 :
1013 :
1014 :
1015 6173299 : void BoundaryInfo::add_node(const Node * node,
1016 : const std::vector<boundary_id_type> & ids)
1017 : {
1018 6173299 : if (ids.empty())
1019 5977306 : return;
1020 :
1021 1584 : libmesh_assert(node);
1022 :
1023 : // Don't add the same ID twice
1024 1584 : auto bounds = _boundary_node_id.equal_range(node);
1025 :
1026 : // The entries in the ids vector may be non-unique. If we expected
1027 : // *lots* of ids, it might be fastest to construct a std::set from
1028 : // the entries, but for a small number of entries, which is more
1029 : // typical, it is probably faster to copy the vector and do sort+unique.
1030 : // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
1031 197577 : std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
1032 195993 : std::sort(unique_ids.begin(), unique_ids.end());
1033 : std::vector<boundary_id_type>::iterator new_end =
1034 195993 : std::unique(unique_ids.begin(), unique_ids.end());
1035 :
1036 486538 : for (auto & id : as_range(unique_ids.begin(), new_end))
1037 : {
1038 290545 : libmesh_error_msg_if(id == invalid_id,
1039 : "ERROR: You may not set a boundary ID of "
1040 : << invalid_id
1041 : << "\n That is reserved for internal use.");
1042 :
1043 2472 : bool already_inserted = false;
1044 466027 : for (const auto & pr : as_range(bounds))
1045 426409 : if (pr.second == id)
1046 : {
1047 1356 : already_inserted = true;
1048 1356 : break;
1049 : }
1050 290545 : if (already_inserted)
1051 249571 : continue;
1052 :
1053 1116 : _boundary_node_id.emplace(node, id);
1054 38502 : _boundary_ids.insert(id);
1055 38502 : _node_boundary_ids.insert(id); // Also add this ID to the set of node boundary IDs
1056 : }
1057 : }
1058 :
1059 :
1060 :
1061 0 : void BoundaryInfo::clear_boundary_node_ids()
1062 : {
1063 0 : _boundary_node_id.clear();
1064 0 : }
1065 :
1066 40336 : void BoundaryInfo::add_edge(const dof_id_type e,
1067 : const unsigned short int edge,
1068 : const boundary_id_type id)
1069 : {
1070 40336 : this->add_edge (_mesh->elem_ptr(e), edge, id);
1071 40336 : }
1072 :
1073 :
1074 :
1075 44611 : void BoundaryInfo::add_edge(const Elem * elem,
1076 : const unsigned short int edge,
1077 : const boundary_id_type id)
1078 : {
1079 1286 : libmesh_assert(elem);
1080 :
1081 : // Only add BCs for level-0 elements.
1082 1286 : libmesh_assert_equal_to (elem->level(), 0);
1083 :
1084 : // Only add BCs for edges that exist.
1085 1286 : libmesh_assert_less (edge, elem->n_edges());
1086 :
1087 44611 : libmesh_error_msg_if(id == invalid_id,
1088 : "ERROR: You may not set a boundary ID of "
1089 : << invalid_id
1090 : << "\n That is reserved for internal use.");
1091 :
1092 : // Don't add the same ID twice
1093 107572 : for (const auto & pr : as_range(_boundary_edge_id.equal_range(elem)))
1094 74561 : if (pr.second.first == edge &&
1095 11600 : pr.second.second == id)
1096 320 : return;
1097 :
1098 266858 : _boundary_edge_id.emplace(elem, std::make_pair(edge, id));
1099 32045 : _boundary_ids.insert(id);
1100 32045 : _edge_boundary_ids.insert(id); // Also add this ID to the set of edge boundary IDs
1101 : }
1102 :
1103 :
1104 :
1105 26668769 : void BoundaryInfo::add_edge(const Elem * elem,
1106 : const unsigned short int edge,
1107 : const std::vector<boundary_id_type> & ids)
1108 : {
1109 26668769 : if (ids.empty())
1110 26668769 : return;
1111 :
1112 0 : libmesh_assert(elem);
1113 :
1114 : // Only add BCs for level-0 elements.
1115 0 : libmesh_assert_equal_to (elem->level(), 0);
1116 :
1117 : // Only add BCs for edges that exist.
1118 0 : libmesh_assert_less (edge, elem->n_edges());
1119 :
1120 : // Don't add the same ID twice
1121 0 : auto bounds = _boundary_edge_id.equal_range(elem);
1122 :
1123 : // The entries in the ids vector may be non-unique. If we expected
1124 : // *lots* of ids, it might be fastest to construct a std::set from
1125 : // the entries, but for a small number of entries, which is more
1126 : // typical, it is probably faster to copy the vector and do sort+unique.
1127 : // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
1128 0 : std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
1129 0 : std::sort(unique_ids.begin(), unique_ids.end());
1130 : std::vector<boundary_id_type>::iterator new_end =
1131 0 : std::unique(unique_ids.begin(), unique_ids.end());
1132 :
1133 0 : for (auto & id : as_range(unique_ids.begin(), new_end))
1134 : {
1135 0 : libmesh_error_msg_if(id == invalid_id,
1136 : "ERROR: You may not set a boundary ID of "
1137 : << invalid_id
1138 : << "\n That is reserved for internal use.");
1139 :
1140 0 : bool already_inserted = false;
1141 0 : for (const auto & pr : as_range(bounds))
1142 0 : if (pr.second.first == edge &&
1143 0 : pr.second.second == id)
1144 : {
1145 0 : already_inserted = true;
1146 0 : break;
1147 : }
1148 0 : if (already_inserted)
1149 0 : continue;
1150 :
1151 0 : _boundary_edge_id.emplace(elem, std::make_pair(edge, id));
1152 0 : _boundary_ids.insert(id);
1153 0 : _edge_boundary_ids.insert(id); // Also add this ID to the set of edge boundary IDs
1154 : }
1155 : }
1156 :
1157 :
1158 :
1159 39936 : void BoundaryInfo::add_shellface(const dof_id_type e,
1160 : const unsigned short int shellface,
1161 : const boundary_id_type id)
1162 : {
1163 39936 : this->add_shellface (_mesh->elem_ptr(e), shellface, id);
1164 39936 : }
1165 :
1166 :
1167 :
1168 537274 : void BoundaryInfo::add_shellface(const Elem * elem,
1169 : const unsigned short int shellface,
1170 : const boundary_id_type id)
1171 : {
1172 13316 : libmesh_assert(elem);
1173 :
1174 : // Only add BCs for level-0 elements.
1175 13316 : libmesh_assert_equal_to (elem->level(), 0);
1176 :
1177 : // Shells only have 2 faces
1178 13316 : libmesh_assert_less(shellface, 2);
1179 :
1180 537274 : libmesh_error_msg_if(id == invalid_id,
1181 : "ERROR: You may not set a boundary ID of "
1182 : << invalid_id
1183 : << "\n That is reserved for internal use.");
1184 :
1185 : // Don't add the same ID twice
1186 805840 : for (const auto & pr : as_range(_boundary_shellface_id.equal_range(elem)))
1187 314230 : if (pr.second.first == shellface &&
1188 45664 : pr.second.second == id)
1189 0 : return;
1190 :
1191 491610 : _boundary_shellface_id.emplace(elem, std::make_pair(shellface, id));
1192 478294 : _boundary_ids.insert(id);
1193 478294 : _shellface_boundary_ids.insert(id); // Also add this ID to the set of shellface boundary IDs
1194 : }
1195 :
1196 :
1197 :
1198 8707090 : void BoundaryInfo::add_shellface(const Elem * elem,
1199 : const unsigned short int shellface,
1200 : const std::vector<boundary_id_type> & ids)
1201 : {
1202 8707090 : if (ids.empty())
1203 8707090 : return;
1204 :
1205 0 : libmesh_assert(elem);
1206 :
1207 : // Only add BCs for level-0 elements.
1208 0 : libmesh_assert_equal_to (elem->level(), 0);
1209 :
1210 : // Shells only have 2 faces
1211 0 : libmesh_assert_less(shellface, 2);
1212 :
1213 : // Don't add the same ID twice
1214 0 : auto bounds = _boundary_shellface_id.equal_range(elem);
1215 :
1216 : // The entries in the ids vector may be non-unique. If we expected
1217 : // *lots* of ids, it might be fastest to construct a std::set from
1218 : // the entries, but for a small number of entries, which is more
1219 : // typical, it is probably faster to copy the vector and do sort+unique.
1220 : // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
1221 0 : std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
1222 0 : std::sort(unique_ids.begin(), unique_ids.end());
1223 : std::vector<boundary_id_type>::iterator new_end =
1224 0 : std::unique(unique_ids.begin(), unique_ids.end());
1225 :
1226 0 : for (auto & id : as_range(unique_ids.begin(), new_end))
1227 : {
1228 0 : libmesh_error_msg_if(id == invalid_id,
1229 : "ERROR: You may not set a boundary ID of "
1230 : << invalid_id
1231 : << "\n That is reserved for internal use.");
1232 :
1233 0 : bool already_inserted = false;
1234 0 : for (const auto & pr : as_range(bounds))
1235 0 : if (pr.second.first == shellface &&
1236 0 : pr.second.second == id)
1237 : {
1238 0 : already_inserted = true;
1239 0 : break;
1240 : }
1241 0 : if (already_inserted)
1242 0 : continue;
1243 :
1244 0 : _boundary_shellface_id.emplace(elem, std::make_pair(shellface, id));
1245 0 : _boundary_ids.insert(id);
1246 0 : _shellface_boundary_ids.insert(id); // Also add this ID to the set of shellface boundary IDs
1247 : }
1248 : }
1249 :
1250 :
1251 107334 : void BoundaryInfo::add_side(const dof_id_type e,
1252 : const unsigned short int side,
1253 : const boundary_id_type id)
1254 : {
1255 107334 : this->add_side (_mesh->elem_ptr(e), side, id);
1256 107334 : }
1257 :
1258 :
1259 :
1260 25021407 : void BoundaryInfo::add_side(const Elem * elem,
1261 : const unsigned short int side,
1262 : const boundary_id_type id)
1263 : {
1264 223113 : libmesh_assert(elem);
1265 :
1266 : // Only add BCs for sides that exist.
1267 223113 : libmesh_assert_less (side, elem->n_sides());
1268 :
1269 25021407 : libmesh_error_msg_if(id == invalid_id, "ERROR: You may not set a boundary ID of "
1270 : << invalid_id
1271 : << "\n That is reserved for internal use.");
1272 :
1273 : // Don't add the same ID twice
1274 31930381 : for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
1275 23066207 : if (pr.second.first == side &&
1276 16158452 : pr.second.second == id)
1277 3546 : return;
1278 :
1279 : #ifdef LIBMESH_ENABLE_AMR
1280 : // Users try to mark boundary on child elements
1281 : // If this happens, we will allow users to remove
1282 : // side from child elements as well
1283 8864174 : if (elem->level())
1284 : {
1285 632 : _children_on_boundary = true;
1286 :
1287 : // Here we have to stop and check if we already have this boundary defined on the
1288 : // parent (if yes, no need to add)
1289 32 : std::vector<boundary_id_type> bd_ids;
1290 632 : this->boundary_ids(elem,side,bd_ids);
1291 :
1292 632 : if(std::find(bd_ids.begin(), bd_ids.end(), id) != bd_ids.end())
1293 167 : libmesh_not_implemented_msg("Trying to add boundary ID "
1294 : + std::to_string(id)
1295 : + " which already exists on the ancestors.");
1296 : }
1297 : #endif
1298 :
1299 8864139 : _boundary_side_id.emplace(elem, std::make_pair(side, id));
1300 8644574 : _boundary_ids.insert(id);
1301 8644574 : _side_boundary_ids.insert(id); // Also add this ID to the set of side boundary IDs
1302 : }
1303 :
1304 :
1305 :
1306 17996598 : void BoundaryInfo::add_side(const Elem * elem,
1307 : const unsigned short int side,
1308 : const std::vector<boundary_id_type> & ids)
1309 : {
1310 17996598 : if (ids.empty())
1311 16162256 : return;
1312 :
1313 52406 : libmesh_assert(elem);
1314 :
1315 : // Only add BCs for sides that exist.
1316 52406 : libmesh_assert_less (side, elem->n_sides());
1317 :
1318 : #ifdef LIBMESH_ENABLE_AMR
1319 : // Users try to mark boundary on child elements
1320 : // If this happens, we will allow users to remove
1321 : // side from child elements as well
1322 1834342 : if (elem->level())
1323 : {
1324 35 : _children_on_boundary = true;
1325 :
1326 : // Here we have to stop and check if we already have this boundary defined on the
1327 : // parent (if yes, no need to add)
1328 4 : std::vector<boundary_id_type> bd_ids;
1329 35 : this->boundary_ids(elem,side,bd_ids);
1330 :
1331 35 : for (const auto id : ids)
1332 35 : if(std::find(bd_ids.begin(), bd_ids.end(), id) != bd_ids.end())
1333 167 : libmesh_not_implemented_msg("Trying to add boundary ID "
1334 : + std::to_string(id)
1335 : + " which already exists on the ancestors.");
1336 : }
1337 : #endif
1338 :
1339 : // Don't add the same ID twice
1340 52404 : auto bounds = _boundary_side_id.equal_range(elem);
1341 :
1342 : // The entries in the ids vector may be non-unique. If we expected
1343 : // *lots* of ids, it might be fastest to construct a std::set from
1344 : // the entries, but for a small number of entries, which is more
1345 : // typical, it is probably faster to copy the vector and do sort+unique.
1346 : // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
1347 1886711 : std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
1348 1834307 : std::sort(unique_ids.begin(), unique_ids.end());
1349 : std::vector<boundary_id_type>::iterator new_end =
1350 1834307 : std::unique(unique_ids.begin(), unique_ids.end());
1351 :
1352 3668614 : for (auto & id : as_range(unique_ids.begin(), new_end))
1353 : {
1354 1834307 : libmesh_error_msg_if(id == invalid_id,
1355 : "ERROR: You may not set a boundary ID of "
1356 : << invalid_id
1357 : << "\n That is reserved for internal use.");
1358 :
1359 52404 : bool already_inserted = false;
1360 1882865 : for (const auto & pr : as_range(bounds))
1361 48558 : if (pr.second.first == side && pr.second.second == id)
1362 : {
1363 0 : already_inserted = true;
1364 0 : break;
1365 : }
1366 1834307 : if (already_inserted)
1367 0 : continue;
1368 :
1369 1834307 : _boundary_side_id.emplace(elem, std::make_pair(side, id));
1370 1781903 : _boundary_ids.insert(id);
1371 1781903 : _side_boundary_ids.insert(id); // Also add this ID to the set of side boundary IDs
1372 : }
1373 : }
1374 :
1375 :
1376 :
1377 1508619 : bool BoundaryInfo::has_boundary_id(const Node * const node,
1378 : const boundary_id_type id) const
1379 : {
1380 2191204 : for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
1381 2190049 : if (pr.second == id)
1382 61038 : return true;
1383 :
1384 66 : return false;
1385 : }
1386 :
1387 :
1388 :
1389 126005724 : void BoundaryInfo::boundary_ids (const Node * node,
1390 : std::vector<boundary_id_type> & vec_to_fill) const
1391 : {
1392 : // Clear out any previous contents
1393 29846972 : vec_to_fill.clear();
1394 :
1395 145581866 : for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
1396 19576142 : vec_to_fill.push_back(pr.second);
1397 126005724 : }
1398 :
1399 :
1400 :
1401 82614544 : unsigned int BoundaryInfo::n_boundary_ids(const Node * node) const
1402 : {
1403 158970 : auto pos = _boundary_node_id.equal_range(node);
1404 82773514 : return cast_int<unsigned int>(std::distance(pos.first, pos.second));
1405 : }
1406 :
1407 :
1408 :
1409 224689646 : void BoundaryInfo::edge_boundary_ids (const Elem * const elem,
1410 : const unsigned short int edge,
1411 : std::vector<boundary_id_type> & vec_to_fill) const
1412 : {
1413 22863974 : libmesh_assert(elem);
1414 :
1415 : // Clear out any previous contents
1416 22863974 : vec_to_fill.clear();
1417 :
1418 : // Only query BCs for edges that exist.
1419 22863974 : libmesh_assert_less (edge, elem->n_edges());
1420 :
1421 : // Only level-0 elements store BCs. If this is not a level-0
1422 : // element get its level-0 parent and infer the BCs.
1423 224689646 : const Elem * searched_elem = elem;
1424 : #ifdef LIBMESH_ENABLE_AMR
1425 224689646 : if (elem->level() != 0)
1426 : {
1427 : // Find all the sides that contain edge. If one of those is a boundary
1428 : // side, then this must be a boundary edge. In that case, we just use the
1429 : // top-level parent.
1430 11009277 : bool found_boundary_edge = false;
1431 72378393 : for (auto side : elem->side_index_range())
1432 : {
1433 59181020 : if (elem->is_edge_on_side(edge,side))
1434 : {
1435 19135244 : if (elem->neighbor_ptr(side) == nullptr)
1436 : {
1437 871984 : searched_elem = elem->top_parent ();
1438 652828 : found_boundary_edge = true;
1439 652828 : break;
1440 : }
1441 : }
1442 : }
1443 :
1444 11009277 : if (!found_boundary_edge)
1445 : {
1446 : // Child element is not on external edge, but it may have internal
1447 : // "boundary" IDs. We will walk up the tree, at each level checking that
1448 : // the current child is actually on the same edge of the parent that is
1449 : // currently being searched for (i.e. that was passed in as "edge").
1450 22135039 : while (searched_elem->parent() != nullptr)
1451 : {
1452 16305496 : const Elem * parent = searched_elem->parent();
1453 20720305 : if (parent->is_child_on_edge(parent->which_child_am_i(searched_elem), edge) == false)
1454 12218403 : return;
1455 8501902 : searched_elem = parent;
1456 : }
1457 : }
1458 : }
1459 : #endif
1460 :
1461 : // Check each element in the range to see if its edge matches the requested edge.
1462 212495475 : for (const auto & pr : as_range(_boundary_edge_id.equal_range(searched_elem)))
1463 24232 : if (pr.second.first == edge)
1464 2474 : vec_to_fill.push_back(pr.second.second);
1465 : }
1466 :
1467 :
1468 :
1469 97838610 : unsigned int BoundaryInfo::n_edge_boundary_ids (const Elem * const elem,
1470 : const unsigned short int edge) const
1471 : {
1472 1224768 : std::vector<boundary_id_type> ids;
1473 97838610 : this->edge_boundary_ids(elem, edge, ids);
1474 98451042 : return cast_int<unsigned int>(ids.size());
1475 : }
1476 :
1477 :
1478 :
1479 47977711 : void BoundaryInfo::raw_edge_boundary_ids (const Elem * const elem,
1480 : const unsigned short int edge,
1481 : std::vector<boundary_id_type> & vec_to_fill) const
1482 : {
1483 22091666 : libmesh_assert(elem);
1484 :
1485 : // Only query BCs for edges that exist.
1486 22091666 : libmesh_assert_less (edge, elem->n_edges());
1487 :
1488 : // Clear out any previous contents
1489 22091666 : vec_to_fill.clear();
1490 :
1491 : // Only level-0 elements store BCs.
1492 48760435 : if (elem->parent())
1493 10809304 : return;
1494 :
1495 : // Check each element in the range to see if its edge matches the requested edge.
1496 36415622 : for (const auto & pr : as_range(_boundary_edge_id.equal_range(elem)))
1497 3384 : if (pr.second.first == edge)
1498 282 : vec_to_fill.push_back(pr.second.second);
1499 : }
1500 :
1501 :
1502 :
1503 70205602 : void BoundaryInfo::shellface_boundary_ids (const Elem * const elem,
1504 : const unsigned short int shellface,
1505 : std::vector<boundary_id_type> & vec_to_fill) const
1506 : {
1507 8567474 : libmesh_assert(elem);
1508 :
1509 : // Shells only have 2 faces
1510 8567474 : libmesh_assert_less(shellface, 2);
1511 :
1512 : // Clear out any previous contents
1513 8567474 : vec_to_fill.clear();
1514 :
1515 : // Only level-0 elements store BCs. If this is not a level-0
1516 : // element get its level-0 parent and infer the BCs.
1517 70205602 : const Elem * searched_elem = elem;
1518 : #ifdef LIBMESH_ENABLE_AMR
1519 70205602 : if (elem->level() != 0)
1520 : {
1521 27951048 : while (searched_elem->parent() != nullptr)
1522 : {
1523 17328592 : const Elem * parent = searched_elem->parent();
1524 21437532 : searched_elem = parent;
1525 : }
1526 : }
1527 : #endif
1528 :
1529 : // Check each element in the range to see if its shellface matches the requested shellface.
1530 70749802 : for (const auto & pr : as_range(_boundary_shellface_id.equal_range(searched_elem)))
1531 544200 : if (pr.second.first == shellface)
1532 272196 : vec_to_fill.push_back(pr.second.second);
1533 70205602 : }
1534 :
1535 :
1536 :
1537 29921706 : unsigned int BoundaryInfo::n_shellface_boundary_ids (const Elem * const elem,
1538 : const unsigned short int shellface) const
1539 : {
1540 400384 : std::vector<boundary_id_type> ids;
1541 29921706 : this->shellface_boundary_ids(elem, shellface, ids);
1542 30125226 : return cast_int<unsigned int>(ids.size());
1543 : }
1544 :
1545 :
1546 :
1547 16768316 : void BoundaryInfo::raw_shellface_boundary_ids (const Elem * const elem,
1548 : const unsigned short int shellface,
1549 : std::vector<boundary_id_type> & vec_to_fill) const
1550 : {
1551 8327790 : libmesh_assert(elem);
1552 :
1553 : // Shells only have 2 faces
1554 8327790 : libmesh_assert_less(shellface, 2);
1555 :
1556 : // Clear out any previous contents
1557 8327790 : vec_to_fill.clear();
1558 :
1559 : // Only level-0 elements store BCs.
1560 17034880 : if (elem->parent())
1561 4752760 : return;
1562 :
1563 : // Check each element in the range to see if its shellface matches the requested shellface.
1564 11668222 : for (const auto & pr : as_range(_boundary_shellface_id.equal_range(elem)))
1565 66584 : if (pr.second.first == shellface)
1566 33292 : vec_to_fill.push_back(pr.second.second);
1567 : }
1568 :
1569 :
1570 :
1571 7089713 : bool BoundaryInfo::has_boundary_id(const Elem * const elem,
1572 : const unsigned short int side,
1573 : const boundary_id_type id) const
1574 : {
1575 534233 : std::vector<boundary_id_type> ids;
1576 7089713 : this->boundary_ids(elem, side, ids);
1577 7882225 : return (std::find(ids.begin(), ids.end(), id) != ids.end());
1578 : }
1579 :
1580 :
1581 82 : void BoundaryInfo::side_boundary_ids (const Elem * const elem,
1582 : std::vector<std::vector<boundary_id_type>> & vec_to_fill) const
1583 : {
1584 5 : libmesh_assert(elem);
1585 :
1586 : // Clear out any previous contents
1587 5 : vec_to_fill.clear();
1588 82 : auto num_sides = elem->n_sides();
1589 :
1590 : // No sides, no boundary ids
1591 82 : if (!num_sides)
1592 70 : return;
1593 :
1594 : // We are going to gather boundary ids for each side
1595 82 : vec_to_fill.resize(num_sides);
1596 : // In most cases only level-0 elements store BCs.
1597 : // In certain applications (such as time-dependent domains), however, children
1598 : // need to store BCs too. This case is covered with the _children_on_boundary
1599 : // flag.
1600 82 : const Elem * searched_elem = elem;
1601 :
1602 : #ifdef LIBMESH_ENABLE_AMR
1603 82 : if (elem->level() != 0)
1604 : {
1605 : // If we have children on the boundaries, we need to search for boundary IDs on the
1606 : // child and its ancestors too if they share the side.
1607 70 : if (_children_on_boundary)
1608 : {
1609 : // Loop over ancestors to check if they have boundary ids on the same side
1610 35 : std::vector<bool> search_on_side(elem->n_sides(), true);
1611 2 : bool keep_searching = true;
1612 105 : while (searched_elem && keep_searching)
1613 : {
1614 245 : for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
1615 : {
1616 875 : for (const auto side : make_range(elem->n_sides()))
1617 : // Here we need to check if the boundary id already exists
1618 752 : if (search_on_side[side] && pr.second.first == side &&
1619 672 : std::find(vec_to_fill[side].begin(), vec_to_fill[side].end(), pr.second.second) ==
1620 52 : vec_to_fill[side].end())
1621 111 : vec_to_fill[side].push_back(pr.second.second);
1622 : }
1623 :
1624 8 : const Elem * parent = searched_elem->parent();
1625 70 : const auto child_index = parent ? parent->which_child_am_i(searched_elem) : libMesh::invalid_uint;
1626 350 : for (const auto side : make_range(elem->n_sides()))
1627 : // If the parent doesn't exist or if the child is not on the correct side of the
1628 : // parent we are done checking the ancestors
1629 320 : if (search_on_side[side] &&
1630 156 : (!parent || parent->is_child_on_side(child_index, side) == false))
1631 16 : search_on_side[side] = false;
1632 :
1633 70 : searched_elem = parent;
1634 : // if found what we needed on all sides, exit
1635 70 : keep_searching = *std::max_element(search_on_side.begin(), search_on_side.end());
1636 : }
1637 :
1638 2 : return;
1639 : }
1640 :
1641 : // Children not on boundaries case.
1642 : // It could be that a children is interior to the parent (search_on_side = false will handle that)
1643 : // However, since no children on boundaries, we know that it's either the top parent or nothing
1644 35 : std::vector<bool> search_on_side(elem->n_sides(), true);
1645 175 : for (const auto side : make_range(elem->n_sides()))
1646 : {
1647 : // Reset the search level for each side
1648 8 : const Elem * searched_elem_for_side = elem;
1649 :
1650 : // If we don't have children on boundaries and we are on an external boundary,
1651 : // we will just use the top parent. search_on_side[side] = true works
1652 148 : if (elem->neighbor_ptr(side) == nullptr)
1653 66 : continue;
1654 : // Otherwise we loop over the ancestors and check if they have a different BC for us
1655 : else
1656 148 : while (searched_elem_for_side->parent() != nullptr)
1657 : {
1658 4 : const Elem * parent = searched_elem_for_side->parent();
1659 74 : if (search_on_side[side] && parent->is_child_on_side(parent->which_child_am_i(searched_elem_for_side), side) == false)
1660 8 : search_on_side[side] = false;
1661 4 : searched_elem_for_side = parent;
1662 : }
1663 : }
1664 : // Now search on the top parent, only if we need to (element is not deep inside the top parent)
1665 37 : if (*std::max_element(search_on_side.begin(), search_on_side.end()))
1666 175 : for (const auto & pr : as_range(_boundary_side_id.equal_range(elem->top_parent())))
1667 148 : if (search_on_side[pr.second.first])
1668 74 : vec_to_fill[pr.second.first].push_back(pr.second.second);
1669 2 : return;
1670 : }
1671 : #endif
1672 :
1673 : // Check each element in the range to see if its side matches the requested side.
1674 72 : for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
1675 65 : vec_to_fill[pr.second.first].push_back(pr.second.second);
1676 : }
1677 :
1678 35621035 : void BoundaryInfo::boundary_ids (const Elem * const elem,
1679 : const unsigned short int side,
1680 : std::vector<boundary_id_type> & vec_to_fill) const
1681 : {
1682 17569013 : libmesh_assert(elem);
1683 :
1684 : // Only query BCs for sides that exist.
1685 17569013 : libmesh_assert_less (side, elem->n_sides());
1686 :
1687 : // Clear out any previous contents
1688 17569013 : vec_to_fill.clear();
1689 :
1690 : // In most cases only level-0 elements store BCs.
1691 : // In certain applications (such as time-dependent domains), however, children
1692 : // need to store BCs too. This case is covered with the _children_on_boundary
1693 : // flag.
1694 35621035 : const Elem * searched_elem = elem;
1695 :
1696 : #ifdef LIBMESH_ENABLE_AMR
1697 :
1698 35621035 : if (elem->level() != 0)
1699 : {
1700 : // If we have children on the boundaries, we need to search for boundary IDs on the
1701 : // child and its ancestors too if they share the side.
1702 13393998 : if (_children_on_boundary)
1703 : {
1704 : // Loop over ancestors to check if they have boundary ids on the same side
1705 2756 : while (searched_elem)
1706 : {
1707 6567 : for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
1708 : // Here we need to check if the boundary id already exists
1709 4247 : if (pr.second.first == side &&
1710 2565 : std::find(vec_to_fill.begin(), vec_to_fill.end(), pr.second.second) ==
1711 2100 : vec_to_fill.end())
1712 767 : vec_to_fill.push_back(pr.second.second);
1713 :
1714 :
1715 1368 : const Elem * parent = searched_elem->parent();
1716 : // If the parent doesn't exist or if the child is not on the correct side of the
1717 : // parent we are done checking the ancestors
1718 2756 : if (!parent || parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
1719 11190588 : return;
1720 :
1721 1037 : searched_elem = parent;
1722 : }
1723 :
1724 0 : return;
1725 : }
1726 :
1727 : // If we don't have children on boundaries and we are on an external boundary,
1728 : // we just look for the top parent
1729 13740233 : if (elem->neighbor_ptr(side) == nullptr)
1730 623375 : searched_elem = elem->top_parent();
1731 : // Otherwise we loop over the ancestors and check if they have a different BC for us
1732 : else
1733 23273939 : while (searched_elem->parent() != nullptr)
1734 : {
1735 15365792 : const Elem * parent = searched_elem->parent();
1736 21170119 : if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
1737 8089578 : return;
1738 :
1739 9981250 : searched_elem = parent;
1740 : }
1741 : }
1742 :
1743 : #endif
1744 :
1745 : // Check each element in the range to see if its side matches the requested side.
1746 45989970 : for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
1747 21559523 : if (pr.second.first == side)
1748 7370134 : vec_to_fill.push_back(pr.second.second);
1749 : }
1750 :
1751 :
1752 :
1753 :
1754 189486 : unsigned int BoundaryInfo::n_boundary_ids (const Elem * const elem,
1755 : const unsigned short int side) const
1756 : {
1757 10704 : std::vector<boundary_id_type> ids;
1758 189486 : this->boundary_ids(elem, side, ids);
1759 195974 : return cast_int<unsigned int>(ids.size());
1760 : }
1761 :
1762 :
1763 217153339 : unsigned int BoundaryInfo::n_raw_boundary_ids (const Elem * const elem,
1764 : const unsigned short int side) const
1765 : {
1766 1363602 : std::vector<boundary_id_type> ids;
1767 217153339 : this->raw_boundary_ids(elem, side, ids);
1768 217854798 : return cast_int<unsigned int>(ids.size());
1769 : }
1770 :
1771 :
1772 :
1773 278811102 : void BoundaryInfo::raw_boundary_ids (const Elem * const elem,
1774 : const unsigned short int side,
1775 : std::vector<boundary_id_type> & vec_to_fill) const
1776 : {
1777 17421607 : libmesh_assert(elem);
1778 :
1779 : // Only query BCs for sides that exist.
1780 17421607 : libmesh_assert_less (side, elem->n_sides());
1781 :
1782 : // Clear out any previous contents
1783 17421607 : vec_to_fill.clear();
1784 :
1785 : // Only level-0 elements store BCs.
1786 279829761 : if (elem->parent() && !_children_on_boundary)
1787 9286648 : return;
1788 :
1789 : // Check each element in the range to see if its side matches the requested side.
1790 264769181 : for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
1791 94052184 : if (pr.second.first == side)
1792 29958729 : vec_to_fill.push_back(pr.second.second);
1793 : }
1794 :
1795 :
1796 :
1797 4353545 : void BoundaryInfo::copy_boundary_ids (const BoundaryInfo & old_boundary_info,
1798 : const Elem * const old_elem,
1799 : const Elem * const new_elem)
1800 : {
1801 133282 : libmesh_assert_equal_to (old_elem->n_sides(), new_elem->n_sides());
1802 133282 : libmesh_assert_equal_to (old_elem->n_edges(), new_elem->n_edges());
1803 :
1804 266564 : std::vector<boundary_id_type> bndry_ids;
1805 :
1806 22219972 : for (auto s : old_elem->side_index_range())
1807 : {
1808 17866427 : old_boundary_info.raw_boundary_ids (old_elem, s, bndry_ids);
1809 17866427 : this->add_side (new_elem, s, bndry_ids);
1810 : }
1811 :
1812 31022314 : for (auto e : old_elem->edge_index_range())
1813 : {
1814 26668769 : old_boundary_info.raw_edge_boundary_ids (old_elem, e, bndry_ids);
1815 26668769 : this->add_edge (new_elem, e, bndry_ids);
1816 : }
1817 :
1818 13060635 : for (unsigned short sf=0; sf != 2; sf++)
1819 : {
1820 8707090 : old_boundary_info.raw_shellface_boundary_ids (old_elem, sf, bndry_ids);
1821 8707090 : this->add_shellface (new_elem, sf, bndry_ids);
1822 : }
1823 4353545 : }
1824 :
1825 :
1826 :
1827 56199675 : void BoundaryInfo::remove (const Node * node)
1828 : {
1829 251662 : libmesh_assert(node);
1830 :
1831 : // Erase everything associated with node
1832 251662 : _boundary_node_id.erase (node);
1833 56199675 : }
1834 :
1835 :
1836 :
1837 107945 : void BoundaryInfo::remove_node (const Node * node,
1838 : const boundary_id_type id)
1839 : {
1840 3134 : libmesh_assert(node);
1841 :
1842 : // Erase (node, id) entry from map.
1843 107945 : erase_if(_boundary_node_id, node,
1844 181917 : [id](decltype(_boundary_node_id)::mapped_type & val)
1845 176583 : {return val == id;});
1846 107945 : }
1847 :
1848 :
1849 :
1850 43765468 : void BoundaryInfo::remove (const Elem * elem)
1851 : {
1852 264657 : libmesh_assert(elem);
1853 :
1854 : // Erase everything associated with elem
1855 264657 : _boundary_edge_id.erase (elem);
1856 264657 : _boundary_side_id.erase (elem);
1857 264657 : _boundary_shellface_id.erase (elem);
1858 43765468 : }
1859 :
1860 :
1861 :
1862 741036 : void BoundaryInfo::remove_edge (const Elem * elem,
1863 : const unsigned short int edge)
1864 : {
1865 38012 : libmesh_assert(elem);
1866 :
1867 : // Only touch BCs for edges that exist.
1868 38012 : libmesh_assert_less (edge, elem->n_edges());
1869 :
1870 : // Only level 0 elements are stored in BoundaryInfo.
1871 38012 : libmesh_assert_equal_to (elem->level(), 0);
1872 :
1873 : // Erase (elem, edge, *) entries from map.
1874 741036 : erase_if(_boundary_edge_id, elem,
1875 0 : [edge](decltype(_boundary_edge_id)::mapped_type & pr)
1876 0 : {return pr.first == edge;});
1877 741036 : }
1878 :
1879 :
1880 :
1881 110 : void BoundaryInfo::remove_edge (const Elem * elem,
1882 : const unsigned short int edge,
1883 : const boundary_id_type id)
1884 : {
1885 4 : libmesh_assert(elem);
1886 :
1887 : // Only touch BCs for edges that exist.
1888 4 : libmesh_assert_less (edge, elem->n_edges());
1889 :
1890 : // Only level 0 elements are stored in BoundaryInfo.
1891 4 : libmesh_assert_equal_to (elem->level(), 0);
1892 :
1893 : // Erase (elem, edge, id) entries from map.
1894 110 : erase_if(_boundary_edge_id, elem,
1895 330 : [edge, id](decltype(_boundary_edge_id)::mapped_type & pr)
1896 220 : {return pr.first == edge && pr.second == id;});
1897 110 : }
1898 :
1899 :
1900 0 : void BoundaryInfo::remove_shellface (const Elem * elem,
1901 : const unsigned short int shellface)
1902 : {
1903 0 : libmesh_assert(elem);
1904 :
1905 : // Only level 0 elements are stored in BoundaryInfo.
1906 0 : libmesh_assert_equal_to (elem->level(), 0);
1907 :
1908 : // Shells only have 2 faces
1909 0 : libmesh_assert_less(shellface, 2);
1910 :
1911 : // Erase (elem, shellface, *) entries from map.
1912 0 : erase_if(_boundary_shellface_id, elem,
1913 0 : [shellface](decltype(_boundary_shellface_id)::mapped_type & pr)
1914 0 : {return pr.first == shellface;});
1915 0 : }
1916 :
1917 :
1918 :
1919 0 : void BoundaryInfo::remove_shellface (const Elem * elem,
1920 : const unsigned short int shellface,
1921 : const boundary_id_type id)
1922 : {
1923 0 : libmesh_assert(elem);
1924 :
1925 : // Only level 0 elements are stored in BoundaryInfo.
1926 0 : libmesh_assert_equal_to (elem->level(), 0);
1927 :
1928 : // Shells only have 2 faces
1929 0 : libmesh_assert_less(shellface, 2);
1930 :
1931 : // Erase (elem, shellface, id) entries from map.
1932 0 : erase_if(_boundary_shellface_id, elem,
1933 0 : [shellface, id](decltype(_boundary_shellface_id)::mapped_type & pr)
1934 0 : {return pr.first == shellface && pr.second == id;});
1935 0 : }
1936 :
1937 376906 : void BoundaryInfo::remove_side (const Elem * elem,
1938 : const unsigned short int side)
1939 : {
1940 18968 : libmesh_assert(elem);
1941 :
1942 : // Only touch BCs for sides that exist.
1943 18968 : libmesh_assert_less (side, elem->n_sides());
1944 :
1945 : // Erase (elem, side, *) entries from map.
1946 376906 : erase_if(_boundary_side_id, elem,
1947 68806 : [side](decltype(_boundary_side_id)::mapped_type & pr)
1948 63727 : {return pr.first == side;});
1949 376906 : }
1950 :
1951 :
1952 :
1953 6441 : void BoundaryInfo::remove_side (const Elem * elem,
1954 : const unsigned short int side,
1955 : const boundary_id_type id)
1956 : {
1957 228 : libmesh_assert(elem);
1958 :
1959 : // Only touch BCs for sides that exist.
1960 228 : libmesh_assert_less (side, elem->n_sides());
1961 :
1962 : #ifdef LIBMESH_ENABLE_AMR
1963 : // Here we have to stop and check if somebody tries to remove an ancestor's boundary ID
1964 : // through a child
1965 6441 : if (elem->level())
1966 : {
1967 4 : std::vector<boundary_id_type> bd_ids;
1968 35 : this->boundary_ids(elem,side,bd_ids);
1969 35 : if(std::find(bd_ids.begin(), bd_ids.end(), id) != bd_ids.end())
1970 : {
1971 4 : std::vector<boundary_id_type> raw_bd_ids;
1972 35 : this->raw_boundary_ids(elem, side, raw_bd_ids);
1973 35 : if(std::find(raw_bd_ids.begin(), raw_bd_ids.end(), id) == raw_bd_ids.end())
1974 167 : libmesh_not_implemented_msg("We cannot delete boundary ID "
1975 : + std::to_string(id) +
1976 : " using a child because it is inherited from an ancestor.");
1977 : }
1978 : }
1979 : #endif
1980 :
1981 : // Erase (elem, side, id) entries from map.
1982 6406 : erase_if(_boundary_side_id, elem,
1983 25366 : [side, id](decltype(_boundary_side_id)::mapped_type & pr)
1984 18960 : {return pr.first == side && pr.second == id;});
1985 6406 : }
1986 :
1987 :
1988 :
1989 568 : void BoundaryInfo::remove_id (boundary_id_type id, const bool global)
1990 : {
1991 : // Pass global==false to the sub-methods here, so we can avoid
1992 : // unnecessary communications
1993 568 : this->remove_side_id(id, false);
1994 568 : this->remove_edge_id(id, false);
1995 568 : this->remove_shellface_id(id, false);
1996 568 : this->remove_node_id(id, false);
1997 :
1998 568 : if (global)
1999 0 : _global_boundary_ids.erase(id);
2000 568 : }
2001 :
2002 :
2003 :
2004 852 : void BoundaryInfo::remove_side_id (boundary_id_type id, const bool global)
2005 : {
2006 : // Erase id from ids containers
2007 24 : _side_boundary_ids.erase(id);
2008 :
2009 24 : if (!_edge_boundary_ids.count(id) &&
2010 48 : !_shellface_boundary_ids.count(id) &&
2011 24 : !_node_boundary_ids.count(id))
2012 2 : _boundary_ids.erase(id);
2013 :
2014 24 : _ss_id_to_name.erase(id);
2015 852 : if (global)
2016 : {
2017 0 : bool someone_has_it = _boundary_ids.count(id);
2018 0 : this->comm().max(someone_has_it);
2019 0 : if (!someone_has_it)
2020 0 : _global_boundary_ids.erase(id);
2021 : }
2022 :
2023 : // Erase (*, *, id) entries from map.
2024 852 : erase_if(_boundary_side_id,
2025 2240 : [id](decltype(_boundary_side_id)::mapped_type & pr)
2026 2112 : {return pr.second == id;});
2027 852 : }
2028 :
2029 :
2030 :
2031 568 : void BoundaryInfo::remove_edge_id (boundary_id_type id, const bool global)
2032 : {
2033 : // Erase id from ids containers
2034 16 : _edge_boundary_ids.erase(id);
2035 :
2036 16 : if (!_side_boundary_ids.count(id) &&
2037 32 : !_shellface_boundary_ids.count(id) &&
2038 16 : !_node_boundary_ids.count(id))
2039 2 : _boundary_ids.erase(id);
2040 :
2041 16 : _es_id_to_name.erase(id);
2042 568 : if (global)
2043 : {
2044 0 : bool someone_has_it = _boundary_ids.count(id);
2045 0 : this->comm().max(someone_has_it);
2046 0 : if (!someone_has_it)
2047 0 : _global_boundary_ids.erase(id);
2048 : }
2049 :
2050 : // Erase (*, *, id) entries from map.
2051 568 : erase_if(_boundary_edge_id,
2052 0 : [id](decltype(_boundary_edge_id)::mapped_type & pr)
2053 0 : {return pr.second == id;});
2054 568 : }
2055 :
2056 :
2057 :
2058 568 : void BoundaryInfo::remove_shellface_id (boundary_id_type id, const bool global)
2059 : {
2060 : // Erase id from ids containers
2061 16 : _shellface_boundary_ids.erase(id);
2062 :
2063 16 : if (!_side_boundary_ids.count(id) &&
2064 32 : !_edge_boundary_ids.count(id) &&
2065 16 : !_node_boundary_ids.count(id))
2066 2 : _boundary_ids.erase(id);
2067 :
2068 568 : if (global)
2069 : {
2070 0 : bool someone_has_it = _boundary_ids.count(id);
2071 0 : this->comm().max(someone_has_it);
2072 0 : if (!someone_has_it)
2073 0 : _global_boundary_ids.erase(id);
2074 : }
2075 :
2076 : // Erase (*, *, id) entries from map.
2077 568 : erase_if(_boundary_shellface_id,
2078 0 : [id](decltype(_boundary_shellface_id)::mapped_type & pr)
2079 0 : {return pr.second == id;});
2080 568 : }
2081 :
2082 :
2083 :
2084 852 : void BoundaryInfo::remove_node_id (boundary_id_type id, const bool global)
2085 : {
2086 : // Erase id from ids containers
2087 24 : _node_boundary_ids.erase(id);
2088 :
2089 24 : if (!_side_boundary_ids.count(id) &&
2090 40 : !_edge_boundary_ids.count(id) &&
2091 16 : !_shellface_boundary_ids.count(id))
2092 16 : _boundary_ids.erase(id);
2093 :
2094 24 : _ns_id_to_name.erase(id);
2095 852 : if (global)
2096 : {
2097 0 : bool someone_has_it = _boundary_ids.count(id);
2098 0 : this->comm().max(someone_has_it);
2099 0 : if (!someone_has_it)
2100 0 : _global_boundary_ids.erase(id);
2101 : }
2102 :
2103 : // Erase (*, id) entries from map.
2104 852 : erase_if(_boundary_node_id,
2105 3255 : [id](decltype(_boundary_node_id)::mapped_type & val)
2106 3069 : {return val == id;});
2107 852 : }
2108 :
2109 :
2110 :
2111 3692 : void BoundaryInfo::renumber_id (boundary_id_type old_id,
2112 : boundary_id_type new_id)
2113 : {
2114 3692 : if (old_id == new_id)
2115 : {
2116 : // If the IDs are the same, this is a no-op.
2117 48 : return;
2118 : }
2119 :
2120 56 : bool found_node = false;
2121 80132 : for (auto & p : _boundary_node_id)
2122 78144 : if (p.second == old_id)
2123 : {
2124 : // If we already have this id on this node, we don't want to
2125 : // create a duplicate in our multimap
2126 13164 : this->remove_node(p.first, new_id);
2127 13164 : p.second = new_id;
2128 456 : found_node = true;
2129 : }
2130 1988 : if (found_node)
2131 : {
2132 56 : _node_boundary_ids.erase(old_id);
2133 1500 : _node_boundary_ids.insert(new_id);
2134 : }
2135 :
2136 56 : bool found_edge = false;
2137 1988 : for (auto & p : _boundary_edge_id)
2138 0 : if (p.second.second == old_id)
2139 : {
2140 : // If we already have this id on this edge, we don't want to
2141 : // create a duplicate in our multimap
2142 0 : this->remove_edge(p.first, p.second.first, new_id);
2143 0 : p.second.second = new_id;
2144 0 : found_edge = true;
2145 : }
2146 1988 : if (found_edge)
2147 : {
2148 0 : _edge_boundary_ids.erase(old_id);
2149 0 : _edge_boundary_ids.insert(new_id);
2150 : }
2151 :
2152 56 : bool found_shellface = false;
2153 1988 : for (auto & p : _boundary_shellface_id)
2154 0 : if (p.second.second == old_id)
2155 : {
2156 : // If we already have this id on this shellface, we don't want
2157 : // to create a duplicate in our multimap
2158 0 : this->remove_shellface(p.first, p.second.first, new_id);
2159 0 : p.second.second = new_id;
2160 0 : found_shellface = true;
2161 : }
2162 1988 : if (found_shellface)
2163 : {
2164 0 : _shellface_boundary_ids.erase(old_id);
2165 0 : _shellface_boundary_ids.insert(new_id);
2166 : }
2167 :
2168 56 : bool found_side = false;
2169 37092 : for (auto & p : _boundary_side_id)
2170 35104 : if (p.second.second == old_id)
2171 : {
2172 : // If we already have this id on this side, we don't want to
2173 : // create a duplicate in our multimap
2174 5944 : this->remove_side(p.first, p.second.first, new_id);
2175 5944 : p.second.second = new_id;
2176 208 : found_side = true;
2177 : }
2178 1988 : if (found_side)
2179 : {
2180 56 : _side_boundary_ids.erase(old_id);
2181 1500 : _side_boundary_ids.insert(new_id);
2182 : }
2183 :
2184 1988 : if (found_node || found_edge || found_shellface || found_side)
2185 : {
2186 56 : _boundary_ids.erase(old_id);
2187 1500 : _boundary_ids.insert(new_id);
2188 56 : _global_boundary_ids.erase(old_id);
2189 1500 : _global_boundary_ids.insert(new_id);
2190 : }
2191 :
2192 1988 : renumber_name(_ss_id_to_name, old_id, new_id);
2193 1988 : renumber_name(_ns_id_to_name, old_id, new_id);
2194 1988 : renumber_name(_es_id_to_name, old_id, new_id);
2195 :
2196 1988 : this->libmesh_assert_valid_multimaps();
2197 : }
2198 :
2199 :
2200 :
2201 142 : void BoundaryInfo::renumber_side_id (boundary_id_type old_id,
2202 : boundary_id_type new_id)
2203 : {
2204 : // If the IDs are the same, this is a no-op.
2205 142 : if (old_id == new_id)
2206 0 : return;
2207 :
2208 4 : bool found_side = false;
2209 2782 : for (auto & p : _boundary_side_id)
2210 2640 : if (p.second.second == old_id)
2211 : {
2212 : // If we already have this id on this side, we don't want to
2213 : // create a duplicate in our multimap
2214 440 : this->remove_side(p.first, p.second.first, new_id);
2215 440 : p.second.second = new_id;
2216 16 : found_side = true;
2217 : }
2218 4 : _side_boundary_ids.erase(old_id);
2219 :
2220 142 : if (found_side)
2221 : {
2222 106 : _side_boundary_ids.insert(new_id);
2223 :
2224 4 : if (!_shellface_boundary_ids.count(old_id) &&
2225 8 : !_edge_boundary_ids.count(old_id) &&
2226 4 : !_node_boundary_ids.count(old_id))
2227 : {
2228 0 : _boundary_ids.erase(old_id);
2229 0 : _global_boundary_ids.erase(old_id);
2230 : }
2231 106 : _boundary_ids.insert(new_id);
2232 106 : _global_boundary_ids.insert(new_id);
2233 : }
2234 :
2235 142 : renumber_name(_ss_id_to_name, old_id, new_id);
2236 :
2237 142 : this->libmesh_assert_valid_multimaps();
2238 : }
2239 :
2240 :
2241 :
2242 110 : void BoundaryInfo::renumber_edge_id (boundary_id_type old_id,
2243 : boundary_id_type new_id)
2244 : {
2245 : // If the IDs are the same, this is a no-op.
2246 110 : if (old_id == new_id)
2247 0 : return;
2248 :
2249 4 : bool found_edge = false;
2250 330 : for (auto & p : _boundary_edge_id)
2251 220 : if (p.second.second == old_id)
2252 : {
2253 : // If we already have this id on this edge, we don't want to
2254 : // create a duplicate in our multimap
2255 110 : this->remove_edge(p.first, p.second.first, new_id);
2256 110 : p.second.second = new_id;
2257 4 : found_edge = true;
2258 : }
2259 4 : _edge_boundary_ids.erase(old_id);
2260 :
2261 110 : if (found_edge)
2262 : {
2263 106 : _edge_boundary_ids.insert(new_id);
2264 :
2265 4 : if (!_shellface_boundary_ids.count(old_id) &&
2266 4 : !_side_boundary_ids.count(old_id) &&
2267 0 : !_node_boundary_ids.count(old_id))
2268 : {
2269 0 : _boundary_ids.erase(old_id);
2270 0 : _global_boundary_ids.erase(old_id);
2271 : }
2272 106 : _boundary_ids.insert(new_id);
2273 106 : _global_boundary_ids.insert(new_id);
2274 : }
2275 :
2276 110 : renumber_name(_es_id_to_name, old_id, new_id);
2277 :
2278 110 : this->libmesh_assert_valid_multimaps();
2279 : }
2280 :
2281 :
2282 :
2283 0 : void BoundaryInfo::renumber_shellface_id (boundary_id_type old_id,
2284 : boundary_id_type new_id)
2285 : {
2286 : // If the IDs are the same, this is a no-op.
2287 0 : if (old_id == new_id)
2288 0 : return;
2289 :
2290 0 : bool found_shellface = false;
2291 0 : for (auto & p : _boundary_shellface_id)
2292 0 : if (p.second.second == old_id)
2293 : {
2294 : // If we already have this id on this shellface, we don't want
2295 : // to create a duplicate in our multimap
2296 0 : this->remove_shellface(p.first, p.second.first, new_id);
2297 0 : p.second.second = new_id;
2298 0 : found_shellface = true;
2299 : }
2300 0 : _shellface_boundary_ids.erase(old_id);
2301 :
2302 0 : if (found_shellface)
2303 : {
2304 0 : _shellface_boundary_ids.insert(new_id);
2305 :
2306 0 : if (!_edge_boundary_ids.count(old_id) &&
2307 0 : !_side_boundary_ids.count(old_id) &&
2308 0 : !_node_boundary_ids.count(old_id))
2309 : {
2310 0 : _boundary_ids.erase(old_id);
2311 0 : _global_boundary_ids.erase(old_id);
2312 : }
2313 0 : _boundary_ids.insert(new_id);
2314 0 : _global_boundary_ids.insert(new_id);
2315 : }
2316 :
2317 0 : this->libmesh_assert_valid_multimaps();
2318 : }
2319 :
2320 :
2321 :
2322 142 : void BoundaryInfo::renumber_node_id (boundary_id_type old_id,
2323 : boundary_id_type new_id)
2324 : {
2325 : // If the IDs are the same, this is a no-op.
2326 142 : if (old_id == new_id)
2327 0 : return;
2328 :
2329 4 : bool found_node = false;
2330 6082 : for (auto & p : _boundary_node_id)
2331 5940 : if (p.second == old_id)
2332 : {
2333 : // If we already have this id on this node, we don't want to
2334 : // create a duplicate in our multimap
2335 990 : this->remove_node(p.first, new_id);
2336 990 : p.second = new_id;
2337 36 : found_node = true;
2338 : }
2339 4 : _node_boundary_ids.erase(old_id);
2340 :
2341 142 : if (found_node)
2342 : {
2343 106 : _node_boundary_ids.insert(new_id);
2344 :
2345 4 : if (!_shellface_boundary_ids.count(old_id) &&
2346 8 : !_side_boundary_ids.count(old_id) &&
2347 4 : !_edge_boundary_ids.count(old_id))
2348 : {
2349 4 : _boundary_ids.erase(old_id);
2350 4 : _global_boundary_ids.erase(old_id);
2351 : }
2352 106 : _boundary_ids.insert(new_id);
2353 106 : _global_boundary_ids.insert(new_id);
2354 : }
2355 :
2356 142 : renumber_name(_ns_id_to_name, old_id, new_id);
2357 :
2358 142 : this->libmesh_assert_valid_multimaps();
2359 : }
2360 :
2361 :
2362 :
2363 234 : unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
2364 : const boundary_id_type boundary_id_in) const
2365 : {
2366 234 : const Elem * searched_elem = elem;
2367 :
2368 : // If we don't have a time-dependent domain, we can just go ahead and use the top parent
2369 : // (since only those contain boundary conditions). Otherwise, we keep the element
2370 234 : if (elem->level() != 0 && !_children_on_boundary)
2371 0 : searched_elem = elem->top_parent();
2372 :
2373 : // elem may have zero or multiple occurrences
2374 489 : for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
2375 : {
2376 : // if this is true we found the requested boundary_id
2377 : // of the element and want to return the side
2378 378 : if (pr.second.second == boundary_id_in)
2379 : {
2380 123 : unsigned int side = pr.second.first;
2381 :
2382 : // Here we branch out. If we don't allow time-dependent boundary domains,
2383 : // we need to check if our parents are consistent.
2384 123 : if (!_children_on_boundary)
2385 : {
2386 : #ifdef LIBMESH_ENABLE_AMR
2387 : // If we're on this external boundary then we share this
2388 : // external boundary id
2389 91 : if (elem->neighbor_ptr(side) == nullptr)
2390 9 : return side;
2391 :
2392 : // Internal boundary case
2393 7 : const Elem * p = elem;
2394 :
2395 : // If we're on an internal boundary then we need to be sure
2396 : // it's the same internal boundary as our top_parent
2397 168 : while (p != nullptr)
2398 : {
2399 14 : const Elem * parent = p->parent();
2400 154 : if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side))
2401 0 : break;
2402 7 : p = parent;
2403 : }
2404 :
2405 : // We're on that side of our top_parent; return it
2406 14 : if (!p)
2407 84 : return side;
2408 : #else
2409 : // do not forget to return the internal boundary when AMR is disabled
2410 : return side;
2411 : #endif
2412 : }
2413 : // Otherwise we need to check if the child's ancestors have something on
2414 : // the side of the child
2415 : else
2416 39 : return side;
2417 : }
2418 : }
2419 :
2420 : #ifdef LIBMESH_ENABLE_AMR
2421 : // We might have instances (especially with moving boundary domains) when we
2422 : // query the paren't boundary ID on a child. We only do this till we find the
2423 : // the first side, for multiple sides see above.
2424 111 : if (_children_on_boundary && elem->level() != 0)
2425 : {
2426 78 : for (auto side : make_range(elem->n_sides()))
2427 : {
2428 4 : const Elem * p = elem;
2429 164 : while (p->parent() != nullptr)
2430 : {
2431 117 : const Elem * parent = p->parent();
2432 :
2433 : // First we make sure the parent shares this side
2434 117 : if (parent->is_child_on_side(parent->which_child_am_i(p), side))
2435 : {
2436 : // parent may have multiple boundary ids
2437 351 : for (const auto & pr : as_range(_boundary_side_id.equal_range(parent)))
2438 : // if this is true we found the requested boundary_id
2439 : // of the element and want to return the side
2440 273 : if (pr.second.first == side && pr.second.second == boundary_id_in)
2441 2 : return side;
2442 :
2443 8 : p = parent;
2444 : }
2445 : // If the parent is not on the same side, other ancestors won't be on the same side either
2446 : else
2447 0 : break;
2448 : }
2449 : }
2450 : }
2451 : #endif
2452 :
2453 : // if we get here, we found elem in the data structure but not
2454 : // the requested boundary id, so return the default value
2455 6 : return libMesh::invalid_uint;
2456 : }
2457 :
2458 :
2459 : std::vector<unsigned int>
2460 65698 : BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
2461 : const boundary_id_type boundary_id_in) const
2462 : {
2463 25668 : std::vector<unsigned int> returnval;
2464 :
2465 65698 : const Elem * searched_elem = elem;
2466 65698 : if (elem->level() != 0 && !_children_on_boundary)
2467 49065 : searched_elem = elem->top_parent();
2468 :
2469 : // elem may have zero or multiple occurrences
2470 184871 : for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
2471 : {
2472 : // if this is true we found the requested boundary_id
2473 : // of the element and want to return the side
2474 119173 : if (pr.second.second == boundary_id_in)
2475 : {
2476 63480 : unsigned int side = pr.second.first;
2477 :
2478 : // Here we branch out. If we don't allow time-dependent boundary domains,
2479 : // we need to check if our parents are consistent.
2480 63480 : if (!_children_on_boundary)
2481 : {
2482 : // If we're on this external boundary then we share this
2483 : // external boundary id
2484 82498 : if (elem->neighbor_ptr(side) == nullptr)
2485 : {
2486 63176 : returnval.push_back(side);
2487 63176 : continue;
2488 : }
2489 :
2490 : // If we're on an internal boundary then we need to be sure
2491 : // it's the same internal boundary as our top_parent
2492 0 : const Elem * p = elem;
2493 :
2494 : #ifdef LIBMESH_ENABLE_AMR
2495 :
2496 608 : while (p != nullptr)
2497 : {
2498 0 : const Elem * parent = p->parent();
2499 608 : if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side))
2500 0 : break;
2501 0 : p = parent;
2502 : }
2503 : #endif
2504 : // We're on that side of our top_parent; return it
2505 304 : if (!p)
2506 304 : returnval.push_back(side);
2507 : }
2508 : // Otherwise we trust what we got and return the side
2509 : else
2510 0 : returnval.push_back(side);
2511 : }
2512 : }
2513 :
2514 : #ifdef LIBMESH_ENABLE_AMR
2515 : // We might have instances (especially with moving boundary domains) when we
2516 : // query the parent boundary ID on a child.
2517 65698 : if (_children_on_boundary && elem->level() != 0)
2518 : {
2519 255 : for (auto side : make_range(elem->n_sides()))
2520 : {
2521 8 : const Elem * p = elem;
2522 318 : while (p->parent() != nullptr)
2523 : {
2524 306 : const Elem * parent = p->parent();
2525 : // First we make sure the parent shares this side
2526 306 : if (parent->is_child_on_side(parent->which_child_am_i(p), side))
2527 : {
2528 : // parent may have multiple boundary ids
2529 306 : for (const auto & pr : as_range(_boundary_side_id.equal_range(parent)))
2530 : {
2531 : // if this is true we found the requested boundary_id
2532 : // of the element and want to add the side to the vector. We
2533 : // also need to check if the side is already in the vector. This might
2534 : // happen if the child inherits the boundary from the parent.
2535 212 : if (pr.second.first == side && pr.second.second == boundary_id_in &&
2536 208 : std::find(returnval.begin(), returnval.end(), side) == returnval.end())
2537 102 : returnval.push_back(side);
2538 : }
2539 : }
2540 : // If the parent is not on the same side, other ancestors won't be on the same side either
2541 : else
2542 8 : break;
2543 8 : p = parent;
2544 : }
2545 : }
2546 : }
2547 : #endif
2548 :
2549 91366 : return returnval;
2550 : }
2551 :
2552 : void
2553 7310 : BoundaryInfo::build_node_boundary_ids(std::vector<boundary_id_type> & b_ids) const
2554 : {
2555 644 : b_ids.clear();
2556 :
2557 1652032 : for (const auto & pr : _boundary_node_id)
2558 : {
2559 1644722 : boundary_id_type id = pr.second;
2560 :
2561 1644722 : if (std::find(b_ids.begin(),b_ids.end(),id) == b_ids.end())
2562 26818 : b_ids.push_back(id);
2563 : }
2564 7310 : }
2565 :
2566 : void
2567 7310 : BoundaryInfo::build_side_boundary_ids(std::vector<boundary_id_type> & b_ids) const
2568 : {
2569 644 : b_ids.clear();
2570 :
2571 983514 : for (const auto & pr : _boundary_side_id)
2572 : {
2573 976204 : boundary_id_type id = pr.second.second;
2574 :
2575 976204 : if (std::find(b_ids.begin(),b_ids.end(),id) == b_ids.end())
2576 27768 : b_ids.push_back(id);
2577 : }
2578 7310 : }
2579 :
2580 : void
2581 7310 : BoundaryInfo::build_shellface_boundary_ids(std::vector<boundary_id_type> & b_ids) const
2582 : {
2583 644 : b_ids.clear();
2584 :
2585 87182 : for (const auto & pr :_boundary_shellface_id)
2586 : {
2587 79872 : boundary_id_type id = pr.second.second;
2588 :
2589 79872 : if (std::find(b_ids.begin(),b_ids.end(),id) == b_ids.end())
2590 192 : b_ids.push_back(id);
2591 : }
2592 7310 : }
2593 :
2594 : #ifdef LIBMESH_ENABLE_AMR
2595 : void
2596 19234506 : BoundaryInfo::transfer_boundary_ids_from_children(const Elem * const parent)
2597 : {
2598 : // this is only needed when we allow boundary to be associated with children elements
2599 : // also, we only transfer the parent's boundary ids when we are actually coarsen the child element
2600 19234530 : if (!_children_on_boundary ||
2601 386 : !(!parent->active() && parent->refinement_flag() == Elem::COARSEN_INACTIVE))
2602 921948 : return;
2603 :
2604 : // We assume that edges can be divided ito two pieces, while triangles and
2605 : // quads can be divided into four smaller areas. This is double because we'll need
2606 : // to convert the ratio of the children with given boundary id to a double.
2607 82 : const double number_of_sides_on_children = std::pow(2, parent->dim()-1);
2608 :
2609 : // In this case the input argument elem is the parent element. We need to check all of its sides
2610 : // to grab any potential boundary ids.
2611 410 : for (unsigned int side_i = 0; side_i < parent->n_sides(); ++side_i)
2612 : {
2613 : // An temporary storage to count how many times the children's boundaries occur. the general
2614 : // consensus is that if the boundary occurs more than once we propagate upon coarsening. Otherwise,
2615 : // it will get deleted.
2616 32 : std::map<unsigned short int, unsigned short int> boundary_counts;
2617 :
2618 1640 : for (const auto & child_i : make_range(parent->n_children()))
2619 : {
2620 : // We only need to check the children which share the side
2621 1312 : if (parent->is_child_on_side(child_i, side_i))
2622 : {
2623 : // Fetching the boundary tags on the child's side
2624 902 : for (const auto & pr : as_range(_boundary_side_id.equal_range(parent->child_ptr(child_i))))
2625 : {
2626 : // Making sure we are on the same boundary
2627 246 : if (pr.second.first == side_i)
2628 123 : ++boundary_counts[pr.second.second];
2629 : }
2630 : }
2631 : }
2632 :
2633 : // This is where the decision is made. If 50% of the children have the tags,
2634 : // we propagate them upwards upon coarsening. Otherwise, they are deleted.
2635 410 : for (const auto & boundary : boundary_counts)
2636 82 : if (boundary.second / number_of_sides_on_children > 0.5)
2637 41 : this->add_side(parent, side_i, boundary.first);
2638 : }
2639 :
2640 410 : for (const auto & child_i : make_range(parent->n_children()))
2641 328 : this->remove(parent->child_ptr(child_i));
2642 : }
2643 : #endif
2644 :
2645 17923 : std::size_t BoundaryInfo::n_boundary_conds () const
2646 : {
2647 : // in serial we know the number of bcs from the
2648 : // size of the container
2649 17923 : if (_mesh->is_serial())
2650 7088 : return _boundary_side_id.size();
2651 :
2652 : // in parallel we need to sum the number of local bcs
2653 42 : parallel_object_only();
2654 :
2655 10835 : std::size_t nbcs=0;
2656 :
2657 72506 : for (const auto & pr : _boundary_side_id)
2658 61711 : if (pr.first->processor_id() == this->processor_id())
2659 16504 : nbcs++;
2660 :
2661 10835 : this->comm().sum (nbcs);
2662 :
2663 10835 : return nbcs;
2664 : }
2665 :
2666 38082 : std::size_t BoundaryInfo::n_edge_conds () const
2667 : {
2668 : // in serial we know the number of nodesets from the
2669 : // size of the container
2670 38082 : if (_mesh->is_serial())
2671 8232 : return _boundary_edge_id.size();
2672 :
2673 : // in parallel we need to sum the number of local nodesets
2674 26 : parallel_object_only();
2675 :
2676 29850 : std::size_t n_edge_bcs=0;
2677 :
2678 31670 : for (const auto & pr : _boundary_edge_id)
2679 1820 : if (pr.first->processor_id() == this->processor_id())
2680 792 : n_edge_bcs++;
2681 :
2682 29850 : this->comm().sum (n_edge_bcs);
2683 :
2684 29850 : return n_edge_bcs;
2685 : }
2686 :
2687 :
2688 5542 : std::size_t BoundaryInfo::n_shellface_conds () const
2689 : {
2690 : // in serial we know the number of nodesets from the
2691 : // size of the container
2692 5542 : if (_mesh->is_serial())
2693 856 : return _boundary_shellface_id.size();
2694 :
2695 : // in parallel we need to sum the number of local nodesets
2696 4 : parallel_object_only();
2697 :
2698 4686 : std::size_t n_shellface_bcs=0;
2699 :
2700 46320 : for (const auto & pr : _boundary_shellface_id)
2701 41634 : if (pr.first->processor_id() == this->processor_id())
2702 26640 : n_shellface_bcs++;
2703 :
2704 4686 : this->comm().sum (n_shellface_bcs);
2705 :
2706 4686 : return n_shellface_bcs;
2707 : }
2708 :
2709 :
2710 9021 : std::size_t BoundaryInfo::n_nodeset_conds () const
2711 : {
2712 : // in serial we know the number of nodesets from the
2713 : // size of the container
2714 9021 : if (_mesh->is_serial())
2715 4399 : return _boundary_node_id.size();
2716 :
2717 : // in parallel we need to sum the number of local nodesets
2718 4 : parallel_object_only();
2719 :
2720 4622 : std::size_t n_nodesets=0;
2721 :
2722 174325 : for (const auto & pr : _boundary_node_id)
2723 169759 : if (pr.first->processor_id() == this->processor_id())
2724 31584 : n_nodesets++;
2725 :
2726 4622 : this->comm().sum (n_nodesets);
2727 :
2728 4622 : return n_nodesets;
2729 : }
2730 :
2731 :
2732 :
2733 : std::vector<BoundaryInfo::NodeBCTuple>
2734 27965 : BoundaryInfo::build_node_list(NodeBCTupleSortBy sort_by) const
2735 : {
2736 1009 : std::vector<NodeBCTuple> bc_tuples;
2737 27965 : bc_tuples.reserve(_boundary_node_id.size());
2738 :
2739 2140844 : for (const auto & [node, bid] : _boundary_node_id)
2740 2112879 : bc_tuples.emplace_back(node->id(), bid);
2741 :
2742 : // This list is currently in memory address (arbitrary) order, so
2743 : // sort, using the specified ordering, to make it consistent on all procs.
2744 27965 : if (sort_by == NodeBCTupleSortBy::NODE_ID)
2745 27965 : std::sort(bc_tuples.begin(), bc_tuples.end());
2746 0 : else if (sort_by == NodeBCTupleSortBy::BOUNDARY_ID)
2747 0 : std::sort(bc_tuples.begin(), bc_tuples.end(),
2748 0 : [](const NodeBCTuple & left, const NodeBCTuple & right)
2749 0 : {return std::get<1>(left) < std::get<1>(right);});
2750 :
2751 27965 : return bc_tuples;
2752 : }
2753 :
2754 :
2755 : void
2756 142 : BoundaryInfo::build_node_list_from_side_list(const std::set<boundary_id_type> & sideset_list)
2757 : {
2758 : // If we're on a distributed mesh, even the owner of a node is not
2759 : // guaranteed to be able to properly assign its new boundary id(s)!
2760 : // Nodal neighbors are not always ghosted, and a nodal neighbor
2761 : // might have a boundary side.
2762 142 : const bool mesh_is_serial = _mesh->is_serial();
2763 :
2764 : typedef std::set<std::pair<dof_id_type, boundary_id_type>> set_type;
2765 : typedef std::vector<std::pair<dof_id_type, boundary_id_type>> vec_type;
2766 :
2767 8 : const processor_id_type my_proc_id = this->processor_id();
2768 4 : std::unordered_map<processor_id_type, set_type> nodes_to_push;
2769 4 : std::unordered_map<processor_id_type, vec_type> node_vecs_to_push;
2770 :
2771 : // For avoiding extraneous element side construction
2772 4 : ElemSideBuilder side_builder;
2773 : // Pull objects out of the loop to reduce heap operations
2774 : const Elem * side;
2775 :
2776 : // Loop over the side list
2777 937 : for (const auto & [elem, id_pair] : _boundary_side_id)
2778 : {
2779 : // Don't add remote sides
2780 795 : if (elem->is_remote())
2781 140 : continue;
2782 :
2783 795 : auto [sidenum, bcid] = id_pair;
2784 :
2785 795 : if (!sideset_list.empty() && !sideset_list.count(bcid))
2786 132 : continue;
2787 :
2788 : // Need to loop over the sides of any possible children
2789 96 : std::vector<const Elem *> family;
2790 : #ifdef LIBMESH_ENABLE_AMR
2791 655 : elem->active_family_tree_by_side (family, sidenum);
2792 : #else
2793 : family.push_back(elem);
2794 : #endif
2795 :
2796 1310 : for (const auto & cur_elem : family)
2797 : {
2798 655 : side = &side_builder(*cur_elem, sidenum);
2799 :
2800 : // Add each node node on the side with the side's boundary id
2801 1965 : for (auto i : side->node_index_range())
2802 : {
2803 1406 : this->add_node(side->node_ptr(i), bcid);
2804 1310 : if (!mesh_is_serial)
2805 : {
2806 : const processor_id_type proc_id =
2807 974 : side->node_ptr(i)->processor_id();
2808 974 : if (proc_id != my_proc_id)
2809 590 : nodes_to_push[proc_id].emplace(side->node_id(i), bcid);
2810 : }
2811 : }
2812 : }
2813 : }
2814 :
2815 : // If we're on a serial mesh then we're done.
2816 142 : if (mesh_is_serial)
2817 4 : return;
2818 :
2819 : // Otherwise we need to push ghost node bcids to their owners, then
2820 : // pull ghost node bcids from their owners.
2821 :
2822 372 : for (auto & [proc_id, s] : nodes_to_push)
2823 : {
2824 244 : node_vecs_to_push[proc_id].assign(s.begin(), s.end());
2825 0 : s.clear();
2826 : }
2827 :
2828 : auto nodes_action_functor =
2829 244 : [this]
2830 : (processor_id_type,
2831 423 : const vec_type & received_nodes)
2832 : {
2833 667 : for (const auto & [dof_id, bndry_id] : received_nodes)
2834 423 : this->add_node(_mesh->node_ptr(dof_id), bndry_id);
2835 128 : };
2836 :
2837 : Parallel::push_parallel_vector_data
2838 128 : (this->comm(), node_vecs_to_push, nodes_action_functor);
2839 :
2840 : // At this point we should know all the BCs for our own nodes; now
2841 : // we need BCs for ghost nodes.
2842 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2843 0 : node_ids_requested;
2844 :
2845 : // Determine what nodes we need to request
2846 2902 : for (const auto & node : _mesh->node_ptr_range())
2847 : {
2848 1323 : const processor_id_type pid = node->processor_id();
2849 1323 : if (pid != my_proc_id)
2850 963 : node_ids_requested[pid].push_back(node->id());
2851 128 : }
2852 :
2853 : typedef std::vector<boundary_id_type> datum_type;
2854 :
2855 : auto node_bcid_gather_functor =
2856 378 : [this]
2857 : (processor_id_type,
2858 : const std::vector<dof_id_type> & ids,
2859 963 : std::vector<datum_type> & data)
2860 : {
2861 0 : const std::size_t query_size = ids.size();
2862 378 : data.resize(query_size);
2863 :
2864 1341 : for (std::size_t i=0; i != query_size; ++i)
2865 963 : this->boundary_ids(_mesh->node_ptr(ids[i]), data[i]);
2866 506 : };
2867 :
2868 : auto node_bcid_action_functor =
2869 378 : [this]
2870 : (processor_id_type,
2871 : const std::vector<dof_id_type> & ids,
2872 963 : const std::vector<datum_type> & data)
2873 : {
2874 1341 : for (auto i : index_range(ids))
2875 963 : this->add_node(_mesh->node_ptr(ids[i]), data[i]);
2876 128 : };
2877 :
2878 0 : datum_type * datum_type_ex = nullptr;
2879 : Parallel::pull_parallel_vector_data
2880 128 : (this->comm(), node_ids_requested, node_bcid_gather_functor,
2881 : node_bcid_action_functor, datum_type_ex);
2882 : }
2883 :
2884 0 : void BoundaryInfo::parallel_sync_side_ids()
2885 : {
2886 : // we need BCs for ghost elements.
2887 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2888 0 : elem_ids_requested;
2889 :
2890 : // Determine what elements we need to request
2891 0 : for (const auto & elem : _mesh->element_ptr_range())
2892 : {
2893 0 : const processor_id_type pid = elem->processor_id();
2894 0 : if (pid != this->processor_id())
2895 0 : elem_ids_requested[pid].push_back(elem->id());
2896 0 : }
2897 :
2898 : typedef std::vector<std::pair<unsigned short int, boundary_id_type>> datum_type;
2899 :
2900 : // gather the element ID, side, and boundary_id_type for the ghost elements
2901 : auto elem_id_gather_functor =
2902 0 : [this]
2903 : (processor_id_type,
2904 : const std::vector<dof_id_type> & ids,
2905 0 : std::vector<datum_type> & data)
2906 : {
2907 0 : data.resize(ids.size());
2908 0 : for (auto i : index_range(ids))
2909 : {
2910 0 : Elem * elem = _mesh->elem_ptr(ids[i]);
2911 0 : for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
2912 0 : data[i].push_back(std::make_pair(pr.second.first, pr.second.second));
2913 : }
2914 0 : };
2915 : // update the _boundary_side_id on this processor
2916 : auto elem_id_action_functor =
2917 0 : [this]
2918 : (processor_id_type,
2919 : const std::vector<dof_id_type> & ids,
2920 0 : std::vector<datum_type> & data)
2921 : {
2922 0 : for (auto i : index_range(ids))
2923 : {
2924 0 : Elem * elem = _mesh->elem_ptr(ids[i]);
2925 : //clear boundary sides for this element
2926 0 : _boundary_side_id.erase(elem);
2927 : // update boundary sides for it
2928 0 : for (const auto & [side_id, bndry_id] : data[i])
2929 0 : _boundary_side_id.insert(std::make_pair(elem, std::make_pair(side_id, bndry_id)));
2930 : }
2931 0 : };
2932 :
2933 :
2934 0 : datum_type * datum_type_ex = nullptr;
2935 : Parallel::pull_parallel_vector_data
2936 0 : (this->comm(), elem_ids_requested, elem_id_gather_functor,
2937 : elem_id_action_functor, datum_type_ex);
2938 0 : }
2939 :
2940 0 : void BoundaryInfo::parallel_sync_node_ids()
2941 : {
2942 : // we need BCs for ghost nodes.
2943 : std::unordered_map<processor_id_type, std::vector<dof_id_type>>
2944 0 : node_ids_requested;
2945 :
2946 : // Determine what nodes we need to request
2947 0 : for (const auto & node : _mesh->node_ptr_range())
2948 : {
2949 0 : const processor_id_type pid = node->processor_id();
2950 0 : if (pid != this->processor_id())
2951 0 : node_ids_requested[pid].push_back(node->id());
2952 0 : }
2953 :
2954 : typedef std::vector<boundary_id_type> datum_type;
2955 :
2956 : // gather the node ID and boundary_id_type for the ghost nodes
2957 : auto node_id_gather_functor =
2958 0 : [this]
2959 : (processor_id_type,
2960 : const std::vector<dof_id_type> & ids,
2961 0 : std::vector<datum_type> & data)
2962 : {
2963 0 : data.resize(ids.size());
2964 0 : for (auto i : index_range(ids))
2965 : {
2966 0 : Node * node = _mesh->node_ptr(ids[i]);
2967 0 : for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
2968 0 : data[i].push_back(pr.second);
2969 : }
2970 0 : };
2971 :
2972 : // update the _boundary_node_id on this processor
2973 : auto node_id_action_functor =
2974 0 : [this]
2975 : (processor_id_type,
2976 : const std::vector<dof_id_type> & ids,
2977 0 : std::vector<datum_type> & data)
2978 : {
2979 0 : for (auto i : index_range(ids))
2980 : {
2981 0 : Node * node = _mesh->node_ptr(ids[i]);
2982 : //clear boundary node
2983 0 : _boundary_node_id.erase(node);
2984 : // update boundary node
2985 0 : for (const auto & pr : data[i])
2986 0 : _boundary_node_id.insert(std::make_pair(node, pr));
2987 : }
2988 0 : };
2989 :
2990 :
2991 0 : datum_type * datum_type_ex = nullptr;
2992 : Parallel::pull_parallel_vector_data
2993 0 : (this->comm(), node_ids_requested, node_id_gather_functor,
2994 : node_id_action_functor, datum_type_ex);
2995 0 : }
2996 :
2997 284 : void BoundaryInfo::build_side_list_from_node_list(const std::set<boundary_id_type> & nodeset_list,
2998 : bool skip_interior_sides)
2999 : {
3000 : // Skipping interior sides on a distributed mesh means resolving the
3001 : // subdomain id of any neighbor we only see as a RemoteElem, which
3002 : // takes a collective communication below. When we're not doing
3003 : // that, this function stays purely local, the way it has always
3004 : // been: existing callers are not required to call us in parallel
3005 : // (AbaqusIO, our only in-library caller, only reads on processor 0),
3006 : // so we must not add a collective to that path. When we are doing
3007 : // it, every processor has to call us, with the same arguments.
3008 284 : const bool need_parallel_sync = skip_interior_sides && !_mesh->is_serial();
3009 :
3010 : // Check for early return. When we're going to communicate below
3011 : // this has to be a collective decision: on a DistributedMesh it's
3012 : // entirely possible for _boundary_node_id to be empty on some
3013 : // processors (e.g. ones that don't own or ghost any node in the
3014 : // nodeset) while it's non-empty on others. Returning early based
3015 : // only on our own local state would leave us out of step with the
3016 : // processors that don't return early, causing a deadlock. So in
3017 : // that case we only take the early return if *every* processor
3018 : // agrees there's nothing to do.
3019 284 : bool no_boundary_nodes = _boundary_node_id.empty();
3020 284 : if (need_parallel_sync)
3021 135 : this->comm().min(no_boundary_nodes);
3022 284 : if (no_boundary_nodes)
3023 : {
3024 0 : libMesh::out << "No boundary node IDs have been added: cannot build side list!" << std::endl;
3025 149 : return;
3026 : }
3027 :
3028 : // For avoiding extraneous element side construction
3029 8 : ElemSideBuilder side_builder;
3030 : // Pull objects out of the loop to reduce heap operations
3031 : const Elem * side_elem;
3032 8 : std::vector<boundary_id_type> matched_ids;
3033 :
3034 : // Add elem's side to each of the given sidesets, making sure each
3035 : // sideset inherits its nodeset's name, if there is one.
3036 : auto add_matched_side =
3037 1186 : [this]
3038 : (const Elem * elem,
3039 : unsigned short int side,
3040 3172 : const std::vector<boundary_id_type> & bids)
3041 : {
3042 2812 : for (const auto & bid : bids)
3043 : {
3044 1406 : add_side(elem, side, bid);
3045 :
3046 : // Let the sideset inherit any non-empty name from the nodeset
3047 1406 : std::string & nset_name = nodeset_name(bid);
3048 :
3049 1406 : if (nset_name != "")
3050 140 : sideset_name(bid) = nset_name;
3051 : }
3052 334 : };
3053 :
3054 : // Sides whose nodes are all in the nodeset, but whose neighbor is
3055 : // an unresolved RemoteElem: we can't tell locally whether they're
3056 : // truly interior to a subdomain, since RemoteElem doesn't carry a
3057 : // meaningful subdomain_id(). We resolve them below by asking
3058 : // whichever processor owns each element: since that processor owns
3059 : // the element, its neighbor connectivity for it is guaranteed to be
3060 : // fully resolved, unlike ours. For each owning processor we keep
3061 : // the (elem id, side) queries and, in the same order, the sidesets
3062 : // each side would be added to.
3063 : typedef std::pair<dof_id_type, unsigned short int> elem_side_type;
3064 8 : std::unordered_map<processor_id_type, std::vector<elem_side_type>> sides_requested;
3065 : std::unordered_map<processor_id_type, std::vector<std::vector<boundary_id_type>>>
3066 8 : requested_boundary_ids;
3067 :
3068 2572 : for (const auto & elem : _mesh->active_element_ptr_range())
3069 4582 : for (auto side : elem->side_index_range())
3070 : {
3071 : // A side whose neighbor is in the same subdomain is interior to
3072 : // a block, not a true boundary or subdomain interface, even if
3073 : // all of its nodes happen to lie in the nodeset (e.g. on a mesh
3074 : // that is only one element deep in some direction). Skip it
3075 : // if the caller has opted in to that behavior.
3076 3360 : const Elem * neigh = skip_interior_sides ? elem->neighbor_ptr(side) : nullptr;
3077 :
3078 2204 : if (neigh && neigh != remote_elem &&
3079 1346 : neigh->subdomain_id() == elem->subdomain_id())
3080 1954 : continue;
3081 :
3082 2014 : side_elem = &side_builder(*elem, side);
3083 :
3084 : // map from nodeset_id to count for that ID
3085 136 : std::map<boundary_id_type, unsigned> nodesets_node_count;
3086 :
3087 : // For each nodeset that this node is a member of, increment the associated
3088 : // nodeset ID count
3089 5898 : for (const auto & node : side_elem->node_ref_range())
3090 10628 : for (const auto & pr : as_range(_boundary_node_id.equal_range(&node)))
3091 6744 : if (nodeset_list.empty() || nodeset_list.count(pr.second))
3092 3204 : nodesets_node_count[pr.second]++;
3093 :
3094 : // Now check to see what nodeset_counts have the correct
3095 : // number of nodes in them.
3096 : // Note that matched_ids may have been moved from on a prior
3097 : // iteration, so clear() rather than assuming it's empty.
3098 136 : matched_ids.clear();
3099 3768 : for (const auto & pr : nodesets_node_count)
3100 1754 : if (pr.second == side_elem->n_nodes())
3101 1474 : matched_ids.push_back(pr.first);
3102 :
3103 2014 : if (matched_ids.empty())
3104 514 : continue;
3105 :
3106 : // We don't yet know whether this side is truly interior;
3107 : // defer the decision until we've resolved neigh's real
3108 : // subdomain id.
3109 1474 : if (neigh == remote_elem)
3110 : {
3111 68 : const processor_id_type pid = elem->processor_id();
3112 68 : sides_requested[pid].emplace_back(elem->id(), side);
3113 0 : requested_boundary_ids[pid].push_back(std::move(matched_ids));
3114 0 : continue;
3115 68 : }
3116 :
3117 1406 : add_matched_side(elem, side, matched_ids);
3118 216 : } // end for side
3119 :
3120 : // If we're not resolving RemoteElem neighbors then nothing was
3121 : // deferred, and we must not communicate: see need_parallel_sync
3122 : // above.
3123 232 : if (!need_parallel_sync)
3124 : {
3125 6 : libmesh_assert(sides_requested.empty());
3126 6 : return;
3127 : }
3128 :
3129 : // Even if this processor has no deferred sides of its own, it may
3130 : // still own elements that other processors need to query, so every
3131 : // processor has to participate in the collective communication
3132 : // below rather than skipping it when its own request list is empty.
3133 :
3134 : auto neighbor_subdomain_gather_functor =
3135 68 : [this]
3136 : (processor_id_type,
3137 : const std::vector<elem_side_type> & elem_sides,
3138 68 : std::vector<subdomain_id_type> & neighbor_subdomains)
3139 : {
3140 68 : neighbor_subdomains.resize(elem_sides.size());
3141 136 : for (auto i : index_range(elem_sides))
3142 : {
3143 68 : const Elem & elem = _mesh->elem_ref(elem_sides[i].first);
3144 68 : const Elem * neigh = elem.neighbor_ptr(elem_sides[i].second);
3145 :
3146 : // The requester saw a RemoteElem neighbor here, so there is
3147 : // a neighbor, and we own elem, so our connectivity for it is
3148 : // guaranteed to be fully resolved.
3149 0 : libmesh_assert(neigh);
3150 0 : libmesh_assert(neigh != remote_elem);
3151 :
3152 68 : neighbor_subdomains[i] = neigh->subdomain_id();
3153 : }
3154 201 : };
3155 :
3156 : auto neighbor_subdomain_action_functor =
3157 68 : [this, &requested_boundary_ids, &add_matched_side]
3158 : (processor_id_type pid,
3159 : const std::vector<elem_side_type> & elem_sides,
3160 0 : const std::vector<subdomain_id_type> & neighbor_subdomains)
3161 : {
3162 68 : const auto & bids = libmesh_map_find(requested_boundary_ids, pid);
3163 0 : libmesh_assert_equal_to(bids.size(), elem_sides.size());
3164 :
3165 136 : for (auto i : index_range(elem_sides))
3166 : {
3167 68 : const Elem * elem = _mesh->elem_ptr(elem_sides[i].first);
3168 :
3169 : // If the side turned out to be truly interior after all,
3170 : // skip it.
3171 68 : if (neighbor_subdomains[i] != elem->subdomain_id())
3172 0 : add_matched_side(elem, elem_sides[i].second, bids[i]);
3173 : }
3174 133 : };
3175 :
3176 2 : subdomain_id_type * datum_type_ex = nullptr;
3177 : Parallel::pull_parallel_vector_data
3178 135 : (this->comm(), sides_requested, neighbor_subdomain_gather_functor,
3179 : neighbor_subdomain_action_functor, datum_type_ex);
3180 : }
3181 :
3182 :
3183 :
3184 :
3185 : std::vector<BoundaryInfo::BCTuple>
3186 22668 : BoundaryInfo::build_side_list(BCTupleSortBy sort_by) const
3187 : {
3188 861 : std::vector<BCTuple> bc_triples;
3189 22668 : bc_triples.reserve(_boundary_side_id.size());
3190 :
3191 628534 : for (const auto & [elem, id_pair] : _boundary_side_id)
3192 605866 : bc_triples.emplace_back(elem->id(), id_pair.first, id_pair.second);
3193 :
3194 : // bc_triples is currently in whatever order the Elem pointers in
3195 : // the _boundary_side_id multimap are in, and in particular might be
3196 : // in different orders on different processors. To avoid this
3197 : // inconsistency, we'll sort using the default operator< for tuples.
3198 22668 : if (sort_by == BCTupleSortBy::ELEM_ID)
3199 22668 : std::sort(bc_triples.begin(), bc_triples.end());
3200 0 : else if (sort_by == BCTupleSortBy::SIDE_ID)
3201 0 : std::sort(bc_triples.begin(), bc_triples.end(),
3202 0 : [](const BCTuple & left, const BCTuple & right)
3203 0 : {return std::get<1>(left) < std::get<1>(right);});
3204 0 : else if (sort_by == BCTupleSortBy::BOUNDARY_ID)
3205 0 : std::sort(bc_triples.begin(), bc_triples.end(),
3206 0 : [](const BCTuple & left, const BCTuple & right)
3207 0 : {return std::get<2>(left) < std::get<2>(right);});
3208 :
3209 22668 : return bc_triples;
3210 : }
3211 :
3212 :
3213 :
3214 : std::vector<BoundaryInfo::BCTuple>
3215 0 : BoundaryInfo::build_active_side_list () const
3216 : {
3217 0 : std::vector<BCTuple> bc_triples;
3218 0 : bc_triples.reserve(_boundary_side_id.size());
3219 :
3220 0 : for (const auto & [elem, id_pair] : _boundary_side_id)
3221 : {
3222 : // Don't add remote sides
3223 0 : if (elem->is_remote())
3224 0 : continue;
3225 :
3226 : // Loop over the sides of possible children
3227 0 : std::vector<const Elem *> family;
3228 : #ifdef LIBMESH_ENABLE_AMR
3229 0 : elem->active_family_tree_by_side(family, id_pair.first);
3230 : #else
3231 : family.push_back(elem);
3232 : #endif
3233 :
3234 : // Populate the list items
3235 0 : for (const auto & f : family)
3236 0 : bc_triples.emplace_back(f->id(), id_pair.first, id_pair.second);
3237 : }
3238 :
3239 : // This list is currently in memory address (arbitrary) order, so
3240 : // sort to make it consistent on all procs.
3241 0 : std::sort(bc_triples.begin(), bc_triples.end());
3242 :
3243 0 : return bc_triples;
3244 : }
3245 :
3246 :
3247 : std::vector<BoundaryInfo::BCTuple>
3248 5224 : BoundaryInfo::build_edge_list() const
3249 : {
3250 368 : std::vector<BCTuple> bc_triples;
3251 5224 : bc_triples.reserve(_boundary_edge_id.size());
3252 :
3253 23140 : for (const auto & [elem, id_pair] : _boundary_edge_id)
3254 17916 : bc_triples.emplace_back(elem->id(), id_pair.first, id_pair.second);
3255 :
3256 : // This list is currently in memory address (arbitrary) order, so
3257 : // sort to make it consistent on all procs.
3258 5224 : std::sort(bc_triples.begin(), bc_triples.end());
3259 :
3260 5224 : return bc_triples;
3261 : }
3262 :
3263 :
3264 : std::vector<BoundaryInfo::BCTuple>
3265 5146 : BoundaryInfo::build_shellface_list() const
3266 : {
3267 364 : std::vector<BCTuple> bc_triples;
3268 5146 : bc_triples.reserve(_boundary_shellface_id.size());
3269 :
3270 45082 : for (const auto & [elem, id_pair] : _boundary_shellface_id)
3271 39936 : bc_triples.emplace_back(elem->id(), id_pair.first, id_pair.second);
3272 :
3273 : // This list is currently in memory address (arbitrary) order, so
3274 : // sort to make it consistent on all procs.
3275 5146 : std::sort(bc_triples.begin(), bc_triples.end());
3276 :
3277 5146 : return bc_triples;
3278 : }
3279 :
3280 :
3281 0 : void BoundaryInfo::print_info(std::ostream & out_stream) const
3282 : {
3283 : // Print out the nodal BCs
3284 0 : if (!_boundary_node_id.empty())
3285 : {
3286 0 : out_stream << "Nodal Boundary conditions:" << std::endl
3287 0 : << "--------------------------" << std::endl
3288 0 : << " (Node No., ID) " << std::endl;
3289 :
3290 0 : for (const auto & [node, bndry_id] : _boundary_node_id)
3291 0 : out_stream << " (" << node->id()
3292 0 : << ", " << bndry_id
3293 0 : << ")" << std::endl;
3294 : }
3295 :
3296 : // Print out the element edge BCs
3297 0 : if (!_boundary_edge_id.empty())
3298 : {
3299 0 : out_stream << std::endl
3300 0 : << "Edge Boundary conditions:" << std::endl
3301 0 : << "-------------------------" << std::endl
3302 0 : << " (Elem No., Edge No., ID) " << std::endl;
3303 :
3304 0 : for (const auto & [elem, id_pair] : _boundary_edge_id)
3305 0 : out_stream << " (" << elem->id()
3306 0 : << ", " << id_pair.first
3307 0 : << ", " << id_pair.second
3308 0 : << ")" << std::endl;
3309 : }
3310 :
3311 : // Print out the element shellface BCs
3312 0 : if (!_boundary_shellface_id.empty())
3313 : {
3314 0 : out_stream << std::endl
3315 0 : << "Shell-face Boundary conditions:" << std::endl
3316 0 : << "-------------------------" << std::endl
3317 0 : << " (Elem No., Shell-face No., ID) " << std::endl;
3318 :
3319 0 : for (const auto & [elem, id_pair] : _boundary_shellface_id)
3320 0 : out_stream << " (" << elem->id()
3321 0 : << ", " << id_pair.first
3322 0 : << ", " << id_pair.second
3323 0 : << ")" << std::endl;
3324 : }
3325 :
3326 : // Print out the element side BCs
3327 0 : if (!_boundary_side_id.empty())
3328 : {
3329 0 : out_stream << std::endl
3330 0 : << "Side Boundary conditions:" << std::endl
3331 0 : << "-------------------------" << std::endl
3332 0 : << " (Elem No., Side No., ID) " << std::endl;
3333 :
3334 0 : for (const auto & [elem, id_pair] : _boundary_side_id)
3335 0 : out_stream << " (" << elem->id()
3336 0 : << ", " << id_pair.first
3337 0 : << ", " << id_pair.second
3338 0 : << ")" << std::endl;
3339 : }
3340 0 : }
3341 :
3342 :
3343 :
3344 0 : void BoundaryInfo::print_summary(std::ostream & out_stream) const
3345 : {
3346 : // Print out the nodal BCs
3347 0 : if (!_boundary_node_id.empty())
3348 : {
3349 0 : out_stream << "Nodal Boundary conditions:" << std::endl
3350 0 : << "--------------------------" << std::endl
3351 0 : << " (ID, number of nodes) " << std::endl;
3352 :
3353 0 : std::map<boundary_id_type, std::size_t> ID_counts;
3354 :
3355 0 : for (const auto & pr : _boundary_node_id)
3356 0 : ID_counts[pr.second]++;
3357 :
3358 0 : for (const auto & [bndry_id, cnt] : ID_counts)
3359 0 : out_stream << " (" << bndry_id
3360 0 : << ", " << cnt
3361 0 : << ")" << std::endl;
3362 : }
3363 :
3364 : // Print out the element edge BCs
3365 0 : if (!_boundary_edge_id.empty())
3366 : {
3367 0 : out_stream << std::endl
3368 0 : << "Edge Boundary conditions:" << std::endl
3369 0 : << "-------------------------" << std::endl
3370 0 : << " (ID, number of edges) " << std::endl;
3371 :
3372 0 : std::map<boundary_id_type, std::size_t> ID_counts;
3373 :
3374 0 : for (const auto & pr : _boundary_edge_id)
3375 0 : ID_counts[pr.second.second]++;
3376 :
3377 0 : for (const auto & [bndry_id, cnt] : ID_counts)
3378 0 : out_stream << " (" << bndry_id
3379 0 : << ", " << cnt
3380 0 : << ")" << std::endl;
3381 : }
3382 :
3383 :
3384 : // Print out the element edge BCs
3385 0 : if (!_boundary_shellface_id.empty())
3386 : {
3387 0 : out_stream << std::endl
3388 0 : << "Shell-face Boundary conditions:" << std::endl
3389 0 : << "-------------------------" << std::endl
3390 0 : << " (ID, number of shellfaces) " << std::endl;
3391 :
3392 0 : std::map<boundary_id_type, std::size_t> ID_counts;
3393 :
3394 0 : for (const auto & pr : _boundary_shellface_id)
3395 0 : ID_counts[pr.second.second]++;
3396 :
3397 0 : for (const auto & [bndry_id, cnt] : ID_counts)
3398 0 : out_stream << " (" << bndry_id
3399 0 : << ", " << cnt
3400 0 : << ")" << std::endl;
3401 : }
3402 :
3403 : // Print out the element side BCs
3404 0 : if (!_boundary_side_id.empty())
3405 : {
3406 0 : out_stream << std::endl
3407 0 : << "Side Boundary conditions:" << std::endl
3408 0 : << "-------------------------" << std::endl
3409 0 : << " (ID, number of sides) " << std::endl;
3410 :
3411 0 : std::map<boundary_id_type, std::size_t> ID_counts;
3412 :
3413 0 : for (const auto & pr : _boundary_side_id)
3414 0 : ID_counts[pr.second.second]++;
3415 :
3416 0 : for (const auto & [bndry_id, cnt] : ID_counts)
3417 0 : out_stream << " (" << bndry_id
3418 0 : << ", " << cnt
3419 0 : << ")" << std::endl;
3420 : }
3421 0 : }
3422 :
3423 :
3424 31539 : const std::string & BoundaryInfo::get_sideset_name(boundary_id_type id) const
3425 : {
3426 31539 : static const std::string empty_string;
3427 31539 : if (const auto it = _ss_id_to_name.find(id);
3428 1708 : it == _ss_id_to_name.end())
3429 241 : return empty_string;
3430 : else
3431 26448 : return it->second;
3432 : }
3433 :
3434 :
3435 1671701 : std::string & BoundaryInfo::sideset_name(boundary_id_type id)
3436 : {
3437 1671701 : return _ss_id_to_name[id];
3438 : }
3439 :
3440 25643 : const std::string & BoundaryInfo::get_nodeset_name(boundary_id_type id) const
3441 : {
3442 25643 : static const std::string empty_string;
3443 25643 : if (const auto it = _ns_id_to_name.find(id);
3444 1471 : it == _ns_id_to_name.end())
3445 40 : return empty_string;
3446 : else
3447 25170 : return it->second;
3448 : }
3449 :
3450 1669326 : std::string & BoundaryInfo::nodeset_name(boundary_id_type id)
3451 : {
3452 1669326 : return _ns_id_to_name[id];
3453 : }
3454 :
3455 269 : const std::string & BoundaryInfo::get_edgeset_name(boundary_id_type id) const
3456 : {
3457 269 : static const std::string empty_string;
3458 269 : if (const auto it = _es_id_to_name.find(id);
3459 23 : it == _es_id_to_name.end())
3460 21 : return empty_string;
3461 : else
3462 24 : return it->second;
3463 : }
3464 :
3465 :
3466 934 : std::string & BoundaryInfo::edgeset_name(boundary_id_type id)
3467 : {
3468 934 : return _es_id_to_name[id];
3469 : }
3470 :
3471 2382 : boundary_id_type BoundaryInfo::get_id_by_name(std::string_view name) const
3472 : {
3473 : // Search sidesets
3474 5955 : for (const auto & [ss_id, ss_name] : _ss_id_to_name)
3475 5887 : if (ss_name == name)
3476 2382 : return ss_id;
3477 :
3478 : // Search nodesets
3479 0 : for (const auto & [ns_id, ns_name] : _ns_id_to_name)
3480 0 : if (ns_name == name)
3481 0 : return ns_id;
3482 :
3483 : // Search edgesets
3484 0 : for (const auto & [es_id, es_name] : _es_id_to_name)
3485 0 : if (es_name == name)
3486 0 : return es_id;
3487 :
3488 : // If we made it here without returning, we don't have a sideset,
3489 : // nodeset, or edgeset by the requested name, so return invalid_id
3490 0 : return invalid_id;
3491 : }
3492 :
3493 :
3494 :
3495 5893 : void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & requested_boundary_ids,
3496 : dof_id_type first_free_node_id,
3497 : std::map<dof_id_type, dof_id_type> * node_id_map,
3498 : dof_id_type first_free_elem_id,
3499 : std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> * side_id_map,
3500 : const std::set<subdomain_id_type> & subdomains_relative_to)
3501 : {
3502 : // We'll use the standard DistributedMesh trick of dividing up id
3503 : // ranges by processor_id to avoid picking duplicate node ids.
3504 : dof_id_type
3505 6059 : next_node_id = first_free_node_id + this->processor_id();
3506 :
3507 : // We have to be careful how we select different element ids on
3508 : // different processors, because we may be adding sides of refined
3509 : // elements and sides of their parents (which are thereby also
3510 : // parents of their sides), but libMesh convention requires (and
3511 : // some libMesh communication code assumes) that parent elements
3512 : // have lower ids than their children.
3513 : //
3514 : // We'll start with temporary ids in a temporary map stratefied by
3515 : // element level, so we can sort by element level after.
3516 : // Make sure we're sorting by ids here, not anything that might
3517 : // differ from processor to processor, so we can do an
3518 : // embarrassingly parallel sort. This is a serialized data
3519 : // structure, but it's at worst O(N^2/3) so we'll hold off on doing
3520 : // this distributed until someone's benchmark screams at us.
3521 : std::vector
3522 : <std::set<std::pair<dof_id_type, unsigned char>>>
3523 498 : side_id_set_by_level;
3524 :
3525 : // For avoiding extraneous element side construction
3526 332 : ElemSideBuilder side_builder;
3527 : // Pull objects out of the loop to reduce heap operations
3528 : const Elem * side;
3529 :
3530 : // We'll pass through the mesh once first to build
3531 : // the maps and count boundary nodes and elements.
3532 : // To find local boundary nodes, we have to examine all elements
3533 : // here rather than just local elements, because it's possible to
3534 : // have a local boundary node that's not on a local boundary
3535 : // element, e.g. at the tip of a triangle.
3536 :
3537 : // We'll loop through two different ranges here: first all elements,
3538 : // looking for local nodes, and second through unpartitioned
3539 : // elements, looking for all remaining nodes.
3540 6059 : const MeshBase::const_element_iterator end_el = _mesh->elements_end();
3541 166 : bool hit_end_el = false;
3542 : const MeshBase::const_element_iterator end_unpartitioned_el =
3543 6059 : _mesh->pid_elements_end(DofObject::invalid_processor_id);
3544 :
3545 21256 : for (MeshBase::const_element_iterator el = _mesh->elements_begin();
3546 174500 : !hit_end_el || (el != end_unpartitioned_el); ++el)
3547 : {
3548 165196 : if ((el == end_el) && !hit_end_el)
3549 : {
3550 : // Note that we're done with local nodes and just looking
3551 : // for remaining unpartitioned nodes
3552 166 : hit_end_el = true;
3553 :
3554 : // Join up the local results from other processors
3555 5893 : if (side_id_map)
3556 : {
3557 3337 : std::size_t n_levels = side_id_set_by_level.size();
3558 3337 : this->comm().max(n_levels);
3559 3337 : side_id_set_by_level.resize(n_levels);
3560 7384 : for (auto l : make_range(n_levels))
3561 4161 : this->comm().set_union(side_id_set_by_level[l]);
3562 : }
3563 5893 : if (node_id_map)
3564 2556 : this->comm().set_union(*node_id_map);
3565 :
3566 : // Finally we'll pass through any unpartitioned elements to add them
3567 : // to the maps and counts.
3568 5893 : next_node_id = first_free_node_id + this->n_processors();
3569 :
3570 11620 : el = _mesh->pid_elements_begin(DofObject::invalid_processor_id);
3571 5893 : if (el == end_unpartitioned_el)
3572 166 : break;
3573 : }
3574 :
3575 159303 : const Elem * elem = *el;
3576 :
3577 : // If the subdomains_relative_to container has the
3578 : // invalid_subdomain_id, we fall back on the "old" behavior of
3579 : // adding sides regardless of this Elem's subdomain. Otherwise,
3580 : // if the subdomains_relative_to container doesn't contain the
3581 : // current Elem's subdomain_id(), we won't add any sides from
3582 : // it.
3583 19952 : if (!subdomains_relative_to.count(Elem::invalid_subdomain_id) &&
3584 154855 : !subdomains_relative_to.count(elem->subdomain_id()))
3585 10424 : continue;
3586 :
3587 852357 : for (auto s : elem->side_index_range())
3588 : {
3589 41644 : bool add_this_side = false;
3590 :
3591 : // Find all the boundary side ids for this Elem side.
3592 83288 : std::vector<boundary_id_type> bcids;
3593 704294 : this->boundary_ids(elem, s, bcids);
3594 :
3595 752829 : for (const boundary_id_type bcid : bcids)
3596 : {
3597 : // if the user wants this id, we want this side
3598 5046 : if (requested_boundary_ids.count(bcid))
3599 : {
3600 2822 : add_this_side = true;
3601 2822 : break;
3602 : }
3603 : }
3604 :
3605 : // We may still want to add this side if the user called
3606 : // sync() with no requested_boundary_ids. This corresponds
3607 : // to the "old" style of calling sync() in which the entire
3608 : // boundary was copied to the BoundaryMesh, and handles the
3609 : // case where elements on the geometric boundary are not in
3610 : // any sidesets.
3611 449476 : if (requested_boundary_ids.count(invalid_id) &&
3612 407832 : elem->neighbor_ptr(s) == nullptr)
3613 1120 : add_this_side = true;
3614 :
3615 687888 : if (add_this_side)
3616 : {
3617 : // We only assign ids for our own and for
3618 : // unpartitioned elements
3619 47827 : if (side_id_map &&
3620 16952 : ((elem->processor_id() == this->processor_id()) ||
3621 756 : (elem->processor_id() ==
3622 : DofObject::invalid_processor_id)))
3623 : {
3624 1512 : std::pair<dof_id_type, unsigned char> side_pair(elem->id(), s);
3625 9072 : auto level = elem->level();
3626 9828 : if (side_id_set_by_level.size() <= level)
3627 2343 : side_id_set_by_level.resize(level+1);
3628 1512 : auto & level_side_id_set = side_id_set_by_level[level];
3629 756 : libmesh_assert (!level_side_id_set.count(side_pair));
3630 8316 : level_side_id_set.insert(side_pair);
3631 : }
3632 :
3633 46315 : side = &side_builder(*elem, s);
3634 279531 : for (auto n : side->node_index_range())
3635 : {
3636 14410 : const Node & node = side->node_ref(n);
3637 :
3638 : // In parallel we don't know enough to number
3639 : // others' nodes ourselves.
3640 247626 : if (!hit_end_el &&
3641 28820 : (node.processor_id() != this->processor_id()))
3642 146756 : continue;
3643 :
3644 86460 : dof_id_type node_id = node.id();
3645 86460 : if (node_id_map && !node_id_map->count(node_id))
3646 : {
3647 20364 : (*node_id_map)[node_id] = next_node_id;
3648 22061 : next_node_id += this->n_processors() + 1;
3649 : }
3650 : }
3651 : }
3652 : }
3653 : }
3654 :
3655 : // FIXME: should we renumber node_id_map's image to be contiguous
3656 : // here, rather than waiting for renumbering in mesh preparation to
3657 : // do it later?
3658 :
3659 : // We do need to build up side_id_map here, to handle proper sorting
3660 : // by level.
3661 5893 : if (side_id_map)
3662 : {
3663 94 : dof_id_type next_elem_id = first_free_elem_id;
3664 7384 : for (auto level : make_range(side_id_set_by_level.size()))
3665 : {
3666 57723 : for (auto side_pair : side_id_set_by_level[level])
3667 53676 : (*side_id_map)[side_pair] = next_elem_id++;
3668 : }
3669 : }
3670 5893 : }
3671 :
3672 1633 : void BoundaryInfo::clear_stitched_boundary_side_ids (const boundary_id_type sideset_id,
3673 : const boundary_id_type other_sideset_id,
3674 : const bool clear_nodeset_data)
3675 : {
3676 46 : auto end_it = _boundary_side_id.end();
3677 46 : auto it = _boundary_side_id.begin();
3678 :
3679 : // This predicate checks to see whether the pred_pr triplet's boundary ID matches sideset_id
3680 : // (other_sideset_id) *and* whether there is a boundary information triplet on the other side of
3681 : // the face whose boundary ID matches the other_sideset_id (sideset_id). We return a pair where
3682 : // first is a boolean indicating our condition is true or false, and second is an iterator to the
3683 : // neighboring triplet if our condition is true
3684 : auto predicate =
3685 71556 : [sideset_id, other_sideset_id](
3686 : const std::pair<const Elem *, std::pair<unsigned short int, boundary_id_type>> & pred_pr,
3687 : const std::multimap<const Elem *, std::pair<unsigned short int, boundary_id_type>> &
3688 7540 : pred_container) {
3689 75828 : const Elem & elem = *pred_pr.first;
3690 75828 : const auto elem_side = pred_pr.second.first;
3691 75828 : const Elem * const other_elem = elem.neighbor_ptr(elem_side);
3692 75828 : if (!other_elem)
3693 1896 : return std::make_pair(false, pred_container.end());
3694 :
3695 8520 : const auto elem_side_bnd_id = pred_pr.second.second;
3696 240 : auto other_elem_side_bnd_id = BoundaryInfo::invalid_id;
3697 8520 : if (elem_side_bnd_id == sideset_id)
3698 153 : other_elem_side_bnd_id = other_sideset_id;
3699 3514 : else if (elem_side_bnd_id == other_sideset_id)
3700 71 : other_elem_side_bnd_id = sideset_id;
3701 : else
3702 16 : return std::make_pair(false, pred_container.end());
3703 :
3704 7952 : const auto other_elem_side = other_elem->which_neighbor_am_i(&elem);
3705 : const typename std::decay<decltype(pred_container)>::type::value_type other_sideset_info(
3706 448 : other_elem, std::make_pair(other_elem_side, other_elem_side_bnd_id));
3707 448 : auto other_range = pred_container.equal_range(other_elem);
3708 224 : libmesh_assert_msg(
3709 : other_range.first != other_range.second,
3710 : "No matching sideset information for other element in boundary information");
3711 7952 : auto other_it = std::find(other_range.first, other_range.second, other_sideset_info);
3712 224 : libmesh_assert_msg(
3713 : other_it != pred_container.end(),
3714 : "No matching sideset information for other element in boundary information");
3715 224 : return std::make_pair(true, other_it);
3716 1633 : };
3717 :
3718 77461 : for (; it != end_it;)
3719 : {
3720 77964 : auto pred_result = predicate(*it, _boundary_side_id);
3721 75828 : if (pred_result.first)
3722 : {
3723 : // First erase associated nodeset information. Do it from both
3724 : // sides, so we get any higher-order nodes if we're looking at
3725 : // them from a lower-order side, and so we only remove the two
3726 : // boundary ids used for stitching.
3727 7952 : if (clear_nodeset_data)
3728 : {
3729 7952 : const Elem & elem = *it->first;
3730 7952 : const Elem & neigh = *pred_result.second->first;
3731 7952 : const auto elem_side = it->second.first;
3732 7952 : const boundary_id_type neigh_side = pred_result.second->second.first;
3733 7952 : const auto elem_bcid = it->second.second;
3734 7952 : const boundary_id_type neigh_bcid = pred_result.second->second.second;
3735 :
3736 55083 : for (const auto local_node_num : elem.nodes_on_side(elem_side))
3737 48229 : this->remove_node(elem.node_ptr(local_node_num), elem_bcid);
3738 :
3739 55060 : for (const auto local_node_num : neigh.nodes_on_side(neigh_side))
3740 48204 : this->remove_node(neigh.node_ptr(local_node_num), neigh_bcid);
3741 : }
3742 :
3743 : // Now erase the sideset information for our element and its
3744 : // neighbor, together. This is safe since a multimap doesn't
3745 : // invalidate iterators.
3746 7728 : _boundary_side_id.erase(pred_result.second);
3747 7728 : it = _boundary_side_id.erase(it);
3748 : }
3749 : else
3750 1912 : ++it;
3751 : }
3752 :
3753 : // Removing stitched-away boundary ids might have removed an id
3754 : // *entirely*, so we need to recompute boundary id sets to check
3755 : // for that.
3756 1633 : this->regenerate_id_sets();
3757 1633 : this->libmesh_assert_valid_multimaps();
3758 1633 : }
3759 :
3760 : const std::set<boundary_id_type> &
3761 568 : BoundaryInfo::get_global_boundary_ids() const
3762 : {
3763 16 : libmesh_assert(_mesh && _mesh->_preparation.has_boundary_id_sets);
3764 568 : return _global_boundary_ids;
3765 : }
3766 :
3767 : void
3768 4015 : BoundaryInfo::libmesh_assert_valid_multimaps() const
3769 : {
3770 : #ifndef NDEBUG
3771 456 : auto verify_multimap = [](const auto & themap) {
3772 13772 : for (const auto & [key, val] : themap)
3773 : {
3774 13316 : auto range = themap.equal_range(key);
3775 :
3776 13316 : int count = 0;
3777 41436 : for (auto it = range.first; it != range.second; ++it)
3778 28120 : if (it->second == val)
3779 13316 : ++count;
3780 :
3781 13316 : libmesh_assert(count == 1);
3782 : }
3783 456 : };
3784 :
3785 114 : verify_multimap(this->_boundary_node_id);
3786 114 : verify_multimap(this->_boundary_edge_id);
3787 114 : verify_multimap(this->_boundary_shellface_id);
3788 114 : verify_multimap(this->_boundary_side_id);
3789 : #else
3790 3901 : return;
3791 : #endif
3792 114 : }
3793 :
3794 : } // namespace libMesh
|
