Transfer child element's boundary info (add or remove) during coarsening, if the child element has boundary info associated with itself

Author: dewenyushu

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.
    *

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 boundary be associated with children elements,
+   * i.e., `_children_on_boundary = true`. A tpyical child boundary ID may be
+   * lost during coarsening if that child has different IDs from its siblings.
    */
   bool coarsen_elements ();
 

src/mesh/boundary_info.C

@@ -1600,7 +1600,7 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
   // Otherwise, we return boundary information of all its ancestors
   if (elem->level() != 0)
   {
-    Elem * parent = elem->parent();
+    const Elem * parent = elem->parent();
     while (parent != nullptr)
     {
       searched_elem_vec.push_back(parent);
@@ -1757,6 +1757,47 @@ 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
+  if (!_children_on_boundary || !elem->active())
+    return;
+
+  // we only transfer the parent's boundary ids when we are actually coarsen the child element
+  libmesh_assert (elem->level()!=0 && elem->refinement_flag() == Elem::COARSEN);
+
+  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

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.