diff --git a/packages/panzer/adapters-ioss/test/ioss_connmngr/tIOSSConnManager.cpp b/packages/panzer/adapters-ioss/test/ioss_connmngr/tIOSSConnManager.cpp
index 1e009a17a076..dbd838a91860 100644
--- a/packages/panzer/adapters-ioss/test/ioss_connmngr/tIOSSConnManager.cpp
+++ b/packages/panzer/adapters-ioss/test/ioss_connmngr/tIOSSConnManager.cpp
@@ -77,12 +77,6 @@
 #include "Ioss_IOFactory.h"
 #include "Ioss_DatabaseIO.h"
 
-//#ifdef HAVE_MPI
-//   #include "Epetra_MpiComm.h"
-//#else
-//   #include "Epetra_SerialComm.h"
-//#endif
-
 using Teuchos::RCP;
 using Teuchos::rcp;
 
diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/gather_orientation.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/gather_orientation.cpp
index 1b5cb7be4fa2..150381a87a44 100644
--- a/packages/panzer/adapters-stk/test/evaluator_tests/gather_orientation.cpp
+++ b/packages/panzer/adapters-stk/test/evaluator_tests/gather_orientation.cpp
@@ -68,8 +68,6 @@ using Teuchos::rcp;
 #include "Teuchos_DefaultMpiComm.hpp"
 #include "Teuchos_OpaqueWrapper.hpp"
 
-#include "Epetra_MpiComm.h"
-
 #include "user_app_EquationSetFactory.hpp"
 
 #include <cstdio> // for get char
diff --git a/packages/panzer/adapters-stk/test/stk_interface_test/tPointLocationSearch.cpp b/packages/panzer/adapters-stk/test/stk_interface_test/tPointLocationSearch.cpp
index eff6c85dbf76..f143fc08bcd4 100644
--- a/packages/panzer/adapters-stk/test/stk_interface_test/tPointLocationSearch.cpp
+++ b/packages/panzer/adapters-stk/test/stk_interface_test/tPointLocationSearch.cpp
@@ -53,12 +53,6 @@
 
 #include "Shards_BasicTopologies.hpp"
 
-#ifdef HAVE_MPI
-   #include "Epetra_MpiComm.h"
-#else
-   #include "Epetra_SerialComm.h"
-#endif
-
 // STK_search objects
 #include <stk_search/CoarseSearch.hpp>
 #include <stk_search/BoundingBox.hpp>
diff --git a/packages/panzer/adapters-stk/test/stk_interface_test/tUniformRef.cpp b/packages/panzer/adapters-stk/test/stk_interface_test/tUniformRef.cpp
index d8457e4acf20..c49a1eea7100 100644
--- a/packages/panzer/adapters-stk/test/stk_interface_test/tUniformRef.cpp
+++ b/packages/panzer/adapters-stk/test/stk_interface_test/tUniformRef.cpp
@@ -55,12 +55,6 @@
 
 #include "Shards_BasicTopologies.hpp"
 
-#ifdef HAVE_MPI
-   #include "Epetra_MpiComm.h"
-#else
-   #include "Epetra_SerialComm.h"
-#endif
-
 /*
 #include "stk_mesh/base/GetEntities.hpp"
 #include "stk_mesh/base/Selector.hpp"
diff --git a/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory.hpp b/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory.hpp
index 26068011a557..31c31145dea2 100644
--- a/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory.hpp
+++ b/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory.hpp
@@ -82,6 +82,7 @@ template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename Globa
 class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
                                     , public ThyraObjFactory<double> {
 public:
+   typedef TpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> TLOC;
    typedef BlockedTpetraLinearObjContainer<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> BTLOC;
    typedef Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> VectorType;
    typedef Tpetra::CrsMatrix<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> CrsMatrixType;
@@ -166,32 +167,62 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
    //! Use preconstructed scatter evaluators
    template <typename EvalT>
    Teuchos::RCP<panzer::CloneableEvaluator> buildScatter() const
-   { return Teuchos::rcp(new ScatterResidual_BlockedTpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(blockedDOFManager_)); }
+   {
+      if(gidProviders_.front()->containsBlockedDOFManager()) {
+         return Teuchos::rcp(new ScatterResidual_BlockedTpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGlobalIndexer()));
+      }
+      return Teuchos::rcp(new ScatterResidual_Tpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getRangeGlobalIndexer()));
+   }
 
    //! Use preconstructed gather evaluators
    template <typename EvalT>
    Teuchos::RCP<panzer::CloneableEvaluator > buildGather() const
-   { return Teuchos::rcp(new GatherSolution_BlockedTpetra<EvalT,Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(blockedDOFManager_)); }
+   {
+      if(gidProviders_.front()->containsBlockedDOFManager()) {
+         return Teuchos::rcp(new GatherSolution_BlockedTpetra<EvalT,Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGlobalIndexer()));
+      }
+      return Teuchos::rcp(new GatherSolution_Tpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getRangeGlobalIndexer()));
+   }
 
    //! Use preconstructed gather evaluators
    template <typename EvalT>
    Teuchos::RCP<panzer::CloneableEvaluator > buildGatherTangent() const
-   { return Teuchos::rcp(new GatherTangent_BlockedTpetra<EvalT,Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(blockedDOFManager_)); }
+   {
+      if(gidProviders_.front()->containsBlockedDOFManager()) {
+         return Teuchos::rcp(new GatherTangent_BlockedTpetra<EvalT,Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGlobalIndexer()));
+      }
+      return Teuchos::rcp(new GatherTangent_Tpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getRangeGlobalIndexer()));
+   }
 
    //! Use preconstructed gather evaluators
    template <typename EvalT>
    Teuchos::RCP<panzer::CloneableEvaluator > buildGatherDomain() const
-   { return Teuchos::rcp(new GatherSolution_BlockedTpetra<EvalT,Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(blockedDOFManager_)); }
+   {
+      if(gidProviders_.back()->containsBlockedDOFManager()) {
+         return Teuchos::rcp(new GatherSolution_BlockedTpetra<EvalT,Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getDomainBlockedIndexer()));
+      }
+      return Teuchos::rcp(new GatherSolution_Tpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getDomainGlobalIndexer()));
+   }
 
    //! Use preconstructed gather evaluators
    template <typename EvalT>
    Teuchos::RCP<panzer::CloneableEvaluator > buildGatherOrientation() const
-   { return Teuchos::rcp(new GatherOrientation<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT>(nc2c_vector(blockedDOFManager_->getFieldDOFManagers()))); }
+   {
+      if(gidProviders_.front()->containsBlockedDOFManager()) {
+         return Teuchos::rcp(new GatherOrientation<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT>(nc2c_vector(getGlobalIndexer()->getFieldDOFManagers())));
+      }
+      return Teuchos::rcp(new GatherOrientation<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT>(getRangeGlobalIndexer()));
+   }
 
    //! Use preconstructed dirichlet scatter evaluators
    template <typename EvalT>
    Teuchos::RCP<panzer::CloneableEvaluator> buildScatterDirichlet() const
-   { return Teuchos::rcp(new ScatterDirichletResidual_BlockedTpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(blockedDOFManager_)); }
+   {
+      if(gidProviders_.front()->containsBlockedDOFManager()) {
+         return Teuchos::rcp(new ScatterDirichletResidual_BlockedTpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGlobalIndexer()));
+      }
+      return Teuchos::rcp(new ScatterDirichletResidual_Tpetra<EvalT,Traits,LocalOrdinalT,GlobalOrdinalT,NodeT>(getRangeGlobalIndexer()));
+   }
 
 /*************** Generic helper functions for container setup *******************/
 
@@ -267,9 +298,11 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
 
    //! get the map from the matrix
    virtual Teuchos::RCP<const MapType> getMap(int i) const;
+   virtual Teuchos::RCP<const MapType> getColMap(int i) const;
 
    //! get the ghosted map from the matrix
    virtual Teuchos::RCP<const MapType> getGhostedMap(int i) const;
+   virtual Teuchos::RCP<const MapType> getColGhostedMap(int i) const;
 
    //! get the graph of the crs matrix
    virtual Teuchos::RCP<const CrsGraphType> getGraph(int i,int j) const;
@@ -279,9 +312,11 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
 
    //! get importer for converting an overalapped object to a "normal" object
    virtual Teuchos::RCP<const ImportType> getGhostedImport(int i) const;
+   virtual Teuchos::RCP<const ImportType> getColGhostedImport(int i) const;
 
    //! get exporter for converting an overalapped object to a "normal" object
    virtual Teuchos::RCP<const ExportType> getGhostedExport(int j) const;
+   virtual Teuchos::RCP<const ExportType> getColGhostedExport(int j) const;
 
    Teuchos::RCP<CrsMatrixType> getTpetraMatrix(int i,int j) const;
    Teuchos::RCP<CrsMatrixType> getGhostedTpetraMatrix(int i,int j) const;
@@ -308,28 +343,134 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
    virtual void endFill(LinearObjContainer & loc) const;
 
    Teuchos::RCP<const panzer::BlockedDOFManager> getGlobalIndexer() const
-   { return blockedDOFManager_; }
+   {
+      TEUCHOS_ASSERT(!gidProviders_.empty());
+      return gidProviders_.front()->getBlockedIndexer();
+   }
 
    //! Get the domain unique global indexer this factory was created with.
    Teuchos::RCP<const panzer::GlobalIndexer> getDomainGlobalIndexer() const
-   { return blockProvider_; }
+   {
+      TEUCHOS_ASSERT(!gidProviders_.empty());
+      return gidProviders_.back()->getGlobalIndexer();
+   }
+
+   //! Get the domain unique global indexer this factory was created with.
+   Teuchos::RCP<const panzer::BlockedDOFManager> getDomainBlockedIndexer() const
+   {
+      TEUCHOS_ASSERT(!gidProviders_.empty());
+      return gidProviders_.back()->getBlockedIndexer();
+   }
 
    //! Get the range unique global indexer this factory was created with.
    Teuchos::RCP<const panzer::GlobalIndexer> getRangeGlobalIndexer() const
-   { return blockProvider_; }
+   {
+      TEUCHOS_ASSERT(!gidProviders_.empty());
+      return gidProviders_.front()->getGlobalIndexer();
+   }
 
 protected:
