diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt b/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt
index 3b2202563a32..d871d0375cb0 100644
--- a/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt
+++ b/packages/panzer/adapters-stk/test/evaluator_tests/CMakeLists.txt
@@ -43,6 +43,13 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST(
   NUM_MPI_PROCS 2
   )
 
+TRIBITS_ADD_EXECUTABLE_AND_TEST(
+  tGatherSolution_BlockedTpetra
+  SOURCES tpetra_blocked_gather_solution.cpp ${UNIT_TEST_DRIVER}
+  COMM serial mpi
+  NUM_MPI_PROCS 2
+  )
+
 TRIBITS_ADD_EXECUTABLE_AND_TEST(
   tScatterResidual_Tpetra
   SOURCES tpetra_scatter_residual.cpp ${UNIT_TEST_DRIVER}
diff --git a/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_gather_solution.cpp b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_gather_solution.cpp
new file mode 100644
index 000000000000..279956f8d6eb
--- /dev/null
+++ b/packages/panzer/adapters-stk/test/evaluator_tests/tpetra_blocked_gather_solution.cpp
@@ -0,0 +1,721 @@
+// @HEADER
+// *****************************************************************************
+//           Panzer: A partial differential equation assembly
+//       engine for strongly coupled complex multiphysics systems
+//
+// Copyright 2011 NTESS and the Panzer contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+///////////////////////////////////////////////////////////////////////////////
+//
+//  Include Files
+//
+///////////////////////////////////////////////////////////////////////////////
+
+// C++
+#include <cstdio>
+#include <string>
+#include <vector>
+
+// Kokkos
+#include "Kokkos_View_Fad.hpp"
+
+// Panzer
+#include "PanzerAdaptersSTK_config.hpp"
+#include "Panzer_BasisIRLayout.hpp"
+#include "Panzer_BlockedTpetraLinearObjFactory.hpp"
+#include "Panzer_BlockedDOFManager.hpp"
+#include "Panzer_DOFManager.hpp"
+#include "Panzer_Evaluator_WithBaseImpl.hpp"
+#include "Panzer_FieldManagerBuilder.hpp"
+#include "Panzer_GatherOrientation.hpp"
+#include "Panzer_PureBasis.hpp"
+#include "Panzer_STKConnManager.hpp"
+#include "Panzer_STK_Interface.hpp"
+#include "Panzer_STK_SetupUtilities.hpp"
+#include "Panzer_STK_SquareQuadMeshFactory.hpp"
+#include "Panzer_STK_Version.hpp"
+#include "Panzer_Workset.hpp"
+#include "Panzer_LOCPair_GlobalEvaluationData.hpp"
+#include "Panzer_GlobalEvaluationDataContainer.hpp"
+
+// Teuchos
+#include "Teuchos_DefaultMpiComm.hpp"
+#include "Teuchos_GlobalMPISession.hpp"
+#include "Teuchos_OpaqueWrapper.hpp"
+#include "Teuchos_RCP.hpp"
+#include "Teuchos_TimeMonitor.hpp"
+#include "Teuchos_UnitTestHarness.hpp"
+
+// Thyra
+#include "Thyra_ProductVectorBase.hpp"
+#include "Thyra_VectorStdOps.hpp"
+
+// Tpetra
+#include "Tpetra_Vector.hpp"
+
+// user_app
+#include "user_app_EquationSetFactory.hpp"
+
+typedef double ScalarT;
+using LocalOrdinalT = panzer::LocalOrdinal;
+using GlobalOrdinalT = panzer::GlobalOrdinal;
+
+typedef Tpetra::Vector<ScalarT, LocalOrdinalT, GlobalOrdinalT> VectorType;
+typedef Tpetra::Operator<ScalarT, LocalOrdinalT, GlobalOrdinalT> OperatorType;
+typedef Tpetra::CrsMatrix<ScalarT, LocalOrdinalT, GlobalOrdinalT> CrsMatrixType;
+typedef Tpetra::CrsGraph<LocalOrdinalT, GlobalOrdinalT> CrsGraphType;
+typedef Tpetra::Map<LocalOrdinalT, GlobalOrdinalT> MapType;
+typedef Tpetra::Import<LocalOrdinalT, GlobalOrdinalT> ImportType;
+typedef Tpetra::Export<LocalOrdinalT, GlobalOrdinalT> ExportType;
+
+typedef Thyra::TpetraLinearOp<ScalarT, LocalOrdinalT, GlobalOrdinalT> ThyraLinearOp;
+
+typedef panzer::BlockedTpetraLinearObjFactory<panzer::Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT> BlockedTpetraLinObjFactoryType;
+typedef panzer::TpetraLinearObjFactory<panzer::Traits, ScalarT, LocalOrdinalT, GlobalOrdinalT> TpetraLinObjFactoryType;
+typedef panzer::BlockedTpetraLinearObjContainer<ScalarT, LocalOrdinalT, GlobalOrdinalT> BlockedTpetraLinObjContainerType;
+typedef panzer::TpetraLinearObjContainer<ScalarT, LocalOrdinalT, GlobalOrdinalT> TpetraLinObjContainerType;
+
+namespace panzer
+{
+
+  Teuchos::RCP<panzer::PureBasis> buildBasis(std::size_t worksetSize, const std::string &basisName);
+  void testInitialization(const Teuchos::RCP<Teuchos::ParameterList> &ipb);
+  Teuchos::RCP<panzer_stk::STK_Interface> buildMesh(int elemX, int elemY);
+  void testGatherScatter(const bool enable_tangents, Teuchos::FancyOStream &out, bool &success);
+
+  // Test without tangent fields in gather evaluator
+  TEUCHOS_UNIT_TEST(tpetra_assembly, gather_solution_no_tangents)
+  {
+    testGatherScatter(false, out, success);
+  }
+
+  // Test with tangent fields in gather evaluator
+  TEUCHOS_UNIT_TEST(tpetra_assembly, gather_solution_tangents)
+  {
+    testGatherScatter(true, out, success);
+  }
+
+  // enable_tangents determines whether tangent fields dx/dp are added to gather evaluator.
+  // These are used when computing df/dx*dx/dp with the tangent evaluation type
+  void testGatherScatter(const bool enable_tangents, Teuchos::FancyOStream &out, bool &success)
+  {
+#ifdef HAVE_MPI
+    Teuchos::RCP<const Teuchos::MpiComm<int>> tComm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));
+#else
+    Teuchos::RCP<const Teuchos::SerialComm<int>> tComm = Teuchos::rcp(new Teuchos::SerialComm<int>(MPI_COMM_WORLD));
+#endif
+
+    int myRank = tComm->getRank();
+    int numProcs = tComm->getSize();
+
+    const std::size_t workset_size = 4 / numProcs;
+    const std::string fieldName1_q1 = "U";
+    const std::string fieldName2_q1 = "V";
+    const std::string fieldName_qedge1 = "B";
+    const int num_tangent = enable_tangents ? 5 : 0;
+
+    Teuchos::RCP<panzer_stk::STK_Interface> mesh = buildMesh(2, 2);
+
+    // build input physics block
+    Teuchos::RCP<panzer::PureBasis> basis_q1 = buildBasis(workset_size, "Q1");
+    Teuchos::RCP<panzer::PureBasis> basis_qedge1 = buildBasis(workset_size, "QEdge1");
+
+    Teuchos::RCP<Teuchos::ParameterList> ipb = Teuchos::parameterList();
+    testInitialization(ipb);
+
+    const int default_int_order = 1;
+    std::string eBlockID = "eblock-0_0";
+    Teuchos::RCP<user_app::MyFactory> eqset_factory = Teuchos::rcp(new user_app::MyFactory);
+    panzer::CellData cellData(workset_size, mesh->getCellTopology("eblock-0_0"));
+    Teuchos::RCP<panzer::GlobalData> gd = panzer::createGlobalData();
+    Teuchos::RCP<panzer::PhysicsBlock> physicsBlock =
+        Teuchos::rcp(new PhysicsBlock(ipb, eBlockID, default_int_order, cellData, eqset_factory, gd, false));
+
+    Teuchos::RCP<std::vector<panzer::Workset>> work_sets = panzer_stk::buildWorksets(*mesh, physicsBlock->elementBlockID(),
+                                                                                     physicsBlock->getWorksetNeeds());
+    TEST_EQUALITY(work_sets->size(), 1);
+
+    // build connection manager and field manager
+    const Teuchos::RCP<panzer::ConnManager> conn_manager = Teuchos::rcp(new panzer_stk::STKConnManager(mesh));
+    Teuchos::RCP<panzer::BlockedDOFManager> blocked_dofManager = Teuchos::rcp(new panzer::BlockedDOFManager(conn_manager, MPI_COMM_WORLD));
+
+    blocked_dofManager->addField(fieldName1_q1, Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis())));
+    blocked_dofManager->addField(fieldName2_q1, Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis())));
+    blocked_dofManager->addField(fieldName_qedge1, Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_qedge1->getIntrepid2Basis())));
+
+    std::vector<std::vector<std::string> > fieldOrder(3);
+    fieldOrder[0].push_back(fieldName1_q1);
+    fieldOrder[1].push_back(fieldName_qedge1);
+    fieldOrder[2].push_back(fieldName2_q1);
+    blocked_dofManager->setFieldOrder(fieldOrder);
+
+    blocked_dofManager->buildGlobalUnknowns();
+
+    // setup linear object factory
+    /////////////////////////////////////////////////////////////
+
+    Teuchos::RCP<BlockedTpetraLinObjFactoryType> t_lof = Teuchos::rcp(new BlockedTpetraLinObjFactoryType(tComm.getConst(), blocked_dofManager));
+    Teuchos::RCP<LinearObjFactory<panzer::Traits>> lof = t_lof;
+    Teuchos::RCP<LinearObjContainer> loc = t_lof->buildGhostedLinearObjContainer();
+    t_lof->initializeGhostedContainer(LinearObjContainer::X, *loc);
+    loc->initialize();
+
+    Teuchos::RCP<BlockedTpetraLinObjContainerType> t_loc = Teuchos::rcp_dynamic_cast<BlockedTpetraLinObjContainerType>(loc);
+    Teuchos::RCP<Thyra::VectorBase<double>> x_vec = t_loc->get_x_th();
+    Thyra::assign(x_vec.ptr(), 123.0 + myRank);
+
+    // need a place to evaluate the tangent fields, so we create a 
+    // unblocked DOFManager and LOF and set up if needed
+    std::vector<Teuchos::RCP<GlobalEvaluationData>> tangentContainers;
+    Teuchos::RCP<panzer::DOFManager> dofManager = Teuchos::rcp(new panzer::DOFManager(conn_manager, MPI_COMM_WORLD));
+    Teuchos::RCP<TpetraLinObjFactoryType> tangent_lof = Teuchos::rcp(new TpetraLinObjFactoryType(tComm.getConst(), dofManager));
+    Teuchos::RCP<LinearObjFactory<panzer::Traits>> parent_tangent_lof = tangent_lof;
+ 
+    if (enable_tangents)
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      using Teuchos::rcp_dynamic_cast;
+      using Thyra::ProductVectorBase;
+      using LOCPair = panzer::LOCPair_GlobalEvaluationData;
+
+      std::vector<std::string> tangent_fieldOrder;
+      for (int i(0); i < num_tangent; ++i) {
+        std::stringstream ssedge;
+        ssedge << fieldName_qedge1 << " Tangent " << i;
+        std::stringstream ss1, ss2;
+        ss1 << fieldName1_q1 << " Tangent " << i;
+        ss2 << fieldName2_q1 << " Tangent " << i;
+ 
+        dofManager->addField(ss1.str(), Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis())));
+        dofManager->addField(ss2.str(), Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_q1->getIntrepid2Basis())));
+        dofManager->addField(ssedge.str(), Teuchos::rcp(new panzer::Intrepid2FieldPattern(basis_qedge1->getIntrepid2Basis())));
+        tangent_fieldOrder.push_back(ss1.str());
+        tangent_fieldOrder.push_back(ss2.str());
+        tangent_fieldOrder.push_back(ssedge.str());
+      }
+      dofManager->setFieldOrder(tangent_fieldOrder);
+      dofManager->buildGlobalUnknowns();
+
+      // generate and evaluate some fields
+      Teuchos::RCP<LinearObjContainer> tangent_loc = tangent_lof->buildGhostedLinearObjContainer();
+      tangent_lof->initializeGhostedContainer(LinearObjContainer::X, *tangent_loc);
+      tangent_loc->initialize();
+
+      for (int i(0); i < num_tangent; ++i)
+      {
+        auto locPair = Teuchos::rcp(new LOCPair(tangent_lof, panzer::LinearObjContainer::X));
+
+        auto global_t_loc = rcp_dynamic_cast<TpetraLinObjContainerType>(locPair->getGlobalLOC());
+        Teuchos::RCP<Thyra::VectorBase<double>> global_x_vec = global_t_loc->get_x_th();
+        Thyra::assign(global_x_vec.ptr(), 0.123 + myRank + i);
+
+        auto ghosted_t_loc = rcp_dynamic_cast<TpetraLinObjContainerType>(locPair->getGhostedLOC());
+        Teuchos::RCP<Thyra::VectorBase<double>> ghosted_x_vec = ghosted_t_loc->get_x_th();
+        Thyra::assign(ghosted_x_vec.ptr(), 0.123 + myRank + i);
+
+        tangentContainers.push_back(locPair);
+      } // end loop over the tangents
+    }   // end if (enable_tangents)
+
+    // setup field manager, add evaluator under test
+    /////////////////////////////////////////////////////////////
+
+    PHX::FieldManager<panzer::Traits> fm;
+
+    std::vector<PHX::index_size_type> derivative_dimensions;
+    derivative_dimensions.push_back(12);
+    fm.setKokkosExtendedDataTypeDimensions<panzer::Traits::Jacobian>(derivative_dimensions);
+
+    std::vector<PHX::index_size_type> tan_derivative_dimensions;
+    if (enable_tangents)
+      tan_derivative_dimensions.push_back(num_tangent);
+    else
+      tan_derivative_dimensions.push_back(0);
+    fm.setKokkosExtendedDataTypeDimensions<panzer::Traits::Tangent>(tan_derivative_dimensions);
+
+    Teuchos::RCP<PHX::FieldTag> evalField_q1, evalField_qedge1;
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+      names->push_back(fieldName1_q1);
+      names->push_back(fieldName2_q1);
+
+      Teuchos::ParameterList pl;
+      pl.set("Basis", basis_q1);
+      pl.set("DOF Names", names);
+      pl.set("Indexer Names", names);
+
+      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Residual>(pl);
+
+      TEST_EQUALITY(evaluator->evaluatedFields().size(), 2);
+
+      fm.registerEvaluator<panzer::Traits::Residual>(evaluator);
+      fm.requireField<panzer::Traits::Residual>(*evaluator->evaluatedFields()[0]);
+    }
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+      names->push_back(fieldName_qedge1);
+
+      Teuchos::ParameterList pl;
+      pl.set("Basis", basis_qedge1);
+      pl.set("DOF Names", names);
+      pl.set("Indexer Names", names);
+
+      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Residual>(pl);
+
+      TEST_EQUALITY(evaluator->evaluatedFields().size(), 1);
+
+      fm.registerEvaluator<panzer::Traits::Residual>(evaluator);
+      fm.requireField<panzer::Traits::Residual>(*evaluator->evaluatedFields()[0]);
+    }
+
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+      names->push_back(fieldName1_q1);
+      names->push_back(fieldName2_q1);
+
+      Teuchos::ParameterList pl;
+      pl.set("Basis", basis_q1);
+      pl.set("DOF Names", names);
+      pl.set("Indexer Names", names);
+
+      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Jacobian>(pl);
+
+      TEST_EQUALITY(evaluator->evaluatedFields().size(), 2);
+
+      fm.registerEvaluator<panzer::Traits::Jacobian>(evaluator);
+      fm.requireField<panzer::Traits::Jacobian>(*evaluator->evaluatedFields()[0]);
+    }
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+      names->push_back(fieldName_qedge1);
+
+      Teuchos::ParameterList pl;
+      pl.set("Basis", basis_qedge1);
+      pl.set("DOF Names", names);
+      pl.set("Indexer Names", names);
+
+      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Jacobian>(pl);
+
+      TEST_EQUALITY(evaluator->evaluatedFields().size(), 1);
+
+      fm.registerEvaluator<panzer::Traits::Jacobian>(evaluator);
+      fm.requireField<panzer::Traits::Jacobian>(*evaluator->evaluatedFields()[0]);
+    }
+
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+      names->push_back(fieldName1_q1);
+      names->push_back(fieldName2_q1);
+
+      Teuchos::ParameterList pl;
+      pl.set("Basis", basis_q1);
+      pl.set("DOF Names", names);
+      pl.set("Indexer Names", names);
+
+      if (enable_tangents)
+      {
+        RCP<std::vector<std::vector<std::string>>> tangent_names =
+            rcp(new std::vector<std::vector<std::string>>(2));
+        for (int i = 0; i < num_tangent; ++i)
+        {
+          std::stringstream ss1, ss2;
+          ss1 << fieldName1_q1 << " Tangent " << i;
+          ss2 << fieldName2_q1 << " Tangent " << i;
+          (*tangent_names)[0].push_back(ss1.str());
+          (*tangent_names)[1].push_back(ss2.str());
+        }
+        pl.set("Tangent Names", tangent_names);
+      }
+
+      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Tangent>(pl);
+
+      TEST_EQUALITY(evaluator->evaluatedFields().size(), 2);
+
+      fm.registerEvaluator<panzer::Traits::Tangent>(evaluator);
+      fm.requireField<panzer::Traits::Tangent>(*evaluator->evaluatedFields()[0]);
+    }
+    {
+      using Teuchos::RCP;
+      using Teuchos::rcp;
+      RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+      names->push_back(fieldName_qedge1);
+
+      Teuchos::ParameterList pl;
+      pl.set("Basis", basis_qedge1);
+      pl.set("DOF Names", names);
+      pl.set("Indexer Names", names);
+
+      if (enable_tangents)
+      {
+        RCP<std::vector<std::vector<std::string>>> tangent_names =
+            rcp(new std::vector<std::vector<std::string>>(1));
+        for (int i = 0; i < num_tangent; ++i)
+        {
+          std::stringstream ss;
+          ss << fieldName_qedge1 << " Tangent " << i;
+          (*tangent_names)[0].push_back(ss.str());
+        }
+        pl.set("Tangent Names", tangent_names);
+      }
+
+      Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator = lof->buildGather<panzer::Traits::Tangent>(pl);
+
+      TEST_EQUALITY(evaluator->evaluatedFields().size(), 1);
+
+      fm.registerEvaluator<panzer::Traits::Tangent>(evaluator);
+      fm.requireField<panzer::Traits::Tangent>(*evaluator->evaluatedFields()[0]);
+    }
+
+    if (enable_tangents)
+    {
+      for (int i = 0; i < num_tangent; ++i)
+      {
+        using Teuchos::RCP;
+        using Teuchos::rcp;
+        RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+        RCP<std::vector<std::string>> tangent_names = rcp(new std::vector<std::string>);
+        names->push_back(fieldName1_q1);
+        names->push_back(fieldName2_q1);
+        {
+          std::stringstream ss1, ss2;
+          ss1 << fieldName1_q1 << " Tangent " << i;
+          ss2 << fieldName2_q1 << " Tangent " << i;
+          tangent_names->push_back(ss1.str());
+          tangent_names->push_back(ss2.str());
+        }
+
+        Teuchos::ParameterList pl;
+        pl.set("Basis", basis_q1);
+        pl.set("DOF Names", tangent_names);
+        pl.set("Indexer Names", tangent_names);
+
+        {
+          std::stringstream ss;
+          ss << "Tangent Container " << i;
+          pl.set("Global Data Key", ss.str());
+        }
+
+        Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator =
+            parent_tangent_lof->buildGatherTangent<panzer::Traits::Tangent>(pl);
+
+        TEST_EQUALITY(evaluator->evaluatedFields().size(), 2);
+
+        fm.registerEvaluator<panzer::Traits::Tangent>(evaluator);
+      }
+      for (int i = 0; i < num_tangent; ++i)
+      {
+        using Teuchos::RCP;
+        using Teuchos::rcp;
+        RCP<std::vector<std::string>> names = rcp(new std::vector<std::string>);
+        RCP<std::vector<std::string>> tangent_names = rcp(new std::vector<std::string>);
+        names->push_back(fieldName_qedge1);
+        {
+          std::stringstream ss;
+          ss << fieldName_qedge1 << " Tangent " << i;
+          tangent_names->push_back(ss.str());
+        }
+
+        Teuchos::ParameterList pl;
+        pl.set("Basis", basis_qedge1);
+        pl.set("DOF Names", tangent_names);
+        pl.set("Indexer Names", tangent_names);
+
+        {
+          std::stringstream ss;
+          ss << "Tangent Container " << i;
+          pl.set("Global Data Key", ss.str());
+        }
+
+        Teuchos::RCP<PHX::Evaluator<panzer::Traits>> evaluator =
+            parent_tangent_lof->buildGatherTangent<panzer::Traits::Tangent>(pl);
+
+        TEST_EQUALITY(evaluator->evaluatedFields().size(), 1);
+
+        fm.registerEvaluator<panzer::Traits::Tangent>(evaluator);
+      }
+    }
+
+    panzer::Traits::SD sd;
+
+    panzer::Workset &workset = (*work_sets)[0];
+    workset.alpha = 0.0;
+    workset.beta = 2.0; // derivatives multiplied by 2
+    workset.time = 0.0;
+    workset.evaluate_transient_terms = false;
+
+    sd.worksets_ = work_sets;
+
+    fm.postRegistrationSetup(sd);
+
+    panzer::Traits::PED ped;
+    ped.gedc->addDataObject("Solution Gather Container", loc);
+    if (enable_tangents)
+    {
+      for (int i(0); i < num_tangent; ++i)
+      {
+        std::stringstream ss;
+        ss << "Tangent Container " << i;
+        ped.gedc->addDataObject(ss.str(), tangentContainers[i]);
+      }
+    }
+
+    fm.preEvaluate<panzer::Traits::Residual>(ped);
+    fm.evaluateFields<panzer::Traits::Residual>(workset);
+    fm.postEvaluate<panzer::Traits::Residual>(0);
+
+    fm.preEvaluate<panzer::Traits::Jacobian>(ped);
+    fm.evaluateFields<panzer::Traits::Jacobian>(workset);
+    fm.postEvaluate<panzer::Traits::Jacobian>(0);
+
+    fm.preEvaluate<panzer::Traits::Tangent>(ped);
+    fm.evaluateFields<panzer::Traits::Tangent>(workset);
+    fm.postEvaluate<panzer::Traits::Tangent>(0);
+
+    // test Residual fields
+    {
+      PHX::MDField<panzer::Traits::Residual::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData1_q1(fieldName1_q1, basis_q1->functional);
+      PHX::MDField<panzer::Traits::Residual::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData2_q1(fieldName2_q1, basis_qedge1->functional);
+
+      fm.getFieldData<panzer::Traits::Residual>(fieldData1_q1);
+      fm.getFieldData<panzer::Traits::Residual>(fieldData2_q1);
+
+      TEST_EQUALITY(fieldData1_q1.extent(0), Teuchos::as<unsigned int>(4 / numProcs));
+      TEST_EQUALITY(fieldData1_q1.extent(1), 4);
+      TEST_EQUALITY(fieldData2_q1.extent(0), Teuchos::as<unsigned int>(4 / numProcs));
+      TEST_EQUALITY(fieldData2_q1.extent(1), 4);
+      TEST_EQUALITY(fieldData1_q1.size(), Teuchos::as<unsigned int>(4 * 4 / numProcs));
+      TEST_EQUALITY(fieldData2_q1.size(), Teuchos::as<unsigned int>(4 * 4 / numProcs));
+
+      auto fieldData1_q1_h = Kokkos::create_mirror_view(fieldData1_q1.get_static_view());
+      auto fieldData2_q1_h = Kokkos::create_mirror_view(fieldData2_q1.get_static_view());
+      Kokkos::deep_copy(fieldData1_q1_h, fieldData1_q1.get_static_view());
+      Kokkos::deep_copy(fieldData2_q1_h, fieldData2_q1.get_static_view());
+
+      for (unsigned int i = 0; i < fieldData1_q1.extent(0); i++)
+        for (unsigned int j = 0; j < fieldData1_q1.extent(1); j++)
+          TEST_EQUALITY(fieldData1_q1_h(i, j), 123.0 + myRank);
+
+      for (unsigned int i = 0; i < fieldData2_q1.extent(0); i++)
+        for (unsigned int j = 0; j < fieldData2_q1.extent(1); j++)
+          TEST_EQUALITY(fieldData2_q1_h(i, j), 123.0 + myRank);
+    }
+    {
+      PHX::MDField<panzer::Traits::Residual::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData_qedge1(fieldName_qedge1, basis_qedge1->functional);
+
+      fm.getFieldData<panzer::Traits::Residual>(fieldData_qedge1);
+
+      auto fieldData_qedge1_h = Kokkos::create_mirror_view(fieldData_qedge1.get_static_view());
+      Kokkos::deep_copy(fieldData_qedge1_h, fieldData_qedge1.get_static_view());
+
+      TEST_EQUALITY(fieldData_qedge1.extent(0), Teuchos::as<unsigned int>(4 / numProcs));
+      TEST_EQUALITY(fieldData_qedge1.extent(1), 4);
+      TEST_EQUALITY(fieldData_qedge1.size(), Teuchos::as<unsigned int>(4 * 4 / numProcs));
+
+      for (unsigned int cell = 0; cell < fieldData_qedge1.extent(0); ++cell)
+        for (unsigned int pt = 0; pt < fieldData_qedge1.extent(1); pt++)
+          TEST_EQUALITY(fieldData_qedge1_h(cell, pt), 123.0 + myRank);
+    }
+
+    // test Jacobian fields
+    {
+      PHX::MDField<panzer::Traits::Jacobian::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData1_q1(fieldName1_q1, basis_q1->functional);
+      PHX::MDField<panzer::Traits::Jacobian::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData2_q1(fieldName2_q1, basis_qedge1->functional);
+
+      fm.getFieldData<panzer::Traits::Jacobian>(fieldData1_q1);
+      fm.getFieldData<panzer::Traits::Jacobian>(fieldData2_q1);
+
+      auto fieldData1_q1_h = Kokkos::create_mirror_view(fieldData1_q1.get_static_view());
+      auto fieldData2_q1_h = Kokkos::create_mirror_view(fieldData2_q1.get_static_view());
+      Kokkos::deep_copy(fieldData1_q1_h, fieldData1_q1.get_static_view());
+      Kokkos::deep_copy(fieldData2_q1_h, fieldData2_q1.get_static_view());
+
+      for (unsigned int cell = 0; cell < fieldData1_q1.extent(0); ++cell)
+      {
+        for (unsigned int pt = 0; pt < fieldData1_q1.extent(1); pt++)
+        {
+          TEST_EQUALITY(fieldData1_q1_h(cell, pt), 123.0 + myRank);
+          TEST_EQUALITY(fieldData1_q1_h(cell, pt).availableSize(), 12);
+        }
+      }
+      for (unsigned int cell = 0; cell < fieldData2_q1.extent(0); ++cell)
+      {
+        for (unsigned int pt = 0; pt < fieldData2_q1.extent(1); pt++)
+        {
+          TEST_EQUALITY(fieldData2_q1_h(cell, pt), 123.0 + myRank);
+          TEST_EQUALITY(fieldData2_q1_h(cell, pt).availableSize(), 12);
+        }
+      }
+    }
+    {
+      PHX::MDField<panzer::Traits::Jacobian::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData_qedge1(fieldName_qedge1, basis_qedge1->functional);
+
+      fm.getFieldData<panzer::Traits::Jacobian>(fieldData_qedge1);
+
+      auto fieldData_qedge1_h = Kokkos::create_mirror_view(fieldData_qedge1.get_static_view());
+      Kokkos::deep_copy(fieldData_qedge1_h, fieldData_qedge1.get_static_view());
+
+      for (unsigned int cell = 0; cell < fieldData_qedge1.extent(0); ++cell)
+      {
+        for (unsigned int pt = 0; pt < fieldData_qedge1.extent(1); ++pt)
+        {
+          TEST_EQUALITY(fieldData_qedge1_h(cell, pt), 123.0 + myRank);
+          TEST_EQUALITY(fieldData_qedge1_h(cell, pt).availableSize(), 12);
+        }
+      }
+    }
+
+    // test Tangent fields
+    {
+      PHX::MDField<panzer::Traits::Tangent::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData1_q1(fieldName1_q1, basis_q1->functional);
+      PHX::MDField<panzer::Traits::Tangent::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData2_q1(fieldName2_q1, basis_qedge1->functional);
+
+      fm.getFieldData<panzer::Traits::Tangent>(fieldData1_q1);
+      fm.getFieldData<panzer::Traits::Tangent>(fieldData2_q1);
+
+      auto fieldData1_q1_h = Kokkos::create_mirror_view(fieldData1_q1.get_static_view());
+      auto fieldData2_q1_h = Kokkos::create_mirror_view(fieldData2_q1.get_static_view());
+      Kokkos::deep_copy(fieldData1_q1_h, fieldData1_q1.get_static_view());
+      Kokkos::deep_copy(fieldData2_q1_h, fieldData2_q1.get_static_view());
+
+      for (unsigned int cell = 0; cell < fieldData1_q1.extent(0); ++cell)
+      {
+        for (unsigned int pt = 0; pt < fieldData1_q1.extent(1); pt++)
+        {
+          if (enable_tangents)
+          {
+            TEST_EQUALITY(fieldData1_q1_h(cell, pt).val(), 123.0 + myRank);
+            TEST_EQUALITY(fieldData1_q1_h(cell, pt).availableSize(), num_tangent);
+            for (int i = 0; i < num_tangent; ++i)
+              TEST_EQUALITY(fieldData1_q1_h(cell, pt).dx(i), 0.123 + myRank + i);
+          }
+          else
+          {
+            TEST_EQUALITY(fieldData1_q1_h(cell, pt), 123.0 + myRank);
+            TEST_EQUALITY(fieldData1_q1_h(cell, pt).availableSize(), 0);
+          }
+        }
+      }
+      for (unsigned int cell = 0; cell < fieldData2_q1.extent(0); ++cell)
+      {
+        for (unsigned int pt = 0; pt < fieldData2_q1.extent(1); pt++)
+        {
+          if (enable_tangents)
+          {
+            TEST_EQUALITY(fieldData2_q1_h(cell, pt).val(), 123.0 + myRank);
+            TEST_EQUALITY(fieldData2_q1_h(cell, pt).availableSize(), num_tangent);
+            for (int i = 0; i < num_tangent; ++i)
+            {
+              TEST_EQUALITY(fieldData2_q1_h(cell, pt).dx(i), 0.123 + myRank + i);
+              TEST_EQUALITY(fieldData2_q1_h(cell, pt).dx(i), 0.123 + myRank + i);
+            }
+          }
+          else
+          {
+            TEST_EQUALITY(fieldData2_q1_h(cell, pt), 123.0 + myRank);
+            TEST_EQUALITY(fieldData2_q1_h(cell, pt).availableSize(), 0);
+          }
+        }
+      }
+    }
+    {
+      PHX::MDField<panzer::Traits::Tangent::ScalarT, panzer::Cell, panzer::BASIS>
+          fieldData_qedge1(fieldName_qedge1, basis_qedge1->functional);
+
+      fm.getFieldData<panzer::Traits::Tangent>(fieldData_qedge1);
+
+      auto fieldData_qedge1_h = Kokkos::create_mirror_view(fieldData_qedge1.get_static_view());
+      Kokkos::deep_copy(fieldData_qedge1_h, fieldData_qedge1.get_static_view());
+
+      for (unsigned int cell = 0; cell < fieldData_qedge1.extent(0); ++cell)
+      {
+        for (unsigned int pt = 0; pt < fieldData_qedge1.extent(1); ++pt)
+        {
+          if (enable_tangents)
+          {
+            TEST_EQUALITY(fieldData_qedge1_h(cell, pt).val(), 123.0 + myRank);
+            TEST_EQUALITY(fieldData_qedge1_h(cell, pt).availableSize(), num_tangent);
+            for (int i = 0; i < num_tangent; ++i)
+              TEST_EQUALITY(fieldData_qedge1_h(cell, pt).dx(i), 0.123 + myRank + i);
+          }
+          else
+          {
+            TEST_EQUALITY(fieldData_qedge1_h(cell, pt), 123.0 + myRank);
+            TEST_EQUALITY(fieldData_qedge1_h(cell, pt).availableSize(), 0);
+          }
+        }
+      }
+    }
+  }
+
+  Teuchos::RCP<panzer::PureBasis> buildBasis(std::size_t worksetSize, const std::string &basisName)
+  {
+    Teuchos::RCP<shards::CellTopology> topo =
+        Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4>>()));
+
+    panzer::CellData cellData(worksetSize, topo);
+    return Teuchos::rcp(new panzer::PureBasis(basisName, 1, cellData));
+  }
+
+  Teuchos::RCP<panzer_stk::STK_Interface> buildMesh(int elemX, int elemY)
+  {
+    Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList);
+    pl->set("X Blocks", 1);
+    pl->set("Y Blocks", 1);
+    pl->set("X Elements", elemX);
+    pl->set("Y Elements", elemY);
+
+    panzer_stk::SquareQuadMeshFactory factory;
+    factory.setParameterList(pl);
+    Teuchos::RCP<panzer_stk::STK_Interface> mesh = factory.buildUncommitedMesh(MPI_COMM_WORLD);
+    factory.completeMeshConstruction(*mesh, MPI_COMM_WORLD);
+
+    return mesh;
+  }
+
+  void testInitialization(const Teuchos::RCP<Teuchos::ParameterList> &ipb)
+  {
+    // Physics block
+    ipb->setName("test physics");
+    {
+      Teuchos::ParameterList &p = ipb->sublist("a");
+      p.set("Type", "Energy");
+      p.set("Prefix", "");
+      p.set("Model ID", "solid");
+      p.set("Basis Type", "HGrad");
+      p.set("Basis Order", 1);
+      p.set("Integration Order", 1);
+    }
+    {
+      Teuchos::ParameterList &p = ipb->sublist("b");
+      p.set("Type", "Energy");
+      p.set("Prefix", "ION_");
+      p.set("Model ID", "solid");
+      p.set("Basis Type", "HCurl");
+      p.set("Basis Order", 1);
+      p.set("Integration Order", 1);
+    }
+  }
+
+}
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra.hpp
index aec43b41dfbc..6d9bde9d1a3b 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra.hpp
@@ -163,7 +163,7 @@ class GatherSolution_BlockedTpetra<panzer::Traits::Tangent,TRAITS,S,LO,GO,NodeT>
 public:
 
    GatherSolution_BlockedTpetra(const Teuchos::RCP<const BlockedDOFManager> & indexer)
