From 0f276f9a013f5c4ec1be15d35add0e7014034bba Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 15 Jan 2026 13:27:23 -0600 Subject: [PATCH 1/7] Add Elem::operator!= --- include/geom/elem.h | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/include/geom/elem.h b/include/geom/elem.h index 833645d96d..4f4665c80d 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -290,13 +290,19 @@ class Elem : public ReferenceCountedObject, virtual dof_id_type key () const; /** - * \returns \p true if two elements are equivalent, false otherwise. - * This is true if the elements are connected to identical global - * nodes, regardless of how those nodes might be numbered local - * to the elements. + * \returns \p true if two elements are equivalent, \p false + * otherwise. This is true if the elements are connected to + * identical global nodes, regardless of how those nodes might be + * numbered local to the elements. */ bool operator == (const Elem & rhs) const; + /** + * \returns \p false if two elements are equivalent, \p true + * otherwise. + */ + bool operator != (const Elem & rhs) const; + /** * \returns \p true if two elements have equal topologies, false * otherwise. @@ -2579,6 +2585,14 @@ subdomain_id_type & Elem::subdomain_id () +inline +bool Elem::operator != (const Elem & rhs) const +{ + return !(*this == rhs); +} + + + inline const Elem * Elem::neighbor_ptr (unsigned int i) const { From f2034d34cd948758c29ff497fda06b3feaa25e1b Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:04:27 -0600 Subject: [PATCH 2/7] Add Elem::total_family_tree_by_side() --- include/geom/elem.h | 15 +++++++++++++++ include/geom/elem_internal.h | 29 +++++++++++++++++++++++++++++ src/geom/elem.C | 18 ++++++++++++++++++ 3 files changed, 62 insertions(+) diff --git a/include/geom/elem.h b/include/geom/elem.h index 4f4665c80d..0eb92f8202 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -1591,6 +1591,21 @@ class Elem : public ReferenceCountedObject, unsigned int side, bool reset = true); + /** + * Same as the \p total_family_tree() member, but only adds elements + * which are next to \p side. + */ + void total_family_tree_by_side (std::vector & family, + unsigned int side, + bool reset = true) const; + + /** + * Non-const version of function above; fills a vector of non-const pointers. + */ + void total_family_tree_by_side (std::vector & family, + unsigned int side, + bool reset = true); + /** * Same as the \p active_family_tree() member, but only adds elements * which are next to \p side. diff --git a/include/geom/elem_internal.h b/include/geom/elem_internal.h index 16144ac5a4..42138ce687 100644 --- a/include/geom/elem_internal.h +++ b/include/geom/elem_internal.h @@ -144,6 +144,35 @@ family_tree_by_side (T elem, +template +void +total_family_tree_by_side (T elem, + std::vector & family, + unsigned int s, + bool reset) +{ + // Clear the vector if the flag reset tells us to. + if (reset) + family.clear(); + + libmesh_assert_less (s, elem->n_sides()); + + // Add this element to the family tree. + family.push_back(elem); + + // Recurse into the elements children, if it has them. + // Do not clear the vector any more. + if (elem->has_children()) + { + const unsigned int nc = elem->n_children(); + for (unsigned int c = 0; c != nc; c++) + if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, s)) + ElemInternal::total_family_tree_by_side (elem->child_ptr(c), family, s, false); + } +} + + + template void active_family_tree_by_side (T elem, diff --git a/src/geom/elem.C b/src/geom/elem.C index 03ce4cff6f..755dfffd13 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -2134,6 +2134,24 @@ void Elem:: family_tree_by_side (std::vector & family, +void Elem::total_family_tree_by_side (std::vector & family, + unsigned int side, + bool reset) const +{ + ElemInternal::total_family_tree_by_side(this, family, side, reset); +} + + + +void Elem::total_family_tree_by_side (std::vector & family, + unsigned int side, + bool reset) +{ + ElemInternal::total_family_tree_by_side(this, family, side, reset); +} + + + void Elem::active_family_tree_by_side (std::vector & family, unsigned int side, bool reset) const From 9cba5172a5295e23c9e58b8d8c67b6ccd701051c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:06:29 -0600 Subject: [PATCH 3/7] Elem::make_links_to_me_local fixes --- src/geom/elem.C | 62 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 46 insertions(+), 16 deletions(-) diff --git a/src/geom/elem.C b/src/geom/elem.C index 755dfffd13..1d84e6c005 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -1430,8 +1430,11 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn) libmesh_assert(neigh); libmesh_assert(!neigh->is_remote()); + const unsigned int this_level = this->level(); + const unsigned int neigh_level = neigh->level(); + // We never have neighbors more refined than us - libmesh_assert_less_equal (neigh->level(), this->level()); + libmesh_assert_less_equal (neigh_level, this_level); // We never have subactive neighbors of non subactive elements libmesh_assert(!neigh->subactive() || this->subactive()); @@ -1439,16 +1442,10 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn) // If we have a neighbor less refined than us then it must not // have any more refined descendants we could have pointed to // instead. - libmesh_assert((neigh->level() == this->level()) || + libmesh_assert((neigh_level == this_level) || (neigh->active() && !this->subactive()) || (!neigh->has_children() && this->subactive())); - // If neigh is at our level, then its family might have - // remote_elem neighbor links which need to point to us - // instead, but if not, then we're done. - if (neigh->level() != this->level()) - return; - // What side of neigh are we on? nn. // // We can't use the usual Elem method because we're in the middle of @@ -1459,20 +1456,27 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn) // not-technically-neighbors until they're // not-even-geometrically-neighbors. - // Find any elements that ought to point to elem + // Find any elements that might need to point to elem std::vector neigh_family; -#ifdef LIBMESH_ENABLE_AMR - if (this->active()) - neigh->family_tree_by_side(neigh_family, nn); - else -#endif - neigh_family.push_back(neigh); + neigh->total_family_tree_by_side(neigh_family, nn); + + // Pull objects out of the loop to reduce heap operations + std::unique_ptr my_side, neigh_side; // And point them to elem for (auto & neigh_family_member : neigh_family) { // Only subactive elements point to other subactive elements - if (this->subactive() && !neigh_family_member->subactive()) + const bool member_subactive = neigh_family_member->subactive(); + if (this->subactive() && !member_subactive) + continue; + + // neigh (and possibly some of its family) might be at a lower + // level than us, with a neighbor at that lower level. We + // would exit early if neigh was at a lower level, except we do + // want to handle neighbor links from subactive descendents. + const unsigned int member_level = neigh_family_member->level(); + if (member_level < this_level) continue; // Ideally, the neighbor link ought to either be correct @@ -1495,6 +1499,32 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn) (neigh_family_member->neighbor_ptr(nn) == remote_elem)); #endif + // Some of neigh's family might be at a higher (finer) level + // than us, and might need to point back to one of our children + // rather than to us. + if (member_level > this_level) + { + if (this->ancestor()) + continue; + + if (member_subactive && this->has_children()) + continue; + } + + // If neigh is at a coarser level than us, some of neigh's + // subactive family just won't line up with us! + if (neigh_level < this_level && + member_level > neigh_level) + { + libmesh_assert(member_subactive); + + this->side_ptr(my_side, n); + neigh_family_member->side_ptr(neigh_side, nn); + + if (*my_side != *neigh_side) + continue; + } + neigh_family_member->set_neighbor(nn, this); } } From 160bec3129e546cbae88f017bdf0f87dafaacb18 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:06:52 -0600 Subject: [PATCH 4/7] Set remote neighbors when refining elements --- src/geom/elem_refinement.C | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/geom/elem_refinement.C b/src/geom/elem_refinement.C index c81cdcb58e..971b39b747 100644 --- a/src/geom/elem_refinement.C +++ b/src/geom/elem_refinement.C @@ -140,6 +140,23 @@ void Elem::refine (MeshRefinement & mesh_refinement) } } + // Get any new remote neighbor links right; even find_neighbors + // relies on those + for (unsigned int s : this->side_index_range()) + { + if (this->neighbor_ptr(s) != remote_elem) + continue; + + for (unsigned int c = 0; c != nc; c++) + { + Elem * current_child = this->child_ptr(c); + if (current_child != remote_elem && + this->is_child_on_side(c, s)) + current_child->set_neighbor + (s, const_cast(remote_elem)); + } + } + // Un-set my refinement flag now this->set_refinement_flag(Elem::INACTIVE); From 2908212b86d96ab622be77697d9f5ca3a47eb188 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:07:52 -0600 Subject: [PATCH 5/7] Handle subactive elements in Packing unpack --- src/parallel/parallel_elem.C | 134 +++++++++++++++++------------------ 1 file changed, 66 insertions(+), 68 deletions(-) diff --git a/src/parallel/parallel_elem.C b/src/parallel/parallel_elem.C index dbad025296..519eda2db3 100644 --- a/src/parallel/parallel_elem.C +++ b/src/parallel/parallel_elem.C @@ -610,82 +610,79 @@ Packing::unpack (std::vector::const_iterator in, // to update a remote_elem link, and we can check for possible // inconsistencies along the way. // - // For subactive elements, we don't bother keeping neighbor - // links in good shape, so there's nothing we need to set or can - // safely assert here. - if (!elem->subactive()) - for (auto n : elem->side_index_range()) - { - const dof_id_type neighbor_id = - cast_int(*in++); + // Even for subactive elements, we'll try to keep neighbor links + // in good shape now, if only so any future find_neighbors() is + // idempotent. + for (auto n : elem->side_index_range()) + { + const dof_id_type neighbor_id = + cast_int(*in++); - const dof_id_type neighbor_side = - cast_int(*in++); + const dof_id_type neighbor_side = + cast_int(*in++); - // If the sending processor sees a domain boundary here, - // we'd better agree ... unless all we see is a - // remote_elem? In that case maybe we just couldn't keep - // up with a user's delete_elem. Let's trust them. - if (neighbor_id == DofObject::invalid_id) - { - const Elem * my_neigh = elem->neighbor_ptr(n); - if (my_neigh == remote_elem) - elem->set_neighbor(n, nullptr); - else - libmesh_assert (!my_neigh); - continue; - } + // If the sending processor sees a domain boundary here, + // we'd better agree ... unless all we see is a remote_elem? + // In that case maybe we just couldn't keep up with a user's + // delete_elem. Let's trust them. + if (neighbor_id == DofObject::invalid_id) + { + const Elem * my_neigh = elem->neighbor_ptr(n); + if (my_neigh == remote_elem) + elem->set_neighbor(n, nullptr); + else + libmesh_assert (!my_neigh); + continue; + } - // If the sending processor has a remote_elem neighbor here, - // then all we know is that we'd better *not* have a domain - // boundary ... except that maybe it's the *sending* - // processor who missed a delete_elem we saw. - if (neighbor_id == remote_elem->id()) - { - // At this level of the code we can't even assert in - // cases where the neighbor should know what they're - // talking about, so skip it. + // If the sending processor has a remote_elem neighbor here, + // then all we know is that we'd better *not* have a domain + // boundary ... except that maybe it's the *sending* + // processor who missed a delete_elem we saw. + if (neighbor_id == remote_elem->id()) + { + // At this level of the code we can't even assert in + // cases where the neighbor should know what they're + // talking about, so skip it. - // libmesh_assert(elem->neighbor_ptr(n)); - continue; - } + // libmesh_assert(elem->neighbor_ptr(n)); + continue; + } - Elem * neigh = mesh->query_elem_ptr(neighbor_id); + Elem * neigh = mesh->query_elem_ptr(neighbor_id); - // The sending processor sees a neighbor here, so if we - // don't have that neighboring element, then we'd better - // have a remote_elem signifying that fact. - if (!neigh) - { - libmesh_assert_equal_to (elem->neighbor_ptr(n), remote_elem); - continue; - } + // The sending processor sees a neighbor here, so if we + // don't have that neighboring element, then we'd better + // have a remote_elem signifying that fact. + if (!neigh) + { + libmesh_assert_equal_to (elem->neighbor_ptr(n), remote_elem); + continue; + } - // The sending processor has a neighbor here, and we have - // that element, but that does *NOT* mean we're already - // linking to it. Perhaps we initially received both elem - // and neigh from processors on which their mutual link was - // remote? - libmesh_assert(elem->neighbor_ptr(n) == neigh || - elem->neighbor_ptr(n) == remote_elem); - - // If the link was originally remote, we should update it, - // and make sure the appropriate parts of its family link - // back to us. - if (elem->neighbor_ptr(n) == remote_elem) - { - elem->set_neighbor(n, neigh); + // The sending processor has a neighbor here, and we have + // that element, but that does *NOT* mean we're already + // linking to it. Perhaps we initially received both elem + // and neigh from processors on which their mutual link was + // remote? + libmesh_assert(elem->neighbor_ptr(n) == neigh || + elem->neighbor_ptr(n) == remote_elem); + + // If the link was originally remote, we should update it, + // and make sure the appropriate parts of its family link + // back to us. + if (elem->neighbor_ptr(n) == remote_elem) + { + elem->set_neighbor(n, neigh); + if (neighbor_side != libMesh::invalid_uint) elem->make_links_to_me_local(n, neighbor_side); - } - else - libmesh_assert(neigh->level() < elem->level() || - neigh->neighbor_ptr(neighbor_side) == elem); - } - else - // We skip these to go to the boundary information if the element is - // actually subactive - in += 2*elem->n_sides(); + } + else + libmesh_assert(elem->subactive() || + neigh->level() < elem->level() || + neigh->neighbor_ptr(neighbor_side) == elem); + } // Our p level and refinement flags should be "close to" correct // if we're not an active element - we might have a p level @@ -841,7 +838,8 @@ Packing::unpack (std::vector::const_iterator in, // to us. elem->set_neighbor(n, neigh); - elem->make_links_to_me_local(n, neighbor_side); + if (neighbor_side != libMesh::invalid_uint) + elem->make_links_to_me_local(n, neighbor_side); } elem->unpack_indexing(in); From 61d43512f4ce730e0286c4348d4928e21296daf7 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 28 Jan 2026 09:46:30 -0600 Subject: [PATCH 6/7] Weaken overzealous make_links_to_me_local assert --- src/geom/elem.C | 1 + 1 file changed, 1 insertion(+) diff --git a/src/geom/elem.C b/src/geom/elem.C index 1d84e6c005..5a8ed25f4d 100644 --- a/src/geom/elem.C +++ b/src/geom/elem.C @@ -1489,6 +1489,7 @@ void Elem::make_links_to_me_local(unsigned int n, unsigned int nn) #ifdef LIBMESH_ENABLE_AMR libmesh_assert((neigh_family_member->neighbor_ptr(nn) && (neigh_family_member->neighbor_ptr(nn)->active() || + this->is_ancestor_of(neigh_family_member->neighbor_ptr(nn)) || neigh_family_member->neighbor_ptr(nn)->is_ancestor_of(this))) || (neigh_family_member->neighbor_ptr(nn) == remote_elem) || ((this->refinement_flag() == JUST_REFINED) && From 6b8b8701837ecae1fd4757e0ffa8137b97698512 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 28 Jan 2026 10:22:37 -0600 Subject: [PATCH 7/7] Make more remote links local when unpacking This fixed a corner case for me, IIRC with unpacking of element sets where we had a non-remote neighbor link whose backlinks were fine but whose descendants' backlinks weren't. --- src/parallel/parallel_elem.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parallel/parallel_elem.C b/src/parallel/parallel_elem.C index 519eda2db3..5043bf2116 100644 --- a/src/parallel/parallel_elem.C +++ b/src/parallel/parallel_elem.C @@ -674,14 +674,14 @@ Packing::unpack (std::vector::const_iterator in, if (elem->neighbor_ptr(n) == remote_elem) { elem->set_neighbor(n, neigh); - - if (neighbor_side != libMesh::invalid_uint) - elem->make_links_to_me_local(n, neighbor_side); } else libmesh_assert(elem->subactive() || neigh->level() < elem->level() || neigh->neighbor_ptr(neighbor_side) == elem); + + if (neighbor_side != libMesh::invalid_uint) + elem->make_links_to_me_local(n, neighbor_side); } // Our p level and refinement flags should be "close to" correct