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..7afeadb1b1c 100644
--- a/examples/fem_system/fem_system_ex3/fem_system_ex3.C
+++ b/examples/fem_system/fem_system_ex3/fem_system_ex3.C
@@ -173,6 +173,10 @@ int main (int argc, char ** argv)
             mesh.get_boundary_info().add_edge(elem, e, edge_boundary_id);
     }
 
+  // We're all done adding to the boundary_info; let's make sure its
+  // caches know about it.
+  mesh.get_boundary_info().regenerate_id_sets();
+
   // Create an equation systems object.
   EquationSystems equation_systems (mesh);
 
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..ba93de15483 100644
--- a/examples/fem_system/fem_system_ex4/fem_system_ex4.C
+++ b/examples/fem_system/fem_system_ex4/fem_system_ex4.C
@@ -157,6 +157,9 @@ int main (int argc, char ** argv)
     if (elem->dim() < dim)
       elem->subdomain_id() = 1;
 
+  // Make sure the mesh knows we added new subdomains.
+  mesh.cache_elem_data();
+
   mesh_refinement.uniformly_refine(coarserefinements);
 
   // Print information about the mesh to the screen.
diff --git a/examples/fem_system/fem_system_ex5/fem_system_ex5.C b/examples/fem_system/fem_system_ex5/fem_system_ex5.C
index fc6e33faa6b..d2e566b5a28 100644
--- a/examples/fem_system/fem_system_ex5/fem_system_ex5.C
+++ b/examples/fem_system/fem_system_ex5/fem_system_ex5.C
@@ -225,6 +225,9 @@ int main (int argc, char ** argv)
         }
     }
 
+  // We're done modifying the BoundaryInfo; get its caches up to date.
+  mesh.get_boundary_info().regenerate_id_sets();
+
   // Create an equation systems object.
   EquationSystems equation_systems (mesh);
 
diff --git a/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h b/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h
index 6b2e1ca5d60..229c6c6efb8 100644
--- a/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h
+++ b/examples/miscellaneous/miscellaneous_ex9/augment_sparsity_on_interface.h
@@ -56,6 +56,13 @@ class AugmentSparsityOnInterface : public GhostingFunctor
                              boundary_id_type crack_boundary_lower,
                              boundary_id_type crack_boundary_upper);
 
+
+  virtual std::unique_ptr<GhostingFunctor> clone () const override
+  {
+    return std::make_unique<AugmentSparsityOnInterface>
+      (_mesh, _crack_boundary_lower, _crack_boundary_upper);
+  }
+
   /**
    * @return a const reference to the lower-to-upper element ID map.
    */
diff --git a/examples/subdomains/subdomains_ex1/subdomains_ex1.C b/examples/subdomains/subdomains_ex1/subdomains_ex1.C
index 8b5e121f4ac..10e85dfe585 100644
--- a/examples/subdomains/subdomains_ex1/subdomains_ex1.C
+++ b/examples/subdomains/subdomains_ex1/subdomains_ex1.C
@@ -243,6 +243,9 @@ int main (int argc, char ** argv)
         elem->subdomain_id() = 1;
     }
 
+  // Make sure the mesh knows we added new subdomains.
+  mesh.cache_elem_data();
+
   // Create an equation systems object.
   EquationSystems equation_systems (mesh);
 
diff --git a/examples/subdomains/subdomains_ex2/subdomains_ex2.C b/examples/subdomains/subdomains_ex2/subdomains_ex2.C
index abe827d038a..47966c45126 100644
--- a/examples/subdomains/subdomains_ex2/subdomains_ex2.C
+++ b/examples/subdomains/subdomains_ex2/subdomains_ex2.C
@@ -195,6 +195,9 @@ int main (int argc, char ** argv)
         elem->subdomain_id() = 1;
     }
 
+  // Make sure the mesh knows we added new subdomains.
+  mesh.cache_elem_data();
+
   // Print information about the mesh to the screen.
   mesh.print_info();
 
diff --git a/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C b/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
index f0fc5c9ba80..61cc9073804 100644
--- a/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
+++ b/examples/systems_of_equations/systems_of_equations_ex6/systems_of_equations_ex6.C
@@ -507,6 +507,10 @@ int main (int argc, char ** argv)
             mesh.get_boundary_info().add_edge(elem, e, EDGE_BOUNDARY_ID);
     }
 
+  // We're all done adding to the boundary_info; let's make sure its
+  // caches know about it.
+  mesh.get_boundary_info().regenerate_id_sets();
+
   // Create an equation systems object.
   EquationSystems equation_systems (mesh);
 
diff --git a/include/base/libmesh_common.h b/include/base/libmesh_common.h
index e2fe143a244..4d61e95b3d7 100644
--- a/include/base/libmesh_common.h
+++ b/include/base/libmesh_common.h
@@ -271,7 +271,7 @@ extern bool warned_about_auto_ptr;
   } while (0)
 
 // The libmesh_dbg_var() macro indicates that an argument to a function
