From aa69f30fd62929d574225944f46e5b40aebfadc9 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Jan 2024 12:59:24 -0600 Subject: [PATCH 01/38] 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 ca5010032d75f7f752239054d817e9035983c326 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 22 Jan 2024 10:44:24 -0600 Subject: [PATCH 02/38] 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 b2b95e0cb3482ac2eeed90b87e57208537f044fa Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 1 Dec 2025 11:51:52 -0600 Subject: [PATCH 03/38] 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 9b3e41c261bae5c0e841a5fddc311819d3eb3f44 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 1 Dec 2025 14:40:57 -0600 Subject: [PATCH 04/38] 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 95432e8d7bc1e51270f1f77cf4136fa2493108ea Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 2 Dec 2025 22:11:11 -0600 Subject: [PATCH 05/38] 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 de39e08ef4f84ddb93e864ff29f2029be4444d3e Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:13:42 -0600 Subject: [PATCH 06/38] 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 1814a3573517bc170936fc4a13cdc9a40f8db2fc Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:14:07 -0600 Subject: [PATCH 07/38] 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 64b5b0f8349ad5a35f656ff31659323e421a0a92 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:15:17 -0600 Subject: [PATCH 08/38] 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 cc28c2228718d6953296c788e433c342cbc6b1ef Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:16:33 -0600 Subject: [PATCH 09/38] 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 f29a4622bc27e761713f3b572a0c160bb35dcb05 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:17:40 -0600 Subject: [PATCH 10/38] 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 a4fa81cc25808308ef983e3a8e2d9f12b0c52af8 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 3 Dec 2025 13:32:10 -0600 Subject: [PATCH 11/38] 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 d19bcb2e3e925255ff77becbf115b0bd9a5615f3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 13:26:54 -0600 Subject: [PATCH 12/38] 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 97725fee1bdf825af41cdf8fd2ebba99a76a6424 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 13:27:21 -0600 Subject: [PATCH 13/38] 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 aa107d4971f1423b394c9ea4d77d8594af9676c8 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 20:55:01 -0600 Subject: [PATCH 14/38] 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 e39fc7537685bb9896e8c8cc6e92478dab8b4588 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 4 Dec 2025 20:56:43 -0600 Subject: [PATCH 15/38] 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 8000dc5a961dd4e54d7777865f25d22ec7626ec3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 5 Dec 2025 17:23:17 -0600 Subject: [PATCH 16/38] 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 8b0d7cb3305776e2f31204455d148fc7f8958257 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 8 Dec 2025 17:55:02 -0600 Subject: [PATCH 17/38] 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 8e20d3b4dbd73cf6b77f851259a8c0d416c5bc02 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 13:38:03 -0600 Subject: [PATCH 18/38] 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 2e00e2451457cd642648e93566c16c0b6f737163 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 13:38:31 -0600 Subject: [PATCH 19/38] 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 a8a3184acea49c261286d31f47e9ec321d099e04 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 21:46:30 -0600 Subject: [PATCH 20/38] 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 a1905040a88e3fda3e6e99c36492376ae79f21ba Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 21:50:23 -0600 Subject: [PATCH 21/38] 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 e72b98b3131aff27e526f5a2b8f2eb456a221f41 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 9 Dec 2025 21:53:36 -0600 Subject: [PATCH 22/38] 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 b3d674a455492d39491219935788a7ef4316e4ec Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:16:47 -0600 Subject: [PATCH 23/38] 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 6da6844cfcdae28ad47fc2c379c32080669ae12a Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:24:40 -0600 Subject: [PATCH 24/38] 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 8c4020447588b8a912005ad4a2c262352c4701cd Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:25:37 -0600 Subject: [PATCH 25/38] 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 223ee3c2c60d57d3ba95c0446bb04c4a000dbded Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:26:19 -0600 Subject: [PATCH 26/38] 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 ad56caae700f65e6b7a1a3d322729b12751c4645 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 13:26:37 -0600 Subject: [PATCH 27/38] 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 d75d3a5aed2a032a7d63e9db82652c0933736dd3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 15:31:20 -0600 Subject: [PATCH 28/38] 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 91559ea356b7278a8b1a883c62e38208926aff36 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 18:09:34 -0600 Subject: [PATCH 29/38] 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 ab86a36caebb763da5f2f362ddb61d326efb0d79 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:19:23 -0600 Subject: [PATCH 30/38] 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 eaa0c2c9a94c4768aa89e8aec047856ac39936f3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:19:53 -0600 Subject: [PATCH 31/38] 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 8ed38bb02dfccbc8ae08f8f2315d825f0ad4e94c Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:20:44 -0600 Subject: [PATCH 32/38] 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 5d237d7cbcab3dba6c520b388c07346d69f2d8af Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 10 Dec 2025 23:24:54 -0600 Subject: [PATCH 33/38] 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 3bb4e11f4311de40db85646d906ba03317d77849 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Thu, 11 Dec 2025 11:49:45 -0600 Subject: [PATCH 34/38] 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 15f48ed4af2669d121a22bf528b03518b06ff9b3 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 12 Dec 2025 11:14:40 -0600 Subject: [PATCH 35/38] 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 bf8fdf3d4a1e78029605fbe003cb49d87bb7bb8f Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 12 Dec 2025 11:26:49 -0600 Subject: [PATCH 36/38] 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 d130b95183949164e70917bfa162a07856f7d6d0 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Fri, 12 Dec 2025 11:28:19 -0600 Subject: [PATCH 37/38] 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 624c8d772873c2ed3561b0f3f91bc46b1012b0ed Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 15 Dec 2025 12:50:02 -0600 Subject: [PATCH 38/38] 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