-     : gidIndexer_(indexer) {}
+     : globalIndexer_(indexer) {}
 
    GatherSolution_BlockedTpetra(const Teuchos::RCP<const BlockedDOFManager> & indexer,
                                 const Teuchos::ParameterList& p);
@@ -176,13 +176,13 @@ class GatherSolution_BlockedTpetra<panzer::Traits::Tangent,TRAITS,S,LO,GO,NodeT>
   void evaluateFields(typename TRAITS::EvalData d);
 
   virtual Teuchos::RCP<CloneableEvaluator> clone(const Teuchos::ParameterList & pl) const
-  { return Teuchos::rcp(new GatherSolution_BlockedTpetra<panzer::Traits::Tangent,TRAITS,S,LO,GO>(gidIndexer_,pl)); }
+  { return Teuchos::rcp(new GatherSolution_BlockedTpetra<panzer::Traits::Tangent,TRAITS,S,LO,GO>(globalIndexer_,pl)); }
 
 
 private:
   typedef typename panzer::Traits::Tangent EvalT;
   typedef typename panzer::Traits::Tangent::ScalarT ScalarT;
-  //typedef typename panzer::Traits::RealType RealT;
+  typedef typename panzer::Traits::RealType RealT;
 
   typedef BlockedTpetraLinearObjContainer<S,LO,GO,NodeT> ContainerType;
   typedef Tpetra::Vector<S,LO,GO,NodeT> VectorType;