+   /*************** Utility class for handling blocked and nonblocked DOF managers *******************/
+
+   /** This classes is mean to abstract away the different global indexer types and hide
+    * if this is a blocked data structure or a unblocked one.
+    */
+   class DOFManagerContainer
+   {
+   public:
+      DOFManagerContainer() {}
+      DOFManagerContainer(const Teuchos::RCP<const GlobalIndexer> &ugi)
+      {
+         setGlobalIndexer(ugi);
+      }
+
+      void setGlobalIndexer(const Teuchos::RCP<const GlobalIndexer> &ugi)
+      {
+         using Teuchos::RCP;
+         using Teuchos::rcp_dynamic_cast;
+
+         auto blockedDOFManager = rcp_dynamic_cast<const BlockedDOFManager>(ugi);
+         auto flatDOFManager = rcp_dynamic_cast<const GlobalIndexer>(ugi);
+
+         if (blockedDOFManager != Teuchos::null)
+         {
+            // set BlockedDOFManager
+            blockedDOFManager_ = blockedDOFManager;
+
+            // get all GID providers
+            auto dofManagers = blockedDOFManager_->getFieldDOFManagers();
+            for (auto itr = dofManagers.begin(); itr != dofManagers.end(); ++itr)
+               gidProviders_.push_back(*itr);
+         }
+         else if (flatDOFManager != Teuchos::null)
+         {
+            // for absolute clarity, nullify the blockedDOFManager_
+            blockedDOFManager_ = Teuchos::null;
+
+            // you have only a single GID provider
+            gidProviders_.push_back(flatDOFManager);
+         }
+         else
+         {
+            TEUCHOS_ASSERT(false);
+         }
+      }
+
+      //! Get the number of global indexers (not including the blocked one) contained.
+      int getFieldBlocks() const
+      {
+         return Teuchos::as<int>(gidProviders_.size());
+      }
+
+      /** Return true if this contains a blocked DOFManager as opposed to only a single DOFManager.
+       * If this returns true then <code>getGlobalIndexer</code> will return a <code>BlockedDOFManager<int,GO></code>,
+       * other wise it will return a <code>GlobalIndexer<int,GO></code>.
+       */
+      bool containsBlockedDOFManager() const
+      {
+         return blockedDOFManager_ != Teuchos::null;
+      }
+
+      //! Get the "parent" global indexer (if <code>containsBlockedDOFManager()==false</code> this will throw)
+      Teuchos::RCP<const BlockedDOFManager> getBlockedIndexer() const
+      {
+         TEUCHOS_ASSERT(containsBlockedDOFManager());
+         return blockedDOFManager_;
+      }
+
+      //! Get the "parent" global indexer (if <code>getFieldBlocks()>1</code> this will be blocked, otherwise it may be either)
+      Teuchos::RCP<const GlobalIndexer> getGlobalIndexer() const
+      {
+         if (blockedDOFManager_ != Teuchos::null)
+            return blockedDOFManager_;
+
+         TEUCHOS_ASSERT(gidProviders_.size() == 1);
+         return gidProviders_[0];
+      }
+
+      //! Get DOFManagers associated with the blocks
+      const std::vector<Teuchos::RCP<const GlobalIndexer>> &getFieldDOFManagers() const
+      {
+         return gidProviders_;
+      }
+
+   private:
+      Teuchos::RCP<const BlockedDOFManager> blockedDOFManager_;
+      std::vector<Teuchos::RCP<const GlobalIndexer>> gidProviders_;
+   };
+
 /*************** Generic methods/members *******************/
 
    // Get the global indexer associated with a particular block
    Teuchos::RCP<const GlobalIndexer> getGlobalIndexer(int i) const;
+   Teuchos::RCP<const GlobalIndexer> getColGlobalIndexer(int i) const;
 
    //! Allocate the space in the std::vector objects so we can fill with appropriate Tpetra data
-   void makeRoomForBlocks(std::size_t blockCnt);
+   void makeRoomForBlocks(std::size_t blockCnt, std::size_t colBlockCnt = 0);
 
-   Teuchos::RCP<const GlobalIndexer> blockProvider_;
-   Teuchos::RCP<const BlockedDOFManager> blockedDOFManager_;
-   std::vector<Teuchos::RCP<const GlobalIndexer> > gidProviders_;
+   std::vector<Teuchos::RCP<const DOFManagerContainer> > gidProviders_;
+
+   bool useColGidProviders_;
 
    // which block entries are ignored
   std::unordered_set<std::pair<int,int>,panzer::pair_hash> excludedPairs_;
@@ -337,10 +478,10 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
 /*************** Thyra based methods/members *******************/
 
    void ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
-                                 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const;
+                                 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out, bool col) const;
    void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase<ScalarT> & in,Thyra::LinearOpBase<ScalarT> & out) const;
    void globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
-                                 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const;
+                                 const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out, bool col) const;
 
    mutable Teuchos::RCP<Thyra::ProductVectorSpaceBase<ScalarT> > rangeSpace_;
    mutable Teuchos::RCP<Thyra::ProductVectorSpaceBase<ScalarT> > domainSpace_;
@@ -356,13 +497,16 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
                                      const Teuchos::Ptr<CrsMatrixType> & A,
                                      bool zeroVectorRows) const;
 
-   void ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const;
+   void ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out, bool col) const;
    void ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType & out) const;
-   void globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const;
+   void globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out, bool col) const;
 
    // get the map from the matrix
    virtual Teuchos::RCP<const MapType> buildTpetraMap(int i) const;
+   virtual Teuchos::RCP<const MapType> buildColTpetraMap(int i) const;
+
    virtual Teuchos::RCP<const MapType> buildTpetraGhostedMap(int i) const;
+   virtual Teuchos::RCP<const MapType> buildColTpetraGhostedMap(int i) const;
 
    // get the graph of the crs matrix
    virtual Teuchos::RCP<const CrsGraphType> buildTpetraGraph(int i,int j) const;
@@ -371,12 +515,22 @@ class BlockedTpetraLinearObjFactory : public LinearObjFactory<Traits>
    // storage for Tpetra graphs and maps
    Teuchos::RCP<const Teuchos::MpiComm<int> > comm_;
    mutable std::vector<Teuchos::RCP<const MapType> > maps_;
+   mutable std::vector<Teuchos::RCP<const MapType> > colMaps_;
+
    mutable std::vector<Teuchos::RCP<const MapType> > ghostedMaps_;
+   mutable std::vector<Teuchos::RCP<const MapType> > colGhostedMaps_;
+
    mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> graphs_ ;
+   mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> colGraphs_ ;
+
    mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> ghostedGraphs_;
+   mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<const CrsGraphType>,panzer::pair_hash> colGhostedGraphs_;
 
    mutable std::vector<Teuchos::RCP<const ImportType> > importers_;
+   mutable std::vector<Teuchos::RCP<const ImportType> > colImporters_;
+
    mutable std::vector<Teuchos::RCP<const ExportType> > exporters_;
+   mutable std::vector<Teuchos::RCP<const ExportType> > colExporters_;
 };
 
 }
diff --git a/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp b/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp
index edd818f8ac0b..58a5338e09c1 100644
--- a/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp
+++ b/packages/panzer/disc-fe/src/lof/Panzer_BlockedTpetraLinearObjFactory_impl.hpp
@@ -66,6 +66,7 @@
 namespace panzer {
 
 using Teuchos::RCP;
+using Teuchos::rcp_dynamic_cast;
 
 // ************************************************************
 // class BlockedTpetraLinearObjFactory
@@ -75,12 +76,11 @@ template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename Globa
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
                               const Teuchos::RCP<const BlockedDOFManager> & gidProvider)
-   : blockProvider_(gidProvider), blockedDOFManager_(gidProvider), comm_(comm)
+   : useColGidProviders_(false), comm_(comm)
 {
-  for(std::size_t i=0;i<gidProvider->getFieldDOFManagers().size();i++)
-    gidProviders_.push_back(gidProvider->getFieldDOFManagers()[i]);
+  gidProviders_.push_back(Teuchos::rcp(new DOFManagerContainer(gidProvider)));
 
-  makeRoomForBlocks(gidProviders_.size());
+  makeRoomForBlocks(getBlockRowCount());
 
   // build and register the gather/scatter evaluators with
   // the base class.
@@ -91,9 +91,17 @@ template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename Globa
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 BlockedTpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
                               const std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> & gidProviders)
-  : gidProviders_(gidProviders), comm_(comm)
+  : useColGidProviders_(true), comm_(comm)
 {
-  makeRoomForBlocks(gidProviders_.size());
+   for(std::size_t i=0;i<gidProviders.size();i++) {
+      gidProviders_.push_back(Teuchos::rcp(new DOFManagerContainer(gidProviders[i])));
+   }
+
+  makeRoomForBlocks(getBlockRowCount(), getBlockColCount());
+
+  // build and register the gather/scatter evaluators with
+  // the base class.
+  this->buildGatherScatterEvaluators(*this);
 }
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
@@ -110,7 +118,7 @@ BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>
 buildLinearObjContainer() const
 {
    std::vector<Teuchos::RCP<const MapType> > blockMaps;
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t blockDim = getBlockRowCount();
    for(std::size_t i=0;i<blockDim;i++)
       blockMaps.push_back(getMap(i));
 
@@ -126,7 +134,7 @@ BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>
 buildGhostedLinearObjContainer() const
 {
    std::vector<Teuchos::RCP<const MapType> > blockMaps;
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t blockDim = getBlockRowCount();
    for(std::size_t i=0;i<blockDim;i++)
       blockMaps.push_back(getGhostedMap(i));
 
@@ -150,13 +158,13 @@ globalToGhostContainer(const LinearObjContainer & in,LinearObjContainer & out,in
    // Operations occur if the GLOBAL container has the correct targets!
    // Users set the GLOBAL continer arguments
    if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
-     globalToGhostThyraVector(b_in.get_x(),b_out.get_x());
+     globalToGhostThyraVector(b_in.get_x(),b_out.get_x(),true);
 
    if ( !is_null(b_in.get_dxdt()) && !is_null(b_out.get_dxdt()) && ((mem & LOC::DxDt)==LOC::DxDt))
-     globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt());
+     globalToGhostThyraVector(b_in.get_dxdt(),b_out.get_dxdt(),true);
 
    if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
-      globalToGhostThyraVector(b_in.get_f(),b_out.get_f());
+      globalToGhostThyraVector(b_in.get_f(),b_out.get_f(),false);
 }
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
@@ -173,10 +181,10 @@ ghostToGlobalContainer(const LinearObjContainer & in,LinearObjContainer & out,in
    // Operations occur if the GLOBAL container has the correct targets!
    // Users set the GLOBAL continer arguments
    if ( !is_null(b_in.get_x()) && !is_null(b_out.get_x()) && ((mem & LOC::X)==LOC::X))
-     ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x());
+     ghostToGlobalThyraVector(b_in.get_x(),b_out.get_x(),true);
 
    if ( !is_null(b_in.get_f()) && !is_null(b_out.get_f()) && ((mem & LOC::F)==LOC::F))
-     ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f());
+     ghostToGlobalThyraVector(b_in.get_f(),b_out.get_f(),false);
 
    if ( !is_null(b_in.get_A()) && !is_null(b_out.get_A()) && ((mem & LOC::Mat)==LOC::Mat))
      ghostToGlobalThyraMatrix(*b_in.get_A(),*b_out.get_A());
@@ -196,7 +204,8 @@ adjustForDirichletConditions(const LinearObjContainer & localBCRows,
    using Thyra::VectorBase;
    using Thyra::ProductVectorBase;
 
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t rowBlockDim = getBlockRowCount();
+   std::size_t colBlockDim = getBlockColCount();
 
    // first cast to block LOCs
    const BTLOC & b_localBCRows = Teuchos::dyn_cast<const BTLOC>(localBCRows);
@@ -215,13 +224,13 @@ adjustForDirichletConditions(const LinearObjContainer & localBCRows,
    if(adjustX) f = rcp_dynamic_cast<ProductVectorBase<ScalarT> >(b_ghosted.get_x());
 
    // sanity check!
-   if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) blockDim);
-   if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) blockDim);
-   if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) blockDim);
-   TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) blockDim);
-   TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) blockDim);
+   if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productRange()->numBlocks()==(int) rowBlockDim);
+   if(A!=Teuchos::null) TEUCHOS_ASSERT(A->productDomain()->numBlocks()==(int) colBlockDim);
+   if(f!=Teuchos::null) TEUCHOS_ASSERT(f->productSpace()->numBlocks()==(int) rowBlockDim);
+   TEUCHOS_ASSERT(local_bcs->productSpace()->numBlocks()==(int) rowBlockDim);
+   TEUCHOS_ASSERT(global_bcs->productSpace()->numBlocks()==(int) rowBlockDim);
 
