diff --git a/include/base/periodic_boundary_base.h b/include/base/periodic_boundary_base.h index 30b8cfcc87c..01353f0f838 100644 --- a/include/base/periodic_boundary_base.h +++ b/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 & 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> _transformation_matrix; + EnforcementType _enforcement_type; + }; } // namespace libmesh diff --git a/include/geom/elem.h b/include/geom/elem.h index c6f3a70a039..1f3b9c190b8 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -336,6 +336,18 @@ class Elem : public ReferenceCountedObject, 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 + 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. diff --git a/src/base/periodic_boundary_base.C b/src/base/periodic_boundary_base.C index f6d06eb8573..85d9cebcdd4 100644 --- a/src/base/periodic_boundary_base.C +++ b/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 & 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 diff --git a/src/fe/fe_base.C b/src/fe/fe_base.C index 799cf38306f..9abd6e21cf7 100644 --- a/src/fe/fe_base.C +++ b/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); diff --git a/src/geom/elem.C b/src/geom/elem.C index 0018b8ed37b..571c9510669 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -1153,6 +1153,51 @@ const Elem * Elem::topological_neighbor (const unsigned int i, } +const std::pair +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 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 bc_ids; + mesh.get_boundary_info().boundary_ids(this, cast_int(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 neighbor_and_side (neighbor_i, neighbor_side); + return neighbor_and_side; + } + } + const std::pair neighbor_and_side (nullptr, 0); + return neighbor_and_side; +} + + bool Elem::has_topological_neighbor (const Elem * elem, const MeshBase & mesh, const PointLocatorBase & point_locator,