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