From df071fa51236ee9de64ce3c5ccbee4334a9e5df5 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 6 Mar 2026 11:45:15 -0600 Subject: [PATCH 1/7] Add a _mesh_local_subdomains cache --- include/mesh/mesh_base.h | 15 +++++++++++++++ src/mesh/distributed_mesh.C | 2 +- src/mesh/mesh_base.C | 13 +++++++++++++ 3 files changed, 29 insertions(+), 1 deletion(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index cf454353b2..93d865c601 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -2057,6 +2057,15 @@ class MeshBase : public ParallelObject const std::set & get_mesh_subdomains() const { libmesh_assert(this->is_prepared()); return _mesh_subdomains; } + + /** + * \return The cached mesh subdomains. As long as the mesh is prepared, this + * should contain all the subdomain ids across processors. Relies on the mesh + * being prepared + */ + const std::set & get_mesh_local_subdomains() const + { libmesh_assert(this->is_prepared()); return _mesh_local_subdomains; } + #ifdef LIBMESH_ENABLE_PERIODIC /** * Register a pair of boundaries as disjoint neighbor boundary pairs. @@ -2332,6 +2341,12 @@ class MeshBase : public ParallelObject */ std::set _mesh_subdomains; + /** + * We also cache the subdomain ids of the elements owned by this + * processor. + */ + std::set _mesh_local_subdomains; + /** * Map from "element set code" to list of set ids to which that element * belongs (and vice-versa). Remarks: diff --git a/src/mesh/distributed_mesh.C b/src/mesh/distributed_mesh.C index 6087bcada6..24cbaa6475 100644 --- a/src/mesh/distributed_mesh.C +++ b/src/mesh/distributed_mesh.C @@ -1078,7 +1078,7 @@ void DistributedMesh::redistribute () void DistributedMesh::update_post_partitioning () { - // this->recalculate_n_partitions(); + this->UnstructuredMesh::update_post_partitioning(); // Let's do the base class cache clearing first, just in case our // later computations are ever changed to make use of a local diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 72641a3765..54f729b2eb 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -128,6 +128,7 @@ MeshBase::MeshBase (const MeshBase & other_mesh) : _elem_default_orders(other_mesh._elem_default_orders), _supported_nodal_order(other_mesh._supported_nodal_order), _mesh_subdomains(other_mesh._mesh_subdomains), + _mesh_local_subdomains(other_mesh._mesh_local_subdomains), _elemset_codes_inverse_map(other_mesh._elemset_codes_inverse_map), _all_elemset_ids(other_mesh._all_elemset_ids), _spatial_dimension(other_mesh._spatial_dimension), @@ -220,6 +221,7 @@ MeshBase& MeshBase::operator= (MeshBase && other_mesh) _elem_default_orders = std::move(other_mesh.elem_default_orders()); _supported_nodal_order = other_mesh.supported_nodal_order(); _mesh_subdomains = other_mesh._mesh_subdomains; + _mesh_local_subdomains = other_mesh._mesh_local_subdomains; _elemset_codes = std::move(other_mesh._elemset_codes); _elemset_codes_inverse_map = std::move(other_mesh._elemset_codes_inverse_map); _all_elemset_ids = std::move(other_mesh._all_elemset_ids); @@ -338,6 +340,8 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const return false; if (_mesh_subdomains != other_mesh._mesh_subdomains) return false; + if (_mesh_local_subdomains != other_mesh._mesh_local_subdomains) + return false; if (_all_elemset_ids != other_mesh._all_elemset_ids) return false; if (_elem_integer_names != other_mesh._elem_integer_names) @@ -1183,6 +1187,10 @@ void MeshBase::update_post_partitioning() // over local elements is obsolete if our partitioner changed the // definition of "local". _const_active_local_element_stored_range.reset(nullptr); + _mesh_local_subdomains.clear(); + + for (const Elem * elem : this->active_local_element_ptr_range()) + _mesh_local_subdomains.insert(elem->subdomain_id()); } @@ -1992,6 +2000,7 @@ void MeshBase::cache_elem_data() _elem_dims.clear(); _elem_default_orders.clear(); _mesh_subdomains.clear(); + _mesh_local_subdomains.clear(); _supported_nodal_order = MAXIMUM; for (const auto & elem : this->active_element_ptr_range()) @@ -1999,6 +2008,8 @@ void MeshBase::cache_elem_data() _elem_dims.insert(cast_int(elem->dim())); _elem_default_orders.insert(elem->default_order()); _mesh_subdomains.insert(elem->subdomain_id()); + if (elem->processor_id() == this->processor_id()) + _mesh_local_subdomains.insert(elem->subdomain_id()); _supported_nodal_order = static_cast (std::min(static_cast(_supported_nodal_order), @@ -2435,6 +2446,7 @@ MeshBase::copy_cached_data(const MeshBase & other_mesh) this->_elem_default_orders = other_mesh._elem_default_orders; this->_supported_nodal_order = other_mesh._supported_nodal_order; this->_mesh_subdomains = other_mesh._mesh_subdomains; + this->_mesh_local_subdomains = other_mesh._mesh_local_subdomains; } @@ -2689,6 +2701,7 @@ MeshBase::copy_constraint_rows(const SparseMatrix & constraint_operator, (std::min(static_cast(this->_supported_nodal_order), static_cast(added_elem->supported_nodal_order()))); this->_mesh_subdomains.insert(new_sbd_id); + this->_mesh_local_subdomains.insert(new_sbd_id); node_to_elem_ptrs.emplace(n, std::make_pair(added_elem->id(), 0)); existing_unconstrained_columns.emplace(j,n->id()); From 125ed26967fe1c10aab5603ddf45b3e95f74658c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 6 Mar 2026 11:45:56 -0600 Subject: [PATCH 2/7] Use caches when returning or counting subdomains Especially n_subdomains() should be an O(1) call, not O(N). --- include/mesh/mesh_base.h | 4 +++ src/mesh/mesh_base.C | 74 ++++++++++++++++------------------------ 2 files changed, 33 insertions(+), 45 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index 93d865c601..858c500804 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -1503,6 +1503,8 @@ class MeshBase : public ParallelObject * materials in a solid mechanics application, or regions where different * physical processes are important. The subdomain mapping is independent * from the parallel decomposition. + * + * This relies on the mesh cached element data being prepared. */ subdomain_id_type n_subdomains () const; @@ -1512,6 +1514,8 @@ class MeshBase : public ParallelObject * materials in a solid mechanics application, or regions where different * physical processes are important. The subdomain mapping is independent * from the parallel decomposition. + * + * This relies on the mesh cached element data being prepared. */ subdomain_id_type n_local_subdomains () const; diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 54f729b2eb..820906ba10 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -54,6 +54,26 @@ #include "libmesh/periodic_boundaries.h" #include "libmesh/periodic_boundary.h" +namespace { + +struct SBDInserter { + std::set my_ids; + + SBDInserter () {} + SBDInserter (SBDInserter &, Threads::split) {} + + void operator()(const ConstElemRange & range) { + for (const Elem * elem : range) + my_ids.insert(elem->subdomain_id()); + } + + void join(SBDInserter & other) { + my_ids.merge(other.my_ids); + } +}; + +} + namespace libMesh { @@ -1129,42 +1149,12 @@ void MeshBase::remove_ghosting_functor(GhostingFunctor & ghosting_functor) void MeshBase::subdomain_ids (std::set & ids, const bool global /* = true */) const { - // This requires an inspection on every processor - if (global) - parallel_object_only(); - - struct SBDInserter { - std::set my_ids; - - SBDInserter () {} - SBDInserter (SBDInserter &, Threads::split) {} - - void operator()(const ConstElemRange & range) { - for (const Elem * elem : range) - my_ids.insert(elem->subdomain_id()); - } - - void join(SBDInserter & other) { - my_ids.merge(other.my_ids); - } - }; - - SBDInserter inserter; - Threads::parallel_reduce(this->active_local_element_stored_range(), inserter); - - ids.swap(inserter.my_ids); + libmesh_assert(this->preparation().has_cached_elem_data); if (global) - { - // Only include the unpartitioned elements if the user requests the global IDs. - // In the case of the local subdomain IDs, it doesn't make sense to include the - // unpartitioned elements because said elements do not have a sense of locality. - for (const auto & elem : this->active_unpartitioned_element_ptr_range()) - ids.insert(elem->subdomain_id()); - - // Some subdomains may only live on other processors - this->comm().set_union(ids); - } + ids = this->get_mesh_subdomains(); + else + ids = this->get_mesh_local_subdomains(); } @@ -1187,24 +1177,18 @@ void MeshBase::update_post_partitioning() // over local elements is obsolete if our partitioner changed the // definition of "local". _const_active_local_element_stored_range.reset(nullptr); - _mesh_local_subdomains.clear(); - for (const Elem * elem : this->active_local_element_ptr_range()) - _mesh_local_subdomains.insert(elem->subdomain_id()); + SBDInserter inserter; + Threads::parallel_reduce(this->active_local_element_stored_range(), inserter); + + _mesh_local_subdomains.swap(inserter.my_ids); } subdomain_id_type MeshBase::n_subdomains() const { - // This requires an inspection on every processor - parallel_object_only(); - - std::set ids; - - this->subdomain_ids (ids); - - return cast_int(ids.size()); + return cast_int(this->get_mesh_subdomains().size()); } From 865bb1250db6623cf7094a8122dc6f451916cb06 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 16 Mar 2026 11:27:47 -0500 Subject: [PATCH 3/7] update_post_partitioning even on trivial partition This fixes the case where in serial we didn't know how many local subdomains we had, probably among other issues. --- src/partitioning/partitioner.C | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/partitioning/partitioner.C b/src/partitioning/partitioner.C index b05db2d87f..c3afee769e 100644 --- a/src/partitioning/partitioner.C +++ b/src/partitioning/partitioner.C @@ -215,6 +215,11 @@ void Partitioner::partition (MeshBase & mesh, if (n_parts == 1) { this->single_partition (mesh); + + // Give derived Mesh classes a chance to update any cached data + // to reflect the new partitioning + mesh.update_post_partitioning(); + return; } From 77453997944c21d7976f4ef45090bf19dc874da9 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 17 Mar 2026 13:14:58 -0500 Subject: [PATCH 4/7] Avoid false positives in _mesh_local_subdomains If we thought of making a constraining NodeElem local and then changed our minds, that means it's not local! --- src/mesh/mesh_base.C | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 820906ba10..c95a1e8ea5 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -2685,7 +2685,6 @@ MeshBase::copy_constraint_rows(const SparseMatrix & constraint_operator, (std::min(static_cast(this->_supported_nodal_order), static_cast(added_elem->supported_nodal_order()))); this->_mesh_subdomains.insert(new_sbd_id); - this->_mesh_local_subdomains.insert(new_sbd_id); node_to_elem_ptrs.emplace(n, std::make_pair(added_elem->id(), 0)); existing_unconstrained_columns.emplace(j,n->id()); @@ -2693,13 +2692,18 @@ MeshBase::copy_constraint_rows(const SparseMatrix & constraint_operator, // DistributedMesh doesn't get confused and think you're not // adding them on all processors at once. int n_pids = 0; + processor_id_type best_pid = DofObject::invalid_processor_id; for (auto [pid, count] : pids) if (count >= n_pids) { n_pids = count; - added_elem->processor_id() = pid; - n->processor_id() = pid; + best_pid = pid; } + libmesh_assert_not_equal_to(best_pid, DofObject::invalid_processor_id); + added_elem->processor_id() = best_pid; + n->processor_id() = best_pid; + if (this->processor_id() == best_pid) + this->_mesh_local_subdomains.insert(new_sbd_id); } // Calculate constraint rows in an indexed form that's easy for us From 189aee62a699ce2bbedd5b9406d0ce0c134b6714 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 17 Mar 2026 13:15:51 -0500 Subject: [PATCH 5/7] Fix mesh preparation in unit test --- tests/mesh/mesh_elem_test.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/mesh/mesh_elem_test.h b/tests/mesh/mesh_elem_test.h index 9fd75b94c1..f459a2843b 100644 --- a/tests/mesh/mesh_elem_test.h +++ b/tests/mesh/mesh_elem_test.h @@ -108,6 +108,11 @@ class MeshPerElemTest : public PerElemTest #endif } + // Setting all those processor ids to 0 changes our sets of local + // subdomains too. + other_mesh.cache_elem_data(); + this->_mesh->cache_elem_data(); + #ifdef LIBMESH_ENABLE_UNIQUE_ID other_mesh.set_next_unique_id(this->_mesh->parallel_max_unique_id()); this->_mesh->set_next_unique_id(this->_mesh->parallel_max_unique_id()); From ef2fd2f3e623b5e806287d36dbb11ff656a97420 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 24 Mar 2026 11:06:38 -0500 Subject: [PATCH 6/7] Weaken preparation assert when getting subdomains We need the caches here to be ready but we don't need full preparation. --- include/mesh/mesh_base.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index 858c500804..fc83477c52 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -2059,7 +2059,7 @@ class MeshBase : public ParallelObject * being prepared */ const std::set & get_mesh_subdomains() const - { libmesh_assert(this->is_prepared()); return _mesh_subdomains; } + { libmesh_assert(this->preparation().has_cached_elem_data); return _mesh_subdomains; } /** @@ -2068,7 +2068,7 @@ class MeshBase : public ParallelObject * being prepared */ const std::set & get_mesh_local_subdomains() const - { libmesh_assert(this->is_prepared()); return _mesh_local_subdomains; } + { libmesh_assert(this->preparation().has_cached_elem_data); return _mesh_local_subdomains; } #ifdef LIBMESH_ENABLE_PERIODIC /** From 77fab2141123864eda767024751822b7fb678d97 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 24 Mar 2026 13:04:35 -0500 Subject: [PATCH 7/7] cache_elem_data() in Tetgen triangulate() --- src/mesh/mesh_tetgen_interface.C | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/mesh/mesh_tetgen_interface.C b/src/mesh/mesh_tetgen_interface.C index 081ee58499..8923a7a19d 100644 --- a/src/mesh/mesh_tetgen_interface.C +++ b/src/mesh/mesh_tetgen_interface.C @@ -125,6 +125,12 @@ void TetGenMeshInterface::triangulate_pointset () // We don't do this by default. if (this->_smooth_after_generating) LaplaceMeshSmoother(this->_mesh, 2).smooth(); + + // We've added a bunch of elements. A full prepare_for_use() would + // be expensive here and user code hasn't been expecting it, but we + // do have code downstream expecting element caches to be up to + // date. + this->_mesh.cache_elem_data(); }