@@ -194,10 +194,14 @@ class GatherSolution_BlockedTpetra<panzer::Traits::Tangent,TRAITS,S,LO,GO,NodeT>
 
   // maps the local (field,element,basis) triplet to a global ID
   // for scattering
-  Teuchos::RCP<const BlockedDOFManager> gidIndexer_;
+  Teuchos::RCP<const BlockedDOFManager> globalIndexer_;
 
   std::vector<int> fieldIds_; // field IDs needing mapping
 
+  //! Returns the index into the Thyra ProductVector sub-block. Size
+  //! of number of fields to scatter
+  std::vector<int> productVectorBlockIndex_;
+
   std::vector< PHX::MDField<ScalarT,Cell,NODE> > gatherFields_;
 
   std::vector<std::string> indexerNames_;
@@ -206,9 +210,16 @@ class GatherSolution_BlockedTpetra<panzer::Traits::Tangent,TRAITS,S,LO,GO,NodeT>
 
   Teuchos::RCP<const BlockedTpetraLinearObjContainer<S,LO,GO,NodeT> > blockedContainer_;
 
+  //! Local indices for unknowns
+  PHX::View<LO**> worksetLIDs_;
+
+  //! Offset into the cell lids for each field. Size of number of fields to scatter.
+  std::vector<PHX::View<int*>> fieldOffsets_;
+
   // Fields for storing tangent components dx/dp of solution vector x
   bool has_tangent_fields_;
