Allow boundary be associated with child elements

include/mesh/boundary_info.h

@@ -550,6 +550,20 @@ class BoundaryInfo : public ParallelObject
    */
   void build_shellface_boundary_ids(std::vector<boundary_id_type> & b_ids) const;
 
+  /**
+  * Update parent's boundary id list so that this information is consistent with
+  * its children elements
+  *
+  * This is useful when `_children_on_boundary = true`, and is used when the
+  * element is about to get coarsened i.e., in MeshRefinement::_coarsen_elements()
+  *
+  * Specifically, when coarsen a child element who is the last child on that
+  * boundary, we remove that boundary on the parent's side accordingly.
+  *
+  * Otherwise, we add the parent's side to the boundary.
+  */
+  void transfer_boundary_ids_to_parent(const Elem * const elem);
+
   /**
    * \returns The number of element-side-based boundary conditions.
    *
@@ -877,6 +891,18 @@ class BoundaryInfo : public ParallelObject
   const std::multimap<const Elem *, std::pair<unsigned short int, boundary_id_type>> & get_sideset_map() const
   { return _boundary_side_id; }
 
+  /**
+   * \returns Whether or not there may be child elements directly assigned boundary sides
+   */
+  bool is_children_on_boundary_side() const
+  { return _children_on_boundary; }
+
+  /**
+   * Whether or not to allow directly setting boundary sides on child elements
+   */
+  void allow_children_on_boundary_side(const bool children_on_boundary)
+  { _children_on_boundary = children_on_boundary; }
+
 private:
 
   /**
@@ -927,6 +953,13 @@ class BoundaryInfo : public ParallelObject
                 std::pair<unsigned short int, boundary_id_type>>
   _boundary_side_id;
 
+  /*
+   * Whether or not children elements are associated to any boundary
+   * It is false by default. The flag will be turned on if add_side
+   * function is called with a child element
+   */
+  bool _children_on_boundary;
+
   /**
    * A collection of user-specified boundary ids for sides, edges, nodes,
    * and shell faces.

include/mesh/mesh_refinement.h

@@ -241,6 +241,10 @@ class MeshRefinement : public ParallelObject
    *
    * \note This function used to take an argument, \p maintain_level_one,
    * new code should use face_level_mismatch_limit() instead.
+   *
+   * \note When we allow boundaries to be directly associated with child elements,
+   * i.e., `_children_on_boundary = true`. A child's boundary ID may be
+   * lost during coarsening if it differs from its siblings on that parent side.
    */
   bool coarsen_elements ();
 

src/mesh/boundary_info.C

@@ -101,7 +101,8 @@ const boundary_id_type BoundaryInfo::invalid_id = -123;
 // BoundaryInfo functions
 BoundaryInfo::BoundaryInfo(MeshBase & m) :
   ParallelObject(m.comm()),
-  _mesh (&m)
+  _mesh (&m),
+  _children_on_boundary(false)
 {
 }
 
