diff --git a/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C b/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C index 558f91947c..c290431cf3 100644 --- a/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C +++ b/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C @@ -221,6 +221,10 @@ int main (int argc, char** argv) sync_obj); } + // We'll need to update caches of subdomain ids too, to keep the + // mesh prepared. + mesh.cache_elem_data(); + // Now we set the sources of the field: prism-shaped objects that are // determined here by containing certain points: auto p_pt_lctr = mesh.sub_point_locator(); diff --git a/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C b/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C index b8b4b07dce..e5c825c14f 100644 --- a/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C +++ b/examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C @@ -150,6 +150,9 @@ int main (int argc, char ** argv) const Real phi = infile("phi", 0.), theta = infile("theta", 0.), psi = infile("psi", 0.); GradDivExactSolution::RM(MeshTools::Modification::rotate(mesh, phi, theta, psi)); + // Mesh rotation invalidates some caches + mesh.complete_preparation(); + // Print information about the mesh to the screen. mesh.print_info(); diff --git a/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C b/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C index c12beefbb9..fc9c1e6120 100644 --- a/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C +++ b/examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C @@ -113,6 +113,9 @@ int main (int argc, char ** argv) const Real phi = infile("phi", 0.), theta = infile("theta", 0.), psi = infile("psi", 0.); CurlCurlExactSolution::RM(MeshTools::Modification::rotate(mesh, phi, theta, psi)); + // Mesh rotation invalidates some caches + mesh.complete_preparation(); + // Print information about the mesh to the screen. mesh.print_info(); diff --git a/include/mesh/distributed_mesh.h b/include/mesh/distributed_mesh.h index ec1febb473..cec39aeb4b 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 59d727a428..45ded45b42 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 @@ -239,8 +247,7 @@ class MeshBase : public ParallelObject * consider ourself prepared. This is a very coarse setting; it is * generally more efficient to mark finer-grained settings instead. */ - void unset_is_prepared() - { _preparation = false; } + void unset_is_prepared(); /** * Tells this we have done some operation creating unpartitioned @@ -423,23 +430,20 @@ class MeshBase : public ParallelObject * \returns A const reference to a std::set of element dimensions * present in the mesh. */ - const std::set & elem_dimensions() const - { return _elem_dims; } + const std::set & elem_dimensions() const; /** * \returns A const reference to a std::set of element default * orders present in the mesh. */ - const std::set & elem_default_orders() const - { return _elem_default_orders; } + const std::set & elem_default_orders() const; /** * \returns The smallest supported_nodal_order() of any element * present in the mesh, which is thus the maximum supported nodal * order on the mesh as a whole. */ - Order supported_nodal_order() const - { return _supported_nodal_order; } + Order supported_nodal_order() const; /** * Most of the time you should not need to call this, as the element @@ -982,10 +986,17 @@ class MeshBase : public ParallelObject * * If \p assert_valid is left as true, then in dbg mode extensive * consistency checking is performed before returning. + * + * If \p check_non_remote is set to false, then only sides which + * currently have remote neighbors are checked for possible local + * neighbors. This is intended to handle a corner case where + * ancestor neighbors are redistributed to a processor only by other + * processors who do not see that neighbor link. */ virtual void find_neighbors (const bool reset_remote_elements = false, const bool reset_current_list = true, - const bool assert_valid = true) = 0; + const bool assert_valid = true, + const bool check_non_remote = true) = 0; /** * Removes any orphaned nodes, nodes not connected to any elements. @@ -1997,8 +2008,7 @@ class MeshBase : public ParallelObject * should contain all the subdomain ids across processors. Relies on the mesh * being prepared */ - const std::set & get_mesh_subdomains() const - { libmesh_assert(this->is_prepared()); return _mesh_subdomains; } + const std::set & get_mesh_subdomains() const; #ifdef LIBMESH_ENABLE_PERIODIC /** @@ -2043,6 +2053,9 @@ class MeshBase : public ParallelObject bool operator== (const Preparation & other) const; bool operator!= (const Preparation & other) const; + // Assert that a Preparation object is identical across processors + void libmesh_assert_consistent (const Parallel::Communicator & libmesh_dbg_var(comm)); + bool is_partitioned; bool has_synched_id_counts; bool has_neighbor_ptrs; @@ -2092,7 +2105,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 @@ -2100,6 +2113,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. */ @@ -2546,6 +2565,48 @@ MeshBase::const_node_iterator : MeshBase::const_node_filter_iter }; +// ------------------------------------------------------------ +// Elem class member functions +inline +const std::set & MeshBase::elem_dimensions() const +{ + libmesh_assert(_preparation.has_cached_elem_data); + return _elem_dims; +} + + +inline +const std::set & MeshBase::elem_default_orders() const +{ + libmesh_assert(_preparation.has_cached_elem_data); + return _elem_default_orders; +} + + +inline +Order MeshBase::supported_nodal_order() const +{ + libmesh_assert(_preparation.has_cached_elem_data); + return _supported_nodal_order; +} + + +inline +const std::set & MeshBase::get_mesh_subdomains() const +{ + libmesh_assert(_preparation.has_cached_elem_data); + return _mesh_subdomains; +} + + +inline +unsigned int MeshBase::spatial_dimension () const +{ + libmesh_assert(_preparation.has_cached_elem_data); + + return cast_int(_spatial_dimension); +} + template inline unsigned int MeshBase::add_elem_datum(const std::string & name, diff --git a/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h index 50b4c6dd07..c1a6d3002c 100644 --- a/include/mesh/mesh_tools.h +++ b/include/mesh/mesh_tools.h @@ -394,24 +394,35 @@ 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); ///@{ /** - * 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 */ #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 @@ -644,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 + +inline void libmesh_assert_valid_is_prepared (const MeshBase &) {} + +inline void libmesh_assert_equal_n_systems (const MeshBase &) {} + +inline void libmesh_assert_old_dof_objects (const MeshBase &) {} + +inline void libmesh_assert_valid_node_pointers (const MeshBase &) {} + +inline void libmesh_assert_valid_remote_elems (const MeshBase &) {} + +inline void libmesh_assert_valid_elem_ids (const MeshBase &) {} + +inline void libmesh_assert_valid_amr_elem_ids (const MeshBase &) {} + +inline void libmesh_assert_valid_amr_interior_parents (const MeshBase &) {} + +inline void libmesh_assert_contiguous_dof_ids (const MeshBase &, unsigned int) {} + +template +inline void libmesh_assert_topology_consistent_procids (const MeshBase &) {} + +inline void libmesh_assert_canonical_node_procids (const MeshBase &) {} + +inline void libmesh_assert_valid_refinement_tree (const MeshBase &) {} + +#endif // NDEBUG + + + } // end namespace MeshTools } // namespace libMesh diff --git a/include/mesh/replicated_mesh.h b/include/mesh/replicated_mesh.h index 6cb5e60a96..fbd8d4d44b 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/include/mesh/unstructured_mesh.h b/include/mesh/unstructured_mesh.h index 8d838a67e1..c931534192 100644 --- a/include/mesh/unstructured_mesh.h +++ b/include/mesh/unstructured_mesh.h @@ -295,7 +295,8 @@ class UnstructuredMesh : public MeshBase */ virtual void find_neighbors (const bool reset_remote_elements = false, const bool reset_current_list = true, - const bool assert_valid = true) override; + const bool assert_valid = true, + const bool check_non_remote = true) override; #ifdef LIBMESH_ENABLE_AMR /** diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C index 085a1f46bc..9b6b09a2f8 100644 --- a/src/mesh/boundary_info.C +++ b/src/mesh/boundary_info.C @@ -2144,6 +2144,8 @@ void BoundaryInfo::renumber_id (boundary_id_type old_id, { _boundary_ids.erase(old_id); _boundary_ids.insert(new_id); + _global_boundary_ids.erase(old_id); + _global_boundary_ids.insert(new_id); } renumber_name(_ss_id_to_name, old_id, new_id); @@ -2183,8 +2185,10 @@ void BoundaryInfo::renumber_side_id (boundary_id_type old_id, !_node_boundary_ids.count(old_id)) { _boundary_ids.erase(old_id); + _global_boundary_ids.erase(old_id); } _boundary_ids.insert(new_id); + _global_boundary_ids.insert(new_id); } renumber_name(_ss_id_to_name, old_id, new_id); @@ -2222,8 +2226,10 @@ void BoundaryInfo::renumber_edge_id (boundary_id_type old_id, !_node_boundary_ids.count(old_id)) { _boundary_ids.erase(old_id); + _global_boundary_ids.erase(old_id); } _boundary_ids.insert(new_id); + _global_boundary_ids.insert(new_id); } renumber_name(_es_id_to_name, old_id, new_id); @@ -2261,8 +2267,10 @@ void BoundaryInfo::renumber_shellface_id (boundary_id_type old_id, !_node_boundary_ids.count(old_id)) { _boundary_ids.erase(old_id); + _global_boundary_ids.erase(old_id); } _boundary_ids.insert(new_id); + _global_boundary_ids.insert(new_id); } this->libmesh_assert_valid_multimaps(); @@ -2298,8 +2306,10 @@ void BoundaryInfo::renumber_node_id (boundary_id_type old_id, !_edge_boundary_ids.count(old_id)) { _boundary_ids.erase(old_id); + _global_boundary_ids.erase(old_id); } _boundary_ids.insert(new_id); + _global_boundary_ids.insert(new_id); } renumber_name(_ns_id_to_name, old_id, new_id); diff --git a/src/mesh/distributed_mesh.C b/src/mesh/distributed_mesh.C index f78ec505bc..eb03428a2c 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 () @@ -1038,9 +1042,13 @@ void DistributedMesh::redistribute () this->update_parallel_id_counts(); - // We ought to still have valid neighbor links; we communicate - // them for newly-redistributed elements - // this->find_neighbors(); + // We communicate valid neighbor links for newly-redistributed + // elements, but we may still have remote_elem links that should + // be corrected but whose correction wasn't communicated. + this->find_neighbors(/*reset_remote_elements*/ false, + /*reset_current_list*/ false /*non-default*/, + /*assert_valid*/ false, + /*check_non_remote*/ false /*non-default!*/); // Is this necessary? If we are called from prepare_for_use(), this will be called // anyway... but users can always call partition directly, in which case we do need diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C index 7dc7a6f31c..04812b636d 100644 --- a/src/mesh/mesh_base.C +++ b/src/mesh/mesh_base.C @@ -270,7 +270,57 @@ 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("assert_equal_to()", "MeshBase"); + + std::string_view local_diff = first_difference_from(other_mesh); + + bool diff_found = !local_diff.empty(); + this->comm().max(diff_found); + + if (diff_found) + { + // Construct a user-friendly message to throw on pid 0 + std::set unique_diffs; + if (!local_diff.empty()) + 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 += "\nMeshes 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. + } +#else + libmesh_ignore(other_mesh, failure_context); +#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 // @@ -278,21 +328,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::libmesh_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 @@ -300,61 +348,43 @@ 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; + 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(_skip_detect_interior_parents); + 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 (_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 (_skip_detect_interior_parents != other_mesh._skip_detect_interior_parents) - 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; 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(); @@ -363,15 +393,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)) { @@ -384,14 +414,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* @@ -400,13 +430,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); } @@ -590,11 +620,6 @@ void MeshBase::change_elemset_id(elemset_id_type old_id, elemset_id_type new_id) _all_elemset_ids.insert(new_id); } -unsigned int MeshBase::spatial_dimension () const -{ - return cast_int(_spatial_dimension); -} - void MeshBase::set_spatial_dimension(unsigned char d) @@ -865,10 +890,16 @@ void MeshBase::complete_preparation() libmesh_assert(this->comm().verify(this->is_serial())); + _preparation.libmesh_assert_consistent(this->comm()); + +#ifdef DEBUG // If we don't go into this method with valid constraint rows, we're // only going to be able to make that worse. -#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 @@ -993,6 +1024,10 @@ void MeshBase::complete_preparation() #endif #ifdef DEBUG + // The if() here avoids both unnecessary work *and* stack overflow + if (was_partly_prepared) + MeshTools::libmesh_assert_valid_is_prepared(*this); + MeshTools::libmesh_assert_valid_boundary_ids(*this); #ifdef LIBMESH_ENABLE_UNIQUE_ID MeshTools::libmesh_assert_valid_unique_ids(*this); @@ -1044,6 +1079,13 @@ bool MeshBase::is_prepared() const return static_cast(_preparation); } +void MeshBase::unset_is_prepared() +{ + _preparation = false; +} + + + void MeshBase::add_ghosting_functor(GhostingFunctor & ghosting_functor) { // We used to implicitly support duplicate inserts to std::set @@ -1228,7 +1270,10 @@ std::string MeshBase::get_info(const unsigned int verbosity /* = 0 */, const boo --_elem_dims.end(), // --end() is valid if the set is non-empty std::ostream_iterator(oss, ", ")); oss << cast_int(*_elem_dims.rbegin()); - oss << "}\n"; + oss << "}"; + if (!this->preparation().has_cached_elem_data) + oss << " (may be out of date)"; + oss << '\n'; } if (!_elem_default_orders.empty()) @@ -1240,12 +1285,23 @@ std::string MeshBase::get_info(const unsigned int verbosity /* = 0 */, const boo [](Order o) { return Utility::enum_to_string(o); }); oss << Utility::enum_to_string(*_elem_default_orders.rbegin()); - oss << "}\n"; + oss << "}"; + if (!this->preparation().has_cached_elem_data) + oss << " (may be out of date)"; + oss << '\n'; } - oss << " supported_nodal_order()=" << this->supported_nodal_order() << '\n' - << " spatial_dimension()=" << this->spatial_dimension() << '\n' - << " n_nodes()=" << this->n_nodes() << '\n' + oss << " supported_nodal_order()=" << this->_supported_nodal_order; + if (!this->preparation().has_cached_elem_data) + oss << " (may be out of date)"; + oss << '\n'; + + oss << " spatial_dimension()=" << int(this->_spatial_dimension); + if (!this->preparation().has_cached_elem_data) + oss << " (may be out of date)"; + oss << '\n'; + + oss << " n_nodes()=" << this->n_nodes() << '\n' << " n_local_nodes()=" << this->n_local_nodes() << '\n' << " n_elem()=" << this->n_elem() << '\n' << " n_local_elem()=" << this->n_local_elem() << '\n'; @@ -1965,19 +2021,31 @@ 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); + // This requires up-to-date mesh dimensions, but if we don't have + // them cached then we can't update them without changing the mesh + // in unexpected ways that interfere with our tests of + // partially-prepared meshes in MeshTools::*valid_is_prepared + std::set elem_dims_copy; + if (_preparation.has_cached_elem_data) + elem_dims_copy = this->elem_dimensions(); + else + { + for (const auto & elem : this->active_element_ptr_range()) + elem_dims_copy.insert(cast_int(elem->dim())); + if (!this->is_serial()) + this->comm().set_union(elem_dims_copy); + } // Early return if the mesh is empty or has elements of a single spatial dimension. - if (this->elem_dimensions().size() <= 1) + if (elem_dims_copy.size() <= 1) { _preparation.has_interior_parent_ptrs = true; return; } // Convenient elem_dimensions iterators - const auto dim_start = this->elem_dimensions().begin(); - const auto dim_end = this->elem_dimensions().end(); + const auto dim_start = elem_dims_copy.begin(); + const auto dim_end = elem_dims_copy.end(); // In this function we find only +1 dimensional interior parents, // (so, for a given element el, the interior parent p must satisfy p.dim() == el.dim() + 1). @@ -2779,6 +2847,21 @@ MeshBase::Preparation::operator!= (const Preparation & other) const } +void +MeshBase::Preparation::libmesh_assert_consistent (const Parallel::Communicator & libmesh_dbg_var(comm)) +{ + libmesh_assert(comm.verify(is_partitioned)); + libmesh_assert(comm.verify(has_synched_id_counts)); + libmesh_assert(comm.verify(has_neighbor_ptrs)); + libmesh_assert(comm.verify(has_cached_elem_data)); + libmesh_assert(comm.verify(has_interior_parent_ptrs)); + libmesh_assert(comm.verify(has_removed_remote_elements)); + libmesh_assert(comm.verify(has_removed_orphaned_nodes)); + libmesh_assert(comm.verify(has_reinit_ghosting_functors)); + libmesh_assert(comm.verify(has_boundary_id_sets)); +} + + // Explicit instantiations for our template function template LIBMESH_EXPORT void MeshBase::copy_constraint_rows(const SparseMatrix & constraint_operator, diff --git a/src/mesh/mesh_modification.C b/src/mesh/mesh_modification.C index 5c27728245..2f11e52162 100644 --- a/src/mesh/mesh_modification.C +++ b/src/mesh/mesh_modification.C @@ -307,6 +307,10 @@ void MeshTools::Modification::redistribute (MeshBase & mesh, #endif } + // If we just moved a mesh in or out out of the X axis or XY plane + // then we might have changed its spatial_dimension() + mesh.unset_has_cached_elem_data(); + // We haven't changed any topology, but just changing geometry could // have invalidated a point locator. mesh.clear_point_locator(); @@ -324,6 +328,10 @@ void MeshTools::Modification::translate (MeshBase & mesh, for (auto & node : mesh.node_ptr_range()) *node += p; + // If we just moved a mesh in or out out of the X axis or XY plane + // then we might have changed its spatial_dimension() + mesh.unset_has_cached_elem_data(); + // We haven't changed any topology, but just changing geometry could // have invalidated a point locator. mesh.clear_point_locator(); @@ -365,15 +373,16 @@ MeshTools::Modification::rotate (MeshBase & mesh, #if LIBMESH_DIM == 3 const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi); - if (theta) - mesh.set_spatial_dimension(3); - for (auto & node : mesh.node_ptr_range()) { Point & pt = *node; pt = R * pt; } + // If we just moved a mesh in or out out of the X axis or XY plane + // then we might have changed its spatial_dimension() + mesh.unset_has_cached_elem_data(); + return R; #else @@ -419,6 +428,10 @@ void MeshTools::Modification::scale (MeshBase & mesh, for (auto & node : mesh.node_ptr_range()) (*node)(2) *= z_scale; + // If we just collapsed a manifold onto the X axis or XY plane + // then we might have changed its spatial_dimension() + mesh.unset_has_cached_elem_data(); + // We haven't changed any topology, but just changing geometry could // have invalidated a point locator. mesh.clear_point_locator(); diff --git a/src/mesh/mesh_smoother_vsmoother.C b/src/mesh/mesh_smoother_vsmoother.C index 3d1c2ce589..caef0ecfdf 100644 --- a/src/mesh/mesh_smoother_vsmoother.C +++ b/src/mesh/mesh_smoother_vsmoother.C @@ -63,6 +63,9 @@ VariationalMeshSmoother::VariationalMeshSmoother( void VariationalMeshSmoother::setup() { // Check for multiple dimensions + if (!_mesh.preparation().has_cached_elem_data) + _mesh.cache_elem_data(); + if (_mesh.elem_dimensions().size() > 1) libmesh_not_implemented_msg("Meshes containing elements of differing dimension are not yet supported."); diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index a99396654e..57ac0caf34 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/mesh/mesh_tools.C b/src/mesh/mesh_tools.C index cc239f759d..1d9450edf8 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,70 @@ 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 (mesh.allow_detect_interior_parents() && + 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(); + } + + // Restore original flag values + 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; +} + + } @@ -1229,37 +1293,66 @@ 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; - std::unique_ptr mesh_clone = mesh.clone(); + // 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; - // 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); + // 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); - bool old_skip_partitioning = mesh_clone->skip_partitioning(); - mesh_clone->skip_partitioning(true); + // Check whether the original and clone compare equal + return (mesh == *mesh_clone); +} - bool old_allow_remote_element_removal = mesh_clone->allow_remote_element_removal(); - mesh_clone->allow_remote_element_removal(false); - // Call prepare_for_use() on the clone - mesh_clone->prepare_for_use(); - // Restore original flag values - mesh_clone->allow_renumbering(old_allow_renumbering); - mesh_clone->skip_partitioning(old_skip_partitioning); - mesh_clone->allow_remote_element_removal(old_allow_remote_element_removal); +#ifndef NDEBUG - // Check whether the original and clone compare equal - return (mesh == *mesh_clone); + +void libmesh_assert_valid_is_prepared (const MeshBase & mesh) +{ + LOG_SCOPE("libmesh_assert_valid_is_prepared()", "MeshTools"); + + 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; + + // 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 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().\n" + "Efficiently-prepared mesh != exhaustively-prepared clone."); } -#ifndef NDEBUG + void libmesh_assert_equal_n_systems (const MeshBase & mesh) { diff --git a/src/mesh/replicated_mesh.C b/src/mesh/replicated_mesh.C index 7310c287ef..9641e75232 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 ""; } diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C index 700ea19444..251b922e5f 100644 --- a/src/mesh/unstructured_mesh.C +++ b/src/mesh/unstructured_mesh.C @@ -900,7 +900,7 @@ void UnstructuredMesh::copy_nodes_and_elements(const MeshBase & other_mesh, const bool skipped_partitioning = this->skip_partitioning(); this->skip_partitioning(true); - const bool was_prepared = this->is_prepared(); + const Preparation old_preparation = this->preparation(); this->prepare_for_use(); //But in the long term, don't change our policies. @@ -913,9 +913,19 @@ void UnstructuredMesh::copy_nodes_and_elements(const MeshBase & other_mesh, // That prepare_for_use() call marked us as prepared, but we // specifically avoided some important preparation, so we might not // actually be prepared now. - if (skip_find_neighbors || - !was_prepared || !other_mesh.is_prepared()) - this->unset_is_prepared(); + if (skip_find_neighbors) + this->unset_has_neighbor_ptrs(); + + const Preparation other_preparation = other_mesh.preparation(); + if (!old_preparation.is_partitioned || + !other_preparation.is_partitioned) + this->unset_is_partitioned(); + if (!old_preparation.has_removed_orphaned_nodes || + !other_preparation.has_removed_orphaned_nodes) + this->unset_has_removed_orphaned_nodes(); + if (!old_preparation.has_removed_remote_elements || + !other_preparation.has_removed_remote_elements) + this->unset_has_removed_remote_elements(); } // In general we've just invalidated just about everything, and we'd @@ -946,7 +956,8 @@ UnstructuredMesh::~UnstructuredMesh () void UnstructuredMesh::find_neighbors (const bool reset_remote_elements, const bool reset_current_list, - const bool assert_valid) + const bool assert_valid, + const bool check_non_remote) { // We might actually want to run this on an empty mesh // (e.g. the boundary mesh for a nonexistent bcid!) @@ -988,7 +999,10 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements, // If we haven't yet found a neighbor on this side, try. // Even if we think our neighbor is remote, that // information may be out of date. - if (element->neighbor_ptr(ms) == nullptr || + // + // If we're only checking remote neighbors, after a + // redistribution, then we'll skip the non-remote ones + if ((element->neighbor_ptr(ms) == nullptr && check_non_remote) || element->neighbor_ptr(ms) == remote_elem) { // Get the key for the side of this element. Use the diff --git a/src/systems/equation_systems.C b/src/systems/equation_systems.C index f498bbb18a..38abf2f357 100644 --- a/src/systems/equation_systems.C +++ b/src/systems/equation_systems.C @@ -105,6 +105,14 @@ void EquationSystems::reinit_mesh () libmesh_assert_not_equal_to (n_sys, 0); + // Our DofMaps are going to expect a prepared mesh later, but our + // older codes might have used some utilities that mark a mesh + // unprepared without re-preparing it. At *this* point hopefully + // everybody's done changing the mesh and we can make sure it's + // prepared. + if (!_mesh.is_prepared()) + _mesh.complete_preparation(); + // Tell all the \p DofObject entities how many systems // there are. for (auto & node : _mesh.node_ptr_range()) diff --git a/src/systems/fem_system.C b/src/systems/fem_system.C index 6587b2d552..c1ef31b846 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()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_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()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_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()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_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 f4f764e307..0d41502b9a 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" @@ -251,6 +252,11 @@ void NonlinearImplicitSystem::assembly(bool get_residual, bool /*apply_heterogeneous_constraints*/, bool /*apply_no_constraints*/) { + libmesh_assert(this->get_mesh().is_prepared()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_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 a26c45e06f..ec314ac3e9 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 @@ -210,6 +211,11 @@ void System::init_data () MeshBase & mesh = this->get_mesh(); + // Hopefully the user or the EquationSystems prepared the mesh, but + // if not then we'd better do it now. + if (!mesh.is_prepared()) + mesh.complete_preparation(); + // Distribute the degrees of freedom on the mesh auto total_dofs = _dof_map->distribute_dofs (mesh); @@ -556,6 +562,11 @@ void System::assemble () // Log how long the user's assembly code takes LOG_SCOPE("assemble()", "System"); + libmesh_assert(this->get_mesh().is_prepared()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); +#endif + // Call the user-specified assembly function this->user_assembly(); } @@ -567,6 +578,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()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_valid_is_prepared(this->get_mesh()); +#endif + // Call the user-specified quantity of interest function this->user_QOI(qoi_indices); } @@ -580,6 +596,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()); +#if defined(DEBUG) && !defined(LIBMESH_ENABLE_DEPRECATED) + MeshTools::libmesh_assert_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); diff --git a/tests/mesh/mesh_smoother_test.C b/tests/mesh/mesh_smoother_test.C index e828e9f653..7f499c9761 100644 --- a/tests/mesh/mesh_smoother_test.C +++ b/tests/mesh/mesh_smoother_test.C @@ -592,12 +592,15 @@ 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(); } + // We may have just invalidated the get_mesh_subdomains() cache by + // adding a new one, and libMesh marked our spatial_dimension() + // cache as invalidated in redistribute(). Just rebuild our + // caches before we start asking for anything like + // elem_default_orders(). + mesh.cache_elem_data(); + // Get the mesh order const auto & elem_orders = mesh.elem_default_orders(); libmesh_error_msg_if(elem_orders.size() != 1,