-  std::vector< std::vector< PHX::MDField<const ScalarT,Cell,NODE> > > tangentFields_;
+  std::vector< std::vector< PHX::MDField<const RealT,Cell,NODE> > > tangentFields_;
+  PHX::ViewOfViews<2,PHX::View<const RealT**>> tangentFieldsVoV_;
 
   GatherSolution_BlockedTpetra();
 };
diff --git a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp
index b0ef54fdd70b..52488585d37e 100644
--- a/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp
+++ b/packages/panzer/disc-fe/src/evaluators/Panzer_GatherSolution_BlockedTpetra_impl.hpp
@@ -8,8 +8,8 @@
 // *****************************************************************************
 // @HEADER
 
-#ifndef PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
-#define PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
+#ifndef PANZER_GATHER_SOLUTION_BLOCKED_TPETRA_IMPL_HPP
+#define PANZER_GATHER_SOLUTION_BLOCKED_TPETRA_IMPL_HPP
 
 #include "Teuchos_Assert.hpp"
 #include "Phalanx_DataLayout.hpp"
@@ -216,7 +216,7 @@ panzer::GatherSolution_BlockedTpetra<panzer::Traits::Tangent, TRAITS,S,LO,GO,Nod
 GatherSolution_BlockedTpetra(
   const Teuchos::RCP<const BlockedDOFManager> & indexer,
   const Teuchos::ParameterList& p)