-   for(std::size_t i=0;i<blockDim;i++) {
+   for(std::size_t i=0;i<rowBlockDim;i++) {
       // grab epetra vector
       RCP<const VectorType> t_local_bcs = rcp_dynamic_cast<const ThyraVector>(local_bcs->getVectorBlock(i),true)->getConstTpetraVector();
       RCP<const VectorType> t_global_bcs = rcp_dynamic_cast<const ThyraVector>(global_bcs->getVectorBlock(i),true)->getConstTpetraVector();
@@ -234,9 +243,9 @@ adjustForDirichletConditions(const LinearObjContainer & localBCRows,
       else
         t_f = rcp_dynamic_cast<ThyraVector>(th_f,true)->getTpetraVector();
 
-      for(std::size_t j=0;j<blockDim;j++) {
-        RCP<const MapType> map_i = getGhostedMap(i);
-        RCP<const MapType> map_j = getGhostedMap(j);
+      for(std::size_t j=0;j<colBlockDim;j++) {
+      //   RCP<const MapType> map_i = getGhostedMap(i);
+      //   RCP<const MapType> map_j = getGhostedMap(j);
 
          // pull out epetra values
          RCP<LinearOpBase<ScalarT> > th_A = (A== Teuchos::null)? Teuchos::null : A->getNonconstBlock(i,j);
@@ -359,10 +368,11 @@ buildReadOnlyDomainContainer() const
   using TVROGED = panzer::TpetraVector_ReadOnly_GlobalEvaluationData<ScalarT,
     LocalOrdinalT, GlobalOrdinalT, NodeT>;
   vector<RCP<ReadOnlyVector_GlobalEvaluationData>> gedBlocks;
+
   for (int i(0); i < getBlockColCount(); ++i)
   {
     auto tvroged = rcp(new TVROGED);
-    tvroged->initialize(getGhostedImport(i), getGhostedMap(i), getMap(i));
+    tvroged->initialize(getColGhostedImport(i), getColGhostedMap(i), getColMap(i));
     gedBlocks.push_back(tvroged);
   }
   auto ged = rcp(new BVROGED);
@@ -485,17 +495,31 @@ Teuchos::RCP<const GlobalIndexer>
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 getGlobalIndexer(int i) const
 {
-   return gidProviders_[i];
+   return gidProviders_[0]->getFieldDOFManagers()[i];
+}
+
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const GlobalIndexer>
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+getColGlobalIndexer(int i) const
+{
+   return gidProviders_.back()->getFieldDOFManagers()[i];
 }
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 void BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
-makeRoomForBlocks(std::size_t blockCnt)
+makeRoomForBlocks(std::size_t blockCnt, std::size_t colBlockCnt)
 {
    maps_.resize(blockCnt);
    ghostedMaps_.resize(blockCnt);
    importers_.resize(blockCnt);
    exporters_.resize(blockCnt);
+   if(colBlockCnt > 0) {
+      colMaps_.resize(blockCnt);
+      colGhostedMaps_.resize(blockCnt);
+      colImporters_.resize(blockCnt);
+      colExporters_.resize(blockCnt);
+   }
 }
 
 // Thyra methods
@@ -509,8 +533,8 @@ getThyraDomainSpace() const
    if(domainSpace_==Teuchos::null) {
       // loop over all vectors and build the vector space
       std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
-      for(std::size_t i=0;i<gidProviders_.size();i++)
-         vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
+      for(std::size_t i=0;i<getBlockColCount();i++)
+         vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getColMap(i)));
 
       domainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
    }
@@ -526,7 +550,7 @@ getThyraRangeSpace() const
    if(rangeSpace_==Teuchos::null) {
       // loop over all vectors and build the vector space
       std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
-      for(std::size_t i=0;i<gidProviders_.size();i++)
+      for(std::size_t i=0;i<getBlockRowCount();i++)
          vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getMap(i)));
 
       rangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
@@ -569,9 +593,9 @@ getThyraDomainVector() const
    Thyra::assign(vec.ptr(),0.0);
 
    Teuchos::RCP<Thyra::ProductVectorBase<ScalarT> > p_vec = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<ScalarT> >(vec);
-   for(std::size_t i=0;i<gidProviders_.size();i++) {
+   for(std::size_t i=0;i<getBlockColCount();i++) {
       TEUCHOS_ASSERT(Teuchos::rcp_dynamic_cast<Thyra::SpmdVectorBase<ScalarT> >(p_vec->getNonconstVectorBlock(i))->spmdSpace()->localSubDim()==
-                     Teuchos::as<int>(getMap(i)->getLocalNumElements()));
+                     Teuchos::as<int>(getColMap(i)->getLocalNumElements()));
    }
 
    return vec;
@@ -597,14 +621,15 @@ getThyraMatrix() const
    Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
 
    // get the block dimension
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t rowBlockDim = getBlockRowCount();
+   std::size_t colBlockDim = getBlockColCount();
 
    // this operator will be square
-   blockedOp->beginBlockFill(blockDim,blockDim);
+   blockedOp->beginBlockFill(rowBlockDim,colBlockDim);
 
    // loop over each block
-   for(std::size_t i=0;i<blockDim;i++) {
-      for(std::size_t j=0;j<blockDim;j++) {
+   for(std::size_t i=0;i<rowBlockDim;i++) {
+      for(std::size_t j=0;j<colBlockDim;j++) {
          if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
             // build (i,j) block matrix and add it to blocked operator
             Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getTpetraMatrix(i,j));
@@ -627,8 +652,8 @@ getGhostedThyraDomainSpace() const
    if(ghostedDomainSpace_==Teuchos::null) {
       // loop over all vectors and build the vector space
       std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
-      for(std::size_t i=0;i<gidProviders_.size();i++)
-         vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
+      for(std::size_t i=0;i<getBlockColCount();i++)
+         vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getColGhostedMap(i)));
 
       ghostedDomainSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
    }
@@ -644,7 +669,7 @@ getGhostedThyraRangeSpace() const
    if(ghostedRangeSpace_==Teuchos::null) {
       // loop over all vectors and build the vector space
       std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<ScalarT> > > vsArray;
-      for(std::size_t i=0;i<gidProviders_.size();i++)
+      for(std::size_t i=0;i<getBlockRowCount();i++)
          vsArray.push_back(Thyra::createVectorSpace<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedMap(i)));
 
       ghostedRangeSpace_ = Thyra::productVectorSpace<ScalarT>(vsArray);
@@ -685,14 +710,15 @@ getGhostedThyraMatrix() const
    Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<ScalarT> > blockedOp = Thyra::defaultBlockedLinearOp<ScalarT>();
 
    // get the block dimension
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t rowBlockDim = getBlockRowCount();
+   std::size_t colBlockDim = getBlockColCount();
 
    // this operator will be square
-   blockedOp->beginBlockFill(blockDim,blockDim);
+   blockedOp->beginBlockFill(rowBlockDim,colBlockDim);
 
    // loop over each block
-   for(std::size_t i=0;i<blockDim;i++) {
-      for(std::size_t j=0;j<blockDim;j++) {
+   for(std::size_t i=0;i<rowBlockDim;i++) {
+      for(std::size_t j=0;j<colBlockDim;j++) {
          if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
             // build (i,j) block matrix and add it to blocked operator
             Teuchos::RCP<Thyra::LinearOpBase<ScalarT> > block = Thyra::createLinearOp<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>(getGhostedTpetraMatrix(i,j));
@@ -710,13 +736,13 @@ getGhostedThyraMatrix() const
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 void BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
-                         const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
+                         const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out, bool col) const
 {
    using Teuchos::RCP;
    using Teuchos::rcp_dynamic_cast;
    using Thyra::ProductVectorBase;
 
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
 
    // get product vectors
    RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
@@ -731,7 +757,7 @@ ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > &
       RCP<VectorType> tp_out      = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
 
       // use Tpetra to do global communication
-      ghostToGlobalTpetraVector(i,*tp_in,*tp_out);
+      ghostToGlobalTpetraVector(i, *tp_in, *tp_out, col);
    }
 }
 
@@ -745,19 +771,20 @@ ghostToGlobalThyraMatrix(const Thyra::LinearOpBase<ScalarT> & in,Thyra::LinearOp
    using Thyra::LinearOpBase;
    using Thyra::PhysicallyBlockedLinearOpBase;
 
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t rowBlockDim = getBlockRowCount();
+   std::size_t colBlockDim = getBlockColCount();
 
    // get product vectors
    const PhysicallyBlockedLinearOpBase<ScalarT> & prod_in = dyn_cast<const PhysicallyBlockedLinearOpBase<ScalarT> >(in);
    PhysicallyBlockedLinearOpBase<ScalarT> & prod_out      = dyn_cast<PhysicallyBlockedLinearOpBase<ScalarT> >(out);
 
-   TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) blockDim);
-   TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) blockDim);
-   TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) blockDim);
-   TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) blockDim);
+   TEUCHOS_ASSERT(prod_in.productRange()->numBlocks()==(int) rowBlockDim);
+   TEUCHOS_ASSERT(prod_in.productDomain()->numBlocks()==(int) colBlockDim);
+   TEUCHOS_ASSERT(prod_out.productRange()->numBlocks()==(int) rowBlockDim);
+   TEUCHOS_ASSERT(prod_out.productDomain()->numBlocks()==(int) colBlockDim);
 