-// is used only in debug mode (i.e., when NDEBUG is not defined).
+// is used only in debug and devel modes (i.e., when NDEBUG is not defined).
 #ifndef NDEBUG
 #define libmesh_dbg_var(var) var
 #else
diff --git a/src/geom/face_inf_quad.C b/src/geom/face_inf_quad.C
index 6070eed4fa9..6abfb1f21a6 100644
--- a/src/geom/face_inf_quad.C
+++ b/src/geom/face_inf_quad.C
@@ -184,7 +184,9 @@ bool InfQuad::is_child_on_side(const unsigned int c,
   libmesh_assert_less (c, this->n_children());
   libmesh_assert_less (s, this->n_sides());
 
-  return (s == 0 || s == c+1);
+  return (s == 0 ||
+          (c == 0 && s == 2) ||
+          (c == 1 && s == 1));
 }
 
 
diff --git a/src/mesh/boundary_info.C b/src/mesh/boundary_info.C
index e8561d333cc..6801a86cf86 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;
@@ -511,9 +513,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())
@@ -3300,11 +3301,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;
@@ -3337,14 +3355,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)
@@ -3401,9 +3424,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);
@@ -3428,8 +3454,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/unstructured_mesh.C b/src/mesh/unstructured_mesh.C
index dfa28eb5047..935780b5ddb 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)
diff --git a/tests/base/overlapping_coupling_test.C b/tests/base/overlapping_coupling_test.C
index 2b71345a077..29f29802324 100644
--- a/tests/base/overlapping_coupling_test.C
+++ b/tests/base/overlapping_coupling_test.C
@@ -60,6 +60,17 @@ class OverlappingCouplingFunctor : public GhostingFunctor
 
   virtual ~OverlappingCouplingFunctor(){};
 
+  virtual std::unique_ptr<GhostingFunctor> clone () const override
+  {
+    auto clone_functor =
+      std::make_unique<OverlappingCouplingFunctor>(_system);
+
+    auto coupling_matrix =
+      std::make_unique<CouplingMatrix>(*_coupling_matrix);
+    clone_functor->set_coupling_matrix(coupling_matrix);
+    return clone_functor;
+  }
+
   void set_coupling_matrix (std::unique_ptr<CouplingMatrix> & coupling_matrix)
   { _coupling_matrix = std::move(coupling_matrix); }
 
diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C
index 578ff186835..23d23480f64 100644
--- a/tests/geom/elem_test.C
+++ b/tests/geom/elem_test.C
@@ -679,6 +679,21 @@ public:
                   for (const Node & node : child_side->node_ref_range())
                     CPPUNIT_ASSERT(parent_side->contains_point(node));
 
+                // We can at least check for sharing some node with
+                // a side, on both finite and infinite elements, for
+                // those children that aren't just the middle elements
+                // of triangular sides.
+                bool shares_a_side_node = false;
+                bool shares_a_vertex = false;
+                for (const Node & node : child_side->node_ref_range())
+                  if (parent_side->local_node(node.id()) != invalid_uint)
+                    shares_a_side_node = true;
+                for (const Node & node : elem->node_ref_range())
+                  if (parent->local_node(node.id()) < parent->n_vertices())
+                    shares_a_vertex = true;
+
+                CPPUNIT_ASSERT(shares_a_side_node || !shares_a_vertex);
+
                 if (elem->neighbor_ptr(s) && !elem->neighbor_ptr(s)->is_remote())
                   CPPUNIT_ASSERT_EQUAL(parent->child_neighbor(elem->neighbor_ptr(s)), elem);
               }
@@ -699,6 +714,15 @@ public:
                 if (!parent_edge->infinite())
                   for (const Node & node : child_edge->node_ref_range())
                     CPPUNIT_ASSERT(parent_edge->contains_point(node));
+
+                // We can at least check for sharing some node with
+                // an edge, on both finite and infinite elements
+                bool shares_an_edge_node = false;
+                for (const Node & node : child_edge->node_ref_range())
+                  if (parent_edge->local_node(node.id()) != invalid_uint)
+                    shares_an_edge_node = true;
+
+                CPPUNIT_ASSERT(shares_an_edge_node);
               }
           }
 
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);
diff --git a/tests/systems/systems_test.C b/tests/systems/systems_test.C
index a181cd628b6..edd31f31ecc 100644
--- a/tests/systems/systems_test.C
+++ b/tests/systems/systems_test.C
@@ -62,6 +62,11 @@ public:
   _mesh(mesh)
   {}
 
+  virtual std::unique_ptr<GhostingFunctor> clone () const override
+  {
+    return std::make_unique<AugmentSparsityOnNodes>(_mesh);
+  }
+
   /**
    * User-defined function to augment the sparsity pattern.
    */