Merge pull request #4182 from roystgnr/interior_parent_parallel

Track a MeshBase::interior_mesh() to enable BoundaryMesh serialization

Author: roystgnr

include/mesh/mesh_base.h

@@ -1787,6 +1787,25 @@ class MeshBase : public ParallelObject
    */
   void detect_interior_parents();
 
+  /**
+   * \return A mesh that may own interior parents of elements in this
+   * mesh.  In most cases this mesh includes its own interior parents,
+   * but in cases where a separate "interior" mesh was used to create
+   * this mesh as a distinct lower-dimensional boundary (or boundary
+   * subset) mesh, the original mesh will be returned here.
+   */
+  const MeshBase & interior_mesh() const { return *_interior_mesh; }
+
+  /**
+   * \return A writeable reference to the interior mesh.
+   */
+  MeshBase & interior_mesh() { return *_interior_mesh; }
+
+  /**
+   * Sets the interior mesh.  For advanced use only.
+   */
+  void set_interior_mesh(MeshBase & int_mesh) { _interior_mesh = &int_mesh; }
+
   /**
    * \return The cached mesh subdomains. As long as the mesh is prepared, this
    * should contain all the subdomain ids across processors. Relies on the mesh
@@ -1906,6 +1925,12 @@ class MeshBase : public ParallelObject
   unique_id_type _next_unique_id;
 #endif
 
+  /**
+   * Defaulting to \p this, a pointer to the mesh used to generate
+   * boundary elements on \p this.
+   */
+  MeshBase *_interior_mesh;
+
   /**
    * If this is true then no partitioning should be done with the
    * possible exception of orphaned nodes.

src/mesh/boundary_info.C

@@ -478,12 +478,18 @@ void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_i
    * into dangling pointers, and it won't be easy to tell which.  So
    * we'll keep around a distributed copy for that case, and query it
    * to fix up interior_parent() links as necessary.
+   *
+   * We'll also need to make sure to unserialize the mesh *before* we
+   * prepare the boundary mesh for use, or the prepare_for_use() call
+   * on a refined boundary mesh will happily notice that it can find
+   * and restore some refined elements' interior_parent pointers, not
+   * knowing that those interior parents are about to go remote.
    */
   std::unique_ptr<MeshBase> mesh_copy;
   if (boundary_mesh.is_serial() && !_mesh->is_serial())
     mesh_copy = _mesh->clone();
 
