Ghost periodic point neighbors in GhostPointNeighbors

This fixes my distributed periodic boundary issue. Currently our
periodic constraint algorithm relies on identifying a primary element in
charge of constraining (vertex) dofs. That designation is determined by
looking at point neighbors and comparing things like ids. If not all the
point neighbors are available, then a process may make an invalid
determination that it owns a primary elem when in fact the primary elem
is off-process. In order to to validly determine primary elems, we must
have all the periodic point neighbors available

include/base/default_coupling.h

@@ -86,7 +86,7 @@ class DefaultCoupling : public GhostingFunctor
 
 #ifdef LIBMESH_ENABLE_PERIODIC
   // Set PeriodicBoundaries to couple
-  void set_periodic_boundaries(const PeriodicBoundaries * periodic_bcs)
+  void set_periodic_boundaries(const PeriodicBoundaries * periodic_bcs) override
   { _periodic_bcs = periodic_bcs; }
 #endif
 

include/base/ghost_point_neighbors.h

@@ -26,6 +26,9 @@
 
 namespace libMesh
 {
+#ifdef LIBMESH_ENABLE_PERIODIC
+class PeriodicBoundaries;
+#endif
 
 /**
  * This class implements the original default geometry ghosting
@@ -42,12 +45,24 @@ class GhostPointNeighbors : public GhostingFunctor
   /**
    * Constructor.
    */
-  GhostPointNeighbors(const MeshBase & mesh) : GhostingFunctor(mesh) {}
+  GhostPointNeighbors(const MeshBase & mesh) :
+      GhostingFunctor(mesh)
+#ifdef LIBMESH_ENABLE_PERIODIC
+      ,
+      _periodic_bcs(nullptr)
+#endif
+    {}
 
   /**
    * Constructor.
    */
-  GhostPointNeighbors(const GhostPointNeighbors & other) : GhostingFunctor(other){}
+  GhostPointNeighbors(const GhostPointNeighbors & other) :
+      GhostingFunctor(other)
+#ifdef LIBMESH_ENABLE_PERIODIC
+        ,
+      _periodic_bcs(other._periodic_bcs)
+#endif
+    {}
 
   /**
    * A clone() is needed because GhostingFunctor can not be shared between
@@ -56,6 +71,24 @@ class GhostPointNeighbors : public GhostingFunctor
   virtual std::unique_ptr<GhostingFunctor> clone () const override
   { return libmesh_make_unique<GhostPointNeighbors>(*this); }
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  // Set PeriodicBoundaries to couple
+  void set_periodic_boundaries(const PeriodicBoundaries * periodic_bcs) override
+  { _periodic_bcs = periodic_bcs; }
+#endif
+
+  /**
+   * If we have periodic boundaries, then we'll need the mesh to have
+   * an updated point locator whenever we're about to query them.
+   */
+  virtual void mesh_reinit () override;
+
+  virtual void redistribute () override
+  { this->mesh_reinit(); }
+
+  virtual void delete_remote_elements() override
+  { this->mesh_reinit(); }
+
   /**
    * For the specified range of active elements, find their point
    * neighbors and interior_parent elements, ignoring those on
@@ -65,6 +98,11 @@ class GhostPointNeighbors : public GhostingFunctor
                            const MeshBase::const_element_iterator & range_end,
                            processor_id_type p,
                            map_type & coupled_elements) override;
+
+private:
+#ifdef LIBMESH_ENABLE_PERIODIC
+  const PeriodicBoundaries * _periodic_bcs;
+#endif
 };
 
 } // namespace libMesh

include/base/ghosting_functor.h

@@ -36,7 +36,9 @@ namespace libMesh
 // Forward Declarations
 class CouplingMatrix;
 class Elem;
-
+#ifdef LIBMESH_ENABLE_PERIODIC
+class PeriodicBoundaries;
+#endif
 
 /**
  * This abstract base class defines the interface by which library
@@ -194,6 +196,11 @@ class GhostingFunctor : public ReferenceCountedObject<GhostingFunctor>
    */
   virtual void set_mesh(const MeshBase * mesh) { _mesh = mesh; }
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  // Set PeriodicBoundaries to couple
+  virtual void set_periodic_boundaries(const PeriodicBoundaries *) {}
+#endif
+
   /**
    *  Return the mesh associated with ghosting functor
    */

include/base/point_neighbor_coupling.h

@@ -88,7 +88,7 @@ class PointNeighborCoupling : public GhostingFunctor
   // Set PeriodicBoundaries to couple.
   //
   // FIXME: This capability is not currently implemented.
-  void set_periodic_boundaries(const PeriodicBoundaries * periodic_bcs)
+  void set_periodic_boundaries(const PeriodicBoundaries * periodic_bcs) override
   { _periodic_bcs = periodic_bcs; }
 #endif
 

src/base/ghost_point_neighbors.C

@@ -21,9 +21,15 @@
 
 #include "libmesh/elem.h"
 #include "libmesh/remote_elem.h"
+#ifdef LIBMESH_ENABLE_PERIODIC
+#include "libmesh/periodic_boundaries.h"
+#include "libmesh/boundary_info.h"
+#endif
 
 // C++ Includes
 #include <unordered_set>
+#include <set>
+#include <vector>
 
 namespace libMesh
 {
@@ -36,6 +42,22 @@ void GhostPointNeighbors::operator()
 {
   libmesh_assert(_mesh);
 
+#ifdef LIBMESH_ENABLE_PERIODIC
+  bool check_periodic_bcs =
+    (_periodic_bcs && !_periodic_bcs->empty());
+
+  std::unique_ptr<PointLocatorBase> point_locator;
+  if (check_periodic_bcs)
+    {
+      libmesh_assert(_mesh);
+      point_locator = _mesh->sub_point_locator();
+    }
+
+  std::set<const Elem *> periodic_elems_examined;
+  const BoundaryInfo & binfo = _mesh->get_boundary_info();
+  std::vector<boundary_id_type> ppn_bcids;
+#endif
+
   // Using a connected_nodes set rather than point_neighbors() would
   // give us correct results even in corner cases, such as where two
   // elements meet only at a corner.  ;-)
@@ -67,30 +89,97 @@ void GhostPointNeighbors::operator()
   CouplingMatrix * nullcm = nullptr;
 
   for (const auto & elem : as_range(range_begin, range_end))
+  {
+    libmesh_assert(_mesh->query_elem_ptr(elem->id()) == elem);
+
+    if (elem->processor_id() != p)
+      coupled_elements.emplace(elem, nullcm);
+
+    std::set<const Elem *> elem_point_neighbors;
+    elem->find_point_neighbors(elem_point_neighbors);
+
+    for (const auto & neigh : elem_point_neighbors)
+      coupled_elements.emplace(neigh, nullcm);
+
+    // An interior_parent isn't on the same manifold so won't be
+    // found as a point neighbor, and it may not share nodes so we
+    // can't use a connected_nodes test.
+    //
+    // Trying to preserve interior_parent() only works if it's on
+    // the same Mesh, which is *not* guaranteed!  So we'll
+    // double-check.
+    const Elem * ip = elem->interior_parent();
+    if (ip && ip->processor_id() != p &&
+        _mesh->query_elem_ptr(ip->id()) == ip)
+      coupled_elements.emplace(ip, nullcm);
+
+#ifdef LIBMESH_ENABLE_PERIODIC
+    if (check_periodic_bcs)
     {
-      libmesh_assert(_mesh->query_elem_ptr(elem->id()) == elem);
-
-      if (elem->processor_id() != p)
-        coupled_elements.emplace(elem, nullcm);
-
-      std::set<const Elem *> elem_point_neighbors;
-      elem->find_point_neighbors(elem_point_neighbors);
-
-      for (const auto & neigh : elem_point_neighbors)
-        coupled_elements.emplace(neigh, nullcm);
-
-      // An interior_parent isn't on the same manifold so won't be
-      // found as a point neighbor, and it may not share nodes so we
-      // can't use a connected_nodes test.
-      //
-      // Trying to preserve interior_parent() only works if it's on
-      // the same Mesh, which is *not* guaranteed!  So we'll
-      // double-check.
-      const Elem * ip = elem->interior_parent();
-      if (ip && ip->processor_id() != p &&
-          _mesh->query_elem_ptr(ip->id()) == ip)
-        coupled_elements.emplace(ip, nullcm);
+      for (const auto s : elem->side_index_range())
+      {
+        if (elem->neighbor_ptr(s))
+          continue;
+
+        const Elem * const periodic_neigh = elem->topological_neighbor
+          (s, *_mesh, *point_locator, _periodic_bcs);
+
+        if (periodic_neigh && periodic_neigh != remote_elem)
+        {
+          std::set <const Elem *> periodic_point_neighbors;
+
+          // This fills point neighbors *including* periodic_neigh
+          periodic_neigh->find_point_neighbors(periodic_point_neighbors);
+
+          for (const Elem * const ppn : periodic_point_neighbors)
+          {
+            // Don't need to ghost RemoteElem or an element we already own or an
+            // element we've already examined
+            if (ppn == remote_elem || ppn->processor_id() == _mesh->processor_id() ||
+                periodic_elems_examined.count(ppn))
+              continue;
+
+            // We only need to keep point neighbors that are along the periodic boundaries
+            bool on_periodic_boundary = false;
+            for (const auto ppn_s : ppn->side_index_range())
+            {
+              binfo.boundary_ids(ppn, ppn_s, ppn_bcids);
+              for (const auto & pb_pair : *_periodic_bcs)
+                if (std::find(ppn_bcids.begin(), ppn_bcids.end(), pb_pair.first) != ppn_bcids.end())
+                {
+                  on_periodic_boundary = true;
+                  goto jump;
+                }
+            }
+          jump:
+            if (on_periodic_boundary)
+              coupled_elements.emplace(ppn, nullcm);
+
+            periodic_elems_examined.insert(ppn);
+          }
+        }
+      }
     }
+#endif // LIBMESH_ENABLE_PERIODIC
+  }
+}
+
+void GhostPointNeighbors::mesh_reinit()
+{
+  // Unless we have periodic boundary conditions, we don't need
+  // anything precomputed.
+#ifdef LIBMESH_ENABLE_PERIODIC
+  if (!_periodic_bcs || _periodic_bcs->empty())
+    return;
+#endif
+
+  // If we do have periodic boundary conditions, we'll need a master
+  // point locator, so we'd better have a mesh to build it on.
+  libmesh_assert(_mesh);
+
+  // Make sure an up-to-date master point locator has been
+  // constructed; we'll need to grab sub-locators soon.
+  _mesh->sub_point_locator();
 }
 
 } // namespace libMesh