diff --git a/include/numerics/distributed_vector.h b/include/numerics/distributed_vector.h
index c3f88eb2faf..0ac01cb21ce 100644
--- a/include/numerics/distributed_vector.h
+++ b/include/numerics/distributed_vector.h
@@ -57,14 +57,16 @@ class DistributedVector : public NumericVector<T>
    *  Dummy-Constructor. Dimension=0
    */
   explicit
-  DistributedVector (const ParallelType = AUTOMATIC);
+  DistributedVector (const ParallelType = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set dimension to \p n and initialize all elements with zero.
    */
   explicit
   DistributedVector (const numeric_index_type n,
-                     const ParallelType ptype = AUTOMATIC);
+                     const ParallelType ptype = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global dimension
@@ -72,7 +74,8 @@ class DistributedVector : public NumericVector<T>
    */
   DistributedVector (const numeric_index_type n,
 		     const numeric_index_type n_local,
-                     const ParallelType ptype = AUTOMATIC);
+                     const ParallelType ptype = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global
@@ -82,7 +85,8 @@ class DistributedVector : public NumericVector<T>
   DistributedVector (const numeric_index_type N,
 		     const numeric_index_type n_local,
 		     const std::vector<numeric_index_type>& ghost,
-                     const ParallelType ptype = AUTOMATIC);
+                     const ParallelType ptype = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Destructor, deallocates memory. Made virtual to allow
@@ -476,7 +480,8 @@ class DistributedVector : public NumericVector<T>
 // DistributedVector inline methods
 template <typename T>
 inline
-DistributedVector<T>::DistributedVector (const ParallelType ptype) :
+DistributedVector<T>::DistributedVector (const ParallelType ptype, const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm),
   _global_size      (0),
   _local_size       (0),
   _first_local_index(0),
@@ -490,7 +495,9 @@ DistributedVector<T>::DistributedVector (const ParallelType ptype) :
 template <typename T>
 inline
 DistributedVector<T>::DistributedVector (const numeric_index_type n,
-                                         const ParallelType ptype)
+                                         const ParallelType ptype,
+                                         const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n, false, ptype);
 }
@@ -501,7 +508,9 @@ template <typename T>
 inline
 DistributedVector<T>::DistributedVector (const numeric_index_type n,
 					 const numeric_index_type n_local,
-                                         const ParallelType ptype)
+                                         const ParallelType ptype,
+                                         const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n_local, false, ptype);
 }
@@ -513,7 +522,9 @@ inline
 DistributedVector<T>::DistributedVector (const numeric_index_type n,
 					 const numeric_index_type n_local,
 		                         const std::vector<numeric_index_type>& ghost,
-                                         const ParallelType ptype)
+                                         const ParallelType ptype,
+                                         const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n_local, ghost, false, ptype);
 }
@@ -571,7 +582,7 @@ void DistributedVector<T>::init (const numeric_index_type n,
 
   local_sizes[libMesh::processor_id()] = n_local;
 
-  CommWorld.sum(local_sizes);
+  this->communicator().sum(local_sizes);
 
   // _first_local_index is the sum of _local_size
   // for all processor ids less than ours
@@ -834,7 +845,7 @@ Real DistributedVector<T>::min () const
   for (numeric_index_type i = 1; i < _values.size(); ++i)
     local_min = std::min(libmesh_real(_values[i]), local_min);
 
-  CommWorld.min(local_min);
+  this->communicator().min(local_min);
 
   return local_min;
 }
@@ -857,7 +868,7 @@ Real DistributedVector<T>::max() const
   for (numeric_index_type i = 1; i < _values.size(); ++i)
     local_max = std::max(libmesh_real(_values[i]), local_max);
 
-  CommWorld.max(local_max);
+  this->communicator().max(local_max);
 
   return local_max;
 }
diff --git a/include/numerics/eigen_sparse_vector.h b/include/numerics/eigen_sparse_vector.h
index ade6413f12d..7613a255f3a 100644
--- a/include/numerics/eigen_sparse_vector.h
+++ b/include/numerics/eigen_sparse_vector.h
@@ -61,14 +61,16 @@ class EigenSparseVector : public NumericVector<T>
    *  Dummy-Constructor. Dimension=0
    */
   explicit
-  EigenSparseVector (const ParallelType = AUTOMATIC);
+  EigenSparseVector (const ParallelType = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set dimension to \p n and initialize all elements with zero.
    */
   explicit
   EigenSparseVector (const numeric_index_type n,
-                 const ParallelType = AUTOMATIC);
+                     const ParallelType = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global dimension
@@ -76,7 +78,8 @@ class EigenSparseVector : public NumericVector<T>
    */
   EigenSparseVector (const numeric_index_type n,
 		     const numeric_index_type n_local,
-		     const ParallelType = AUTOMATIC);
+		     const ParallelType = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global
@@ -86,7 +89,8 @@ class EigenSparseVector : public NumericVector<T>
   EigenSparseVector (const numeric_index_type N,
 		     const numeric_index_type n_local,
 		     const std::vector<numeric_index_type>& ghost,
-		     const ParallelType = AUTOMATIC);
+		     const ParallelType = AUTOMATIC,
+                     const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Destructor, deallocates memory. Made virtual to allow
@@ -474,7 +478,8 @@ class EigenSparseVector : public NumericVector<T>
 // EigenSparseVector inline methods
 template <typename T>
 inline
-EigenSparseVector<T>::EigenSparseVector (const ParallelType ptype)
+EigenSparseVector<T>::EigenSparseVector (const ParallelType ptype, const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->_type = ptype;
 }
@@ -484,7 +489,9 @@ EigenSparseVector<T>::EigenSparseVector (const ParallelType ptype)
 template <typename T>
 inline
 EigenSparseVector<T>::EigenSparseVector (const numeric_index_type n,
-					 const ParallelType ptype)
+					 const ParallelType ptype,
+                                         const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n, false, ptype);
 }
@@ -495,7 +502,9 @@ template <typename T>
 inline
 EigenSparseVector<T>::EigenSparseVector (const numeric_index_type n,
 					 const numeric_index_type n_local,
-					 const ParallelType ptype)
+					 const ParallelType ptype,
+                                         const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n_local, false, ptype);
 }
@@ -507,7 +516,9 @@ inline
 EigenSparseVector<T>::EigenSparseVector (const numeric_index_type N,
 					 const numeric_index_type n_local,
 					 const std::vector<numeric_index_type>& ghost,
-					 const ParallelType ptype)
+					 const ParallelType ptype,
+                                         const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(N, n_local, ghost, false, ptype);
 }
diff --git a/include/numerics/laspack_vector.h b/include/numerics/laspack_vector.h
index fbc2fb8c413..44ed55857f6 100644
--- a/include/numerics/laspack_vector.h
+++ b/include/numerics/laspack_vector.h
@@ -62,14 +62,16 @@ class LaspackVector : public NumericVector<T>
    *  Dummy-Constructor. Dimension=0
    */
   explicit
-  LaspackVector (const ParallelType = AUTOMATIC);
+  LaspackVector (const ParallelType = AUTOMATIC,
+                 const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set dimension to \p n and initialize all elements with zero.
    */
   explicit
   LaspackVector (const numeric_index_type n,
-                 const ParallelType = AUTOMATIC);
+                 const ParallelType = AUTOMATIC,
+                 const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global dimension
@@ -77,7 +79,8 @@ class LaspackVector : public NumericVector<T>
    */
   LaspackVector (const numeric_index_type n,
 		 const numeric_index_type n_local,
-                 const ParallelType = AUTOMATIC);
+                 const ParallelType = AUTOMATIC,
+                 const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global
@@ -87,7 +90,8 @@ class LaspackVector : public NumericVector<T>
   LaspackVector (const numeric_index_type N,
 		 const numeric_index_type n_local,
 		 const std::vector<numeric_index_type>& ghost,
-                 const ParallelType = AUTOMATIC);
+                 const ParallelType = AUTOMATIC,
+                 const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Destructor, deallocates memory. Made virtual to allow
@@ -463,7 +467,9 @@ class LaspackVector : public NumericVector<T>
 // LaspackVector inline methods
 template <typename T>
 inline
-LaspackVector<T>::LaspackVector (const ParallelType ptype)
+LaspackVector<T>::LaspackVector (const ParallelType ptype,
+                                 const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->_type = ptype;
 }
@@ -473,7 +479,9 @@ LaspackVector<T>::LaspackVector (const ParallelType ptype)
 template <typename T>
 inline
 LaspackVector<T>::LaspackVector (const numeric_index_type n,
-                                 const ParallelType ptype)
+                                 const ParallelType ptype,
+                                 const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n, false, ptype);
 }
@@ -484,7 +492,9 @@ template <typename T>
 inline
 LaspackVector<T>::LaspackVector (const numeric_index_type n,
 				 const numeric_index_type n_local,
-                                 const ParallelType ptype)
+                                 const ParallelType ptype,
+                                 const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(n, n_local, false, ptype);
 }
@@ -496,7 +506,9 @@ inline
 LaspackVector<T>::LaspackVector (const numeric_index_type N,
 	                         const numeric_index_type n_local,
 	                         const std::vector<numeric_index_type>& ghost,
-                                 const ParallelType ptype)
+                                 const ParallelType ptype,
+                                 const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm)
 {
   this->init(N, n_local, ghost, false, ptype);
 }
diff --git a/include/numerics/numeric_vector.h b/include/numerics/numeric_vector.h
index 377e06abb95..f4066503cb9 100644
--- a/include/numerics/numeric_vector.h
+++ b/include/numerics/numeric_vector.h
@@ -28,6 +28,7 @@
 #include "libmesh/id_types.h"
 #include "libmesh/reference_counted_object.h"
 #include "libmesh/libmesh.h"
+#include "libmesh/parallel_object.h"
 
 // C++ includes
 #include <cstddef>
@@ -54,7 +55,8 @@ template <typename T> class ShellMatrix;
  * @author Benjamin S. Kirk, 2003
  */
 template <typename T>
-class NumericVector : public ReferenceCountedObject<NumericVector<T> >
+class NumericVector : public ReferenceCountedObject<NumericVector<T> >,
+                      public ParallelObject
 {
 public:
 
@@ -62,14 +64,16 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T> >
    *  Dummy-Constructor. Dimension=0
    */
   explicit
-  NumericVector (const ParallelType ptype = AUTOMATIC);
+  NumericVector (const ParallelType ptype = AUTOMATIC,
+                 const Parallel::Communicator &comm=libMesh::CommWorld);
 
   /**
    * Constructor. Set dimension to \p n and initialize all elements with zero.
    */
   explicit
   NumericVector (const numeric_index_type n,
-                 const ParallelType ptype = AUTOMATIC);
+                 const ParallelType ptype = AUTOMATIC,
+                 const Parallel::Communicator &comm=libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global dimension
@@ -77,7 +81,8 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T> >
    */
   NumericVector (const numeric_index_type n,
 		 const numeric_index_type n_local,
-                 const ParallelType ptype = AUTOMATIC);
+                 const ParallelType ptype = AUTOMATIC,
+                 const Parallel::Communicator &comm=libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global
@@ -87,7 +92,8 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T> >
   NumericVector (const numeric_index_type N,
 		 const numeric_index_type n_local,
 		 const std::vector<numeric_index_type>& ghost,
-                 const ParallelType ptype = AUTOMATIC);
+                 const ParallelType ptype = AUTOMATIC,
+                 const Parallel::Communicator &comm=libMesh::CommWorld);
 
 public:
 
@@ -102,7 +108,8 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T> >
    * \p solver_package
    */
   static AutoPtr<NumericVector<T> >