-   for(std::size_t i=0;i<blockDim;i++) {
-      for(std::size_t j=0;j<blockDim;j++) {
+   for(std::size_t i=0;i<rowBlockDim;i++) {
+      for(std::size_t j=0;j<colBlockDim;j++) {
          if(excludedPairs_.find(std::make_pair(i,j))==excludedPairs_.end()) {
             // extract the blocks
             RCP<const LinearOpBase<ScalarT> > th_in = prod_in.getBlock(i,j);
@@ -784,13 +811,13 @@ ghostToGlobalThyraMatrix(const Thyra::LinearOpBase<ScalarT> & in,Thyra::LinearOp
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 void BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > & in,
-                         const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out) const
+                         const Teuchos::RCP<Thyra::VectorBase<ScalarT> > & out, bool col) const
 {
    using Teuchos::RCP;
    using Teuchos::rcp_dynamic_cast;
    using Thyra::ProductVectorBase;
 
-   std::size_t blockDim = gidProviders_.size();
+   std::size_t blockDim = col ? getBlockColCount() : getBlockRowCount();
 
    // get product vectors
    RCP<const ProductVectorBase<ScalarT> > prod_in = rcp_dynamic_cast<const ProductVectorBase<ScalarT> >(in,true);
@@ -805,7 +832,7 @@ globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > &
       RCP<VectorType> tp_out      = rcp_dynamic_cast<ThyraVector>(prod_out->getNonconstVectorBlock(i),true)->getTpetraVector();
 
       // use Tpetra to do global communication
-      globalToGhostTpetraVector(i,*tp_in,*tp_out);
+      globalToGhostTpetraVector(i, *tp_in, *tp_out, col);
    }
 }
 
@@ -814,12 +841,12 @@ globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<ScalarT> > &
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 void BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
-ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out) const
+ghostToGlobalTpetraVector(int i,const VectorType & in,VectorType & out, bool col) const
 {
    using Teuchos::RCP;
 
    // do the global distribution
-   RCP<const ExportType> exporter = getGhostedExport(i);
+   RCP<const ExportType> exporter = col ? getColGhostedExport(i) : getGhostedExport(i);
    out.putScalar(0.0);
    out.doExport(in,*exporter,Tpetra::ADD);
 }
@@ -844,12 +871,12 @@ ghostToGlobalTpetraMatrix(int blockRow,const CrsMatrixType & in,CrsMatrixType &
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 void BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
-globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out) const
+globalToGhostTpetraVector(int i,const VectorType & in,VectorType & out, bool col) const
 {
    using Teuchos::RCP;
 
    // do the global distribution
-   RCP<const ImportType> importer = getGhostedImport(i);
+   RCP<const ImportType> importer = col ? getColGhostedImport(i) : getGhostedImport(i);
    out.putScalar(0.0);
    out.doImport(in,*importer,Tpetra::INSERT);
 }
@@ -866,6 +893,21 @@ getMap(int i) const
    return maps_[i];
 }
 
+// get the col map from the matrix
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+getColMap(int i) const
+{
+   if(not useColGidProviders_)
+     return getMap(i);
+
+   if(colMaps_[i]==Teuchos::null)
+      colMaps_[i] = buildColTpetraMap(i);
+
+   return colMaps_[i];
+}
+
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
@@ -877,6 +919,20 @@ getGhostedMap(int i) const
    return ghostedMaps_[i];
 }
 
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+getColGhostedMap(int i) const
+{
+   if(not useColGidProviders_)
+     return getGhostedMap(i);
+
+   if(colGhostedMaps_[i]==Teuchos::null)
+      colGhostedMaps_[i] = buildColTpetraGhostedMap(i);
+
+   return colGhostedMaps_[i];
+}
+
 // get the graph of the crs matrix
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
@@ -929,6 +985,20 @@ getGhostedImport(int i) const
    return importers_[i];
 }
 
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const  Tpetra::Import<LocalOrdinalT,GlobalOrdinalT,NodeT> >
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+getColGhostedImport(int i) const
+{
+   if(not useColGidProviders_)
+     return getGhostedImport(i);
+
+   if(colImporters_[i]==Teuchos::null)
+      colImporters_[i] = Teuchos::rcp(new ImportType(getColMap(i),getColGhostedMap(i)));
+
+   return importers_[i];
+}
+
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 Teuchos::RCP<const  Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
@@ -940,6 +1010,20 @@ getGhostedExport(int i) const
    return exporters_[i];
 }
 
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const  Tpetra::Export<LocalOrdinalT,GlobalOrdinalT,NodeT> >
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+getColGhostedExport(int i) const
+{
+   if(not useColGidProviders_)
+     return getGhostedExport(i);
+
+   if(colExporters_[i]==Teuchos::null)
+      colExporters_[i] = Teuchos::rcp(new ExportType(getColGhostedMap(i),getColMap(i)));
+
+   return exporters_[i];
+}
+
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
@@ -953,6 +1037,19 @@ buildTpetraMap(int i) const
    return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
 }
 
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+buildColTpetraMap(int i) const
+{
+   std::vector<GlobalOrdinalT> indices;
+
+   // get the global indices
+   getColGlobalIndexer(i)->getOwnedIndices(indices);
+
+   return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
+}
+
 // build the ghosted map
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
@@ -967,6 +1064,20 @@ buildTpetraGhostedMap(int i) const
    return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
 }
 
+// build the ghosted col map
+template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
+Teuchos::RCP<const Tpetra::Map<LocalOrdinalT,GlobalOrdinalT,NodeT> >
+BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
+buildColTpetraGhostedMap(int i) const
+{
+   std::vector<GlobalOrdinalT> indices;
+
+   // get the global indices
+   getColGlobalIndexer(i)->getOwnedAndGhostedIndices(indices);
+
+   return Teuchos::rcp(new MapType(Teuchos::OrdinalTraits<GlobalOrdinalT>::invalid(),indices,0,comm_));
+}
+
 // get the graph of the crs matrix
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
 Teuchos::RCP<const Tpetra::CrsGraph<LocalOrdinalT,GlobalOrdinalT,NodeT> >
@@ -979,7 +1090,10 @@ buildTpetraGraph(int i,int j) const
    // build the map and allocate the space for the graph and
    // grab the ghosted graph
    RCP<const MapType> map_i = getMap(i);
-   RCP<const MapType> map_j = getMap(j);
+   RCP<const MapType> map_j = getColMap(j);
+
+   TEUCHOS_ASSERT(map_i!=Teuchos::null);
+   TEUCHOS_ASSERT(map_j!=Teuchos::null);
 
    RCP<CrsGraphType> graph  = rcp(new CrsGraphType(map_i,0));
    RCP<const CrsGraphType> oGraph = getGhostedGraph(i,j);
@@ -1003,17 +1117,17 @@ buildTpetraGhostedGraph(int i,int j) const
    // build the map and allocate the space for the graph and
    // grab the ghosted graph
    RCP<const MapType> map_i = getGhostedMap(i);
-   RCP<const MapType> map_j = getGhostedMap(j);
+   RCP<const MapType> map_j = getColGhostedMap(j);
 
    std::vector<std::string> elementBlockIds;
 
    Teuchos::RCP<const GlobalIndexer> rowProvider, colProvider;
 
    rowProvider = getGlobalIndexer(i);
-   colProvider = getGlobalIndexer(j);
+   colProvider = getColGlobalIndexer(j);
 
-   gidProviders_[0]->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
-                                                          // same element blocks
+   rowProvider->getElementBlockIds(elementBlockIds); // each sub provider "should" have the
+                                                     // same element blocks
 
    // Count number of entries in each row of graph; needed for graph constructor
    std::vector<size_t> nEntriesPerRow(map_i->getLocalNumElements(), 0);
@@ -1021,8 +1135,8 @@ buildTpetraGhostedGraph(int i,int j) const
    for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
       std::string blockId = *blockItr;
       // grab elements for this block
-      const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
-                                                                                                // same elements in each element block
+      const std::vector<LocalOrdinalT> & elements = rowProvider->getElementBlock(blockId); // each sub provider "should" have the
+                                                                                           // same elements in each element block
 
       // get information about number of indicies
       std::vector<GlobalOrdinalT> row_gids;
@@ -1030,27 +1144,26 @@ buildTpetraGhostedGraph(int i,int j) const
 
       // loop over the elemnts
       for(std::size_t elmt=0;elmt<elements.size();elmt++) {
-
          rowProvider->getElementGIDs(elements[elmt],row_gids);
          colProvider->getElementGIDs(elements[elmt],col_gids);
+
          for(std::size_t row=0;row<row_gids.size();row++) {
             LocalOrdinalT lid = map_i->getLocalElement(row_gids[row]);
             nEntriesPerRow[lid] += col_gids.size();
          }
       }
    }
+
    Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
    RCP<CrsGraphType> graph  = rcp(new CrsGraphType(map_i,map_j, nEntriesPerRowView));
 
