Modify parallel packing for elements.

Author: grmnptr

include/mesh/boundary_info.h

@@ -480,6 +480,13 @@ class BoundaryInfo : public ParallelObject
   unsigned int n_boundary_ids (const Elem * const elem,
                                const unsigned short int side) const;
 
+  /**
+   * \returns The number of raw (excludes ancestors) boundary ids associated with the \p side
+   * side of element \p elem.
+   */
+  unsigned int n_raw_boundary_ids (const Elem * const elem,
+                                   const unsigned short int side) const;
+
   /**
    * \returns The list of boundary ids associated with the \p side side of
    * element \p elem.
@@ -552,15 +559,15 @@ class BoundaryInfo : public ParallelObject
 
   /**
   * Update parent's boundary id list so that this information is consistent with
-  * its children elements
+  * its children
   *
   * 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.
+  * Specifically, when we coarsen an element whose children have different boundary ids.
+  * In such scenarios, the parent will inherit the children's boundaries if at
+  * least two of them own a boundary while sharing the side of the parent.
+  * Otherwise, we delete the boundary from the children and the parent as well.
   */
   void transfer_boundary_ids_to_parent(const Elem * const elem);
 
@@ -954,8 +961,8 @@ class BoundaryInfo : public ParallelObject
   _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
+   * Whether or not children elements are associated with 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;

src/mesh/boundary_info.C

@@ -1268,21 +1268,22 @@ void BoundaryInfo::boundary_ids (const Elem * const elem,
                       vec_to_fill.end())
                 vec_to_fill.push_back(pr.second.second);
 
+      // Loop over ancestors to check if they have boundary ids on the same side
       while (searched_elem->parent() != nullptr)
       {
-              const Elem * parent = searched_elem->parent();
-              if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
-                return;
-
-              searched_elem = parent;
-
-              if (_children_on_boundary)
-                for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
-                      // Here we need to check if the boundary id already exists
-                      if (pr.second.first == side &&
-                          std::find(vec_to_fill.begin(), vec_to_fill.end(), pr.second.second) ==
-                              vec_to_fill.end())
-                          vec_to_fill.push_back(pr.second.second);
+        const Elem * parent = searched_elem->parent();
+        if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
+          return;
+
+        searched_elem = parent;
+
+        if (_children_on_boundary)
+          for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
+            // Here we need to check if the boundary id already exists
+            if (pr.second.first == side &&
+                std::find(vec_to_fill.begin(), vec_to_fill.end(), pr.second.second) ==
+                vec_to_fill.end())
+              vec_to_fill.push_back(pr.second.second);
       }
 
       return;
@@ -1296,11 +1297,11 @@ void BoundaryInfo::boundary_ids (const Elem * const elem,
     else
       while (searched_elem->parent() != nullptr)
       {
-              const Elem * parent = searched_elem->parent();
-              if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
-                return;
+        const Elem * parent = searched_elem->parent();
+        if (parent->is_child_on_side(parent->which_child_am_i(searched_elem), side) == false)
+          return;
 
-              searched_elem = parent;
+        searched_elem = parent;
       }
   }
 
@@ -1324,6 +1325,15 @@ unsigned int BoundaryInfo::n_boundary_ids (const Elem * const elem,
 }
 
 
+unsigned int BoundaryInfo::n_raw_boundary_ids (const Elem * const elem,
+                                               const unsigned short int side) const
+{
+  std::vector<boundary_id_type> ids;
+  this->raw_boundary_ids(elem, side, ids);
+  return cast_int<unsigned int>(ids.size());
+}
+
+
 
 void BoundaryInfo::raw_boundary_ids (const Elem * const elem,
                                      const unsigned short int side,
@@ -1499,9 +1509,6 @@ 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)
@@ -1661,7 +1668,7 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
             p = parent;
           }
 #endif
-              // We're on that side of our top_parent; return it
+          // We're on that side of our top_parent; return it
           if (!p)
             return side;
         }
@@ -1719,7 +1726,7 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
 
   const Elem * searched_elem = elem;
   if (elem->level() != 0 && !_children_on_boundary)
