Merge pull request #2725 from lindsayad/ghost-periodic-point-neighbors

Ghost periodic point neighbors in GhostPointNeighbors

Author: lindsayad

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,27 @@ 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
+        ,
+      // We do not simply want to copy over the other's periodic bcs because
+      // they may very well correspond to periodic bcs from a \p DofMap entirely
+      // unrelated to any \p DofMaps that depend on this objects ghosting
+      _periodic_bcs(nullptr)
+#endif
+    {}
 
   /**
    * A clone() is needed because GhostingFunctor can not be shared between
@@ -56,6 +74,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 +101,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,20 @@ 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)
+      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> appn_bcids;
+  std::vector<const Elem *> active_periodic_neighbors;
+#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.  ;-)
@@ -90,7 +110,87 @@ void GhostPointNeighbors::operator()
       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)
+        {
+          for (const auto s : elem->side_index_range())
+            {
+              if (elem->neighbor_ptr(s))
+                continue;
+
+              const Elem * const equal_level_periodic_neigh = elem->topological_neighbor
+                (s, *_mesh, *point_locator, _periodic_bcs);
+
+              if (!equal_level_periodic_neigh || equal_level_periodic_neigh == remote_elem)
+                continue;
+
+              equal_level_periodic_neigh->active_family_tree_by_topological_neighbor(
+                active_periodic_neighbors,
+                elem,
+                *_mesh,
+                *point_locator,
+                _periodic_bcs,
+                /*reset=*/true);
+
+              for (const Elem * const active_periodic_neigh : active_periodic_neighbors)
+                {
+                  std::set <const Elem *> active_periodic_point_neighbors;
+
+                  // This fills point neighbors *including*
+                  // active_periodic_neigh. The documentation for this method
+                  // states that this will return *active* point neighbors
+                  active_periodic_neigh->find_point_neighbors(active_periodic_point_neighbors);
+
+                  for (const Elem * const appn : active_periodic_point_neighbors)
+                    {
+                      // Don't need to ghost RemoteElem or an element we already own or an
+                      // element we've already examined
+                      if (appn == remote_elem || appn->processor_id() == p ||
+                          periodic_elems_examined.count(appn))
+                        continue;
+
+                      // We only need to keep point neighbors that are along the periodic boundaries
+                      bool on_periodic_boundary = false;
+                      for (const auto appn_s : appn->side_index_range())
+                        {
+                          binfo.boundary_ids(appn, appn_s, appn_bcids);
+                          for (const auto appn_bcid : appn_bcids)
+                            if (_periodic_bcs->find(appn_bcid) != _periodic_bcs->end())
+                              {
+                                on_periodic_boundary = true;
+                                goto jump;
+                              }
+                        }
+                      jump:
+                        if (on_periodic_boundary)
+                          coupled_elements.emplace(appn, nullcm);
+
+                        periodic_elems_examined.insert(appn);
+                    }
+                }
+            }
+        }
+#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

src/mesh/mesh_base.C

@@ -104,39 +104,39 @@ MeshBase::MeshBase (const MeshBase & other_mesh) :
   _default_ghosting(libmesh_make_unique<GhostPointNeighbors>(*this)),
   _point_locator_close_to_point_tol(other_mesh._point_locator_close_to_point_tol)
 {
-   for (const auto & gf : other_mesh._ghosting_functors )
-   {
-     std::shared_ptr<GhostingFunctor> clone_gf = gf->clone();
-     // Some subclasses of GhostingFunctor might not override the
-     // clone function yet. If this is the case, GhostingFunctor will
-     // return nullptr by default. The clone function should be overridden
-     // in all derived classes. This following code ("else") is written
-     // for API upgrade. That will allow users gradually to update their code.
-     // Once the API upgrade is done, we will come back and delete "else."
-     if (clone_gf)
-     {
-       clone_gf->set_mesh(this);
-       add_ghosting_functor(clone_gf);
-     }
-     else
-     {
-       libmesh_deprecated();
-       add_ghosting_functor(*gf);
-     }
-   }
-
-  // Make sure we don't accidentally delete the other mesh's default
-  // ghosting functor; we'll use our own if that's needed.
-  if (other_mesh._ghosting_functors.count(other_mesh._default_ghosting.get()))
+  const GhostingFunctor * const other_default_ghosting = other_mesh._default_ghosting.get();
+
+  for (GhostingFunctor * const gf : other_mesh._ghosting_functors)
     {
-      _ghosting_functors.erase(other_mesh._default_ghosting.get());
-      _ghosting_functors.insert(_default_ghosting.get());
+      // If the other mesh is using default ghosting, then we will use our own
+      // default ghosting
+      if (gf == other_default_ghosting)
+        {
+          _ghosting_functors.insert(_default_ghosting.get());
+          continue;
+        }
+
+      std::shared_ptr<GhostingFunctor> clone_gf = gf->clone();
+      // Some subclasses of GhostingFunctor might not override the
+      // clone function yet. If this is the case, GhostingFunctor will
+      // return nullptr by default. The clone function should be overridden
+      // in all derived classes. This following code ("else") is written
+      // for API upgrade. That will allow users gradually to update their code.
+      // Once the API upgrade is done, we will come back and delete "else."
+      if (clone_gf)
+        {
+          clone_gf->set_mesh(this);
+          add_ghosting_functor(clone_gf);
+        }
+      else
+        {
+          libmesh_deprecated();
+          add_ghosting_functor(*gf);
+        }
     }
 
   if (other_mesh._partitioner.get())
-    {
-      _partitioner = other_mesh._partitioner->clone();
-    }
+    _partitioner = other_mesh._partitioner->clone();
 }