From cdd8204e02eb694fc55cfcd8b1ccdb401c1abce8 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 1 Oct 2020 18:34:57 -0700 Subject: [PATCH 1/5] 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 | 2 +- include/base/ghost_point_neighbors.h | 42 ++++++++++++- include/base/ghosting_functor.h | 9 ++- include/base/point_neighbor_coupling.h | 2 +- src/base/ghost_point_neighbors.C | 86 ++++++++++++++++++++++++++ 5 files changed, 136 insertions(+), 5 deletions(-) diff --git a/include/base/default_coupling.h b/include/base/default_coupling.h index 9e3a27ea9d5..17bdda5e3e7 100644 --- a/include/base/default_coupling.h +++ b/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 diff --git a/include/base/ghost_point_neighbors.h b/include/base/ghost_point_neighbors.h index c8ade660ec7..1b520b1d525 100644 --- a/include/base/ghost_point_neighbors.h +++ b/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 clone () const override { return libmesh_make_unique(*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 diff --git a/include/base/ghosting_functor.h b/include/base/ghosting_functor.h index 4c979a23da3..49c57435f6f 100644 --- a/include/base/ghosting_functor.h +++ b/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 */ 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 */ diff --git a/include/base/point_neighbor_coupling.h b/include/base/point_neighbor_coupling.h index cea5f6821de..7d51bf011c8 100644 --- a/include/base/point_neighbor_coupling.h +++ b/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 diff --git a/src/base/ghost_point_neighbors.C b/src/base/ghost_point_neighbors.C index 8d50d3fdd32..bb0959ad696 100644 --- a/src/base/ghost_point_neighbors.C +++ b/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 +#include +#include namespace libMesh { @@ -36,6 +42,19 @@ void GhostPointNeighbors::operator() { libmesh_assert(_mesh); +#ifdef LIBMESH_ENABLE_PERIODIC + bool check_periodic_bcs = + (_periodic_bcs && !_periodic_bcs->empty()); + + std::unique_ptr point_locator; + if (check_periodic_bcs) + point_locator = _mesh->sub_point_locator(); + + std::set periodic_elems_examined; + const BoundaryInfo & binfo = _mesh->get_boundary_info(); + std::vector 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. ;-) @@ -90,7 +109,74 @@ 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 periodic_neigh = elem->topological_neighbor + (s, *_mesh, *point_locator, _periodic_bcs); + + if (periodic_neigh && periodic_neigh != remote_elem) + { + std::set 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 ppn_bcid : ppn_bcids) + if (_periodic_bcs->find(ppn_bcid) != _periodic_bcs->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 From b18512d2f9226f0bc756eb5eeb6d85c4b000962a Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 2 Oct 2020 12:15:20 -0700 Subject: [PATCH 2/5] Don't copy others periodic bcs --- include/base/ghost_point_neighbors.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/base/ghost_point_neighbors.h b/include/base/ghost_point_neighbors.h index 1b520b1d525..c87076106db 100644 --- a/include/base/ghost_point_neighbors.h +++ b/include/base/ghost_point_neighbors.h @@ -60,7 +60,10 @@ class GhostPointNeighbors : public GhostingFunctor GhostingFunctor(other) #ifdef LIBMESH_ENABLE_PERIODIC , - _periodic_bcs(other._periodic_bcs) + // 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 {} From ed92819f63be6cb4022c572a84597766e3c46637 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 2 Oct 2020 12:15:55 -0700 Subject: [PATCH 3/5] Make sure to ghost active point neighbors --- src/base/ghost_point_neighbors.C | 48 +++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/src/base/ghost_point_neighbors.C b/src/base/ghost_point_neighbors.C index bb0959ad696..5654c616283 100644 --- a/src/base/ghost_point_neighbors.C +++ b/src/base/ghost_point_neighbors.C @@ -52,7 +52,8 @@ void GhostPointNeighbors::operator() std::set periodic_elems_examined; const BoundaryInfo & binfo = _mesh->get_boundary_info(); - std::vector ppn_bcids; + std::vector appn_bcids; + std::vector active_periodic_neighbors; #endif // Using a connected_nodes set rather than point_neighbors() would @@ -118,41 +119,54 @@ void GhostPointNeighbors::operator() if (elem->neighbor_ptr(s)) continue; - const Elem * const periodic_neigh = elem->topological_neighbor + const Elem * const equal_level_periodic_neigh = elem->topological_neighbor (s, *_mesh, *point_locator, _periodic_bcs); - if (periodic_neigh && periodic_neigh != remote_elem) + 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 periodic_point_neighbors; + std::set active_periodic_point_neighbors; - // This fills point neighbors *including* periodic_neigh - periodic_neigh->find_point_neighbors(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 ppn : 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 (ppn == remote_elem || ppn->processor_id() == _mesh->processor_id() || - periodic_elems_examined.count(ppn)) + if (appn == remote_elem || appn->processor_id() == _mesh->processor_id() || + 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 ppn_s : ppn->side_index_range()) + for (const auto appn_s : appn->side_index_range()) { - binfo.boundary_ids(ppn, ppn_s, ppn_bcids); - for (const auto ppn_bcid : ppn_bcids) - if (_periodic_bcs->find(ppn_bcid) != _periodic_bcs->end()) + 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(ppn, nullcm); + jump: + if (on_periodic_boundary) + coupled_elements.emplace(appn, nullcm); - periodic_elems_examined.insert(ppn); + periodic_elems_examined.insert(appn); } } } From b7c6704b913ec965602c49b12fbf688b16bbcab8 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 2 Oct 2020 12:16:07 -0700 Subject: [PATCH 4/5] Don't duplicate default ghosting when copy constructing new mesh Since we now clone all the ghosting functors coming from the other mesh, the old code added the new object's default ghosting functor as well as a clone of the other object's default ghosting functor. With the new code path we only add one default ghosting-like object --- src/mesh/mesh_base.C | 58 ++++++++++++++++++++++---------------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 3ce888d6edf..8293cd54ae5 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -104,39 +104,39 @@ MeshBase::MeshBase (const MeshBase & other_mesh) : _default_ghosting(libmesh_make_unique(*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 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 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(); } From 14490ef36b1da13baa6462c492184be053cbd8dd Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Sat, 3 Oct 2020 18:12:49 -0700 Subject: [PATCH 5/5] Don't use _mesh->processor_id(); use p like I'm meant to --- src/base/ghost_point_neighbors.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/base/ghost_point_neighbors.C b/src/base/ghost_point_neighbors.C index 5654c616283..58b944d75d4 100644 --- a/src/base/ghost_point_neighbors.C +++ b/src/base/ghost_point_neighbors.C @@ -146,7 +146,7 @@ void GhostPointNeighbors::operator() { // 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() == _mesh->processor_id() || + if (appn == remote_elem || appn->processor_id() == p || periodic_elems_examined.count(appn)) continue;