@@ -962,8 +963,11 @@ void BoundaryInfo::add_side(const Elem * elem,
 {
   libmesh_assert(elem);
 
-  // Only add BCs for level-0 elements.
-  libmesh_assert_equal_to (elem->level(), 0);
+  // Users try to mark boundary on child elements
+  // If this happens, we will allow users to remove
+  // side from child elements as well
+  if (elem->level())
+    _children_on_boundary = true;
 
   libmesh_error_msg_if(id == invalid_id, "ERROR: You may not set a boundary ID of "
                        << invalid_id
@@ -991,8 +995,11 @@ void BoundaryInfo::add_side(const Elem * elem,
 
   libmesh_assert(elem);
 
-  // Only add BCs for level-0 elements.
-  libmesh_assert_equal_to (elem->level(), 0);
+  // Users try to mark boundary on child elements
+  // If this happens, we will allow users to remove
+  // side from child elements as well
+  if (elem->level())
+    _children_on_boundary = true;
 
   // Don't add the same ID twice
   auto bounds = _boundary_side_id.equal_range(elem);
@@ -1240,8 +1247,15 @@ void BoundaryInfo::boundary_ids (const Elem * const elem,
   // Clear out any previous contents
   vec_to_fill.clear();
 
-  // Only level-0 elements store BCs.  If this is not a level-0
-  // element get its level-0 parent and infer the BCs.
+  // Search BC on the current element. If we find anything, we should return
+  for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
+    if (pr.second.first == side)
+      vec_to_fill.push_back(pr.second.second);
+
+  if (vec_to_fill.size())
+    return;
+
+  // We check the top parent now
   const Elem * searched_elem = elem;
   if (elem->level() != 0)
     {
@@ -1288,7 +1302,7 @@ void BoundaryInfo::raw_boundary_ids (const Elem * const elem,
   vec_to_fill.clear();
 
   // Only level-0 elements store BCs.
-  if (elem->parent())
+  if (elem->parent() && !_children_on_boundary)
     return;
 
   // Check each element in the range to see if its side matches the requested side.
@@ -1438,8 +1452,8 @@ void BoundaryInfo::remove_side (const Elem * elem,
 {
   libmesh_assert(elem);
 
-  // Only level 0 elements are stored in BoundaryInfo.
-  libmesh_assert_equal_to (elem->level(), 0);
+  // Only level 0 elements unless the flag "_children_on_boundary" is on.
+  libmesh_assert(elem->level() == 0 || _children_on_boundary);
 
   // Erase (elem, side, *) entries from map.
   erase_if(_boundary_side_id, elem,
@@ -1455,6 +1469,9 @@ void BoundaryInfo::remove_side (const Elem * elem,
 {
   libmesh_assert(elem);
 
+  // Only level 0 elements unless the flag "_children_on_boundary" is on.
+  libmesh_assert(elem->level() == 0 || _children_on_boundary);
+
   // Erase (elem, side, id) entries from map.
   erase_if(_boundary_side_id, elem,
            [side, id](decltype(_boundary_side_id)::mapped_type & pr)
@@ -1575,12 +1592,28 @@ void BoundaryInfo::renumber_id (boundary_id_type old_id,
 unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
                                                  const boundary_id_type boundary_id_in) const
 {
-  const Elem * searched_elem = elem;
+  std::vector<const Elem *>  searched_elem_vec;
+  // If elem has boundary information, we return that when
+  // the flag "_children_on_boundary" is on
+  if (_children_on_boundary)
+    searched_elem_vec.push_back(elem);
+  // Otherwise, we return boundary information of all its ancestors
   if (elem->level() != 0)
-    searched_elem = elem->top_parent();
+  {
+    const Elem * parent = elem->parent();
+    while (parent != nullptr)
+    {
+      searched_elem_vec.push_back(parent);
+      parent = parent->parent();
+    }
+  }
+  else if (!_children_on_boundary)
+    searched_elem_vec.push_back(elem);
 
-  // elem may have zero or multiple occurrences
-  for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
+  for (const Elem * searched_elem : searched_elem_vec)
+  {
+    // elem may have zero or multiple occurrences
+    for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
     {
       // if this is true we found the requested boundary_id
       // of the element and want to return the side
@@ -1611,7 +1644,8 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
           if (!p)
             return side;
         }
-    }
+     }
+  }
 
   // if we get here, we found elem in the data structure but not
   // the requested boundary id, so return the default value
@@ -1625,12 +1659,22 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
 {
   std::vector<unsigned int> returnval;
 
-  const Elem * searched_elem = elem;
+  std::vector<const Elem *>  searched_elem_vec;
+  // If elem has boundary information, that is part of return when
+  // the flag "_children_on_boundary" is on
+  if (_children_on_boundary)
+    searched_elem_vec.push_back(elem);
+  // Return boundary information of its parent as well
   if (elem->level() != 0)
-    searched_elem = elem->top_parent();
+    searched_elem_vec.push_back(elem->top_parent());
+  else if (!_children_on_boundary)
+    searched_elem_vec.push_back(elem);
 
-  // elem may have zero or multiple occurrences
-  for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
+  for (auto it = searched_elem_vec.begin(); it != searched_elem_vec.end(); ++it)
+  {
+    const Elem * searched_elem = *it;
+    // elem may have zero or multiple occurrences
+    for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
     {
       // if this is true we found the requested boundary_id
       // of the element and want to return the side
@@ -1665,7 +1709,7 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
             returnval.push_back(side);
         }
     }
-
+  }
   return returnval;
 }
 
@@ -1712,6 +1756,48 @@ BoundaryInfo::build_shellface_boundary_ids(std::vector<boundary_id_type> & b_ids
     }
 }
 
+void
+BoundaryInfo::transfer_boundary_ids_to_parent(const Elem * const elem)
+{
+  // this is only needed when we allow boundary to be associated with children elements
+  // also, we only transfer the parent's boundary ids when we are actually coarsen the child element
+  if (!_children_on_boundary ||
+      !elem->active() ||
+      elem->level()==0 ||
+      elem->refinement_flag() != Elem::COARSEN)
+    return;
+
+  const Elem * parent = elem->parent();
+
+  for (const auto & pr : as_range(_boundary_side_id.equal_range(elem)))
+  {
+    auto side = pr.second.first;
+    auto bnd_id = pr.second.second;
+    // Track if any of the sibling elements is on this boundary.
+    // If yes, we make sure that the corresponding parent side is added to the boundary.
+    // Otherwise, we remove the parent side from the boundary.
+    std::vector<const Elem *> family_on_side;
+    elem->family_tree_by_side(family_on_side, side);
+    for (auto relative : family_on_side)
+    {
+      if (relative != elem && relative->level() == elem->level()) // check only siblings
+      {
+        for (unsigned int i = 0; i < relative->n_sides(); ++i)
+          if (this->has_boundary_id(relative, i, bnd_id))
+          {
+            // Do not worry, `add_side` will avoid adding duplicate sides on the same boundary
+            // Note: it is assumed that the child and parent elements' side numberings are identical.
+            // I.e., a child's ith side is encompassed in the parent's jth side, where i=j.
+            this->add_side(parent, side, bnd_id);
+            return;
+          }
+      }
+    }
+    // No relative share the same boundary, therefore, we remove it from the parent, if any.
+    this->remove_side(parent, side, bnd_id);
+  }
+}
+
 std::size_t BoundaryInfo::n_boundary_conds () const
 {
   // in serial we know the number of bcs from the
@@ -1880,7 +1966,8 @@ BoundaryInfo::build_node_list_from_side_list()
       // Need to loop over the sides of any possible children
       std::vector<const Elem *> family;
 #ifdef LIBMESH_ENABLE_AMR
-      elem->active_family_tree_by_side (family, id_pair.first);
+      if (!elem->subactive())
+        elem->active_family_tree_by_side (family, id_pair.first);
 #else
       family.push_back(elem);
 #endif
@@ -2242,7 +2329,8 @@ BoundaryInfo::build_active_side_list () const
       // Loop over the sides of possible children
       std::vector<const Elem *> family;
 #ifdef LIBMESH_ENABLE_AMR
-      elem->active_family_tree_by_side(family, id_pair.first);
+      if (!elem->subactive())
+        elem->active_family_tree_by_side(family, id_pair.first);
 #else
       family.push_back(elem);
 #endif

src/mesh/mesh_refinement.C

@@ -1371,6 +1371,11 @@ bool MeshRefinement::_coarsen_elements ()
 
   for (auto & elem : _mesh.element_ptr_range())
     {
+      // Make sure we transfer the element's boundary id(s)
+      // up to its parent when necessary before coarsening.
+      // This can be adding or removing the corresonding boundary info.
+      _mesh.get_boundary_info().transfer_boundary_ids_to_parent(elem);
+
       // active elements flagged for coarsening will
       // no longer be deleted until MeshRefinement::contract()
       if (elem->refinement_flag() == Elem::COARSEN)
@@ -1384,8 +1389,11 @@ bool MeshRefinement::_coarsen_elements ()
           elem->nullify_neighbors();
 
           // Remove any boundary information associated
-          // with this element
-          _mesh.get_boundary_info().remove (elem);
+          // with this element if we do not allow children to have boundary info
+          // otherwise we will have trouble in boundary info consistency among
+          // parent and children elements
+          if (!_mesh.get_boundary_info().is_children_on_boundary_side())
+            _mesh.get_boundary_info().remove (elem);
 
           // Add this iterator to the _unused_elements
           // data structure so we might fill it.