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();
 
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);
 
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);
 
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 b1c40a44ee8..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
@@ -165,7 +173,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> & partitioner() { return _partitioner; }
 
@@ -194,7 +202,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 +210,134 @@ 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 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 _is_prepared; }
+  { return _preparation; }
+
+  /**
+   * \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
+   * Tells this we have done some operation where we should no longer
+   * consider ourself prepared.  This is a very coarse setting; it is
+   * generally more efficient 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.
+   *
+   * 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 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, 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; }
+
+  /**
+   * 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.
+   *
+   * 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; }
+
+  /**
+   * Tells this we have done some operation (e.g. refining elements
+   * with interior parents) which requires interior parent pointers to
+   * 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; }
+
+  /**
+   * 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.  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; }
+
+  /**
+   * 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.  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; }
+
+  /**
+   * 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 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; }
 
   /**
    * \returns \p true if all elements and nodes of the mesh
@@ -260,7 +385,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()
@@ -400,16 +527,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;
 
@@ -742,7 +870,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 +889,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 +908,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 +976,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 +1250,26 @@ class MeshBase : public ParallelObject
                                           const std::vector<T> * 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
+   *
+   * 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
@@ -1144,6 +1286,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 +1995,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:
 
@@ -1883,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
@@ -1891,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.
    */
@@ -1923,9 +2132,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/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h
index 04ea94fba0b..016b94bee2a 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
+
+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 <typename DofObjectSubclass>
+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
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/boundary_info.C b/src/mesh/boundary_info.C
index 29a3f530e4d..1e39b3d5755 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;
@@ -433,6 +435,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;
 }
 
 
@@ -511,9 +515,8 @@ void BoundaryInfo::sync (const std::set<boundary_id_type> & requested_boundary_i
   boundary_mesh.set_n_partitions() = _mesh->n_partitions();
 
   std::map<dof_id_type, dof_id_type> node_id_map;
-  std::map<std::pair<dof_id_type, unsigned char>, 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())
@@ -3297,11 +3300,28 @@ void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & requested_bo
                                  std::map<std::pair<dof_id_type, unsigned char>, dof_id_type> * side_id_map,
                                  const std::set<subdomain_id_type> & 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
+    <std::set<std::pair<dof_id_type, unsigned char>>>
+    side_id_set_by_level;
 
   // For avoiding extraneous element side construction
   ElemSideBuilder side_builder;
@@ -3334,14 +3354,19 @@ void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & 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)
@@ -3398,9 +3423,12 @@ void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & requested_bo
                     DofObject::invalid_processor_id)))
                 {
                   std::pair<dof_id_type, unsigned char> 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);
@@ -3425,8 +3453,21 @@ void BoundaryInfo::_find_id_maps(const std::set<boundary_id_type> & 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,
diff --git a/src/mesh/distributed_mesh.C b/src/mesh/distributed_mesh.C
index 24e352322d2..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<const DistributedMesh *>(&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<Elem *>(other_elem)))
-        return false;
+        return "_extra_ghost_elems entry";
     }
 
-  return true;
+  return "";
 }
 
 DistributedMesh::~DistributedMesh ()
@@ -205,7 +209,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 +295,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 +1617,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 +1784,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..cad02e90497 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
@@ -265,7 +265,54 @@ 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<std::string> 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 += "\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
   //
@@ -273,21 +320,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 (_is_prepared != other_mesh._is_prepared)
-    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
@@ -295,59 +340,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<bool>(_default_ghosting) != static_cast<bool>(other_mesh._default_ghosting))
-    return false;
+    return "_default_ghosting";
   if (static_cast<bool>(_partitioner) != static_cast<bool>(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<bool>(_disjoint_neighbor_boundary_pairs) !=
       static_cast<bool>(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 +383,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 +404,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 +420,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);
 }
 
 
@@ -585,6 +612,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<unsigned int>(_spatial_dimension);
 }
 
@@ -794,6 +823,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,18 +865,36 @@ 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();
 
   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
@@ -871,15 +920,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 +947,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,41 +965,47 @@ 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
+  // 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);
@@ -958,6 +1021,8 @@ MeshBase::reinit_ghosting_functors()
       libmesh_assert(gf);
       gf->mesh_reinit();
     }
+
+  _preparation.has_reinit_ghosting_functors = true;
 }
 
 void MeshBase::clear ()
@@ -965,8 +1030,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 +1748,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 +1953,8 @@ void MeshBase::cache_elem_data()
 #endif
         }
     }
+
+  _preparation.has_cached_elem_data = true;
 }
 
 
@@ -1903,10 +1972,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 +2077,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/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();
 }
 
 
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<DistributedMesh>(_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<EquationSystems>(*_mesh_copy);
   _system =
       &(_equation_systems->add_system<VariationalSmootherSystem>("variational_smoother_system"));
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/mesh/mesh_tools.C b/src/mesh/mesh_tools.C
index 7e98732dcb9..1ab3d632240 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<MeshBase> 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<MeshBase> 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;
+}
+
+
 }
 
 
@@ -1227,26 +1289,26 @@ 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())
+  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<MeshBase> 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);
-  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);
-  mesh_clone->allow_remote_element_removal(false);
-  mesh_clone->prepare_for_use();
-  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<MeshBase> mesh_clone = reprepared_mesh_clone(mesh);
 
   return (mesh == *mesh_clone);
 }
@@ -1255,6 +1317,40 @@ bool valid_is_prepared (const MeshBase & mesh)
 
 #ifndef NDEBUG
 
+
+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<MeshBase> 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.");
+}
+
+
+
+
+
 void libmesh_assert_equal_n_systems (const MeshBase & mesh)
 {
   LOG_SCOPE("libmesh_assert_equal_n_systems()", "MeshTools");
diff --git a/src/mesh/replicated_mesh.C b/src/mesh/replicated_mesh.C
index 059d3fcbe78..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<const ReplicatedMesh *>(&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 "";
 }
 
 
@@ -115,7 +119,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 +647,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 +826,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..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<MeshBase *>(&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 *>(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)
@@ -1256,15 +1298,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;
 }
 
 
diff --git a/src/systems/fem_system.C b/src/systems/fem_system.C
index 70cd06b6fac..1398ad6348e 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
+  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());
+#ifdef DEBUG
+  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());
+#ifdef DEBUG
+  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 f8273a11974..2897e00d891 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
+  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 e2a7f0f30d1..ac8c5fbba54 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
+  MeshTools::libmesh_assert_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
+  MeshTools::libmesh_assert_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
+  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 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
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);