diff --git a/include/Makefile.in b/include/Makefile.in
index b24d74be7f..2456d7c8e9 100644
--- a/include/Makefile.in
+++ b/include/Makefile.in
@@ -693,6 +693,7 @@ include_HEADERS = \
         geom/cell_tet.h \
         geom/cell_tet10.h \
         geom/cell_tet4.h \
+        geom/compare_elems_by_level.h \
         geom/edge.h \
         geom/edge_edge2.h \
         geom/edge_edge3.h \
diff --git a/include/geom/compare_elems_by_level.h b/include/geom/compare_elems_by_level.h
new file mode 100644
index 0000000000..e15ef5e5fa
--- /dev/null
+++ b/include/geom/compare_elems_by_level.h
@@ -0,0 +1,61 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+
+#ifndef LIBMESH_COMPARE_ELEMS_BY_LEVEL_H
+#define LIBMESH_COMPARE_ELEMS_BY_LEVEL_H
+
+
+// Local Includes
+#include "libmesh/elem.h"
+
+// C++ Includes
+#include <numeric>
+#include <set>
+
+
+
+
+namespace libMesh {
+
+/**
+ * Specific weak ordering for Elem *'s to be used in a set.
+ * We use the id, but first sort by level.  This guarantees
+ * when traversing the set from beginning to end the lower
+ * level (parent) elements are encountered first.
+ */
+struct CompareElemIdsByLevel
+{
+  bool operator()(const Elem * a,
+                  const Elem * b) const
+  {
+    libmesh_assert (a);
+    libmesh_assert (b);
+    const unsigned int
+      al = a->level(), bl = b->level();
+    const dof_id_type
+      aid = a->id(),   bid = b->id();
+
+    return (al == bl) ? aid < bid : al < bl;
+  }
+};
+
+
+} // namespace libMesh
+
+#endif // LIBMESH_COMPARE_ELEMS_BY_LEVEL_H
diff --git a/include/include_HEADERS b/include/include_HEADERS
index 15952c17bc..a3d2aedf6f 100644
--- a/include/include_HEADERS
+++ b/include/include_HEADERS
@@ -125,6 +125,7 @@ include_HEADERS =  \
         geom/cell_tet.h \
         geom/cell_tet10.h \
         geom/cell_tet4.h \
+        geom/compare_elems_by_level.h \
         geom/edge.h \
         geom/edge_edge2.h \
         geom/edge_edge3.h \
diff --git a/include/libmesh/Makefile.am b/include/libmesh/Makefile.am
index 7d90864412..420e3ed62b 100644
--- a/include/libmesh/Makefile.am
+++ b/include/libmesh/Makefile.am
@@ -114,6 +114,7 @@ BUILT_SOURCES = \
         cell_tet.h \
         cell_tet10.h \
         cell_tet4.h \
+        compare_elems_by_level.h \
         edge.h \
         edge_edge2.h \
         edge_edge3.h \
@@ -852,6 +853,9 @@ cell_tet10.h: $(top_srcdir)/include/geom/cell_tet10.h
 cell_tet4.h: $(top_srcdir)/include/geom/cell_tet4.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
+compare_elems_by_level.h: $(top_srcdir)/include/geom/compare_elems_by_level.h
+	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
+
 edge.h: $(top_srcdir)/include/geom/edge.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
diff --git a/include/libmesh/Makefile.in b/include/libmesh/Makefile.in
index fd14f4e735..e6831f8532 100644
--- a/include/libmesh/Makefile.in
+++ b/include/libmesh/Makefile.in
@@ -522,17 +522,17 @@ BUILT_SOURCES = auto_ptr.h default_coupling.h dirichlet_boundaries.h \
 	cell_inf_hex8.h cell_inf_prism.h cell_inf_prism12.h \
 	cell_inf_prism6.h cell_prism.h cell_prism15.h cell_prism18.h \
 	cell_prism6.h cell_pyramid.h cell_pyramid13.h cell_pyramid14.h \
-	cell_pyramid5.h cell_tet.h cell_tet10.h cell_tet4.h edge.h \
-	edge_edge2.h edge_edge3.h edge_edge4.h edge_inf_edge2.h elem.h \
-	elem_cutter.h elem_hash.h elem_quality.h elem_range.h face.h \
-	face_inf_quad.h face_inf_quad4.h face_inf_quad6.h face_quad.h \
-	face_quad4.h face_quad4_shell.h face_quad8.h face_quad9.h \
-	face_tri.h face_tri3.h face_tri3_shell.h \
-	face_tri3_subdivision.h face_tri6.h node.h node_elem.h \
-	node_range.h plane.h point.h reference_elem.h remote_elem.h \
-	side.h sphere.h stored_range.h surface.h abaqus_io.h \
-	boundary_info.h boundary_mesh.h checkpoint_io.h \
-	distributed_mesh.h ensight_io.h exodusII_io.h \
+	cell_pyramid5.h cell_tet.h cell_tet10.h cell_tet4.h \
+	compare_elems_by_level.h edge.h edge_edge2.h edge_edge3.h \
+	edge_edge4.h edge_inf_edge2.h elem.h elem_cutter.h elem_hash.h \
+	elem_quality.h elem_range.h face.h face_inf_quad.h \
+	face_inf_quad4.h face_inf_quad6.h face_quad.h face_quad4.h \
+	face_quad4_shell.h face_quad8.h face_quad9.h face_tri.h \
+	face_tri3.h face_tri3_shell.h face_tri3_subdivision.h \
+	face_tri6.h node.h node_elem.h node_range.h plane.h point.h \
+	reference_elem.h remote_elem.h side.h sphere.h stored_range.h \
+	surface.h abaqus_io.h boundary_info.h boundary_mesh.h \
+	checkpoint_io.h distributed_mesh.h ensight_io.h exodusII_io.h \
 	exodusII_io_helper.h fro_io.h gmsh_io.h gmv_io.h gnuplot_io.h \
 	inf_elem_builder.h matlab_io.h medit_io.h mesh.h mesh_base.h \
 	mesh_communication.h mesh_function.h mesh_generation.h \
@@ -1198,6 +1198,9 @@ cell_tet10.h: $(top_srcdir)/include/geom/cell_tet10.h
 cell_tet4.h: $(top_srcdir)/include/geom/cell_tet4.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
+compare_elems_by_level.h: $(top_srcdir)/include/geom/compare_elems_by_level.h
+	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
+
 edge.h: $(top_srcdir)/include/geom/edge.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
diff --git a/include/mesh/checkpoint_io.h b/include/mesh/checkpoint_io.h
index a3d9cdb5c2..57554ddf44 100644
--- a/include/mesh/checkpoint_io.h
+++ b/include/mesh/checkpoint_io.h
@@ -22,7 +22,7 @@
 
 
 // Local includes
-#include "libmesh/libmesh.h"
+#include "libmesh/compare_elems_by_level.h"
 #include "libmesh/mesh_input.h"
 #include "libmesh/mesh_output.h"
 #include "libmesh/parallel_object.h"
@@ -39,13 +39,13 @@ class Xdr;
 /**
  * The CheckpointIO class can be used to write simplified restart
  * files that can be used to restart simulations that have
- * crashed. Only N-to-N (procs) restart is supported with CheckpointIO
- * files.
+ * crashed.
  *
  * \author Benjamin Kirk
  * \author John Peterson
  * \author Derek Gaston
- * \date 2013
+ * \author Roy Stogner
+ * \date 2017
  */
 class CheckpointIO : public MeshInput<MeshBase>,
                      public MeshOutput<MeshBase>,
@@ -112,18 +112,24 @@ class CheckpointIO : public MeshInput<MeshBase>,
   std::string &       version ()       { return _version; }
 
   /**
-   * Get/Set the processor_id to use
+   * Get/Set the processor id or processor ids to use.
+   *
+   * The default processor_id to use is the processor_id() of the
+   * mesh.
    *
    * This is used for m->n parallel checkpoint file writing:
-   * You can force CheckpointIO to view the world as if it is on a particular
-   * processor_id by setting it here
+   * You can force CheckpointIO to write out different partitions of a
+   * mesh by setting which partitions to write from each processor here.
    */
-  const processor_id_type & current_processor_id() const { return _my_processor_id; }
-  processor_id_type & current_processor_id() { return _my_processor_id; }
+  const std::vector<processor_id_type> & current_processor_ids() const { return _my_processor_ids; }
+  std::vector<processor_id_type> & current_processor_ids() { return _my_processor_ids; }
 
   /**
    * Get/Set the n_processors to use
    *
+   * The default n_processors to use is the n_processors() of the
+   * mesh.
+   *
    * This is used for m->n parallel checkpoint file writing:
    * You can force CheckpointIO to view the world as if it contains this number of
    * processors by setting it here
@@ -136,39 +142,39 @@ class CheckpointIO : public MeshInput<MeshBase>,
   // Write Implementation
 
   /**
-   * Build up the elem list
+   * Write subdomain name information - NEW in 0.9.2 format
    */