-  MeshSerializer serializer
+  auto serializer = std::make_unique<MeshSerializer>
     (const_cast<MeshBase &>(*_mesh), boundary_mesh.is_serial());
 
   /**
@@ -566,6 +572,11 @@ void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_i
   // interior partitioning.
   boundary_mesh.partitioner().reset(nullptr);
 
+  // Deserialize the interior mesh before the boundary mesh
+  // prepare_for_use() can come to erroneous conclusions about which
+  // of its elements are semilocal
+  serializer.reset();
+
   // Make boundary_mesh nodes and elements contiguous
   boundary_mesh.prepare_for_use();
 
@@ -650,10 +661,22 @@ void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_bou
   LOG_SCOPE("add_elements()", "BoundaryInfo");
 
   // We're not prepared to mix serial and distributed meshes in this
-  // method, so make sure they match from the start.
+  // method, so make sure their statuses match from the start.
+  //
+  // Specifically test *is_serial* here - we can handle a mix of
+  // ReplicatedMesh and serialized DistributedMesh.
   libmesh_assert_equal_to(_mesh->is_serial(),
                           boundary_mesh.is_serial());
 
+  // If the boundary mesh already has interior pointers pointing at
+  // elements in a third mesh then we're in trouble
+  libmesh_assert(&boundary_mesh.interior_mesh() == &boundary_mesh ||
+                 &boundary_mesh.interior_mesh() == _mesh);
+
+  // And now we're going to add interior pointers to elements from
+  // this mesh
+  boundary_mesh.set_interior_mesh(*_mesh);
+
   std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> side_id_map;
   this->_find_id_maps(requested_boundary_ids,
                       0,
@@ -682,25 +705,18 @@ void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_bou
           !subdomains_relative_to.count(elem->subdomain_id()))
         continue;
 
-      // Get the top-level parent for this element
-      const Elem * top_parent = elem->top_parent();
-
-      // Find all the boundary side ids for this Elem.
-      auto bounds = _boundary_side_id.equal_range(top_parent);
-
       for (auto s : elem->side_index_range())
         {
           bool add_this_side = false;
-          boundary_id_type this_bcid = invalid_id;
 
-          for (const auto & pr : as_range(bounds))
-            {
-              this_bcid = pr.second.second;
+          // Find all the boundary side ids for this Elem side.
+          std::vector<boundary_id_type> bcids;
+          this->boundary_ids(elem, s, bcids);
 
-              // if this side is flagged with a boundary condition
-              // and the user wants this id
-              if ((pr.second.first == s) &&
-                  (requested_boundary_ids.count(this_bcid)))
+          for (const boundary_id_type bcid : bcids)
+            {
+              // if the user wants this id, we want this side
+              if (requested_boundary_ids.count(bcid))
                 {
                   add_this_side = true;
                   break;
@@ -713,8 +729,7 @@ void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_bou
           // boundary was copied to the BoundaryMesh, and handles the
           // case where elements on the geometric boundary are not in
           // any sidesets.
-          if (bounds.first == bounds.second            &&
-              requested_boundary_ids.count(invalid_id) &&
+          if (requested_boundary_ids.count(invalid_id) &&
               elem->neighbor_ptr(s) == nullptr)
             add_this_side = true;
 
@@ -765,6 +780,9 @@ void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_bou
         new_elem->set_extra_integer(parent_side_index_tag, s);
 
 #ifdef LIBMESH_ENABLE_AMR
+      new_elem->set_refinement_flag(elem->refinement_flag());
+      new_elem->set_p_refinement_flag(elem->p_refinement_flag());
+
       // Set parent links
       if (elem->parent())
         {
@@ -778,8 +796,6 @@ void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_bou
 
           new_elem->set_parent(side_parent);
 
-          side_parent->set_refinement_flag(Elem::INACTIVE);
-
           // Figuring out which child we are of our parent
           // is a trick.  Due to libMesh child numbering
           // conventions, if we are an element on a vertex,
@@ -802,6 +818,22 @@ void BoundaryInfo::add_elements(const std::set<boundary_id_type> & requested_bou
               side_parent->add_child(new_elem, 3);
             }
         }
+
+      // Set remote_elem child links if necessary.  Rather than
+      // worrying about which interior child corresponds to which side
+      // child we'll just set all null links to be remote and we'll
+      // rely on our detection of actual semilocal children to
+      // overwrite the links that shouldn't be remote.
+      if (elem->has_children())
+        for (auto c : make_range(elem->n_children()))
+          if (elem->child_ptr(c) == remote_elem &&
+              elem->is_child_on_side(c, s))
+            {
+              for (auto sc : make_range(new_elem->n_children()))
+                if (!new_elem->raw_child_ptr(sc))
+                  new_elem->add_child
+                    (const_cast<RemoteElem*>(remote_elem), sc);
+            }
 #endif
 
       new_elem->set_interior_parent (elem);
@@ -3136,26 +3168,18 @@ void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & requested_bo
           !subdomains_relative_to.count(elem->subdomain_id()))
         continue;
 
-      // Get the top-level parent for this element. This is used for
-      // searching for boundary sides on this element.
-      const Elem * top_parent = elem->top_parent();
-
-      // Find all the boundary side ids for this Elem.
-      auto bounds = _boundary_side_id.equal_range(top_parent);
-
       for (auto s : elem->side_index_range())
         {
           bool add_this_side = false;
-          boundary_id_type this_bcid = invalid_id;
 
-          for (const auto & pr : as_range(bounds))
-            {
-              this_bcid = pr.second.second;
+          // Find all the boundary side ids for this Elem side.
+          std::vector<boundary_id_type> bcids;
+          this->boundary_ids(elem, s, bcids);
 
-              // if this side is flagged with a boundary condition
-              // and the user wants this id
-              if ((pr.second.first == s) &&
-                  (requested_boundary_ids.count(this_bcid)))
+          for (const boundary_id_type bcid : bcids)
+            {
+              // if the user wants this id, we want this side
+              if (requested_boundary_ids.count(bcid))
                 {
                   add_this_side = true;
                   break;
@@ -3168,8 +3192,7 @@ void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & requested_bo
           // boundary was copied to the BoundaryMesh, and handles the
           // case where elements on the geometric boundary are not in
           // any sidesets.
-          if (bounds.first == bounds.second            &&
-              requested_boundary_ids.count(invalid_id) &&
+          if (requested_boundary_ids.count(invalid_id) &&
               elem->neighbor_ptr(s) == nullptr)
             add_this_side = true;
 

src/mesh/mesh_base.C

@@ -70,6 +70,7 @@ MeshBase::MeshBase (const Parallel::Communicator & comm_in,
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
   _next_unique_id(DofObject::invalid_unique_id),
 #endif
+  _interior_mesh(this),
   _skip_noncritical_partitioning(false),
   _skip_all_partitioning(libMesh::on_command_line("--skip-partitioning")),
   _skip_renumber_nodes_and_elements(false),
@@ -101,6 +102,10 @@ MeshBase::MeshBase (const MeshBase & other_mesh) :
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
   _next_unique_id(other_mesh._next_unique_id),
 #endif
+  // If the other mesh interior_parent pointers just go back to
+  // itself, so should we
+  _interior_mesh((other_mesh._interior_mesh == &other_mesh) ?
+                 this : other_mesh._interior_mesh),
   _skip_noncritical_partitioning(other_mesh._skip_noncritical_partitioning),
   _skip_all_partitioning(other_mesh._skip_all_partitioning),
   _skip_renumber_nodes_and_elements(other_mesh._skip_renumber_nodes_and_elements),
@@ -168,9 +173,13 @@ MeshBase& MeshBase::operator= (MeshBase && other_mesh)
   _is_prepared = other_mesh.is_prepared();
   _point_locator = std::move(other_mesh._point_locator);
   _count_lower_dim_elems_in_point_locator = other_mesh.get_count_lower_dim_elems_in_point_locator();
-  #ifdef LIBMESH_ENABLE_UNIQUE_ID
-    _next_unique_id = other_mesh.next_unique_id();
-  #endif
+#ifdef LIBMESH_ENABLE_UNIQUE_ID
+  _next_unique_id = other_mesh.next_unique_id();
+#endif
+  // If the other mesh interior_parent pointers just go back to
+  // itself, so should we
+  _interior_mesh = (other_mesh._interior_mesh == &other_mesh) ?
+                   this : other_mesh._interior_mesh;
   _skip_noncritical_partitioning = other_mesh.skip_noncritical_partitioning();
   _skip_all_partitioning = other_mesh.skip_partitioning();
   _skip_renumber_nodes_and_elements = !(other_mesh.allow_renumbering());
@@ -247,6 +256,15 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const
   if (_count_lower_dim_elems_in_point_locator !=
         other_mesh._count_lower_dim_elems_in_point_locator)
     return false;
+
+  // We should either both have our own interior parents or both not;
+  // but if we both don't then we can't really assert anything else
+  // because pointing at the same interior mesh is fair but so is
+  // pointing at two different copies of "the same" interior mesh.
+  if ((_interior_mesh == this) !=
+      (other_mesh._interior_mesh == &other_mesh))
+    return false;
+
   if (_skip_noncritical_partitioning != other_mesh._skip_noncritical_partitioning)
     return false;
   if (_skip_all_partitioning != other_mesh._skip_all_partitioning)
@@ -1794,11 +1812,17 @@ void MeshBase::detect_interior_parents()
   // This requires an inspection on every processor
   parallel_object_only();
 
-  // Check if the mesh contains mixed dimensions. If so, then set interior parents, otherwise return.
+  // Check if the mesh contains mixed dimensions. If so, then we may
+  // have interior parents to set.  Otherwise return.
   if (this->elem_dimensions().size() == 1)
     return;
 
-  //This map will be used to set interior parents
+  // Do we have interior parent pointers going to a different mesh?
+  // If so then we'll still check to make sure that's the only place
+  // they go, so we can libmesh_not_implemented() if not.
+  const bool separate_interior_mesh = (&(this->interior_mesh()) != this);
+
+  // This map will be used to set interior parents
   std::unordered_map<dof_id_type, std::vector<dof_id_type>> node_to_elem;
 
   for (const auto & elem : this->element_ptr_range())
@@ -1877,6 +1901,14 @@ void MeshBase::detect_interior_parents()
                   break;
                 }
             }
+
+          // Do we have a mixed dimensional mesh that contains some of
+          // its own interior parents, but we already expect to have
+          // interior parents on a different mesh?  That's going to
+          // take some work to support if anyone needs it.
+          if (separate_interior_mesh)
+            libmesh_not_implemented_msg
+              ("interior_parent() values in multiple meshes are unsupported.");
         }
     }
 }

src/parallel/parallel_elem.C

@@ -579,7 +579,8 @@ Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
           }
         else
           {
-            Elem * ip = mesh->query_elem_ptr(interior_parent_id);
+            Elem * ip =
+              mesh->interior_mesh().query_elem_ptr(interior_parent_id);
 
             // The sending processor sees an interior parent here, so
             // if we don't have that interior element, then we'd
@@ -794,7 +795,8 @@ Packing<Elem *>::unpack (std::vector<largest_id_type>::const_iterator in,
           {
             // If we don't have the interior parent element, then it's
             // a remote_elem until we get it.
-            Elem * ip = mesh->query_elem_ptr(interior_parent_id);
+            Elem * ip =
+              mesh->interior_mesh().query_elem_ptr(interior_parent_id);
             if (!ip )
               elem->set_interior_parent
                 (const_cast<RemoteElem *>(remote_elem));