-  build(const SolverPackage solver_package = libMesh::default_solver_package());
+  build(const SolverPackage solver_package = libMesh::default_solver_package(),
+        const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * @returns true if the vector has been initialized,
@@ -646,7 +653,8 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T> >
 
 template <typename T>
 inline
-NumericVector<T>::NumericVector (const ParallelType ptype) :
+NumericVector<T>::NumericVector (const ParallelType ptype, const Parallel::Communicator &comm) :
+  ParallelObject(comm),
   _is_closed(false),
   _is_initialized(false),
   _type(ptype)
@@ -658,7 +666,9 @@ NumericVector<T>::NumericVector (const ParallelType ptype) :
 template <typename T>
 inline
 NumericVector<T>::NumericVector (const numeric_index_type /*n*/,
-                                 const ParallelType ptype) :
+                                 const ParallelType ptype,
+                                 const Parallel::Communicator &comm) :
+  ParallelObject(comm),
   _is_closed(false),
   _is_initialized(false),
   _type(ptype)
@@ -673,7 +683,9 @@ template <typename T>
 inline
 NumericVector<T>::NumericVector (const numeric_index_type /*n*/,
 				 const numeric_index_type /*n_local*/,
-                                 const ParallelType ptype) :
+                                 const ParallelType ptype,
+                                 const Parallel::Communicator &comm) :
+  ParallelObject(comm),
   _is_closed(false),
   _is_initialized(false),
   _type(ptype)
@@ -689,7 +701,9 @@ inline
 NumericVector<T>::NumericVector (const numeric_index_type /*n*/,
 				 const numeric_index_type /*n_local*/,
 				 const std::vector<numeric_index_type>& /*ghost*/,
-                                 const ParallelType ptype) :
+                                 const ParallelType ptype,
+                                 const Parallel::Communicator &comm) :
+  ParallelObject(comm),
   _is_closed(false),
   _is_initialized(false),
   _type(ptype)
diff --git a/include/numerics/petsc_macro.h b/include/numerics/petsc_macro.h
index 8dc1a70eb30..546c8378a82 100644
--- a/include/numerics/petsc_macro.h
+++ b/include/numerics/petsc_macro.h
@@ -124,6 +124,8 @@ typedef enum { PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER} PetscCop
 #  define ISCreateLibMesh(comm,n,idx,mode,is) ISCreateGeneral((comm),(n),(idx),(mode),(is))
 #endif
 
+#define LIBMESH_CHKERRABORT(ierr) CHKERRABORT(this->communicator().get(), ierr);
+
 #else // LIBMESH_HAVE_PETSC
 
 #define PETSC_VERSION_LESS_THAN(major,minor,subminor) 1
diff --git a/include/numerics/petsc_vector.h b/include/numerics/petsc_vector.h
index 1f31e20e670..f32526352e8 100644
--- a/include/numerics/petsc_vector.h
+++ b/include/numerics/petsc_vector.h
@@ -66,14 +66,16 @@ class PetscVector : public NumericVector<T>
    *  Dummy-Constructor. Dimension=0
    */
   explicit
-  PetscVector (const ParallelType type = AUTOMATIC);
+  PetscVector (const ParallelType type = AUTOMATIC,
+               const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set dimension to \p n and initialize all elements with zero.
    */
   explicit
   PetscVector (const numeric_index_type n,
-               const ParallelType type = AUTOMATIC);
+               const ParallelType type = AUTOMATIC,
+               const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global dimension
@@ -81,7 +83,8 @@ class PetscVector : public NumericVector<T>
    */
   PetscVector (const numeric_index_type n,
 	       const numeric_index_type n_local,
-               const ParallelType type = AUTOMATIC);
+               const ParallelType type = AUTOMATIC,
+               const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global
@@ -91,7 +94,8 @@ class PetscVector : public NumericVector<T>
   PetscVector (const numeric_index_type N,
 	       const numeric_index_type n_local,
 	       const std::vector<numeric_index_type>& ghost,
-               const ParallelType type = AUTOMATIC);
+               const ParallelType type = AUTOMATIC,
+               const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor.  Creates a PetscVector assuming you already have a
@@ -100,7 +104,7 @@ class PetscVector : public NumericVector<T>
    * This allows ownership of \p v to remain with the original creator,
    * and to simply provide additional functionality with the PetscVector.
    */
-  PetscVector(Vec v);
+  PetscVector(Vec v, const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Destructor, deallocates memory. Made virtual to allow
@@ -598,12 +602,13 @@ class PetscVector : public NumericVector<T>
 
 template <typename T>
 inline
-PetscVector<T>::PetscVector (const ParallelType ptype)
-  : _array_is_present(false),
-    _local_form(NULL),
-    _values(NULL),
-    _global_to_local_map(),
-    _destroy_vec_on_exit(true)
+PetscVector<T>::PetscVector (const ParallelType ptype, const Parallel::Communicator &comm)
+    : NumericVector<T>(ptype, comm),
+      _array_is_present(false),
+      _local_form(NULL),
+      _values(NULL),
+      _global_to_local_map(),
+      _destroy_vec_on_exit(true)
 {
   this->_type = ptype;
 }
@@ -613,8 +618,10 @@ PetscVector<T>::PetscVector (const ParallelType ptype)
 template <typename T>
 inline
 PetscVector<T>::PetscVector (const numeric_index_type n,
-                             const ParallelType ptype)
-  : _array_is_present(false),
+                             const ParallelType ptype,
+                             const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm),
+    _array_is_present(false),
     _local_form(NULL),
     _values(NULL),
     _global_to_local_map(),
@@ -629,8 +636,10 @@ template <typename T>
 inline
 PetscVector<T>::PetscVector (const numeric_index_type n,
 			     const numeric_index_type n_local,
-                             const ParallelType ptype)
-  : _array_is_present(false),
+                             const ParallelType ptype,
+                             const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm),
+    _array_is_present(false),
     _local_form(NULL),
     _values(NULL),
     _global_to_local_map(),
@@ -646,8 +655,10 @@ inline
 PetscVector<T>::PetscVector (const numeric_index_type n,
 			     const numeric_index_type n_local,
 			     const std::vector<numeric_index_type>& ghost,
-                             const ParallelType ptype)
-  : _array_is_present(false),
+                             const ParallelType ptype,
+                             const Parallel::Communicator &comm)
+  : NumericVector<T>(ptype, comm),
+    _array_is_present(false),
     _local_form(NULL),
     _values(NULL),
     _global_to_local_map(),
