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