Enable disconnected neighbor detection across disconnected boundary pairs in libMesh

include/mesh/mesh_base.h

@@ -36,6 +36,16 @@
 #include <string>
 #include <memory>
 
+#include "libmesh/vector_value.h"
+
+// periodic boundary condition support
+// Use forward declarations inside the libMesh namespace
+namespace libMesh
+{
+  class PeriodicBoundary;
+  class PeriodicBoundaries;
+}
+
 namespace libMesh
 {
 
@@ -1820,8 +1830,30 @@ class MeshBase : public ParallelObject
   const std::set<subdomain_id_type> & get_mesh_subdomains() const
   { libmesh_assert(this->is_prepared()); return _mesh_subdomains; }
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  /**
+   * Register a pair of boundaries as disjoint neighbor boundary pairs.
+   */
+  void add_disjoint_neighbor_boundary_pairs(const boundary_id_type b1,
+                                   const boundary_id_type b2,
+                                   const RealVectorValue & translation);
+
+  PeriodicBoundaries * get_disjoint_neighbor_boundary_pairs();
+
+  const PeriodicBoundaries * get_disjoint_neighbor_boundary_pairs() const;
+
+  void remove_disjoint_boundary_pair(const boundary_id_type b1,
+                                     const boundary_id_type b2);
+#endif
+
+
 protected:
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  /// @brief The disjoint boundary id pairs.
+  std::unique_ptr<PeriodicBoundaries> _disjoint_boundary_pairs;
+#endif
+
   /**
    * This class holds the boundary information.  It can store nodes, edges,
    * and faces with a corresponding id that facilitates setting boundary

src/ghosting/ghost_point_neighbors.C

@@ -46,8 +46,11 @@ void GhostPointNeighbors::operator()
   bool check_periodic_bcs =
     (_periodic_bcs && !_periodic_bcs->empty());
 
+  auto * db = _mesh->get_disjoint_neighbor_boundary_pairs();
+  bool check_disjoint_bcs = (db && !db->empty());
+
   std::unique_ptr<PointLocatorBase> point_locator;
-  if (check_periodic_bcs)
+  if (check_periodic_bcs || check_disjoint_bcs)
       point_locator = _mesh->sub_point_locator();
 
   std::set<const Elem *> periodic_elems_examined;
@@ -176,6 +179,30 @@ void GhostPointNeighbors::operator()
                 }
             }
         }
+
+      if (check_disjoint_bcs)
+        {
+          // Also ghost their disjoint neighbors
+          for (auto s : elem->side_index_range())
+          {
+            for (const auto & [id, boundary_ptr] : *db)
+            {
+              if (!_mesh->get_boundary_info().has_boundary_id(elem, s, id))
+                continue;
+
+              unsigned int neigh_side = invalid_uint;
+              const Elem * neigh =
+                  db->neighbor(id, *point_locator, elem, s, &neigh_side);
+
+              if (!neigh || neigh == remote_elem)
+                continue;
+
+              if (neigh->processor_id() != p)
+                coupled_elements.emplace(neigh, nullcm);
+            }
+          }
+        }
+
 #endif // LIBMESH_ENABLE_PERIODIC
     }
 }

src/mesh/boundary_info.C

@@ -2325,19 +2325,22 @@ unsigned int BoundaryInfo::side_with_boundary_id(const Elem * const elem,
           if (elem->neighbor_ptr(side) == nullptr)
             return side;
 
-          // If we're on an internal boundary then we need to be sure
-          // it's the same internal boundary as our top_parent
+          // Internal boundary case
           const Elem * p = elem;
 
 #ifdef LIBMESH_ENABLE_AMR
-
+          // If we're on an internal boundary then we need to be sure
+          // it's the same internal boundary as our top_parent
           while (p != nullptr)
           {
             const Elem * parent = p->parent();
             if (parent && !parent->is_child_on_side(parent->which_child_am_i(p), side))
               break;
             p = parent;
           }
+#else
+          // do not forget to return the internal boundary when AMR is disabled
+          return side;
 #endif
           // We're on that side of our top_parent; return it
           if (!p)

src/mesh/mesh_base.C

@@ -49,6 +49,9 @@
 #include <sstream>   // for std::ostringstream
 #include <unordered_map>
 
+#include "libmesh/periodic_boundaries.h"
+#include "libmesh/periodic_boundary.h"
+
 namespace libMesh
 {
 
@@ -154,6 +157,18 @@ MeshBase::MeshBase (const MeshBase & other_mesh) :
   if (other_mesh._partitioner.get())
     _partitioner = other_mesh._partitioner->clone();
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  // Deep copy of all periodic boundaries
+  if (other_mesh._disjoint_boundary_pairs)
+    {
+      _disjoint_boundary_pairs = std::make_unique<PeriodicBoundaries>();
+
+      for (const auto & [id, pb] : *other_mesh._disjoint_boundary_pairs)
+        if (pb)
+          (*_disjoint_boundary_pairs)[id] = pb->clone();
+    }
+#endif
+
   // _elemset_codes stores pointers to entries in _elemset_codes_inverse_map,
   // so it is not possible to simply copy it directly from other_mesh
   for (const auto & [set, code] : _elemset_codes_inverse_map)
@@ -199,6 +214,20 @@ MeshBase& MeshBase::operator= (MeshBase && other_mesh)
   _node_integer_default_values = std::move(other_mesh._node_integer_default_values);
   _point_locator_close_to_point_tol = other_mesh.get_point_locator_close_to_point_tol();
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  // Deep copy of all periodic boundaries:
+  // We must clone each PeriodicBoundaryBase in the source map,
+  // since unique_ptr cannot be copied and we need independent instances
+  if (other_mesh._disjoint_boundary_pairs)
+    {
+      _disjoint_boundary_pairs = std::make_unique<PeriodicBoundaries>();
+
+      for (const auto & [id, pb] : *other_mesh._disjoint_boundary_pairs)
+        if (pb)
+          (*_disjoint_boundary_pairs)[id] = pb->clone();
+    }
+#endif
+
   // This relies on our subclasses *not* invalidating pointers when we
   // do their portion of the move assignment later!
   boundary_info = std::move(other_mesh.boundary_info);
@@ -306,6 +335,14 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const
   if (*boundary_info != *other_mesh.boundary_info)
     return false;
 
+  // First check whether the "existence" of the two pointers differs (one present, one absent)
+  if ((bool)_disjoint_boundary_pairs != (bool)other_mesh._disjoint_boundary_pairs)
+    return false;
+  // If both exist, compare the contents (Weak Test: just compare sizes like `_ghosting_functors`)
+  if (_disjoint_boundary_pairs &&
+      (_disjoint_boundary_pairs->size() != other_mesh._disjoint_boundary_pairs->size()))
+    return false;
+
   const constraint_rows_type & other_rows =
     other_mesh.get_constraint_rows();
   for (const auto & [node, row] : this->_constraint_rows)
@@ -1953,6 +1990,78 @@ void MeshBase::detect_interior_parents()
 
 
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  /**
+   * Register a pair of boundaries as disjoint neighbor boundary pairs.
+   */
+  void MeshBase::add_disjoint_neighbor_boundary_pairs(const boundary_id_type b1,
+                                             const boundary_id_type b2,
+                                             const RealVectorValue & translation)
+    {
+      // Lazily allocate the container the first time it’s needed
+      if (!_disjoint_boundary_pairs)
+        _disjoint_boundary_pairs = std::make_unique<PeriodicBoundaries>();
+
+      PeriodicBoundaries & db = *_disjoint_boundary_pairs;
+
+      // Create forward and inverse boundary mappings
+      PeriodicBoundary forward(translation);
+      PeriodicBoundary inverse(translation * -1.0);
+
+      forward.myboundary       = b1;
+      forward.pairedboundary   = b2;
+      inverse.myboundary       = b2;
+      inverse.pairedboundary   = b1;
+
+      // Add both directions into the container
+      db.emplace(b1, forward.clone());
+      db.emplace(b2, inverse.clone());
+    }
+
+  PeriodicBoundaries * MeshBase::get_disjoint_neighbor_boundary_pairs()
+    {
+      return _disjoint_boundary_pairs.get();
+    }
+
+  const PeriodicBoundaries * MeshBase::get_disjoint_neighbor_boundary_pairs() const
+    {
+      return _disjoint_boundary_pairs.get();
+    }
+
+    void MeshBase::remove_disjoint_boundary_pair(const boundary_id_type b1,
+                                                   const boundary_id_type b2)
+    {
+      // Nothing to remove if not allocated or empty
+      if (!_disjoint_boundary_pairs || _disjoint_boundary_pairs->empty())
+        return;
+
+      auto & pairs = *_disjoint_boundary_pairs;
+
+      // Helper to check and erase both directions
+      auto erase_if_match = [](boundary_id_type key,
+                              boundary_id_type pair,
+                              PeriodicBoundaries & pb_map)
+        {
+          auto it = pb_map.find(key);
+          if (it != pb_map.end())
+            {
+              const auto & pb = *(it->second);
+              // Check both directions
+              if ((pb.myboundary == key && pb.pairedboundary == pair) ||
+                  (pb.pairedboundary == key && pb.myboundary == pair))
+                pb_map.erase(it);
+            }
+        };
+
+      erase_if_match(b1, b2, pairs);
+      erase_if_match(b2, b1, pairs);
+    }
+
+
+#endif
+
+
+
 void MeshBase::set_point_locator_close_to_point_tol(Real val)
 {
   _point_locator_close_to_point_tol = val;

src/mesh/mesh_refinement.C

@@ -684,6 +684,22 @@ bool MeshRefinement::refine_elements ()
   // This function must be run on all processors at once
   parallel_object_only();
 
+  // Prevent refinement when the mesh has disjoint neighbor boundary pairs.
+  //
+  // Disjoint boundary interfaces violate the topological continuity
+  // assumed by the refinement algorithm. Refinement on such meshes can
+  // produce invalid neighbor relationships (for example reverse indices
+  // equal to invalid_uint), which may trigger assertions like
+  // `rev < neigh->n_neighbors()` or corrupt neighbor data.
+  //
+  // For safety, refinement is disabled while disjoint neighbor boundary pairs are
+  // present.
+  if (_mesh.get_disjoint_neighbor_boundary_pairs())
+    {
+      libmesh_not_implemented_msg(
+        "Mesh refinement is not yet implemented for meshes with disjoint boundary interfaces.\n");
+    }
+
   if (_face_level_mismatch_limit)
     libmesh_assert(test_level_one(true));