-  void build_elem_list();
+  void write_subdomain_names(Xdr & io) const;
 
   /**
-   * Build up the node list
+   * Write the connectivity for part of a mesh
    */
-  void build_node_list();
+  void write_connectivity (Xdr & io,
+                           const std::set<const Elem *, CompareElemIdsByLevel> & elements) const;
 
   /**
-   * Write subdomain name information - NEW in 0.9.2 format
+   * Write the remote_elem neighbor links for part of a mesh
    */
-  void write_subdomain_names(Xdr & io) const;
+  void write_remote_elem (Xdr & io,
+                          const std::set<const Elem *, CompareElemIdsByLevel> & elements) const;
 
   /**
-   * Write the connectivity for a parallel, distributed mesh
+   * Write the nodal locations for part of a mesh
    */
-  void write_connectivity (Xdr & io) const;
+  void write_nodes (Xdr & io,
+                    const std::set<const Node *> & nodeset) const;
 
   /**
-   * Write the nodal locations for a parallel, distributed mesh
+   * Write the side boundary conditions for part of a mesh
    */
-  void write_nodes (Xdr & io) const;
+  void write_bcs (Xdr & io,
+                  const std::set<const Elem *, CompareElemIdsByLevel> & elements) const;
 
   /**
-   * Write the boundary conditions for a parallel, distributed mesh
+   * Write the nodal boundary conditions for part of a mesh
    */
-  void write_bcs (Xdr & io) const;
-
-  /**
-   * Write the boundary conditions for a parallel, distributed mesh
-   */
-  void write_nodesets (Xdr & io) const;
+  void write_nodesets (Xdr & io,
+                       const std::set<const Node *> & nodeset) const;
 
   /**
    * Write boundary names information (sideset and nodeset) - NEW in 0.9.2 format
@@ -188,6 +194,11 @@ class CheckpointIO : public MeshInput<MeshBase>,
    */
   void read_connectivity (Xdr & io);
 
+  /**
+   * Read the remote_elem neighbor links for a parallel, distributed mesh
+   */
+  void read_remote_elem (Xdr & io);
+
   /**
    * Read the nodal locations for a parallel, distributed mesh
    */
@@ -214,21 +225,19 @@ class CheckpointIO : public MeshInput<MeshBase>,
    * Implemented by looping over all the active elements and finding
    * the maximum level.
    */
-  unsigned int n_active_levels_on_processor(const MeshBase & mesh) const;
+  unsigned int n_active_levels_in(MeshBase::const_element_iterator begin,
+                                  MeshBase::const_element_iterator end) const;
 
   bool _binary;
   bool _parallel;
   std::string _version;
   unsigned int _mesh_dimension;
+  unsigned int _n_active_levels;
 
-  /// These are sets of IDs to make the lookup for boundary conditions simpler
-  std::set<largest_id_type> _local_elements;
-  std::set<largest_id_type> _nodes_connected_to_local_elements;
-
-  /// The processor_id to use
-  processor_id_type _my_processor_id;
+  // The processor ids to write
+  std::vector<processor_id_type> _my_processor_ids;
 
-  /// The number of processors to use
+  // The largest processor id to write
   processor_id_type _my_n_processors;
 };
 
diff --git a/include/mesh/distributed_mesh.h b/include/mesh/distributed_mesh.h
index e2e2dd2dd4..a2c2adafa2 100644
--- a/include/mesh/distributed_mesh.h
+++ b/include/mesh/distributed_mesh.h
@@ -118,6 +118,13 @@ class DistributedMesh : public UnstructuredMesh
   virtual bool is_serial () const libmesh_override
   { return _is_serial; }
 
