Weak enforcement of periodic boundary conditions

include/base/periodic_boundary_base.h

@@ -38,6 +38,13 @@ namespace libMesh
 class Elem;
 class MeshBase;
 
+/**
+ * The type of enforcement of condition; strong(default) = constraint matrix
+ */
+enum class EnforcementType
+  { STRONG_ENFORCEMENT=0,
+    WEAK_ENFORCEMENT=1 };
+
 /**
  * The base class for defining periodic boundaries.
  *
@@ -126,6 +133,16 @@ class PeriodicBoundaryBase
    */
   const std::set<unsigned int> & get_variables() const;
 
+  /**
+   * Get the enforcement type.
+   */
+  const EnforcementType & get_enforcement_type() const;
+
+  /**
+   * Set the enforcement type.
+   */
+  void set_enforcement_type(const EnforcementType & e_type);
+
 protected:
 
   /**
@@ -145,6 +162,8 @@ class PeriodicBoundaryBase
    */
   std::unique_ptr<DenseMatrix<Real>> _transformation_matrix;
 
+  EnforcementType _enforcement_type;
+
 };
 
 } // namespace libmesh

include/geom/elem.h

@@ -336,6 +336,18 @@ class Elem : public ReferenceCountedObject<Elem>,
                                const PointLocatorBase & point_locator,
                                const PeriodicBoundaries * pb);
 
+  /**
+   * \returns A pointer to the \f$ i^{th} \f$ neighbor of this element
+   * for interior elements.  If an element is on a periodic
+   * boundary, it will return a corresponding element on the opposite
+   * side, along with the side neigh_side.
+   */
+  const std::pair<const Elem *, unsigned int>
+  topological_neighbor_side(const unsigned int i,
+                            const MeshBase & mesh,
+                            const PointLocatorBase & point_locator,
+                            const PeriodicBoundaries * pb) const;
+
   /**
    * \returns \p true if the element \p elem in question is a neighbor or
    * topological neighbor of this element, \p false otherwise.

src/base/periodic_boundary_base.C

@@ -33,7 +33,8 @@ namespace libMesh
 
 PeriodicBoundaryBase::PeriodicBoundaryBase() :
   myboundary(BoundaryInfo::invalid_id),
-  pairedboundary(BoundaryInfo::invalid_id)
+  pairedboundary(BoundaryInfo::invalid_id),
+  _enforcement_type(EnforcementType::STRONG_ENFORCEMENT)
 {
 }
 
@@ -42,7 +43,8 @@ PeriodicBoundaryBase::PeriodicBoundaryBase() :
 PeriodicBoundaryBase::PeriodicBoundaryBase(const PeriodicBoundaryBase & o) :
   myboundary(o.myboundary),
   pairedboundary(o.pairedboundary),
-  variables(o.variables)
+  variables(o.variables),
+  _enforcement_type(EnforcementType::STRONG_ENFORCEMENT)
 {
   // Make a deep copy of _transformation_matrix, if it's not null
   if(o._transformation_matrix)
@@ -111,6 +113,18 @@ const std::set<unsigned int> & PeriodicBoundaryBase::get_variables() const
   return variables;
 }
 
+
+const EnforcementType & PeriodicBoundaryBase::get_enforcement_type() const
+{
+  return _enforcement_type;
+}
+
+
+void PeriodicBoundaryBase::set_enforcement_type(const EnforcementType & e_type)
+{
+  this->_enforcement_type = e_type;
+}
+
 } // namespace libMesh
 
 #endif // LIBMESH_ENABLE_PERIODIC

src/fe/fe_base.C

@@ -1928,7 +1928,7 @@ compute_periodic_constraints (DofConstraints & constraints,
       for (const auto & boundary_id : bc_ids)
         {
           const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
-          if (!periodic || !periodic->is_my_variable(variable_number))
+          if (!periodic || !periodic->is_my_variable(variable_number) || (periodic->get_enforcement_type()==EnforcementType::WEAK_ENFORCEMENT))
             continue;
 
           libmesh_assert(point_locator);

src/geom/elem.C

@@ -1153,6 +1153,51 @@ const Elem * Elem::topological_neighbor (const unsigned int i,
 }
 
 
+const std::pair<const Elem *, unsigned int>
+Elem::topological_neighbor_side(const unsigned int i,
+                                const MeshBase & mesh,
+                                const PointLocatorBase & point_locator,
+                                const PeriodicBoundaries * pb) const
+{
+  libmesh_assert_less (i, this->n_neighbors());
+
+  const Elem * neighbor_i = this->neighbor_ptr(i);
+  if (neighbor_i != nullptr)
+  {
+    unsigned int neighbor_side = neighbor_i->which_neighbor_am_i(this);
+    const std::pair<const Elem *, unsigned int> neighbor_and_side (neighbor_i, neighbor_side);
+    return neighbor_and_side;
+  }
+
+  if (pb)
+    {
+      // Since the neighbor is nullptr it must be on a boundary. We need
+      // see if this is a periodic boundary in which case it will have a
+      // topological neighbor
+      std::vector<boundary_id_type> bc_ids;
+      mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
+      for (const auto & id : bc_ids)
+        if (pb->boundary(id))
+          {
+            unsigned int neighbor_side = 0;
+            unsigned int * neigh_side = &neighbor_side;
+            neighbor_i = pb->neighbor(id, point_locator, this, i, neigh_side);
+
+            // Since coarse elements do not have more refined
+            // neighbors we need to make sure that we don't return one
+            // of these types of neighbors.
+            if (neighbor_i)
+              while (level() < neighbor_i->level())
+                neighbor_i = neighbor_i->parent();
+            const std::pair<const Elem *, unsigned int> neighbor_and_side (neighbor_i, neighbor_side);
+            return neighbor_and_side;
+          }
+    }
+  const std::pair<const Elem *, unsigned int> neighbor_and_side (nullptr, 0);
+  return neighbor_and_side;
+}
+
+
 bool Elem::has_topological_neighbor (const Elem * elem,
                                      const MeshBase & mesh,
                                      const PointLocatorBase & point_locator,