-
-
    // graph information about the mesh
    for(blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
       std::string blockId = *blockItr;
 
       // grab elements for this block
-      const std::vector<LocalOrdinalT> & elements = gidProviders_[0]->getElementBlock(blockId); // each sub provider "should" have the
-                                                                                                // same elements in each element block
+      const std::vector<LocalOrdinalT> & elements = rowProvider->getElementBlock(blockId); // each sub provider "should" have the
+                                                                                           // same elements in each element block
 
       // get information about number of indicies
       std::vector<GlobalOrdinalT> row_gids;
@@ -1058,7 +1171,6 @@ buildTpetraGhostedGraph(int i,int j) const
 
       // loop over the elemnts
       for(std::size_t elmt=0;elmt<elements.size();elmt++) {
-
          rowProvider->getElementGIDs(elements[elmt],row_gids);
          colProvider->getElementGIDs(elements[elmt],col_gids);
          for(std::size_t row=0;row<row_gids.size();row++)
@@ -1068,7 +1180,7 @@ buildTpetraGhostedGraph(int i,int j) const
 
    // finish filling the graph: Make sure the colmap and row maps coincide to
    //                           minimize calls to LID lookups
-   graph->fillComplete(getMap(j),getMap(i));
+   graph->fillComplete(map_j, map_i);
 
    return graph;
 }
@@ -1079,7 +1191,7 @@ BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>
 getTpetraMatrix(int i,int j) const
 {
    Teuchos::RCP<const MapType> map_i = getMap(i);
-   Teuchos::RCP<const MapType> map_j = getMap(j);
+   Teuchos::RCP<const MapType> map_j = getColMap(j);
 
    Teuchos::RCP<const CrsGraphType> tGraph = getGraph(i,j);
    Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
@@ -1094,11 +1206,11 @@ BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>
 getGhostedTpetraMatrix(int i,int j) const
 {
    Teuchos::RCP<const MapType> map_i = getGhostedMap(i);
-   Teuchos::RCP<const MapType> map_j = getGhostedMap(j);
+   Teuchos::RCP<const MapType> map_j = getColGhostedMap(j);
 
    Teuchos::RCP<const CrsGraphType> tGraph = getGhostedGraph(i,j);
    Teuchos::RCP<CrsMatrixType> mat = Teuchos::rcp(new CrsMatrixType(tGraph));
-   mat->fillComplete(getMap(j),getMap(i));
+   mat->fillComplete(map_j,map_i);
 
    return mat;
 }
@@ -1108,7 +1220,7 @@ Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 getTpetraDomainVector(int i) const
 {
-   Teuchos::RCP<const MapType> tMap = getMap(i);
+   Teuchos::RCP<const MapType> tMap = getColMap(i);
    return Teuchos::rcp(new VectorType(tMap));
 }
 
@@ -1117,7 +1229,7 @@ Teuchos::RCP<Tpetra::Vector<ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT> >
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 getGhostedTpetraDomainVector(int i) const
 {
-   Teuchos::RCP<const MapType> tMap = getGhostedMap(i);
+   Teuchos::RCP<const MapType> tMap = getColGhostedMap(i);
    return Teuchos::rcp(new VectorType(tMap));
 }
 
@@ -1144,7 +1256,7 @@ int
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 getBlockRowCount() const
 {
-   return gidProviders_.size();
+   return gidProviders_[0]->getFieldBlocks();
 }
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
@@ -1152,7 +1264,11 @@ int
 BlockedTpetraLinearObjFactory<Traits,ScalarT,LocalOrdinalT,GlobalOrdinalT,NodeT>::
 getBlockColCount() const
 {
-   return gidProviders_.size();
+   if(not useColGidProviders_) {
+      return gidProviders_[0]->getFieldBlocks();
+   }
+
+   return gidProviders_.back()->getFieldBlocks();
 }
 
 template <typename Traits,typename ScalarT,typename LocalOrdinalT,typename GlobalOrdinalT,typename NodeT>
diff --git a/packages/panzer/disc-fe/src/lof/Panzer_LinearObjFactory_Utilities.cpp b/packages/panzer/disc-fe/src/lof/Panzer_LinearObjFactory_Utilities.cpp
index 32729261eb49..39a9ff95887b 100644
--- a/packages/panzer/disc-fe/src/lof/Panzer_LinearObjFactory_Utilities.cpp
+++ b/packages/panzer/disc-fe/src/lof/Panzer_LinearObjFactory_Utilities.cpp
@@ -86,9 +86,10 @@ Teuchos::RCP<const LinearObjFactory<panzer::Traits> > cloneWithNewRangeAndDomain
 
   Ptr<const BlockedTpetraLOF> blk_tpetra_lof = ptr_dynamic_cast<const BlockedTpetraLOF>(ptrFromRef(lof));
   if(blk_tpetra_lof!=null) {
-    TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
-                               "panzer::cloneWithNewRangeAndDomain: Blocked Tpetra LOF does not yet support "
-                               "different range and domain indexers!");
+    auto rangeUGI  = (rUgi==null ? blk_tpetra_lof->getRangeGlobalIndexer() : rUgi);
+    auto domainUGI = (dUgi==null ? blk_tpetra_lof->getDomainGlobalIndexer() : dUgi);
+    auto mpiComm = rcp(new Teuchos::MpiComm<int>(blk_tpetra_lof->getComm()));
+    return rcp(new BlockedTpetraLOF(mpiComm, {rangeUGI, domainUGI}));
   }
 
   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
diff --git a/packages/panzer/disc-fe/test/CMakeLists.txt b/packages/panzer/disc-fe/test/CMakeLists.txt
index ed2d5a256571..7c75953d4d04 100644
--- a/packages/panzer/disc-fe/test/CMakeLists.txt
+++ b/packages/panzer/disc-fe/test/CMakeLists.txt
@@ -4,6 +4,7 @@ ADD_SUBDIRECTORY(closure_model)
 ADD_SUBDIRECTORY(closure_model_composite)
 ADD_SUBDIRECTORY(physics_block)
 ADD_SUBDIRECTORY(epetra_gatherscatter)
+ADD_SUBDIRECTORY(tpetra_gatherscatter)
 ADD_SUBDIRECTORY(la_factory)
 ADD_SUBDIRECTORY(evaluator_tests)
 ADD_SUBDIRECTORY(workset_tests)
diff --git a/packages/panzer/disc-fe/test/la_factory/tTpetra_tCloneLOF.cpp b/packages/panzer/disc-fe/test/la_factory/tTpetra_tCloneLOF.cpp
index 20d2583a0b15..01590fa88c2a 100644
--- a/packages/panzer/disc-fe/test/la_factory/tTpetra_tCloneLOF.cpp
+++ b/packages/panzer/disc-fe/test/la_factory/tTpetra_tCloneLOF.cpp
@@ -139,9 +139,8 @@ namespace panzer
 
    TEUCHOS_UNIT_TEST(tCloneLOF, blocked_tpetra)
    {
-      // TODO: Uncomment when BlockedTpetraLinearObjFactory will be supported in cloneWithNewDomain method
-      // typedef Thyra::ProductVectorBase<double> PVector;
-      // typedef Thyra::BlockedLinearOpBase<double> BLinearOp;
+      typedef Thyra::ProductVectorBase<double> PVector;
+      typedef Thyra::BlockedLinearOpBase<double> BLinearOp;
 
 // build global (or serial communicator)
 #ifdef HAVE_MPI
@@ -160,7 +159,7 @@ namespace panzer
 
       RCP<panzer::BlockedDOFManager> indexer = rcp(new panzer::BlockedDOFManager());
       {
-         std::vector<std::vector<std::string>> fieldOrder(2);
+         std::vector<std::vector<std::string> > fieldOrder(2);
 
          indexer->setConnManager(connManager, MPI_COMM_WORLD);
          indexer->addField("U", patternC1);
@@ -192,49 +191,45 @@ namespace panzer
 
       // setup factory
       RCP<BlockedTpetraLinObjFactoryType> bt_lof = Teuchos::rcp(new BlockedTpetraLinObjFactoryType(tComm, indexer));
-
-      // NOT supported yet
-      TEST_THROW(cloneWithNewDomain(*bt_lof, control_indexer), std::logic_error);
-
-      // TODO: Uncomment when BlockedTpetraLinearObjFactory will be supported in cloneWithNewDomain method
-      // // this is the member we are testing!
-      // RCP<const LinearObjFactory<Traits>> control_lof = cloneWithNewDomain(*bt_lof, control_indexer);
-      // RCP<const BlockedTpetraLinObjFactoryType> bt_control_lof = rcp_dynamic_cast<const BlockedTpetraLinObjFactoryType>(control_lof);
-
-      // RCP<BLinearOp> mat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getThyraMatrix());
-      // RCP<BLinearOp> gmat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getGhostedThyraMatrix());
-      // RCP<PVector> x = rcp_dynamic_cast<PVector>(bt_control_lof->getThyraDomainVector());
-      // RCP<PVector> gx = rcp_dynamic_cast<PVector>(bt_control_lof->getGhostedThyraDomainVector());
-      // RCP<PVector> f = rcp_dynamic_cast<PVector>(bt_control_lof->getThyraRangeVector());
-      // RCP<PVector> gf = rcp_dynamic_cast<PVector>(bt_control_lof->getGhostedThyraRangeVector());
-
-      // TEST_EQUALITY(x->productSpace()->numBlocks(), 1);
-      // TEST_EQUALITY(x->productSpace()->dim(), 18);
-      // TEST_EQUALITY(gx->productSpace()->numBlocks(), 1);
-      // TEST_EQUALITY(gx->productSpace()->dim(), 10 + 15);
-
-      // TEST_EQUALITY(f->productSpace()->numBlocks(), 2);
-      // TEST_EQUALITY(f->productSpace()->dim(), 36);
-      // TEST_EQUALITY(gf->productSpace()->numBlocks(), 2);
-      // TEST_EQUALITY(gf->productSpace()->dim(), 50);
-
-      // TEST_EQUALITY(mat->productRange()->numBlocks(), 2);
-      // TEST_EQUALITY(mat->productRange()->dim(), 36);
-      // TEST_EQUALITY(mat->productDomain()->numBlocks(), 1);
-      // TEST_EQUALITY(mat->productDomain()->dim(), 18);
-
-      // TEST_EQUALITY(gmat->productRange()->numBlocks(), 2);
-      // TEST_EQUALITY(gmat->productRange()->dim(), 50);
-      // TEST_EQUALITY(gmat->productDomain()->numBlocks(), 1);
-      // TEST_EQUALITY(gmat->productDomain()->dim(), 10 + 15);
+      RCP<const LinearObjFactory<Traits>> control_lof = cloneWithNewDomain(*bt_lof, control_indexer);
+
+      TEST_ASSERT(nullptr != control_lof.getRawPtr());
+
+      RCP<const BlockedTpetraLinObjFactoryType> bt_control_lof = rcp_dynamic_cast<const BlockedTpetraLinObjFactoryType>(control_lof, true);
+
+      RCP<BLinearOp> mat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getThyraMatrix());
+      RCP<BLinearOp> gmat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getGhostedThyraMatrix());
+      RCP<PVector> x = rcp_dynamic_cast<PVector>(bt_control_lof->getThyraDomainVector());
+      RCP<PVector> gx = rcp_dynamic_cast<PVector>(bt_control_lof->getGhostedThyraDomainVector());
+      RCP<PVector> f = rcp_dynamic_cast<PVector>(bt_control_lof->getThyraRangeVector());
+      RCP<PVector> gf = rcp_dynamic_cast<PVector>(bt_control_lof->getGhostedThyraRangeVector());
+
+      TEST_EQUALITY(x->productSpace()->numBlocks(), 1);
+      TEST_EQUALITY(x->productSpace()->dim(), 18);
+      TEST_EQUALITY(gx->productSpace()->numBlocks(), 1);
+      TEST_EQUALITY(gx->productSpace()->dim(), 10 + 15);
+
+      TEST_EQUALITY(f->productSpace()->numBlocks(), 2);
+      TEST_EQUALITY(f->productSpace()->dim(), 36);
+      TEST_EQUALITY(gf->productSpace()->numBlocks(), 2);
+      TEST_EQUALITY(gf->productSpace()->dim(), 50);
+
+      TEST_EQUALITY(mat->productRange()->numBlocks(), 2);
+      TEST_EQUALITY(mat->productRange()->dim(), 36);
+      TEST_EQUALITY(mat->productDomain()->numBlocks(), 1);
+      TEST_EQUALITY(mat->productDomain()->dim(), 18);
+
+      TEST_EQUALITY(gmat->productRange()->numBlocks(), 2);
+      TEST_EQUALITY(gmat->productRange()->dim(), 50);
+      TEST_EQUALITY(gmat->productDomain()->numBlocks(), 1);
+      TEST_EQUALITY(gmat->productDomain()->dim(), 10 + 15);
    }
 
    TEUCHOS_UNIT_TEST(tCloneLOF, blocked_tpetra_nonblocked_domain)
    {
-      // TODO: Uncomment when BlockedTpetraLinearObjFactory will be supported in cloneWithNewDomain method
-      // typedef Thyra::ProductVectorBase<double> PVector;
-      // typedef Thyra::BlockedLinearOpBase<double> BLinearOp;
-      // typedef Thyra::VectorBase<double> Vector;
+      typedef Thyra::ProductVectorBase<double> PVector;
+      typedef Thyra::BlockedLinearOpBase<double> BLinearOp;
+      typedef Thyra::VectorBase<double> Vector;
 
 // build global (or serial communicator)
 #ifdef HAVE_MPI
@@ -275,45 +270,39 @@ namespace panzer
       }
 
       // setup factory
-      out << "build lof" << std::endl;
       RCP<BlockedTpetraLinObjFactoryType> bt_lof = Teuchos::rcp(new BlockedTpetraLinObjFactoryType(tComm, indexer));
 
