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