From d9b13e40365dc031625e18814514f3502ef69219 Mon Sep 17 00:00:00 2001 From: ttruster Date: Thu, 13 Apr 2023 21:58:33 -0400 Subject: [PATCH 1/4] Add enum for strong and weak enforcement --- include/base/periodic_boundary_base.h | 19 +++++++++++++++++++ src/base/periodic_boundary_base.C | 12 ++++++++++++ src/fe/fe_base.C | 2 +- 3 files changed, 32 insertions(+), 1 deletion(-) 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/src/base/periodic_boundary_base.C b/src/base/periodic_boundary_base.C index f6d06eb8573..9154fbc576e 100644 --- a/src/base/periodic_boundary_base.C +++ b/src/base/periodic_boundary_base.C @@ -111,6 +111,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); From 30efee2939a380e0c06b5b9be168fe00c101abc1 Mon Sep 17 00:00:00 2001 From: ttruster Date: Thu, 20 Apr 2023 14:30:30 -0400 Subject: [PATCH 2/4] Add topological neighbor method with return of neighbor side --- include/geom/elem.h | 12 ++++++++++++ src/geom/elem.C | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/include/geom/elem.h b/include/geom/elem.h index c6f3a70a039..1a2cb8a4d21 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -335,6 +335,18 @@ class Elem : public ReferenceCountedObject, MeshBase & mesh, 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 Elem * topological_neighbor_side(const unsigned int i, + const MeshBase & mesh, + const PointLocatorBase & point_locator, + const PeriodicBoundaries * pb, + unsigned int * neigh_side) const; /** * \returns \p true if the element \p elem in question is a neighbor or diff --git a/src/geom/elem.C b/src/geom/elem.C index 0018b8ed37b..f5bb514fb03 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -1153,6 +1153,45 @@ const Elem * Elem::topological_neighbor (const unsigned int i, } +const Elem * +Elem::topological_neighbor_side(const unsigned int i, + const MeshBase & mesh, + const PointLocatorBase & point_locator, + const PeriodicBoundaries * pb, + unsigned int * neigh_side) const +{ + libmesh_assert_less (i, this->n_neighbors()); + + const Elem * neighbor_i = this->neighbor_ptr(i); + if (neighbor_i != nullptr) + return neighbor_i; + + 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)) + { + 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(); + return neighbor_i; + } + } + + return nullptr; +} + + bool Elem::has_topological_neighbor (const Elem * elem, const MeshBase & mesh, const PointLocatorBase & point_locator, From e8c737aa9f7ba066da7194caeb0ba3699bbdad47 Mon Sep 17 00:00:00 2001 From: ttruster Date: Thu, 20 Apr 2023 14:49:14 -0400 Subject: [PATCH 3/4] remove whitespace --- include/geom/elem.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/geom/elem.h b/include/geom/elem.h index 1a2cb8a4d21..513ccf900a6 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -335,7 +335,7 @@ class Elem : public ReferenceCountedObject, MeshBase & mesh, 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 From 481e525216d68c5675ab48a04dd1b5dc9a3ebef4 Mon Sep 17 00:00:00 2001 From: ttruster Date: Thu, 20 Apr 2023 22:46:45 -0400 Subject: [PATCH 4/4] Add constructor for default of enforcement type Return a pair for the topological neighbor --- include/geom/elem.h | 10 +++++----- src/base/periodic_boundary_base.C | 6 ++++-- src/geom/elem.C | 20 +++++++++++++------- 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/include/geom/elem.h b/include/geom/elem.h index 513ccf900a6..1f3b9c190b8 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -342,11 +342,11 @@ class Elem : public ReferenceCountedObject, * boundary, it will return a corresponding element on the opposite * side, along with the side neigh_side. */ - const Elem * topological_neighbor_side(const unsigned int i, - const MeshBase & mesh, - const PointLocatorBase & point_locator, - const PeriodicBoundaries * pb, - unsigned int * neigh_side) const; + 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 diff --git a/src/base/periodic_boundary_base.C b/src/base/periodic_boundary_base.C index 9154fbc576e..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) diff --git a/src/geom/elem.C b/src/geom/elem.C index f5bb514fb03..571c9510669 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -1153,18 +1153,21 @@ const Elem * Elem::topological_neighbor (const unsigned int i, } -const Elem * +const std::pair Elem::topological_neighbor_side(const unsigned int i, const MeshBase & mesh, const PointLocatorBase & point_locator, - const PeriodicBoundaries * pb, - unsigned int * neigh_side) const + const PeriodicBoundaries * pb) const { libmesh_assert_less (i, this->n_neighbors()); const Elem * neighbor_i = this->neighbor_ptr(i); if (neighbor_i != nullptr) - return neighbor_i; + { + 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) { @@ -1176,6 +1179,8 @@ Elem::topological_neighbor_side(const unsigned int i, 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 @@ -1184,11 +1189,12 @@ Elem::topological_neighbor_side(const unsigned int i, if (neighbor_i) while (level() < neighbor_i->level()) neighbor_i = neighbor_i->parent(); - return neighbor_i; + const std::pair neighbor_and_side (neighbor_i, neighbor_side); + return neighbor_and_side; } } - - return nullptr; + const std::pair neighbor_and_side (nullptr, 0); + return neighbor_and_side; }