@@ -662,8 +673,9 @@ PetscVector<T>::PetscVector (const numeric_index_type n,
 
 template <typename T>
 inline
-PetscVector<T>::PetscVector (Vec v)
-  : _array_is_present(false),
+PetscVector<T>::PetscVector (Vec v, const Parallel::Communicator &comm)
+  : NumericVector<T>(AUTOMATIC, comm),
+    _array_is_present(false),
     _local_form(NULL),
     _values(NULL),
     _global_to_local_map(),
@@ -1025,7 +1037,7 @@ template <typename T>
 inline
 AutoPtr<NumericVector<T> > PetscVector<T>::zero_clone () const
 {
-  AutoPtr<NumericVector<T> > cloned_vector (new PetscVector<T>);
+  AutoPtr<NumericVector<T> > cloned_vector (new PetscVector<T>(this->type(), this->communicator()));
 
   cloned_vector->init(*this);
 
@@ -1038,7 +1050,7 @@ template <typename T>
 inline
 AutoPtr<NumericVector<T> > PetscVector<T>::clone () const
 {
-  AutoPtr<NumericVector<T> > cloned_vector (new PetscVector<T>);
+  AutoPtr<NumericVector<T> > cloned_vector (new PetscVector<T>(this->type(), this->communicator()));
 
   cloned_vector->init(*this, true);
 
diff --git a/include/numerics/trilinos_epetra_vector.h b/include/numerics/trilinos_epetra_vector.h
index 3887392c19e..1875aa35a6a 100644
--- a/include/numerics/trilinos_epetra_vector.h
+++ b/include/numerics/trilinos_epetra_vector.h
@@ -65,14 +65,16 @@ class EpetraVector : public NumericVector<T>
    *  Dummy-Constructor. Dimension=0
    */
   explicit
-  EpetraVector (const ParallelType type = AUTOMATIC);
+  EpetraVector (const ParallelType type = AUTOMATIC,
+                const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set dimension to \p n and initialize all elements with zero.
    */
   explicit
   EpetraVector (const numeric_index_type n,
-                const ParallelType type = AUTOMATIC);
+                const ParallelType type = AUTOMATIC,
+                const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global dimension
@@ -80,7 +82,8 @@ class EpetraVector : public NumericVector<T>
    */
   EpetraVector (const numeric_index_type n,
 	        const numeric_index_type n_local,
-                const ParallelType type = AUTOMATIC);
+                const ParallelType type = AUTOMATIC,
+                const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor. Set local dimension to \p n_local, the global
@@ -90,7 +93,8 @@ class EpetraVector : public NumericVector<T>
   EpetraVector (const numeric_index_type N,
 		const numeric_index_type n_local,
 		const std::vector<numeric_index_type>& ghost,
-                const ParallelType type = AUTOMATIC);
+                const ParallelType type = AUTOMATIC,
+                const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Constructor.  Creates a EpetraVector assuming you already have a
@@ -99,7 +103,8 @@ class EpetraVector : public NumericVector<T>
    * This allows ownership of v to remain with the original creator,
    * and to simply provide additional functionality with the EpetraVector.
    */
-  EpetraVector(Epetra_Vector & v);
+  EpetraVector(Epetra_Vector & v,
+               const Parallel::Communicator &comm = libMesh::CommWorld);
 
   /**
    * Destructor, deallocates memory. Made virtual to allow
@@ -620,8 +625,10 @@ class EpetraVector : public NumericVector<T>
 
 template <typename T>
 inline
-EpetraVector<T>::EpetraVector (const ParallelType type)
-: _destroy_vec_on_exit(true),
+EpetraVector<T>::EpetraVector (const ParallelType type,
+                               const Parallel::Communicator &comm)
+: NumericVector<T>(type, comm),
+  _destroy_vec_on_exit(true),
   myFirstID_(0),
   myNumIDs_(0),
   myCoefs_(NULL),
@@ -641,8 +648,10 @@ EpetraVector<T>::EpetraVector (const ParallelType type)
 template <typename T>
 inline
 EpetraVector<T>::EpetraVector (const numeric_index_type n,
-                               const ParallelType type)
-: _destroy_vec_on_exit(true),
+                               const ParallelType type,
+                               const Parallel::Communicator &comm)
+: NumericVector<T>(type, comm),
+  _destroy_vec_on_exit(true),
   myFirstID_(0),
   myNumIDs_(0),
   myCoefs_(NULL),
@@ -664,8 +673,10 @@ template <typename T>
 inline
 EpetraVector<T>::EpetraVector (const numeric_index_type n,
 			       const numeric_index_type n_local,
-                               const ParallelType type)
-: _destroy_vec_on_exit(true),
+                               const ParallelType type,
+                               const Parallel::Communicator &comm)
+: NumericVector<T>(type, comm),
+  _destroy_vec_on_exit(true),
   myFirstID_(0),
   myNumIDs_(0),
   myCoefs_(NULL),
@@ -685,18 +696,20 @@ EpetraVector<T>::EpetraVector (const numeric_index_type n,
 
 template <typename T>
 inline
-EpetraVector<T>::EpetraVector(Epetra_Vector & v)
-: _destroy_vec_on_exit(false),
-  myFirstID_(0),
-  myNumIDs_(0),
-  myCoefs_(NULL),
-  nonlocalIDs_(NULL),
-  nonlocalElementSize_(NULL),
-  numNonlocalIDs_(0),
-  allocatedNonlocalLength_(0),
-  nonlocalCoefs_(NULL),
-  last_edit(0),
-  ignoreNonLocalEntries_(false)
+EpetraVector<T>::EpetraVector(Epetra_Vector & v,
+                              const Parallel::Communicator &comm)
+  : NumericVector<T>(AUTOMATIC, comm),
+    _destroy_vec_on_exit(false),
+    myFirstID_(0),
+    myNumIDs_(0),
+    myCoefs_(NULL),
+    nonlocalIDs_(NULL),
+    nonlocalElementSize_(NULL),
+    numNonlocalIDs_(0),
+    allocatedNonlocalLength_(0),
+    nonlocalCoefs_(NULL),
+    last_edit(0),
+    ignoreNonLocalEntries_(false)
 {
   _vec = &v;
 
@@ -721,8 +734,10 @@ inline
 EpetraVector<T>::EpetraVector (const numeric_index_type n,
 			       const numeric_index_type n_local,
 		               const std::vector<numeric_index_type>& ghost,
-                               const ParallelType type)
-: _destroy_vec_on_exit(true),
+                               const ParallelType type,
+                               const Parallel::Communicator &comm)
+: NumericVector<T>(AUTOMATIC, comm),
+  _destroy_vec_on_exit(true),
   myFirstID_(0),
   myNumIDs_(0),
   myCoefs_(NULL),
@@ -894,7 +909,7 @@ template <typename T>
 inline
 AutoPtr<NumericVector<T> > EpetraVector<T>::zero_clone () const
 {
-  AutoPtr<NumericVector<T> > cloned_vector (new EpetraVector<T>);
+  AutoPtr<NumericVector<T> > cloned_vector (new EpetraVector<T>(AUTOMATIC, this->communicator()));
 
   cloned_vector->init(*this);
 
@@ -907,7 +922,7 @@ template <typename T>
 inline
 AutoPtr<NumericVector<T> > EpetraVector<T>::clone () const
 {
-  AutoPtr<NumericVector<T> > cloned_vector (new EpetraVector<T>);
+  AutoPtr<NumericVector<T> > cloned_vector (new EpetraVector<T>(AUTOMATIC, this->communicator()));
 
   cloned_vector->init(*this, true);
 
diff --git a/src/base/dof_map_constraints.C b/src/base/dof_map_constraints.C
index 2fad1bd2e6b..ee048339536 100644
--- a/src/base/dof_map_constraints.C
+++ b/src/base/dof_map_constraints.C
@@ -1559,7 +1559,7 @@ void DofMap::enforce_constraints_exactly (const System &system,
   AutoPtr<NumericVector<Number> > v_built;
   if (v->type() == SERIAL)
     {
-      v_built = NumericVector<Number>::build();
+      v_built = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
       v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
       v_built->close();
 
@@ -1574,7 +1574,7 @@ void DofMap::enforce_constraints_exactly (const System &system,
     }
   else if (v->type() == PARALLEL)
     {
-      v_built = NumericVector<Number>::build();
+      v_built = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
       v_built->init (v->size(), v->size(), true, SERIAL);
       v->localize(*v_built);
       v_built->close();
diff --git a/src/error_estimation/adjoint_refinement_estimator.C b/src/error_estimation/adjoint_refinement_estimator.C
index 3033d0e6166..a7a57341ea4 100644
--- a/src/error_estimation/adjoint_refinement_estimator.C
+++ b/src/error_estimation/adjoint_refinement_estimator.C
@@ -153,7 +153,7 @@ void AdjointRefinementEstimator::estimate_error (const System& _system,
 
   // Copy the projected coarse grid solutions, which will be
   // overwritten by solve()
-  projected_solution = NumericVector<Number>::build().release();
+  projected_solution = NumericVector<Number>::build(libMesh::default_solver_package(), mesh.communicator()).release();
   projected_solution->init(system.solution->size(), true, SERIAL);
   system.solution->localize(*projected_solution,
 			    system.get_dof_map().get_send_list());
@@ -313,12 +313,12 @@ void AdjointRefinementEstimator::estimate_error (const System& _system,
   // Localize the global rhs and adjoint solution vectors (which might be shared on multiple processsors) onto a
   // local ghosted vector, this ensures each processor has all the dof_indices to compute an error indicator for
   // an element it owns
-  AutoPtr<NumericVector<Number> > localized_projected_residual = NumericVector<Number>::build();
+  AutoPtr<NumericVector<Number> > localized_projected_residual = NumericVector<Number>::build(libMesh::default_solver_package(), system.communicator());
   localized_projected_residual->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
   projected_residual.localize(*localized_projected_residual, system.get_dof_map().get_send_list());
 
   // Each adjoint solution will also require ghosting; for efficiency we'll reuse the same memory
-  AutoPtr<NumericVector<Number> > localized_adjoint_solution = NumericVector<Number>::build();
+  AutoPtr<NumericVector<Number> > localized_adjoint_solution = NumericVector<Number>::build(libMesh::default_solver_package(), system.communicator());
   localized_adjoint_solution->init(system.n_dofs(), system.n_local_dofs(), system.get_dof_map().get_send_list(), false, GHOSTED);
 
   // We will loop over each adjoint solution, localize that adjoint
diff --git a/src/error_estimation/exact_error_estimator.C b/src/error_estimation/exact_error_estimator.C
index dd9e37745d2..0e448b2174e 100644
--- a/src/error_estimation/exact_error_estimator.C
+++ b/src/error_estimation/exact_error_estimator.C
@@ -260,7 +260,7 @@ void ExactErrorEstimator::estimate_error (const System& system,
 
       // Prepare a global solution and a MeshFunction of the fine system if we need one
       AutoPtr<MeshFunction> fine_values;
-      AutoPtr<NumericVector<Number> > fine_soln = NumericVector<Number>::build();
+      AutoPtr<NumericVector<Number> > fine_soln = NumericVector<Number>::build(libMesh::default_solver_package(), mesh.communicator());
       if (_equation_systems_fine)
       {
 	const System& fine_system = _equation_systems_fine->get_system(system.name());
diff --git a/src/error_estimation/exact_solution.C b/src/error_estimation/exact_solution.C
index c929a0b0dc7..51c57f8aa68 100644
--- a/src/error_estimation/exact_solution.C
+++ b/src/error_estimation/exact_solution.C
@@ -556,7 +556,7 @@ void ExactSolution::_compute_error(const std::string& sys_name,
 
   // Prepare a global solution and a MeshFunction of the coarse system if we need one
   AutoPtr<MeshFunction> coarse_values;
-  AutoPtr<NumericVector<Number> > comparison_soln = NumericVector<Number>::build();
+  AutoPtr<NumericVector<Number> > comparison_soln = NumericVector<Number>::build(libMesh::default_solver_package(), _equation_systems.communicator());
   if (_equation_systems_fine)
     {
       const System& comparison_system
diff --git a/src/error_estimation/uniform_refinement_estimator.C b/src/error_estimation/uniform_refinement_estimator.C
index 5fcaf7245f0..b1d95cb2bdf 100644
--- a/src/error_estimation/uniform_refinement_estimator.C
+++ b/src/error_estimation/uniform_refinement_estimator.C
@@ -280,7 +280,7 @@ void UniformRefinementEstimator::_estimate_error (const EquationSystems* _es,
       // Copy the projected coarse grid solutions, which will be
       // overwritten by solve()
 //      projected_solutions[i] = system.solution->clone().release();
-      projected_solutions[i] = NumericVector<Number>::build().release();
+      projected_solutions[i] = NumericVector<Number>::build(libMesh::default_solver_package(), system.communicator()).release();
       projected_solutions[i]->init(system.solution->size(), true, SERIAL);
       system.solution->localize(*projected_solutions[i],
                                 system.get_dof_map().get_send_list());
diff --git a/src/numerics/distributed_vector.C b/src/numerics/distributed_vector.C
index d12bda1182b..4c9602a6178 100644
--- a/src/numerics/distributed_vector.C
+++ b/src/numerics/distributed_vector.C
@@ -53,7 +53,7 @@ T DistributedVector<T>::sum () const
   for (numeric_index_type i=0; i<local_size(); i++)
     local_sum += _values[i];
 
-  CommWorld.sum(local_sum);
+  this->communicator().sum(local_sum);
 
   return local_sum;
 }
@@ -75,7 +75,7 @@ Real DistributedVector<T>::l1_norm () const
   for (numeric_index_type i=0; i<local_size(); i++)
     local_l1 += std::abs(_values[i]);
 
-  CommWorld.sum(local_l1);
+  this->communicator().sum(local_l1);
 
   return local_l1;
 }
@@ -97,7 +97,7 @@ Real DistributedVector<T>::l2_norm () const
   for (numeric_index_type i=0; i<local_size(); i++)
     local_l2 += TensorTools::norm_sq(_values[i]);
 
-  CommWorld.sum(local_l2);
+  this->communicator().sum(local_l2);
 
   return std::sqrt(local_l2);
 }
@@ -123,7 +123,7 @@ Real DistributedVector<T>::linfty_norm () const
                                 // types are the same, as required
                                 // by std::max
 
-  CommWorld.max(local_linfty);
+  this->communicator().max(local_linfty);
 
   return local_linfty;
 }
@@ -365,7 +365,7 @@ T DistributedVector<T>::dot (const NumericVector<T>& V) const
     local_dot += this->_values[i] * v->_values[i];
 
   // The local dot products are now summed via MPI
-  CommWorld.sum(local_dot);
+  this->communicator().sum(local_dot);
 
   return local_dot;
 }
@@ -548,7 +548,7 @@ void DistributedVector<T>::localize (std::vector<T>& v_local) const
 
   v_local = this->_values;
 
-  CommWorld.allgather (v_local);
+  this->communicator().allgather (v_local);
 
 #ifndef LIBMESH_HAVE_MPI
   libmesh_assert_equal_to (local_size(), size());
@@ -570,7 +570,7 @@ void DistributedVector<T>::localize_to_one (std::vector<T>& v_local,
 
   v_local = this->_values;
 
-  CommWorld.gather (pid, v_local);
+  this->communicator().gather (pid, v_local);
 
 #ifndef LIBMESH_HAVE_MPI
   libmesh_assert_equal_to (local_size(), size());
diff --git a/src/numerics/numeric_vector.C b/src/numerics/numeric_vector.C
index 704d863329d..4e0947f41a5 100644
--- a/src/numerics/numeric_vector.C
+++ b/src/numerics/numeric_vector.C
@@ -43,7 +43,7 @@ namespace libMesh
 // Full specialization for Real datatypes
 template <typename T>
 AutoPtr<NumericVector<T> >
-NumericVector<T>::build(const SolverPackage solver_package)
+NumericVector<T>::build(const SolverPackage solver_package, const Parallel::Communicator &comm)
 {
   // Build the appropriate vector
   switch (solver_package)
@@ -53,7 +53,7 @@ NumericVector<T>::build(const SolverPackage solver_package)
 #ifdef LIBMESH_HAVE_LASPACK
     case LASPACK_SOLVERS:
       {
-	AutoPtr<NumericVector<T> > ap(new LaspackVector<T>);
+	AutoPtr<NumericVector<T> > ap(new LaspackVector<T>(AUTOMATIC, comm));
 	return ap;
       }
 #endif
@@ -62,7 +62,7 @@ NumericVector<T>::build(const SolverPackage solver_package)
 #ifdef LIBMESH_HAVE_PETSC
     case PETSC_SOLVERS:
       {
-	AutoPtr<NumericVector<T> > ap(new PetscVector<T>);
+	AutoPtr<NumericVector<T> > ap(new PetscVector<T>(AUTOMATIC, comm));
 	return ap;
       }
 #endif
@@ -71,7 +71,7 @@ NumericVector<T>::build(const SolverPackage solver_package)
 #ifdef LIBMESH_HAVE_TRILINOS
     case TRILINOS_SOLVERS:
       {
-	AutoPtr<NumericVector<T> > ap(new EpetraVector<T>);
+	AutoPtr<NumericVector<T> > ap(new EpetraVector<T>(AUTOMATIC, comm));
 	return ap;
       }
 #endif
@@ -80,15 +80,16 @@ NumericVector<T>::build(const SolverPackage solver_package)
 #ifdef LIBMESH_HAVE_EIGEN
     case EIGEN_SOLVERS:
       {
-	AutoPtr<NumericVector<T> > ap(new EigenSparseVector<T>);
+	AutoPtr<NumericVector<T> > ap(new EigenSparseVector<T>(AUTOMATIC, comm));
 	return ap;
       }
 #endif
 
 
     default:
-      AutoPtr<NumericVector<T> > ap(new DistributedVector<T>);
+      AutoPtr<NumericVector<T> > ap(new DistributedVector<T>(AUTOMATIC, comm));
       return ap;
+
     }
 
   AutoPtr<NumericVector<T> > ap(NULL);
@@ -119,7 +120,7 @@ int NumericVector<T>::compare (const NumericVector<T> &other_vector,
          && i<last_local_index());
 
   // Find the correct first differing index in parallel
-  CommWorld.min(first_different_i);
+  this->communicator().min(first_different_i);
 
   if (first_different_i == std::numeric_limits<int>::max())
     return -1;
@@ -152,7 +153,7 @@ int NumericVector<T>::local_relative_compare (const NumericVector<T> &other_vect
          && i<last_local_index());
 
   // Find the correct first differing index in parallel
-  CommWorld.min(first_different_i);
+  this->communicator().min(first_different_i);
 
   if (first_different_i == std::numeric_limits<int>::max())
     return -1;
@@ -188,7 +189,7 @@ int NumericVector<T>::global_relative_compare (const NumericVector<T> &other_vec
          && i<last_local_index());
 
   // Find the correct first differing index in parallel
-  CommWorld.min(first_different_i);
+  this->communicator().min(first_different_i);
 
   if (first_different_i == std::numeric_limits<int>::max())
     return -1;
@@ -317,7 +318,7 @@ Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indi
   for(; it!=it_end; ++it)
     norm += std::abs(v(*it));
 
-  CommWorld.sum(norm);
+  this->communicator().sum(norm);
 
   return norm;
 }
@@ -335,7 +336,7 @@ Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indi
   for(; it!=it_end; ++it)
     norm += TensorTools::norm_sq(v(*it));
 
-  CommWorld.sum(norm);
+  this->communicator().sum(norm);
 
   return std::sqrt(norm);
 }
@@ -357,7 +358,7 @@ Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> &
         norm = value;
     }
 
-  CommWorld.max(norm);
+  this->communicator().max(norm);
 
   return norm;
 }
diff --git a/src/numerics/petsc_vector.C b/src/numerics/petsc_vector.C
index 4ace34f8652..5e5ef9bc726 100644
--- a/src/numerics/petsc_vector.C
+++ b/src/numerics/petsc_vector.C
@@ -58,7 +58,7 @@ T PetscVector<T>::sum () const
   PetscScalar value=0.;
 
   ierr = VecSum (_vec, &value);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   return static_cast<T>(value);
 }
@@ -74,7 +74,7 @@ Real PetscVector<T>::l1_norm () const
   PetscReal value=0.;
 
   ierr = VecNorm (_vec, NORM_1, &value);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   return static_cast<Real>(value);
 }
@@ -91,7 +91,7 @@ Real PetscVector<T>::l2_norm () const
   PetscReal value=0.;
 
   ierr = VecNorm (_vec, NORM_2, &value);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   return static_cast<Real>(value);
 }
@@ -109,7 +109,7 @@ Real PetscVector<T>::linfty_norm () const
   PetscReal value=0.;
 
   ierr = VecNorm (_vec, NORM_INFINITY, &value);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   return static_cast<Real>(value);
 }
@@ -156,7 +156,7 @@ void PetscVector<T>::set (const numeric_index_type i, const T value)
   PetscScalar petsc_value = static_cast<PetscScalar>(value);
 
   ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   this->_is_closed = false;
 }
@@ -170,7 +170,7 @@ void PetscVector<T>::reciprocal()
 
   // VecReciprocal has been in PETSc since at least 2.3.3 days
   ierr = VecReciprocal(_vec);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 }
 
 
@@ -186,7 +186,7 @@ void PetscVector<T>::add (const numeric_index_type i, const T value)
   PetscScalar petsc_value = static_cast<PetscScalar>(value);
 
   ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   this->_is_closed = false;
 }
@@ -234,7 +234,7 @@ void PetscVector<T>::add_vector (const NumericVector<T>& V_in,
   // The const_cast<> is not elegant, but it is required since PETSc
   // is not const-correct.
   ierr = MatMultAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 }
 
 
@@ -266,7 +266,7 @@ void PetscVector<T>::add_vector_transpose (const NumericVector<T>& V_in,
   // The const_cast<> is not elegant, but it is required since PETSc
   // is not const-correct.
   ierr = MatMultTransposeAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 }
 
 
@@ -302,7 +302,7 @@ void PetscVector<T>::add_vector_conjugate_transpose (const NumericVector<T>& V_i
   // The const_cast<> is not elegant, but it is required since PETSc
   // is not const-correct.
   PetscErrorCode ierr = MatMultHermitianTranspose(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec);
-  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+  LIBMESH_CHKERRABORT(ierr);
 
   // Add the temporary copy to the matvec result
   this->add(1., *this_clone);
@@ -327,17 +327,17 @@ void PetscVector<T>::add (const T v_in)
       for (PetscInt i=0; i<n; i++)
 	{
 	  ierr = VecGetArray (_vec, &values);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 	  PetscInt ig = fli + i;
 
 	  PetscScalar value = (values[i] + v);
 
 	  ierr = VecRestoreArray (_vec, &values);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 	  ierr = VecSetValues (_vec, 1, &ig, &value, INSERT_VALUES);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	}
     }
   else
@@ -346,28 +346,28 @@ void PetscVector<T>::add (const T v_in)
 	 handling.  */
       Vec loc_vec;
       ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 
       PetscInt n=0;
       ierr = VecGetSize(loc_vec, &n);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 
       for (PetscInt i=0; i<n; i++)
 	{
 	  ierr = VecGetArray (loc_vec, &values);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 	  PetscScalar value = (values[i] + v);
 
 	  ierr = VecRestoreArray (loc_vec, &values);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 	  ierr = VecSetValues (loc_vec, 1, &i, &value, INSERT_VALUES);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	}
 
       ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
     }
 
   this->_is_closed = false;
@@ -402,11 +402,11 @@ void PetscVector<T>::add (const T a_in, const NumericVector<T>& v_in)
 #if PETSC_VERSION_LESS_THAN(2,3,0)
       // 2.2.x & earlier style
       ierr = VecAXPY(&a, v->_vec, _vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #else
       // 2.3.x & later style
       ierr = VecAXPY(_vec, a, v->_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #endif
     }
   else
@@ -414,24 +414,24 @@ void PetscVector<T>::add (const T a_in, const NumericVector<T>& v_in)
       Vec loc_vec;
       Vec v_loc_vec;
       ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
       ierr = VecGhostGetLocalForm (v->_vec,&v_loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 
 #if PETSC_VERSION_LESS_THAN(2,3,0)
       // 2.2.x & earlier style
       ierr = VecAXPY(&a, v_loc_vec, loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #else
       // 2.3.x & later style
       ierr = VecAXPY(loc_vec, a, v_loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #endif
 
       ierr = VecGhostRestoreLocalForm (v->_vec,&v_loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
       ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
     }
 }
 
@@ -498,31 +498,31 @@ void PetscVector<T>::scale (const T factor_in)
 #if PETSC_VERSION_LESS_THAN(2,3,0)
       // 2.2.x & earlier style
       ierr = VecScale(&factor, _vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #else
       // 2.3.x & later style
       ierr = VecScale(_vec, factor);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #endif
     }
   else
     {
       Vec loc_vec;
       ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 
 #if PETSC_VERSION_LESS_THAN(2,3,0)
       // 2.2.x & earlier style
       ierr = VecScale(&factor, loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #else
       // 2.3.x & later style
       ierr = VecScale(loc_vec, factor);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 #endif
 
       ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
     }
 }
 
@@ -536,19 +536,19 @@ void PetscVector<T>::abs()
   if(this->type() != GHOSTED)
     {
       ierr = VecAbs(_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
     }
   else
     {
       Vec loc_vec;
       ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 
       ierr = VecAbs(loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
 
       ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
     }
 }
 
@@ -568,7 +568,7 @@ T PetscVector<T>::dot (const NumericVector<T>& V) const
 
   // 2.3.x (at least) style.  Untested for previous versions.
   ierr = VecDot(this->_vec, v->_vec, &value);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   return static_cast<T>(value);
 }
@@ -589,7 +589,7 @@ T PetscVector<T>::indefinite_dot (const NumericVector<T>& V) const
 
   // 2.3.x (at least) style.  Untested for previous versions.
   ierr = VecTDot(this->_vec, v->_vec, &value);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   return static_cast<T>(value);
 }
@@ -612,31 +612,31 @@ PetscVector<T>::operator = (const T s_in)
 #if PETSC_VERSION_LESS_THAN(2,3,0)
 	  // 2.2.x & earlier style
 	  ierr = VecSet(&s, _vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 #else
 	  // 2.3.x & later style
 	  ierr = VecSet(_vec, s);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 #endif
 	}
       else
 	{
 	  Vec loc_vec;
 	  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 #if PETSC_VERSION_LESS_THAN(2,3,0)
 	  // 2.2.x & earlier style
 	  ierr = VecSet(&s, loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 #else
 	  // 2.3.x & later style
 	  ierr = VecSet(loc_vec, s);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 #endif
 
 	  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	}
     }
 
@@ -681,7 +681,7 @@ PetscVector<T>::operator = (const PetscVector<T>& v)
        causes no difficulty.  See discussion in libmesh-devel of
        June 24, 2010.  */
     ierr = VecCopy (v._vec, this->_vec);
-    CHKERRABORT(libMesh::COMM_WORLD,ierr);
+    LIBMESH_CHKERRABORT(ierr);
   }
   else
   {
@@ -695,24 +695,24 @@ PetscVector<T>::operator = (const PetscVector<T>& v)
       if(this->type() != GHOSTED)
       {
         ierr = VecCopy (v._vec, this->_vec);
-        CHKERRABORT(libMesh::COMM_WORLD,ierr);
+        LIBMESH_CHKERRABORT(ierr);
       }
       else
       {
         Vec loc_vec;
         Vec v_loc_vec;
         ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-        CHKERRABORT(libMesh::COMM_WORLD,ierr);
+        LIBMESH_CHKERRABORT(ierr);
         ierr = VecGhostGetLocalForm (v._vec,&v_loc_vec);
-        CHKERRABORT(libMesh::COMM_WORLD,ierr);
+        LIBMESH_CHKERRABORT(ierr);
 
         ierr = VecCopy (v_loc_vec, loc_vec);
-        CHKERRABORT(libMesh::COMM_WORLD,ierr);
+        LIBMESH_CHKERRABORT(ierr);
 
         ierr = VecGhostRestoreLocalForm (v._vec,&v_loc_vec);
-        CHKERRABORT(libMesh::COMM_WORLD,ierr);
+        LIBMESH_CHKERRABORT(ierr);
         ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-        CHKERRABORT(libMesh::COMM_WORLD,ierr);
+        LIBMESH_CHKERRABORT(ierr);
       }
     }
   }
@@ -742,13 +742,13 @@ PetscVector<T>::operator = (const std::vector<T>& v)
   if (this->size() == v.size())
     {
       ierr = VecGetArray (_vec, &values);
- 	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ 	     LIBMESH_CHKERRABORT(ierr);
 
       for (numeric_index_type i=0; i<nl; i++)
 	values[i] =  static_cast<PetscScalar>(v[i+ioff]);
 
       ierr = VecRestoreArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
     }
 
   /**
@@ -760,13 +760,13 @@ PetscVector<T>::operator = (const std::vector<T>& v)
       libmesh_assert_equal_to (this->local_size(), v.size());
 
       ierr = VecGetArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
 
       for (numeric_index_type i=0; i<nl; i++)
 	values[i] = static_cast<PetscScalar>(v[i]);
 
       ierr = VecRestoreArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
     }
 
   // Make sure ghost dofs are up to date
@@ -799,41 +799,41 @@ void PetscVector<T>::localize (NumericVector<T>& v_local_in) const
   std::vector<PetscInt> idx(n); Utility::iota (idx.begin(), idx.end(), 0);
 
   // Create the index set & scatter object
-  ierr = ISCreateLibMesh(libMesh::COMM_WORLD, n, &idx[0], PETSC_USE_POINTER, &is);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+  ierr = ISCreateLibMesh(this->communicator().get(), n, &idx[0], PETSC_USE_POINTER, &is);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterCreate(_vec,          is,
 			  v_local->_vec, is,
 			  &scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   // Perform the scatter
 #if PETSC_VERSION_LESS_THAN(2,3,3)
 
   ierr = VecScatterBegin(_vec, v_local->_vec, INSERT_VALUES,
 			 SCATTER_FORWARD, scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterEnd  (_vec, v_local->_vec, INSERT_VALUES,
 			 SCATTER_FORWARD, scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 #else
   // API argument order change in PETSc 2.3.3
   ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
 			 INSERT_VALUES, SCATTER_FORWARD);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterEnd  (scatter, _vec, v_local->_vec,
                          INSERT_VALUES, SCATTER_FORWARD);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 #endif
 
   // Clean up
   ierr = LibMeshISDestroy (&is);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = LibMeshVecScatterDestroy(&scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   // Make sure ghost dofs are up to date
   if (v_local->type() == GHOSTED)
@@ -882,17 +882,17 @@ void PetscVector<T>::localize (NumericVector<T>& v_local_in,
 
   // Create the index set & scatter object
   if (idx.empty())
-    ierr = ISCreateLibMesh(libMesh::COMM_WORLD,
+    ierr = ISCreateLibMesh(this->communicator().get(),
                            n_sl+this->local_size(), PETSC_NULL, PETSC_USE_POINTER, &is);
   else
-    ierr = ISCreateLibMesh(libMesh::COMM_WORLD,
+    ierr = ISCreateLibMesh(this->communicator().get(),
 			   n_sl+this->local_size(), &idx[0], PETSC_USE_POINTER, &is);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterCreate(_vec,          is,
 			  v_local->_vec, is,
 			  &scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
 
   // Perform the scatter
@@ -900,32 +900,32 @@ void PetscVector<T>::localize (NumericVector<T>& v_local_in,
 
   ierr = VecScatterBegin(_vec, v_local->_vec, INSERT_VALUES,
 			 SCATTER_FORWARD, scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterEnd  (_vec, v_local->_vec, INSERT_VALUES,
 			 SCATTER_FORWARD, scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
 #else
 
   // API argument order change in PETSc 2.3.3
   ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
                          INSERT_VALUES, SCATTER_FORWARD);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterEnd  (scatter, _vec, v_local->_vec,
                          INSERT_VALUES, SCATTER_FORWARD);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
 #endif
 
 
   // Clean up
   ierr = LibMeshISDestroy (&is);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   ierr = LibMeshVecScatterDestroy(&scatter);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   // Make sure ghost dofs are up to date
   if (v_local->type() == GHOSTED)
@@ -957,7 +957,7 @@ void PetscVector<T>::localize (const numeric_index_type first_local_idx,
 
   // Build a parallel vector, initialize it with the local
   // parts of (*this)
-  PetscVector<T> parallel_vec;
+  PetscVector<T> parallel_vec(PARALLEL, this->communicator());
 
   parallel_vec.init (my_size, my_local_size, true, PARALLEL);
 
@@ -972,45 +972,45 @@ void PetscVector<T>::localize (const numeric_index_type first_local_idx,
     Utility::iota (idx.begin(), idx.end(), first_local_idx);
 
     // Create the index set & scatter object
-    ierr = ISCreateLibMesh(libMesh::COMM_WORLD, my_local_size,
+    ierr = ISCreateLibMesh(this->communicator().get(), my_local_size,
                            my_local_size ? &idx[0] : NULL, PETSC_USE_POINTER, &is);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
     ierr = VecScatterCreate(_vec,              is,
 			    parallel_vec._vec, is,
 			    &scatter);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
     // Perform the scatter
 #if PETSC_VERSION_LESS_THAN(2,3,3)
 
     ierr = VecScatterBegin(_vec, parallel_vec._vec, INSERT_VALUES,
 			   SCATTER_FORWARD, scatter);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
     ierr = VecScatterEnd  (_vec, parallel_vec._vec, INSERT_VALUES,
 			   SCATTER_FORWARD, scatter);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
 #else
 
       // API argument order change in PETSc 2.3.3
     ierr = VecScatterBegin(scatter, _vec, parallel_vec._vec,
 			   INSERT_VALUES, SCATTER_FORWARD);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
     ierr = VecScatterEnd  (scatter, _vec, parallel_vec._vec,
 			   INSERT_VALUES, SCATTER_FORWARD);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
 #endif
 
     // Clean up
     ierr = LibMeshISDestroy (&is);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
 
     ierr = LibMeshVecScatterDestroy(&scatter);
-           CHKERRABORT(libMesh::COMM_WORLD,ierr);
+           LIBMESH_CHKERRABORT(ierr);
   }
 
   // localize like normal
@@ -1038,7 +1038,7 @@ void PetscVector<T>::localize (std::vector<T>& v_local) const
   v_local.resize(n, 0.);
 
   ierr = VecGetArray (_vec, &values);
-	 CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	 LIBMESH_CHKERRABORT(ierr);
 
   numeric_index_type ioff = first_local_index();
 
@@ -1046,9 +1046,9 @@ void PetscVector<T>::localize (std::vector<T>& v_local) const
     v_local[i+ioff] = static_cast<T>(values[i]);
 
   ierr = VecRestoreArray (_vec, &values);
-	 CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	 LIBMESH_CHKERRABORT(ierr);
 
-  CommWorld.sum(v_local);
+  this->communicator().sum(v_local);
 }
 
 
@@ -1075,13 +1075,13 @@ void PetscVector<Real>::localize_to_one (std::vector<Real>& v_local,
   if (n == nl)
     {
       ierr = VecGetArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
 
       for (PetscInt i=0; i<n; i++)
 	v_local[i] = static_cast<Real>(values[i]);
 
       ierr = VecRestoreArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
     }
 
   // otherwise multiple processors
@@ -1092,18 +1092,18 @@ void PetscVector<Real>::localize_to_one (std::vector<Real>& v_local,
 
       {
 	ierr = VecGetArray (_vec, &values);
-	       CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	       LIBMESH_CHKERRABORT(ierr);
 
 	for (PetscInt i=0; i<nl; i++)
 	  local_values[i+ioff] = static_cast<Real>(values[i]);
 
 	ierr = VecRestoreArray (_vec, &values);
-	       CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	       LIBMESH_CHKERRABORT(ierr);
       }
 
 
       MPI_Reduce (&local_values[0], &v_local[0], n, MPI_REAL, MPI_SUM,
-		  pid, libMesh::COMM_WORLD);
+		  pid, this->communicator().get());
     }
 }
 
@@ -1135,13 +1135,13 @@ void PetscVector<Complex>::localize_to_one (std::vector<Complex>& v_local,
   if (n == nl)
     {
       ierr = VecGetArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
 
       for (PetscInt i=0; i<n; i++)
 	v_local[i] = static_cast<Complex>(values[i]);
 
       ierr = VecRestoreArray (_vec, &values);
-	     CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	     LIBMESH_CHKERRABORT(ierr);
     }
 
   // otherwise multiple processors
@@ -1157,7 +1157,7 @@ void PetscVector<Complex>::localize_to_one (std::vector<Complex>& v_local,
 
       {
 	ierr = VecGetArray (_vec, &values);
-	       CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	       LIBMESH_CHKERRABORT(ierr);
 
 	// provide my local share to the real and imag buffers
 	for (PetscInt i=0; i<nl; i++)
@@ -1167,7 +1167,7 @@ void PetscVector<Complex>::localize_to_one (std::vector<Complex>& v_local,
 	  }
 
 	ierr = VecRestoreArray (_vec, &values);
-	       CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	       LIBMESH_CHKERRABORT(ierr);
       }
 
       /* have buffers of the real and imaginary part of v_local.
@@ -1181,11 +1181,11 @@ void PetscVector<Complex>::localize_to_one (std::vector<Complex>& v_local,
       // collect entries from other proc's in real_v_local, imag_v_local
       MPI_Reduce (&real_local_values[0], &real_v_local[0], n,
 		  MPI_REAL, MPI_SUM,
-		  pid, libMesh::COMM_WORLD);
+		  pid, this->communicator().get());
 
       MPI_Reduce (&imag_local_values[0], &imag_v_local[0], n,
 		  MPI_REAL, MPI_SUM,
-		  pid, libMesh::COMM_WORLD);
+		  pid, this->communicator().get());
 
       // copy real_v_local and imag_v_local to v_local
       for (PetscInt i=0; i<n; i++)
@@ -1226,7 +1226,7 @@ void PetscVector<T>::pointwise_mult (const NumericVector<T>& vec1,
       ierr = VecPointwiseMult(this->vec(),
 			      const_cast<PetscVector<T>*>(vec1_petsc)->vec(),
 			      const_cast<PetscVector<T>*>(vec2_petsc)->vec());
-      CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      LIBMESH_CHKERRABORT(ierr);
     }
   else
     {
@@ -1234,21 +1234,21 @@ void PetscVector<T>::pointwise_mult (const NumericVector<T>& vec1,
 	  Vec v1_loc_vec;
 	  Vec v2_loc_vec;
 	  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	  ierr = VecGhostGetLocalForm (const_cast<PetscVector<T>*>(vec1_petsc)->vec(),&v1_loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	  ierr = VecGhostGetLocalForm (const_cast<PetscVector<T>*>(vec2_petsc)->vec(),&v2_loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 	  ierr = VecPointwiseMult(loc_vec,v1_loc_vec,v2_loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 
 	  ierr = VecGhostRestoreLocalForm (const_cast<PetscVector<T>*>(vec1_petsc)->vec(),&v1_loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	  ierr = VecGhostRestoreLocalForm (const_cast<PetscVector<T>*>(vec2_petsc)->vec(),&v2_loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
 	  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
-	  CHKERRABORT(libMesh::COMM_WORLD,ierr);
+	  LIBMESH_CHKERRABORT(ierr);
     }
 
 #endif
@@ -1267,9 +1267,9 @@ void PetscVector<T>::print_matlab (const std::string name) const
   PetscViewer petsc_viewer;
 
 
-  ierr = PetscViewerCreate (libMesh::COMM_WORLD,
+  ierr = PetscViewerCreate (this->communicator().get(),
 			    &petsc_viewer);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 
   /**
    * Create an ASCII file containing the matrix
@@ -1277,17 +1277,17 @@ void PetscVector<T>::print_matlab (const std::string name) const
    */
   if (name != "NULL")
     {
-      ierr = PetscViewerASCIIOpen( libMesh::COMM_WORLD,
+      ierr = PetscViewerASCIIOpen( this->communicator().get(),
 				   name.c_str(),
 				   &petsc_viewer);
-             CHKERRABORT(libMesh::COMM_WORLD,ierr);
+             LIBMESH_CHKERRABORT(ierr);
 
       ierr = PetscViewerSetFormat (petsc_viewer,
 				   PETSC_VIEWER_ASCII_MATLAB);
-             CHKERRABORT(libMesh::COMM_WORLD,ierr);
+             LIBMESH_CHKERRABORT(ierr);
 
       ierr = VecView (_vec, petsc_viewer);
-             CHKERRABORT(libMesh::COMM_WORLD,ierr);
+             LIBMESH_CHKERRABORT(ierr);
     }
 
   /**
@@ -1297,10 +1297,10 @@ void PetscVector<T>::print_matlab (const std::string name) const
     {
       ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
 				   PETSC_VIEWER_ASCII_MATLAB);
-             CHKERRABORT(libMesh::COMM_WORLD,ierr);
+             LIBMESH_CHKERRABORT(ierr);
 
       ierr = VecView (_vec, PETSC_VIEWER_STDOUT_WORLD);
-             CHKERRABORT(libMesh::COMM_WORLD,ierr);
+             LIBMESH_CHKERRABORT(ierr);
     }
 
 
@@ -1308,7 +1308,7 @@ void PetscVector<T>::print_matlab (const std::string name) const
    * Destroy the viewer.
    */
   ierr = LibMeshPetscViewerDestroy (&petsc_viewer);
-         CHKERRABORT(libMesh::COMM_WORLD,ierr);
+         LIBMESH_CHKERRABORT(ierr);
 }
 
 
@@ -1340,12 +1340,12 @@ void PetscVector<T>::create_subvector(NumericVector<T>& subvector,
       // entries) is not currently offered by the PetscVector
       // class.  Should we differentiate here between sequential and
       // parallel vector creation based on libMesh::n_processors() ?
-      ierr = VecCreateMPI(libMesh::COMM_WORLD,
+      ierr = VecCreateMPI(this->communicator().get(),
 			  PETSC_DECIDE,          // n_local
 			  rows.size(),           // n_global
-			  &(petsc_subvector->_vec)); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+			  &(petsc_subvector->_vec)); LIBMESH_CHKERRABORT(ierr);
 
-      ierr = VecSetFromOptions (petsc_subvector->_vec); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+      ierr = VecSetFromOptions (petsc_subvector->_vec); LIBMESH_CHKERRABORT(ierr);
 
       // Mark the subvector as initialized
       petsc_subvector->_is_initialized = true;
@@ -1360,24 +1360,24 @@ void PetscVector<T>::create_subvector(NumericVector<T>& subvector,
   Utility::iota (idx.begin(), idx.end(), 0);
 
   // Construct index sets
-  ierr = ISCreateLibMesh(libMesh::COMM_WORLD,
+  ierr = ISCreateLibMesh(this->communicator().get(),
 			 rows.size(),
 			 (PetscInt*) &rows[0],
 			 PETSC_USE_POINTER,
-			 &parent_is); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+			 &parent_is); LIBMESH_CHKERRABORT(ierr);
 
-  ierr = ISCreateLibMesh(libMesh::COMM_WORLD,
+  ierr = ISCreateLibMesh(this->communicator().get(),
 			 rows.size(),
 			 (PetscInt*) &idx[0],
 			 PETSC_USE_POINTER,
-			 &subvector_is); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+			 &subvector_is); LIBMESH_CHKERRABORT(ierr);
 
   // Construct the scatter object
   ierr = VecScatterCreate(this->_vec,
 			  parent_is,
 			  petsc_subvector->_vec,
 			  subvector_is,
-			  &scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+			  &scatter); LIBMESH_CHKERRABORT(ierr);
 
   // Actually perform the scatter
 #if PETSC_VERSION_LESS_THAN(2,3,3)
@@ -1385,32 +1385,32 @@ void PetscVector<T>::create_subvector(NumericVector<T>& subvector,
 			 petsc_subvector->_vec,
 			 INSERT_VALUES,
 			 SCATTER_FORWARD,
-			 scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+			 scatter); LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterEnd(this->_vec,
 		       petsc_subvector->_vec,
 		       INSERT_VALUES,
 		       SCATTER_FORWARD,
-		       scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+		       scatter); LIBMESH_CHKERRABORT(ierr);
 #else
   // API argument order change in PETSc 2.3.3
   ierr = VecScatterBegin(scatter,
 			 this->_vec,
 			 petsc_subvector->_vec,
 			 INSERT_VALUES,
-			 SCATTER_FORWARD); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+			 SCATTER_FORWARD); LIBMESH_CHKERRABORT(ierr);
 
   ierr = VecScatterEnd(scatter,
 		       this->_vec,
 		       petsc_subvector->_vec,
 		       INSERT_VALUES,
-		       SCATTER_FORWARD); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+		       SCATTER_FORWARD); LIBMESH_CHKERRABORT(ierr);
 #endif
 
   // Clean up
-  ierr = LibMeshISDestroy(&parent_is);       CHKERRABORT(libMesh::COMM_WORLD,ierr);
-  ierr = LibMeshISDestroy(&subvector_is);    CHKERRABORT(libMesh::COMM_WORLD,ierr);
-  ierr = LibMeshVecScatterDestroy(&scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);
+  ierr = LibMeshISDestroy(&parent_is);       LIBMESH_CHKERRABORT(ierr);
+  ierr = LibMeshISDestroy(&subvector_is);    LIBMESH_CHKERRABORT(ierr);
+  ierr = LibMeshVecScatterDestroy(&scatter); LIBMESH_CHKERRABORT(ierr);
 
 }
 
diff --git a/src/numerics/trilinos_epetra_vector.C b/src/numerics/trilinos_epetra_vector.C
index c8617a172fa..69a14c88659 100644
--- a/src/numerics/trilinos_epetra_vector.C
+++ b/src/numerics/trilinos_epetra_vector.C
@@ -58,7 +58,7 @@ T EpetraVector<T>::sum () const
   for (unsigned int i=0; i<nl; i++)
     sum += values[i];
 
-  CommWorld.sum<T>(sum);
+  this->communicator().sum(sum);
 
   return sum;
 }
@@ -504,7 +504,7 @@ void EpetraVector<T>::localize (const unsigned int first_local_idx,
 
   // Build a parallel vector, initialize it with the local
   // parts of (*this)
-  EpetraVector<T> parallel_vec;
+  EpetraVector<T> parallel_vec(PARALLEL, this->communicator());
 
   parallel_vec.init (my_size, my_local_size, true, PARALLEL);
 
@@ -538,7 +538,7 @@ void EpetraVector<T>::localize (std::vector<T>& v_local) const
   for (unsigned int i=0; i<nl; i++)
     v_local.push_back((*this->_vec)[i]);
 
-  CommWorld.allgather (v_local);
+  this->communicator().allgather (v_local);
 }
 
 
@@ -564,7 +564,7 @@ void EpetraVector<T>::localize_to_one (std::vector<T>&  v_local,
   for (unsigned int i=0; i<nl; i++)
     v_local.push_back((*this->_vec)[i]);
 
-  CommWorld.gather (pid, v_local);
+  this->communicator().gather (pid, v_local);
 }
 
 
diff --git a/src/reduced_basis/rb_construction.C b/src/reduced_basis/rb_construction.C
index 5a6db118e70..3d8fa23a52b 100644
--- a/src/reduced_basis/rb_construction.C
+++ b/src/reduced_basis/rb_construction.C
@@ -510,7 +510,7 @@ void RBConstruction::allocate_data_structures()
   for(unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
   {
     // Initialize the memory for the vectors
-    Fq_vector[q] = NumericVector<Number>::build().release();
+    Fq_vector[q] = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
     Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
   }
 
@@ -521,7 +521,7 @@ void RBConstruction::allocate_data_structures()
     for(unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
     {
       // Initialize the memory for the vectors
-      non_dirichlet_Fq_vector[q] = NumericVector<Number>::build().release();
+      non_dirichlet_Fq_vector[q] = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
       non_dirichlet_Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
     }
   }
@@ -530,7 +530,7 @@ void RBConstruction::allocate_data_structures()
     for(unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
     {
       // Initialize the memory for the truth output vectors
-      outputs_vector[n][q_l] = (NumericVector<Number>::build().release());
+      outputs_vector[n][q_l] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
       outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
     }
 
@@ -543,7 +543,7 @@ void RBConstruction::allocate_data_structures()
       for(unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
       {
         // Initialize the memory for the truth output vectors
-        non_dirichlet_outputs_vector[n][q_l] = (NumericVector<Number>::build().release());
+        non_dirichlet_outputs_vector[n][q_l] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
         non_dirichlet_outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
       }
     }
@@ -764,7 +764,7 @@ void RBConstruction::truth_assembly()
       matrix->add(get_rb_theta_expansion().eval_A_theta(q_a, mu), *get_Aq(q_a));
     }
 
-    AutoPtr< NumericVector<Number> > temp_vec = NumericVector<Number>::build();
+    AutoPtr< NumericVector<Number> > temp_vec = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
     temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
     for(unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
     {
@@ -1302,12 +1302,12 @@ void RBConstruction::enrich_RB_space()
 {
   START_LOG("enrich_RB_space()", "RBConstruction");
 
-  NumericVector<Number>* new_bf = NumericVector<Number>::build().release();
+  NumericVector<Number>* new_bf = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
   new_bf->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
   *new_bf = *solution;
 
   // compute orthogonalization
-  AutoPtr< NumericVector<Number> > proj_index = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > proj_index = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   proj_index->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
   if(single_matrix_mode)
@@ -1495,7 +1495,7 @@ void RBConstruction::update_RB_system_matrices()
 
   unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
 
-  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   temp->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
   for(unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
@@ -1617,7 +1617,7 @@ void RBConstruction::update_residual_terms(bool compute_inner_products)
       // Initialize the vector in which we'll store the representor
       if(!get_rb_evaluation().Aq_representor[q_a][i])
       {
-        get_rb_evaluation().Aq_representor[q_a][i] = (NumericVector<Number>::build().release());
+        get_rb_evaluation().Aq_representor[q_a][i] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
         get_rb_evaluation().Aq_representor[q_a][i]->init(this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
       }
 
@@ -1802,7 +1802,7 @@ void RBConstruction::compute_output_dual_innerprods()
     std::vector< NumericVector<Number>* > L_q_representor(max_Q_l);
     for(unsigned int q=0; q<max_Q_l; q++)
     {
-      L_q_representor[q] = (NumericVector<Number>::build().release());
+      L_q_representor[q] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
       L_q_representor[q]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
     }
 
@@ -1963,7 +1963,7 @@ void RBConstruction::compute_Fq_representor_innerprods(bool compute_inner_produc
     {
       if(!Fq_representor[q_f])
       {
-        Fq_representor[q_f] = (NumericVector<Number>::build().release());
+        Fq_representor[q_f] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
         Fq_representor[q_f]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
       }
 
@@ -2434,7 +2434,7 @@ void RBConstruction::read_riesz_representors_from_files(const std::string& riesz
 
     read_serialized_data(fqr_data, false);
 
-    Fq_representor[i] = NumericVector<Number>::build().release();
+    Fq_representor[i] = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
     Fq_representor[i]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
     // No need to copy, just swap
@@ -2487,7 +2487,7 @@ void RBConstruction::read_riesz_representors_from_files(const std::string& riesz
 
       read_serialized_data(aqr_data, false);
 
-      get_rb_evaluation().Aq_representor[i][j] = NumericVector<Number>::build().release();
+      get_rb_evaluation().Aq_representor[i][j] = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
       get_rb_evaluation().Aq_representor[i][j]->init (n_dofs(), n_local_dofs(),
                                             false, libMeshEnums::PARALLEL);
 
diff --git a/src/reduced_basis/rb_construction_base.C b/src/reduced_basis/rb_construction_base.C
index 292273dacc8..fe03a4f113d 100644
--- a/src/reduced_basis/rb_construction_base.C
+++ b/src/reduced_basis/rb_construction_base.C
@@ -92,7 +92,7 @@ void RBConstructionBase<Base>::init_data ()
 
   // Initialize the inner product storage vector, which is useful for
   // storing intermediate results when evaluating inner products
-  inner_product_storage_vector = NumericVector<Number>::build();
+  inner_product_storage_vector = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   inner_product_storage_vector->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 }
 
@@ -297,7 +297,7 @@ void RBConstructionBase<Base>::load_training_set(std::map< std::string, std::vec
   it = training_parameters.begin();
   for( ; it != it_end; ++it)
   {
-    it->second = NumericVector<Number>::build().release();
+    it->second = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
     it->second->init(n_global_training_samples, n_local_training_samples, false, libMeshEnums::PARALLEL);
   }
 
@@ -390,7 +390,7 @@ void RBConstructionBase<Base>::generate_training_parameters_random(const Paralle
     for( ; it != it_end; ++it)
     {
       std::string param_name = it->first;
-      training_parameters_in[param_name] = NumericVector<Number>::build().release();
+      training_parameters_in[param_name] = NumericVector<Number>::build(libMesh::default_solver_package(), communicator).release();
 
       if(!serial_training_set)
       {
@@ -523,7 +523,7 @@ void RBConstructionBase<Base>::generate_training_parameters_partially_random(con
     for( ; it != it_end; ++it)
     {
       std::string param_name = it->first;
-      training_parameters_in[param_name] = NumericVector<Number>::build().release();
+      training_parameters_in[param_name] = NumericVector<Number>::build(libMesh::default_solver_package(), communicator).release();
 
       if(!serial_training_set)
       {
@@ -679,7 +679,7 @@ void RBConstructionBase<Base>::generate_training_parameters_deterministic(const
     for( ; it != it_end; ++it)
     {
       std::string param_name = it->first;
-      training_parameters_in[param_name] = NumericVector<Number>::build().release();
+      training_parameters_in[param_name] = NumericVector<Number>::build(libMesh::default_solver_package(), communicator).release();
 
       if(!serial_training_set)
       {
diff --git a/src/reduced_basis/rb_eim_assembly.C b/src/reduced_basis/rb_eim_assembly.C
index 19735bca301..666e817a98a 100644
--- a/src/reduced_basis/rb_eim_assembly.C
+++ b/src/reduced_basis/rb_eim_assembly.C
@@ -36,7 +36,7 @@ RBEIMAssembly::RBEIMAssembly(RBEIMConstruction& rb_eim_con_in,
                              unsigned int basis_function_index_in)
 : _rb_eim_con(rb_eim_con_in),
   _basis_function_index(basis_function_index_in),
-  _ghosted_basis_function(NumericVector<Number>::build())
+  _ghosted_basis_function(NumericVector<Number>::build(libMesh::default_solver_package(), rb_eim_con_in.system().communicator()))
 {
   // localize the vector that stores the basis function for this assembly object,
   // i.e. the vector that is used in evaluate_basis_function_at_quad_pts
diff --git a/src/reduced_basis/rb_eim_construction.C b/src/reduced_basis/rb_eim_construction.C
index ffbeaa23322..f0fbdafd401 100644
--- a/src/reduced_basis/rb_eim_construction.C
+++ b/src/reduced_basis/rb_eim_construction.C
@@ -166,7 +166,7 @@ void RBEIMConstruction::initialize_rb_construction()
   Parent::initialize_rb_construction();
 
   // initialize a serial vector that we will use for MeshFunction evaluations
-  _ghosted_meshfunction_vector = NumericVector<Number>::build();
+  _ghosted_meshfunction_vector = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   _ghosted_meshfunction_vector->init (this->n_dofs(), this->n_local_dofs(),
                                       this->get_dof_map().get_send_list(), false,
                                       GHOSTED);
@@ -350,7 +350,7 @@ void RBEIMConstruction::enrich_RB_space()
     eim_eval.interpolation_points.push_back(optimal_point);
     eim_eval.interpolation_points_var.push_back(optimal_var);
 
-    NumericVector<Number>* new_bf = NumericVector<Number>::build().release();
+    NumericVector<Number>* new_bf = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
     new_bf->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
     *new_bf = *solution;
     get_rb_evaluation().basis_functions.push_back( new_bf );
diff --git a/src/reduced_basis/rb_evaluation.C b/src/reduced_basis/rb_evaluation.C
index dd291cdbe71..65e61c262cf 100644
--- a/src/reduced_basis/rb_evaluation.C
+++ b/src/reduced_basis/rb_evaluation.C
@@ -1064,7 +1064,7 @@ void RBEvaluation::read_in_vectors(System& sys,
 
 	  sys.read_serialized_data(vector_data, false);
 
-	  vectors[i] = NumericVector<Number>::build().release();
+	  vectors[i] = NumericVector<Number>::build(libMesh::default_solver_package(), sys.communicator()).release();
 	  vectors[i]->init (sys.n_dofs(), sys.n_local_dofs(), false, libMeshEnums::PARALLEL);
 
 	  // No need to copy, just swap
@@ -1080,7 +1080,7 @@ void RBEvaluation::read_in_vectors(System& sys,
       // Allocate storage for each vector
       for(unsigned int i=0; i<vectors.size(); i++)
 	{
-	  vectors[i] = NumericVector<Number>::build().release();
+	  vectors[i] = NumericVector<Number>::build(libMesh::default_solver_package(), sys.communicator()).release();
 	  vectors[i]->init (sys.n_dofs(), sys.n_local_dofs(), false, libMeshEnums::PARALLEL);
 	}
 
diff --git a/src/reduced_basis/transient_rb_construction.C b/src/reduced_basis/transient_rb_construction.C
index de70192275c..7e752caa50f 100644
--- a/src/reduced_basis/transient_rb_construction.C
+++ b/src/reduced_basis/transient_rb_construction.C
@@ -261,7 +261,7 @@ void TransientRBConstruction::allocate_data_structures()
 
   for(unsigned int i=0; i<n_time_levels; i++)
   {
-    temporal_data[i] = (NumericVector<Number>::build().release());
+    temporal_data[i] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
     temporal_data[i]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
   }
 
@@ -293,7 +293,7 @@ void TransientRBConstruction::assemble_affine_expansion()
     // Load the initial condition into the solution vector
     initialize_truth();
 
-    AutoPtr< NumericVector<Number> > temp1 = NumericVector<Number>::build();
+    AutoPtr< NumericVector<Number> > temp1 = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
     temp1->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
     // First compute the right-hand side vector for the L2 projection
@@ -435,7 +435,7 @@ void TransientRBConstruction::mass_matrix_scaled_matvec(Number scalar,
 
   const unsigned int Q_m = trans_theta_expansion.get_n_M_terms();
 
-  AutoPtr< NumericVector<Number> > temp_vec = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > temp_vec = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
   for(unsigned int q=0; q<Q_m; q++)
@@ -490,7 +490,7 @@ void TransientRBConstruction::truth_assembly()
     add_scaled_mass_matrix(1./dt, matrix);
     mass_matrix_scaled_matvec(1./dt, *rhs, *current_local_solution);
 
-    AutoPtr< NumericVector<Number> > temp_vec = NumericVector<Number>::build();
+    AutoPtr< NumericVector<Number> > temp_vec = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
     temp_vec->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
     for(unsigned int q_a=0; q_a<Q_a; q_a++)
@@ -914,7 +914,7 @@ Number TransientRBConstruction::set_error_temporal_data()
   {
     unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
 
-    AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build();
+    AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
     temp->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
     // First compute the right-hand side vector for the projection
@@ -1031,7 +1031,7 @@ void TransientRBConstruction::add_IC_to_RB_space()
   initialize_truth();
 
   // load the new basis function into the basis_functions vector.
-  get_rb_evaluation().basis_functions.push_back( NumericVector<Number>::build().release() );
+  get_rb_evaluation().basis_functions.push_back( NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release() );
   NumericVector<Number>& current_bf = get_rb_evaluation().get_basis_function(get_rb_evaluation().get_n_basis_functions()-1);
   current_bf.init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
   current_bf = *solution;
@@ -1173,7 +1173,7 @@ void TransientRBConstruction::enrich_RB_space()
   while (true)
     {
       // load the new basis function into the basis_functions vector.
-      get_rb_evaluation().basis_functions.push_back( NumericVector<Number>::build().release() );
+      get_rb_evaluation().basis_functions.push_back( NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release() );
       NumericVector<Number>& current_bf = get_rb_evaluation().get_basis_function(get_rb_evaluation().get_n_basis_functions()-1);
       current_bf.init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
       current_bf.zero();
@@ -1315,7 +1315,7 @@ void TransientRBConstruction::update_RB_system_matrices()
 
   unsigned int RB_size = get_rb_evaluation().get_n_basis_functions();
 
-  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   temp->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
   for(unsigned int i=(RB_size-delta_N); i<RB_size; i++)
@@ -1427,7 +1427,7 @@ void TransientRBConstruction::update_residual_terms(bool compute_inner_products)
       // Initialize the vectors when we need them
       if(!trans_rb_eval.M_q_representor[q_m][i])
       {
-        trans_rb_eval.M_q_representor[q_m][i] = (NumericVector<Number>::build().release());
+        trans_rb_eval.M_q_representor[q_m][i] = (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release());
         trans_rb_eval.M_q_representor[q_m][i]->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
       }
 
@@ -1618,10 +1618,10 @@ void TransientRBConstruction::update_RB_initial_condition_all_N()
   // Load the initial condition into the solution vector
   initialize_truth();
 
-  AutoPtr< NumericVector<Number> > temp1 = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > temp1 = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   temp1->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
-  AutoPtr< NumericVector<Number> > temp2 = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > temp2 = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   temp2->init (this->n_dofs(), this->n_local_dofs(), false, libMeshEnums::PARALLEL);
 
 
@@ -1892,7 +1892,7 @@ void TransientRBConstruction::read_riesz_representors_from_files(const std::stri
 
 	read_serialized_data(aqr_data, false);
 
-	trans_rb_eval.M_q_representor[i][j] = NumericVector<Number>::build().release();
+	trans_rb_eval.M_q_representor[i][j] = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
 	trans_rb_eval.M_q_representor[i][j]->init (n_dofs(), n_local_dofs(),
 	                 	     false, libMeshEnums::PARALLEL);
 
diff --git a/src/solution_transfer/meshfunction_solution_transfer.C b/src/solution_transfer/meshfunction_solution_transfer.C
index cac810bb154..3fd7a2068a0 100644
--- a/src/solution_transfer/meshfunction_solution_transfer.C
+++ b/src/solution_transfer/meshfunction_solution_transfer.C
@@ -50,7 +50,7 @@ MeshFunctionSolutionTransfer::transfer(const Variable & from_var, const Variable
   EquationSystems & from_es = from_sys->get_equation_systems();
 
   //Create a serialized version of the solution vector
-  NumericVector<Number> * serialized_solution = NumericVector<Number>::build().release();
+  NumericVector<Number> * serialized_solution = NumericVector<Number>::build(libMesh::default_solver_package(), from_sys->get_mesh().communicator()).release();
   serialized_solution->init(from_sys->n_dofs(), false, SERIAL);
 
   // Need to pull down a full copy of this vector on every processor so we can get values in parallel
diff --git a/src/solvers/petsc_diff_solver.C b/src/solvers/petsc_diff_solver.C
index 1eadf869b93..9542880444d 100644
--- a/src/solvers/petsc_diff_solver.C
+++ b/src/solvers/petsc_diff_solver.C
@@ -80,7 +80,7 @@ __libmesh_petsc_diff_solver_residual (SNES, Vec x, Vec r, void *ctx)
     *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
   PetscVector<Number>& R_system =
     *libmesh_cast_ptr<PetscVector<Number>*>(sys.rhs);
-  PetscVector<Number> X_input(x), R_input(r);
+  PetscVector<Number> X_input(x, sys.communicator()), R_input(r, sys.communicator());
 
   // DiffSystem assembles from the solution and into the rhs, so swap
   // those with our input vectors before assembling.  They'll probably
@@ -126,7 +126,7 @@ __libmesh_petsc_diff_solver_jacobian (SNES, Vec x, Mat *libmesh_dbg_var(j), Mat
 
   PetscVector<Number>& X_system =
     *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
-  PetscVector<Number> X_input(x);
+  PetscVector<Number> X_input(x, sys.communicator());
 
   PetscMatrix<Number> J_input(*pc);
   PetscMatrix<Number>& J_system =
diff --git a/src/solvers/petsc_linear_solver.C b/src/solvers/petsc_linear_solver.C
index 9f15623c3e7..b3f427f6419 100644
--- a/src/solvers/petsc_linear_solver.C
+++ b/src/solvers/petsc_linear_solver.C
@@ -63,8 +63,8 @@ extern "C"
   {
     Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
 
-    PetscVector<Number> x_vec(x);
-    PetscVector<Number> y_vec(y);
+    PetscVector<Number> x_vec(x, preconditioner->communicator());
+    PetscVector<Number> y_vec(y, preconditioner->communicator());
 
     preconditioner->apply(x_vec,y_vec);
 
@@ -94,8 +94,8 @@ extern "C"
     PetscErrorCode ierr = PCShellGetContext(pc,&ctx);CHKERRQ(ierr);
     Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
 
-    PetscVector<Number> x_vec(x);
-    PetscVector<Number> y_vec(y);
+    PetscVector<Number> x_vec(x/*, preconditioner->communicator()*/);
+    PetscVector<Number> y_vec(y/*, preconditioner->communicator()*/);
 
     preconditioner->apply(x_vec,y_vec);
 
@@ -1823,8 +1823,8 @@ PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult(Mat mat, Vec arg,
   const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
 
   /* Make \p NumericVector instances around the vectors.  */
-  PetscVector<T> arg_global(arg);
-  PetscVector<T> dest_global(dest);
+  PetscVector<T> arg_global(arg/*, shell_matrix.communicator()*/);
+  PetscVector<T> dest_global(dest/*, shell_matrix.communicator()*/);
 
   /* Call the user function.  */
   shell_matrix.vector_mult(dest_global,arg_global);
@@ -1847,9 +1847,9 @@ PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult_add(Mat mat, Vec a
   const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
 
   /* Make \p NumericVector instances around the vectors.  */
-  PetscVector<T> arg_global(arg);
-  PetscVector<T> dest_global(dest);
-  PetscVector<T> add_global(add);
+  PetscVector<T> arg_global(arg/*, shell_matrix.communicator()*/);
+  PetscVector<T> dest_global(dest/*, shell_matrix.communicator()*/);
+  PetscVector<T> add_global(add/*, shell_matrix.communicator()*/);
 
   if(add!=arg)
     {
@@ -1877,7 +1877,7 @@ PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_get_diagonal(Mat mat, V
   const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
 
   /* Make \p NumericVector instances around the vector.  */
-  PetscVector<T> dest_global(dest);
+  PetscVector<T> dest_global(dest/*, shell_matrix.communicator()*/);
 
   /* Call the user function.  */
   shell_matrix.get_diagonal(dest_global);
diff --git a/src/solvers/petsc_nonlinear_solver.C b/src/solvers/petsc_nonlinear_solver.C
index 4379208c67f..d8fb5931843 100644
--- a/src/solvers/petsc_nonlinear_solver.C
+++ b/src/solvers/petsc_nonlinear_solver.C
@@ -102,7 +102,7 @@ extern "C"
 
     PetscVector<Number>& X_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
     PetscVector<Number>& R_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.rhs);
-    PetscVector<Number> X_global(x), R(r);
+    PetscVector<Number> X_global(x, sys.communicator()), R(r, sys.communicator());
 
     // Use the systems update() to get a good local version of the parallel solution
     X_global.swap(X_sys);
@@ -180,7 +180,7 @@ extern "C"
     PetscMatrix<Number> Jac(*jac);
     PetscVector<Number>& X_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
     PetscMatrix<Number>& Jac_sys = *libmesh_cast_ptr<PetscMatrix<Number>*>(sys.matrix);
-    PetscVector<Number> X_global(x);
+    PetscVector<Number> X_global(x, sys.communicator());
 
     // Set the dof maps
     PC.attach_dof_map(sys.get_dof_map());
diff --git a/src/solvers/petscdmlibmesh.C b/src/solvers/petscdmlibmesh.C
index 7ad92c47487..93cd886a648 100644
--- a/src/solvers/petscdmlibmesh.C
+++ b/src/solvers/petscdmlibmesh.C
@@ -720,7 +720,7 @@ static PetscErrorCode DMlibMeshFunction(DM dm, Vec x, Vec r)
   NonlinearImplicitSystem& sys = *_sys;
   PetscVector<Number>& X_sys = *libmesh_cast_ptr<PetscVector<Number>* >(sys.solution.get());
   PetscVector<Number>& R_sys = *libmesh_cast_ptr<PetscVector<Number>* >(sys.rhs);
-  PetscVector<Number> X_global(x), R(r);
+  PetscVector<Number> X_global(x, _sys->communicator()), R(r, _sys->communicator());
 
   // Use the systems update() to get a good local version of the parallel solution
   X_global.swap(X_sys);
@@ -796,7 +796,7 @@ static PetscErrorCode DMlibMeshJacobian(DM dm, Vec x, Mat jac, Mat pc, MatStruct
   PetscMatrix<Number>  Jac(jac);
   PetscVector<Number>& X_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
   PetscMatrix<Number>& Jac_sys = *libmesh_cast_ptr<PetscMatrix<Number>*>(sys.matrix);
-  PetscVector<Number>  X_global(x);
+  PetscVector<Number>  X_global(x, sys.communicator());
 
   // Set the dof maps
   the_pc.attach_dof_map(sys.get_dof_map());
@@ -872,8 +872,8 @@ static PetscErrorCode DMVariableBounds_libMesh(DM dm, Vec xl, Vec xu)
   NonlinearImplicitSystem* _sys;
   ierr = DMLibMeshGetSystem(dm, _sys); CHKERRQ(ierr);
   NonlinearImplicitSystem& sys = *_sys;
-  PetscVector<Number> XL(xl);
-  PetscVector<Number> XU(xu);
+  PetscVector<Number> XL(xl, sys.communicator());
+  PetscVector<Number> XU(xu, sys.communicator());
   PetscFunctionBegin;
 
   ierr = VecSet(xl, SNES_VI_NINF); CHKERRQ(ierr);
diff --git a/src/solvers/slepc_eigen_solver.C b/src/solvers/slepc_eigen_solver.C
index 09eec8062ad..54a0cd1569a 100644
--- a/src/solvers/slepc_eigen_solver.C
+++ b/src/solvers/slepc_eigen_solver.C
@@ -746,8 +746,8 @@ PetscErrorCode SlepcEigenSolver<T>::_petsc_shell_matrix_mult(Mat mat, Vec arg, V
   const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
 
   /* Make \p NumericVector instances around the vectors.  */
-  PetscVector<T> arg_global(arg);
-  PetscVector<T> dest_global(dest);
+  PetscVector<T> arg_global(arg, shell_matrix.communicator());
+  PetscVector<T> dest_global(dest, shell_matrix.communicator());
 
   /* Call the user function.  */
   shell_matrix.vector_mult(dest_global,arg_global);
@@ -768,7 +768,7 @@ PetscErrorCode SlepcEigenSolver<T>::_petsc_shell_matrix_get_diagonal(Mat mat, Ve
   const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
 
   /* Make \p NumericVector instances around the vector.  */
-  PetscVector<T> dest_global(dest);
+  PetscVector<T> dest_global(dest, shell_matrix.communicator());
 
   /* Call the user function.  */
   shell_matrix.get_diagonal(dest_global);
diff --git a/src/solvers/trilinos_nox_nonlinear_solver.C b/src/solvers/trilinos_nox_nonlinear_solver.C
index 6822b9db22f..df94face0a8 100644
--- a/src/solvers/trilinos_nox_nonlinear_solver.C
+++ b/src/solvers/trilinos_nox_nonlinear_solver.C
@@ -92,7 +92,7 @@ bool Problem_Interface::computeF(const Epetra_Vector& x, Epetra_Vector& r,
 
   NonlinearImplicitSystem &sys = _solver->system();
 
-  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x)), R(r);
+  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.communicator()), R(r, sys.communicator());
   EpetraVector<Number>& X_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.solution.get());
   EpetraVector<Number>& R_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.rhs);
 
@@ -149,7 +149,7 @@ bool Problem_Interface::computeJacobian(const Epetra_Vector & x,
 
   EpetraMatrix<Number> Jac(&dynamic_cast<Epetra_FECrsMatrix &>(jac));
   EpetraVector<Number>& X_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.solution.get());
-  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x));
+  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.communicator());
 
   // Set the dof maps
   Jac.attach_dof_map(sys.get_dof_map());
@@ -209,7 +209,7 @@ bool Problem_Interface::computePreconditioner(const Epetra_Vector & x,
 
   EpetraMatrix<Number> Jac(dynamic_cast<Epetra_FECrsMatrix *>(tpc.mat()));
   EpetraVector<Number>& X_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.solution.get());
-  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x));
+  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.communicator());
 
   // Set the dof maps
   Jac.attach_dof_map(sys.get_dof_map());
diff --git a/src/solvers/unsteady_solver.C b/src/solvers/unsteady_solver.C
index 5da8b42c4e5..67756a8511f 100644
--- a/src/solvers/unsteady_solver.C
+++ b/src/solvers/unsteady_solver.C
@@ -29,7 +29,7 @@ namespace libMesh
 
 UnsteadySolver::UnsteadySolver (sys_type& s)
   : TimeSolver(s),
-    old_local_nonlinear_solution (NumericVector<Number>::build()),
+    old_local_nonlinear_solution (NumericVector<Number>::build(libMesh::default_solver_package(), s.communicator())),
     first_solve                  (true),
     first_adjoint_step (true)
 {
diff --git a/src/systems/condensed_eigen_system.C b/src/systems/condensed_eigen_system.C
index 02f4d965d7f..87e0e5ab080 100644
--- a/src/systems/condensed_eigen_system.C
+++ b/src/systems/condensed_eigen_system.C
@@ -188,7 +188,7 @@ std::pair<Real, Real> CondensedEigenSystem::get_eigenpair(unsigned int i)
 
   // This function assumes that condensed_solve has just been called.
   // If this is not the case, then we will trip an asset in get_eigenpair
-  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build();
+  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   unsigned int n_local = local_non_condensed_dofs_vector.size();
   unsigned int n       = n_local;
   this->communicator().sum(n);
diff --git a/src/systems/system.C b/src/systems/system.C
index 59586e3601d..048df622048 100644
--- a/src/systems/system.C
+++ b/src/systems/system.C
@@ -60,8 +60,8 @@ System::System (EquationSystems& es,
   assemble_before_solve             (true),
   use_fixed_solution                (false),
   extra_quadrature_order            (0),
-  solution                          (NumericVector<Number>::build()),
-  current_local_solution            (NumericVector<Number>::build()),
+  solution                          (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator())),
+  current_local_solution            (NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator())),
   time                              (0.),
   qoi                               (0),
   _init_system_function             (NULL),
@@ -681,7 +681,7 @@ NumericVector<Number> & System::add_vector (const std::string& vec_name,
     return *(_vectors[vec_name]);
 
   // Otherwise build the vector
-  NumericVector<Number>* buf = NumericVector<Number>::build().release();
+  NumericVector<Number>* buf = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator()).release();
   _vectors.insert (std::make_pair (vec_name, buf));
   _vector_projections.insert (std::make_pair (vec_name, projections));
 
@@ -1444,7 +1444,7 @@ Real System::calculate_norm(const NumericVector<Number>& v,
     }
 
   // Localize the potentially parallel vector
-  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
+  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
   local_v->init(v.size(), true, SERIAL);
   v.localize (*local_v, _dof_map->get_send_list());
 
diff --git a/src/systems/system_projection.C b/src/systems/system_projection.C
index 125712e6777..49d1c17835f 100644
--- a/src/systems/system_projection.C
+++ b/src/systems/system_projection.C
@@ -317,8 +317,8 @@ void System::project_vector (const NumericVector<Number>& old_v,
       projection_list.unique();
 
       new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
-      new_vector_built = NumericVector<Number>::build();
-      local_old_vector_built = NumericVector<Number>::build();
+      new_vector_built = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
+      local_old_vector_built = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
       new_vector_ptr = new_vector_built.get();
       local_old_vector = local_old_vector_built.get();
       new_vector_ptr->init(this->n_dofs(), false, SERIAL);
@@ -340,7 +340,7 @@ void System::project_vector (const NumericVector<Number>& old_v,
       new_v.init (this->n_dofs(), this->n_local_dofs(),
                   this->get_dof_map().get_send_list(), false, GHOSTED);
 
-      local_old_vector_built = NumericVector<Number>::build();
+      local_old_vector_built = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
       new_vector_ptr = &new_v;
       local_old_vector = local_old_vector_built.get();
       local_old_vector->init(old_v.size(), old_v.local_size(),
@@ -403,7 +403,7 @@ void System::project_vector (const NumericVector<Number>& old_v,
   // creating a temporary parallel vector to use localize! - RHS
   if (old_v.type() == SERIAL)
     {
-      AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build();
+      AutoPtr<NumericVector<Number> > dist_v = NumericVector<Number>::build(libMesh::default_solver_package(), this->communicator());
       dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
       dist_v->close();