-  : gidIndexer_(indexer)
+  : globalIndexer_(indexer)
   , has_tangent_fields_(false)
 {
   typedef std::vector< std::vector<std::string> > vvstring;
@@ -250,7 +250,7 @@ GatherSolution_BlockedTpetra(
       tangentFields_[fd].resize(tangent_field_names[fd].size());
       for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
         tangentFields_[fd][i] =
-          PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
+          PHX::MDField<const RealT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
         this->addDependentField(tangentFields_[fd][i]);
       }
     }
@@ -268,17 +268,60 @@ GatherSolution_BlockedTpetra(
 // **********************************************************************
 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
 void panzer::GatherSolution_BlockedTpetra<panzer::Traits::Tangent, TRAITS,S,LO,GO,NodeT>::
-postRegistrationSetup(typename TRAITS::SetupData /* d */,
+postRegistrationSetup(typename TRAITS::SetupData d,
                       PHX::FieldManager<TRAITS>& /* fm */)
 {
   TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
 
-  fieldIds_.resize(gatherFields_.size());
+  const Workset & workset_0 = (*d.worksets_)[0];
+  const std::string blockId = this->wda(workset_0).block_id;
 
+  fieldIds_.resize(gatherFields_.size());
+  fieldOffsets_.resize(gatherFields_.size());
+  productVectorBlockIndex_.resize(gatherFields_.size());
+  int maxElementBlockGIDCount = -1;
   for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
-    // get field ID from DOF manager
-    const std::string& fieldName = indexerNames_[fd];
-    fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName);
+
+    const std::string fieldName = indexerNames_[fd];
+    const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
+    productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
+    const auto& subGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
+    fieldIds_[fd] = subGlobalIndexer->getFieldNum(fieldName); // Field number in the sub-global-indexer
+
+    const std::vector<int>& offsets = subGlobalIndexer->getGIDFieldOffsets(blockId,fieldIds_[fd]);
+    fieldOffsets_[fd] = PHX::View<int*>("GatherSolution_BlockedTpetra(Tangent):fieldOffsets",offsets.size());
+    auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
+    for (std::size_t i=0; i < offsets.size(); ++i)
+      hostOffsets(i) = offsets[i];
+    Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
+    maxElementBlockGIDCount = std::max(subGlobalIndexer->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
+  }
+
+  // We will use one workset lid view for all fields, but has to be
+  // sized big enough to hold the largest elementBlockGIDCount in the
+  // ProductVector.
+  worksetLIDs_ = PHX::View<LO**>("ScatterResidual_BlockedTpetra(Tangent):worksetLIDs",
+                                                gatherFields_[0].extent(0),
+                                                maxElementBlockGIDCount);
+
+  // Set up storage for tangentFields using view of views
+  // We also need storage for the number of tangent fields associated with
+  // each gatherField
+
+  if (has_tangent_fields_) {
+
+    size_t inner_vector_max_size = 0;
+    for (std::size_t fd = 0; fd < tangentFields_.size(); ++fd)
+      inner_vector_max_size = std::max(inner_vector_max_size,tangentFields_[fd].size());
+    tangentFieldsVoV_.initialize("GatherSolution_BlockedTpetra<Tangent>::tangentFieldsVoV_",gatherFields_.size(),inner_vector_max_size);
+
+    for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
+      for (std::size_t i=0; i<tangentFields_[fd].size(); ++i) {
+        tangentFieldsVoV_.addView(tangentFields_[fd][i].get_static_view(),fd,i);
+      }
+    }
+
+    tangentFieldsVoV_.syncHostToDevice();
   }
 
   indexerNames_.clear();  // Don't need this anymore
@@ -298,75 +341,72 @@ template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
 void panzer::GatherSolution_BlockedTpetra<panzer::Traits::Tangent, TRAITS,S,LO,GO,NodeT>::
 evaluateFields(typename TRAITS::EvalData workset)
 {
-   using Teuchos::RCP;
-   using Teuchos::ArrayRCP;
-   using Teuchos::ptrFromRef;
-   using Teuchos::rcp_dynamic_cast;
-
-   using Thyra::VectorBase;
-   using Thyra::SpmdVectorBase;
-   using Thyra::ProductVectorBase;
+  using Teuchos::RCP;
+  using Teuchos::ArrayRCP;
+  using Teuchos::ptrFromRef;
+  using Teuchos::rcp_dynamic_cast;
 
-   Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
-   out.setShowProcRank(true);
-   out.setOutputToRootOnly(-1);
+  using Thyra::VectorBase;
+  using Thyra::SpmdVectorBase;
+  using Thyra::ProductVectorBase;
 
-   std::vector<std::pair<int,GO> > GIDs;
-   std::vector<LO> LIDs;
+  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
+  out.setShowProcRank(true);
+  out.setOutputToRootOnly(-1);
 
-   // for convenience pull out some objects from workset
-   std::string blockId = this->wda(workset).block_id;
-   const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
+  const PHX::View<const int*> & localCellIds = this->wda(workset).getLocalCellIDs();
 
-   Teuchos::RCP<ProductVectorBase<double> > x;
-   if (useTimeDerivativeSolutionVector_)
-     x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
-   else
-     x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
+  Teuchos::RCP<ProductVectorBase<double> > blockedSolution;
+  if (useTimeDerivativeSolutionVector_)
+    blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
+  else
+    blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
 
-   // gather operation for each cell in workset
-   for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
-      LO cellLocalId = localCellIds[worksetCellIndex];
+  // Loop over fields to gather
+  int currentWorksetLIDSubBlock = -1;
+  for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
+    // workset LIDs only change if in different sub blocks 
+    if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
+      const auto& blockIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]];
+      const std::string blockId = this->wda(workset).block_id;
+      const int num_dofs = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]]->getElementBlockGIDCount(blockId);
+      blockIndexer->getElementLIDs(localCellIds,worksetLIDs_,num_dofs); 
+      currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
+    }
 