-      // NOT supported yet
-      TEST_THROW(cloneWithNewDomain(*bt_lof, control_indexer), std::logic_error);
-
-      // TODO: Uncomment when BlockedTpetraLinearObjFactory will be supported in cloneWithNewDomain method
-      // // this is the member we are testing!
-      // out << "cloning lof" << std::endl;
-      // RCP<const LinearObjFactory<Traits>> control_lof = cloneWithNewDomain(*bt_lof, control_indexer);
-
-      // out << "casting lof" << std::endl;
-      // RCP<const BlockedTpetraLinObjFactoryType> bt_control_lof = rcp_dynamic_cast<const BlockedTpetraLinObjFactoryType>(control_lof, true);
-
-      // out << "using casted lof" << std::endl;
-      // RCP<BLinearOp> mat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getThyraMatrix(), true);
-      // RCP<BLinearOp> gmat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getGhostedThyraMatrix(), true);
-      // RCP<Vector> x = bt_control_lof->getThyraDomainVector();
-      // RCP<Vector> gx = bt_control_lof->getGhostedThyraDomainVector();
-      // RCP<PVector> f = rcp_dynamic_cast<PVector>(bt_control_lof->getThyraRangeVector(), true);
-      // RCP<PVector> gf = rcp_dynamic_cast<PVector>(bt_control_lof->getGhostedThyraRangeVector(), true);
-
-      // TEST_EQUALITY(x->space()->dim(), 18);
-      // TEST_EQUALITY(gx->space()->dim(), 10 + 15);
-
-      // TEST_EQUALITY(f->productSpace()->numBlocks(), 2);
-      // TEST_EQUALITY(f->productSpace()->dim(), 36);
-      // TEST_EQUALITY(gf->productSpace()->numBlocks(), 2);
-      // TEST_EQUALITY(gf->productSpace()->dim(), 50);
-
-      // TEST_EQUALITY(mat->productRange()->numBlocks(), 2);
-      // TEST_EQUALITY(mat->productRange()->dim(), 36);
-      // TEST_EQUALITY(mat->productDomain()->numBlocks(), 1);
-      // TEST_EQUALITY(mat->productDomain()->dim(), 18);
-
-      // TEST_EQUALITY(gmat->productRange()->numBlocks(), 2);
-      // TEST_EQUALITY(gmat->productRange()->dim(), 50);
-      // TEST_EQUALITY(gmat->productDomain()->numBlocks(), 1);
-      // TEST_EQUALITY(gmat->productDomain()->dim(), 10 + 15);
+      // this is the member we are testing!
+      RCP<const LinearObjFactory<Traits>> control_lof = cloneWithNewDomain(*bt_lof, control_indexer);
+      TEST_ASSERT(nullptr != control_lof.getRawPtr());
+
+      RCP<const BlockedTpetraLinObjFactoryType> bt_control_lof = rcp_dynamic_cast<const BlockedTpetraLinObjFactoryType>(control_lof, true);
+      TEST_ASSERT(nullptr != bt_control_lof.getRawPtr());
+
+      RCP<BLinearOp> mat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getThyraMatrix(), true);
+      RCP<BLinearOp> gmat = rcp_dynamic_cast<BLinearOp>(bt_control_lof->getGhostedThyraMatrix(), true);
+      RCP<Vector> x = bt_control_lof->getThyraDomainVector();
+      RCP<Vector> gx = bt_control_lof->getGhostedThyraDomainVector();
+      RCP<PVector> f = rcp_dynamic_cast<PVector>(bt_control_lof->getThyraRangeVector(), true);
+      RCP<PVector> gf = rcp_dynamic_cast<PVector>(bt_control_lof->getGhostedThyraRangeVector(), true);
+
+      TEST_EQUALITY(x->space()->dim(), 18);
+      TEST_EQUALITY(gx->space()->dim(), 10 + 15);
+
+      TEST_EQUALITY(f->productSpace()->numBlocks(), 2);
+      TEST_EQUALITY(f->productSpace()->dim(), 36);
+      TEST_EQUALITY(gf->productSpace()->numBlocks(), 2);
+      TEST_EQUALITY(gf->productSpace()->dim(), 50);
+
+      TEST_EQUALITY(mat->productRange()->numBlocks(), 2);
+      TEST_EQUALITY(mat->productRange()->dim(), 36);
+      TEST_EQUALITY(mat->productDomain()->numBlocks(), 1);
+      TEST_EQUALITY(mat->productDomain()->dim(), 18);
+
+      TEST_EQUALITY(gmat->productRange()->numBlocks(), 2);
+      TEST_EQUALITY(gmat->productRange()->dim(), 50);
+      TEST_EQUALITY(gmat->productDomain()->numBlocks(), 1);
+      TEST_EQUALITY(gmat->productDomain()->dim(), 10 + 15);
    }
 
 }
diff --git a/packages/panzer/disc-fe/test/physics_block/physics_block.cpp b/packages/panzer/disc-fe/test/physics_block/physics_block.cpp
index 2dd66ac9256f..63a50df98c9d 100644
--- a/packages/panzer/disc-fe/test/physics_block/physics_block.cpp
+++ b/packages/panzer/disc-fe/test/physics_block/physics_block.cpp
@@ -55,6 +55,7 @@
 #ifdef PANZER_HAVE_EPETRA_STACK
 #include "Panzer_BlockedEpetraLinearObjFactory.hpp"
 #endif
+#include "Panzer_BlockedTpetraLinearObjFactory.hpp"
 
 #include "user_app_EquationSetFactory.hpp"
 #include "Panzer_ClosureModel_Factory_TemplateManager.hpp"
@@ -238,6 +239,32 @@ namespace panzer {
   }
 #endif
 
+  TEUCHOS_UNIT_TEST(physics_block, nontemplate_tpetra_evaluator_builders)
+  {
+    Teuchos::RCP<const Teuchos::MpiComm<int>> comm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));
+    Teuchos::RCP<const BlockedDOFManager> dof = Teuchos::rcp(new BlockedDOFManager());
+    auto btlof = panzer::BlockedTpetraLinearObjFactory<panzer::Traits, double, panzer::LocalOrdinal, panzer::GlobalOrdinal>(comm, dof);
+
+    Teuchos::RCP<panzer::PhysicsBlock> physics_block =
+        panzer_test_utils::createPhysicsBlock();
+
+    PHX::FieldManager<panzer::Traits> fm;
+
+    Teuchos::ParameterList user_data("User Data");
+
+    physics_block->buildAndRegisterEquationSetEvaluators(fm, user_data);
+    physics_block->buildAndRegisterGatherAndOrientationEvaluators(fm, btlof, user_data);
+    physics_block->buildAndRegisterDOFProjectionsToIPEvaluators(fm, Teuchos::null, user_data);
+    physics_block->buildAndRegisterScatterEvaluators(fm, btlof, user_data);
+
+    Teuchos::RCP<panzer::ClosureModelFactory_TemplateManager<panzer::Traits>> factory =
+        panzer_test_utils::buildModelFactory();
+
+    Teuchos::RCP<Teuchos::ParameterList> models = panzer_test_utils::buildModelDescriptors();
+
+    physics_block->buildAndRegisterClosureModelEvaluators(fm, *factory, *models, user_data);
+  }
+
   TEUCHOS_UNIT_TEST(physics_block, elementBlockID)
   {
 
@@ -287,4 +314,38 @@ namespace panzer {
   }
 #endif
 
+  TEUCHOS_UNIT_TEST(physics_block, templated_tpetra_evaluator_builders)
+  {
+    Teuchos::RCP<const Teuchos::MpiComm<int>> comm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));
+    Teuchos::RCP<const BlockedDOFManager> dof = Teuchos::rcp(new BlockedDOFManager());
+    auto btlof = panzer::BlockedTpetraLinearObjFactory<panzer::Traits, double, panzer::LocalOrdinal, panzer::GlobalOrdinal>(comm, dof);
+
+    Teuchos::RCP<panzer::PhysicsBlock> physics_block =
+      panzer_test_utils::createPhysicsBlock();
+
+    Teuchos::ParameterList user_data("User Data");
+
+    PHX::FieldManager<panzer::Traits> fm;
+
+    physics_block->buildAndRegisterEquationSetEvaluatorsForType<panzer::Traits::Residual>(fm, user_data);
+    physics_block->buildAndRegisterEquationSetEvaluatorsForType<panzer::Traits::Jacobian>(fm, user_data);
+
+    physics_block->buildAndRegisterGatherAndOrientationEvaluatorsForType<panzer::Traits::Residual>(fm, btlof, user_data);
+    physics_block->buildAndRegisterGatherAndOrientationEvaluatorsForType<panzer::Traits::Jacobian>(fm, btlof, user_data);
+
+    physics_block->buildAndRegisterDOFProjectionsToIPEvaluatorsForType<panzer::Traits::Residual>(fm, Teuchos::null, user_data);
+    physics_block->buildAndRegisterDOFProjectionsToIPEvaluatorsForType<panzer::Traits::Jacobian>(fm, Teuchos::null, user_data);
+
+    physics_block->buildAndRegisterScatterEvaluatorsForType<panzer::Traits::Residual>(fm, btlof, user_data);
+    physics_block->buildAndRegisterScatterEvaluatorsForType<panzer::Traits::Jacobian>(fm, btlof, user_data);
+
+    Teuchos::RCP<panzer::ClosureModelFactory_TemplateManager<panzer::Traits> > factory =
+      panzer_test_utils::buildModelFactory();
+
+    Teuchos::RCP<Teuchos::ParameterList> models = panzer_test_utils::buildModelDescriptors();
+
+    physics_block->buildAndRegisterClosureModelEvaluatorsForType<panzer::Traits::Residual>(fm, *factory, *models, user_data);
+    physics_block->buildAndRegisterClosureModelEvaluatorsForType<panzer::Traits::Jacobian>(fm, *factory, *models, user_data);
+  }
+
 }
