diff --git a/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C b/examples/miscellaneous/miscellaneous_ex15/miscellaneous_ex15.C
index 6394f16648..f338719f03 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/include/mesh/distributed_mesh.h b/include/mesh/distributed_mesh.h
index 8f610c9b13..cf145508fd 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 eaee2a3b2d..e8b7d42e1b 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
@@ -983,10 +991,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.
@@ -2061,6 +2076,18 @@ class MeshBase : public ParallelObject
              has_reinit_ghosting_functors == other.has_reinit_ghosting_functors &&
              has_boundary_id_sets == other.has_boundary_id_sets;
     }
+
+    void libmesh_assert_consistent (const Parallel::Communicator & libmesh_dbg_var(comm)) {
+      libmesh_assert(comm.verify(is_partitioned));
+      libmesh_assert(comm.verify(has_synched_id_counts));
+      libmesh_assert(comm.verify(has_neighbor_ptrs));
+      libmesh_assert(comm.verify(has_cached_elem_data));
+      libmesh_assert(comm.verify(has_interior_parent_ptrs));
+      libmesh_assert(comm.verify(has_removed_remote_elements));
+      libmesh_assert(comm.verify(has_removed_orphaned_nodes));
+      libmesh_assert(comm.verify(has_reinit_ghosting_functors));
+      libmesh_assert(comm.verify(has_boundary_id_sets));
+    }
   };
 
 protected:
@@ -2101,7 +2128,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
@@ -2109,6 +2136,12 @@ class MeshBase : public ParallelObject
    */
   bool nodes_and_elements_equal(const MeshBase & other_mesh) const;
 
+  /**
+   *
+   */
+  std::string_view first_difference_from(const MeshBase & other_mesh) const;
+
+
   /**
    * \returns A writable reference to the number of partitions.
    */
diff --git a/include/mesh/mesh_tools.h b/include/mesh/mesh_tools.h
index 04ea94fba0..2f9f89b114 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 <typename DofObjectSubclass>
+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 9e2d5e950b..28ddb7eb4d 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 39c50c0cb0..f2f7e8f36e 100644
--- a/include/mesh/unstructured_mesh.h
+++ b/include/mesh/unstructured_mesh.h
@@ -294,7 +294,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/distributed_mesh.C b/src/mesh/distributed_mesh.C
index 9fdc2e3fae..825ec3bee0 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 ()
@@ -1037,9 +1041,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 9282a104d3..46f7895dbe 100644
--- a/src/mesh/mesh_base.C
+++ b/src/mesh/mesh_base.C
@@ -267,7 +267,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("operator==()", "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<std::string> 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
   //
@@ -275,21 +325,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
@@ -297,59 +345,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();
@@ -358,15 +388,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))
         {
@@ -379,14 +409,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*
@@ -395,13 +425,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);
 }
 
 
@@ -587,6 +617,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);
 }
 
@@ -860,10 +892,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
@@ -979,6 +1017,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);
diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C
index 0d72d8aaa5..08c5175a52 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 7fdcf9406a..1ab3d63224 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;
+}
+
+
 }
 
 
@@ -1229,24 +1291,24 @@ 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<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 2e8a4cf48f..8e26947adb 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 "";
 }
 
 
diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C
index ca5a90500b..5a1a734986 100644
--- a/src/mesh/unstructured_mesh.C
+++ b/src/mesh/unstructured_mesh.C
@@ -943,7 +943,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!)
@@ -985,7 +986,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/fem_system.C b/src/systems/fem_system.C
index 70cd06b6fa..51eda44382 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 f8273a1197..7f21f91ace 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());
+#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 a34a71983b..b9b018a080 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());
+#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();
 }
@@ -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());
+#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);
 }
@@ -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());
+#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);