-      gidIndexer_->getElementGIDsPair(cellLocalId,GIDs,blockId);
+    const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
+    const auto& subblockSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<RealT,LO,GO,NodeT>>(blockedSolution->getNonconstVectorBlock(blockRowIndex),true))->getTpetraVector());
+    const auto kokkosSolution = subblockSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
 
-      // caculate the local IDs for this element
-      LIDs.resize(GIDs.size());
-      for(std::size_t i=0;i<GIDs.size();i++) {
-         // used for doing local ID lookups
-         RCP<const MapType> x_map = blockedContainer_->getMapForBlock(GIDs[i].first);
+    // Class data fields for lambda capture
+    const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
+    const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
+    const PHX::View<ScalarT**> fieldValues = gatherFields_[fieldIndex].get_static_view();        
 
-         LIDs[i] = x_map->getLocalElement(GIDs[i].second);
-      }
+    if (has_tangent_fields_) { 
+      const int numTangents = tangentFields_[fieldIndex].size();
+      const auto tangentFieldsDevice = tangentFieldsVoV_.getViewDevice();
+      const auto kokkosTangents = Kokkos::subview(tangentFieldsDevice,fieldIndex,Kokkos::ALL());
+      Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {  
+        for (int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
+          const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
+	       fieldValues(cell,basis).zero();
+          fieldValues(cell,basis).val() = kokkosSolution(rowLID,0);
+          for (int i_tangent=0; i_tangent<numTangents; ++i_tangent)
+            fieldValues(cell,basis).fastAccessDx(i_tangent) = kokkosTangents(i_tangent)(cell,basis);
+        }
+      });
+    } else {
+      Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {  
+        for (int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
+          const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
+	       fieldValues(cell,basis).zero();
+          fieldValues(cell,basis) = kokkosSolution(rowLID,0);
+        }
+      });
+    }
+  }
 
-      // loop over the fields to be gathered
-      Teuchos::ArrayRCP<const double> local_x;
-      for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
-         int fieldNum = fieldIds_[fieldIndex];
-         int indexerId = gidIndexer_->getFieldBlock(fieldNum);
-
-         // grab local data for inputing
-         RCP<SpmdVectorBase<double> > block_x = rcp_dynamic_cast<SpmdVectorBase<double> >(x->getNonconstVectorBlock(indexerId));
-         block_x->getLocalData(ptrFromRef(local_x));
-
-         const std::vector<int> & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum);
-
-         // loop over basis functions and fill the fields
-         for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
-            int offset = elmtOffset[basis];
-            int lid = LIDs[offset];
-
-            if (!has_tangent_fields_)
-              (gatherFields_[fieldIndex])(worksetCellIndex,basis) = local_x[lid];
-            else {
-              (gatherFields_[fieldIndex])(worksetCellIndex,basis).val() = local_x[lid];
-              for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
-                (gatherFields_[fieldIndex])(worksetCellIndex,basis).fastAccessDx(i) =
-                  tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
-            }
-         }
-      }
-   }
 }
 
 // **********************************************************************