+  /**
+   * Asserts that not all elements and nodes of the mesh necessarily
+   * exist on the current processor.
+   */
+  virtual void set_distributed ()
+  { _is_serial = false; }
+
   /**
    * @returns \p true if new elements and nodes can and should be
    * created in synchronization on all processors, \p false otherwise
diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h
index ff7b841068..11669090f2 100644
--- a/include/mesh/mesh_base.h
+++ b/include/mesh/mesh_base.h
@@ -134,6 +134,14 @@ class MeshBase : public ParallelObject
   virtual bool is_serial () const
   { return true; }
 
+  /**
+   * Asserts that not all elements and nodes of the mesh necessarily
+   * exist on the current processor.  Only valid to call on classes
+   * which can be created in a distributed form.
+   */
+  virtual void set_distributed ()
+  { libmesh_error(); }
+
   /**
    * @returns \p true if new elements and nodes can and should be
    * created in synchronization on all processors, \p false otherwise
diff --git a/include/mesh/mesh_communication.h b/include/mesh/mesh_communication.h
index 26f6baf55a..7161801c38 100644
--- a/include/mesh/mesh_communication.h
+++ b/include/mesh/mesh_communication.h
@@ -21,6 +21,7 @@
 #define LIBMESH_MESH_COMMUNICATION_H
 
 // Local Includes
+#include "libmesh/compare_elems_by_level.h"
 #include "libmesh/libmesh_common.h"
 #include "libmesh/mesh_tools.h"
 
@@ -231,6 +232,31 @@ class MeshCommunication
 };
 
 
+// Related utilities
+
+// Ask a mesh's ghosting functors to insert into a set all elements
+// that are either on or connected to processor id \p pid.  Ask only
+// for newly-coarsened elements if \p newly_coarsened_only is true.
+void query_ghosting_functors(const MeshBase & mesh,
+                             processor_id_type pid,
+                             bool newly_coarsened_only,
+                             std::set<const Elem *, CompareElemIdsByLevel> & connected_elements);
+
+// Take a set of elements and insert all elements' immediate
+// children as well.
+void connect_children(const MeshBase & mesh,
+                      processor_id_type pid,
+                      std::set<const Elem *, CompareElemIdsByLevel> & connected_elements);
+
+// Take a set of elements and insert all elements' ancestors and
+// subactive descendants as well.
+void connect_families(std::set<const Elem *, CompareElemIdsByLevel> & connected_elements);
+
+// Take a set of elements and create a set of connected nodes.
+void reconnect_nodes (const std::set<const Elem *, CompareElemIdsByLevel> & connected_elements,
+                      std::set<const Node *> & connected_nodes);
+
+
 
 //--------------------------------------------------------------
 // MeshCommunication inline members
diff --git a/src/apps/splitter.C b/src/apps/splitter.C
index 4e3e881443..e7d98a4f2a 100644
--- a/src/apps/splitter.C
+++ b/src/apps/splitter.C
@@ -19,7 +19,7 @@
 // Read in a Mesh file and write out partitionings of it that are suitable
 // for reading into a DistributedMesh
 #include "libmesh/libmesh.h"
-#include "libmesh/mesh.h"
+#include "libmesh/replicated_mesh.h"
 #include "libmesh/checkpoint_io.h"
 #include "libmesh/metis_partitioner.h"
 #include "libmesh/getpot.h"
@@ -59,7 +59,7 @@ int main (int argc, char ** argv)
 
   Parallel::Communicator & comm = init.comm();
 
-  Mesh mesh(init.comm());
+  ReplicatedMesh mesh(init.comm());
 
   libMesh::out << "Reading " << filename << std::endl;
 
@@ -142,19 +142,14 @@ int main (int argc, char ** argv)
         {
           libMesh::out << "Writing " << my_num_chunks << " Files" << std::endl;
 
+          CheckpointIO cpr(mesh);
+          cpr.current_processor_ids().clear();
           for (unsigned int i = my_first_chunk; i < my_first_chunk + my_num_chunks; i++)
-            {
-              libMesh::out << " " << 100 * (static_cast<double>(i) / static_cast<double>(my_num_chunks)) << "% Complete" << std::endl;
-
-              CheckpointIO cpr(mesh);
-              cpr.current_processor_id() = i;
-              cpr.current_n_processors() = n_procs;
-              cpr.binary() = true;
-              cpr.parallel() = true;
-              cpr.write(remove_extension(filename) + ".cpr");
-            }
-
-          libMesh::out << " 100% Complete" << std::endl;
+            cpr.current_processor_ids().push_back(i);
+          cpr.current_n_processors() = n_procs;
+          cpr.binary() = true;
+          cpr.parallel() = true;
+          cpr.write(remove_extension(filename) + ".cpr");
         }
     }
 
diff --git a/src/mesh/checkpoint_io.C b/src/mesh/checkpoint_io.C
index cc27381ccd..b5549153cf 100644
--- a/src/mesh/checkpoint_io.C
+++ b/src/mesh/checkpoint_io.C
@@ -27,19 +27,20 @@
 #include <sstream> // for ostringstream
 
 // Local includes
-#include "libmesh/xdr_io.h"
-#include "libmesh/xdr_cxx.h"
+#include "libmesh/boundary_info.h"
+#include "libmesh/distributed_mesh.h"
+#include "libmesh/elem.h"
 #include "libmesh/enum_xdr_mode.h"
+#include "libmesh/libmesh_logging.h"
 #include "libmesh/mesh_base.h"
+#include "libmesh/mesh_communication.h"
+#include "libmesh/mesh_tools.h"
 #include "libmesh/node.h"
-#include "libmesh/elem.h"
-#include "libmesh/boundary_info.h"
 #include "libmesh/parallel.h"
-#include "libmesh/mesh_tools.h"
 #include "libmesh/partitioner.h"
-#include "libmesh/libmesh_logging.h"
-#include "libmesh/mesh_communication.h"
-#include "libmesh/distributed_mesh.h"
+#include "libmesh/remote_elem.h"
+#include "libmesh/xdr_io.h"
+#include "libmesh/xdr_cxx.h"
 
 namespace libMesh
 {
@@ -53,7 +54,7 @@ CheckpointIO::CheckpointIO (MeshBase & mesh, const bool binary_in) :
   _binary             (binary_in),
   _parallel           (false),
   _version            ("checkpoint-1.1"),
-  _my_processor_id    (processor_id()),
+  _my_processor_ids   (1, processor_id()),
   _my_n_processors    (n_processors())
 {
 }
@@ -65,7 +66,7 @@ CheckpointIO::CheckpointIO (const MeshBase & mesh, const bool binary_in) :
   ParallelObject      (mesh),
   _binary             (binary_in),
   _parallel           (false),
-  _my_processor_id    (processor_id()),
+  _my_processor_ids   (1, processor_id()),
   _my_n_processors    (n_processors())
 {
 }
@@ -86,33 +87,28 @@ void CheckpointIO::write (const std::string & name)
   // convenient reference to our mesh
   const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
 
-  // Try to dynamic cast the mesh to see if it's a DistributedMesh object
-  // Note: Just using is_serial() is not good enough because the Mesh won't
-  // have been prepared yet when is when that flag gets set to false... sigh.
-  _parallel = _parallel || !mesh.is_replicated();
+  // FIXME: For backwards compatibility, we'll assume for now that we
+  // only want to write distributed meshes in parallel.  Later we can
+  // do a gather_to_zero() and support that case too.
+  _parallel = _parallel || !mesh.is_serial();
 
-  // If this is a serial mesh then we're only going to write it on processor 0
-  if (_parallel || _my_processor_id == 0)
+  // We'll write a header file from processor 0 to make it easier to do unambiguous
+  // restarts later:
+  if (this->processor_id() == 0)
     {
-      std::ostringstream file_name_stream;
-
-      file_name_stream << name;
-
-      if (_parallel)
-        file_name_stream << "-" << _my_n_processors << "-" << _my_processor_id;
-
-      Xdr io (file_name_stream.str(), this->binary() ? ENCODE : WRITE);
+      Xdr io (name, this->binary() ? ENCODE : WRITE);
 
       // write the version
       io.data(_version, "# version");
 
-      // Write out whether or not this is a serial mesh (helps with error checking on read)
+      // Write out the max mesh dimension for backwards compatibility
+      // with code that sets it independently of element dimensions
       {
         unsigned int mesh_dimension = mesh.mesh_dimension();
         io.data(mesh_dimension, "# dimensions");
       }
 
-      // Write out whether or not this is a serial mesh (helps with error checking on read)
+      // Write out whether or not this is serial output
       {
         unsigned int parallel = _parallel;
         io.data(parallel, "# parallel");
@@ -122,91 +118,87 @@ void CheckpointIO::write (const std::string & name)
       // so we can check it upon reading the file
       if (_parallel)
         {
-          largest_id_type n_procs = _my_n_processors;
+          processor_id_type n_procs = _my_n_processors;
           io.data(n_procs, "# n_procs");
         }
 
-      // Build the list of local elements
-      this->build_elem_list();
-
-      // Build the list of nodes connected to local elements
-      this->build_node_list();
-
       // write subdomain names
       this->write_subdomain_names(io);
-
-      // write the nodal locations
-      this->write_nodes (io);
-
-      // write connectivity
-      this->write_connectivity (io);
-
-      // write the boundary condition information
-      this->write_bcs (io);
-
-      // write the nodeset information
-      this->write_nodesets (io);
-
-      // close it up
-      io.close();
     }
 
-  // this->comm().barrier();
-}
-
-void CheckpointIO::build_elem_list()
-{
-  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
-
-  MeshBase::const_element_iterator it = mesh.elements_begin();
-  MeshBase::const_element_iterator end = mesh.elements_end();
+  // If this is a serial mesh written to a serial file then we're only
+  // going to write local data from processor 0.  If this is a mesh being
+  // written in parallel then we're going to write from every
+  // processor.
+  std::vector<processor_id_type> ids_to_write;
 
   if (_parallel)
     {
-      it  = mesh.pid_elements_begin(_my_processor_id);
-      end = mesh.pid_elements_end(_my_processor_id);
+      ids_to_write = _my_processor_ids;
+    }
+  else if (mesh.is_serial())
+    {
+      if (mesh.processor_id() == 0)
+        {
+          // placeholder
+          ids_to_write.push_back(0);
+        }
+    }
+  else
+    {
+      libmesh_error_msg("Cannot write serial checkpoint from distributed mesh");
     }
 
-  std::set<const Elem *> neighbors;
-
-  for (; it != end; ++it)
+  for (std::vector<processor_id_type>::const_iterator
+         id_it = ids_to_write.begin(), id_end = ids_to_write.end();
+         id_it != id_end; ++id_it)
     {
-      Elem * elem = *it;
+      const processor_id_type my_pid = *id_it;
 
-      _local_elements.insert(elem->id());
+      std::ostringstream file_name_stream;
 
-      // Make sure the neighbors vector is cleared out in case this Elem
-      // isn't active and we skip calling find_point_neighbors().
-      neighbors.clear();
+      file_name_stream << name << "-" << (_parallel ? _my_n_processors : 1) << "-" << my_pid;
 
-      // Also need to add in all the point neighbors of this element.
-      // Point neighbors can only be found for active elements...
-      if (elem->active())
-        elem->find_point_neighbors(neighbors);
+      Xdr io (file_name_stream.str(), this->binary() ? ENCODE : WRITE);
 
-      std::set<const Elem *>::iterator
-        set_it = neighbors.begin(),
-        set_end = neighbors.end();
+      std::set<const Elem *, CompareElemIdsByLevel> elements;
 
-      for (; set_it != set_end; ++set_it)
-        _local_elements.insert((*set_it)->id());
-    }
-}
+      // For serial files we write everything
+      if (!_parallel)
+        elements.insert(mesh.elements_begin(), mesh.elements_end());
+      // For parallel files we write what we're asked plus what the
+      // ghosting functors and mesh structure say is associated with
+      // what we're asked.
+      else
+        {
+          query_ghosting_functors(mesh, my_pid, false, elements);
+          connect_children(mesh, my_pid, elements);
+          connect_families(elements);
+        }
 
-void CheckpointIO::build_node_list()
-{
-  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
+      std::set<const Node *> connected_nodes;
+      reconnect_nodes(elements, connected_nodes);
 
-  std::set<largest_id_type>::iterator it  = _local_elements.begin();
-  const std::set<largest_id_type>::iterator end  = _local_elements.end();
+      // write the nodal locations
+      this->write_nodes (io, connected_nodes);
 
-  for (; it != end; ++it)
-    {
-      const Elem * elem = mesh.elem_ptr(*it);
+      // write connectivity
+      this->write_connectivity (io, elements);
 
-      for (unsigned int n = 0; n < elem->n_nodes(); n++)
-        _nodes_connected_to_local_elements.insert(elem->node_id(n));
+      // write remote_elem connectivity
+      this->write_remote_elem (io, elements);
+
+      // write the boundary condition information
+      this->write_bcs (io, elements);
+
+      // write the nodeset information
+      this->write_nodesets (io, connected_nodes);
+
+      // close it up
+      io.close();
     }
+
+  // this->comm().barrier();
 }
 
 void CheckpointIO::write_subdomain_names(Xdr & io) const
@@ -246,12 +238,10 @@ void CheckpointIO::write_subdomain_names(Xdr & io) const
 
 
 
-void CheckpointIO::write_nodes (Xdr & io) const
+void CheckpointIO::write_nodes (Xdr & io,
+                                const std::set<const Node *> & nodeset) const
 {
-  // convenient reference to our mesh
-  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
-
-  unsigned int n_nodes_here = _nodes_connected_to_local_elements.size();
+  dof_id_type n_nodes_here = nodeset.size();
 
   io.data(n_nodes_here, "# n_nodes on proc");
 
@@ -261,13 +251,10 @@ void CheckpointIO::write_nodes (Xdr & io) const
   // For the coordinates
   std::vector<Real> coords(LIBMESH_DIM);
 
-  std::set<largest_id_type>::iterator
-    it  = _nodes_connected_to_local_elements.begin(),
-    end = _nodes_connected_to_local_elements.end();
-
-  for (; it != end; ++it)
+  for (std::set<const Node *>::iterator it = nodeset.begin(),
+         end = nodeset.end(); it != end; ++it)
     {
-      const Node & node = mesh.node_ref(*it);
+      const Node & node = **it;
 
       id_pid[0] = node.id();
       id_pid[1] = node.processor_id();
@@ -296,94 +283,105 @@ void CheckpointIO::write_nodes (Xdr & io) const
 
 
 
-void CheckpointIO::write_connectivity (Xdr & io) const
+void CheckpointIO::write_connectivity (Xdr & io,
+                                       const std::set<const Elem *, CompareElemIdsByLevel> & elements) const
 {
   libmesh_assert (io.writing());
 
-  // convenient reference to our mesh
-  const MeshBase & mesh = MeshOutput<MeshBase>::mesh();
-
-  // We will only write active elements and their parents.
-  unsigned int n_active_levels = n_active_levels_on_processor(mesh);
-
-  std::vector<xdr_id_type> n_elem_at_level(n_active_levels, 0);
-
-  // Find the number of elements at each level
-  {
-    std::set<largest_id_type>::iterator it = _local_elements.begin();
-    const std::set<largest_id_type>::iterator end = _local_elements.end();
-
-    for (; it != end; ++it)
-      n_elem_at_level[mesh.elem_ptr(*it)->level()]++;
-  }
-
-  io.data(n_active_levels, "# n_active_levels");
-
   // Put these out here to reduce memory churn
   // id type pid subdomain_id parent_id
   std::vector<largest_id_type> elem_data(5);
   std::vector<largest_id_type> conn_data;
 
-  for (unsigned int level=0; level < n_active_levels; level++)
-    {
-      std::ostringstream comment;
-      comment << "# n_elem at level ";
-      comment << level ;
-      io.data (n_elem_at_level[level], comment.str().c_str());
+  dof_id_type n_elems_here = elements.size();
 
-      std::set<largest_id_type>::iterator it = _local_elements.begin();
-      const std::set<largest_id_type>::iterator end = _local_elements.end();
+  io.data(n_elems_here, "# number of elements");
 
-      for (; it != end; ++it)
-        {
-          const Elem & elem = mesh.elem_ref(*it);
-
-          // Only write elements of the current level.
-          if (elem.level() != level)
-            continue;
+  for (std::set<const Elem *, CompareElemIdsByLevel>::const_iterator it = elements.begin(),
+         end = elements.end(); it != end; ++it)
+    {
+      const Elem & elem = **it;
 
-          unsigned int n_nodes = elem.n_nodes();
+      unsigned int n_nodes = elem.n_nodes();
 
-          elem_data[0] = elem.id();
-          elem_data[1] = elem.type();
-          elem_data[2] = elem.processor_id();
-          elem_data[3] = elem.subdomain_id();
+      elem_data[0] = elem.id();
+      elem_data[1] = elem.type();
+      elem_data[2] = elem.processor_id();
+      elem_data[3] = elem.subdomain_id();
 
-          if (elem.parent() != libmesh_nullptr)
-            elem_data[4] = elem.parent()->id();
-          else
-            elem_data[4] = DofObject::invalid_processor_id;
+      if (elem.parent() != libmesh_nullptr)
+        elem_data[4] = elem.parent()->id();
+      else
+        elem_data[4] = DofObject::invalid_processor_id;
 
-          conn_data.resize(n_nodes);
+      conn_data.resize(n_nodes);
 
-          for (unsigned int i=0; i<n_nodes; i++)
-            conn_data[i] = elem.node_id(i);
+      for (unsigned int i=0; i<n_nodes; i++)
+        conn_data[i] = elem.node_id(i);
 
-          io.data_stream(&elem_data[0],
-                         cast_int<unsigned int>(elem_data.size()),
-                         cast_int<unsigned int>(elem_data.size()));
+      io.data_stream(&elem_data[0],
+                     cast_int<unsigned int>(elem_data.size()),
+                     cast_int<unsigned int>(elem_data.size()));
 
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
-          largest_id_type unique_id = elem.unique_id();
+      largest_id_type unique_id = elem.unique_id();
 
-          io.data(unique_id, "# unique id");
+      io.data(unique_id, "# unique id");
 #endif
 
 #ifdef LIBMESH_ENABLE_AMR
-          unsigned int p_level = elem.p_level();
+      unsigned int p_level = elem.p_level();
 
-          io.data(p_level, "# p_level");
+      io.data(p_level, "# p_level");
 #endif
-          io.data_stream(&conn_data[0],
-                         cast_int<unsigned int>(conn_data.size()),
-                         cast_int<unsigned int>(conn_data.size()));
+      io.data_stream(&conn_data[0],
+                     cast_int<unsigned int>(conn_data.size()),
+                     cast_int<unsigned int>(conn_data.size()));
+    }
+}
+
+
+void CheckpointIO::write_remote_elem (Xdr & io,
+                                      const std::set<const Elem *, CompareElemIdsByLevel> & elements) const
+{
+  libmesh_assert (io.writing());
+
+  // Find the number of remote_elem links
+  dof_id_type n_remote_elem = 0;
+  std::vector<dof_id_type> elem_ids;
+  std::vector<unsigned short int> elem_sides;
+
+  for (std::set<const Elem *, CompareElemIdsByLevel>::const_iterator it = elements.begin(),
+         end = elements.end(); it != end; ++it)
+    {
+      const Elem & elem = **it;
+
+      for (unsigned int n=0; n != elem.n_neighbors(); ++n)
+        {
+          const Elem * neigh = elem.neighbor_ptr(n);
+          if (neigh == remote_elem ||
+              (neigh && !elements.count(neigh)))
+            {
+              elem_ids.push_back(elem.id());
+              elem_sides.push_back(n);
+              n_remote_elem++;
+            }
         }
     }
+
+  io.data(n_remote_elem, "# n_remote_elem");
+
+  if (n_remote_elem)
+    {
+      io.data_stream(&elem_ids[0], n_remote_elem);
+      io.data_stream(&elem_sides[0], n_remote_elem);
+    }
 }
 
 
 
-void CheckpointIO::write_bcs (Xdr & io) const
+void CheckpointIO::write_bcs (Xdr & io,
+                              const std::set<const Elem *, CompareElemIdsByLevel> & elements) const
 {
   libmesh_assert (io.writing());
 
@@ -397,37 +395,32 @@ void CheckpointIO::write_bcs (Xdr & io) const
   write_bc_names(io, boundary_info, true);  // sideset names
   write_bc_names(io, boundary_info, false); // nodeset names
 
-  std::vector<dof_id_type> orig_element_id_list;
-  std::vector<unsigned short int> orig_side_list;
-  std::vector<boundary_id_type> orig_bc_id_list;
+  std::vector<dof_id_type> full_element_id_list;
+  std::vector<unsigned short int> full_side_list;
+  std::vector<boundary_id_type> full_bc_id_list;
 
-  boundary_info.build_side_list(orig_element_id_list, orig_side_list, orig_bc_id_list);
+  boundary_info.build_side_list(full_element_id_list, full_side_list, full_bc_id_list);
 
-  size_t num_orig_elements = orig_element_id_list.size();
+  std::size_t bc_size = full_element_id_list.size();
+  libmesh_assert_equal_to(bc_size, full_side_list.size());
+  libmesh_assert_equal_to(bc_size, full_bc_id_list.size());
 
   std::vector<dof_id_type> element_id_list;
   std::vector<unsigned short int> side_list;
   std::vector<boundary_id_type> bc_id_list;
 
-  element_id_list.reserve(num_orig_elements);
-  side_list.reserve(num_orig_elements);
-  bc_id_list.reserve(num_orig_elements);
+  element_id_list.reserve(bc_size);
+  side_list.reserve(bc_size);
+  bc_id_list.reserve(bc_size);
+
+  for (std::size_t i = 0; i != bc_size; ++i)
+    if (elements.count(mesh.elem_ptr(full_element_id_list[i])))
+      {
+        element_id_list.push_back(full_element_id_list[i]);
+        side_list.push_back(full_side_list[i]);
+        bc_id_list.push_back(full_bc_id_list[i]);
+      }
 
-  // Filter these lists to only include local elements
-  for (unsigned int i = 0; i < num_orig_elements; i++)
-    {
-      if (_local_elements.find(orig_element_id_list[i]) != _local_elements.end())
-        {
-          // Only write out BCs for truly local elements
-          if (_parallel &&
-              mesh.elem_ptr(orig_element_id_list[i])->processor_id() != _my_processor_id)
-            continue;
-
-          element_id_list.push_back(orig_element_id_list[i]);
-          side_list.push_back(orig_side_list[i]);
-          bc_id_list.push_back(orig_bc_id_list[i]);
-        }
-    }
 
   io.data(element_id_list, "# element ids for bcs");
   io.data(side_list, "# sides of elements for bcs");
@@ -436,7 +429,8 @@ void CheckpointIO::write_bcs (Xdr & io) const
 
 
 
-void CheckpointIO::write_nodesets (Xdr & io) const
+void CheckpointIO::write_nodesets (Xdr & io,
+                                   const std::set<const Node *> & nodeset) const
 {
   libmesh_assert (io.writing());
 
@@ -446,33 +440,26 @@ void CheckpointIO::write_nodesets (Xdr & io) const
   // and our boundary info object
   const BoundaryInfo & boundary_info = mesh.get_boundary_info();
 
-  std::vector<dof_id_type> orig_node_id_list;
-  std::vector<boundary_id_type> orig_bc_id_list;
+  std::vector<dof_id_type> full_node_id_list;
+  std::vector<boundary_id_type> full_bc_id_list;
 
-  boundary_info.build_node_list(orig_node_id_list, orig_bc_id_list);
+  boundary_info.build_node_list(full_node_id_list, full_bc_id_list);
 
-  size_t num_orig_nodes = orig_node_id_list.size();
+  std::size_t nodeset_size = full_node_id_list.size();
+  libmesh_assert_equal_to(nodeset_size, full_bc_id_list.size());
 
   std::vector<dof_id_type> node_id_list;
   std::vector<boundary_id_type> bc_id_list;
 
-  node_id_list.reserve(num_orig_nodes);
-  bc_id_list.reserve(num_orig_nodes);
+  node_id_list.reserve(nodeset_size);
+  bc_id_list.reserve(nodeset_size);
 
-  // Filter these lists to only include nodes connected to local elements
-  for (unsigned int i = 0; i < num_orig_nodes; i++)
-    {
-      if (_nodes_connected_to_local_elements.find(orig_node_id_list[i]) != _nodes_connected_to_local_elements.end())
-        {
-          // Only write out BCs for truly local nodes
-          if (_parallel &&
-              mesh.node_ptr(orig_node_id_list[i])->processor_id() != _my_processor_id)
-            continue;
-
-          node_id_list.push_back(orig_node_id_list[i]);
-          bc_id_list.push_back(orig_bc_id_list[i]);
-        }
-    }
+  for (std::size_t i = 0; i != nodeset_size; ++i)
+    if (nodeset.count(mesh.node_ptr(full_node_id_list[i])))
+      {
+        node_id_list.push_back(full_node_id_list[i]);
+        bc_id_list.push_back(full_bc_id_list[i]);
+      }
 
   io.data(node_id_list, "# node id list");
   io.data(bc_id_list, "# nodeset bc id list");
@@ -523,82 +510,129 @@ void CheckpointIO::read (const std::string & name)
 
   MeshBase & mesh = MeshInput<MeshBase>::mesh();
 
-  // Try to dynamic cast the mesh to see if it's a DistributedMesh object
-  // Note: Just using is_serial() is not good enough because the Mesh won't
-  // have been prepared yet when is when that flag gets set to false... sigh.
-  _parallel = _parallel || !mesh.is_replicated();
+  // Will this be a parallel input file?  With how many processors?  Stay tuned!
+  unsigned int input_parallel;
+  processor_id_type input_n_procs;
 
-  // If this is a serial mesh then we're going to only read it on processor 0 and broadcast it
-  if (_parallel || _my_processor_id == 0)
+  // We'll write a header file from processor 0 and broadcast.
+  if (this->processor_id() == 0)
     {
-      std::ostringstream file_name_stream;
-
-      file_name_stream << name;
-
-      if (_parallel)
-        file_name_stream << "-" << _my_n_processors << "-" << _my_processor_id;
-
       {
-        std::ifstream in (file_name_stream.str().c_str());
+        std::ifstream in (name.c_str());
 
         if (!in.good())
-          libmesh_error_msg("ERROR: cannot locate specified file:\n\t" << file_name_stream.str());
+          libmesh_error_msg("ERROR: cannot locate header file:\n\t" << name);
       }
 
-      Xdr io (file_name_stream.str(), this->binary() ? DECODE : READ);
+      Xdr io (name, this->binary() ? DECODE : READ);
 
-      // read the version
-      io.data (_version);
+      // read the version, but don't care about it
+      std::string input_version;
+      io.data(input_version);
 
       // read the dimension
       io.data (_mesh_dimension);
+      this->comm().broadcast(_mesh_dimension);
       mesh.set_mesh_dimension(_mesh_dimension);
 
-      // Check if the mesh we're reading is the same as the one that was written
-      {
-        unsigned int parallel;
-        io.data(parallel, "# parallel");
-
-        if (_parallel != parallel)
-          libmesh_error_msg("Attempted to read a " <<
-                            (parallel ? "parallel" : "non-parallel")
-                            << " (" << parallel << ')'
-                            << " checkpoint file with a " <<
-                            (_parallel ? "parallel" : "non-parallel")
-                            << " checkpoint object!");
-      }
+      // Read whether or not this is a parallel file
+      io.data(input_parallel);
+      this->comm().broadcast(input_parallel);
 
-      // If this is a parallel mesh then we need to check to ensure we're reading this on the same number of procs
-      if (_parallel)
+      // With how many processors?
+      if (input_parallel)
         {
-          largest_id_type n_procs;
-          io.data(n_procs, "# n_procs");
-
-          if (n_procs != _my_n_processors)
-            libmesh_error_msg("Attempted to utilize a checkpoint file on " << _my_n_processors << " processors but it was written using " << n_procs << "!!");
+          io.data(input_n_procs);
+          this->comm().broadcast(input_n_procs);
         }
+      else
+        input_n_procs = 1;
 
-      // read subdomain names
+      // read and broadcast subdomain names
       this->read_subdomain_names(io);
 
-      // read the nodal locations
-      this->read_nodes (io);
+      std::map<subdomain_id_type, std::string> & subdomain_map =
+        mesh.set_subdomain_name_map();
+
+      this->comm().broadcast(subdomain_map);
+    }
+  // We'll receive the header broadcast everywhere else.
+  else
+    {
+      this->comm().broadcast(_mesh_dimension);
+      mesh.set_mesh_dimension(_mesh_dimension);
+
+      this->comm().broadcast(input_parallel);
+      if (input_parallel)
+        this->comm().broadcast(input_n_procs);
+      else
+        input_n_procs = 1;
+
+      std::map<subdomain_id_type, std::string> & subdomain_map =
+        mesh.set_subdomain_name_map();
+      this->comm().broadcast(subdomain_map);
+    }
 
-      // read connectivity
-      this->read_connectivity (io);
 
-      // read the boundary conditions
-      this->read_bcs (io);
+  // If this is a serial read then we're going to only read the mesh
+  // on processor 0, then broadcast it
+  if ((input_parallel && !mesh.is_replicated()) || mesh.processor_id() == 0)
+    {
+      // If we're trying to read a parallel checkpoint file on a
+      // replicated mesh, we'll read every file on processor 0 so we
+      // can broadcast it later.  If we're on a distributed mesh then
+      // we'll read every id to it's own processor and we'll "wrap
+      // around" with any ids that exceed our processor count.
+      const processor_id_type begin_proc_id =
+        (input_parallel && !mesh.is_replicated()) ?
+        mesh.processor_id() : 0;
+      const processor_id_type stride =
+        (input_parallel && !mesh.is_replicated()) ?
+        mesh.n_processors() : 1;
+
+      for (processor_id_type proc_id = begin_proc_id; proc_id < input_n_procs; proc_id += stride)
+        {
+          std::ostringstream file_name_stream;
 
-      // read the nodesets
-      this->read_nodesets (io);
+          file_name_stream << name;
 
-      io.close();
+          file_name_stream << "-" << (input_parallel ? input_n_procs : 1) << "-" << proc_id;
+
+          {
+            std::ifstream in (file_name_stream.str().c_str());
+
+            if (!in.good())
+              libmesh_error_msg("ERROR: cannot locate specified file:\n\t" << file_name_stream.str());
+          }
+
+          Xdr io (file_name_stream.str(), this->binary() ? DECODE : READ);
+
+          // read the nodal locations
+          this->read_nodes (io);
+
+          // read connectivity
+          this->read_connectivity (io);
+
+          // read remote_elem connectivity
+          this->read_remote_elem (io);
+
+          // read the boundary conditions
+          this->read_bcs (io);
+
+          // read the nodesets
+          this->read_nodesets (io);
+
+          io.close();
+        }
     }
 
-  // If the mesh is serial then we only read it on processor 0 so we need to broadcast it
-  if (!_parallel)
+  // If the mesh was only read on processor 0 then we need to broadcast it
+  if (mesh.is_replicated())
     MeshCommunication().broadcast(mesh);
+  // If the mesh is really distributed then we need to make sure it
+  // knows that
+  else if (mesh.n_processors() > 1)
+    mesh.set_distributed();
 }
 
 
@@ -668,15 +702,37 @@ void CheckpointIO::read_nodes (Xdr & io)
       p(2) = coords[2];
 #endif
 
+      const dof_id_type id = cast_int<dof_id_type>(id_pid[0]);
+
+      // "Wrap around" if we see more processors than we're using.
+      processor_id_type pid =
+        cast_int<processor_id_type>(id_pid[1] % mesh.n_processors());
+
+      // If we already have this node (e.g. from another file, when
+      // reading multiple distributed CheckpointIO files into a
+      // ReplicatedMesh) then we don't want to add it again (because
+      // ReplicatedMesh can't handle that) but we do want to assert
+      // consistency between what we're reading and what we have.
+      const Node * old_node = mesh.query_node_ptr(id);
+
+      if (old_node)
+        {
+          libmesh_assert_equal_to(pid, old_node->processor_id());
+#ifdef LIBMESH_ENABLE_UNIQUE_ID
+          libmesh_assert_equal_to(unique_id, old_node->unique_id());
+#endif
+        }
+      else
+        {
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
-      Node * node =
+          Node * node =
 #endif
-        mesh.add_point(p, cast_int<dof_id_type>(id_pid[0]),
-                       cast_int<processor_id_type>(id_pid[1]));
+            mesh.add_point(p, id, pid);
 
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
-      node->set_unique_id() = unique_id;
+          node->set_unique_id() = unique_id;
 #endif
+        }
     }
 }
 
@@ -687,59 +743,82 @@ void CheckpointIO::read_connectivity (Xdr & io)
   // convenient reference to our mesh
   MeshBase & mesh = MeshInput<MeshBase>::mesh();
 
-  unsigned int n_active_levels;
-  io.data(n_active_levels, "# n_active_levels");
+  dof_id_type n_elems_here;
+  io.data(n_elems_here);
 
   // Keep track of the highest dimensional element we've added to the mesh
   unsigned int highest_elem_dim = 1;
 
-  for (unsigned int level=0; level < n_active_levels; level++)
+  for (unsigned int i=0; i<n_elems_here; i++)
     {
-      xdr_id_type n_elem_at_level = 0;
-      io.data (n_elem_at_level, "");
-
-      for (unsigned int i=0; i<n_elem_at_level; i++)
-        {
-          // id type pid subdomain_id parent_id
-          std::vector<largest_id_type> elem_data(5);
-          io.data_stream
-            (&elem_data[0], cast_int<unsigned int>(elem_data.size()),
-             cast_int<unsigned int>(elem_data.size()));
+      // id type pid subdomain_id parent_id
+      std::vector<largest_id_type> elem_data(5);
+      io.data_stream
+        (&elem_data[0], cast_int<unsigned int>(elem_data.size()),
+         cast_int<unsigned int>(elem_data.size()));
 
 #ifdef LIBMESH_ENABLE_UNIQUE_ID
-          largest_id_type unique_id = 0;
-          io.data(unique_id, "# unique id");
+      largest_id_type unique_id = 0;
+      io.data(unique_id, "# unique id");
 #endif
 
 #ifdef LIBMESH_ENABLE_AMR
-          unsigned int p_level = 0;
+      unsigned int p_level = 0;
 
-          io.data(p_level, "# p_level");
+      io.data(p_level, "# p_level");
 #endif
 
-          unsigned int n_nodes = Elem::type_to_n_nodes_map[elem_data[1]];
-
-          // Snag the node ids this element was connected to
-          std::vector<largest_id_type> conn_data(n_nodes);
-          io.data_stream
-            (&conn_data[0], cast_int<unsigned int>(conn_data.size()),
-             cast_int<unsigned int>(conn_data.size()));
-
-          const dof_id_type id                 =
-            cast_int<dof_id_type>      (elem_data[0]);
-          const ElemType elem_type             =
-            static_cast<ElemType>      (elem_data[1]);
-          const processor_id_type proc_id      =
-            cast_int<processor_id_type>(elem_data[2]);
-          const subdomain_id_type subdomain_id =
-            cast_int<subdomain_id_type>(elem_data[3]);
-          const dof_id_type parent_id          =
-            cast_int<dof_id_type>      (elem_data[4]);
-
-          Elem * parent =
-            (parent_id == DofObject::invalid_processor_id) ?
-            libmesh_nullptr : mesh.elem_ptr(parent_id);
+      unsigned int n_nodes = Elem::type_to_n_nodes_map[elem_data[1]];
+
+      // Snag the node ids this element was connected to
+      std::vector<largest_id_type> conn_data(n_nodes);
+      io.data_stream
+        (&conn_data[0], cast_int<unsigned int>(conn_data.size()),
+         cast_int<unsigned int>(conn_data.size()));
+
+      const dof_id_type id                 =
+        cast_int<dof_id_type>      (elem_data[0]);
+      const ElemType elem_type             =
+        static_cast<ElemType>      (elem_data[1]);
+      const processor_id_type proc_id      =
+        cast_int<processor_id_type>
+        (elem_data[2] % mesh.n_processors());
+      const subdomain_id_type subdomain_id =
+        cast_int<subdomain_id_type>(elem_data[3]);
+      const dof_id_type parent_id          =
+        cast_int<dof_id_type>      (elem_data[4]);
+
+      Elem * parent =
+        (parent_id == DofObject::invalid_processor_id) ?
+        libmesh_nullptr : mesh.elem_ptr(parent_id);
+
+      Elem * old_elem = mesh.query_elem_ptr(id);
+
+      // If we already have this element (e.g. from another file,
+      // when reading multiple distributed CheckpointIO files into
+      // a ReplicatedMesh) then we don't want to add it again
+      // (because ReplicatedMesh can't handle that) but we do want
+      // to assert consistency between what we're reading and what
+      // we have.
+      if (old_elem)
+        {
+          libmesh_assert_equal_to(elem_type, old_elem->type());
+          libmesh_assert_equal_to(proc_id, old_elem->processor_id());
+          libmesh_assert_equal_to(subdomain_id, old_elem->subdomain_id());
+          if (parent)
+            libmesh_assert_equal_to(parent, old_elem->parent());
+          else
+            libmesh_assert(!old_elem->parent());
+
+          libmesh_assert_equal_to(old_elem->n_nodes(), conn_data.size());
 
+          for (std::size_t n=0; n != conn_data.size(); ++n)
+            libmesh_assert_equal_to
+              (old_elem->node_id(n),
+               cast_int<dof_id_type>(conn_data[n]));
+        }
+      else
+        {
           // Create the element
           Elem * elem = Elem::build(elem_type, parent).release();
 
@@ -781,6 +860,34 @@ void CheckpointIO::read_connectivity (Xdr & io)
 }
 
 
+void CheckpointIO::read_remote_elem (Xdr & io)
+{
+  // convenient reference to our mesh
+  MeshBase & mesh = MeshInput<MeshBase>::mesh();
+
+  // Find the number of remote_elem links
+  dof_id_type n_remote_elem = 0;
+
+  io.data(n_remote_elem, "# n_remote_elem");
+
+  std::vector<dof_id_type> elem_ids(n_remote_elem);
+  std::vector<unsigned short int> elem_sides(n_remote_elem);
+
+  if (n_remote_elem)
+    {
+      io.data_stream(&elem_ids[0], n_remote_elem);
+      io.data_stream(&elem_sides[0], n_remote_elem);
+
+      for (std::size_t i=0; i != n_remote_elem; ++i)
+        {
+          Elem & elem = mesh.elem_ref(elem_ids[i]);
+          elem.set_neighbor(elem_sides[i],
+                            const_cast<RemoteElem *>(remote_elem));
+        }
+    }
+}
+
+
 
 void CheckpointIO::read_bcs (Xdr & io)
 {
@@ -855,15 +962,14 @@ void CheckpointIO::read_bc_names(Xdr & io, BoundaryInfo & info, bool is_sideset)
 }
 
 
-unsigned int CheckpointIO::n_active_levels_on_processor(const MeshBase & mesh) const
+unsigned int CheckpointIO::n_active_levels_in(MeshBase::const_element_iterator begin,
+                                              MeshBase::const_element_iterator end) const
 {
   unsigned int max_level = 0;
 
-  std::set<largest_id_type>::iterator it = _local_elements.begin();
-  const std::set<largest_id_type>::iterator end = _local_elements.end();
-
-  for (; it != end; ++it)
-    max_level = std::max(mesh.elem_ptr((*it))->level(), max_level);
+  for (MeshBase::const_element_iterator it = begin;
+         it != end; ++it)
+    max_level = std::max((*it)->level(), max_level);
 
   return max_level + 1;
 }
diff --git a/src/mesh/mesh_communication.C b/src/mesh/mesh_communication.C
index c722447dd2..7c7bf4f54b 100644
--- a/src/mesh/mesh_communication.C
+++ b/src/mesh/mesh_communication.C
@@ -49,28 +49,6 @@ namespace {
 
 using namespace libMesh;
 
-/**
- * Specific weak ordering for Elem *'s to be used in a set.
- * We use the id, but first sort by level.  This guarantees
- * when traversing the set from beginning to end the lower
- * level (parent) elements are encountered first.
- */
-struct CompareElemIdsByLevel
-{
-  bool operator()(const Elem * a,
-                  const Elem * b) const
-  {
-    libmesh_assert (a);
-    libmesh_assert (b);
-    const unsigned int
-      al = a->level(), bl = b->level();
-    const dof_id_type
-      aid = a->id(),   bid = b->id();
-
-    return (al == bl) ? aid < bid : al < bl;
-  }
-};
-
 struct SyncNeighbors
 {
   typedef std::vector<dof_id_type> datum;
@@ -147,7 +125,15 @@ struct SyncNeighbors
 };
 
 
