Integrate disconnected neighbor detection via PeriodicBoundaries

- Introduced PeriodicBoundaries::disconnected_neighbor() to locate geometrically coincident elements across disconnected boundaries using PointLocator.
- Replaced legacy geometry-based neighbor search in UnstructuredMesh::find_neighbors() with unified PeriodicBoundaries-based logic.

Author: ChengHauYang

.DS_Store

@@ -68,7 +68,8 @@ class PeriodicBoundaries : public std::map<boundary_id_type, std::unique_ptr<Per
                         const PointLocatorBase & point_locator,
                         const Elem * e,
                         unsigned int side,
-                        unsigned int * neigh_side = nullptr) const;
+                        unsigned int * neigh_side = nullptr,
+                        bool skip_finding_check = false) const;
 };
 
 } // namespace libMesh

include/base/periodic_boundaries.h

@@ -36,6 +36,10 @@
 #include <string>
 #include <memory>
 
+// periodic boundary condition support
+#include "libmesh/periodic_boundaries.h"
+#include "libmesh/periodic_boundary.h"
+
 namespace libMesh
 {
 
@@ -1820,27 +1824,49 @@ 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 disconnected boundaries.
    */
-  void add_disconnected_boundaries (const boundary_id_type b1,
-                                    const boundary_id_type b2)
-  {
-    if (b1 == b2)
-      _boundary_id_pairs.emplace(b1, b2);
-    else
+  void add_disconnected_boundaries(const boundary_id_type b1,
+                                  const boundary_id_type b2)
     {
-      _boundary_id_pairs.emplace(b1, b2);
-      _boundary_id_pairs.emplace(b2, b1);
+      // Lazily allocate the container the first time it’s needed
+      if (!_disconnected_boundary_pairs)
+        _disconnected_boundary_pairs = std::make_unique<PeriodicBoundaries>();
+
+      // Create forward and inverse boundary mappings
+      PeriodicBoundary forward(RealVectorValue(0., 0., 0.));
+      PeriodicBoundary inverse(RealVectorValue(0., 0., 0.));
+
+      forward.myboundary       = b1;
+      forward.pairedboundary   = b2;
+      inverse.myboundary       = b2;
+      inverse.pairedboundary   = b1;
+
+      // Add both directions into the container
+      _disconnected_boundary_pairs->emplace(b1, forward.clone());
+      _disconnected_boundary_pairs->emplace(b2, inverse.clone());
+    }
+
+  PeriodicBoundaries * get_disconnected_boundaries()
+    {
+      return _disconnected_boundary_pairs.get();
     }
-  }
 
+  const PeriodicBoundaries * get_disconnected_boundaries() const
+    {
+      return _disconnected_boundary_pairs.get();
+    }
+#endif
 
 
 protected:
 
-  /// @brief a map of pairs of boundary ids between which new boundary sides are created
-  std::unordered_map<boundary_id_type, boundary_id_type> _boundary_id_pairs;
+ #ifdef LIBMESH_ENABLE_PERIODIC
+  /// @brief The disconnected boundary id pairs.
+  std::unique_ptr<PeriodicBoundaries> _disconnected_boundary_pairs;
+#endif
 
   /**
    * This class holds the boundary information.  It can store nodes, edges,

include/mesh/mesh_base.h

@@ -60,7 +60,8 @@ const Elem * PeriodicBoundaries::neighbor(boundary_id_type boundary_id,
                                           const PointLocatorBase & point_locator,
                                           const Elem * e,
                                           unsigned int side,
-                                          unsigned int * neigh_side) const
+                                          unsigned int * neigh_side,
+                                          bool skip_finding_check) const
 {
   std::unique_ptr<const Elem> neigh_side_proxy;
 
@@ -80,6 +81,10 @@ const Elem * PeriodicBoundaries::neighbor(boundary_id_type boundary_id,
   const MeshBase & mesh = point_locator.get_mesh();
   for(const Elem * elem_it : candidate_elements)
     {
+
+      if (elem_it == e) // skip self
+        continue;
+
       std::vector<unsigned int> neigh_sides =
         mesh.get_boundary_info().sides_with_boundary_id(elem_it, b->pairedboundary);
 
@@ -101,18 +106,18 @@ const Elem * PeriodicBoundaries::neighbor(boundary_id_type boundary_id,
         }
     }
 
-  // If we should have found a periodic neighbor but didn't then
-  // either we're on a ghosted element with a remote periodic neighbor
-  // or we're on a mesh with an inconsistent periodic boundary.
-  libmesh_error_msg_if(mesh.is_serial() ||
-                       (e->processor_id() == mesh.processor_id()),
-                       "Periodic boundary neighbor not found");
+  if (!skip_finding_check)
+    // If we should have found a periodic neighbor but didn't then
+    // either we're on a ghosted element with a remote periodic neighbor
+    // or we're on a mesh with an inconsistent periodic boundary.
+    libmesh_error_msg_if(mesh.is_serial() ||
+                        (e->processor_id() == mesh.processor_id()),
+                        "Periodic boundary neighbor not found");
 
   if (neigh_side)
     *neigh_side = libMesh::invalid_uint;
   return remote_elem;
 }
-
 } // namespace libMesh
 
 

src/base/periodic_boundaries.C

@@ -833,41 +833,6 @@ bool Elem::topologically_equal (const Elem & rhs) const
 
 
 
-bool Elem::geometrically_equal(const Elem &rhs, const Real tol_in) const
-{
-  // 1. Basic checks: dimension & node count
-  if (this->dim() != rhs.dim())
-    return false;
-  if (this->n_nodes() != rhs.n_nodes())
-    return false;
-
-  // 2. Compare centroids within tolerance
-  const Real tol = tol_in > 0 ? tol_in : 1e-12;
-  const Point c1 = this->vertex_average();
-  const Point c2 = rhs.vertex_average();
-  if ((c1 - c2).norm() > tol)
-    return false;
-
-  // 3. Unordered comparison within tolerance
-  const unsigned int n = this->n_nodes();
-  std::vector<Point> pts1(n), pts2(n);
-  for (unsigned int i = 0; i < n; ++i)
-    pts1[i] = this->point(i);
-  for (unsigned int i = 0; i < n; ++i)
-    pts2[i] = rhs.point(i);
-
-  for (const auto &p1 : pts1)
-    {
-      auto found = std::any_of(pts2.begin(), pts2.end(),
-                              [&](const Point &p2)
-                              { return (p1 - p2).norm() <= tol; });
-      if (!found)
-        return false;
-    }
-
-  return true;
-}
-
 bool Elem::is_semilocal(const processor_id_type my_pid) const
 {
   std::set<const Elem *> point_neighbors;

src/geom/elem.C

@@ -47,9 +47,6 @@
 #include <sstream>
 #include <unordered_map>
 
-// TIMPI includes
-#include "timpi/parallel_implementation.h"
-#include "timpi/parallel_sync.h"
 namespace {
 
 using namespace libMesh;
@@ -996,111 +993,40 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
       }
   }
 
-{
-  // Build or obtain a point locator (sub-locator) and initialize it
-  std::unique_ptr<PointLocatorBase> point_locator = this->sub_point_locator();
-
-  // Get a reference to the boundary info object for convenience
-  auto & bi = this->get_boundary_info();
+#ifdef LIBMESH_ENABLE_PERIODIC
+  auto * db = this->get_disconnected_boundaries();
 
-  // Check whether the element side belongs to a paired boundary
-  auto pick_paired_id = [this, &bi](const Elem * e, unsigned int side,
-                                    boundary_id_type &current_id,
-                                    boundary_id_type &paired_id) -> bool
-  {
-    std::vector<boundary_id_type> ids;
-    bi.boundary_ids(e, side, ids);
-    for (auto id : ids)
+  if (db)
     {
-      auto it = _boundary_id_pairs.find(id);
-      if (it != _boundary_id_pairs.end())
-      {
-        current_id = id;
-        paired_id  = it->second;
-        return true;
-      }
-    }
-    return false;
-  };
-
-  // Main loop: iterate over all elements and their sides to find disconnected neighbors
-  for (const auto & element : this->element_ptr_range())
-  {
-    for (auto ms : element->side_index_range())
-    {
-      // Skip if this side already has a valid neighbor (including remote neighbors)
-      if (element->neighbor_ptr(ms) != nullptr &&
-          element->neighbor_ptr(ms) != remote_elem)
-        continue;
-
-      // Determine whether this side belongs to a paired (disconnected) boundary
-      boundary_id_type current_id;
-      boundary_id_type paired_id;
-      if (!pick_paired_id(element, ms, current_id, paired_id))
-        continue; // No valid pair found, skip this side
-
-      // Build the geometry of this side
-      std::unique_ptr<const Elem> this_side;
-      element->build_side_ptr(this_side, ms);
-
-      // Compute centroid and approximate element size
-      const Point c_mid  = this_side->vertex_average();
-
-      // Use the point locator to find candidate elements near the side centroid
-      std::set<const Elem *> const_candidate_elements;
-      point_locator->operator()(c_mid, const_candidate_elements);
+      // Obtain a point locator
+      std::unique_ptr<PointLocatorBase> point_locator = this->sub_point_locator();
 
-      std::set<Elem *> candidate_elements;
-      for (auto ce : const_candidate_elements)
-          candidate_elements.insert(this->elem_ptr(ce->id()));
+      // Get the disconnected boundaries object (from periodic BCs)
+      auto * db = this->get_disconnected_boundaries();
 
-      bool paired_found = false; // Flag to mark if a valid match was found
-
-      // Loop through candidate elements to look for potential paired neighbors
-      for (auto * cand_elem : candidate_elements)
+      for (const auto & element : this->element_ptr_range())
         {
-          if (!cand_elem || cand_elem == element)
-            continue;
-
-          // If paired_id == current_id, it means that there is no paired boundary on the other side element.
-          // In this case, we should consider all sides of the candidate element.
-          std::vector<unsigned int> cand_neigh_sides =
-            (paired_id == current_id) ?
-              std::vector<unsigned int>(cand_elem->side_index_range().begin(),
-                                    cand_elem->side_index_range().end()) :
-              bi.sides_with_boundary_id(cand_elem, paired_id);
-
-          if (cand_neigh_sides.empty())
-            continue;
-
-          // Compare each candidate side to the current one
-          for (auto ns : cand_neigh_sides)
+          for (auto ms : element->side_index_range())
             {
-              // Skip if this side already has a valid (non-remote) neighbor
-              if (cand_elem->neighbor_ptr(ns) != nullptr &&
-                  cand_elem->neighbor_ptr(ns) != remote_elem)
+              // Skip if this side already has a valid neighbor (including remote neighbors)
+              if (element->neighbor_ptr(ms) != nullptr &&
+                  element->neighbor_ptr(ms) != remote_elem)
                 continue;
 
-              // Build the geometry of the candidate side
-              std::unique_ptr<const Elem> cand_side;
-              cand_elem->build_side_ptr(cand_side, ns);
-
-              // Check geometric equality
-              if (this_side->geometrically_equal(*cand_side))
+                for (const auto & [id, boundary_ptr] : *db)
                 {
-                  element->set_neighbor(ms, cand_elem);
-                  cand_elem->set_neighbor(ns, element);
-                  paired_found = true;
+                  unsigned int neigh_side;
+                  const Elem * neigh = db->neighbor(id, *point_locator, element, ms, &neigh_side, true /*skip_finding_check*/);
+                  if (neigh && neigh != remote_elem && neigh != element)
+                    {
+                      element->set_neighbor(ms, this->elem_ptr(neigh->id()));
+                      this->elem_ptr(neigh->id())->set_neighbor(neigh_side, element);
+                    }
                 }
             }
-
-          if (paired_found)
-            break; // Stop searching through other candidates
         }
     }
-  }
-}
-
+#endif // LIBMESH_ENABLE_PERIODIC
 
 #ifdef LIBMESH_ENABLE_AMR
 
@@ -1325,6 +1251,7 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
 
 #endif // AMR
 
+
 #ifdef DEBUG
   MeshTools::libmesh_assert_valid_neighbors(*this,
                                             !reset_remote_elements);