-      searched_elem = elem->top_parent();
+    searched_elem = elem->top_parent();
 
   // elem may have zero or multiple occurrences
   for (const auto & pr : as_range(_boundary_side_id.equal_range(searched_elem)))
@@ -1737,10 +1744,10 @@ BoundaryInfo::sides_with_boundary_id(const Elem * const elem,
           // If we're on this external boundary then we share this
           // external boundary id
           if (elem->neighbor_ptr(side) == nullptr)
-          {
-            returnval.push_back(side);
-            continue;
-          }
+            {
+              returnval.push_back(side);
+              continue;
+            }
 
           // If we're on an internal boundary then we need to be sure
           // it's the same internal boundary as our top_parent
@@ -2069,7 +2076,6 @@ 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
-      // if (!elem->subactive())
       elem->active_family_tree_by_side (family, id_pair.first);
 #else
       family.push_back(elem);
@@ -2432,8 +2438,7 @@ BoundaryInfo::build_active_side_list () const
       // Loop over the sides of possible children
       std::vector<const Elem *> family;
 #ifdef LIBMESH_ENABLE_AMR
-      // if (!elem->subactive())
-        elem->active_family_tree_by_side(family, id_pair.first);
+      elem->active_family_tree_by_side(family, id_pair.first);
 #else
       family.push_back(elem);
 #endif

src/mesh/mesh_refinement.C

@@ -659,7 +659,6 @@ bool MeshRefinement::coarsen_elements ()
   if (_face_level_mismatch_limit)
     libmesh_assert(test_level_one(true));
   libmesh_assert(this->make_coarsening_compatible());
-
   // FIXME: This won't pass unless we add a redundant find_neighbors()
   // call or replace find_neighbors() with on-the-fly neighbor updating
   // libmesh_assert(!this->eliminate_unrefined_patches());
@@ -706,8 +705,6 @@ bool MeshRefinement::refine_elements ()
         elem->set_refinement_flag(Elem::DO_NOTHING);
     }
 
-
-
   // Parallel consistency has to come first, or coarsening
   // along processor boundaries might occasionally be falsely
   // prevented
@@ -736,6 +733,7 @@ bool MeshRefinement::refine_elements ()
   if (_face_level_mismatch_limit)
     libmesh_assert(test_level_one(true));
   libmesh_assert(this->make_refinement_compatible());
+
   // FIXME: This won't pass unless we add a redundant find_neighbors()
   // call or replace find_neighbors() with on-the-fly neighbor updating
   // libmesh_assert(!this->eliminate_unrefined_patches());
@@ -1393,8 +1391,8 @@ bool MeshRefinement::_coarsen_elements ()
           // 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);
+          // 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.
@@ -1406,8 +1404,6 @@ bool MeshRefinement::_coarsen_elements ()
           // _mesh.delete_elem(elem);
         }
 
-
-
       // inactive elements flagged for coarsening
       // will become active
       else if (elem->refinement_flag() == Elem::COARSEN_INACTIVE)
@@ -1417,8 +1413,6 @@ bool MeshRefinement::_coarsen_elements ()
 
           // the mesh has certainly changed
           mesh_changed = true;
-
-
         }
       if (elem->p_refinement_flag() == Elem::COARSEN)
         {
@@ -1433,7 +1427,6 @@ bool MeshRefinement::_coarsen_elements ()
               elem->set_p_refinement_flag(Elem::DO_NOTHING);
             }
         }
-
     }
 
   this->comm().max(mesh_p_changed);
@@ -1460,8 +1453,6 @@ bool MeshRefinement::_coarsen_elements ()
 #endif
     }
 
-
-
   // If p levels changed all we need to do is make sure that parent p
   // levels changed in sync
   if (mesh_p_changed && !_mesh.is_serial())
@@ -1470,8 +1461,6 @@ bool MeshRefinement::_coarsen_elements ()
     }
 
   return (mesh_changed || mesh_p_changed);
-
-
 }