-void query_ghosting_functors(MeshBase & mesh,
+}
+
+
+
+namespace libMesh
+{
+
+
+void query_ghosting_functors(const MeshBase & mesh,
                              processor_id_type pid,
                              bool newly_coarsened_only,
                              std::set<const Elem *, CompareElemIdsByLevel> & connected_elements)
@@ -197,7 +183,7 @@ void query_ghosting_functors(MeshBase & mesh,
 }
 
 
-void connect_children(MeshBase & mesh,
+void connect_children(const MeshBase & mesh,
                       processor_id_type pid,
                       std::set<const Elem *, CompareElemIdsByLevel> & connected_elements)
 {
@@ -311,11 +297,7 @@ void reconnect_nodes (const std::set<const Elem *, CompareElemIdsByLevel> & conn
     }
 }
 
-}
-
 
-namespace libMesh
-{
 
 
 // ------------------------------------------------------------
diff --git a/src/partitioning/parmetis_partitioner.C b/src/partitioning/parmetis_partitioner.C
index 12c4b2bf59..6ed7021d0c 100644
--- a/src/partitioning/parmetis_partitioner.C
+++ b/src/partitioning/parmetis_partitioner.C
@@ -116,6 +116,9 @@ void ParmetisPartitioner::_do_repartition (MeshBase & mesh,
   // Revert to METIS on one processor.
   if (mesh.n_processors() == 1)
     {
+      // Make sure the mesh knows it's serial
+      mesh.allgather();
+
       MetisPartitioner mp;
       mp.partition (mesh, n_sbdmns);
       return;
diff --git a/src/systems/equation_systems.C b/src/systems/equation_systems.C
index e76896a98e..5a29155c0f 100644
--- a/src/systems/equation_systems.C
+++ b/src/systems/equation_systems.C
@@ -783,7 +783,7 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
         // We used to simply call non_const_sys.solution->close()
         // here, but that is not allowed when the solution vector is
         // locked read-only, for example when printing the solution
-        // during during the middle of a solve...  So try to be a bit
+        // during the middle of a solve...  So try to be a bit
         // more careful about calling close() unnecessarily.
         libmesh_assert(this->comm().verify(non_const_sys.solution->closed()));
         if (!non_const_sys.solution->closed())
diff --git a/src/utils/xdr_cxx.C b/src/utils/xdr_cxx.C
index adc809634f..dde36e86f8 100644
--- a/src/utils/xdr_cxx.C
+++ b/src/utils/xdr_cxx.C
@@ -564,6 +564,9 @@ xdrproc_t xdr_translator<long int>() { return (xdrproc_t)(xdr_long); }
 template <>
 xdrproc_t xdr_translator<unsigned long int>() { return (xdrproc_t)(xdr_u_long); }
 
+template <>
+xdrproc_t xdr_translator<long long>() { return (xdrproc_t)(xdr_longlong_t); }
+
 template <>
 xdrproc_t xdr_translator<unsigned long long>() { return (xdrproc_t)(xdr_u_longlong_t); }
 
@@ -779,25 +782,12 @@ void Xdr::data_stream (T * val, const unsigned int len, const unsigned int line_
 
         unsigned int size_of_type = cast_int<unsigned int>(sizeof(T));
 
-        if (size_of_type <= 4) // 32-bit types
-          {
-            xdr_vector(xdrs.get(),
-                       (char *) val,
-                       len,
-                       size_of_type,
-                       (xdrproc_t) xdr_u_int);
-          }
-        else // 64-bit types
-          {
-            xdr_vector(xdrs.get(),
-                       (char *) val,
-                       len,
-                       size_of_type,
-                       (xdrproc_t) xdr_u_hyper);
-          }
-
+        xdr_vector(xdrs.get(),
+                   (char *) val,
+                   len,
+                   size_of_type,
+                   xdr_translator<T>());
 #else
-
         libmesh_error_msg("ERROR: Functionality is not available.\n"    \
                           << "Make sure LIBMESH_HAVE_XDR is defined at build time\n" \
                           << "The XDR interface is not available in this installation");
@@ -814,28 +804,13 @@ void Xdr::data_stream (T * val, const unsigned int len, const unsigned int line_
 
         unsigned int size_of_type = cast_int<unsigned int>(sizeof(T));
 
-        if (size_of_type <= 4) // 32-bit types
-          {
-            if (len > 0)
-              xdr_vector(xdrs.get(),
-                         (char *) val,
-                         len,
-                         size_of_type,
-                         (xdrproc_t) xdr_u_int);
-          }
-        else // 64-bit types
-          {
-            if (len > 0)
-              xdr_vector(xdrs.get(),
-                         (char *) val,
-                         len,
-                         size_of_type,
-                         (xdrproc_t) xdr_u_hyper);
-
-          }
-
+        if (len > 0)
+          xdr_vector(xdrs.get(),
+                     (char *) val,
+                     len,
+                     size_of_type,
+                     xdr_translator<T>());
 #else
-
         libmesh_error_msg("ERROR: Functionality is not available.\n"    \
                           << "Make sure LIBMESH_HAVE_XDR is defined at build time\n" \
                           << "The XDR interface is not available in this installation");
@@ -1532,6 +1507,7 @@ template void Xdr::data<short int>                        (short int &,
 template void Xdr::data<unsigned long int>                (unsigned long int &,               const char *);
 template void Xdr::data<unsigned long long>               (unsigned long long &,              const char *);
 template void Xdr::data<long int>                         (long int &,                        const char *);
+template void Xdr::data<long long>                        (long long &,                       const char *);
 template void Xdr::data<char>                             (char &,                            const char *);
 template void Xdr::data<signed char>                      (signed char &,                     const char *);
 template void Xdr::data<unsigned char>                    (unsigned char &,                   const char *);
@@ -1547,6 +1523,7 @@ template void Xdr::data<std::vector<unsigned int> >       (std::vector<unsigned
 template void Xdr::data<std::vector<short int> >          (std::vector<short int> &,          const char *);
 template void Xdr::data<std::vector<unsigned short int> > (std::vector<unsigned short int> &, const char *);
 template void Xdr::data<std::vector<long int> >           (std::vector<long int> &,           const char *);
+template void Xdr::data<std::vector<long long> >          (std::vector<long long> &,          const char *);
 template void Xdr::data<std::vector<unsigned long int> >  (std::vector<unsigned long int> &,  const char *);
 template void Xdr::data<std::vector<unsigned long long> > (std::vector<unsigned long long> &, const char *);
 template void Xdr::data<std::vector<char> >               (std::vector<char> &,               const char *);
@@ -1559,7 +1536,10 @@ template void Xdr::data<std::vector<std::complex<float> > >  (std::vector<std::c
 template void Xdr::data<std::vector<std::complex<double> > > (std::vector<std::complex<double> > &, const char *);
 template void Xdr::data<std::vector<std::complex<long double> > > (std::vector<std::complex<long double> > &, const char *);
 template void Xdr::data<std::vector<std::string> >        (std::vector<std::string> &,        const char *);
+template void Xdr::data_stream<unsigned char>      (unsigned char * val,      const unsigned int len, const unsigned int line_break);
+template void Xdr::data_stream<short int>          (short int * val,          const unsigned int len, const unsigned int line_break);
 template void Xdr::data_stream<int>                (int * val,                const unsigned int len, const unsigned int line_break);
+template void Xdr::data_stream<long long>          (long long * val,          const unsigned int len, const unsigned int line_break);
 template void Xdr::data_stream<unsigned short int> (unsigned short int * val, const unsigned int len, const unsigned int line_break);
 template void Xdr::data_stream<unsigned int>       (unsigned int * val,       const unsigned int len, const unsigned int line_break);
 template void Xdr::data_stream<unsigned long int>  (unsigned long int * val,  const unsigned int len, const unsigned int line_break);
diff --git a/tests/mesh/checkpoint.C b/tests/mesh/checkpoint.C
index f6d0409589..8c6c0326f6 100644
--- a/tests/mesh/checkpoint.C
+++ b/tests/mesh/checkpoint.C
@@ -4,11 +4,11 @@
 #include <cppunit/TestCase.h>
 #include <libmesh/restore_warnings.h>
 
-#include "libmesh/mesh.h"
+#include "libmesh/distributed_mesh.h"
 #include "libmesh/replicated_mesh.h"
 #include "libmesh/checkpoint_io.h"
 #include "libmesh/mesh_generation.h"
-#include "libmesh/metis_partitioner.h"
+#include "libmesh/partitioner.h"
 
 #include "test_comm.h"
 
@@ -31,7 +31,14 @@ class CheckpointIOTest : public CppUnit::TestCase {
 public:
   CPPUNIT_TEST_SUITE( CheckpointIOTest );
 
-  CPPUNIT_TEST( testSplitter );
+  CPPUNIT_TEST( testAsciiDistRepSplitter );
+  CPPUNIT_TEST( testBinaryDistRepSplitter );
+  CPPUNIT_TEST( testAsciiRepDistSplitter );
+  CPPUNIT_TEST( testBinaryRepDistSplitter );
+  CPPUNIT_TEST( testAsciiRepRepSplitter );
+  CPPUNIT_TEST( testBinaryRepRepSplitter );
+  CPPUNIT_TEST( testAsciiDistDistSplitter );
+  CPPUNIT_TEST( testBinaryDistDistSplitter );
 
   CPPUNIT_TEST_SUITE_END();
 
@@ -47,18 +54,25 @@ public:
   }
 
   // Test that we can write multiple checkpoint files from a single processor.
-  void testSplitter()
+  template <typename MeshA, typename MeshB>
+  void testSplitter(bool binary, bool using_distmesh)
   {
-    // In this test, we partition the mesh into n_procs parts.
-    const unsigned int n_procs = 2;
+    // In this test, we partition the mesh into n_procs parts.  Don't
+    // try to partition a DistributedMesh into more parts than we have
+    // processors, though.
+    const unsigned int n_procs = using_distmesh ?
+      std::min(static_cast<unsigned int>(2), TestCommWorld->size()) :
+      2;
 
     // The number of elements in the original mesh.  For verification
     // later.
     dof_id_type original_n_elem = 0;
 
+    const std::string filename =
+      std::string("checkpoint_splitter.cp") + (binary ? "r" : "a");
+
     {
-      MetisPartitioner partitioner;
-      ReplicatedMesh mesh(*TestCommWorld);
+      MeshA mesh(*TestCommWorld);
 
       MeshTools::Generation::build_square(mesh,
                                           4,  4,
@@ -70,18 +84,20 @@ public:
       original_n_elem = mesh.n_elem();
 
       // Partition the mesh into n_procs pieces
-      partitioner.partition(mesh, n_procs);
-
-      // Write out checkpoint files for each processor.
-      for (unsigned i=0; i<n_procs; ++i)
-        {
-          CheckpointIO cpr(mesh);
-          cpr.current_processor_id() = i;
-          cpr.current_n_processors() = n_procs;
-          cpr.binary() = true;
-          cpr.parallel() = true;
-          cpr.write("checkpoint_splitter.cpr");
-        }
+      mesh.partition(n_procs);
+
+      // Write out checkpoint files for each piece.  Since on a
+      // ReplicatedMesh we might have more pieces than we do
+      // processors, some processors may have to write out more than
+      // one piece.
+      CheckpointIO cpr(mesh);
+      cpr.current_processor_ids().clear();
+      for (processor_id_type pid = mesh.processor_id(); pid < n_procs; pid += mesh.n_processors())
+        cpr.current_processor_ids().push_back(pid);
+      cpr.current_n_processors() = n_procs;
+      cpr.binary() = binary;
+      cpr.parallel() = true;
+      cpr.write(filename);
     }
 
     TestCommWorld->barrier();
@@ -89,26 +105,65 @@ public:
     // Test that we can read in the files we wrote and sum up to the
     // same total number of elements.
     {
-      unsigned int read_in_elements = 0;
+      MeshB mesh(*TestCommWorld);
+      CheckpointIO cpr(mesh);
+      cpr.binary() = binary;
+      cpr.read(filename);
 
-      for (unsigned i=0; i<n_procs; ++i)
+      std::size_t read_in_elements = 0;
+
+      for (unsigned pid=mesh.processor_id(); pid<n_procs; pid += mesh.n_processors())
         {
-          Mesh mesh(*TestCommWorld);
-          CheckpointIO cpr(mesh);
-          cpr.current_processor_id() = i;
-          cpr.current_n_processors() = n_procs;
-          cpr.binary() = true;
-          cpr.parallel() = true;
-          cpr.read("checkpoint_splitter.cpr");
-          read_in_elements += std::distance(mesh.pid_elements_begin(i),
-                                            mesh.pid_elements_end(i));
+          read_in_elements += std::distance(mesh.pid_elements_begin(pid),
+                                            mesh.pid_elements_end(pid));
         }
+      mesh.comm().sum(read_in_elements);
 
-      // Verify that we read in exactly as many elements on each proc as we started with.
+      // Verify that we read in exactly as many elements as we started with.
       CPPUNIT_ASSERT_EQUAL(static_cast<dof_id_type>(read_in_elements), original_n_elem);
     }
   }
 
+  void testAsciiDistRepSplitter()
+  {
+    testSplitter<DistributedMesh, ReplicatedMesh>(false, true);
+  }
+
+  void testBinaryDistRepSplitter()
+  {
+    testSplitter<DistributedMesh, ReplicatedMesh>(true, true);
+  }
+
+  void testAsciiRepDistSplitter()
+  {
+    testSplitter<ReplicatedMesh, DistributedMesh>(false, true);
+  }
+
+  void testBinaryRepDistSplitter()
+  {
+    testSplitter<ReplicatedMesh, DistributedMesh>(true, true);
+  }
+
+  void testAsciiRepRepSplitter()
+  {
+    testSplitter<ReplicatedMesh, ReplicatedMesh>(false, false);
+  }
+
+  void testBinaryRepRepSplitter()
+  {
+    testSplitter<ReplicatedMesh, ReplicatedMesh>(true, false);
+  }
+
+  void testAsciiDistDistSplitter()
+  {
+    testSplitter<DistributedMesh, DistributedMesh>(false, true);
+  }
+
+  void testBinaryDistDistSplitter()
+  {
+    testSplitter<DistributedMesh, DistributedMesh>(true, true);
+  }
+
 };
 
 CPPUNIT_TEST_SUITE_REGISTRATION( CheckpointIOTest );