boundary_info.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // C++ includes
21 #include <iterator> // std::distance
22 
23 // Local includes
24 #include "libmesh/libmesh_config.h"
25 
26 #include "libmesh/boundary_info.h"
28 #include "libmesh/elem.h"
31 #include "libmesh/parallel.h"
32 #include "libmesh/partitioner.h"
33 #include "libmesh/remote_elem.h"
35 
36 namespace libMesh
37 {
38 
39 
40 
41 //------------------------------------------------------
42 // BoundaryInfo static member initializations
44 
45 
46 
47 //------------------------------------------------------
48 // BoundaryInfo functions
50  ParallelObject(m.comm()),
51  _mesh (m)
52 {
53 }
54 
55 BoundaryInfo & BoundaryInfo::operator=(const BoundaryInfo & other_boundary_info)
56 {
57  // Overwrite any preexisting boundary info
58  this->clear();
59 
70  // Copy node boundary info
71  for (const auto & pr : other_boundary_info._boundary_node_id)
72  _boundary_node_id.insert(std::make_pair(_mesh.node_ptr(pr.first->id()),
73  pr.second));
74 
75  // Copy edge boundary info
76  for (const auto & pr : other_boundary_info._boundary_edge_id)
77  _boundary_edge_id.insert(std::make_pair(_mesh.elem_ptr(pr.first->id()),
78  pr.second));
79 
80  // Copy shellface boundary info
81  for (const auto & pr : other_boundary_info._boundary_shellface_id)
82  _boundary_shellface_id.insert(std::make_pair(_mesh.elem_ptr(pr.first->id()),
83  pr.second));
84 
85  // Copy side boundary info
86  for (const auto & pr : other_boundary_info._boundary_side_id)
87  _boundary_side_id.insert(std::make_pair(_mesh.elem_ptr(pr.first->id()),
88  pr.second));
89 
90  _boundary_ids = other_boundary_info._boundary_ids;
91  _side_boundary_ids = other_boundary_info._side_boundary_ids;
92  _node_boundary_ids = other_boundary_info._node_boundary_ids;
93  _edge_boundary_ids = other_boundary_info._edge_boundary_ids;
95 
96  return *this;
97 }
98 
99 
101 {
102  this->clear();
103 }
104 
105 
106 
108 {
109  _boundary_node_id.clear();
110  _boundary_side_id.clear();
111  _boundary_edge_id.clear();
112  _boundary_shellface_id.clear();
113  _boundary_ids.clear();
114  _side_boundary_ids.clear();
115  _node_boundary_ids.clear();
116  _edge_boundary_ids.clear();
117  _shellface_boundary_ids.clear();
118 }
119 
120 
121 
123 {
124  // Clear the old caches
125  _boundary_ids.clear();
126  _side_boundary_ids.clear();
127  _node_boundary_ids.clear();
128  _edge_boundary_ids.clear();
129  _shellface_boundary_ids.clear();
130 
131  // Loop over id maps to regenerate each set.
132  for (const auto & pr : _boundary_node_id)
133  {
134  const boundary_id_type id = pr.second;
135  _boundary_ids.insert(id);
136  _node_boundary_ids.insert(id);
137  }
138 
139  for (const auto & pr : _boundary_edge_id)
140  {
141  const boundary_id_type id = pr.second.second;
142  _boundary_ids.insert(id);
143  _edge_boundary_ids.insert(id);
144  }
145 
146  for (const auto & pr : _boundary_side_id)
147  {
148  const boundary_id_type id = pr.second.second;
149  _boundary_ids.insert(id);
150  _side_boundary_ids.insert(id);
151  }
152 
153  for (const auto & pr : _boundary_shellface_id)
154  {
155  const boundary_id_type id = pr.second.second;
156  _boundary_ids.insert(id);
157  _shellface_boundary_ids.insert(id);
158  }
159 }
160 
161 
162 
163 void BoundaryInfo::sync (UnstructuredMesh & boundary_mesh)
164 {
165  std::set<boundary_id_type> request_boundary_ids(_boundary_ids);
166  request_boundary_ids.insert(invalid_id);
167  if (!_mesh.is_serial())
168  this->comm().set_union(request_boundary_ids);
169 
170  this->sync(request_boundary_ids,
171  boundary_mesh);
172 }
173 
174 
175 
176 void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_ids,
177  UnstructuredMesh & boundary_mesh)
178 {
179  // Call the 3 argument version of this function with a dummy value for the third set.
180  std::set<subdomain_id_type> subdomains_relative_to;
181  subdomains_relative_to.insert(Elem::invalid_subdomain_id);
182 
183  this->sync(requested_boundary_ids,
184  boundary_mesh,
185  subdomains_relative_to);
186 }
187 
188 
189 
190 void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_ids,
191  UnstructuredMesh & boundary_mesh,
192  const std::set<subdomain_id_type> & subdomains_relative_to)
193 {
194  LOG_SCOPE("sync()", "BoundaryInfo");
195 
196  boundary_mesh.clear();
197 
203  if (!_mesh.is_serial())
204  boundary_mesh.delete_remote_elements();
205 
212  MeshSerializer serializer
213  (const_cast<MeshBase &>(_mesh), boundary_mesh.is_serial());
214 
219  boundary_mesh.set_n_partitions() = _mesh.n_partitions();
220 
221  std::map<dof_id_type, dof_id_type> node_id_map;
222  std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> side_id_map;
223 
224  this->_find_id_maps(requested_boundary_ids, 0, &node_id_map, 0, &side_id_map, subdomains_relative_to);
225 
226  // Let's add all the boundary nodes we found to the boundary mesh
227  for (const auto & node : _mesh.node_ptr_range())
228  {
229  dof_id_type node_id = node->id();
230  if (node_id_map.count(node_id))
231  {
232  boundary_mesh.add_point(*node, node_id_map[node_id], node->processor_id());
233 
234  // Copy over all the node's boundary IDs to boundary_mesh
235  std::vector<boundary_id_type> node_boundary_ids;
236  this->boundary_ids(node, node_boundary_ids);
237  for (const auto & node_bid : node_boundary_ids)
238  boundary_mesh.get_boundary_info().add_node(node_id_map[node_id], node_bid);
239  }
240  }
241 
242  // Let's add the elements
243  this->add_elements (requested_boundary_ids, boundary_mesh, subdomains_relative_to);
244 
245  // The new elements are currently using the interior mesh's nodes;
246  // we want them to use the boundary mesh's nodes instead.
247 
248  // This side's Node pointers still point to the nodes of the original mesh.
249  // We need to re-point them to the boundary mesh's nodes! Since we copied *ALL* of
250  // the original mesh's nodes over, we should be guaranteed to have the same ordering.
251  for (auto & new_elem : boundary_mesh.element_ptr_range())
252  {
253  for (auto nn : new_elem->node_index_range())
254  {
255  // Get the correct node pointer, based on the id()
256  Node * new_node =
257  boundary_mesh.node_ptr(node_id_map[new_elem->node_id(nn)]);
258 
259  // sanity check: be sure that the new Node exists and its
260  // global id really matches
261  libmesh_assert (new_node);
262  libmesh_assert_equal_to (new_node->id(),
263  node_id_map[new_elem->node_id(nn)]);
264 
265  // Assign the new node pointer
266  new_elem->set_node(nn) = new_node;
267  }
268  }
269 
270  // Don't repartition this mesh; we want it to stay in sync with the
271  // interior partitioning.
272  boundary_mesh.partitioner().reset(nullptr);
273 
274  // Make boundary_mesh nodes and elements contiguous
275  boundary_mesh.prepare_for_use(/*skip_renumber =*/ false);
276 
277  // and finally distribute element partitioning to the nodes
279 }
280 
281 
283  std::map<dof_id_type, dof_id_type> & node_id_map,
284  std::map<dof_id_type, unsigned char> & side_id_map,
285  Real tolerance)
286 {
287  LOG_SCOPE("get_side_and_node_maps()", "BoundaryInfo");
288 
289  node_id_map.clear();
290  side_id_map.clear();
291 
292  // Pull objects out of the loop to reduce heap operations
293  std::unique_ptr<const Elem> interior_parent_side;
294 
295  for (const auto & boundary_elem : boundary_mesh.active_element_ptr_range())
296  {
297  const Elem * interior_parent = boundary_elem->interior_parent();
298 
299  // Find out which side of interior_parent boundary_elem corresponds to.
300  // Use centroid comparison as a way to check.
301  unsigned char interior_parent_side_index = 0;
302  bool found_matching_sides = false;
303  for (auto side : interior_parent->side_index_range())
304  {
305  interior_parent->build_side_ptr(interior_parent_side, side);
306  Real centroid_distance = (boundary_elem->centroid() - interior_parent_side->centroid()).norm();
307 
308  if (centroid_distance < (tolerance * boundary_elem->hmin()))
309  {
310  interior_parent_side_index = cast_int<unsigned char>(side);
311  found_matching_sides = true;
312  break;
313  }
314  }
315 
316  if (!found_matching_sides)
317  libmesh_error_msg("No matching side found within the specified tolerance");
318 
319  side_id_map[boundary_elem->id()] = interior_parent_side_index;
320 
321  for (auto local_node_index : boundary_elem->node_index_range())
322  {
323  dof_id_type boundary_node_id = boundary_elem->node_id(local_node_index);
324  dof_id_type interior_node_id = interior_parent_side->node_id(local_node_index);
325 
326  node_id_map[interior_node_id] = boundary_node_id;
327  }
328  }
329 }
330 
331 
332 
333 void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_boundary_ids,
334  UnstructuredMesh & boundary_mesh)
335 {
336  // Call the 3 argument version of this function with a dummy value for the third arg.
337  std::set<subdomain_id_type> subdomains_relative_to;
338  subdomains_relative_to.insert(Elem::invalid_subdomain_id);
339 
340  this->add_elements(requested_boundary_ids,
341  boundary_mesh,
342  subdomains_relative_to);
343 }
344 
345 
346 
347 void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_boundary_ids,
348  UnstructuredMesh & boundary_mesh,
349  const std::set<subdomain_id_type> & subdomains_relative_to)
350 {
351  LOG_SCOPE("add_elements()", "BoundaryInfo");
352 
353  // We're not prepared to mix serial and distributed meshes in this
354  // method, so make sure they match from the start.
355  libmesh_assert_equal_to(_mesh.is_serial(),
356  boundary_mesh.is_serial());
357 
358  std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> side_id_map;
359  this->_find_id_maps(requested_boundary_ids,
360  0,
361  nullptr,
362  boundary_mesh.max_elem_id(),
363  &side_id_map,
364  subdomains_relative_to);
365 
366  // We have to add sides *outside* any element loop, because if
367  // boundary_mesh and _mesh are the same then those additions can
368  // invalidate our element iterators. So we just use the element
369  // loop to make a list of sides to add.
370  typedef std::vector<std::pair<dof_id_type, unsigned char>>
371  side_container;
372  side_container sides_to_add;
373 
374  for (const auto & elem : _mesh.element_ptr_range())
375  {
376  // If the subdomains_relative_to container has the
377  // invalid_subdomain_id, we fall back on the "old" behavior of
378  // adding sides regardless of this Elem's subdomain. Otherwise,
379  // if the subdomains_relative_to container doesn't contain the
380  // current Elem's subdomain_id(), we won't add any sides from
381  // it.
382  if (!subdomains_relative_to.count(Elem::invalid_subdomain_id) &&
383  !subdomains_relative_to.count(elem->subdomain_id()))
384  continue;
385 
386  // Get the top-level parent for this element
387  const Elem * top_parent = elem->top_parent();
388 
389  // Find all the boundary side ids for this Elem.
390  auto bounds = _boundary_side_id.equal_range(top_parent);
391 
392  for (auto s : elem->side_index_range())
393  {
394  bool add_this_side = false;
395  boundary_id_type this_bcid = invalid_id;
396 
397  for (const auto & pr : as_range(bounds))
398  {
399  this_bcid = pr.second.second;
400 
401  // if this side is flagged with a boundary condition
402  // and the user wants this id
403  if ((pr.second.first == s) &&
404  (requested_boundary_ids.count(this_bcid)))
405  {
406  add_this_side = true;
407  break;
408  }
409  }
410 
411  // We may still want to add this side if the user called
412  // sync() with no requested_boundary_ids. This corresponds
413  // to the "old" style of calling sync() in which the entire
414  // boundary was copied to the BoundaryMesh, and handles the
415  // case where elements on the geometric boundary are not in
416  // any sidesets.
417  if (bounds.first == bounds.second &&
418  requested_boundary_ids.count(invalid_id) &&
419  elem->neighbor_ptr(s) == nullptr)
420  add_this_side = true;
421 
422  if (add_this_side)
423  sides_to_add.push_back(std::make_pair(elem->id(), s));
424  }
425  }
426 
427 #ifdef LIBMESH_ENABLE_UNIQUE_ID
428  unique_id_type old_max_unique_id = boundary_mesh.parallel_max_unique_id();
429 #endif
430 
431  for (const auto & pr : sides_to_add)
432  {
433  const dof_id_type elem_id = pr.first;
434  const unsigned char s = pr.second;
435  Elem * elem = _mesh.elem_ptr(elem_id);
436 
437  // Build the side - do not use a "proxy" element here:
438  // This will be going into the boundary_mesh and needs to
439  // stand on its own.
440  std::unique_ptr<Elem> side (elem->build_side_ptr(s, false));
441 
442  side->processor_id() = elem->processor_id();
443 
444  const std::pair<dof_id_type, unsigned char> side_pair(elem_id, s);
445 
446  libmesh_assert(side_id_map.count(side_pair));
447 
448  const dof_id_type new_side_id = side_id_map[side_pair];
449 
450  side->set_id(new_side_id);
451 
452 #ifdef LIBMESH_ENABLE_UNIQUE_ID
453  side->set_unique_id() = old_max_unique_id + new_side_id;
454 #endif
455 
456  // Add the side
457  Elem * new_elem = boundary_mesh.add_elem(side.release());
458 
459 #ifdef LIBMESH_ENABLE_AMR
460  // Set parent links
461  if (elem->parent())
462  {
463  const std::pair<dof_id_type, unsigned char> parent_side_pair(elem->parent()->id(), s);
464 
465  libmesh_assert(side_id_map.count(parent_side_pair));
466 
467  Elem * side_parent = boundary_mesh.elem_ptr(side_id_map[parent_side_pair]);
468 
469  libmesh_assert(side_parent);
470 
471  new_elem->set_parent(side_parent);
472 
473  side_parent->set_refinement_flag(Elem::INACTIVE);
474 
475  // Figuring out which child we are of our parent
476  // is a trick. Due to libMesh child numbering
477  // conventions, if we are an element on a vertex,
478  // then we share that vertex with our parent, with
479  // the same local index.
480  bool found_child = false;
481  for (unsigned int v=0; v != new_elem->n_vertices(); ++v)
482  if (new_elem->node_ptr(v) == side_parent->node_ptr(v))
483  {
484  side_parent->add_child(new_elem, v);
485  found_child = true;
486  }
487 
488  // If we don't share any vertex with our parent,
489  // then we're the fourth child (index 3) of a
490  // triangle.
491  if (!found_child)
492  {
493  libmesh_assert_equal_to (new_elem->n_vertices(), 3);
494  side_parent->add_child(new_elem, 3);
495  }
496  }
497 #endif
498 
499  new_elem->set_interior_parent (const_cast<Elem *>(elem));
500 
501  // On non-local elements on DistributedMesh we might have
502  // RemoteElem neighbor links to construct
503  if (!_mesh.is_serial() &&
504  (elem->processor_id() != this->processor_id()))
505  {
506  const unsigned short n_nodes = elem->n_nodes();
507 
508  const unsigned short bdy_n_sides = new_elem->n_sides();
509  const unsigned short bdy_n_nodes = new_elem->n_nodes();
510 
511  // Check every interior side for a RemoteElem
512  for (auto interior_side : elem->side_index_range())
513  {
514  // Might this interior side have a RemoteElem that
515  // needs a corresponding Remote on a boundary side?
516  if (elem->neighbor_ptr(interior_side) != remote_elem)
517  continue;
518 
519  // Which boundary side?
520  for (unsigned short boundary_side = 0;
521  boundary_side != bdy_n_sides; ++boundary_side)
522  {
523  // Look for matching node points. This is safe in
524  // *this* context.
525  bool found_all_nodes = true;
526  for (unsigned short boundary_node = 0;
527  boundary_node != bdy_n_nodes; ++boundary_node)
528  {
529  if (!new_elem->is_node_on_side(boundary_node,
530  boundary_side))
531  continue;
532 
533  bool found_this_node = false;
534  for (unsigned short interior_node = 0;
535  interior_node != n_nodes; ++interior_node)
536  {
537  if (!elem->is_node_on_side(interior_node,
538  interior_side))
539  continue;
540 
541  if (new_elem->point(boundary_node) ==
542  elem->point(interior_node))
543  {
544  found_this_node = true;
545  break;
546  }
547  }
548  if (!found_this_node)
549  {
550  found_all_nodes = false;
551  break;
552  }
553  }
554 
555  if (found_all_nodes)
556  {
557  new_elem->set_neighbor
558  (boundary_side,
559  const_cast<RemoteElem *>(remote_elem));
560  break;
561  }
562  }
563  }
564  }
565  }
566 
567  // We haven't been bothering to keep unique ids consistent on ghost
568  // elements
569  if (!boundary_mesh.is_serial())
571 
572  // Make sure we didn't add ids inconsistently
573 #ifdef DEBUG
574 # ifdef LIBMESH_HAVE_RTTI
575  DistributedMesh * parmesh = dynamic_cast<DistributedMesh *>(&boundary_mesh);
576  if (parmesh)
578 # endif
579 #endif
580 }
581 
582 
583 
585  const boundary_id_type id)
586 {
587  const Node * node_ptr = _mesh.query_node_ptr(node_id);
588 
589  // The user could easily ask for an invalid node id, so let's throw
590  // an easy-to-understand error message when this happens.
591  if (!node_ptr)
592  libmesh_error_msg("BoundaryInfo::add_node(): Could not retrieve pointer for node " << node_id << ", no boundary id was added.");
593 
594  this->add_node (node_ptr, id);
595 }
596 
597 
598 
599 void BoundaryInfo::add_node(const Node * node,
600  const boundary_id_type id)
601 {
602  if (id == invalid_id)
603  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
604  << invalid_id \
605  << "\n That is reserved for internal use.");
606 
607  // Don't add the same ID twice
608  for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
609  if (pr.second == id)
610  return;
611 
612  _boundary_node_id.insert(std::make_pair(node, id));
613  _boundary_ids.insert(id);
614  _node_boundary_ids.insert(id); // Also add this ID to the set of node boundary IDs
615 }
616 
617 
618 
619 void BoundaryInfo::add_node(const Node * node,
620  const std::vector<boundary_id_type> & ids)
621 {
622  if (ids.empty())
623  return;
624 
625  libmesh_assert(node);
626 
627  // Don't add the same ID twice
628  auto bounds = _boundary_node_id.equal_range(node);
629 
630  // The entries in the ids vector may be non-unique. If we expected
631  // *lots* of ids, it might be fastest to construct a std::set from
632  // the entries, but for a small number of entries, which is more
633  // typical, it is probably faster to copy the vector and do sort+unique.
634  // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
635  std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
636  std::sort(unique_ids.begin(), unique_ids.end());
637  std::vector<boundary_id_type>::iterator new_end =
638  std::unique(unique_ids.begin(), unique_ids.end());
639 
640  for (auto & id : as_range(unique_ids.begin(), new_end))
641  {
642  if (id == invalid_id)
643  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
644  << invalid_id \
645  << "\n That is reserved for internal use.");
646 
647  bool already_inserted = false;
648  for (const auto & pr : as_range(bounds))
649  if (pr.second == id)
650  {
651  already_inserted = true;
652  break;
653  }
654  if (already_inserted)
655  continue;
656 
657  _boundary_node_id.insert(std::make_pair(node,id));
658  _boundary_ids.insert(id);
659  _node_boundary_ids.insert(id); // Also add this ID to the set of node boundary IDs
660  }
661 }
662 
663 
664 
666 {
667  _boundary_node_id.clear();
668 }
669 
671  const unsigned short int edge,
672  const boundary_id_type id)
673 {
674  this->add_edge (_mesh.elem_ptr(e), edge, id);
675 }
676 
677 
678 
679 void BoundaryInfo::add_edge(const Elem * elem,
680  const unsigned short int edge,
681  const boundary_id_type id)
682 {
683  libmesh_assert(elem);
684 
685  // Only add BCs for level-0 elements.
686  libmesh_assert_equal_to (elem->level(), 0);
687 
688  if (id == invalid_id)
689  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
690  << invalid_id \
691  << "\n That is reserved for internal use.");
692 
693  // Don't add the same ID twice
694  for (const auto & pr : as_range(_boundary_edge_id.equal_range(elem)))
695  if (pr.second.first == edge &&
696  pr.second.second == id)
697  return;
698 
699  _boundary_edge_id.insert(std::make_pair(elem, std::make_pair(edge, id)));
700  _boundary_ids.insert(id);
701  _edge_boundary_ids.insert(id); // Also add this ID to the set of edge boundary IDs
702 }
703 
704 
705 
706 void BoundaryInfo::add_edge(const Elem * elem,
707  const unsigned short int edge,
708  const std::vector<boundary_id_type> & ids)
709 {
710  if (ids.empty())
711  return;
712 
713  libmesh_assert(elem);
714 
715  // Only add BCs for level-0 elements.
716  libmesh_assert_equal_to (elem->level(), 0);
717 
718  // Don't add the same ID twice
719  auto bounds = _boundary_edge_id.equal_range(elem);
720 
721  // The entries in the ids vector may be non-unique. If we expected
722  // *lots* of ids, it might be fastest to construct a std::set from
723  // the entries, but for a small number of entries, which is more
724  // typical, it is probably faster to copy the vector and do sort+unique.
725  // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
726  std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
727  std::sort(unique_ids.begin(), unique_ids.end());
728  std::vector<boundary_id_type>::iterator new_end =
729  std::unique(unique_ids.begin(), unique_ids.end());
730 
731  for (auto & id : as_range(unique_ids.begin(), new_end))
732  {
733  if (id == invalid_id)
734  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
735  << invalid_id \
736  << "\n That is reserved for internal use.");
737 
738  bool already_inserted = false;
739  for (const auto & pr : as_range(bounds))
740  if (pr.second.first == edge &&
741  pr.second.second == id)
742  {
743  already_inserted = true;
744  break;
745  }
746  if (already_inserted)
747  continue;
748 
749  _boundary_edge_id.insert(std::make_pair(elem, std::make_pair(edge, id)));
750  _boundary_ids.insert(id);
751  _edge_boundary_ids.insert(id); // Also add this ID to the set of edge boundary IDs
752  }
753 }
754 
755 
756 
758  const unsigned short int shellface,
759  const boundary_id_type id)
760 {
761  this->add_shellface (_mesh.elem_ptr(e), shellface, id);
762 }
763 
764 
765 
767  const unsigned short int shellface,
768  const boundary_id_type id)
769 {
770  libmesh_assert(elem);
771 
772  // Only add BCs for level-0 elements.
773  libmesh_assert_equal_to (elem->level(), 0);
774 
775  // Shells only have 2 faces
776  libmesh_assert_less(shellface, 2);
777 
778  if (id == invalid_id)
779  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
780  << invalid_id \
781  << "\n That is reserved for internal use.");
782 
783  // Don't add the same ID twice
784  for (const auto & pr : as_range(_boundary_shellface_id.equal_range(elem)))
785  if (pr.second.first == shellface &&
786  pr.second.second == id)
787  return;
788 
789  _boundary_shellface_id.insert(std::make_pair(elem, std::make_pair(shellface, id)));
790  _boundary_ids.insert(id);
791  _shellface_boundary_ids.insert(id); // Also add this ID to the set of shellface boundary IDs
792 }
793 
794 
795 
797  const unsigned short int shellface,
798  const std::vector<boundary_id_type> & ids)
799 {
800  if (ids.empty())
801  return;
802 
803  libmesh_assert(elem);
804 
805  // Only add BCs for level-0 elements.
806  libmesh_assert_equal_to (elem->level(), 0);
807 
808  // Shells only have 2 faces
809  libmesh_assert_less(shellface, 2);
810 
811  // Don't add the same ID twice
812  auto bounds = _boundary_shellface_id.equal_range(elem);
813 
814  // The entries in the ids vector may be non-unique. If we expected
815  // *lots* of ids, it might be fastest to construct a std::set from
816  // the entries, but for a small number of entries, which is more
817  // typical, it is probably faster to copy the vector and do sort+unique.
818  // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
819  std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
820  std::sort(unique_ids.begin(), unique_ids.end());
821  std::vector<boundary_id_type>::iterator new_end =
822  std::unique(unique_ids.begin(), unique_ids.end());
823 
824  for (auto & id : as_range(unique_ids.begin(), new_end))
825  {
826  if (id == invalid_id)
827  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
828  << invalid_id \
829  << "\n That is reserved for internal use.");
830 
831  bool already_inserted = false;
832  for (const auto & pr : as_range(bounds))
833  if (pr.second.first == shellface &&
834  pr.second.second == id)
835  {
836  already_inserted = true;
837  break;
838  }
839  if (already_inserted)
840  continue;
841 
842  _boundary_shellface_id.insert(std::make_pair(elem, std::make_pair(shellface, id)));
843  _boundary_ids.insert(id);
844  _shellface_boundary_ids.insert(id); // Also add this ID to the set of shellface boundary IDs
845  }
846 }
847 
848 
850  const unsigned short int side,
851  const boundary_id_type id)
852 {
853  this->add_side (_mesh.elem_ptr(e), side, id);
854 }
855 
856 
857 
858 void BoundaryInfo::add_side(const Elem * elem,
859  const unsigned short int side,
860  const boundary_id_type id)
861 {
862  libmesh_assert(elem);
863 
864  // Only add BCs for level-0 elements.
865  libmesh_assert_equal_to (elem->level(), 0);
866 
867  if (id == invalid_id)
868  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
869  << invalid_id \
870  << "\n That is reserved for internal use.");
871 
872  // Don't add the same ID twice
873  for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
874  if (pr.second.first == side &&
875  pr.second.second == id)
876  return;
877 
878  _boundary_side_id.insert(std::make_pair(elem, std::make_pair(side, id)));
879  _boundary_ids.insert(id);
880  _side_boundary_ids.insert(id); // Also add this ID to the set of side boundary IDs
881 }
882 
883 
884 
885 void BoundaryInfo::add_side(const Elem * elem,
886  const unsigned short int side,
887  const std::vector<boundary_id_type> & ids)
888 {
889  if (ids.empty())
890  return;
891 
892  libmesh_assert(elem);
893 
894  // Only add BCs for level-0 elements.
895  libmesh_assert_equal_to (elem->level(), 0);
896 
897  // Don't add the same ID twice
898  auto bounds = _boundary_side_id.equal_range(elem);
899 
900  // The entries in the ids vector may be non-unique. If we expected
901  // *lots* of ids, it might be fastest to construct a std::set from
902  // the entries, but for a small number of entries, which is more
903  // typical, it is probably faster to copy the vector and do sort+unique.
904  // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
905  std::vector<boundary_id_type> unique_ids(ids.begin(), ids.end());
906  std::sort(unique_ids.begin(), unique_ids.end());
907  std::vector<boundary_id_type>::iterator new_end =
908  std::unique(unique_ids.begin(), unique_ids.end());
909 
910  for (auto & id : as_range(unique_ids.begin(), new_end))
911  {
912  if (id == invalid_id)
913  libmesh_error_msg("ERROR: You may not set a boundary ID of " \
914  << invalid_id \
915  << "\n That is reserved for internal use.");
916 
917  bool already_inserted = false;
918  for (const auto & pr : as_range(bounds))
919  if (pr.second.first == side && pr.second.second == id)
920  {
921  already_inserted = true;
922  break;
923  }
924  if (already_inserted)
925  continue;
926 
927  _boundary_side_id.insert(std::make_pair(elem, std::make_pair(side, id)));
928  _boundary_ids.insert(id);
929  _side_boundary_ids.insert(id); // Also add this ID to the set of side boundary IDs
930  }
931 }
932 
933 
934 
935 bool BoundaryInfo::has_boundary_id(const Node * const node,
936  const boundary_id_type id) const
937 {
938  for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
939  if (pr.second == id)
940  return true;
941 
942  return false;
943 }
944 
945 
946 
947 #ifdef LIBMESH_ENABLE_DEPRECATED
948 std::vector<boundary_id_type> BoundaryInfo::boundary_ids(const Node * node) const
949 {
950  libmesh_deprecated();
951 
952  std::vector<boundary_id_type> ids;
953  this->boundary_ids(node, ids);
954  return ids;
955 }
956 #endif
957 
958 
959 
960 void BoundaryInfo::boundary_ids (const Node * node,
961  std::vector<boundary_id_type> & vec_to_fill) const
962 {
963  // Clear out any previous contents
964  vec_to_fill.clear();
965 
966  for (const auto & pr : as_range(_boundary_node_id.equal_range(node)))
967  vec_to_fill.push_back(pr.second);
968 }
969 
970 
971 
972 unsigned int BoundaryInfo::n_boundary_ids(const Node * node) const
973 {
974  auto pos = _boundary_node_id.equal_range(node);
975  return cast_int<unsigned int>(std::distance(pos.first, pos.second));
976 }
977 
978 
979 
980 #ifdef LIBMESH_ENABLE_DEPRECATED
981 std::vector<boundary_id_type> BoundaryInfo::edge_boundary_ids (const Elem * const elem,
982  const unsigned short int edge) const
983 {
984  libmesh_deprecated();
985 
986  std::vector<boundary_id_type> ids;
987  this->edge_boundary_ids(elem, edge, ids);
988  return ids;
989 }
990 #endif
991 
992 
993 
994 void BoundaryInfo::edge_boundary_ids (const Elem * const elem,
995  const unsigned short int edge,
996  std::vector<boundary_id_type> & vec_to_fill) const
997 {
998  libmesh_assert(elem);
999 
1000  // Clear out any previous contents
1001  vec_to_fill.clear();
1002 
1003  // Only level-0 elements store BCs. If this is not a level-0
1004  // element get its level-0 parent and infer the BCs.
1005  const Elem * searched_elem = elem;
1006 #ifdef LIBMESH_ENABLE_AMR
1007  if (elem->level() != 0)
1008  {
1009  // Find all the sides that contain edge. If one of those is a boundary
1010  // side, then this must be a boundary edge. In that case, we just use the
1011  // top-level parent.
1012  bool found_boundary_edge = false;
1013  for (auto side : elem->side_index_range())
1014  {
1015  if (elem->is_edge_on_side(edge,side))
1016  {
1017  if (elem->neighbor_ptr(side) == nullptr)
1018  {
1019  searched_elem = elem->top_parent ();
1020  found_boundary_edge = true;
1021  break;
1022  }
1023  }
1024  }
1025 
1026  if (!found_boundary_edge)
1027  {
1028  // Child element is not on external edge, but it may have internal
1029  // "boundary" IDs. We will walk up the tree, at each level checking that
1030  // the current child is actually on the same edge of the parent that is
1031  // currently being searched for (i.e. that was passed in as "edge").
1032  while (searched_elem->parent() != nullptr)
1033  {
1034  const Elem * parent = searched_elem->parent();
1035  if (parent->is_child_on_edge(parent->which_child_am_i(searched_elem), edge) == false)
1036  return;
1037  searched_elem = parent;
1038  }
1039  }
1040  }
1041 #endif
1042 
1043  // Check each element in the range to see if its edge matches the requested edge.
1044  for (const auto & pr : as_range(_boundary_edge_id.equal_range(searched_elem)))
1045  if (pr.second.first == edge)
1046  vec_to_fill.push_back(pr.second.second);
1047 }
1048 
1049 
1050 
1051 unsigned int BoundaryInfo::n_edge_boundary_ids (const Elem * const elem,
1052  const unsigned short int edge) const
1053 {
1054  std::vector<boundary_id_type> ids;
1055  this->edge_boundary_ids(elem, edge, ids);
1056  return cast_int<unsigned int>(ids.size());
1057 }
1058 
1059 
1060 
1061 #ifdef LIBMESH_ENABLE_DEPRECATED
1062 std::vector<boundary_id_type> BoundaryInfo::raw_edge_boundary_ids (const Elem * const elem,
1063  const unsigned short int edge) const
1064 {
1065  libmesh_deprecated();
1066 
1067  std::vector<boundary_id_type> ids;
1068  this->raw_edge_boundary_ids(elem, edge, ids);
1069  return ids;
1070 }
1071 #endif
1072 
1073 
1074 
1076  const unsigned short int edge,
1077  std::vector<boundary_id_type> & vec_to_fill) const
1078 {
1079  libmesh_assert(elem);
1080 
1081  // Clear out any previous contents
1082  vec_to_fill.clear();
1083 
1084  // Only level-0 elements store BCs.
1085  if (elem->parent())
1086  return;
1087 
1088  // Check each element in the range to see if its edge matches the requested edge.
1089  for (const auto & pr : as_range(_boundary_edge_id.equal_range(elem)))
1090  if (pr.second.first == edge)
1091  vec_to_fill.push_back(pr.second.second);
1092 }
1093 
1094 
1095 
1097  const unsigned short int shellface,
1098  std::vector<boundary_id_type> & vec_to_fill) const
1099 {
1100  libmesh_assert(elem);
1101 
1102  // Shells only have 2 faces
1103  libmesh_assert_less(shellface, 2);
1104 
1105  // Clear out any previous contents
1106  vec_to_fill.clear();
1107 
1108  // Only level-0 elements store BCs. If this is not a level-0
1109  // element get its level-0 parent and infer the BCs.
1110  const Elem * searched_elem = elem;
1111 #ifdef LIBMESH_ENABLE_AMR
1112  if (elem->level() != 0)
1113  {
1114  while (searched_elem->parent() != nullptr)
1115  {
1116  const Elem * parent = searched_elem->parent();
1117  searched_elem = parent;
1118  }
1119  }
1120 #endif
1121 
1122  // Check each element in the range to see if its shellface matches the requested shellface.
1123  for (const auto & pr : as_range(_boundary_shellface_id.equal_range(searched_elem)))
1124  if (pr.second.first == shellface)
1125  vec_to_fill.push_back(pr.second.second);
1126 }
1127 
1128 
1129 
1130 unsigned int BoundaryInfo::n_shellface_boundary_ids (const Elem * const elem,
1131  const unsigned short int shellface) const
1132 {
1133  std::vector<boundary_id_type> ids;
1134  this->shellface_boundary_ids(elem, shellface, ids);
1135  return cast_int<unsigned int>(ids.size());
1136 }
1137 
1138 
1139 
1141  const unsigned short int shellface,
1142  std::vector<boundary_id_type> & vec_to_fill) const
1143 {
1144  libmesh_assert(elem);
1145 
1146  // Shells only have 2 faces
1147  libmesh_assert_less(shellface, 2);
1148 
1149  // Clear out any previous contents
1150  vec_to_fill.clear();
1151 
1152  // Only level-0 elements store BCs.
1153  if (elem->parent())
1154  return;
1155 
1156  // Check each element in the range to see if its shellface matches the requested shellface.
1157  for (const auto & pr : as_range(_boundary_shellface_id.equal_range(elem)))
1158  if (pr.second.first == shellface)
1159  vec_to_fill.push_back(pr.second.second);
1160 }
1161 
1162 
1163 #ifdef LIBMESH_ENABLE_DEPRECATED
1165  const unsigned short int side) const
1166 {
1167  libmesh_deprecated();
1168 
1169  std::vector<boundary_id_type> ids;
1170  this->boundary_ids(elem, side, ids);
1171 
1172  // If the set is empty, return invalid_id
1173  if (ids.empty())
1174  return invalid_id;
1175 
1176  // Otherwise, just return the first id we came across for this
1177  // element on this side.
1178  return *(ids.begin());
1179 }
1180 #endif
1181 
1182 
1183 
1184 bool BoundaryInfo::has_boundary_id(const Elem * const elem,
1185  const unsigned short int side,
1186  const boundary_id_type id) const
1187 {
1188  std::vector<boundary_id_type> ids;
1189  this->boundary_ids(elem, side, ids);
1190  return (std::find(ids.begin(), ids.end(), id) != ids.end());
1191 }
1192 
1193 
1194 
1195 #ifdef LIBMESH_ENABLE_DEPRECATED
1196 std::vector<boundary_id_type> BoundaryInfo::boundary_ids (const Elem * const elem,
1197  const unsigned short int side) const
1198 {
1199  libmesh_deprecated();
1200 
1201  std::vector<boundary_id_type> ids;
1202  this->boundary_ids(elem, side, ids);
1203  return ids;
1204 }
1205 #endif
1206 
1207 
1208 
1209 void BoundaryInfo::boundary_ids (const Elem * const elem,
1210  const unsigned short int side,
1211  std::vector<boundary_id_type> & vec_to_fill) const
1212 {
1213  libmesh_assert(elem);
1214 
1215  // Clear out any previous contents
1216  vec_to_fill.clear();
1217 
1218  // Only level-0 elements store BCs. If this is not a level-0
1219  // element get its level-0 parent and infer the BCs.
1220  const Elem * searched_elem = elem;
1221  if (elem->level() != 0)
1222  {
1223  if (elem->neighbor_ptr(side) == nullptr)
1224  searched_elem = elem->top_parent ();
1225 #ifdef LIBMESH_ENABLE_AMR
1226  else
1227  while (searched_elem->parent() != nullptr)
1228  {
1229  const Elem * parent = searched_elem->parent();
1230  if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
1231  return;
1232  searched_elem = parent;
1233  }
1234 #endif
1235  }
1236 
1237  // Check each element in the range to see if its side matches the requested side.
1238  for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
1239  if (pr.second.first == side)
1240  vec_to_fill.push_back(pr.second.second);
1241 }
1242 
1243 
1244 
1245 
1246 unsigned int BoundaryInfo::n_boundary_ids (const Elem * const elem,
1247  const unsigned short int side) const
1248 {
1249  std::vector<boundary_id_type> ids;
1250  this->boundary_ids(elem, side, ids);
1251  return cast_int<unsigned int>(ids.size());
1252 }
1253 
1254 
1255 
1256 #ifdef LIBMESH_ENABLE_DEPRECATED
1257 std::vector<boundary_id_type> BoundaryInfo::raw_boundary_ids (const Elem * const elem,
1258  const unsigned short int side) const
1259 {
1260  libmesh_deprecated();
1261 
1262  std::vector<boundary_id_type> ids;
1263  this->raw_boundary_ids(elem, side, ids);
1264  return ids;
1265 }
1266 #endif
1267 
1268 
1269 
1270 void BoundaryInfo::raw_boundary_ids (const Elem * const elem,
1271  const unsigned short int side,
1272  std::vector<boundary_id_type> & vec_to_fill) const
1273 {
1274  libmesh_assert(elem);
1275 
1276  // Clear out any previous contents
1277  vec_to_fill.clear();
1278 
1279  // Only level-0 elements store BCs.
1280  if (elem->parent())
1281  return;
1282 
1283  // Check each element in the range to see if its side matches the requested side.
1284  for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
1285  if (pr.second.first == side)
1286  vec_to_fill.push_back(pr.second.second);
1287 }
1288 
1289 
1290 
1291 void BoundaryInfo::copy_boundary_ids (const BoundaryInfo & old_boundary_info,
1292  const Elem * const old_elem,
1293  const Elem * const new_elem)
1294 {
1295  libmesh_assert_equal_to (old_elem->n_sides(), new_elem->n_sides());
1296  libmesh_assert_equal_to (old_elem->n_edges(), new_elem->n_edges());
1297 
1298  std::vector<boundary_id_type> bndry_ids;
1299 
1300  for (auto s : old_elem->side_index_range())
1301  {
1302  old_boundary_info.raw_boundary_ids (old_elem, s, bndry_ids);
1303  this->add_side (new_elem, s, bndry_ids);
1304  }
1305 
1306  for (auto e : old_elem->edge_index_range())
1307  {
1308  old_boundary_info.raw_edge_boundary_ids (old_elem, e, bndry_ids);
1309  this->add_edge (new_elem, e, bndry_ids);
1310  }
1311 
1312  for (unsigned short sf=0; sf != 2; sf++)
1313  {
1314  old_boundary_info.raw_shellface_boundary_ids (old_elem, sf, bndry_ids);
1315  this->add_shellface (new_elem, sf, bndry_ids);
1316  }
1317 }
1318 
1319 
1320 
1321 void BoundaryInfo::remove (const Node * node)
1322 {
1323  libmesh_assert(node);
1324 
1325  // Erase everything associated with node
1326  _boundary_node_id.erase (node);
1327 }
1328 
1329 
1330 
1331 void BoundaryInfo::remove (const Elem * elem)
1332 {
1333  libmesh_assert(elem);
1334 
1335  // Erase everything associated with elem
1336  _boundary_edge_id.erase (elem);
1337  _boundary_side_id.erase (elem);
1338  _boundary_shellface_id.erase (elem);
1339 }
1340 
1341 
1342 
1343 void BoundaryInfo::remove_edge (const Elem * elem,
1344  const unsigned short int edge)
1345 {
1346  libmesh_assert(elem);
1347 
1348  // The user shouldn't be trying to remove only one child's boundary
1349  // id
1350  libmesh_assert_equal_to (elem->level(), 0);
1351 
1352  // Some older compilers don't support erasing from a map with
1353  // const_iterators, so we explicitly use non-const iterators here.
1354  auto e = _boundary_edge_id.equal_range(elem);
1355 
1356  // elem may be there, maybe multiple occurrences
1357  while (e.first != e.second)
1358  {
1359  // if this is true we found the requested edge
1360  // of the element and want to erase the id
1361  if (e.first->second.first == edge)
1362  e.first = _boundary_edge_id.erase(e.first);
1363  else
1364  ++e.first;
1365  }
1366 }
1367 
1368 
1369 
1370 void BoundaryInfo::remove_edge (const Elem * elem,
1371  const unsigned short int edge,
1372  const boundary_id_type id)
1373 {
1374  libmesh_assert(elem);
1375 
1376  // The user shouldn't be trying to remove only one child's boundary
1377  // id
1378  libmesh_assert_equal_to (elem->level(), 0);
1379 
1380  // Some older compilers don't support erasing from a map with
1381  // const_iterators, so we explicitly use non-const iterators here.
1382  auto e = _boundary_edge_id.equal_range(elem);
1383 
1384  // elem may be there, maybe multiple occurrences
1385  while (e.first != e.second)
1386  {
1387  // if this is true we found the requested edge
1388  // of the element and want to erase the requested id
1389  if (e.first->second.first == edge && e.first->second.second == id)
1390  e.first = _boundary_edge_id.erase(e.first);
1391  else
1392  ++e.first;
1393  }
1394 }
1395 
1396 
1398  const unsigned short int shellface)
1399 {
1400  libmesh_assert(elem);
1401 
1402  // The user shouldn't be trying to remove only one child's boundary
1403  // id
1404  libmesh_assert_equal_to (elem->level(), 0);
1405 
1406  // Shells only have 2 faces
1407  libmesh_assert_less(shellface, 2);
1408 
1409  // Some older compilers don't support erasing from a map with
1410  // const_iterators, so we explicitly use non-const iterators here.
1411  auto e = _boundary_shellface_id.equal_range(elem);
1412 
1413  // elem may be there, maybe multiple occurrences
1414  while (e.first != e.second)
1415  {
1416  // if this is true we found the requested shellface
1417  // of the element and want to erase the id
1418  if (e.first->second.first == shellface)
1419  e.first = _boundary_shellface_id.erase(e.first);
1420  else
1421  ++e.first;
1422  }
1423 }
1424 
1425 
1426 
1428  const unsigned short int shellface,
1429  const boundary_id_type id)
1430 {
1431  libmesh_assert(elem);
1432 
1433  // The user shouldn't be trying to remove only one child's boundary
1434  // id
1435  libmesh_assert_equal_to (elem->level(), 0);
1436 
1437  // Shells only have 2 faces
1438  libmesh_assert_less(shellface, 2);
1439 
1440  // Some older compilers don't support erasing from a map with
1441  // const_iterators, so we explicitly use non-const iterators here.
1442  auto e = _boundary_shellface_id.equal_range(elem);
1443 
1444  // elem may be there, maybe multiple occurrences
1445  while (e.first != e.second)
1446  {
1447  // if this is true we found the requested shellface
1448  // of the element and want to erase the requested id
1449  if (e.first->second.first == shellface &&
1450  e.first->second.second == id)
1451  {
1452  // (postfix++ - increment the iterator before it's invalid)
1453  e.first = _boundary_shellface_id.erase(e.first);
1454  }
1455  else
1456  ++e.first;
1457  }
1458 }
1459 
1460 void BoundaryInfo::remove_side (const Elem * elem,
1461  const unsigned short int side)
1462 {
1463  libmesh_assert(elem);
1464 
1465  // The user shouldn't be trying to remove only one child's boundary
1466  // id
1467  libmesh_assert_equal_to (elem->level(), 0);
1468 
1469  // Some older compilers don't support erasing from a map with
1470  // const_iterators, so we explicitly use non-const iterators here.
1471  auto e = _boundary_side_id.equal_range(elem);
1472 
1473  // elem may be there, maybe multiple occurrences
1474  while (e.first != e.second)
1475  {
1476  // if this is true we found the requested side
1477  // of the element and want to erase the id
1478  if (e.first->second.first == side)
1479  e.first = _boundary_side_id.erase(e.first);
1480  else
1481  ++e.first;
1482  }
1483 }
1484 
1485 
1486 
1487 void BoundaryInfo::remove_side (const Elem * elem,
1488  const unsigned short int side,
1489  const boundary_id_type id)
1490 {
1491  libmesh_assert(elem);
1492 
1493  // Some older compilers don't support erasing from a map with
1494  // const_iterators, so we explicitly use non-const iterators here.
1495  auto e = _boundary_side_id.equal_range(elem);
1496 
1497  // elem may be there, maybe multiple occurrences
1498  while (e.first != e.second)
1499  {
1500  // if this is true we found the requested side
1501  // of the element and want to erase the requested id
1502  if (e.first->second.first == side && e.first->second.second == id)
1503  e.first = _boundary_side_id.erase(e.first);
1504  else
1505  ++e.first;
1506  }
1507 }
1508 
1509 
1510 
1512 {
1513  // Erase id from ids containers
1514  _boundary_ids.erase(id);
1515  _side_boundary_ids.erase(id);
1516  _edge_boundary_ids.erase(id);
1517  _shellface_boundary_ids.erase(id);
1518  _node_boundary_ids.erase(id);
1519  _ss_id_to_name.erase(id);
1520  _ns_id_to_name.erase(id);
1521 
1522  // Erase pointers to geometric entities with this id.
1523  for (auto it = _boundary_node_id.begin(); it != _boundary_node_id.end(); /*below*/)
1524  {
1525  if (it->second == id)
1526  it = _boundary_node_id.erase(it);
1527  else
1528  ++it;
1529  }
1530 
1531  for (auto it = _boundary_edge_id.begin(); it != _boundary_edge_id.end(); /*below*/)
1532  {
1533  if (it->second.second == id)
1534  it = _boundary_edge_id.erase(it);
1535  else
1536  ++it;
1537  }
1538 
1539  for (auto it = _boundary_shellface_id.begin(); it != _boundary_shellface_id.end(); /*below*/)
1540  {
1541  if (it->second.second == id)
1542  it = _boundary_shellface_id.erase(it);
1543  else
1544  ++it;
1545  }
1546 
1547  for (auto it = _boundary_side_id.begin(); it != _boundary_side_id.end(); /*below*/)
1548  {
1549  if (it->second.second == id)
1550  it = _boundary_side_id.erase(it);
1551  else
1552  ++it;
1553  }
1554 }
1555 
1556 
1557 
1558 unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
1559  const boundary_id_type boundary_id_in) const
1560 {
1561  const Elem * searched_elem = elem;
1562  if (elem->level() != 0)
1563  searched_elem = elem->top_parent();
1564 
1565  // elem may have zero or multiple occurrences
1566  for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
1567  {
1568  // if this is true we found the requested boundary_id
1569  // of the element and want to return the side
1570  if (pr.second.second == boundary_id_in)
1571  {
1572  unsigned int side = pr.second.first;
1573 
1574  // If we're on this external boundary then we share this
1575  // external boundary id
1576  if (elem->neighbor_ptr(side) == nullptr)
1577  return side;
1578 
1579  // If we're on an internal boundary then we need to be sure
1580  // it's the same internal boundary as our top_parent
1581  const Elem * p = elem;
1582 
1583 #ifdef LIBMESH_ENABLE_AMR
1584 
1585  while (p != nullptr)
1586  {
1587  const Elem * parent = p->parent();
1588  if (!parent->is_child_on_side(parent->which_child_am_i(p), side))
1589  break;
1590  p = parent;
1591  }
1592 #endif
1593  // We're on that side of our top_parent; return it
1594  if (!p)
1595  return side;
1596  }
1597  }
1598 
1599  // if we get here, we found elem in the data structure but not
1600  // the requested boundary id, so return the default value
1601  return libMesh::invalid_uint;
1602 }
1603 
1604 void
1605 BoundaryInfo::build_node_boundary_ids(std::vector<boundary_id_type> & b_ids) const
1606 {
1607  b_ids.clear();
1608 
1609  for (const auto & pr : _boundary_node_id)
1610  {
1611  boundary_id_type id = pr.second;
1612 
1613  if (std::find(b_ids.begin(),b_ids.end(),id) == b_ids.end())
1614  b_ids.push_back(id);
1615  }
1616 }
1617 
1618 void
1619 BoundaryInfo::build_side_boundary_ids(std::vector<boundary_id_type> & b_ids) const
1620 {
1621  b_ids.clear();
1622 
1623  for (const auto & pr : _boundary_side_id)
1624  {
1625  boundary_id_type id = pr.second.second;
1626 
1627  if (std::find(b_ids.begin(),b_ids.end(),id) == b_ids.end())
1628  b_ids.push_back(id);
1629  }
1630 }
1631 
1632 void
1633 BoundaryInfo::build_shellface_boundary_ids(std::vector<boundary_id_type> & b_ids) const
1634 {
1635  b_ids.clear();
1636 
1637  for (const auto & pr :_boundary_shellface_id)
1638  {
1639  boundary_id_type id = pr.second.second;
1640 
1641  if (std::find(b_ids.begin(),b_ids.end(),id) == b_ids.end())
1642  b_ids.push_back(id);
1643  }
1644 }
1645 
1647 {
1648  // in serial we know the number of bcs from the
1649  // size of the container
1650  if (_mesh.is_serial())
1651  return _boundary_side_id.size();
1652 
1653  // in parallel we need to sum the number of local bcs
1654  parallel_object_only();
1655 
1656  std::size_t nbcs=0;
1657 
1658  for (const auto & pr : _boundary_side_id)
1659  if (pr.first->processor_id() == this->processor_id())
1660  nbcs++;
1661 
1662  this->comm().sum (nbcs);
1663 
1664  return nbcs;
1665 }
1666 
1667 std::size_t BoundaryInfo::n_edge_conds () const
1668 {
1669  // in serial we know the number of nodesets from the
1670  // size of the container
1671  if (_mesh.is_serial())
1672  return _boundary_edge_id.size();
1673 
1674  // in parallel we need to sum the number of local nodesets
1675  parallel_object_only();
1676 
1677  std::size_t n_edge_bcs=0;
1678 
1679  for (const auto & pr : _boundary_edge_id)
1680  if (pr.first->processor_id() == this->processor_id())
1681  n_edge_bcs++;
1682 
1683  this->comm().sum (n_edge_bcs);
1684 
1685  return n_edge_bcs;
1686 }
1687 
1688 
1690 {
1691  // in serial we know the number of nodesets from the
1692  // size of the container
1693  if (_mesh.is_serial())
1694  return _boundary_shellface_id.size();
1695 
1696  // in parallel we need to sum the number of local nodesets
1697  parallel_object_only();
1698 
1699  std::size_t n_shellface_bcs=0;
1700 
1701  for (const auto & pr : _boundary_shellface_id)
1702  if (pr.first->processor_id() == this->processor_id())
1703  n_shellface_bcs++;
1704 
1705  this->comm().sum (n_shellface_bcs);
1706 
1707  return n_shellface_bcs;
1708 }
1709 
1710 
1711 std::size_t BoundaryInfo::n_nodeset_conds () const
1712 {
1713  // in serial we know the number of nodesets from the
1714  // size of the container
1715  if (_mesh.is_serial())
1716  return _boundary_node_id.size();
1717 
1718  // in parallel we need to sum the number of local nodesets
1719  parallel_object_only();
1720 
1721  std::size_t n_nodesets=0;
1722 
1723  for (const auto & pr : _boundary_node_id)
1724  if (pr.first->processor_id() == this->processor_id())
1725  n_nodesets++;
1726 
1727  this->comm().sum (n_nodesets);
1728 
1729  return n_nodesets;
1730 }
1731 
1732 
1733 
1734 #ifdef LIBMESH_ENABLE_DEPRECATED
1735 void BoundaryInfo::build_node_list (std::vector<dof_id_type> & nl,
1736  std::vector<boundary_id_type> & il) const
1737 {
1738  libmesh_deprecated();
1739 
1740  // Call the non-deprecated version of this function.
1741  auto bc_tuples = this->build_node_list();
1742 
1743  // Clear the input vectors, just in case they were used for
1744  // something else recently...
1745  nl.clear();
1746  il.clear();
1747 
1748  // Reserve the size, then use push_back
1749  nl.reserve (bc_tuples.size());
1750  il.reserve (bc_tuples.size());
1751 
1752  for (const auto & t : bc_tuples)
1753  {
1754  nl.push_back(std::get<0>(t));
1755  il.push_back(std::get<1>(t));
1756  }
1757 }
1758 #endif
1759 
1760 
1761 std::vector<std::tuple<dof_id_type, boundary_id_type>>
1763 {
1764  std::vector<std::tuple<dof_id_type, boundary_id_type>> bc_tuples;
1765  bc_tuples.reserve(_boundary_node_id.size());
1766 
1767  for (const auto & pr : _boundary_node_id)
1768  bc_tuples.emplace_back(pr.first->id(), pr.second);
1769 
1770  // This list is currently in memory address (arbitrary) order, so
1771  // sort to make it consistent on all procs.
1772  std::sort(bc_tuples.begin(), bc_tuples.end());
1773 
1774  return bc_tuples;
1775 }
1776 
1777 
1778 void
1780 {
1781  // If we're on a distributed mesh, even the owner of a node is not
1782  // guaranteed to be able to properly assign its new boundary id(s)!
1783  // Nodal neighbors are not always ghosted, and a nodal neighbor
1784  // might have a boundary side.
1785  const bool mesh_is_serial = _mesh.is_serial();
1786 
1787  typedef std::set<std::pair<dof_id_type, boundary_id_type>> set_type;
1788 
1789  const processor_id_type n_proc = this->n_processors();
1790  const processor_id_type my_proc_id = this->processor_id();
1791  std::vector<set_type> nodes_to_push(n_proc);
1792 
1793  // Pull objects out of the loop to reduce heap operations
1794  std::unique_ptr<const Elem> side;
1795 
1796  // Loop over the side list
1797  for (const auto & pr : _boundary_side_id)
1798  {
1799  // Don't add remote sides
1800  if (pr.first->is_remote())
1801  continue;
1802 
1803  // Need to loop over the sides of any possible children
1804  std::vector<const Elem *> family;
1805 #ifdef LIBMESH_ENABLE_AMR
1806  pr.first->active_family_tree_by_side (family, pr.second.first);
1807 #else
1808  family.push_back(pr.first);
1809 #endif
1810 
1811  for (const auto & cur_elem : family)
1812  {
1813  cur_elem->build_side_ptr(side, pr.second.first);
1814 
1815  // Add each node node on the side with the side's boundary id
1816  for (auto i : side->node_index_range())
1817  {
1818  const boundary_id_type bcid = pr.second.second;
1819  this->add_node(side->node_ptr(i), bcid);
1820  if (!mesh_is_serial)
1821  {
1822  const processor_id_type proc_id =
1823  side->node_ptr(i)->processor_id();
1824  if (proc_id != my_proc_id)
1825  nodes_to_push[proc_id].insert
1826  (std::make_pair(side->node_id(i), bcid));
1827  }
1828  }
1829  }
1830  }
1831 
1832  // If we're on a serial mesh then we're done.
1833  if (mesh_is_serial)
1834  return;
1835 
1836  // Otherwise we need to push ghost node bcids to their owners, then
1837  // pull ghost node bcids from their owners.
1839  node_pushes_tag = this->comm().get_unique_tag(31337),
1840  node_pulls_tag = this->comm().get_unique_tag(31338),
1841  node_responses_tag = this->comm().get_unique_tag(31339);
1842 
1843  std::vector<Parallel::Request> node_push_requests(n_proc-1);
1844 
1845  for (processor_id_type p = 0; p != n_proc; ++p)
1846  {
1847  if (p == my_proc_id)
1848  continue;
1849 
1851  node_push_requests[p - (p > my_proc_id)];
1852 
1853  this->comm().send
1854  (p, nodes_to_push[p], request, node_pushes_tag);
1855  }
1856 
1857  for (processor_id_type p = 1; p != n_proc; ++p)
1858  {
1859  set_type received_nodes;
1860 
1861  this->comm().receive
1862  (Parallel::any_source, received_nodes, node_pushes_tag);
1863 
1864  for (const auto & pr : received_nodes)
1865  this->add_node(_mesh.node_ptr(pr.first), pr.second);
1866  }
1867 
1868  // At this point we should know all the BCs for our own nodes; now
1869  // we need BCs for ghost nodes.
1870  //
1871  // FIXME - parallel_ghost_sync.h doesn't work here because it
1872  // assumes a fixed size datum on each node.
1873  std::vector<std::vector<dof_id_type>> node_ids_requested(n_proc);
1874 
1875  // Determine what nodes we need to request
1876  for (const auto & node : _mesh.node_ptr_range())
1877  {
1878  const processor_id_type pid = node->processor_id();
1879  if (pid != my_proc_id)
1880  node_ids_requested[pid].push_back(node->id());
1881  }
1882 
1883  typedef std::vector<std::pair<dof_id_type, boundary_id_type>> vec_type;
1884 
1885  std::vector<Parallel::Request>
1886  node_pull_requests(n_proc-1),
1887  node_response_requests(n_proc-1);
1888 
1889  // Make all requests
1890  for (processor_id_type p = 0; p != n_proc; ++p)
1891  {
1892  if (p == my_proc_id)
1893  continue;
1894 
1896  node_pull_requests[p - (p > my_proc_id)];
1897 
1898  this->comm().send
1899  (p, node_ids_requested[p], request, node_pulls_tag);
1900  }
1901 
1902  // Process all incoming requests
1903  std::vector<vec_type> responses(n_proc-1);
1904 
1905  for (processor_id_type p = 1; p != n_proc; ++p)
1906  {
1907  std::vector<dof_id_type> requested_nodes;
1908 
1910  status(this->comm().probe (Parallel::any_source, node_pulls_tag));
1911  const processor_id_type
1912  source_pid = cast_int<processor_id_type>(status.source());
1913 
1914  this->comm().receive
1915  (source_pid, requested_nodes, node_pulls_tag);
1916 
1918  node_response_requests[p-1];
1919 
1920  std::vector<boundary_id_type> bcids;
1921 
1922  for (const auto & id : requested_nodes)
1923  {
1924  this->boundary_ids(_mesh.node_ptr(id), bcids);
1925 
1926  for (const auto & b : bcids)
1927  responses[p-1].push_back(std::make_pair(id, b));
1928  }
1929 
1930  this->comm().send
1931  (source_pid, responses[p-1], request, node_responses_tag);
1932  }
1933 
1934  // Process all incoming responses
1935  for (processor_id_type p = 1; p != n_proc; ++p)
1936  {
1938  status(this->comm().probe (Parallel::any_source, node_responses_tag));
1939  const processor_id_type
1940  source_pid = cast_int<processor_id_type>(status.source());
1941 
1942  vec_type response;
1943 
1944  this->comm().receive
1945  (source_pid, response, node_responses_tag);
1946 
1947  for (const auto & pr : response)
1948  this->add_node(_mesh.node_ptr(pr.first), pr.second);
1949  }
1950 
1951  Parallel::wait (node_push_requests);
1952  Parallel::wait (node_pull_requests);
1953  Parallel::wait (node_response_requests);
1954 }
1955 
1956 
1957 
1958 
1960 {
1961  // Check for early return
1962  if (_boundary_node_id.empty())
1963  {
1964  libMesh::out << "No boundary node IDs have been added: cannot build side list!" << std::endl;
1965  return;
1966  }
1967 
1968  // Pull objects out of the loop to reduce heap operations
1969  std::unique_ptr<const Elem> side_elem;
1970 
1971  for (const auto & elem : _mesh.active_element_ptr_range())
1972  for (auto side : elem->side_index_range())
1973  {
1974  elem->build_side_ptr(side_elem, side);
1975 
1976  // map from nodeset_id to count for that ID
1977  std::map<boundary_id_type, unsigned> nodesets_node_count;
1978 
1979  // For each nodeset that this node is a member of, increment the associated
1980  // nodeset ID count
1981  for (const auto & node : side_elem->node_ref_range())
1982  for (const auto & pr : as_range(_boundary_node_id.equal_range(&node)))
1983  nodesets_node_count[pr.second]++;
1984 
1985  // Now check to see what nodeset_counts have the correct
1986  // number of nodes in them. For any that do, add this side to
1987  // the sideset, making sure the sideset inherits the
1988  // nodeset's name, if there is one.
1989  for (const auto & pr : nodesets_node_count)
1990  if (pr.second == side_elem->n_nodes())
1991  {
1992  add_side(elem, side, pr.first);
1993 
1994  // Let the sideset inherit any non-empty name from the nodeset
1995  std::string & nset_name = nodeset_name(pr.first);
1996 
1997  if (nset_name != "")
1998  sideset_name(pr.first) = nset_name;
1999  }
2000  } // end for side
2001 }
2002 
2003 
2004 
2005 
2006 #ifdef LIBMESH_ENABLE_DEPRECATED
2007 void BoundaryInfo::build_side_list (std::vector<dof_id_type> & el,
2008  std::vector<unsigned short int> & sl,
2009  std::vector<boundary_id_type> & il) const
2010 {
2011  libmesh_deprecated();
2012 
2013  // Call the non-deprecated version of this function.
2014  auto bc_tuples = this->build_side_list();
2015 
2016  // Clear the input vectors, just in case they were used for
2017  // something else recently...
2018  el.clear();
2019  sl.clear();
2020  il.clear();
2021 
2022  // Reserve the size, then use push_back
2023  el.reserve (bc_tuples.size());
2024  sl.reserve (bc_tuples.size());
2025  il.reserve (bc_tuples.size());
2026 
2027  for (const auto & t : bc_tuples)
2028  {
2029  el.push_back(std::get<0>(t));
2030  sl.push_back(std::get<1>(t));
2031  il.push_back(std::get<2>(t));
2032  }
2033 }
2034 #endif
2035 
2036 
2037 std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
2039 {
2040  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> bc_triples;
2041  bc_triples.reserve(_boundary_side_id.size());
2042 
2043  for (const auto & pr : _boundary_side_id)
2044  bc_triples.emplace_back(pr.first->id(), pr.second.first, pr.second.second);
2045 
2046  // bc_triples is currently in whatever order the Elem pointers in
2047  // the _boundary_side_id multimap are in, and in particular might be
2048  // in different orders on different processors. To avoid this
2049  // inconsistency, we'll sort using the default operator< for tuples.
2050  std::sort(bc_triples.begin(), bc_triples.end());
2051 
2052  return bc_triples;
2053 }
2054 
2055 
2056 
2057 #ifdef LIBMESH_ENABLE_DEPRECATED
2058 void BoundaryInfo::build_active_side_list (std::vector<dof_id_type> & el,
2059  std::vector<unsigned short int> & sl,
2060  std::vector<boundary_id_type> & il) const
2061 {
2062  libmesh_deprecated();
2063 
2064  // Call the non-deprecated version of this function.
2065  auto bc_tuples = this->build_active_side_list();
2066 
2067  // Clear the input vectors, just in case they were used for
2068  // something else recently...
2069  el.clear();
2070  sl.clear();
2071  il.clear();
2072 
2073  // Reserve the size, then use push_back
2074  el.reserve (bc_tuples.size());
2075  sl.reserve (bc_tuples.size());
2076  il.reserve (bc_tuples.size());
2077 
2078  for (const auto & t : bc_tuples)
2079  {
2080  el.push_back(std::get<0>(t));
2081  sl.push_back(std::get<1>(t));
2082  il.push_back(std::get<2>(t));
2083  }
2084 }
2085 #endif
2086 
2087 
2088 std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
2090 {
2091  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> bc_triples;
2092  bc_triples.reserve(_boundary_side_id.size());
2093 
2094  for (const auto & pr : _boundary_side_id)
2095  {
2096  // Don't add remote sides
2097  if (pr.first->is_remote())
2098  continue;
2099 
2100  // Loop over the sides of possible children
2101  std::vector<const Elem *> family;
2102 #ifdef LIBMESH_ENABLE_AMR
2103  pr.first->active_family_tree_by_side(family, pr.second.first);
2104 #else
2105  family.push_back(pr.first);
2106 #endif
2107 
2108  // Populate the list items
2109  for (const auto & elem : family)
2110  bc_triples.emplace_back(elem->id(), pr.second.first, pr.second.second);
2111  }
2112 
2113  // This list is currently in memory address (arbitrary) order, so
2114  // sort to make it consistent on all procs.
2115  std::sort(bc_triples.begin(), bc_triples.end());
2116 
2117  return bc_triples;
2118 }
2119 
2120 
2121 #ifdef LIBMESH_ENABLE_DEPRECATED
2122 void BoundaryInfo::build_edge_list (std::vector<dof_id_type> & el,
2123  std::vector<unsigned short int> & sl,
2124  std::vector<boundary_id_type> & il) const
2125 {
2126  libmesh_deprecated();
2127 
2128  // Call the non-deprecated version of this function.
2129  auto bc_tuples = this->build_edge_list();
2130 
2131  // Clear the input vectors, just in case they were used for
2132  // something else recently...
2133  el.clear();
2134  sl.clear();
2135  il.clear();
2136 
2137  // Reserve the size, then use push_back
2138  el.reserve (bc_tuples.size());
2139  sl.reserve (bc_tuples.size());
2140  il.reserve (bc_tuples.size());
2141 
2142  for (const auto & t : bc_tuples)
2143  {
2144  el.push_back(std::get<0>(t));
2145  sl.push_back(std::get<1>(t));
2146  il.push_back(std::get<2>(t));
2147  }
2148 }
2149 #endif
2150 
2151 
2152 std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
2154 {
2155  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> bc_triples;
2156  bc_triples.reserve(_boundary_edge_id.size());
2157 
2158  for (const auto & pr : _boundary_edge_id)
2159  bc_triples.emplace_back(pr.first->id(), pr.second.first, pr.second.second);
2160 
2161  // This list is currently in memory address (arbitrary) order, so
2162  // sort to make it consistent on all procs.
2163  std::sort(bc_triples.begin(), bc_triples.end());
2164 
2165  return bc_triples;
2166 }
2167 
2168 
2169 #ifdef LIBMESH_ENABLE_DEPRECATED
2170 void BoundaryInfo::build_shellface_list (std::vector<dof_id_type> & el,
2171  std::vector<unsigned short int> & sl,
2172  std::vector<boundary_id_type> & il) const
2173 {
2174  libmesh_deprecated();
2175 
2176  // Call the non-deprecated version of this function.
2177  auto bc_tuples = this->build_shellface_list();
2178 
2179  // Clear the input vectors, just in case they were used for
2180  // something else recently...
2181  el.clear();
2182  sl.clear();
2183  il.clear();
2184 
2185  // Reserve the size, then use push_back
2186  el.reserve (bc_tuples.size());
2187  sl.reserve (bc_tuples.size());
2188  il.reserve (bc_tuples.size());
2189 
2190  for (const auto & t : bc_tuples)
2191  {
2192  el.push_back(std::get<0>(t));
2193  sl.push_back(std::get<1>(t));
2194  il.push_back(std::get<2>(t));
2195  }
2196 }
2197 #endif
2198 
2199 
2200 std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>>
2202 {
2203  std::vector<std::tuple<dof_id_type, unsigned short int, boundary_id_type>> bc_triples;
2204  bc_triples.reserve(_boundary_shellface_id.size());
2205 
2206  for (const auto & pr : _boundary_shellface_id)
2207  bc_triples.emplace_back(pr.first->id(), pr.second.first, pr.second.second);
2208 
2209  // This list is currently in memory address (arbitrary) order, so
2210  // sort to make it consistent on all procs.
2211  std::sort(bc_triples.begin(), bc_triples.end());
2212 
2213  return bc_triples;
2214 }
2215 
2216 
2217 void BoundaryInfo::print_info(std::ostream & out_stream) const
2218 {
2219  // Print out the nodal BCs
2220  if (!_boundary_node_id.empty())
2221  {
2222  out_stream << "Nodal Boundary conditions:" << std::endl
2223  << "--------------------------" << std::endl
2224  << " (Node No., ID) " << std::endl;
2225 
2226  for (const auto & pr : _boundary_node_id)
2227  out_stream << " (" << pr.first->id()
2228  << ", " << pr.second
2229  << ")" << std::endl;
2230  }
2231 
2232  // Print out the element edge BCs
2233  if (!_boundary_edge_id.empty())
2234  {
2235  out_stream << std::endl
2236  << "Edge Boundary conditions:" << std::endl
2237  << "-------------------------" << std::endl
2238  << " (Elem No., Edge No., ID) " << std::endl;
2239 
2240  for (const auto & pr : _boundary_edge_id)
2241  out_stream << " (" << pr.first->id()
2242  << ", " << pr.second.first
2243  << ", " << pr.second.second
2244  << ")" << std::endl;
2245  }
2246 
2247  // Print out the element shellface BCs
2248  if (!_boundary_shellface_id.empty())
2249  {
2250  out_stream << std::endl
2251  << "Shell-face Boundary conditions:" << std::endl
2252  << "-------------------------" << std::endl
2253  << " (Elem No., Shell-face No., ID) " << std::endl;
2254 
2255  for (const auto & pr : _boundary_shellface_id)
2256  out_stream << " (" << pr.first->id()
2257  << ", " << pr.second.first
2258  << ", " << pr.second.second
2259  << ")" << std::endl;
2260  }
2261 
2262  // Print out the element side BCs
2263  if (!_boundary_side_id.empty())
2264  {
2265  out_stream << std::endl
2266  << "Side Boundary conditions:" << std::endl
2267  << "-------------------------" << std::endl
2268  << " (Elem No., Side No., ID) " << std::endl;
2269 
2270  for (const auto & pr : _boundary_side_id)
2271  out_stream << " (" << pr.first->id()
2272  << ", " << pr.second.first
2273  << ", " << pr.second.second
2274  << ")" << std::endl;
2275  }
2276 }
2277 
2278 
2279 
2280 void BoundaryInfo::print_summary(std::ostream & out_stream) const
2281 {
2282  // Print out the nodal BCs
2283  if (!_boundary_node_id.empty())
2284  {
2285  out_stream << "Nodal Boundary conditions:" << std::endl
2286  << "--------------------------" << std::endl
2287  << " (ID, number of nodes) " << std::endl;
2288 
2289  std::map<boundary_id_type, std::size_t> ID_counts;
2290 
2291  for (const auto & pr : _boundary_node_id)
2292  ID_counts[pr.second]++;
2293 
2294  for (const auto & pr : ID_counts)
2295  out_stream << " (" << pr.first
2296  << ", " << pr.second
2297  << ")" << std::endl;
2298  }
2299 
2300  // Print out the element edge BCs
2301  if (!_boundary_edge_id.empty())
2302  {
2303  out_stream << std::endl
2304  << "Edge Boundary conditions:" << std::endl
2305  << "-------------------------" << std::endl
2306  << " (ID, number of edges) " << std::endl;
2307 
2308  std::map<boundary_id_type, std::size_t> ID_counts;
2309 
2310  for (const auto & pr : _boundary_edge_id)
2311  ID_counts[pr.second.second]++;
2312 
2313  for (const auto & pr : ID_counts)
2314  out_stream << " (" << pr.first
2315  << ", " << pr.second
2316  << ")" << std::endl;
2317  }
2318 
2319 
2320  // Print out the element edge BCs
2321  if (!_boundary_shellface_id.empty())
2322  {
2323  out_stream << std::endl
2324  << "Shell-face Boundary conditions:" << std::endl
2325  << "-------------------------" << std::endl
2326  << " (ID, number of shellfaces) " << std::endl;
2327 
2328  std::map<boundary_id_type, std::size_t> ID_counts;
2329 
2330  for (const auto & pr : _boundary_shellface_id)
2331  ID_counts[pr.second.second]++;
2332 
2333  for (const auto & pr : ID_counts)
2334  out_stream << " (" << pr.first
2335  << ", " << pr.second
2336  << ")" << std::endl;
2337  }
2338 
2339  // Print out the element side BCs
2340  if (!_boundary_side_id.empty())
2341  {
2342  out_stream << std::endl
2343  << "Side Boundary conditions:" << std::endl
2344  << "-------------------------" << std::endl
2345  << " (ID, number of sides) " << std::endl;
2346 
2347  std::map<boundary_id_type, std::size_t> ID_counts;
2348 
2349  for (const auto & pr : _boundary_side_id)
2350  ID_counts[pr.second.second]++;
2351 
2352  for (const auto & pr : ID_counts)
2353  out_stream << " (" << pr.first
2354  << ", " << pr.second
2355  << ")" << std::endl;
2356  }
2357 }
2358 
2359 
2361 {
2362  static const std::string empty_string;
2363  std::map<boundary_id_type, std::string>::const_iterator it =
2364  _ss_id_to_name.find(id);
2365  if (it == _ss_id_to_name.end())
2366  return empty_string;
2367  else
2368  return it->second;
2369 }
2370 
2371 
2373 {
2374  return _ss_id_to_name[id];
2375 }
2376 
2378 {
2379  static const std::string empty_string;
2380  std::map<boundary_id_type, std::string>::const_iterator it =
2381  _ns_id_to_name.find(id);
2382  if (it == _ns_id_to_name.end())
2383  return empty_string;
2384  else
2385  return it->second;
2386 }
2387 
2389 {
2390  return _ns_id_to_name[id];
2391 }
2392 
2394 {
2395  // Search sidesets
2396  for (const auto & pr : _ss_id_to_name)
2397  if (pr.second == name)
2398  return pr.first;
2399 
2400  // Search nodesets
2401  for (const auto & pr : _ns_id_to_name)
2402  if (pr.second == name)
2403  return pr.first;
2404 
2405  // If we made it here without returning, we don't have a sideset or
2406  // nodeset by the requested name, so return invalid_id
2407  return invalid_id;
2408 }
2409 
2410 
2411 
2412 void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & requested_boundary_ids,
2413  dof_id_type first_free_node_id,
2414  std::map<dof_id_type, dof_id_type> * node_id_map,
2415  dof_id_type first_free_elem_id,
2416  std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> * side_id_map,
2417  const std::set<subdomain_id_type> & subdomains_relative_to)
2418 {
2419  // We'll do the same modulus trick that DistributedMesh uses to avoid
2420  // id conflicts between different processors
2421  dof_id_type
2422  next_node_id = first_free_node_id + this->processor_id(),
2423  next_elem_id = first_free_elem_id + this->processor_id();
2424 
2425  // Pull objects out of the loop to reduce heap operations
2426  std::unique_ptr<const Elem> side;
2427 
2428  // We'll pass through the mesh once first to build
2429  // the maps and count boundary nodes and elements.
2430  // To find local boundary nodes, we have to examine all elements
2431  // here rather than just local elements, because it's possible to
2432  // have a local boundary node that's not on a local boundary
2433  // element, e.g. at the tip of a triangle.
2434 
2435  // We'll loop through two different ranges here: first all elements,
2436  // looking for local nodes, and second through unpartitioned
2437  // elements, looking for all remaining nodes.
2439  bool hit_end_el = false;
2440  const MeshBase::const_element_iterator end_unpartitioned_el =
2442 
2444  !hit_end_el || (el != end_unpartitioned_el); ++el)
2445  {
2446  if ((el == end_el) && !hit_end_el)
2447  {
2448  // Note that we're done with local nodes and just looking
2449  // for remaining unpartitioned nodes
2450  hit_end_el = true;
2451 
2452  // Join up the local results from other processors
2453  if (side_id_map)
2454  this->comm().set_union(*side_id_map);
2455  if (node_id_map)
2456  this->comm().set_union(*node_id_map);
2457 
2458  // Finally we'll pass through any unpartitioned elements to add them
2459  // to the maps and counts.
2460  next_node_id = first_free_node_id + this->n_processors();
2461  next_elem_id = first_free_elem_id + this->n_processors();
2462 
2464  if (el == end_unpartitioned_el)
2465  break;
2466  }
2467 
2468  const Elem * elem = *el;
2469 
2470  // If the subdomains_relative_to container has the
2471  // invalid_subdomain_id, we fall back on the "old" behavior of
2472  // adding sides regardless of this Elem's subdomain. Otherwise,
2473  // if the subdomains_relative_to container doesn't contain the
2474  // current Elem's subdomain_id(), we won't add any sides from
2475  // it.
2476  if (!subdomains_relative_to.count(Elem::invalid_subdomain_id) &&
2477  !subdomains_relative_to.count(elem->subdomain_id()))
2478  continue;
2479 
2480  // Get the top-level parent for this element. This is used for
2481  // searching for boundary sides on this element.
2482  const Elem * top_parent = elem->top_parent();
2483 
2484  // Find all the boundary side ids for this Elem.
2485  auto bounds = _boundary_side_id.equal_range(top_parent);
2486 
2487  for (auto s : elem->side_index_range())
2488  {
2489  bool add_this_side = false;
2490  boundary_id_type this_bcid = invalid_id;
2491 
2492  for (const auto & pr : as_range(bounds))
2493  {
2494  this_bcid = pr.second.second;
2495 
2496  // if this side is flagged with a boundary condition
2497  // and the user wants this id
2498  if ((pr.second.first == s) &&
2499  (requested_boundary_ids.count(this_bcid)))
2500  {
2501  add_this_side = true;
2502  break;
2503  }
2504  }
2505 
2506  // We may still want to add this side if the user called
2507  // sync() with no requested_boundary_ids. This corresponds
2508  // to the "old" style of calling sync() in which the entire
2509  // boundary was copied to the BoundaryMesh, and handles the
2510  // case where elements on the geometric boundary are not in
2511  // any sidesets.
2512  if (bounds.first == bounds.second &&
2513  requested_boundary_ids.count(invalid_id) &&
2514  elem->neighbor_ptr(s) == nullptr)
2515  add_this_side = true;
2516 
2517  if (add_this_side)
2518  {
2519  // We only assign ids for our own and for
2520  // unpartitioned elements
2521  if (side_id_map &&
2522  ((elem->processor_id() == this->processor_id()) ||
2523  (elem->processor_id() ==
2525  {
2526  std::pair<dof_id_type, unsigned char> side_pair(elem->id(), s);
2527  libmesh_assert (!side_id_map->count(side_pair));
2528  (*side_id_map)[side_pair] = next_elem_id;
2529  next_elem_id += this->n_processors() + 1;
2530  }
2531 
2532  elem->build_side_ptr(side, s);
2533  for (auto n : side->node_index_range())
2534  {
2535  const Node & node = side->node_ref(n);
2536 
2537  // In parallel we don't know enough to number
2538  // others' nodes ourselves.
2539  if (!hit_end_el &&
2540  (node.processor_id() != this->processor_id()))
2541  continue;
2542 
2543  dof_id_type node_id = node.id();
2544  if (node_id_map && !node_id_map->count(node_id))
2545  {
2546  (*node_id_map)[node_id] = next_node_id;
2547  next_node_id += this->n_processors() + 1;
2548  }
2549  }
2550  }
2551  }
2552  }
2553 
2554  // FIXME: ought to renumber side/node_id_map image to be contiguous
2555  // to save memory, also ought to reserve memory
2556 }
2557 
2558 
2559 } // namespace libMesh
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
void set_union(T &data, const unsigned int root_id) const
void add_elements(const std::set< boundary_id_type > &requested_boundary_ids, UnstructuredMesh &boundary_mesh)
std::size_t n_boundary_conds() const
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
std::map< boundary_id_type, std::string > _ns_id_to_name
void wait(std::vector< Request > &r)
Definition: request.C:213
std::set< boundary_id_type > _node_boundary_ids
const Elem * parent() const
Definition: elem.h:2479
void make_node_unique_ids_parallel_consistent(MeshBase &)
void remove_edge(const Elem *elem, const unsigned short int edge)
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
virtual unique_id_type parallel_max_unique_id() const =0
std::string & nodeset_name(boundary_id_type id)
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
const unsigned int invalid_uint
Definition: libmesh.h:245
void set_parent(Elem *p)
Definition: elem.h:2495
std::set< boundary_id_type > _edge_boundary_ids
virtual void libmesh_assert_valid_parallel_ids() const override
void sync(UnstructuredMesh &boundary_mesh)
void build_node_list_from_side_list()
const Elem * interior_parent() const
Definition: elem.C:804
const unsigned int any_source
Definition: communicator.h:70
std::size_t n_edge_conds() const
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2166
std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > build_edge_list() const
void remove_id(boundary_id_type id)
BoundaryInfo(MeshBase &m)
Definition: boundary_info.C:49
const Elem * top_parent() const
Definition: elem.h:2503
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const =0
void remove_shellface(const Elem *elem, const unsigned short int shellface)
static void set_node_processor_ids(MeshBase &mesh)
Definition: partitioner.C:679
std::size_t n_shellface_conds() const
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
unsigned int n_shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface) const
std::set< boundary_id_type > _boundary_ids
void remove(const Node *node)
unsigned short int side
Definition: xdr_io.C:50
unsigned int side_with_boundary_id(const Elem *const elem, const boundary_id_type boundary_id) const
The base class for all geometric element types.
Definition: elem.h:100
void print_info(std::ostream &out=libMesh::out) const
uint8_t processor_id_type
Definition: id_types.h:99
MPI_Request request
Definition: request.h:40
virtual std::unique_ptr< Partitioner > & partitioner()
Definition: mesh_base.h:126
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
const Parallel::Communicator & comm() const
std::set< boundary_id_type > _side_boundary_ids
void build_side_boundary_ids(std::vector< boundary_id_type > &b_ids) const
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
MessageTag get_unique_tag(int tagvalue) const
Definition: communicator.C:201
std::vector< boundary_id_type > raw_edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:131
std::multimap< const Elem *, std::pair< unsigned short int, boundary_id_type > > _boundary_side_id
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
void set_interior_parent(Elem *p)
Definition: elem.C:856
Base class for Mesh.
Definition: mesh_base.h:77
std::size_t n_boundary_ids() const
std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > build_active_side_list() const
std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > build_shellface_list() const
IntRange< unsigned short > edge_index_range() const
Definition: elem.h:2157
std::map< boundary_id_type, std::string > _ss_id_to_name
std::vector< boundary_id_type > boundary_ids(const Node *node) const
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
virtual element_iterator elements_begin()=0
boundary_id_type get_id_by_name(const std::string &name) const
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:154
static const subdomain_id_type invalid_subdomain_id
Definition: elem.h:262
void add_node(const Node *node, const boundary_id_type id)
const dof_id_type n_nodes
Definition: tecplot_io.C:68
int8_t boundary_id_type
Definition: id_types.h:51
static const boundary_id_type invalid_id
virtual SimpleRange< element_iterator > element_ptr_range()=0
dof_id_type id() const
Definition: dof_object.h:655
virtual unsigned int n_nodes() const =0
static const processor_id_type invalid_processor_id
Definition: dof_object.h:358
std::multimap< const Node *, boundary_id_type > _boundary_node_id
Base class for Replicated and Distributed meshes.
void build_side_list_from_node_list()
virtual Elem * add_elem(Elem *e)=0
void raw_shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
virtual element_iterator elements_end()=0
virtual const Node * query_node_ptr(const dof_id_type i) const =0
virtual SimpleRange< node_iterator > node_ptr_range()=0
virtual dof_id_type max_elem_id() const =0
Used by the Mesh to keep track of boundary nodes and elements.
Definition: boundary_info.h:57
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
MPI_Status status
Definition: status.h:41
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Definition: mesh_base.C:152
virtual unsigned int n_edges() const =0
boundary_id_type boundary_id(const Elem *const elem, const unsigned short int side) const
void set_neighbor(const unsigned int i, Elem *n)
Definition: elem.h:2083
Mesh data structure which is distributed across all processors.
An object whose state is distributed along a set of processors.
const std::string & get_nodeset_name(boundary_id_type id) const
virtual void clear()
Definition: mesh_base.C:260
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:2620
std::string & sideset_name(boundary_id_type id)
BoundaryInfo & operator=(const BoundaryInfo &other_boundary_info)
Definition: boundary_info.C:55
unsigned int n_partitions() const
Definition: mesh_base.h:875
virtual element_iterator pid_elements_begin(processor_id_type proc_id)=0
virtual const Elem * elem_ptr(const dof_id_type i) const =0
std::vector< boundary_id_type > raw_boundary_ids(const Elem *const elem, const unsigned short int side) const
void build_node_boundary_ids(std::vector< boundary_id_type > &b_ids) const
virtual unsigned int n_sides() const =0
void print_summary(std::ostream &out=libMesh::out) const
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2050
void remove_side(const Elem *elem, const unsigned short int side)
std::vector< std::tuple< dof_id_type, unsigned short int, boundary_id_type > > build_side_list() const
unsigned int level() const
Definition: elem.h:2521
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int n_edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
void copy_boundary_ids(const BoundaryInfo &old_boundary_info, const Elem *const old_elem, const Elem *const new_elem)
Temporarily serializes a DistributedMesh for output.
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1957
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
const std::string & get_sideset_name(boundary_id_type id) const
void add_shellface(const dof_id_type elem, const unsigned short int shellface, const boundary_id_type id)
std::vector< std::tuple< dof_id_type, boundary_id_type > > build_node_list() const
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
std::size_t n_nodeset_conds() const
std::set< boundary_id_type > _shellface_boundary_ids
void _find_id_maps(const std::set< boundary_id_type > &requested_boundary_ids, dof_id_type first_free_node_id, std::map< dof_id_type, dof_id_type > *node_id_map, dof_id_type first_free_elem_id, std::map< std::pair< dof_id_type, unsigned char >, dof_id_type > *side_id_map, const std::set< subdomain_id_type > &subdomains_relative_to)
std::multimap< const Elem *, std::pair< unsigned short int, boundary_id_type > > _boundary_shellface_id
virtual void delete_remote_elements()
Definition: mesh_base.h:196
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
OStreamProxy out(std::cout)
processor_id_type processor_id() const
Definition: dof_object.h:717
const Point & point(const unsigned int i) const
Definition: elem.h:1892
uint8_t unique_id_type
Definition: id_types.h:79
unsigned int & set_n_partitions()
Definition: mesh_base.h:1371
void build_shellface_boundary_ids(std::vector< boundary_id_type > &b_ids) const
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition: elem.C:1518
std::multimap< const Elem *, std::pair< unsigned short int, boundary_id_type > > _boundary_edge_id
void get_side_and_node_maps(UnstructuredMesh &boundary_mesh, std::map< dof_id_type, dof_id_type > &node_id_map, std::map< dof_id_type, unsigned char > &side_id_map, Real tolerance=1.e-6)
uint8_t dof_id_type
Definition: id_types.h:64
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
std::vector< boundary_id_type > edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
const RemoteElem * remote_elem
Definition: remote_elem.C:57
virtual element_iterator pid_elements_end(processor_id_type proc_id)=0