From 7ca5d6f08d9871b9762073a32ff95277a2a32f52 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Jan 2024 12:59:24 -0600 Subject: [PATCH 01/48] Make is_prepared more fine-grained --- include/mesh/mesh_base.h | 200 +++++++++++++++++++++++++++++++---- src/mesh/boundary_info.C | 2 + src/mesh/distributed_mesh.C | 8 +- src/mesh/mesh_base.C | 90 +++++++++++----- src/mesh/replicated_mesh.C | 6 +- src/mesh/unstructured_mesh.C | 4 +- 6 files changed, 260 insertions(+), 50 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index b1c40a44ee8..c30a559a494 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -165,7 +165,7 @@ class MeshBase : public ParallelObject virtual ~MeshBase (); /** - * A partitioner to use at each prepare_for_use() + * A partitioner to use at each partitioning */ virtual std::unique_ptr & partitioner() { return _partitioner; } @@ -194,7 +194,7 @@ class MeshBase : public ParallelObject * No Node is removed from the mesh, however even NodeElem elements * are deleted, so the remaining Nodes will be considered "unused" * and cleared unless they are reconnected to new elements before - * the next \p prepare_for_use() + * the next preparation step. * * This does not affect BoundaryInfo data; any boundary information * associated elements should already be cleared. @@ -202,17 +202,105 @@ class MeshBase : public ParallelObject virtual void clear_elems () = 0; /** - * \returns \p true if the mesh has been prepared via a call - * to \p prepare_for_use, \p false otherwise. + * \returns \p true if the mesh has undergone all the preparation + * done in a call to \p prepare_for_use, \p false otherwise. */ bool is_prepared () const - { return _is_prepared; } + { return _preparation; } /** - * Tells this we have done some operation where we should no longer consider ourself prepared + * \returns the \p Preparation structure with details about in what + * ways \p this mesh is currently prepared or unprepared. This + * structure may change in the future when cache designs change. + */ + struct Preparation; + + Preparation preparation () const + { return _preparation; } + + /** + * Tells this we have done some operation where we should no longer + * consider ourself prepared. This is a very coarse setting; it may + * be preferable to mark finer-grained settings instead. */ void set_isnt_prepared() - { _is_prepared = false; } + { _preparation = false; } + + /** + * Tells this we have done some operation creating unpartitioned + * elements. + */ + void set_isnt_partitioned() + { _preparation.is_partitioned = false; } + + /** + * Tells this we have done some operation (e.g. adding elements to a + * distributed mesh on one processor only) which can lose + * synchronization of id counts. + */ + void set_hasnt_synched_id_counts() + { _preparation.has_synched_id_counts = false; } + + /** + * Tells this we have done some operation (e.g. adding elements + * without setting their neighbor pointers) which requires neighbor + * pointers to be found later + */ + void set_hasnt_neighbor_ptrs() + { _preparation.has_neighbor_ptrs = false; } + + /** + * Tells this we have done some operation (e.g. adding elements with + * a new dimension or subdomain value) which may invalidate cached + * summaries of element data + */ + void set_hasnt_cached_elem_data() + { _preparation.has_cached_elem_data = false; } + + /** + * Tells this we have done some operation (e.g. refining elements + * with interior parents) which requires interior parent pointers to + * be found later + */ + void set_hasnt_interior_parent_ptrs() + { _preparation.has_interior_parent_ptrs = false; } + + /** + * Tells this we have done some operation (e.g. repartitioning) + * which may have left elements as ghosted which on a distributed + * mesh should be remote. + * + * User code should probably never need to use this; we can set it + * in Partitioner. + */ + void set_hasnt_removed_remote_elements() + { _preparation.has_removed_remote_elements = false; } + + /** + * Tells this we have done some operation (e.g. coarsening) + * which may have left orphaned nodes in need of removal. + * + * Most user code should probably never need to use this; we can set + * it in MeshRefinement. + */ + void set_hasnt_removed_orphaned_nodes() + { _preparation.has_removed_orphaned_nodes = false; } + + /** + * Tells this we have done some operation (e.g. adding or removing + * elements) which may require a reinit() of custom ghosting + * functors. + */ + void hasnt_reinit_ghosting_functors() + { _preparation.has_reinit_ghosting_functors = false; } + + /** + * Tells this we have done some operation (e.g. removing elements or + * adding new boundary entries) which may have invalidated our + * cached boundary id sets. + */ + void set_hasnt_boundary_id_sets() + { _preparation.has_boundary_id_sets = false; } /** * \returns \p true if all elements and nodes of the mesh @@ -260,7 +348,9 @@ class MeshBase : public ParallelObject * except for "ghosts" which touch a local element, and deletes * all nodes which are not part of a local or ghost element */ - virtual void delete_remote_elements () {} + virtual void delete_remote_elements () { + _preparation.has_removed_remote_elements = true; + } /** * Loops over ghosting functors and calls mesh_reinit() @@ -742,7 +832,7 @@ class MeshBase : public ParallelObject * To ensure a specific element id, call e->set_id() before adding it; * only do this in parallel if you are manually keeping ids consistent. * - * Users should call MeshBase::prepare_for_use() after elements are + * Users should call MeshBase::complete_preparation() after elements are * added to and/or deleted from the mesh. */ virtual Elem * add_elem (Elem * e) = 0; @@ -761,7 +851,7 @@ class MeshBase : public ParallelObject * Insert elem \p e to the element array, preserving its id * and replacing/deleting any existing element with the same id. * - * Users should call MeshBase::prepare_for_use() after elements are + * Users should call MeshBase::complete_preparation() after elements are * added to and/or deleted from the mesh. */ virtual Elem * insert_elem (Elem * e) = 0; @@ -780,8 +870,8 @@ class MeshBase : public ParallelObject * Removes element \p e from the mesh. This method must be * implemented in derived classes in such a way that it does not * invalidate element iterators. Users should call - * MeshBase::prepare_for_use() after elements are added to and/or - * deleted from the mesh. + * MeshBase::complete_preparation() after elements are added to + * and/or deleted from the mesh. * * \note Calling this method may produce isolated nodes, i.e. nodes * not connected to any element. @@ -848,7 +938,7 @@ class MeshBase : public ParallelObject /** * Removes any orphaned nodes, nodes not connected to any elements. - * Typically done automatically in prepare_for_use + * Typically done automatically in a preparation step */ void remove_orphaned_nodes (); @@ -1122,12 +1212,19 @@ class MeshBase : public ParallelObject const std::vector * default_values = nullptr); /** - * Prepare a newly ecreated (or read) mesh for use. - * This involves 4 steps: - * 1.) call \p find_neighbors() - * 2.) call \p partition() - * 3.) call \p renumber_nodes_and_elements() - * 4.) call \p cache_elem_data() + * Prepare a newly created (or read) mesh for use. + * This involves several steps: + * 1.) renumbering (if enabled) + * 2.) removing any orphaned nodes + * 3.) updating parallel id counts + * 4.) finding neighbor links + * 5.) caching summarized element data + * 6.) finding interior parent links + * 7.) clearing any old point locator + * 8.) calling reinit() on ghosting functors + * 9.) repartitioning (if enabled) + * 10.) removing any remote elements (if enabled) + * 11.) regenerating summarized boundary id sets * * The argument to skip renumbering is now deprecated - to prevent a * mesh from being renumbered, set allow_renumbering(false). The argument to skip @@ -1144,6 +1241,15 @@ class MeshBase : public ParallelObject #endif // LIBMESH_ENABLE_DEPRECATED void prepare_for_use (); + /* + * Prepare a newly created or modified mesh for use. + * + * Unlike \p prepare_for_use(), \p complete_preparation() performs + * *only* those preparatory steps that have been marked as + * necessary in the MeshBase::Preparation state. + */ + void complete_preparation(); + /** * Call the default partitioner (currently \p metis_partition()). */ @@ -1844,6 +1950,58 @@ class MeshBase : public ParallelObject const boundary_id_type b2); #endif + /** + * Flags indicating in what ways a mesh has been prepared for use. + */ + struct Preparation { + bool is_partitioned = false, + has_synched_id_counts = false, + has_neighbor_ptrs = false, + has_cached_elem_data = false, + has_interior_parent_ptrs = false, + has_removed_remote_elements = false, + has_removed_orphaned_nodes = false, + has_boundary_id_sets = false, + has_reinit_ghosting_functors = false; + + operator bool() const { + return is_partitioned && + has_synched_id_counts && + has_neighbor_ptrs && + has_cached_elem_data && + has_interior_parent_ptrs && + has_removed_remote_elements && + has_removed_orphaned_nodes && + has_reinit_ghosting_functors && + has_boundary_id_sets; + } + + Preparation & operator= (bool set_all) { + is_partitioned = set_all; + has_synched_id_counts = set_all; + has_neighbor_ptrs = set_all; + has_cached_elem_data = set_all; + has_interior_parent_ptrs = set_all; + has_removed_remote_elements = set_all; + has_removed_orphaned_nodes = set_all; + has_reinit_ghosting_functors = set_all; + has_boundary_id_sets = set_all; + + return *this; + } + + bool operator== (const Preparation & other) { + return is_partitioned == other.is_partitioned && + has_synched_id_counts == other.has_synched_id_counts && + has_neighbor_ptrs == other.has_neighbor_ptrs && + has_cached_elem_data == other.has_cached_elem_data && + has_interior_parent_ptrs == other.has_interior_parent_ptrs && + has_removed_remote_elements == other.has_removed_remote_elements && + has_removed_orphaned_nodes == other.has_removed_orphaned_nodes && + has_reinit_ghosting_functors == other.has_reinit_ghosting_functors && + has_boundary_id_sets == other.has_boundary_id_sets; + } + }; protected: @@ -1923,9 +2081,9 @@ class MeshBase : public ParallelObject unsigned char _default_mapping_data; /** - * Flag indicating if the mesh has been prepared for use. + * Flags indicating in what ways \p this mesh has been prepared. */ - bool _is_prepared; + Preparation _preparation; /** * A \p PointLocator class for this mesh. diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C index 29a3f530e4d..cda463c554d 100644 --- a/src/mesh/boundary_info.C +++ b/src/mesh/boundary_info.C @@ -433,6 +433,8 @@ void BoundaryInfo::synchronize_global_id_set() libmesh_assert(_mesh); if (!_mesh->is_serial()) _communicator.set_union(_global_boundary_ids); + + _mesh->_preparation.has_boundary_id_sets = true; } diff --git a/src/mesh/distributed_mesh.C b/src/mesh/distributed_mesh.C index 24e352322d2..5a15a04814b 100644 --- a/src/mesh/distributed_mesh.C +++ b/src/mesh/distributed_mesh.C @@ -205,7 +205,7 @@ DistributedMesh::DistributedMesh (const MeshBase & other_mesh) : this->copy_constraint_rows(other_mesh); - this->_is_prepared = other_mesh.is_prepared(); + this->_preparation = other_mesh.preparation(); auto & this_boundary_info = this->get_boundary_info(); const auto & other_boundary_info = other_mesh.get_boundary_info(); @@ -291,6 +291,8 @@ void DistributedMesh::update_parallel_id_counts() ((_next_unique_id + this->n_processors() - 1) / (this->n_processors() + 1) + 1) * (this->n_processors() + 1) + this->processor_id(); #endif + + this->_preparation.has_synched_id_counts = true; } @@ -1611,6 +1613,8 @@ void DistributedMesh::renumber_nodes_and_elements () } } + this->_preparation.has_removed_orphaned_nodes = true; + if (_skip_renumber_nodes_and_elements) { this->update_parallel_id_counts(); @@ -1776,6 +1780,8 @@ void DistributedMesh::delete_remote_elements() this->libmesh_assert_valid_parallel_ids(); this->libmesh_assert_valid_parallel_flags(); #endif + + this->_preparation.has_removed_remote_elements = true; } diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 535117af6b5..28ee49c5514 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -66,7 +66,7 @@ MeshBase::MeshBase (const Parallel::Communicator & comm_in, _n_parts (1), _default_mapping_type(LAGRANGE_MAP), _default_mapping_data(0), - _is_prepared (false), + _preparation (), _point_locator (), _count_lower_dim_elems_in_point_locator(true), _partitioner (), @@ -99,7 +99,7 @@ MeshBase::MeshBase (const MeshBase & other_mesh) : _n_parts (other_mesh._n_parts), _default_mapping_type(other_mesh._default_mapping_type), _default_mapping_data(other_mesh._default_mapping_data), - _is_prepared (other_mesh._is_prepared), + _preparation (other_mesh._preparation), _point_locator (), _count_lower_dim_elems_in_point_locator(other_mesh._count_lower_dim_elems_in_point_locator), _partitioner (), @@ -187,7 +187,7 @@ MeshBase& MeshBase::operator= (MeshBase && other_mesh) _n_parts = other_mesh.n_partitions(); _default_mapping_type = other_mesh.default_mapping_type(); _default_mapping_data = other_mesh.default_mapping_data(); - _is_prepared = other_mesh.is_prepared(); + _preparation = other_mesh._preparation; _point_locator = std::move(other_mesh._point_locator); _count_lower_dim_elems_in_point_locator = other_mesh.get_count_lower_dim_elems_in_point_locator(); #ifdef LIBMESH_ENABLE_UNIQUE_ID @@ -283,7 +283,7 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const return false; if (_default_mapping_data != other_mesh._default_mapping_data) return false; - if (_is_prepared != other_mesh._is_prepared) + if (_preparation != other_mesh._preparation) return false; if (_count_lower_dim_elems_in_point_locator != other_mesh._count_lower_dim_elems_in_point_locator) @@ -794,6 +794,8 @@ void MeshBase::remove_orphaned_nodes () for (const auto & node : this->node_ptr_range()) if (!connected_nodes.count(node)) this->delete_node(node); + + _preparation.has_removed_orphaned_nodes = true; } @@ -834,9 +836,23 @@ void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements) + void MeshBase::prepare_for_use () { - LOG_SCOPE("prepare_for_use()", "MeshBase"); + // Mark everything as unprepared, except for those things we've been + // told we don't need to prepare, for backwards compatibility + this->clear_point_locator(); + _preparation = false; + _preparation.has_neighbor_ptrs = _skip_find_neighbors; + _preparation.has_removed_remote_elements = !_allow_remote_element_removal; + + this->complete_preparation(); +} + + +void MeshBase::complete_preparation() +{ + LOG_SCOPE("complete_preparation()", "MeshBase"); parallel_object_only(); @@ -871,15 +887,21 @@ void MeshBase::prepare_for_use () // using, but our partitioner might need that consistency and/or // might be confused by orphaned nodes. if (!_skip_renumber_nodes_and_elements) - this->renumber_nodes_and_elements(); + { + if (!_preparation.has_removed_orphaned_nodes || + !_preparation.has_synched_id_counts) + this->renumber_nodes_and_elements(); + } else { - this->remove_orphaned_nodes(); - this->update_parallel_id_counts(); + if (!_preparation.has_removed_orphaned_nodes) + this->remove_orphaned_nodes(); + if (!_preparation.has_synched_id_counts) + this->update_parallel_id_counts(); } // Let all the elements find their neighbors - if (!_skip_find_neighbors) + if (!_skip_find_neighbors && !_preparation.has_neighbor_ptrs) this->find_neighbors(); // The user may have set boundary conditions. We require that the @@ -892,11 +914,13 @@ void MeshBase::prepare_for_use () // Search the mesh for all the dimensions of the elements // and cache them. - this->cache_elem_data(); + if (!_preparation.has_cached_elem_data) + this->cache_elem_data(); // Search the mesh for elements that have a neighboring element // of dim+1 and set that element as the interior parent - this->detect_interior_parents(); + if (!_preparation.has_interior_parent_ptrs) + this->detect_interior_parents(); // Fix up node unique ids in case mesh generation code didn't take // exceptional care to do so. @@ -908,39 +932,41 @@ void MeshBase::prepare_for_use () MeshTools::libmesh_assert_valid_unique_ids(*this); #endif - // Reset our PointLocator. Any old locator is invalidated any time - // the elements in the underlying elements in the mesh have changed, - // so we clear it here. - this->clear_point_locator(); - // Allow our GhostingFunctor objects to reinit if necessary. // Do this before partitioning and redistributing, and before // deleting remote elements. - this->reinit_ghosting_functors(); + if (!_preparation.has_reinit_ghosting_functors) + this->reinit_ghosting_functors(); // Partition the mesh unless *all* partitioning is to be skipped. // If only noncritical partitioning is to be skipped, the // partition() call will still check for orphaned nodes. - if (!skip_partitioning()) + if (!skip_partitioning() && !_preparation.is_partitioned) this->partition(); + else + _preparation.is_partitioned = true; // If we're using DistributedMesh, we'll probably want it // parallelized. - if (this->_allow_remote_element_removal) + if (this->_allow_remote_element_removal && + !_preparation.has_removed_remote_elements) this->delete_remote_elements(); + else + _preparation.has_removed_remote_elements = true; // Much of our boundary info may have been for now-remote parts of the mesh, // in which case we don't want to keep local copies of data meant to be // local. On the other hand we may have deleted, or the user may have added in // a distributed fashion, boundary data that is meant to be global. So we // handle both of those scenarios here - this->get_boundary_info().regenerate_id_sets(); + if (!_preparation.has_boundary_id_sets) + this->get_boundary_info().regenerate_id_sets(); if (!_skip_renumber_nodes_and_elements) this->renumber_nodes_and_elements(); - // The mesh is now prepared for use. - _is_prepared = true; + // The mesh is now prepared for use, and it should know it. + libmesh_assert(_preparation); #ifdef DEBUG MeshTools::libmesh_assert_valid_boundary_ids(*this); @@ -958,6 +984,8 @@ MeshBase::reinit_ghosting_functors() libmesh_assert(gf); gf->mesh_reinit(); } + + _preparation.has_reinit_ghosting_functors = true; } void MeshBase::clear () @@ -965,8 +993,8 @@ void MeshBase::clear () // Reset the number of partitions _n_parts = 1; - // Reset the _is_prepared flag - _is_prepared = false; + // Reset the preparation flags + _preparation = false; // Clear boundary information if (boundary_info) @@ -1683,6 +1711,8 @@ void MeshBase::partition (const unsigned int n_parts) // Make sure any other locally cached data is correct this->update_post_partitioning(); } + + _preparation.is_partitioned = true; } void MeshBase::all_second_order (const bool full_ordered) @@ -1886,6 +1916,8 @@ void MeshBase::cache_elem_data() #endif } } + + _preparation.has_cached_elem_data = true; } @@ -1903,10 +1935,16 @@ void MeshBase::detect_interior_parents() // This requires an inspection on every processor parallel_object_only(); + // This requires up-to-date mesh dimensions in cache + libmesh_assert(_preparation.has_cached_elem_data); + // Check if the mesh contains mixed dimensions. If so, then we may // have interior parents to set. Otherwise return. if (this->elem_dimensions().size() == 1) - return; + { + _preparation.has_interior_parent_ptrs = true; + return; + } // Do we have interior parent pointers going to a different mesh? // If so then we'll still check to make sure that's the only place @@ -2002,6 +2040,8 @@ void MeshBase::detect_interior_parents() ("interior_parent() values in multiple meshes are unsupported."); } } + + _preparation.has_interior_parent_ptrs = true; } diff --git a/src/mesh/replicated_mesh.C b/src/mesh/replicated_mesh.C index 059d3fcbe78..b58d941156b 100644 --- a/src/mesh/replicated_mesh.C +++ b/src/mesh/replicated_mesh.C @@ -115,7 +115,7 @@ ReplicatedMesh::ReplicatedMesh (const MeshBase & other_mesh) : this->copy_constraint_rows(other_mesh); - this->_is_prepared = other_mesh.is_prepared(); + this->_preparation = other_mesh.preparation(); auto & this_boundary_info = this->get_boundary_info(); const auto & other_boundary_info = other_mesh.get_boundary_info(); @@ -643,6 +643,8 @@ void ReplicatedMesh::update_parallel_id_counts() #ifdef LIBMESH_ENABLE_UNIQUE_ID _next_unique_id = this->parallel_max_unique_id(); #endif + + this->_preparation.has_synched_id_counts = true; } @@ -820,6 +822,8 @@ void ReplicatedMesh::renumber_nodes_and_elements () } } + this->_preparation.has_removed_orphaned_nodes = true; + libmesh_assert_equal_to (next_free_elem, _elements.size()); libmesh_assert_equal_to (next_free_node, _nodes.size()); diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C index dfa28eb5047..d3a09d63cec 100644 --- a/src/mesh/unstructured_mesh.C +++ b/src/mesh/unstructured_mesh.C @@ -1256,15 +1256,15 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements, libmesh_assert(current_elem->interior_parent()); } } - #endif // AMR - #ifdef DEBUG MeshTools::libmesh_assert_valid_neighbors(*this, !reset_remote_elements); MeshTools::libmesh_assert_valid_amr_interior_parents(*this); #endif + + this->_preparation.has_neighbor_ptrs = true; } From 1de5774f8c0d5d1c4e4e80fd4fb4b0c771f55284 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 22 Jan 2024 10:44:24 -0600 Subject: [PATCH 02/48] Assert our cache is still good --- src/mesh/mesh_base.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 28ee49c5514..9a4301d3dee 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -585,6 +585,8 @@ void MeshBase::change_elemset_id(elemset_id_type old_id, elemset_id_type new_id) unsigned int MeshBase::spatial_dimension () const { + libmesh_assert(_preparation.has_cached_elem_data); + return cast_int(_spatial_dimension); } From bc61dcd7aa82eb3ac68c7a826ee467e0882b6ca2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 1 Dec 2025 11:51:52 -0600 Subject: [PATCH 03/48] Make PerfLog name accurate --- src/mesh/mesh_tools.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C index 7e98732dcb9..7fdcf9406a9 100644 --- a/src/mesh/mesh_tools.C +++ b/src/mesh/mesh_tools.C @@ -1227,7 +1227,7 @@ void clear_spline_nodes(MeshBase & mesh) bool valid_is_prepared (const MeshBase & mesh) { - LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools"); + LOG_SCOPE("valid_is_prepared()", "MeshTools"); if (!mesh.is_prepared()) return true; From 25e1a6efde0f7aacbd2ba20b90b53e8cf12da0b3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 1 Dec 2025 14:40:57 -0600 Subject: [PATCH 04/48] More detailed preparation comments --- include/mesh/mesh_base.h | 57 ++++++++++++++++++++++++++++++---------- 1 file changed, 43 insertions(+), 14 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index c30a559a494..8aaa3f944cd 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -202,8 +202,9 @@ class MeshBase : public ParallelObject virtual void clear_elems () = 0; /** - * \returns \p true if the mesh has undergone all the preparation - * done in a call to \p prepare_for_use, \p false otherwise. + * \returns \p true if the mesh is marked as having undergone all of + * the preparation done in a call to \p prepare_for_use, \p false + * otherwise. */ bool is_prepared () const { return _preparation; } @@ -220,8 +221,8 @@ class MeshBase : public ParallelObject /** * Tells this we have done some operation where we should no longer - * consider ourself prepared. This is a very coarse setting; it may - * be preferable to mark finer-grained settings instead. + * consider ourself prepared. This is a very coarse setting; it is + * generally more efficient to mark finer-grained settings instead. */ void set_isnt_prepared() { _preparation = false; } @@ -229,22 +230,33 @@ class MeshBase : public ParallelObject /** * Tells this we have done some operation creating unpartitioned * elements. + * + * User code which adds elements to this mesh must either partition + * them too or call this method. */ void set_isnt_partitioned() { _preparation.is_partitioned = false; } /** - * Tells this we have done some operation (e.g. adding elements to a + * Tells this we have done some operation (e.g. adding objects to a * distributed mesh on one processor only) which can lose * synchronization of id counts. + * + * User code which does distributed additions of nodes or elements + * must call either this method or \p update_parallel_id_counts(). */ void set_hasnt_synched_id_counts() { _preparation.has_synched_id_counts = false; } /** * Tells this we have done some operation (e.g. adding elements - * without setting their neighbor pointers) which requires neighbor - * pointers to be found later + * without setting their neighbor pointers, or adding disjoint + * neighbor boundary pairs) which requires neighbor pointers to be + * determined later. + * + * User code which adds new elements to this mesh must call this + * function or manually set neighbor pointer from and to those + * elements. */ void set_hasnt_neighbor_ptrs() { _preparation.has_neighbor_ptrs = false; } @@ -252,7 +264,10 @@ class MeshBase : public ParallelObject /** * Tells this we have done some operation (e.g. adding elements with * a new dimension or subdomain value) which may invalidate cached - * summaries of element data + * summaries of element data. + * + * User code which adds new elements to this mesh must call this + * function. */ void set_hasnt_cached_elem_data() { _preparation.has_cached_elem_data = false; } @@ -260,7 +275,11 @@ class MeshBase : public ParallelObject /** * Tells this we have done some operation (e.g. refining elements * with interior parents) which requires interior parent pointers to - * be found later + * be found later. + * + * Most user code will not need to call this method; any user code + * that manipulates interior parents or their boundary elements may + * be an exception. */ void set_hasnt_interior_parent_ptrs() { _preparation.has_interior_parent_ptrs = false; } @@ -271,7 +290,10 @@ class MeshBase : public ParallelObject * mesh should be remote. * * User code should probably never need to use this; we can set it - * in Partitioner. + * in Partitioner. Any user code which manually repartitions + * elements on distributed meshes may need to call this manually, in + * addition to manually communicating elements with newly-created + * ghosting requirements. */ void set_hasnt_removed_remote_elements() { _preparation.has_removed_remote_elements = false; } @@ -281,7 +303,8 @@ class MeshBase : public ParallelObject * which may have left orphaned nodes in need of removal. * * Most user code should probably never need to use this; we can set - * it in MeshRefinement. + * it in MeshRefinement. User code which deletes elements without + * carefully deleting orphaned nodes should call this manually. */ void set_hasnt_removed_orphaned_nodes() { _preparation.has_removed_orphaned_nodes = false; } @@ -290,14 +313,20 @@ class MeshBase : public ParallelObject * Tells this we have done some operation (e.g. adding or removing * elements) which may require a reinit() of custom ghosting * functors. + * + * User code which adds or removes elements should call this method. + * User code which moves nodes ... should probably call this method, + * in case ghosting functors depending on position exist? */ void hasnt_reinit_ghosting_functors() { _preparation.has_reinit_ghosting_functors = false; } /** - * Tells this we have done some operation (e.g. removing elements or - * adding new boundary entries) which may have invalidated our - * cached boundary id sets. + * Tells this we have done some operation which may have invalidated + * our cached boundary id sets. + * + * User code which removes elements, or which adds or removes + * boundary entries, should call this method. */ void set_hasnt_boundary_id_sets() { _preparation.has_boundary_id_sets = false; } From fca90e5d963a9a1d80aea78a8890bef05606b02f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 2 Dec 2025 22:11:11 -0600 Subject: [PATCH 05/48] Test partially-prepared status validity too --- src/mesh/mesh_tools.C | 65 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 56 insertions(+), 9 deletions(-) diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C index 7fdcf9406a9..0650d3f991f 100644 --- a/src/mesh/mesh_tools.C +++ b/src/mesh/mesh_tools.C @@ -1229,21 +1229,68 @@ bool valid_is_prepared (const MeshBase & mesh) { LOG_SCOPE("valid_is_prepared()", "MeshTools"); - if (!mesh.is_prepared()) + const MeshBase::Preparation prep = mesh.preparation(); + + // If the mesh doesn't think *anything* has been prepared, we have + // nothing to check. + if (prep == MeshBase::Preparation()) return true; + // If the mesh thinks it's partitioned, check. These are counts, + // not caches, so it's a real check. + if (prep.is_partitioned) + if (mesh.n_unpartitioned_elem() || + mesh.n_unpartitioned_nodes()) + return false; + + // If the mesh thinks it's prepared in some way, *re*-preparing in + // that way shouldn't change a clone of it, as long as we disallow + // repartitioning or renumbering or remote element removal. std::unique_ptr mesh_clone = mesh.clone(); - // Try preparing (without allowing repartitioning or renumbering, to - // avoid false assertion failures) - bool old_allow_renumbering = mesh_clone->allow_renumbering(); - mesh_clone->allow_renumbering(false); - bool old_allow_remote_element_removal = + const bool old_allow_renumbering = mesh_clone->allow_renumbering(); + const bool old_allow_remote_element_removal = mesh_clone->allow_remote_element_removal(); - bool old_skip_partitioning = mesh_clone->skip_partitioning(); - mesh_clone->skip_partitioning(true); + const bool old_skip_partitioning = mesh_clone->skip_partitioning(); + mesh_clone->allow_renumbering(false); mesh_clone->allow_remote_element_removal(false); - mesh_clone->prepare_for_use(); + mesh_clone->skip_partitioning(true); + + // If the mesh thinks it's already completely prepared, test that + if (prep) + mesh_clone->prepare_for_use(); + // If the mesh thinks it's somewhat prepared, test each way it + // thinks so. + else + { + if (prep.has_synched_id_counts) + mesh_clone->update_parallel_id_counts(); + + if (prep.has_neighbor_ptrs) + mesh_clone->find_neighbors(); + + if (prep.has_cached_elem_data) + mesh_clone->cache_elem_data(); + + if (prep.has_interior_parent_ptrs) + mesh_clone->detect_interior_parents(); + + if (old_allow_remote_element_removal && + prep.has_removed_remote_elements) + mesh_clone->delete_remote_elements(); + + if (prep.has_removed_orphaned_nodes) + mesh_clone->remove_orphaned_nodes(); + + if (prep.has_boundary_id_sets) + mesh_clone->get_boundary_info().regenerate_id_sets(); + + // I don't know how we'll tell if this changes anything, but + // we'll do it for completeness + if (prep.has_reinit_ghosting_functors) + mesh_clone->reinit_ghosting_functors(); + } + mesh_clone->allow_renumbering(old_allow_renumbering); mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal); mesh_clone->skip_partitioning(old_skip_partitioning); From 528a29800b0971ea12ac03485f37e12dc370824a Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:13:42 -0600 Subject: [PATCH 06/48] Update comments --- include/mesh/mesh_base.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index 8aaa3f944cd..c3736034453 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -519,16 +519,17 @@ class MeshBase : public ParallelObject * higher dimensions is checked. Also, x-z and y-z planar meshes are * considered to have spatial dimension == 3. * - * The spatial dimension is updated during prepare_for_use() based + * The spatial dimension is updated during mesh preparation based * on the dimensions of the various elements present in the Mesh, - * but is *never automatically decreased* by this function. + * but is *never automatically decreased*. * * For example, if the user calls set_spatial_dimension(2) and then * later inserts 3D elements into the mesh, * Mesh::spatial_dimension() will return 3 after the next call to - * prepare_for_use(). On the other hand, if the user calls - * set_spatial_dimension(3) and then inserts only x-aligned 1D - * elements into the Mesh, mesh.spatial_dimension() will remain 3. + * prepare_for_use() or complete_preparation(). On the other hand, + * if the user calls set_spatial_dimension(3) and then inserts only + * x-aligned 1D elements into the Mesh, mesh.spatial_dimension() + * will remain 3. */ unsigned int spatial_dimension () const; From 88b85b0aef5b59619ad4109111a50822ce80a2d4 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:14:07 -0600 Subject: [PATCH 07/48] Don't try to use an unprepared mesh in testing This wasn't failing the test, because the kernel here didn't need the unprepared data, but we'd like to simply assert that we're not assembling or solving on unprepared meshes. --- tests/systems/periodic_bc_test.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/systems/periodic_bc_test.C b/tests/systems/periodic_bc_test.C index f7c4b18051b..f4705a8d016 100644 --- a/tests/systems/periodic_bc_test.C +++ b/tests/systems/periodic_bc_test.C @@ -215,7 +215,7 @@ private: // Okay, *now* the mesh knows to save periodic neighbors mesh.allow_remote_element_removal(true); - mesh.delete_remote_elements(); + mesh.prepare_for_use(); const unsigned int u_var = sys.add_variable("u", fe_type); sys.attach_assemble_function (periodic_bc_test_poisson); From 54c20a47cf1627690ba200d42ec621d5f252e29a Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:15:17 -0600 Subject: [PATCH 08/48] cache_elem_data() after modifying mesh subdomains Otherwise the mesh still thinks its fully prepared even when it isn't and we can't assert valid_is_prepared() on it. --- tests/mesh/mesh_smoother_test.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index cfe220b2820..e828e9f6533 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -592,6 +592,10 @@ public: mesh.comm().sum(distorted_subdomain_volumes[sub_id]); } + + // We've just invalidated the get_mesh_subdomains() cache by + // adding a new one; fix it. + mesh.cache_elem_data(); } // Get the mesh order From bff2ec6a0dc230fd7195532a4ebd3b913ed9ec9f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:16:33 -0600 Subject: [PATCH 09/48] Prepare mesh_copy before smoother solves on it We need things like neighbor pointers, and in general we'd like to assert that we're never assembling or solving on unprepared meshes. --- src/mesh/mesh_smoother_vsmoother.C | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/mesh/mesh_smoother_vsmoother.C b/src/mesh/mesh_smoother_vsmoother.C index 5cbc8429962..5efb135c9ed 100644 --- a/src/mesh/mesh_smoother_vsmoother.C +++ b/src/mesh/mesh_smoother_vsmoother.C @@ -88,6 +88,16 @@ void VariationalMeshSmoother::setup() // Create a new mesh, EquationSystems, and System _mesh_copy = std::make_unique(_mesh); + + // If the _mesh wasn't prepared, that's fine (we'll just be moving + // its nodes), but we do need the copy to be prepared before our + // solve does things like looking at neighbors. We'll disable + // repartitioning and renumbering first to make sure that we can + // transfer our geometry changes back to the original mesh. + _mesh_copy->allow_renumbering(false); + _mesh_copy->skip_partitioning(true); + _mesh_copy->complete_preparation(); + _equation_systems = std::make_unique(*_mesh_copy); _system = &(_equation_systems->add_system("variational_smoother_system")); From 2e7cbeab6405b535aeeed24c5ebce710c6bb4a45 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:17:40 -0600 Subject: [PATCH 10/48] Assert that we only assemble on prepared meshes We're trying to allow meshes to keep track of "half-prepared" status, for the sake of efficiency in large chains or webs of mesh generators, but at the very least user assembly kernels should be able to rely on a mesh being prepared for use - they're what we mean by "use"!. --- src/systems/fem_system.C | 16 ++++++++++++++++ src/systems/nonlinear_implicit_system.C | 6 ++++++ src/systems/system.C | 16 ++++++++++++++++ 3 files changed, 38 insertions(+) diff --git a/src/systems/fem_system.C b/src/systems/fem_system.C index 70cd06b6fac..90f4292bda9 100644 --- a/src/systems/fem_system.C +++ b/src/systems/fem_system.C @@ -25,6 +25,7 @@ #include "libmesh/fem_system.h" #include "libmesh/libmesh_logging.h" #include "libmesh/mesh_base.h" +#include "libmesh/mesh_tools.h" #include "libmesh/numeric_vector.h" #include "libmesh/parallel_algebra.h" #include "libmesh/parallel_ghost_sync.h" @@ -898,6 +899,11 @@ void FEMSystem::assembly (bool get_residual, bool get_jacobian, const MeshBase & mesh = this->get_mesh(); + libmesh_assert(mesh.is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(mesh)); +#endif + if (print_solution_norms) { this->solution->close(); @@ -1150,6 +1156,11 @@ void FEMSystem::assemble_qoi (const QoISet & qoi_indices) const MeshBase & mesh = this->get_mesh(); + libmesh_assert(mesh.is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(mesh)); +#endif + this->update(); const unsigned int Nq = this->n_qois(); @@ -1184,6 +1195,11 @@ void FEMSystem::assemble_qoi_derivative (const QoISet & qoi_indices, const MeshBase & mesh = this->get_mesh(); + libmesh_assert(mesh.is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(mesh)); +#endif + this->update(); // The quantity of interest derivative assembly accumulates on diff --git a/src/systems/nonlinear_implicit_system.C b/src/systems/nonlinear_implicit_system.C index f8273a11974..a4b078d735a 100644 --- a/src/systems/nonlinear_implicit_system.C +++ b/src/systems/nonlinear_implicit_system.C @@ -22,6 +22,7 @@ #include "libmesh/diff_solver.h" #include "libmesh/equation_systems.h" #include "libmesh/libmesh_logging.h" +#include "libmesh/mesh_tools.h" #include "libmesh/nonlinear_solver.h" #include "libmesh/sparse_matrix.h" #include "libmesh/static_condensation.h" @@ -247,6 +248,11 @@ void NonlinearImplicitSystem::assembly(bool get_residual, bool /*apply_heterogeneous_constraints*/, bool /*apply_no_constraints*/) { + libmesh_assert(this->get_mesh().is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); +#endif + // Get current_local_solution in sync this->update(); diff --git a/src/systems/system.C b/src/systems/system.C index e2a7f0f30d1..64e81413533 100644 --- a/src/systems/system.C +++ b/src/systems/system.C @@ -23,6 +23,7 @@ #include "libmesh/int_range.h" #include "libmesh/libmesh_logging.h" #include "libmesh/mesh_base.h" +#include "libmesh/mesh_tools.h" #include "libmesh/numeric_vector.h" #include "libmesh/parameter_vector.h" #include "libmesh/point.h" // For point_value @@ -551,6 +552,11 @@ void System::assemble () // Log how long the user's assembly code takes LOG_SCOPE("assemble()", "System"); + libmesh_assert(this->get_mesh().is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); +#endif + // Call the user-specified assembly function this->user_assembly(); } @@ -562,6 +568,11 @@ void System::assemble_qoi (const QoISet & qoi_indices) // Log how long the user's assembly code takes LOG_SCOPE("assemble_qoi()", "System"); + libmesh_assert(this->get_mesh().is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); +#endif + // Call the user-specified quantity of interest function this->user_QOI(qoi_indices); } @@ -575,6 +586,11 @@ void System::assemble_qoi_derivative(const QoISet & qoi_indices, // Log how long the user's assembly code takes LOG_SCOPE("assemble_qoi_derivative()", "System"); + libmesh_assert(this->get_mesh().is_prepared()); +#ifdef DEBUG + libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); +#endif + // Call the user-specified quantity of interest function this->user_QOI_derivative(qoi_indices, include_liftfunc, apply_constraints); From 9c49cca877c60e5b26f1afe1c5d6edfd34b0593b Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:32:10 -0600 Subject: [PATCH 11/48] MeshModification shouldn't leave invalid caches We can rebuild point locators or element caches, but we don't want to use an already-built invalid cache by accident. --- src/mesh/mesh_modification.C | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh_modification.C b/src/mesh/mesh_modification.C index a7429475e58..6015b813854 100644 --- a/src/mesh/mesh_modification.C +++ b/src/mesh/mesh_modification.C @@ -215,6 +215,10 @@ void MeshTools::Modification::distort (MeshBase & mesh, #endif } } + + // We haven't changed any topology, but just changing geometry could + // have invalidated a point locator. + mesh.clear_point_locator(); } @@ -302,6 +306,10 @@ void MeshTools::Modification::redistribute (MeshBase & mesh, (*node)(2) = output_vec(2); #endif } + + // We haven't changed any topology, but just changing geometry could + // have invalidated a point locator. + mesh.clear_point_locator(); } @@ -315,6 +323,10 @@ void MeshTools::Modification::translate (MeshBase & mesh, for (auto & node : mesh.node_ptr_range()) *node += p; + + // We haven't changed any topology, but just changing geometry could + // have invalidated a point locator. + mesh.clear_point_locator(); } @@ -346,6 +358,10 @@ MeshTools::Modification::rotate (MeshBase & mesh, const Real theta, const Real psi) { + // We won't change any topology, but just changing geometry could + // invalidate a point locator. + mesh.clear_point_locator(); + #if LIBMESH_DIM == 3 const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi); @@ -402,6 +418,10 @@ void MeshTools::Modification::scale (MeshBase & mesh, for (auto & node : mesh.node_ptr_range()) (*node)(2) *= z_scale; + + // We haven't changed any topology, but just changing geometry could + // have invalidated a point locator. + mesh.clear_point_locator(); } @@ -1425,8 +1445,6 @@ void MeshTools::Modification::all_tri (MeshBase & mesh) MeshCommunication().make_nodes_parallel_consistent (mesh); } - - // Prepare the newly created mesh for use. mesh.prepare_for_use(); @@ -1585,6 +1603,10 @@ void MeshTools::Modification::smooth (MeshBase & mesh, } } // refinement_level loop } // end iteration + + // We haven't changed any topology, but just changing geometry could + // have invalidated a point locator. + mesh.clear_point_locator(); } @@ -1741,6 +1763,10 @@ void MeshTools::Modification::change_subdomain_id (MeshBase & mesh, if (elem->subdomain_id() == old_id) elem->subdomain_id() = new_id; } + + // We just invalidated mesh.get_subdomain_ids(), but it might not be + // efficient to fix that here. + mesh.set_hasnt_cached_elem_data(); } From 2c17ef283cfa8f82fb57e18ac45a68c2ad893938 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 13:26:54 -0600 Subject: [PATCH 12/48] Assign _children_on_boundary in BoundaryInfo = --- src/mesh/boundary_info.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C index cda463c554d..d7cfe80281c 100644 --- a/src/mesh/boundary_info.C +++ b/src/mesh/boundary_info.C @@ -138,6 +138,8 @@ BoundaryInfo & BoundaryInfo::operator=(const BoundaryInfo & other_boundary_info) for (const auto & [elem, id_pair] : other_boundary_info._boundary_side_id) _boundary_side_id.emplace(_mesh->elem_ptr(elem->id()), id_pair); + _children_on_boundary = other_boundary_info._children_on_boundary; + _boundary_ids = other_boundary_info._boundary_ids; _global_boundary_ids = other_boundary_info._global_boundary_ids; _side_boundary_ids = other_boundary_info._side_boundary_ids; From b6a7832fac53328d673a4249f38db1b5eec966a0 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 13:27:21 -0600 Subject: [PATCH 13/48] Clarify prepare_for_use() doxygen --- include/mesh/mesh_base.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index c3736034453..f0e659c640d 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -1256,6 +1256,13 @@ class MeshBase : public ParallelObject * 10.) removing any remote elements (if enabled) * 11.) regenerating summarized boundary id sets * + * For backwards compatibility, prepare_for_use() performs *all* those + * steps, regardless of the official preparation() state of the + * mesh. In codes which have maintained a valid preparation() state + * via methods such as set_hasnt_synched_id_counts(), calling + * complete_preparation() will result in a fully-prepared mesh at + * less cost. + * * The argument to skip renumbering is now deprecated - to prevent a * mesh from being renumbered, set allow_renumbering(false). The argument to skip * finding neighbors is also deprecated. To prevent find_neighbors, set From a315156a3525ebbd788dfbaafa50677905b4bc12 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 20:55:01 -0600 Subject: [PATCH 14/48] Get interior parent pointers *right* when copying We've added some tricky interior_mesh options that our copy_nodes_and_elements code wasn't handling correctly, and it was possible for a clone() of such a mesh to fail its operator== assertion. --- src/mesh/unstructured_mesh.C | 50 +++++++++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C index d3a09d63cec..ac9da962c8e 100644 --- a/src/mesh/unstructured_mesh.C +++ b/src/mesh/unstructured_mesh.C @@ -735,7 +735,7 @@ void UnstructuredMesh::copy_nodes_and_elements(const MeshBase & other_mesh, // Hold off on trying to set the interior parent because we may actually // add lower dimensional elements before their interior parents - if (el->dim() < other_mesh.mesh_dimension() && old->interior_parent()) + if (old->interior_parent()) ip_map[old] = el.get(); #ifdef LIBMESH_ENABLE_AMR @@ -797,9 +797,51 @@ void UnstructuredMesh::copy_nodes_and_elements(const MeshBase & other_mesh, } } - for (auto & elem_pair : ip_map) - elem_pair.second->set_interior_parent( - this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset)); + // If the other_mesh had some interior parents, we may need to + // copy those pointers (if they're to elements in a third mesh), + // or create new equivalent pointers (if they're to elements we + // just copied), or scream and die (if the other mesh had interior + // parents from a third mesh but we already have interior parents + // that aren't to that same third mesh. + if (!ip_map.empty()) + { + bool existing_interior_parents = false; + for (const auto & elem : this->element_ptr_range()) + if (elem->interior_parent()) + { + existing_interior_parents = true; + break; + } + + MeshBase * other_interior_mesh = + const_cast(&other_mesh.interior_mesh()); + + // If we don't already have interior parents, then we can just + // use whatever interior_mesh we need for the incoming + // elements. + if (!existing_interior_parents) + { + if (other_interior_mesh == &other_mesh) + this->set_interior_mesh(*this); + else + this->set_interior_mesh(*other_interior_mesh); + } + + if (other_interior_mesh == &other_mesh && + _interior_mesh == this) + for (auto & elem_pair : ip_map) + elem_pair.second->set_interior_parent( + this->elem_ptr(elem_pair.first->interior_parent()->id() + element_id_offset)); + else if (other_interior_mesh == _interior_mesh) + for (auto & elem_pair : ip_map) + { + Elem * ip = const_cast(elem_pair.first->interior_parent()); + libmesh_assert(ip == other_interior_mesh->elem_ptr(ip->id())); + elem_pair.second->set_interior_parent(ip); + } + else + libmesh_error_msg("Cannot copy boundary elements between meshes with different interior meshes"); + } // Loop (again) over the elements to fill in the neighbors if (skip_find_neighbors) From b79b7093ca95eebfd01ebf591c139847a536fe3e Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 20:56:43 -0600 Subject: [PATCH 15/48] More testing of preparation() validity in dbg mode --- src/mesh/mesh_base.C | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 9a4301d3dee..610a96eb828 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -860,10 +860,14 @@ void MeshBase::complete_preparation() libmesh_assert(this->comm().verify(this->is_serial())); +#ifdef DEBUG // If we don't go into this method with valid constraint rows, we're // only going to be able to make that worse. -#ifdef DEBUG MeshTools::libmesh_assert_valid_constraint_rows(*this); + + // If this mesh thinks it's already partially prepared, then in + // optimized builds we'll trust it, but in debug builds we'll check. + const bool was_partly_prepared = (_preparation == Preparation()); #endif // A distributed mesh may have processors with no elements (or @@ -971,6 +975,10 @@ void MeshBase::complete_preparation() libmesh_assert(_preparation); #ifdef DEBUG + // The if() here avoids both unnecessary work *and* stack overflow + if (was_partly_prepared) + libmesh_assert(MeshTools::valid_is_prepared(*this)); + MeshTools::libmesh_assert_valid_boundary_ids(*this); #ifdef LIBMESH_ENABLE_UNIQUE_ID MeshTools::libmesh_assert_valid_unique_ids(*this); From a63efcbe7124414994a4c4d926c6853d50548877 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 5 Dec 2025 17:23:17 -0600 Subject: [PATCH 16/48] Update bcid sets with new bcids in eigen ex3 --- examples/eigenproblems/eigenproblems_ex3/eigenproblems_ex3.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/eigenproblems/eigenproblems_ex3/eigenproblems_ex3.C b/examples/eigenproblems/eigenproblems_ex3/eigenproblems_ex3.C index 75f68ba8eed..65a0dca5f2a 100644 --- a/examples/eigenproblems/eigenproblems_ex3/eigenproblems_ex3.C +++ b/examples/eigenproblems/eigenproblems_ex3/eigenproblems_ex3.C @@ -170,6 +170,8 @@ int main (int argc, char ** argv) if (elem->neighbor_ptr (side) == nullptr) mesh.get_boundary_info().add_side(elem, side, BOUNDARY_ID); + mesh.get_boundary_info().regenerate_id_sets(); + // Print information about the mesh to the screen. mesh.print_info(); From ee0d4076bd2cd0740a6ffe131dc4efa56a9aff5c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 8 Dec 2025 17:55:02 -0600 Subject: [PATCH 17/48] Update boundary id sets after boundary update --- examples/fem_system/fem_system_ex3/fem_system_ex3.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/fem_system/fem_system_ex3/fem_system_ex3.C b/examples/fem_system/fem_system_ex3/fem_system_ex3.C index 820d829232a..1ca8f370592 100644 --- a/examples/fem_system/fem_system_ex3/fem_system_ex3.C +++ b/examples/fem_system/fem_system_ex3/fem_system_ex3.C @@ -173,6 +173,8 @@ int main (int argc, char ** argv) mesh.get_boundary_info().add_edge(elem, e, edge_boundary_id); } + mesh.get_boundary_info().regenerate_id_sets(); + // Create an equation systems object. EquationSystems equation_systems (mesh); From a2ba524736fe13f3e63c5ca8103c23499883b28a Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 13:38:03 -0600 Subject: [PATCH 18/48] Remove unnecessary side_id_map computation This must have been an atavism from a previous refactoring. --- src/mesh/boundary_info.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C index d7cfe80281c..0eb7511e660 100644 --- a/src/mesh/boundary_info.C +++ b/src/mesh/boundary_info.C @@ -515,9 +515,8 @@ void BoundaryInfo::sync (const std::set & requested_boundary_i boundary_mesh.set_n_partitions() = _mesh->n_partitions(); std::map node_id_map; - std::map, dof_id_type> side_id_map; - this->_find_id_maps(requested_boundary_ids, 0, &node_id_map, 0, &side_id_map, subdomains_relative_to); + this->_find_id_maps(requested_boundary_ids, 0, &node_id_map, 0, nullptr, subdomains_relative_to); // Let's add all the boundary nodes we found to the boundary mesh for (const auto & node : _mesh->node_ptr_range()) From 970142a90aef80545f46df8dc18b1500695469c1 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 13:38:31 -0600 Subject: [PATCH 19/48] Fix side_id_map generation on refined meshes We never caught this before because we were never doing much subsequent mesh modification to the output, but as soon as we start heavier testing we run into code that screams about the inconsistent parent/child id ordering. --- src/mesh/boundary_info.C | 60 ++++++++++++++++++++++++++++++++-------- 1 file changed, 49 insertions(+), 11 deletions(-) diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C index 0eb7511e660..1e39b3d5755 100644 --- a/src/mesh/boundary_info.C +++ b/src/mesh/boundary_info.C @@ -3300,11 +3300,28 @@ void BoundaryInfo::_find_id_maps(const std::set & requested_bo std::map, dof_id_type> * side_id_map, const std::set & subdomains_relative_to) { - // We'll do the same modulus trick that DistributedMesh uses to avoid - // id conflicts between different processors + // We'll use the standard DistributedMesh trick of dividing up id + // ranges by processor_id to avoid picking duplicate node ids. dof_id_type - next_node_id = first_free_node_id + this->processor_id(), - next_elem_id = first_free_elem_id + this->processor_id(); + next_node_id = first_free_node_id + this->processor_id(); + + // We have to be careful how we select different element ids on + // different processors, because we may be adding sides of refined + // elements and sides of their parents (which are thereby also + // parents of their sides), but libMesh convention requires (and + // some libMesh communication code assumes) that parent elements + // have lower ids than their children. + // + // We'll start with temporary ids in a temporary map stratefied by + // element level, so we can sort by element level after. + // Make sure we're sorting by ids here, not anything that might + // differ from processor to processor, so we can do an + // embarrassingly parallel sort. This is a serialized data + // structure, but it's at worst O(N^2/3) so we'll hold off on doing + // this distributed until someone's benchmark screams at us. + std::vector + >> + side_id_set_by_level; // For avoiding extraneous element side construction ElemSideBuilder side_builder; @@ -3337,14 +3354,19 @@ void BoundaryInfo::_find_id_maps(const std::set & requested_bo // Join up the local results from other processors if (side_id_map) - this->comm().set_union(*side_id_map); + { + std::size_t n_levels = side_id_set_by_level.size(); + this->comm().max(n_levels); + side_id_set_by_level.resize(n_levels); + for (auto l : make_range(n_levels)) + this->comm().set_union(side_id_set_by_level[l]); + } if (node_id_map) this->comm().set_union(*node_id_map); // Finally we'll pass through any unpartitioned elements to add them // to the maps and counts. next_node_id = first_free_node_id + this->n_processors(); - next_elem_id = first_free_elem_id + this->n_processors(); el = _mesh->pid_elements_begin(DofObject::invalid_processor_id); if (el == end_unpartitioned_el) @@ -3401,9 +3423,12 @@ void BoundaryInfo::_find_id_maps(const std::set & requested_bo DofObject::invalid_processor_id))) { std::pair side_pair(elem->id(), s); - libmesh_assert (!side_id_map->count(side_pair)); - (*side_id_map)[side_pair] = next_elem_id; - next_elem_id += this->n_processors() + 1; + auto level = elem->level(); + if (side_id_set_by_level.size() <= level) + side_id_set_by_level.resize(level+1); + auto & level_side_id_set = side_id_set_by_level[level]; + libmesh_assert (!level_side_id_set.count(side_pair)); + level_side_id_set.insert(side_pair); } side = &side_builder(*elem, s); @@ -3428,8 +3453,21 @@ void BoundaryInfo::_find_id_maps(const std::set & requested_bo } } - // FIXME: ought to renumber side/node_id_map image to be contiguous - // to save memory, also ought to reserve memory + // FIXME: should we renumber node_id_map's image to be contiguous + // here, rather than waiting for renumbering in mesh preparation to + // do it later? + + // We do need to build up side_id_map here, to handle proper sorting + // by level. + if (side_id_map) + { + dof_id_type next_elem_id = first_free_elem_id; + for (auto level : make_range(side_id_set_by_level.size())) + { + for (auto side_pair : side_id_set_by_level[level]) + (*side_id_map)[side_pair] = next_elem_id++; + } + } } void BoundaryInfo::clear_stitched_boundary_side_ids (const boundary_id_type sideset_id, From 4aab53dedc7b8d45b97d07f213e4dc06ad02df32 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 21:46:30 -0600 Subject: [PATCH 20/48] Added MeshBase::assert_equal_to It's quite annoying to me to have to hunt down *why* a valid_is_prepared() call fails, which means it's got to be utterly intractable to too many users. So, instead of just asserting the final boolean results of an operator==, let's provide some way to get assertion failures that deliver more information about what wasn't equal. --- include/mesh/distributed_mesh.h | 2 +- include/mesh/mesh_base.h | 18 +++- include/mesh/replicated_mesh.h | 2 +- src/mesh/distributed_mesh.C | 46 +++++----- src/mesh/mesh_base.C | 156 +++++++++++++++++++------------- src/mesh/replicated_mesh.C | 20 ++-- 6 files changed, 146 insertions(+), 98 deletions(-) diff --git a/include/mesh/distributed_mesh.h b/include/mesh/distributed_mesh.h index 8f610c9b139..cf145508fdc 100644 --- a/include/mesh/distributed_mesh.h +++ b/include/mesh/distributed_mesh.h @@ -111,7 +111,7 @@ class DistributedMesh : public UnstructuredMesh * Shim to allow operator == (&) to behave like a virtual function * without having to be one. */ - virtual bool subclass_locally_equals (const MeshBase & other_mesh) const override; + virtual std::string_view subclass_first_difference_from (const MeshBase & other_mesh_base) const override; /** * Virtual copy-constructor, creates a copy of this mesh diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index f0e659c640d..883ff282348 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -135,7 +135,7 @@ class MeshBase : public ParallelObject * ids). * * Though this method is non-virtual, its implementation calls the - * virtual function \p subclass_locally_equals() to test for + * virtual function \p subclass_first_difference_from() to test for * equality of subclass-specific data as well. */ bool operator== (const MeshBase & other_mesh) const; @@ -145,6 +145,14 @@ class MeshBase : public ParallelObject return !(*this == other_mesh); } + /** + * This behaves like libmesh_assert(*this == other_mesh), but gives + * a more useful accounting of the first difference found, if the + * assertion fails. + */ + void assert_equal_to (const MeshBase & other_mesh, + std::string_view failure_context) const; + /** * This behaves the same as operator==, but only for the local and * ghosted aspects of the mesh; i.e. operator== is true iff local @@ -2078,7 +2086,7 @@ class MeshBase : public ParallelObject * Shim to allow operator == (&) to behave like a virtual function * without having to be one. */ - virtual bool subclass_locally_equals (const MeshBase & other_mesh) const = 0; + virtual std::string_view subclass_first_difference_from (const MeshBase & other_mesh) const = 0; /** * Tests for equality of all elements and nodes in the mesh. Helper @@ -2086,6 +2094,12 @@ class MeshBase : public ParallelObject */ bool nodes_and_elements_equal(const MeshBase & other_mesh) const; + /** + * + */ + std::string_view first_difference_from(const MeshBase & other_mesh) const; + + /** * \returns A writable reference to the number of partitions. */ diff --git a/include/mesh/replicated_mesh.h b/include/mesh/replicated_mesh.h index 9e2d5e950b6..28ddb7eb4dd 100644 --- a/include/mesh/replicated_mesh.h +++ b/include/mesh/replicated_mesh.h @@ -97,7 +97,7 @@ class ReplicatedMesh : public UnstructuredMesh * Shim to allow operator == (&) to behave like a virtual function * without having to be one. */ - virtual bool subclass_locally_equals (const MeshBase & other_mesh) const override; + virtual std::string_view subclass_first_difference_from (const MeshBase & other_mesh_base) const override; /** * Virtual copy-constructor, creates a copy of this mesh diff --git a/src/mesh/distributed_mesh.C b/src/mesh/distributed_mesh.C index 5a15a04814b..4989a71a2a7 100644 --- a/src/mesh/distributed_mesh.C +++ b/src/mesh/distributed_mesh.C @@ -95,48 +95,52 @@ MeshBase & DistributedMesh::assign(MeshBase && other_mesh) return *this; } -bool DistributedMesh::subclass_locally_equals(const MeshBase & other_mesh_base) const +std::string_view DistributedMesh::subclass_first_difference_from (const MeshBase & other_mesh_base) const { const DistributedMesh * dist_mesh_ptr = dynamic_cast(&other_mesh_base); if (!dist_mesh_ptr) - return false; + return "DistributedMesh class"; const DistributedMesh & other_mesh = *dist_mesh_ptr; - if (_is_serial != other_mesh._is_serial || - _is_serial_on_proc_0 != other_mesh._is_serial_on_proc_0 || - _deleted_coarse_elements != other_mesh._deleted_coarse_elements || - _n_nodes != other_mesh._n_nodes || - _n_elem != other_mesh._n_elem || - _max_node_id != other_mesh._max_node_id || - _max_elem_id != other_mesh._max_elem_id || - // We expect these things to change in a prepare_for_use(); - // they're conceptually "mutable"... +#define CHECK_MEMBER(member_name) \ + if (member_name != other_mesh.member_name) \ + return #member_name; + + CHECK_MEMBER(_is_serial); + CHECK_MEMBER(_is_serial_on_proc_0); + CHECK_MEMBER(_deleted_coarse_elements); + CHECK_MEMBER(_n_nodes); + CHECK_MEMBER(_n_elem); + CHECK_MEMBER(_max_node_id); + CHECK_MEMBER(_max_elem_id); + // We expect these things to change in a prepare_for_use(); + // they're conceptually "mutable"... /* - _next_free_local_node_id != other_mesh._next_free_local_node_id || - _next_free_local_elem_id != other_mesh._next_free_local_elem_id || - _next_free_unpartitioned_node_id != other_mesh._next_free_unpartitioned_node_id || - _next_free_unpartitioned_elem_id != other_mesh._next_free_unpartitioned_elem_id || + CHECK_MEMBER(_next_free_local_node_id); + CHECK_MEMBER(_next_free_local_elem_id); + CHECK_MEMBER(_next_free_unpartitioned_node_id); + CHECK_MEMBER(_next_free_unpartitioned_elem_id); #ifdef LIBMESH_ENABLE_UNIQUE_ID - _next_unpartitioned_unique_id != other_mesh._next_unpartitioned_unique_id || + CHECK_MEMBER(_next_unpartitioned_unique_id); #endif */ - !this->nodes_and_elements_equal(other_mesh)) - return false; + if (!this->nodes_and_elements_equal(other_mesh)) + return "nodes and/or elements"; if (_extra_ghost_elems.size() != other_mesh._extra_ghost_elems.size()) - return false; + return "_extra_ghost_elems size"; for (auto & elem : _extra_ghost_elems) { libmesh_assert(this->query_elem_ptr(elem->id()) == elem); const Elem * other_elem = other_mesh.query_elem_ptr(elem->id()); if (!other_elem || !other_mesh._extra_ghost_elems.count(const_cast(other_elem))) - return false; + return "_extra_ghost_elems entry"; } - return true; + return ""; } DistributedMesh::~DistributedMesh () diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 610a96eb828..dea9747a920 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -265,7 +265,53 @@ bool MeshBase::operator== (const MeshBase & other_mesh) const } +void MeshBase::assert_equal_to (const MeshBase & other_mesh, + std::string_view failure_context) const +{ +#ifndef NDEBUG + LOG_SCOPE("operator==()", "MeshBase"); + + std::string_view local_diff = first_difference_from(other_mesh); + + bool diff_found = !local_diff.empty(); + if (diff_found) + { + // Construct a user-friendly message to throw on pid 0 + std::set unique_diffs; + unique_diffs.insert(std::string(local_diff)); + this->comm().set_union(unique_diffs); + + if (!this->processor_id()) + { + std::string error_msg {failure_context}; + error_msg += '\n' + + "Meshes failed asserted equality in at least these aspects:\n"; + for (auto & diff : unique_diffs) + { + error_msg += diff; + error_msg += '\n'; + } + libmesh_assert_msg(!diff_found, error_msg); + } + + // We're not going to throw on other processors because we don't + // want to accidentally preempt pid 0's error message. We're + // not even going to exit on other processors because for all we + // know user code is going to catch that error and sync up with + // us later. + } +#endif // NDEBUG +} + + bool MeshBase::locally_equals (const MeshBase & other_mesh) const +{ + const std::string_view diff = first_difference_from(other_mesh); + return diff.empty(); +} + + +std::string_view MeshBase::first_difference_from(const MeshBase & other_mesh) const { // Check whether (almost) everything in the base is equal // @@ -273,21 +319,19 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const // change in a DistributedMesh prepare_for_use(); it's conceptually // "mutable". // - // We use separate if statements instead of logical operators here, - // to make it easy to see the failing condition when using a - // debugger to figure out why a MeshTools::valid_is_prepared(mesh) - // is failing. - if (_n_parts != other_mesh._n_parts) - return false; - if (_default_mapping_type != other_mesh._default_mapping_type) - return false; - if (_default_mapping_data != other_mesh._default_mapping_data) - return false; - if (_preparation != other_mesh._preparation) - return false; - if (_count_lower_dim_elems_in_point_locator != - other_mesh._count_lower_dim_elems_in_point_locator) - return false; + // We use separate tests here and return strings for each test, + // to make it easy to see the failing condition a + // MeshTools::valid_is_prepared(mesh) is failing. + +#define CHECK_MEMBER(member_name) \ + if (member_name != other_mesh.member_name) \ + return #member_name; + + CHECK_MEMBER(_n_parts); + CHECK_MEMBER(_default_mapping_type); + CHECK_MEMBER(_default_mapping_data); + CHECK_MEMBER(_preparation); + CHECK_MEMBER(_count_lower_dim_elems_in_point_locator); // We should either both have our own interior parents or both not; // but if we both don't then we can't really assert anything else @@ -295,59 +339,41 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const // pointing at two different copies of "the same" interior mesh. if ((_interior_mesh == this) != (other_mesh._interior_mesh == &other_mesh)) - return false; - - if (_skip_noncritical_partitioning != other_mesh._skip_noncritical_partitioning) - return false; - if (_skip_all_partitioning != other_mesh._skip_all_partitioning) - return false; - if (_skip_renumber_nodes_and_elements != other_mesh._skip_renumber_nodes_and_elements) - return false; - if (_skip_find_neighbors != other_mesh._skip_find_neighbors) - return false; - if (_allow_remote_element_removal != other_mesh._allow_remote_element_removal) - return false; - if (_allow_node_and_elem_unique_id_overlap != other_mesh._allow_node_and_elem_unique_id_overlap) - return false; - if (_spatial_dimension != other_mesh._spatial_dimension) - return false; - if (_point_locator_close_to_point_tol != other_mesh._point_locator_close_to_point_tol) - return false; - if (_block_id_to_name != other_mesh._block_id_to_name) - return false; - if (_elem_dims != other_mesh._elem_dims) - return false; - if (_elem_default_orders != other_mesh._elem_default_orders) - return false; - if (_supported_nodal_order != other_mesh._supported_nodal_order) - return false; - if (_mesh_subdomains != other_mesh._mesh_subdomains) - return false; - if (_all_elemset_ids != other_mesh._all_elemset_ids) - return false; - if (_elem_integer_names != other_mesh._elem_integer_names) - return false; - if (_elem_integer_default_values != other_mesh._elem_integer_default_values) - return false; - if (_node_integer_names != other_mesh._node_integer_names) - return false; - if (_node_integer_default_values != other_mesh._node_integer_default_values) - return false; + return "_interior_mesh"; + + CHECK_MEMBER(_skip_noncritical_partitioning); + CHECK_MEMBER(_skip_all_partitioning); + CHECK_MEMBER(_skip_renumber_nodes_and_elements); + CHECK_MEMBER(_skip_find_neighbors); + CHECK_MEMBER(_allow_remote_element_removal); + CHECK_MEMBER(_allow_node_and_elem_unique_id_overlap); + CHECK_MEMBER(_spatial_dimension); + CHECK_MEMBER(_point_locator_close_to_point_tol); + CHECK_MEMBER(_block_id_to_name); + CHECK_MEMBER(_elem_dims); + CHECK_MEMBER(_elem_default_orders); + CHECK_MEMBER(_supported_nodal_order); + CHECK_MEMBER(_mesh_subdomains); + CHECK_MEMBER(_all_elemset_ids); + CHECK_MEMBER(_elem_integer_names); + CHECK_MEMBER(_elem_integer_default_values); + CHECK_MEMBER(_node_integer_names); + CHECK_MEMBER(_node_integer_default_values); if (static_cast(_default_ghosting) != static_cast(other_mesh._default_ghosting)) - return false; + return "_default_ghosting"; if (static_cast(_partitioner) != static_cast(other_mesh._partitioner)) - return false; + return "_partitioner"; if (*boundary_info != *other_mesh.boundary_info) - return false; + return "boundary_info"; // First check whether the "existence" of the two pointers differs (one present, one absent) if (static_cast(_disjoint_neighbor_boundary_pairs) != static_cast(other_mesh._disjoint_neighbor_boundary_pairs)) - return false; + return "_disjoint_neighbor_boundary_pairs existence"; // If both exist, compare the contents (Weak Test: just compare sizes like `_ghosting_functors`) if (_disjoint_neighbor_boundary_pairs && (_disjoint_neighbor_boundary_pairs->size() != other_mesh._disjoint_neighbor_boundary_pairs->size())) - return false; + return "_disjoint_neighbor_boundary_pairs size"; const constraint_rows_type & other_rows = other_mesh.get_constraint_rows(); @@ -356,15 +382,15 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const const dof_id_type node_id = node->id(); const Node * other_node = other_mesh.query_node_ptr(node_id); if (!other_node) - return false; + return "_constraint_rows node presence"; auto it = other_rows.find(other_node); if (it == other_rows.end()) - return false; + return "_constraint_rows row presence"; const auto & other_row = it->second; if (row.size() != other_row.size()) - return false; + return "_constraint_rows row size"; for (auto i : index_range(row)) { @@ -377,14 +403,14 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const elem_pair.second != other_elem_pair.second || coef != other_coef) - return false; + return "_constraint_rows row entry"; } } for (const auto & [elemset_code, elemset_ptr] : this->_elemset_codes) if (const auto it = other_mesh._elemset_codes.find(elemset_code); it == other_mesh._elemset_codes.end() || *elemset_ptr != *it->second) - return false; + return "_elemset_codes"; // FIXME: we have no good way to compare ghosting functors, since // they're in a vector of pointers, and we have no way *at all* @@ -393,13 +419,13 @@ bool MeshBase::locally_equals (const MeshBase & other_mesh) const // we have the same number, is all. if (_ghosting_functors.size() != other_mesh._ghosting_functors.size()) - return false; + return "_ghosting_functors size"; // Same deal for partitioners. We tested that we both have one or // both don't, but are they equivalent? Let's guess "yes". // Now let the subclasses decide whether everything else is equal - return this->subclass_locally_equals(other_mesh); + return this->subclass_first_difference_from(other_mesh); } diff --git a/src/mesh/replicated_mesh.C b/src/mesh/replicated_mesh.C index b58d941156b..b54477ca1b6 100644 --- a/src/mesh/replicated_mesh.C +++ b/src/mesh/replicated_mesh.C @@ -59,23 +59,27 @@ ReplicatedMesh::ReplicatedMesh (const Parallel::Communicator & comm_in, } -bool ReplicatedMesh::subclass_locally_equals(const MeshBase & other_mesh_base) const +std::string_view ReplicatedMesh::subclass_first_difference_from (const MeshBase & other_mesh_base) const { const ReplicatedMesh * rep_mesh_ptr = dynamic_cast(&other_mesh_base); if (!rep_mesh_ptr) - return false; + return "ReplicatedMesh class"; const ReplicatedMesh & other_mesh = *rep_mesh_ptr; - if (_n_nodes != other_mesh._n_nodes || - _n_elem != other_mesh._n_elem || +#define CHECK_MEMBER(member_name) \ + if (member_name != other_mesh.member_name) \ + return #member_name; + + CHECK_MEMBER(_n_nodes); + CHECK_MEMBER(_n_elem); #ifdef LIBMESH_ENABLE_UNIQUE_ID - _next_unique_id != other_mesh._next_unique_id || + CHECK_MEMBER(_next_unique_id); #endif - !this->nodes_and_elements_equal(other_mesh)) - return false; + if (!this->nodes_and_elements_equal(other_mesh)) + return "nodes and/or elements"; - return true; + return ""; } From b5ebe71c6cbffb3cfe221e14b039122c15298ee2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 21:50:23 -0600 Subject: [PATCH 21/48] Add MeshTools::libmesh_assert_valid_is_prepared() This should give a more informative assertion failure message when needed. --- include/mesh/mesh_tools.h | 21 +++++-- src/mesh/mesh_tools.C | 128 ++++++++++++++++++++++++++------------ 2 files changed, 103 insertions(+), 46 deletions(-) diff --git a/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h index 04ea94fba0b..c73645f2391 100644 --- a/include/mesh/mesh_tools.h +++ b/include/mesh/mesh_tools.h @@ -394,11 +394,15 @@ void clear_spline_nodes(MeshBase &); /** * A function for testing whether a mesh's cached is_prepared() setting - * is not a false positive. If the mesh is marked as not prepared, or - * if preparing the already-partitioned mesh (without any - * repartitioning or renumbering) does not change it, then we return - * true. If the mesh believes it is prepared but prepare_for_use() - * would change it, we return false. + * is not a false positive, and that its preparation() values are not + * false positives. If the mesh is marked as completely not prepared, or + * if preparing a already-prepared mesh (without any repartitioning or + * renumbering) or preparing after a complete_preparation() (likewise) + * does not change the mesh, then we return true. If the mesh believes it + * is prepared but prepare_for_use() would change it, or if the mesh + * believes it is partially prepared in a certain way but + * complete_preparation() does not completely prepare it, we return + * false. */ bool valid_is_prepared (const MeshBase & mesh); @@ -412,6 +416,13 @@ bool valid_is_prepared (const MeshBase & mesh); #ifndef NDEBUG +/** + * Like libmesh_assert(MeshTools::valid_is_prepared(mesh)), but + * provides more information on preparation incompleteness in case of + * an error. + */ +void libmesh_assert_valid_is_prepared (const MeshBase & mesh); + /** * A function for testing that all DofObjects within a mesh * have the same n_systems count diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C index 0650d3f991f..d5a9ce9f2eb 100644 --- a/src/mesh/mesh_tools.C +++ b/src/mesh/mesh_tools.C @@ -51,7 +51,7 @@ // ------------------------------------------------------------ -// anonymous namespace for helper classes +// anonymous namespace for helper classes and subroutines namespace { using namespace libMesh; @@ -406,6 +406,68 @@ void find_nodal_neighbors_helper(const dof_id_type global_id, neighbors.assign(neighbor_set.begin(), neighbor_set.end()); } + + +std::unique_ptr reprepared_mesh_clone (const MeshBase & mesh) +{ + const MeshBase::Preparation prep = mesh.preparation(); + + // If the mesh thinks it's prepared in some way, *re*-preparing in + // that way shouldn't change a clone of it, as long as we disallow + // repartitioning or renumbering or remote element removal. + std::unique_ptr mesh_clone = mesh.clone(); + + const bool old_allow_renumbering = mesh_clone->allow_renumbering(); + const bool old_allow_remote_element_removal = + mesh_clone->allow_remote_element_removal(); + const bool old_skip_partitioning = mesh_clone->skip_partitioning(); + mesh_clone->allow_renumbering(false); + mesh_clone->allow_remote_element_removal(false); + mesh_clone->skip_partitioning(true); + + // If the mesh thinks it's already completely prepared, test that + if (prep) + mesh_clone->prepare_for_use(); + // If the mesh thinks it's somewhat prepared, test each way it + // thinks so. + else + { + if (prep.has_synched_id_counts) + mesh_clone->update_parallel_id_counts(); + + if (prep.has_neighbor_ptrs) + mesh_clone->find_neighbors(); + + if (prep.has_cached_elem_data) + mesh_clone->cache_elem_data(); + + if (prep.has_interior_parent_ptrs) + mesh_clone->detect_interior_parents(); + + if (old_allow_remote_element_removal && + prep.has_removed_remote_elements) + mesh_clone->delete_remote_elements(); + + if (prep.has_removed_orphaned_nodes) + mesh_clone->remove_orphaned_nodes(); + + if (prep.has_boundary_id_sets) + mesh_clone->get_boundary_info().regenerate_id_sets(); + + // I don't know how we'll tell if this changes anything, but + // we'll do it for completeness + if (prep.has_reinit_ghosting_functors) + mesh_clone->reinit_ghosting_functors(); + } + + mesh_clone->allow_renumbering(old_allow_renumbering); + mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal); + mesh_clone->skip_partitioning(old_skip_partitioning); + + return mesh_clone; +} + + } @@ -1246,61 +1308,45 @@ bool valid_is_prepared (const MeshBase & mesh) // If the mesh thinks it's prepared in some way, *re*-preparing in // that way shouldn't change a clone of it, as long as we disallow // repartitioning or renumbering or remote element removal. - std::unique_ptr mesh_clone = mesh.clone(); + std::unique_ptr mesh_clone = reprepared_mesh_clone(mesh); - const bool old_allow_renumbering = mesh_clone->allow_renumbering(); - const bool old_allow_remote_element_removal = - mesh_clone->allow_remote_element_removal(); - const bool old_skip_partitioning = mesh_clone->skip_partitioning(); - mesh_clone->allow_renumbering(false); - mesh_clone->allow_remote_element_removal(false); - mesh_clone->skip_partitioning(true); + return (mesh == *mesh_clone); +} - // If the mesh thinks it's already completely prepared, test that - if (prep) - mesh_clone->prepare_for_use(); - // If the mesh thinks it's somewhat prepared, test each way it - // thinks so. - else - { - if (prep.has_synched_id_counts) - mesh_clone->update_parallel_id_counts(); - if (prep.has_neighbor_ptrs) - mesh_clone->find_neighbors(); - if (prep.has_cached_elem_data) - mesh_clone->cache_elem_data(); +#ifndef NDEBUG - if (prep.has_interior_parent_ptrs) - mesh_clone->detect_interior_parents(); - if (old_allow_remote_element_removal && - prep.has_removed_remote_elements) - mesh_clone->delete_remote_elements(); +void libmesh_assert_valid_is_prepared (const MeshBase & mesh) +{ + LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools"); - if (prep.has_removed_orphaned_nodes) - mesh_clone->remove_orphaned_nodes(); + const MeshBase::Preparation prep = mesh.preparation(); - if (prep.has_boundary_id_sets) - mesh_clone->get_boundary_info().regenerate_id_sets(); + // If the mesh doesn't think *anything* has been prepared, we have + // nothing to check. + if (prep == MeshBase::Preparation()) + return; - // I don't know how we'll tell if this changes anything, but - // we'll do it for completeness - if (prep.has_reinit_ghosting_functors) - mesh_clone->reinit_ghosting_functors(); - } + // If the mesh thinks it's partitioned, check. These are counts, + // not caches, so it's a real check. + if (prep.is_partitioned) + libmesh_assert_msg + (!mesh.n_unpartitioned_elem() && !mesh.n_unpartitioned_nodes(), + "Mesh data does not match mesh preparation().is_partitioned"); - mesh_clone->allow_renumbering(old_allow_renumbering); - mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal); - mesh_clone->skip_partitioning(old_skip_partitioning); + // If the mesh thinks it's prepared in some way, *re*-preparing in + // that way shouldn't change a clone of it, as long as we disallow + // repartitioning or renumbering or remote element removal. + std::unique_ptr mesh_clone = reprepared_mesh_clone(mesh); - return (mesh == *mesh_clone); + mesh.assert_equal_to(*mesh_clone, "Mesh data does not match mesh preparation()."); } -#ifndef NDEBUG + void libmesh_assert_equal_n_systems (const MeshBase & mesh) { From eba83897489ef1f8576b3b1c072ecd548ac3066d Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 21:53:36 -0600 Subject: [PATCH 22/48] Enable !NDEBUG libmesh_assert_X as libmesh_assert When NDEBUG isn't defined, these should be no-ops, not unavailable, so we don't have to wrap them in ifdefs elsewhere. --- include/mesh/mesh_tools.h | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h index c73645f2391..016b94bee2a 100644 --- a/include/mesh/mesh_tools.h +++ b/include/mesh/mesh_tools.h @@ -409,7 +409,7 @@ bool valid_is_prepared (const MeshBase & mesh); ///@{ /** - * The following functions, only available in builds with NDEBUG + * The following functions, no-ops except in builds with NDEBUG * undefined, are for asserting internal consistency that we hope * should never be broken in opt */ @@ -655,6 +655,39 @@ namespace Private { void globally_renumber_nodes_and_elements (MeshBase &); } // end namespace Private + +// Declare inline no-ops for assertions with NDEBUG defined +#ifdef NDEBUG + +void libmesh_assert_valid_is_prepared (const MeshBase &) {} + +void libmesh_assert_equal_n_systems (const MeshBase &) {} + +void libmesh_assert_old_dof_objects (const MeshBase &) {} + +void libmesh_assert_valid_node_pointers (const MeshBase &) {} + +void libmesh_assert_valid_remote_elems (const MeshBase &) {} + +void libmesh_assert_valid_elem_ids (const MeshBase &) {} + +void libmesh_assert_valid_amr_elem_ids (const MeshBase &) {} + +void libmesh_assert_valid_amr_interior_parents (const MeshBase &) {} + +void libmesh_assert_contiguous_dof_ids (const MeshBase &, unsigned int) {} + +template +void libmesh_assert_topology_consistent_procids (const MeshBase &) {} + +void libmesh_assert_canonical_node_procids (const MeshBase &) {} + +void libmesh_assert_valid_refinement_tree (const MeshBase &) {} + +#endif // NDEBUG + + + } // end namespace MeshTools } // namespace libMesh From 1937cade08fb77d5b01905d7a1d0830f522d8da7 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:16:47 -0600 Subject: [PATCH 23/48] Better libmesh_assert_valid_is_prepared messages This should be a little less cryptic, anyway. --- src/mesh/mesh_tools.C | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C index d5a9ce9f2eb..1ab3d632240 100644 --- a/src/mesh/mesh_tools.C +++ b/src/mesh/mesh_tools.C @@ -1334,14 +1334,17 @@ void libmesh_assert_valid_is_prepared (const MeshBase & mesh) if (prep.is_partitioned) libmesh_assert_msg (!mesh.n_unpartitioned_elem() && !mesh.n_unpartitioned_nodes(), - "Mesh data does not match mesh preparation().is_partitioned"); + "Mesh preparation().is_partitioned does not reflect mesh data."); // If the mesh thinks it's prepared in some way, *re*-preparing in // that way shouldn't change a clone of it, as long as we disallow // repartitioning or renumbering or remote element removal. std::unique_ptr mesh_clone = reprepared_mesh_clone(mesh); - mesh.assert_equal_to(*mesh_clone, "Mesh data does not match mesh preparation()."); + mesh.assert_equal_to + (*mesh_clone, + "Mesh data does not match mesh preparation().\n" + "Efficiently-prepared mesh != exhaustively-prepared clone."); } From 0c74e78fd3fd358ecf9fabc414b8d43ccb4c0c5f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:24:40 -0600 Subject: [PATCH 24/48] Update subdomain sets after subdomains --- examples/fem_system/fem_system_ex4/fem_system_ex4.C | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/fem_system/fem_system_ex4/fem_system_ex4.C b/examples/fem_system/fem_system_ex4/fem_system_ex4.C index 90782ac7539..c8c4716a73c 100644 --- a/examples/fem_system/fem_system_ex4/fem_system_ex4.C +++ b/examples/fem_system/fem_system_ex4/fem_system_ex4.C @@ -156,6 +156,7 @@ int main (int argc, char ** argv) for (auto & elem : mesh.element_ptr_range()) if (elem->dim() < dim) elem->subdomain_id() = 1; + mesh.cache_elem_data(); mesh_refinement.uniformly_refine(coarserefinements); From 50a8ec40203647e303d283e811ccba486b8fd8f5 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:25:37 -0600 Subject: [PATCH 25/48] Fix broken error message --- src/mesh/mesh_base.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index dea9747a920..af2792c0684 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -284,8 +284,7 @@ void MeshBase::assert_equal_to (const MeshBase & other_mesh, if (!this->processor_id()) { std::string error_msg {failure_context}; - error_msg += '\n' + - "Meshes failed asserted equality in at least these aspects:\n"; + error_msg += "\nMeshes failed asserted equality in at least these aspects:\n"; for (auto & diff : unique_diffs) { error_msg += diff; From 727d738ca78ce9ef29e6c749f83672789a2a6e68 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:26:19 -0600 Subject: [PATCH 26/48] Fix unused-parameter warnings --- src/mesh/mesh_base.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index af2792c0684..584d9c69db7 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -299,6 +299,8 @@ void MeshBase::assert_equal_to (const MeshBase & other_mesh, // know user code is going to catch that error and sync up with // us later. } +#else + libmesh_ignore(other_mesh, failure_context); #endif // NDEBUG } From 044f12260f45d818b8d19a79b3d90e04a707baa2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:26:37 -0600 Subject: [PATCH 27/48] Use libmesh_assert_valid_is_prepared() This should be much better for debugging cases where the assert fails. --- src/mesh/mesh_base.C | 4 ++-- src/mesh/mesh_tet_interface.C | 2 +- src/systems/fem_system.C | 6 +++--- src/systems/nonlinear_implicit_system.C | 2 +- src/systems/system.C | 6 +++--- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 584d9c69db7..cad02e90497 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -322,7 +322,7 @@ std::string_view MeshBase::first_difference_from(const MeshBase & other_mesh) co // // We use separate tests here and return strings for each test, // to make it easy to see the failing condition a - // MeshTools::valid_is_prepared(mesh) is failing. + // MeshTools::libmesh_valid_is_prepared(mesh) is failing. #define CHECK_MEMBER(member_name) \ if (member_name != other_mesh.member_name) \ @@ -1004,7 +1004,7 @@ void MeshBase::complete_preparation() #ifdef DEBUG // The if() here avoids both unnecessary work *and* stack overflow if (was_partly_prepared) - libmesh_assert(MeshTools::valid_is_prepared(*this)); + MeshTools::libmesh_assert_valid_is_prepared(*this); MeshTools::libmesh_assert_valid_boundary_ids(*this); #ifdef LIBMESH_ENABLE_UNIQUE_ID diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index 84d84bc724f..962fa141fda 100644 --- a/src/mesh/mesh_tet_interface.C +++ b/src/mesh/mesh_tet_interface.C @@ -132,7 +132,7 @@ BoundingBox MeshTetInterface::volume_to_surface_mesh(UnstructuredMesh & mesh) { // If we've been handed an unprepared mesh then we need to be made // aware of that and fix that; we're relying on neighbor pointers. - libmesh_assert(MeshTools::valid_is_prepared(mesh)); + MeshTools::libmesh_assert_valid_is_prepared(mesh); if (!mesh.is_prepared()) mesh.prepare_for_use(); diff --git a/src/systems/fem_system.C b/src/systems/fem_system.C index 90f4292bda9..1398ad6348e 100644 --- a/src/systems/fem_system.C +++ b/src/systems/fem_system.C @@ -901,7 +901,7 @@ void FEMSystem::assembly (bool get_residual, bool get_jacobian, libmesh_assert(mesh.is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(mesh)); + MeshTools::libmesh_assert_valid_is_prepared(mesh); #endif if (print_solution_norms) @@ -1158,7 +1158,7 @@ void FEMSystem::assemble_qoi (const QoISet & qoi_indices) libmesh_assert(mesh.is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(mesh)); + MeshTools::libmesh_assert_valid_is_prepared(mesh); #endif this->update(); @@ -1197,7 +1197,7 @@ void FEMSystem::assemble_qoi_derivative (const QoISet & qoi_indices, libmesh_assert(mesh.is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(mesh)); + MeshTools::libmesh_assert_valid_is_prepared(mesh); #endif this->update(); diff --git a/src/systems/nonlinear_implicit_system.C b/src/systems/nonlinear_implicit_system.C index a4b078d735a..2897e00d891 100644 --- a/src/systems/nonlinear_implicit_system.C +++ b/src/systems/nonlinear_implicit_system.C @@ -250,7 +250,7 @@ void NonlinearImplicitSystem::assembly(bool get_residual, { libmesh_assert(this->get_mesh().is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); + MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif // Get current_local_solution in sync diff --git a/src/systems/system.C b/src/systems/system.C index 64e81413533..ac8c5fbba54 100644 --- a/src/systems/system.C +++ b/src/systems/system.C @@ -554,7 +554,7 @@ void System::assemble () libmesh_assert(this->get_mesh().is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); + MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif // Call the user-specified assembly function @@ -570,7 +570,7 @@ void System::assemble_qoi (const QoISet & qoi_indices) libmesh_assert(this->get_mesh().is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); + MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif // Call the user-specified quantity of interest function @@ -588,7 +588,7 @@ void System::assemble_qoi_derivative(const QoISet & qoi_indices, libmesh_assert(this->get_mesh().is_prepared()); #ifdef DEBUG - libmesh_assert(MeshTools::valid_is_prepared(this->get_mesh())); + MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif // Call the user-specified quantity of interest function From 0867deb926b24dd4610a5e48f0842f51410abc1b Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 15:31:20 -0600 Subject: [PATCH 28/48] Mark inline alternatives as inline --- include/mesh/mesh_tools.h | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h index 016b94bee2a..2f9f89b1143 100644 --- a/include/mesh/mesh_tools.h +++ b/include/mesh/mesh_tools.h @@ -659,30 +659,30 @@ void globally_renumber_nodes_and_elements (MeshBase &); // Declare inline no-ops for assertions with NDEBUG defined #ifdef NDEBUG -void libmesh_assert_valid_is_prepared (const MeshBase &) {} +inline void libmesh_assert_valid_is_prepared (const MeshBase &) {} -void libmesh_assert_equal_n_systems (const MeshBase &) {} +inline void libmesh_assert_equal_n_systems (const MeshBase &) {} -void libmesh_assert_old_dof_objects (const MeshBase &) {} +inline void libmesh_assert_old_dof_objects (const MeshBase &) {} -void libmesh_assert_valid_node_pointers (const MeshBase &) {} +inline void libmesh_assert_valid_node_pointers (const MeshBase &) {} -void libmesh_assert_valid_remote_elems (const MeshBase &) {} +inline void libmesh_assert_valid_remote_elems (const MeshBase &) {} -void libmesh_assert_valid_elem_ids (const MeshBase &) {} +inline void libmesh_assert_valid_elem_ids (const MeshBase &) {} -void libmesh_assert_valid_amr_elem_ids (const MeshBase &) {} +inline void libmesh_assert_valid_amr_elem_ids (const MeshBase &) {} -void libmesh_assert_valid_amr_interior_parents (const MeshBase &) {} +inline void libmesh_assert_valid_amr_interior_parents (const MeshBase &) {} -void libmesh_assert_contiguous_dof_ids (const MeshBase &, unsigned int) {} +inline void libmesh_assert_contiguous_dof_ids (const MeshBase &, unsigned int) {} template -void libmesh_assert_topology_consistent_procids (const MeshBase &) {} +inline void libmesh_assert_topology_consistent_procids (const MeshBase &) {} -void libmesh_assert_canonical_node_procids (const MeshBase &) {} +inline void libmesh_assert_canonical_node_procids (const MeshBase &) {} -void libmesh_assert_valid_refinement_tree (const MeshBase &) {} +inline void libmesh_assert_valid_refinement_tree (const MeshBase &) {} #endif // NDEBUG From 8d2459bad6b6e7f1a5962f7322512e8ef0cdfdb4 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 18:09:34 -0600 Subject: [PATCH 29/48] Add clone() to our test GhostingFunctor subclasses This was supposed to be a pure virtual function eventually, and it's time to start moving that way. --- tests/base/overlapping_coupling_test.C | 11 +++++++++++ tests/systems/systems_test.C | 5 +++++ 2 files changed, 16 insertions(+) diff --git a/tests/base/overlapping_coupling_test.C b/tests/base/overlapping_coupling_test.C index 2b71345a077..23d9f71a718 100644 --- a/tests/base/overlapping_coupling_test.C +++ b/tests/base/overlapping_coupling_test.C @@ -60,6 +60,17 @@ class OverlappingCouplingFunctor : public GhostingFunctor virtual ~OverlappingCouplingFunctor(){}; + virtual std::unique_ptr clone () const + { + auto clone_functor = + std::make_unique(_system); + + auto coupling_matrix = + std::make_unique(*_coupling_matrix); + clone_functor->set_coupling_matrix(coupling_matrix); + return clone_functor; + } + void set_coupling_matrix (std::unique_ptr & coupling_matrix) { _coupling_matrix = std::move(coupling_matrix); } diff --git a/tests/systems/systems_test.C b/tests/systems/systems_test.C index a181cd628b6..275ce9360e2 100644 --- a/tests/systems/systems_test.C +++ b/tests/systems/systems_test.C @@ -62,6 +62,11 @@ public: _mesh(mesh) {} + virtual std::unique_ptr clone () const + { + return std::make_unique(_mesh); + } + /** * User-defined function to augment the sparsity pattern. */ From 44f1e9d4956b83cd937697d7a63565743efcafa3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:19:23 -0600 Subject: [PATCH 30/48] More updates of mesh caches --- examples/fem_system/fem_system_ex3/fem_system_ex3.C | 2 ++ examples/fem_system/fem_system_ex4/fem_system_ex4.C | 2 ++ examples/fem_system/fem_system_ex5/fem_system_ex5.C | 3 +++ examples/subdomains/subdomains_ex1/subdomains_ex1.C | 3 +++ examples/subdomains/subdomains_ex2/subdomains_ex2.C | 3 +++ .../systems_of_equations_ex6/systems_of_equations_ex6.C | 4 ++++ 6 files changed, 17 insertions(+) diff --git a/examples/fem_system/fem_system_ex3/fem_system_ex3.C b/examples/fem_system/fem_system_ex3/fem_system_ex3.C index 1ca8f370592..7afeadb1b1c 100644 --- a/examples/fem_system/fem_system_ex3/fem_system_ex3.C +++ b/examples/fem_system/fem_system_ex3/fem_system_ex3.C @@ -173,6 +173,8 @@ int main (int argc, char ** argv) mesh.get_boundary_info().add_edge(elem, e, edge_boundary_id); } + // We're all done adding to the boundary_info; let's make sure its + // caches know about it. mesh.get_boundary_info().regenerate_id_sets(); // Create an equation systems object. diff --git a/examples/fem_system/fem_system_ex4/fem_system_ex4.C b/examples/fem_system/fem_system_ex4/fem_system_ex4.C index c8c4716a73c..ba93de15483 100644 --- a/examples/fem_system/fem_system_ex4/fem_system_ex4.C +++ b/examples/fem_system/fem_system_ex4/fem_system_ex4.C @@ -156,6 +156,8 @@ int main (int argc, char ** argv) for (auto & elem : mesh.element_ptr_range()) if (elem->dim() < dim) elem->subdomain_id() = 1; + + // Make sure the mesh knows we added new subdomains. mesh.cache_elem_data(); mesh_refinement.uniformly_refine(coarserefinements); diff --git a/examples/fem_system/fem_system_ex5/fem_system_ex5.C b/examples/fem_system/fem_system_ex5/fem_system_ex5.C index fc6e33faa6b..d2e566b5a28 100644 --- a/examples/fem_system/fem_system_ex5/fem_system_ex5.C +++ b/examples/fem_system/fem_system_ex5/fem_system_ex5.C @@ -225,6 +225,9 @@ int main (int argc, char ** argv) } } + // We're done modifying the BoundaryInfo; get its caches up to date. + mesh.get_boundary_info().regenerate_id_sets(); + // Create an equation systems object. EquationSystems equation_systems (mesh); diff --git a/examples/subdomains/subdomains_ex1/subdomains_ex1.C b/examples/subdomains/subdomains_ex1/subdomains_ex1.C index 8b5e121f4ac..10e85dfe585 100644 --- a/examples/subdomains/subdomains_ex1/subdomains_ex1.C +++ b/examples/subdomains/subdomains_ex1/subdomains_ex1.C @@ -243,6 +243,9 @@ int main (int argc, char ** argv) elem->subdomain_id() = 1; } + // Make sure the mesh knows we added new subdomains. + mesh.cache_elem_data(); + // Create an equation systems object. EquationSystems equation_systems (mesh); diff --git a/examples/subdomains/subdomains_ex2/subdomains_ex2.C b/examples/subdomains/subdomains_ex2/subdomains_ex2.C index abe827d038a..47966c45126 100644 --- a/examples/subdomains/subdomains_ex2/subdomains_ex2.C +++ b/examples/subdomains/subdomains_ex2/subdomains_ex2.C @@ -195,6 +195,9 @@ int main (int argc, char ** argv) elem->subdomain_id() = 1; } + // Make sure the mesh knows we added new subdomains. + mesh.cache_elem_data(); + // Print information about the mesh to the screen. mesh.print_info(); diff --git a/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C b/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C index f0fc5c9ba80..61cc9073804 100644 --- a/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C +++ b/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C @@ -507,6 +507,10 @@ int main (int argc, char ** argv) mesh.get_boundary_info().add_edge(elem, e, EDGE_BOUNDARY_ID); } + // We're all done adding to the boundary_info; let's make sure its + // caches know about it. + mesh.get_boundary_info().regenerate_id_sets(); + // Create an equation systems object. EquationSystems equation_systems (mesh); From 88ac485bd365928c848f865326c29b42d7e20eb6 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:19:53 -0600 Subject: [PATCH 31/48] Implement clone() for misc_ex9 ghosting functor We're making that pure virtual in non-deprecated builds now. --- .../miscellaneous_ex9/augment_sparsity_on_interface.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h b/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h index 6b2e1ca5d60..08fe3d91518 100644 --- a/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h +++ b/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h @@ -56,6 +56,13 @@ class AugmentSparsityOnInterface : public GhostingFunctor boundary_id_type crack_boundary_lower, boundary_id_type crack_boundary_upper); + + virtual std::unique_ptr clone () const + { + return std::make_unique + (_mesh, _crack_boundary_lower, _crack_boundary_upper); + } + /** * @return a const reference to the lower-to-upper element ID map. */ From ab9305b85629ad1098eaf0c38739f06f528b9de1 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:20:44 -0600 Subject: [PATCH 32/48] GhostingFunctor clone is pure virtual (at least in non-deprecated builds) --- include/ghosting/ghosting_functor.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/include/ghosting/ghosting_functor.h b/include/ghosting/ghosting_functor.h index fa1fb62cee6..b9e5458d3a4 100644 --- a/include/ghosting/ghosting_functor.h +++ b/include/ghosting/ghosting_functor.h @@ -213,10 +213,12 @@ class GhostingFunctor : public ReferenceCountedObject * different meshes. The operations in GhostingFunctor are mesh dependent. */ virtual std::unique_ptr clone () const - // Let us return nullptr for backward compatibility. - // We will come back to mark this function as pure virtual - // once the API upgrade is done. - { return nullptr; } +#ifndef LIBMESH_ENABLE_DEPRECATED + = 0; +#else + // Return nullptr for backward compatibility. + { libmesh_deprecated(); return nullptr; } +#endif /** * It should be called after cloning a ghosting functor. From 63df2c253bfccc4de0a65a434b37f15343a20a25 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:24:54 -0600 Subject: [PATCH 33/48] Only enforce preparation with --disable-deprecated This is actually a big behavior change, considering it led to failed assertions in a half dozen of our own examples. Let's give users some time to get up to speed. --- src/systems/fem_system.C | 6 +++--- src/systems/nonlinear_implicit_system.C | 2 +- src/systems/system.C | 6 +++--- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/systems/fem_system.C b/src/systems/fem_system.C index 1398ad6348e..51eda44382d 100644 --- a/src/systems/fem_system.C +++ b/src/systems/fem_system.C @@ -900,7 +900,7 @@ void FEMSystem::assembly (bool get_residual, bool get_jacobian, const MeshBase & mesh = this->get_mesh(); libmesh_assert(mesh.is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(mesh); #endif @@ -1157,7 +1157,7 @@ void FEMSystem::assemble_qoi (const QoISet & qoi_indices) const MeshBase & mesh = this->get_mesh(); libmesh_assert(mesh.is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(mesh); #endif @@ -1196,7 +1196,7 @@ void FEMSystem::assemble_qoi_derivative (const QoISet & qoi_indices, const MeshBase & mesh = this->get_mesh(); libmesh_assert(mesh.is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(mesh); #endif diff --git a/src/systems/nonlinear_implicit_system.C b/src/systems/nonlinear_implicit_system.C index 2897e00d891..7f21f91acea 100644 --- a/src/systems/nonlinear_implicit_system.C +++ b/src/systems/nonlinear_implicit_system.C @@ -249,7 +249,7 @@ void NonlinearImplicitSystem::assembly(bool get_residual, bool /*apply_no_constraints*/) { libmesh_assert(this->get_mesh().is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif diff --git a/src/systems/system.C b/src/systems/system.C index ac8c5fbba54..9fced7d925a 100644 --- a/src/systems/system.C +++ b/src/systems/system.C @@ -553,7 +553,7 @@ void System::assemble () LOG_SCOPE("assemble()", "System"); libmesh_assert(this->get_mesh().is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif @@ -569,7 +569,7 @@ void System::assemble_qoi (const QoISet & qoi_indices) LOG_SCOPE("assemble_qoi()", "System"); libmesh_assert(this->get_mesh().is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif @@ -587,7 +587,7 @@ void System::assemble_qoi_derivative(const QoISet & qoi_indices, LOG_SCOPE("assemble_qoi_derivative()", "System"); libmesh_assert(this->get_mesh().is_prepared()); -#ifdef DEBUG +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); #endif From 20876b1c5fe79e76b14e6952c487d6d07f3a5b8f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 11 Dec 2025 11:49:45 -0600 Subject: [PATCH 34/48] Mark overriding functions with override --- .../miscellaneous_ex9/augment_sparsity_on_interface.h | 2 +- tests/base/overlapping_coupling_test.C | 2 +- tests/systems/systems_test.C | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h b/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h index 08fe3d91518..229c6c6efb8 100644 --- a/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h +++ b/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h @@ -57,7 +57,7 @@ class AugmentSparsityOnInterface : public GhostingFunctor boundary_id_type crack_boundary_upper); - virtual std::unique_ptr clone () const + virtual std::unique_ptr clone () const override { return std::make_unique (_mesh, _crack_boundary_lower, _crack_boundary_upper); diff --git a/tests/base/overlapping_coupling_test.C b/tests/base/overlapping_coupling_test.C index 23d9f71a718..29f29802324 100644 --- a/tests/base/overlapping_coupling_test.C +++ b/tests/base/overlapping_coupling_test.C @@ -60,7 +60,7 @@ class OverlappingCouplingFunctor : public GhostingFunctor virtual ~OverlappingCouplingFunctor(){}; - virtual std::unique_ptr clone () const + virtual std::unique_ptr clone () const override { auto clone_functor = std::make_unique(_system); diff --git a/tests/systems/systems_test.C b/tests/systems/systems_test.C index 275ce9360e2..edd31f31ecc 100644 --- a/tests/systems/systems_test.C +++ b/tests/systems/systems_test.C @@ -62,7 +62,7 @@ public: _mesh(mesh) {} - virtual std::unique_ptr clone () const + virtual std::unique_ptr clone () const override { return std::make_unique(_mesh); } From de9593ca8702d144e35cd72fddc5d607d80fc145 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 12 Dec 2025 11:14:40 -0600 Subject: [PATCH 35/48] Allow for remote_elem interior_parent() values We can have those on a ghosted boundary element --- src/mesh/unstructured_mesh.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C index ac9da962c8e..38bd469cc51 100644 --- a/src/mesh/unstructured_mesh.C +++ b/src/mesh/unstructured_mesh.C @@ -836,7 +836,8 @@ void UnstructuredMesh::copy_nodes_and_elements(const MeshBase & other_mesh, for (auto & elem_pair : ip_map) { Elem * ip = const_cast(elem_pair.first->interior_parent()); - libmesh_assert(ip == other_interior_mesh->elem_ptr(ip->id())); + libmesh_assert(ip == remote_elem || + ip == other_interior_mesh->elem_ptr(ip->id())); elem_pair.second->set_interior_parent(ip); } else From 039fcfc98ae7ebbeb9ce6a8364c6470352d20213 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 12 Dec 2025 11:26:49 -0600 Subject: [PATCH 36/48] Fix InfQuad::is_child_on_side() 2008. This has been buggy since I wrote it in 2008. And yet it took a double-refinement test on a distributed mesh plus an extra layer of dbg-mode internal consistency checks to catch the bug. --- src/geom/face_inf_quad.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/geom/face_inf_quad.C b/src/geom/face_inf_quad.C index 6070eed4fa9..6abfb1f21a6 100644 --- a/src/geom/face_inf_quad.C +++ b/src/geom/face_inf_quad.C @@ -184,7 +184,9 @@ bool InfQuad::is_child_on_side(const unsigned int c, libmesh_assert_less (c, this->n_children()); libmesh_assert_less (s, this->n_sides()); - return (s == 0 || s == c+1); + return (s == 0 || + (c == 0 && s == 2) || + (c == 1 && s == 1)); } From 10979c6e2a85fe1c198456ca536ca386cb2eb77e Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 12 Dec 2025 11:28:19 -0600 Subject: [PATCH 37/48] More InfElem test coverage for is_child_on_* --- tests/geom/elem_test.C | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C index 578ff186835..23d23480f64 100644 --- a/tests/geom/elem_test.C +++ b/tests/geom/elem_test.C @@ -679,6 +679,21 @@ public: for (const Node & node : child_side->node_ref_range()) CPPUNIT_ASSERT(parent_side->contains_point(node)); + // We can at least check for sharing some node with + // a side, on both finite and infinite elements, for + // those children that aren't just the middle elements + // of triangular sides. + bool shares_a_side_node = false; + bool shares_a_vertex = false; + for (const Node & node : child_side->node_ref_range()) + if (parent_side->local_node(node.id()) != invalid_uint) + shares_a_side_node = true; + for (const Node & node : elem->node_ref_range()) + if (parent->local_node(node.id()) < parent->n_vertices()) + shares_a_vertex = true; + + CPPUNIT_ASSERT(shares_a_side_node || !shares_a_vertex); + if (elem->neighbor_ptr(s) && !elem->neighbor_ptr(s)->is_remote()) CPPUNIT_ASSERT_EQUAL(parent->child_neighbor(elem->neighbor_ptr(s)), elem); } @@ -699,6 +714,15 @@ public: if (!parent_edge->infinite()) for (const Node & node : child_edge->node_ref_range()) CPPUNIT_ASSERT(parent_edge->contains_point(node)); + + // We can at least check for sharing some node with + // an edge, on both finite and infinite elements + bool shares_an_edge_node = false; + for (const Node & node : child_edge->node_ref_range()) + if (parent_edge->local_node(node.id()) != invalid_uint) + shares_an_edge_node = true; + + CPPUNIT_ASSERT(shares_an_edge_node); } } From 77156027148353025aa94e729d42b388a85d9dd9 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 15 Dec 2025 12:50:02 -0600 Subject: [PATCH 38/48] Handle skip_partitioning() in assertions If we're skipping partitioning, we shouldn't mark ourselves as partitioned if we're not, and we shouldn't expect that we're marked as partitioned. --- src/mesh/mesh_base.C | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index cad02e90497..0e35a636151 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -976,7 +976,8 @@ void MeshBase::complete_preparation() // partition() call will still check for orphaned nodes. if (!skip_partitioning() && !_preparation.is_partitioned) this->partition(); - else + else if (!this->n_unpartitioned_elem() && + !this->n_unpartitioned_nodes()) _preparation.is_partitioned = true; // If we're using DistributedMesh, we'll probably want it @@ -998,8 +999,15 @@ void MeshBase::complete_preparation() if (!_skip_renumber_nodes_and_elements) this->renumber_nodes_and_elements(); - // The mesh is now prepared for use, and it should know it. - libmesh_assert(_preparation); + // The mesh is now prepared for use, with the possible exception of + // partitioning that was supposed to be skipped, and it should know + // it. +#ifndef NDEBUG + Preparation completed_preparation = _preparation; + if (skip_partitioning()) + completed_preparation.is_partitioned = true; + libmesh_assert(completed_preparation); +#endif #ifdef DEBUG // The if() here avoids both unnecessary work *and* stack overflow From 96dba18dbdb3b857595b5594d266d70ac0e04ff6 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 18 Dec 2025 14:25:09 -0600 Subject: [PATCH 39/48] Preparation::libmesh_assert_consistent() I haven't actually seen this fail yet but it does seem like the sort of thing we should be checking. --- include/mesh/mesh_base.h | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h index 883ff282348..7c753401fd9 100644 --- a/include/mesh/mesh_base.h +++ b/include/mesh/mesh_base.h @@ -2046,6 +2046,18 @@ class MeshBase : public ParallelObject has_reinit_ghosting_functors == other.has_reinit_ghosting_functors && has_boundary_id_sets == other.has_boundary_id_sets; } + + void libmesh_assert_consistent (const Parallel::Communicator & libmesh_dbg_var(comm)) { + libmesh_assert(comm.verify(is_partitioned)); + libmesh_assert(comm.verify(has_synched_id_counts)); + libmesh_assert(comm.verify(has_neighbor_ptrs)); + libmesh_assert(comm.verify(has_cached_elem_data)); + libmesh_assert(comm.verify(has_interior_parent_ptrs)); + libmesh_assert(comm.verify(has_removed_remote_elements)); + libmesh_assert(comm.verify(has_removed_orphaned_nodes)); + libmesh_assert(comm.verify(has_reinit_ghosting_functors)); + libmesh_assert(comm.verify(has_boundary_id_sets)); + } }; protected: From 858a657e3d636d1cb0857607a854304ce4d4f921 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 18 Dec 2025 14:25:40 -0600 Subject: [PATCH 40/48] Fix libmesh_dbg_var() comment I went to double-check that this was the macro I wanted, and the behavior did match my recollection but didn't match my comment --- include/base/libmesh_common.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/base/libmesh_common.h b/include/base/libmesh_common.h index 7ee03d06a68..f0a93a67a98 100644 --- a/include/base/libmesh_common.h +++ b/include/base/libmesh_common.h @@ -271,7 +271,7 @@ extern bool warned_about_auto_ptr; } while (0) // The libmesh_dbg_var() macro indicates that an argument to a function -// is used only in debug mode (i.e., when NDEBUG is not defined). +// is used only in debug and devel modes (i.e., when NDEBUG is not defined). #ifndef NDEBUG #define libmesh_dbg_var(var) var #else From d79bf9107802736dcef75973962bbc5e42ccf296 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 18 Dec 2025 14:26:34 -0600 Subject: [PATCH 41/48] Fix MeshBase::assert_equal_to error handling --- src/mesh/mesh_base.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 0e35a636151..a9cb7b70206 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -274,6 +274,8 @@ void MeshBase::assert_equal_to (const MeshBase & other_mesh, std::string_view local_diff = first_difference_from(other_mesh); bool diff_found = !local_diff.empty(); + this->comm().max(diff_found); + if (diff_found) { // Construct a user-friendly message to throw on pid 0 From a32107e6e0f79f1496fea5af21af9b5abf28325a Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 18 Dec 2025 14:26:52 -0600 Subject: [PATCH 42/48] Clean up MeshBase::assert_equal_to error message --- src/mesh/mesh_base.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index a9cb7b70206..39ca84ef8e9 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -280,7 +280,8 @@ void MeshBase::assert_equal_to (const MeshBase & other_mesh, { // Construct a user-friendly message to throw on pid 0 std::set unique_diffs; - unique_diffs.insert(std::string(local_diff)); + if (!local_diff.empty()) + unique_diffs.insert(std::string(local_diff)); this->comm().set_union(unique_diffs); if (!this->processor_id()) From d8a4658770da808b27344ad7d57fa262f0cba302 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 18 Dec 2025 14:27:12 -0600 Subject: [PATCH 43/48] Test for parallel-consistent preparation --- src/mesh/mesh_base.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 39ca84ef8e9..41e03e67414 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -890,6 +890,8 @@ void MeshBase::complete_preparation() libmesh_assert(this->comm().verify(this->is_serial())); + _preparation.libmesh_assert_consistent(this->comm()); + #ifdef DEBUG // If we don't go into this method with valid constraint rows, we're // only going to be able to make that worse. From 7e131f11a8cb3e44f73f0a5790c647f297cbbfb6 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:04:27 -0600 Subject: [PATCH 44/48] 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 833645d96d5..d6dc64a703b 100644 --- a/include/geom/elem.h +++ b/include/geom/elem.h @@ -1585,6 +1585,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 16144ac5a40..42138ce6870 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 03ce4cff6fa..755dfffd138 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 9cf5648ebc319b91afe6076df696e7ea29a64af1 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 15 Jan 2026 13:27:23 -0600 Subject: [PATCH 45/48] 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 d6dc64a703b..0eb92f8202d 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. @@ -2594,6 +2600,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 5ef677ce84908dc034a6277faa2c673422728f3c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:06:29 -0600 Subject: [PATCH 46/48] 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 755dfffd138..1d84e6c0053 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 4bbac14d3743cfdf6d80e51e09322568154e4fdb Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:06:52 -0600 Subject: [PATCH 47/48] 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 c81cdcb58ed..971b39b747d 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 9e6f0dbd26eae7902aba29f99f381884e6156f4c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 13 Jan 2026 17:07:52 -0600 Subject: [PATCH 48/48] 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 dbad0252965..519eda2db3a 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);