diff --git a/packages/panzer/disc-fe/test/tpetra_gatherscatter/CMakeLists.txt b/packages/panzer/disc-fe/test/tpetra_gatherscatter/CMakeLists.txt
new file mode 100644
index 000000000000..6a7e5abf72cb
--- /dev/null
+++ b/packages/panzer/disc-fe/test/tpetra_gatherscatter/CMakeLists.txt
@@ -0,0 +1,24 @@
+
+
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
+
+SET(UNIT_TEST_DRIVER
+  ${PANZER_UNIT_TEST_MAIN})
+
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+  tTpetraGather
+  SOURCES tTpetraGather.cpp ${UNIT_TEST_DRIVER}
+  COMM serial mpi
+  )
+
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+  tTpetraScatter
+  SOURCES tTpetraScatter.cpp ${UNIT_TEST_DRIVER}
+  COMM serial mpi
+  )
+
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+  tTpetraScatterDirichlet
+  SOURCES tTpetraScatterDirichlet.cpp ${UNIT_TEST_DRIVER}
+  COMM serial mpi
+  )
diff --git a/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraGather.cpp b/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraGather.cpp
new file mode 100644
index 000000000000..1d404224ed95
--- /dev/null
+++ b/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraGather.cpp
@@ -0,0 +1,125 @@
+// @HEADER
+// ***********************************************************************
+//
+//           Panzer: A partial differential equation assembly
+//       engine for strongly coupled complex multiphysics systems
+//                 Copyright (2011) Sandia Corporation
+//
+// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+// the U.S. Government retains certain rights in this software.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the Corporation nor the names of the
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
+// Eric C. Cyr (eccyr@sandia.gov)
+// ***********************************************************************
+// @HEADER
+
+#include <Teuchos_ConfigDefs.hpp>
+#include <Teuchos_UnitTestHarness.hpp>
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_TimeMonitor.hpp>
+
+#include <string>
+#include <iostream>
+#include <vector>
+
+#include "Panzer_GatherSolution_Tpetra.hpp"
+#include "Panzer_PureBasis.hpp"
+#include "Panzer_IntegrationRule.hpp"
+#include "Panzer_CellData.hpp"
+
+namespace panzer {
+
+TEUCHOS_UNIT_TEST(tTpetraGather, constructor)
+{
+
+   using Teuchos::RCP;
+   using Teuchos::rcp;
+
+   typedef panzer::Traits::Residual Residual;
+   typedef panzer::Traits::Jacobian Jacobian;
+
+   Teuchos::RCP<shards::CellTopology> topo 
+    = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));
+
+   // auxiliary information needed to construct basis object
+   std::size_t numCells = 10;
+   std::string basisType = "Q1";
+   panzer::CellData cellData(numCells,topo);
+
+   // build DOF names
+   RCP<std::vector<std::string> > dofNames = rcp(new std::vector<std::string>); 
+   dofNames->push_back("ux"); // in practice these probably would not be gathered together!
+   dofNames->push_back("p");
+
+   // build basis
+   RCP<panzer::PureBasis> basis = rcp(new panzer::PureBasis(basisType,1,cellData));
+
+   // build gather parameter list
+   Teuchos::ParameterList gatherParams;
+   gatherParams.set<RCP<std::vector<std::string> > >("DOF Names",dofNames);
+   gatherParams.set<RCP<std::vector<std::string> > >("Indexer Names",dofNames);
+   gatherParams.set<RCP<panzer::PureBasis> >("Basis",basis);
+
+   // test residual gather evaluator
+   {
+      panzer::GatherSolution_Tpetra<Residual,panzer::Traits,int,panzer::GlobalOrdinal> gatherResidual(Teuchos::null,gatherParams);
+
+      const std::vector<RCP<PHX::FieldTag> > & fields = gatherResidual.evaluatedFields();
+      TEST_EQUALITY(fields.size(),2);
+ 
+      TEST_EQUALITY(fields[0]->name(),"ux");
+      TEST_EQUALITY(fields[1]->name(),"p");
+
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+   }
+
+   // test jacobian gather evaluator
+   {
+      panzer::GatherSolution_Tpetra<Jacobian,panzer::Traits,int,panzer::GlobalOrdinal> gatherJacobian(Teuchos::null,gatherParams);
+
+      const std::vector<RCP<PHX::FieldTag> > & fields = gatherJacobian.evaluatedFields();
+      TEST_EQUALITY(fields.size(),2);
+ 
+      TEST_EQUALITY(fields[0]->name(),"ux");
+      TEST_EQUALITY(fields[1]->name(),"p");
+
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+   }
+}
+
+}
diff --git a/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraScatter.cpp b/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraScatter.cpp
new file mode 100644
index 000000000000..1676fcc8a83d
--- /dev/null
+++ b/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraScatter.cpp
@@ -0,0 +1,140 @@
+// @HEADER
+// ***********************************************************************
+//
+//           Panzer: A partial differential equation assembly
+//       engine for strongly coupled complex multiphysics systems
+//                 Copyright (2011) Sandia Corporation
+//
+// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+// the U.S. Government retains certain rights in this software.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the Corporation nor the names of the
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
+// Eric C. Cyr (eccyr@sandia.gov)
+// ***********************************************************************
+// @HEADER
+
+#include <Teuchos_ConfigDefs.hpp>
+#include <Teuchos_UnitTestHarness.hpp>
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_TimeMonitor.hpp>
+
+#include <string>
+#include <iostream>
+#include <vector>
+#include <map>
+
+#include "Panzer_ScatterResidual_Tpetra.hpp"
+#include "Panzer_PureBasis.hpp"
+#include "Panzer_IntegrationRule.hpp"
+#include "Panzer_CellData.hpp"
+
+namespace panzer {
+
+TEUCHOS_UNIT_TEST(tTpetraScatter, constructor)
+{
+
+   using Teuchos::RCP;
+   using Teuchos::rcp;
+
+   typedef panzer::Traits::Residual Residual;
+   typedef panzer::Traits::Jacobian Jacobian;
+
+   Teuchos::RCP<shards::CellTopology> topo
+    = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));
+
+   // auxiliary information needed to construct basis object
+   std::size_t numCells = 10;
+   std::string basisType = "Q1";
+   panzer::CellData cellData(numCells,topo);
+
+   // build basis
+   RCP<const panzer::PureBasis> basis = rcp(new panzer::PureBasis(basisType,1,cellData));
+
+   std::string scatterName = "Residual_NS";
+
+   // build DOF names
+   RCP<std::vector<std::string> > evaluatedNames = rcp(new std::vector<std::string>); 
+   evaluatedNames->push_back("Residual_ux"); // in practice these probably would not be scattered together!
+   evaluatedNames->push_back("Residual_p");
+
+   // build evaluated map
+   RCP<std::map<std::string,std::string> > evaluatedMap = rcp(new std::map<std::string,std::string>); 
+   evaluatedMap->insert(std::make_pair("Residual_ux","ux")); // in practice these probably would not be scattered together!
+   evaluatedMap->insert(std::make_pair("Residual_p","p"));
+
+   // build scatter parameter list
+   Teuchos::ParameterList scatterParams;
+   scatterParams.set<std::string>("Scatter Name",scatterName);
+   scatterParams.set<RCP<std::vector<std::string> > >("Dependent Names",evaluatedNames);
+   scatterParams.set<RCP<std::map<std::string,std::string> > >("Dependent Map",evaluatedMap);
+   scatterParams.set("Basis",basis);
+
+   // test residual scatter evaluator
+   {
+      panzer::ScatterResidual_Tpetra<Residual,panzer::Traits,int,panzer::GlobalOrdinal> scatterResidual(Teuchos::null,scatterParams);
+
+      const std::vector<RCP<PHX::FieldTag> > & evalFields = scatterResidual.evaluatedFields();
+      TEST_EQUALITY(evalFields.size(),1); // this is a dummy holder for the sake of the field manager
+
+      const std::vector<RCP<PHX::FieldTag> > & fields = scatterResidual.dependentFields();
+      TEST_EQUALITY(fields.size(),2); // these store the residual values
+
+      TEST_EQUALITY(fields[0]->name(),"Residual_ux");
+      TEST_EQUALITY(fields[1]->name(),"Residual_p");
+
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+   }
+
+   // test jacobian scatter evaluator
+   {
+      panzer::ScatterResidual_Tpetra<Jacobian,panzer::Traits,int,panzer::GlobalOrdinal> scatterJacobian(Teuchos::null,scatterParams);
+
+      const std::vector<RCP<PHX::FieldTag> > & evalFields = scatterJacobian.evaluatedFields();
+      TEST_EQUALITY(evalFields.size(),1); // this is a dummy holder for the sake of the field manager
+
+      const std::vector<RCP<PHX::FieldTag> > & fields = scatterJacobian.dependentFields();
+      TEST_EQUALITY(fields.size(),2); // these store the residual values
+
+      TEST_EQUALITY(fields[0]->name(),"Residual_ux");
+      TEST_EQUALITY(fields[1]->name(),"Residual_p");
+
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+   }
+}
+
+}
diff --git a/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraScatterDirichlet.cpp b/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraScatterDirichlet.cpp
new file mode 100644
index 000000000000..6804a03c6632
--- /dev/null
+++ b/packages/panzer/disc-fe/test/tpetra_gatherscatter/tTpetraScatterDirichlet.cpp
@@ -0,0 +1,143 @@
+// @HEADER
+// ***********************************************************************
+//
+//           Panzer: A partial differential equation assembly
+//       engine for strongly coupled complex multiphysics systems
+//                 Copyright (2011) Sandia Corporation
+//
+// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+// the U.S. Government retains certain rights in this software.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the Corporation nor the names of the
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
+// Eric C. Cyr (eccyr@sandia.gov)
+// ***********************************************************************
+// @HEADER
+
+#include <Teuchos_ConfigDefs.hpp>
+#include <Teuchos_UnitTestHarness.hpp>
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_TimeMonitor.hpp>
+
+#include <string>
+#include <iostream>
+#include <vector>
+#include <map>
+
+#include "Panzer_ScatterDirichletResidual_Tpetra.hpp"
+#include "Panzer_PureBasis.hpp"
+#include "Panzer_IntegrationRule.hpp"
+#include "Panzer_CellData.hpp"
+
+namespace panzer {
+
+TEUCHOS_UNIT_TEST(tTpetraScatter, constructor)
+{
+
+   using Teuchos::RCP;
+   using Teuchos::rcp;
+
+   typedef panzer::Traits::Residual Residual;
+   typedef panzer::Traits::Jacobian Jacobian;
+
+   Teuchos::RCP<shards::CellTopology> topo
+     = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));
+
+   // auxiliary information needed to construct basis object
+   std::size_t numCells = 10;
+   std::string basisType = "Q1";
+   panzer::CellData cellData(numCells,topo);
+
+   // build basis
+   RCP<panzer::PureBasis> basis = rcp(new panzer::PureBasis(basisType,1,cellData));
+
+   std::string scatterName = "Residual_NS";
+
+   // build DOF names
+   RCP<std::vector<std::string> > evaluatedNames = rcp(new std::vector<std::string>); 
+   evaluatedNames->push_back("Residual_ux"); // in practice these probably would not be scattered together!
+   evaluatedNames->push_back("Residual_p");
+
+   // build evaluated map
+   RCP<std::map<std::string,std::string> > evaluatedMap = rcp(new std::map<std::string,std::string>); 
+   evaluatedMap->insert(std::make_pair("Residual_ux","ux")); // in practice these probably would not be scattered together!
+   evaluatedMap->insert(std::make_pair("Residual_p","p"));
+
+   // build scatter parameter list
+   Teuchos::ParameterList scatterParams;
+   scatterParams.set<std::string>("Scatter Name",scatterName);
+   scatterParams.set<RCP<std::vector<std::string> > >("Dependent Names",evaluatedNames);
+   scatterParams.set<RCP<std::map<std::string,std::string> > >("Dependent Map",evaluatedMap);
+   scatterParams.set<RCP<panzer::PureBasis> >("Basis",basis);
+   scatterParams.set<int>("Side Subcell Dimension",1);
+   scatterParams.set<int>("Local Side ID",2);
+   scatterParams.set<bool>("Check Apply BC",false);
+
+   // test residual scatter evaluator
+   {
+      panzer::ScatterDirichletResidual_Tpetra<Residual,panzer::Traits,int,panzer::GlobalOrdinal> scatterResidual(Teuchos::null, scatterParams);
+
+      const std::vector<RCP<PHX::FieldTag> > & evalFields = scatterResidual.evaluatedFields();
+      TEST_EQUALITY(evalFields.size(),1); // this is a dummy holder for the sake of the field manager
+
+      const std::vector<RCP<PHX::FieldTag> > & fields = scatterResidual.dependentFields();
+      TEST_EQUALITY(fields.size(),2); // these store the residual values
+
+      TEST_EQUALITY(fields[0]->name(),"Residual_ux");
+      TEST_EQUALITY(fields[1]->name(),"Residual_p");
+
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+   }
+
+   // test jacobian scatter evaluator
+   {
+      panzer::ScatterDirichletResidual_Tpetra<Jacobian,panzer::Traits,int,panzer::GlobalOrdinal> scatterJacobian(Teuchos::null, scatterParams);
+
+      const std::vector<RCP<PHX::FieldTag> > & evalFields = scatterJacobian.evaluatedFields();
+      TEST_EQUALITY(evalFields.size(),1); // this is a dummy holder for the sake of the field manager
+
+      const std::vector<RCP<PHX::FieldTag> > & fields = scatterJacobian.dependentFields();
+      TEST_EQUALITY(fields.size(),2); // these store the residual values
+
+      TEST_EQUALITY(fields[0]->name(),"Residual_ux");
+      TEST_EQUALITY(fields[1]->name(),"Residual_p");
+
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[0]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(0),Teuchos::as<int>(numCells));
+      TEST_EQUALITY(fields[1]->dataLayout().extent_int(1),Teuchos::as<int>(4)); // for Q1
+   }
+}
+
+}
diff --git a/packages/panzer/dof-mgr/test/dofmngr_test/CMakeLists.txt b/packages/panzer/dof-mgr/test/dofmngr_test/CMakeLists.txt
index eeb7aea0e3dc..cf316254c578 100644
--- a/packages/panzer/dof-mgr/test/dofmngr_test/CMakeLists.txt
+++ b/packages/panzer/dof-mgr/test/dofmngr_test/CMakeLists.txt
@@ -40,7 +40,6 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST(
 #  COMM serial mpi
 #  )
 
-IF(PANZER_HAVE_EPETRA)
 TRIBITS_ADD_EXECUTABLE_AND_TEST(
   tUniqueGlobalIndexerUtilities
   SOURCES tGlobalIndexerUtilities.cpp UnitTest_GlobalIndexer.cpp ${UNIT_TEST_DRIVER}
@@ -54,7 +53,6 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST(
   NUM_MPI_PROCS 2
   COMM mpi
   )
-ENDIF()
 
 #TRIBITS_ADD_EXECUTABLE_AND_TEST(
 #  tBlockedDOFManager_fei
diff --git a/packages/panzer/dof-mgr/test/dofmngr_test/tBlockedDOFManager.cpp b/packages/panzer/dof-mgr/test/dofmngr_test/tBlockedDOFManager.cpp
index f7dace4a2a40..a23179877179 100644
--- a/packages/panzer/dof-mgr/test/dofmngr_test/tBlockedDOFManager.cpp
+++ b/packages/panzer/dof-mgr/test/dofmngr_test/tBlockedDOFManager.cpp
@@ -44,6 +44,7 @@
 #include <Teuchos_UnitTestHarness.hpp>
 #include <Teuchos_RCP.hpp>
 #include <Teuchos_TimeMonitor.hpp>
+#include <Teuchos_DefaultComm.hpp>
 
 #include <string>
 #include <iostream>
@@ -62,9 +63,6 @@
 
 #include "Kokkos_DynRankView.hpp"
 
-#include "Epetra_MpiComm.h"
-#include "Epetra_SerialComm.h"
-
 using Teuchos::rcp;
 using Teuchos::rcp_dynamic_cast;
 using Teuchos::RCP;
@@ -86,22 +84,11 @@ Teuchos::RCP<const panzer::FieldPattern> buildFieldPattern()
 // this just excercises a bunch of functions
 TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,assortedTests)
 {
-
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
-   #endif
-
-   // panzer::pauseToAttach();
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   using Teuchos::RCP;
-   using Teuchos::rcp;
-   using Teuchos::rcp_dynamic_cast;
-
-   int myRank = eComm->MyPID();
-   int numProc = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProc = comm->getSize();
 
    RCP<ConnManager> connManager = rcp(new unit_test::ConnManager(myRank,numProc));
    BlockedDOFManager dofManager; 
@@ -134,25 +121,15 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,assortedTests)
    TEST_ASSERT(dofManager.getElementBlock("block_0")==connManager->getElementBlock("block_0"));
    TEST_ASSERT(dofManager.getElementBlock("block_1")==connManager->getElementBlock("block_1"));
    TEST_ASSERT(dofManager.getElementBlock("block_2")==connManager->getElementBlock("block_2"));
-
 }
 
 TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,registerFields)
 {
-
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
-   #endif
-
-   using Teuchos::RCP;
-   using Teuchos::rcp;
-   using Teuchos::rcp_dynamic_cast;
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   int myRank = eComm->MyPID();
-   int numProc = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProc = comm->getSize();
 
    RCP<ConnManager> connManger = rcp(new unit_test::ConnManager(myRank,numProc));
    BlockedDOFManager dofManager; 
@@ -264,29 +241,17 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,registerFields)
 
    TEST_EQUALITY(blk2fn[0],4);
    TEST_EQUALITY(blk2fn[1],5);
