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