-
 }
 
 /* NOTE: disabled as Blocked manager doesn't handle gids anymore.
 
 TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,buildGlobalUnknowns)
 {
-
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
-   #endif
-
-   // panzer::pauseToAttach();
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   using Teuchos::RCP;
-   using Teuchos::rcp;
-   using Teuchos::rcp_dynamic_cast;
-
-   int myRank = eComm->MyPID();
-   int numProc = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProc = comm->getSize();
 
    RCP<ConnManager> connManger = rcp(new unit_test::ConnManager(myRank,numProc));
    BlockedDOFManager dofManager; 
@@ -389,22 +354,11 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,buildGlobalUnknowns)
 
 TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
 {
-
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
-   #endif
-
-   // panzer::pauseToAttach();
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   using Teuchos::RCP;
-   using Teuchos::rcp;
-   using Teuchos::rcp_dynamic_cast;
-
-   int myRank = eComm->MyPID();
-   int numProc = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProc = comm->getSize();
 
    RCP<ConnManager> connManger = rcp(new unit_test::ConnManager(myRank,numProc));
    BlockedDOFManager dofManager; 
@@ -433,13 +387,11 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
 
    dofManager.buildGlobalUnknowns();
 
-   /* NOTE: disabled since BlockManager can't do this anymore...
-
    // check from element block 0
    std::vector<std::pair<int,int> > gids;
    if(myRank==0) {
       // std::vector<std::pair<int,int> > gids;
-      dofManager.getElementGIDs(1,gids);
+      dofManager.getElementGIDsPair(1,gids);
 
       TEST_EQUALITY(gids.size(),16);
       for(std::size_t i=0; i<8;i++) TEST_EQUALITY(gids[i].first,0);
@@ -448,7 +400,7 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
    }
    else if(myRank==1) {
       // std::vector<std::pair<int,int> > gids;
-      dofManager.getElementGIDs(1,gids);
+      dofManager.getElementGIDsPair(1,gids);
 
       TEST_EQUALITY(gids.size(),16);
       for(std::size_t i=0; i<8;i++) TEST_EQUALITY(gids[i].first,0);
@@ -459,7 +411,7 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
    // check from element block 1 
    if(myRank==0) {
       // std::vector<std::pair<int,int> > gids;
-      dofManager.getElementGIDs(4,gids);
+      dofManager.getElementGIDsPair(4,gids);
 
       TEST_EQUALITY(gids.size(),4);
       for(std::size_t i=0;i<4;i++)
@@ -467,7 +419,7 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
    }
    else if(myRank==1) {
       // std::vector<std::pair<int,int> > gids;
-      dofManager.getElementGIDs(3,gids);
+      dofManager.getElementGIDsPair(3,gids);
 
       TEST_EQUALITY(gids.size(),4);
       for(std::size_t i=0;i<4;i++)
@@ -477,7 +429,7 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
    // check from element block 2 
    if(myRank==0) {
       // std::vector<std::pair<int,int> > gids;
-      dofManager.getElementGIDs(3,gids);
+      dofManager.getElementGIDsPair(3,gids);
 
       TEST_EQUALITY(gids.size(),8);
       for(std::size_t i=0;i<4;i++) {
@@ -563,12 +515,10 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,getElement_gids_fieldoffsets)
       for(std::size_t i=0;i<vec->first.size();i++) 
          TEST_EQUALITY(vec->first[i],dofManager.getBlockGIDOffset("block_2",2)+sub_vec->first[i]);
    }
-   */
 }
 
 TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,validFieldOrder)
 {
-
    BlockedDOFManager dofManager; 
    dofManager.setUseDOFManagerFEI(false);
 
@@ -658,25 +608,15 @@ TEUCHOS_UNIT_TEST(tBlockedDOFManager_SimpleTests,validFieldOrder)
  
       TEST_ASSERT(!dofManager.validFieldOrder(order,validFields));
    }
-
 }
 
 TEUCHOS_UNIT_TEST(tBlockedDOFManager,mergetests)
 {
-
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      Teuchos::RCP<Epetra_Comm> eComm = Teuchos::rcp(new Epetra_SerialComm());
-   #endif
-
-   using Teuchos::RCP;
-   using Teuchos::rcp;
-   using Teuchos::rcp_dynamic_cast;
-
-   int myRank = eComm->MyPID();
-   int numProc = eComm->NumProc();
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
+
+   int myRank = comm->getRank();
+   int numProc = comm->getSize();
 
    RCP<ConnManager> connManager = rcp(new unit_test::ConnManager(myRank,numProc));
 
diff --git a/packages/panzer/dof-mgr/test/dofmngr_test/tFilteredUGI.cpp b/packages/panzer/dof-mgr/test/dofmngr_test/tFilteredUGI.cpp
index 337d2f3f0c3f..6b811062dd27 100644
--- a/packages/panzer/dof-mgr/test/dofmngr_test/tFilteredUGI.cpp
+++ b/packages/panzer/dof-mgr/test/dofmngr_test/tFilteredUGI.cpp
@@ -65,9 +65,6 @@
 
 #include "Kokkos_DynRankView.hpp"
 
-// #include "Epetra_MpiComm.h"
-// #include "Epetra_SerialComm.h"
-
 using Teuchos::rcp;
 using Teuchos::rcp_dynamic_cast;
 using Teuchos::RCP;
diff --git a/packages/panzer/dof-mgr/test/dofmngr_test/tGlobalIndexerUtilities.cpp b/packages/panzer/dof-mgr/test/dofmngr_test/tGlobalIndexerUtilities.cpp
index 81ab33363025..e383e3c2fa7b 100644
--- a/packages/panzer/dof-mgr/test/dofmngr_test/tGlobalIndexerUtilities.cpp
+++ b/packages/panzer/dof-mgr/test/dofmngr_test/tGlobalIndexerUtilities.cpp
@@ -45,6 +45,7 @@
 #include <Teuchos_RCP.hpp>
 #include <Teuchos_TimeMonitor.hpp>
 #include <Teuchos_ParameterList.hpp>
+#include <Teuchos_DefaultComm.hpp>
 
 #include <string>
 #include <iostream>
@@ -57,13 +58,6 @@
 
 #include "UnitTest_GlobalIndexer.hpp"
 
-#ifdef HAVE_MPI
-   #include "Epetra_MpiComm.h"
-   #include "mpi.h"
-#else
-   #include "Epetra_SerialComm.h"
-#endif
-
 #include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
 
 using Teuchos::rcp;
@@ -75,16 +69,11 @@ namespace panzer {
 
 TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,GhostedFieldVector)
 {
-
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
-   #endif
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   int myRank = eComm->MyPID();
-   int numProcs = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProcs = comm->getSize();
 
    TEUCHOS_ASSERT(numProcs==2);
 
@@ -179,18 +168,13 @@ void fillFieldContainer(int fieldNum,const std::string & blockId,
 
 TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,updateGhostedDataVector)
 {
-
    typedef Kokkos::DynRankView<int,PHX::Device> IntFieldContainer;
 
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
-   #endif
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   int myRank = eComm->MyPID();
-   int numProcs = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProcs = comm->getSize();
 
    TEUCHOS_ASSERT(numProcs==2);
 
@@ -263,18 +247,13 @@ TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,updateGhostedDataVector)
 
 TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,ArrayToFieldVector_ghost)
 {
-
    typedef Kokkos::DynRankView<int,PHX::Device> IntFieldContainer;
 
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
-   #endif
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   int myRank = eComm->MyPID();
-   int numProcs = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProcs = comm->getSize();
 
    TEUCHOS_ASSERT(numProcs==2);
 
@@ -341,18 +320,13 @@ TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,ArrayToFieldVector_ghost)
 
 TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,ArrayToFieldVector)
 {
-
    typedef Kokkos::DynRankView<int,PHX::Device> IntFieldContainer;
 
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
-   #endif
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
 
-   int myRank = eComm->MyPID();
-   int numProcs = eComm->NumProc();
+   int myRank = comm->getRank();
+   int numProcs = comm->getSize();
 
    TEUCHOS_ASSERT(numProcs==2);
 
@@ -409,18 +383,13 @@ TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,ArrayToFieldVector)
 
 TEUCHOS_UNIT_TEST(tGlobalIndexer_Utilities,ArrayToFieldVector_multicol)
 {
-
    typedef Kokkos::DynRankView<int,PHX::Device> IntFieldContainer;
 
    // build global (or serial communicator)
-   #ifdef HAVE_MPI
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_MpiComm(MPI_COMM_WORLD));
-   #else
-      RCP<Epetra_Comm> eComm = rcp(new Epetra_SerialComm());
-   #endif
-
-   int myRank = eComm->MyPID();
-   int numProcs = eComm->NumProc();
+   RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
+
+   int myRank = comm->getRank();
+   int numProcs = comm->getSize();
 
    TEUCHOS_ASSERT(numProcs==2);