diff --git a/example/PDE-OPT/0ld/stefan-boltzmann/CMakeLists.txt b/example/PDE-OPT/0ld/stefan-boltzmann/CMakeLists.txt
index aaf59a309..8a7abd471 100644
--- a/example/PDE-OPT/0ld/stefan-boltzmann/CMakeLists.txt
+++ b/example/PDE-OPT/0ld/stefan-boltzmann/CMakeLists.txt
@@ -1,5 +1,5 @@
 
-IF(${PROJECT_NAME}_ENABLE_Intrepid AND 
+IF(${PROJECT_NAME}_ENABLE_Intrepid AND
    ${PROJECT_NAME}_ENABLE_Epetra   AND
    ${PROJECT_NAME}_ENABLE_Amesos   AND
    ${PROJECT_NAME}_ENABLE_Amesos2  AND
@@ -17,7 +17,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND
 #    ARGS PrintItAll
 #    NUM_MPI_PROCS 4
 #    NUM_TOTAL_CORES_USED 4
-#    CATEGORIES CONTINUOUS
 #    PASS_REGULAR_EXPRESSION "TEST PASSED"
 #    ADD_DIR_TO_NAME
 #  )
@@ -28,7 +27,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND
 #    ARGS PrintItAll
 #    NUM_MPI_PROCS 4
 #    NUM_TOTAL_CORES_USED 4
-#    CATEGORIES CONTINUOUS
 #    PASS_REGULAR_EXPRESSION "TEST PASSED"
 #    ADD_DIR_TO_NAME
 #  )
diff --git a/example/PDE-OPT/TEST/CMakeLists.txt b/example/PDE-OPT/TEST/CMakeLists.txt
index 818440948..549809721 100644
--- a/example/PDE-OPT/TEST/CMakeLists.txt
+++ b/example/PDE-OPT/TEST/CMakeLists.txt
@@ -1,6 +1,5 @@
 
-IF(${PROJECT_NAME}_ENABLE_Intrepid AND 
-   ${PROJECT_NAME}_ENABLE_Epetra )
+IF(${PROJECT_NAME}_ENABLE_Intrepid)
 
   ROL_ADD_EXECUTABLE_AND_TEST(
     test_01
@@ -92,6 +91,98 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND
 
 ENDIF()
 
+IF(${PROJECT_NAME}_ENABLE_Intrepid2)
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_01_Kokkos
+    SOURCES test_01K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_02_Kokkos
+    SOURCES test_02K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_03_Kokkos
+    SOURCES test_03K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_04_Kokkos
+    SOURCES test_04K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_05_Kokkos
+    SOURCES test_05K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_07_Kokkos
+    SOURCES test_07K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_08_Kokkos
+    SOURCES test_08K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_ADD_EXECUTABLE_AND_TEST(
+    test_09_Kokkos
+    SOURCES test_09K.cpp
+    ARGS PrintItAll
+    NUM_MPI_PROCS 1
+    NUM_TOTAL_CORES_USED 1
+    PASS_REGULAR_EXPRESSION "TEST PASSED"
+    ADD_DIR_TO_NAME
+  )
+
+  ROL_COPY_FILES_TO_BINARY_DIR(
+    TestDataCopyK
+    SOURCE_FILES
+      input_01.xml input_02.xml input_03.xml input_04.xml input_05.xml input_07.xml input_08.xml input_09.xml cube.txt L-cubes.txt
+    SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}"
+    DEST_DIR "${CMAKE_CURRENT_BINARY_DIR}"
+  )
+
+ENDIF()
+
 
 #IF(${PROJECT_NAME}_ENABLE_Sacado)
 #  ROL_ADD_EXECUTABLE_AND_TEST(
@@ -104,6 +195,3 @@ ENDIF()
 #    ADD_DIR_TO_NAME
 #  )
 #ENDIF()
-
-
-
diff --git a/example/PDE-OPT/TEST/test_01K.cpp b/example/PDE-OPT/TEST/test_01K.cpp
new file mode 100644
index 000000000..b347a1539
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_01K.cpp
@@ -0,0 +1,166 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_01.cpp
+    \brief Unit test for the mesh manager and the degree-of-freedom manager.
+           Mesh type: RECTANGLE with QUAD CELLS and HGRAD SPACE.
+*/
+
+#include "ROL_Algorithm.hpp"
+#include "ROL_BoundConstraint_SimOpt.hpp"
+#include "ROL_Vector_SimOpt.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
+#include "Intrepid2_HGRAD_QUAD_C2_FEM.hpp"
+
+#include <Kokkos_Core.hpp>
+
+#include <iostream>
+#include <algorithm>
+
+#include "../TOOLS/meshmanagerK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+
+using scalar_view = typename MeshManager<RealT, DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT, DeviceT>::int_view;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_01.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshManager_Rectangle<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+    int_view cellToEdgeMap = meshmgr.getCellToEdgeMap();
+    ROL::Ptr<std::vector<std::vector<std::vector<int> > > > sideSetsPtr = meshmgr.getSideSets();
+
+    std::vector<std::vector<std::vector<int> > >  &sideSets = *sideSetsPtr;
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl; // << nodes;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl; // << cellToNodeMap;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl; // << cellToEdgeMap;
+    // Print mesh info to file.
+    std::ofstream meshfile;
+    meshfile.open("cells.txt");
+    for (int i=0; i<cellToNodeMap.extent_int(0); ++i) {
+      meshfile << cellToNodeMap(i,0) << "  "
+               << cellToNodeMap(i,1) << "  "
+               << cellToNodeMap(i,2) << "  "
+               << cellToNodeMap(i,3) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("edges.txt");
+    for (int i=0; i<cellToEdgeMap.extent_int(0); ++i) {
+      meshfile << cellToEdgeMap(i,0) << "  "
+               << cellToEdgeMap(i,1) << "  "
+               << cellToEdgeMap(i,2) << "  "
+               << cellToEdgeMap(i,3) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("nodes.txt");
+    for (int i=0; i<nodes.extent_int(0); ++i) {
+      meshfile << nodes(i,0) << "  "
+               << nodes(i,1) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("sideset.txt");
+    for (int i=0; i<static_cast<int>(sideSets.size()); ++i) {
+      for (int j=0; j<static_cast<int>(sideSets[i].size()); ++j) {
+        if (sideSets[i][j].size() > 0) {
+          for (int k=0; k<static_cast<int>(sideSets[i][j].size()); ++k) {
+            meshfile << sideSets[i][j][k] << std::endl;
+          }
+        }
+        meshfile << std::endl << std::endl;
+      }
+    }
+    meshfile.close();
+
+    basis_ptr basisPtrQ1 = ROL::makePtr<Intrepid2::Basis_HGRAD_QUAD_C1_FEM<DeviceT, RealT, RealT>>();
+    basis_ptr basisPtrQ2 = ROL::makePtr<Intrepid2::Basis_HGRAD_QUAD_C2_FEM<DeviceT, RealT, RealT>>();
+
+    std::vector<basis_ptr> basisPtrs(3, ROL::nullPtr);
+    basisPtrs[0] = basisPtrQ2;
+    basisPtrs[1] = basisPtrQ1;
+    basisPtrs[2] = basisPtrQ2;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT> > meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    DofManager<RealT, DeviceT> dofmgr(meshmgrPtr, basisPtrs);
+
+    *outStream << "Number of node dofs = " << dofmgr.getNumNodeDofs() << std::endl; // << dofmgr.getNodeDofs();
+    *outStream << "Number of edge dofs = " << dofmgr.getNumEdgeDofs() << std::endl; // << dofmgr.getEdgeDofs();
+    *outStream << "Number of face dofs = " << dofmgr.getNumFaceDofs() << std::endl; // << dofmgr.getFaceDofs();
+    *outStream << "Number of void dofs = " << dofmgr.getNumVoidDofs() << std::endl; // << dofmgr.getVoidDofs();
+    *outStream << "Total number of dofs = " << dofmgr.getNumDofs() << std::endl; // << dofmgr.getCellDofs();
+
+    std::vector<std::vector<int>> fieldPattern = dofmgr.getFieldPattern();
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField " << i << " pattern:   ";
+      for (int j=0; j<dofmgr.getLocalFieldSize(i); ++j) {
+        *outStream << fieldPattern[i][j] << " ";
+      }
+      *outStream << std::endl;
+    }
+
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField  " << i << std::endl;
+      // *outStream << dofmgr.getFieldDofs(i);
+    }
+
+    bool correct = true;
+    static const int checkDofs[] = {20, 23, 35, 32, 55, 63, 69, 61, 81};
+    for (int i=0; i<dofmgr.getLocalFieldSize(2); ++i) {
+      correct = correct && ( (dofmgr.getFieldDofs(2))(5,i) == checkDofs[i] );
+    }
+    if (!correct) {
+      errorFlag = -1;
+    }
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_02K.cpp b/example/PDE-OPT/TEST/test_02K.cpp
new file mode 100644
index 000000000..a387506d4
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_02K.cpp
@@ -0,0 +1,164 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_02.cpp
+    \brief Unit test for the mesh manager and the degree-of-freedom manager.
+           Mesh type: BACKWARD-FACING STEP with QUAD CELLS and HGRAD SPACE.
+*/
+
+#include "ROL_Algorithm.hpp"
+#include "ROL_BoundConstraint_SimOpt.hpp"
+#include "ROL_Vector_SimOpt.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
+#include "Intrepid2_HGRAD_QUAD_C2_FEM.hpp"
+
+#include <iostream>
+#include <algorithm>
+
+#include "../TOOLS/meshmanagerK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+
+using scalar_view = typename MeshManager<RealT, DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT, DeviceT>::int_view;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_02.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshManager_BackwardFacingStepChannel<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+    int_view cellToEdgeMap = meshmgr.getCellToEdgeMap();
+    ROL::Ptr<std::vector<std::vector<std::vector<int> > > > sideSetsPtr = meshmgr.getSideSets();
+
+    std::vector<std::vector<std::vector<int> > >  &sideSets = *sideSetsPtr;
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl; // << nodes;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl; // << cellToNodeMap;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl; // << cellToEdgeMap;
+    // Print mesh info to file.
+    std::ofstream meshfile;
+    meshfile.open("cells.txt");
+    for (int i=0; i<cellToNodeMap.extent_int(0); ++i) {
+      meshfile << cellToNodeMap(i,0) << "  "
+               << cellToNodeMap(i,1) << "  "
+               << cellToNodeMap(i,2) << "  "
+               << cellToNodeMap(i,3) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("edges.txt");
+    for (int i=0; i<cellToEdgeMap.extent_int(0); ++i) {
+      meshfile << cellToEdgeMap(i,0) << "  "
+               << cellToEdgeMap(i,1) << "  "
+               << cellToEdgeMap(i,2) << "  "
+               << cellToEdgeMap(i,3) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("nodes.txt");
+    for (int i=0; i<nodes.extent_int(0); ++i) {
+      meshfile << nodes(i,0) << "  "
+               << nodes(i,1) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("sideset.txt");
+    for (int i=0; i<static_cast<int>(sideSets.size()); ++i) {
+      for (int j=0; j<static_cast<int>(sideSets[i].size()); ++j) {
+        if (sideSets[i][j].size() > 0) {
+          for (int k=0; k<static_cast<int>(sideSets[i][j].size()); ++k) {
+            meshfile << sideSets[i][j][k] << std::endl;
+          }
+        }
+        meshfile << std::endl << std::endl;
+      }
+    }
+    meshfile.close();
+
+    basis_ptr basisPtrQ1 = ROL::makePtr<Intrepid2::Basis_HGRAD_QUAD_C1_FEM<DeviceT, RealT, RealT>>();;
+    basis_ptr basisPtrQ2 = ROL::makePtr<Intrepid2::Basis_HGRAD_QUAD_C2_FEM<DeviceT, RealT, RealT>>();;
+
+    std::vector<basis_ptr> basisPtrs(3, ROL::nullPtr);
+    basisPtrs[0] = basisPtrQ2;
+    basisPtrs[1] = basisPtrQ1;
+    basisPtrs[2] = basisPtrQ2;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT>> meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    DofManager<RealT, DeviceT> dofmgr(meshmgrPtr, basisPtrs);
+
+    *outStream << "Number of node dofs = " << dofmgr.getNumNodeDofs() << std::endl; // << *(dofmgr.getNodeDofs());
+    *outStream << "Number of edge dofs = " << dofmgr.getNumEdgeDofs() << std::endl; // << *(dofmgr.getEdgeDofs());
+    *outStream << "Number of face dofs = " << dofmgr.getNumFaceDofs() << std::endl; // << *(dofmgr.getFaceDofs());
+    *outStream << "Number of void dofs = " << dofmgr.getNumVoidDofs() << std::endl; // << *(dofmgr.getVoidDofs());
+    *outStream << "Total number of dofs = " << dofmgr.getNumDofs() << std::endl;    // << *(dofmgr.getCellDofs());
+
+    std::vector<std::vector<int>> fieldPattern = dofmgr.getFieldPattern();
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField " << i << " pattern:   ";
+      for (int j=0; j<dofmgr.getLocalFieldSize(i); ++j) {
+        *outStream << fieldPattern[i][j] << " ";
+      }
+      *outStream << std::endl;
+    }
+
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField  " << i << std::endl;
+      // *outStream << *(dofmgr.getFieldDofs(i));
+    }
+
+    bool correct = true;
+    static const int checkDofs[] = {149, 152, 176, 173, 350, 366, 380, 364, 468};
+    for (int i=0; i<dofmgr.getLocalFieldSize(2); ++i) {
+      correct = correct && ( (dofmgr.getFieldDofs(2))(43,i) == checkDofs[i] );
+    }
+    if (!correct) {
+      errorFlag = -1;
+    }
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_03K.cpp b/example/PDE-OPT/TEST/test_03K.cpp
new file mode 100644
index 000000000..ff259287a
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_03K.cpp
@@ -0,0 +1,172 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_03.cpp
+    \brief Unit test for the mesh manager and the degree-of-freedom manager.
+           Mesh type: BRICK with HEX CELLS and HCURL SPACE.
+*/
+
+#include "ROL_Algorithm.hpp"
+#include "ROL_BoundConstraint_SimOpt.hpp"
+#include "ROL_Vector_SimOpt.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HCURL_HEX_I1_FEM.hpp"
+
+#include <iostream>
+#include <algorithm>
+
+#include "../TOOLS/meshmanagerK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+using scalar_view = typename MeshManager<RealT,DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT,DeviceT>::int_view;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_03.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshManager_Brick<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+    int_view cellToEdgeMap = meshmgr.getCellToEdgeMap();
+    ROL::Ptr<std::vector<std::vector<std::vector<int> > > > sideSetsPtr = meshmgr.getSideSets();
+
+    std::vector<std::vector<std::vector<int> > >  &sideSets = *sideSetsPtr;
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl;  // << nodes;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl;  // << cellToNodeMap;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl;  // << cellToEdgeMap;
+    // Print mesh info to file.
+    std::ofstream meshfile;
+    meshfile.open("cells.txt");
+    for (int i=0; i<cellToNodeMap.extent_int(0); ++i) {
+      meshfile << cellToNodeMap(i,0) << "  "
+               << cellToNodeMap(i,1) << "  "
+               << cellToNodeMap(i,2) << "  "
+               << cellToNodeMap(i,3) << "  "
+               << cellToNodeMap(i,4) << "  "
+               << cellToNodeMap(i,5) << "  "
+               << cellToNodeMap(i,6) << "  "
+               << cellToNodeMap(i,7) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("edges.txt");
+    for (int i=0; i<cellToEdgeMap.extent_int(0); ++i) {
+      meshfile << cellToEdgeMap(i,0)  << "  "
+               << cellToEdgeMap(i,1)  << "  "
+               << cellToEdgeMap(i,2)  << "  "
+               << cellToEdgeMap(i,3)  << "  "
+               << cellToEdgeMap(i,4)  << "  "
+               << cellToEdgeMap(i,5)  << "  "
+               << cellToEdgeMap(i,6)  << "  "
+               << cellToEdgeMap(i,7)  << "  "
+               << cellToEdgeMap(i,8)  << "  "
+               << cellToEdgeMap(i,9)  << "  "
+               << cellToEdgeMap(i,10) << "  "
+               << cellToEdgeMap(i,11) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("nodes.txt");
+    for (int i=0; i<nodes.extent_int(0); ++i) {
+      meshfile << nodes(i,0) << "  "
+               << nodes(i,1) << "  "
+               << nodes(i,2) << std::endl;
+    }
+    meshfile.close();
+    meshfile.open("sideset.txt");
+    for (int i=0; i<static_cast<int>(sideSets.size()); ++i) {
+      for (int j=0; j<static_cast<int>(sideSets[i].size()); ++j) {
+        if (sideSets[i][j].size() > 0) {
+          for (int k=0; k<static_cast<int>(sideSets[i][j].size()); ++k) {
+            meshfile << sideSets[i][j][k] << std::endl;
+          }
+        }
+        meshfile << std::endl << std::endl;
+      }
+    }
+    meshfile.close();
+
+    basis_ptr basisPtrNedelec1 = ROL::makePtr<Intrepid2::Basis_HCURL_HEX_I1_FEM<DeviceT, RealT, RealT>>();
+
+    std::vector<basis_ptr> basisPtrs(1, ROL::nullPtr);
+    basisPtrs[0] = basisPtrNedelec1;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT>> meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    DofManager<RealT, DeviceT> dofmgr(meshmgrPtr, basisPtrs);
+
+    *outStream << "Number of node dofs = " << dofmgr.getNumNodeDofs() << std::endl;  // << *(dofmgr.getNodeDofs());
+    *outStream << "Number of edge dofs = " << dofmgr.getNumEdgeDofs() << std::endl;  // << *(dofmgr.getEdgeDofs());
+    *outStream << "Number of face dofs = " << dofmgr.getNumFaceDofs() << std::endl;  // << *(dofmgr.getFaceDofs());
+    *outStream << "Number of void dofs = " << dofmgr.getNumVoidDofs() << std::endl;  // << *(dofmgr.getVoidDofs());
+    *outStream << "Total number of dofs = " << dofmgr.getNumDofs() << std::endl;     // << *(dofmgr.getCellDofs());
+
+    std::vector<std::vector<int> > fieldPattern = dofmgr.getFieldPattern();
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField " << i << " pattern:   ";
+      for (int j=0; j<dofmgr.getLocalFieldSize(i); ++j) {
+        *outStream << fieldPattern[i][j] << " ";
+      }
+      *outStream << std::endl;
+    }
+
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField  " << i << std::endl;
+      // *outStream << dofmgr.getFieldDofs(i);
+    }
+
+    bool correct = true;
+    static const int checkDofs[] = {38, 42, 45, 41, 67, 71, 74, 70, 52, 53, 57, 56};
+    for (int i=0; i<dofmgr.getLocalFieldSize(0); ++i) {
+      correct = correct && ( (dofmgr.getFieldDofs(0))(11,i) == checkDofs[i] );
+    }
+    if (!correct) {
+      errorFlag = -1;
+    }
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_04K.cpp b/example/PDE-OPT/TEST/test_04K.cpp
new file mode 100644
index 000000000..aa02427a0
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_04K.cpp
@@ -0,0 +1,145 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_04.cpp
+    \brief Unit test for the FE class on quads.
+*/
+
+#include "ROL_Types.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HGRAD_QUAD_C1_FEM.hpp"
+#include "Intrepid2_HGRAD_QUAD_C2_FEM.hpp"
+#include "Intrepid2_DefaultCubatureFactory.hpp"
+
+#include "../TOOLS/meshmanagerK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+#include "../TOOLS/feK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+using scalar_view = typename MeshManager<RealT,DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT,DeviceT>::int_view;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_04.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshManager_Rectangle<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT>> meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    // Basis.
+    basis_ptr basis = ROL::makePtr<Intrepid2::Basis_HGRAD_QUAD_C2_FEM<DeviceT, RealT, RealT>>();
+    // Cubature.
+    ROL::Ptr<Intrepid2::Cubature<DeviceT, RealT, RealT>> cubature;
+    shards::CellTopology cellType = basis->getBaseCellTopology();                        // get the cell type from any basis
+    Intrepid2::DefaultCubatureFactory cubFactory;                                         // create cubature factory
+    int cubDegree = 4;                                                                   // set cubature degree, e.g., 2
+    cubature = cubFactory.create<DeviceT, RealT, RealT>(cellType, cubDegree);            // create default cubature
+    // Cell nodes.
+    scalar_view cellNodes("cellNodes", meshmgr.getNumCells(), cellType.getNodeCount(), cellType.getDimension());
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      for (unsigned j=0; j<cellType.getNodeCount(); ++j) {
+        for (unsigned k=0; k<cellType.getDimension(); ++k) {
+          cellNodes(i,j,k) = nodes(cellToNodeMap(i,j), k);
+        }
+      }
+    }
+    // FE object.
+    FE<RealT, DeviceT> fe(cellNodes, basis, cubature);
+
+    // Check integration.
+    scalar_view feVals1("feVals1", meshmgr.getNumCells(), cubature->getNumPoints());
+    scalar_view feVals2("feVals2", meshmgr.getNumCells(), cubature->getNumPoints());
+    scalar_view feGrads1("feGrads1", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feGrads2("feGrads2", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feCoeffs1("feCoeffs1", meshmgr.getNumCells(), basis->getCardinality());
+    scalar_view feCoeffs2("feCoeffs2", meshmgr.getNumCells(), basis->getCardinality());
+    scalar_view feIntegral("feIntegral", meshmgr.getNumCells());
+    // values
+    Kokkos::deep_copy(feCoeffs1, static_cast<RealT>(1));
+    Kokkos::deep_copy(feCoeffs2, static_cast<RealT>(2));
+    fe.evaluateValue(feVals1, feCoeffs1);
+    fe.evaluateValue(feVals2, feCoeffs2);
+    fe.computeIntegral(feIntegral, feVals1, feVals2);
+    RealT valval(0);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      valval += feIntegral(i);
+    }
+    if (std::abs(valval - static_cast<RealT>(1000)) > std::sqrt(ROL::ROL_EPSILON<RealT>())) {
+      errorFlag = -1;
+    }
+    // gradients
+    scalar_view dofCoords("dofCoords", meshmgr.getNumCells(), basis->getCardinality(), cubature->getDimension());
+    fe.computeDofCoords(dofCoords, cellNodes);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      for (int j=0; j<basis->getCardinality(); ++j) {
+        RealT x = dofCoords(i,j,0);
+        RealT y = dofCoords(i,j,1);
+        feCoeffs1(i,j) = x + y;
+        feCoeffs2(i,j) = static_cast<RealT>(2)*x + static_cast<RealT>(3)*y;
+      }
+    }
+    fe.evaluateGradient(feGrads1, feCoeffs1);
+    fe.evaluateGradient(feGrads2, feCoeffs2);
+    fe.computeIntegral(feIntegral, feGrads1, feGrads2);
+    RealT gradgrad(0);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      gradgrad += feIntegral(i);
+    }
+    if (std::abs(gradgrad - static_cast<RealT>(2500)) > std::sqrt(ROL::ROL_EPSILON<RealT>())) {
+      errorFlag = -1;
+    }
+    *outStream << std::setprecision(14) << "\nvalvalIntegral=" << valval << "  " << "gradgradIntegral=" << gradgrad << std::endl;
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_05K.cpp b/example/PDE-OPT/TEST/test_05K.cpp
new file mode 100644
index 000000000..6aee95b49
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_05K.cpp
@@ -0,0 +1,139 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_04.cpp
+    \brief Unit test for the FE_CURL class.
+*/
+
+#include "ROL_Types.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HCURL_HEX_I1_FEM.hpp"
+#include "Intrepid2_DefaultCubatureFactory.hpp"
+
+#include "../TOOLS/meshmanagerK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+#include "../TOOLS/feK_curl.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+using scalar_view = typename MeshManager<RealT,DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT,DeviceT>::int_view;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_05.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshManager_Brick<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT>> meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    // Basis.
+    basis_ptr basis = ROL::makePtr<Intrepid2::Basis_HCURL_HEX_I1_FEM<DeviceT, RealT, RealT>>();
+    // Cubature.
+    ROL::Ptr<Intrepid2::Cubature<DeviceT, RealT, RealT>> cubature;
+    shards::CellTopology cellType = basis->getBaseCellTopology();                        // get the cell type from any basis
+    Intrepid2::DefaultCubatureFactory cubFactory;                                        // create cubature factory
+    int cubDegree = 4;                                                                   // set cubature degree, e.g., 2
+    cubature = cubFactory.create<DeviceT, RealT, RealT>(cellType, cubDegree);            // create default cubature
+    // Cell nodes.
+    scalar_view cellNodes("cellNodes", meshmgr.getNumCells(), cellType.getNodeCount(), cellType.getDimension());
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      for (unsigned j=0; j<cellType.getNodeCount(); ++j) {
+        for (unsigned k=0; k<cellType.getDimension(); ++k) {
+          cellNodes(i,j,k) = nodes(cellToNodeMap(i,j), k);
+        }
+      }
+    }
+    // FE object.
+    FE_CURL<RealT, DeviceT> fe(cellNodes, basis, cubature);
+
+    // Check integration.
+    scalar_view feVals1("feVals1", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feVals2("feVals2", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feCurls1("feCurls1", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feCurls2("feCurls2", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feCoeffs1("feCoeffs1", meshmgr.getNumCells(), basis->getCardinality());
+    scalar_view feCoeffs2("feCoeffs2", meshmgr.getNumCells(), basis->getCardinality());
+    scalar_view feIntegral("feIntegral", meshmgr.getNumCells());
+    // values
+    Kokkos::deep_copy(feCoeffs1,static_cast<RealT>(0));
+    Kokkos::deep_copy(feCoeffs2,static_cast<RealT>(0));
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      feCoeffs1(i, 0) = static_cast<RealT>(1);
+      feCoeffs2(i, 0) = static_cast<RealT>(1);
+    }
+    fe.evaluateValue(feVals1, feCoeffs1);
+    fe.evaluateValue(feVals2, feCoeffs2);
+    fe.computeIntegral(feIntegral, feVals1, feVals2);
+    RealT valval(0);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      valval += feIntegral(i);
+    }
+    if (std::abs(valval - static_cast<RealT>(4)/static_cast<RealT>(9)) > std::sqrt(ROL::ROL_EPSILON<RealT>())) {
+      errorFlag = -1;
+    }
+    // curls
+    fe.evaluateCurl(feCurls1, feCoeffs1);
+    fe.evaluateCurl(feCurls2, feCoeffs2);
+    fe.computeIntegral(feIntegral, feCurls1, feCurls2);
+    RealT curlcurl(0);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      curlcurl += feIntegral(i);
+    }
+    if (std::abs(curlcurl - static_cast<RealT>(8)/static_cast<RealT>(3)) > std::sqrt(ROL::ROL_EPSILON<RealT>())) {
+      errorFlag = -1;
+    }
+    *outStream << std::setprecision(14) << "\nvalvalIntegral=" << valval << "  "
+                                        << "curlcurlIntegral=" << curlcurl << std::endl;
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_07K.cpp b/example/PDE-OPT/TEST/test_07K.cpp
new file mode 100644
index 000000000..c137c7cec
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_07K.cpp
@@ -0,0 +1,137 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_07.cpp
+    \brief Unit test for the MeshReader mesh manager.
+*/
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HGRAD_HEX_C1_FEM.hpp"
+#include "Intrepid2_HGRAD_HEX_C2_FEM.hpp"
+
+#include "../TOOLS/meshreaderK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+using scalar_view = typename MeshManager<RealT,DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT,DeviceT>::int_view;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_07.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshReader<RealT, DeviceT> meshmgr(*parlist);
+
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+    int_view cellToEdgeMap = meshmgr.getCellToEdgeMap();
+    int_view cellToFaceMap = meshmgr.getCellToFaceMap();
+    ROL::Ptr<std::vector<std::vector<std::vector<int> > > > sideSetsPtr = meshmgr.getSideSets();
+
+    std::vector<std::vector<std::vector<int> > >  &sideSets = *sideSetsPtr;
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl;  // << nodes;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl;  // << cellToNodeMap;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl;  // << cellToEdgeMap;
+    *outStream << "Number of faces = " << meshmgr.getNumFaces() << std::endl;  // << cellToFaceMap;
+    // Print sideset info to file.
+    std::ofstream meshfile;
+    meshfile.open("sideset.txt");
+    for (int i=0; i<static_cast<int>(sideSets.size()); ++i) {
+      for (int j=0; j<static_cast<int>(sideSets[i].size()); ++j) {
+        if (sideSets[i][j].size() > 0) {
+          for (int k=0; k<static_cast<int>(sideSets[i][j].size()); ++k) {
+            meshfile << i << " " << j << " " << k << " " << sideSets[i][j][k] << std::endl;
+          }
+        }
+        meshfile << std::endl;
+      }
+      meshfile << std::endl << std::endl;
+    }
+    meshfile.close();
+
+    basis_ptr basisPtrQ1 = ROL::makePtr<Intrepid2::Basis_HGRAD_HEX_C1_FEM<DeviceT, RealT, RealT>>();
+    basis_ptr basisPtrQ2 = ROL::makePtr<Intrepid2::Basis_HGRAD_HEX_C2_FEM<DeviceT, RealT, RealT>>();
+
+    std::vector<basis_ptr> basisPtrs(3, ROL::nullPtr);
+    basisPtrs[0] = basisPtrQ2;
+    basisPtrs[1] = basisPtrQ1;
+    basisPtrs[2] = basisPtrQ2;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT>> meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    DofManager<RealT, DeviceT> dofmgr(meshmgrPtr, basisPtrs);
+
+    *outStream << "Number of node dofs = " << dofmgr.getNumNodeDofs() << std::endl;  // << *(dofmgr.getNodeDofs());
+    *outStream << "Number of edge dofs = " << dofmgr.getNumEdgeDofs() << std::endl;  // << *(dofmgr.getEdgeDofs());
+    *outStream << "Number of face dofs = " << dofmgr.getNumFaceDofs() << std::endl;  // << *(dofmgr.getFaceDofs());
+    *outStream << "Number of void dofs = " << dofmgr.getNumVoidDofs() << std::endl;  // << *(dofmgr.getVoidDofs());
+    *outStream << "Total number of dofs = " << dofmgr.getNumDofs() << std::endl;     // << *(dofmgr.getCellDofs());
+
+    std::vector<std::vector<int> > fieldPattern = dofmgr.getFieldPattern();
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField " << i << " pattern:   ";
+      for (int j=0; j<dofmgr.getLocalFieldSize(i); ++j) {
+        *outStream << fieldPattern[i][j] << " ";
+      }
+      *outStream << std::endl;
+    }
+
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField  " << i << std::endl;
+      // *outStream << *(dofmgr.getFieldDofs(i));
+    }
+
+    bool correct = true;
+    static const int checkDofs[] = {11, 8, 38, 41, 23, 20, 44, 47, 53, 89, 91, 93, 61, 95, 97, 99, 71, 69, 101, 103, 109, 127, 129, 131, 133, 135, 141};
+    for (int i=0; i<dofmgr.getLocalFieldSize(2); ++i) {
+      correct = correct && ( (dofmgr.getFieldDofs(2))(2,i) == checkDofs[i] );
+    }
+    if (!correct) {
+      errorFlag = -1;
+    }
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_08K.cpp b/example/PDE-OPT/TEST/test_08K.cpp
new file mode 100644
index 000000000..4ad6e8861
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_08K.cpp
@@ -0,0 +1,147 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_08.cpp
+    \brief Unit test for the FE class on hexes.
+*/
+
+#include "ROL_Types.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HGRAD_HEX_C1_FEM.hpp"
+#include "Intrepid2_HGRAD_HEX_C2_FEM.hpp"
+#include "Intrepid2_DefaultCubatureFactory.hpp"
+
+#include "../TOOLS/meshreaderK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+#include "../TOOLS/feK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+using scalar_view = typename MeshManager<RealT,DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT,DeviceT>::int_view;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_08.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshReader<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view  cellToNodeMap = meshmgr.getCellToNodeMap();
+
+    *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl;
+    *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl;
+    *outStream << "Number of edges = " << meshmgr.getNumEdges() << std::endl;
+    *outStream << "Number of edges = " << meshmgr.getNumFaces() << std::endl;
+
+    ROL::Ptr<MeshManager<RealT, DeviceT>> meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    // Basis.
+    basis_ptr basis = ROL::makePtr<Intrepid2::Basis_HGRAD_HEX_C2_FEM<DeviceT, RealT, RealT >>();
+    // Cubature.
+    ROL::Ptr<Intrepid2::Cubature<DeviceT, RealT, RealT>> cubature;
+    shards::CellTopology cellType = basis->getBaseCellTopology();              // get the cell type from any basis
+    Intrepid2::DefaultCubatureFactory cubFactory;                              // create cubature factory
+    int cubDegree = 6;                                                         // set cubature degree, e.g., 2
+    cubature = cubFactory.create<DeviceT, RealT, RealT>(cellType, cubDegree);  // create default cubature
+    // Cell nodes.
+    scalar_view cellNodes("cellNodes", meshmgr.getNumCells(), cellType.getNodeCount(), cellType.getDimension());
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      for (unsigned j=0; j<cellType.getNodeCount(); ++j) {
+        for (unsigned k=0; k<cellType.getDimension(); ++k) {
+          cellNodes(i,j,k) = nodes(cellToNodeMap(i,j), k);
+        }
+      }
+    }
+    // FE object.
+    FE<RealT, DeviceT> fe(cellNodes, basis, cubature);
+
+    // Check integration.
+    scalar_view feVals1("feVals1", meshmgr.getNumCells(), cubature->getNumPoints());
+    scalar_view feVals2("feVals2", meshmgr.getNumCells(), cubature->getNumPoints());
+    scalar_view feGrads1("feGrads1", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feGrads2("feGrads2", meshmgr.getNumCells(), cubature->getNumPoints(), cubature->getDimension());
+    scalar_view feCoeffs1("feCoeffs1", meshmgr.getNumCells(), basis->getCardinality());
+    scalar_view feCoeffs2("feCoeffs2", meshmgr.getNumCells(), basis->getCardinality());
+    scalar_view feIntegral("feIntegral", meshmgr.getNumCells());
+    // values
+    Kokkos::deep_copy(feCoeffs1,static_cast<RealT>(1));
+    Kokkos::deep_copy(feCoeffs2,static_cast<RealT>(2));
+    fe.evaluateValue(feVals1, feCoeffs1);
+    fe.evaluateValue(feVals2, feCoeffs2);
+    fe.computeIntegral(feIntegral, feVals1, feVals2);
+    RealT valval(0);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      valval += feIntegral(i);
+    }
+    if (std::abs(valval - static_cast<RealT>(6)) > std::sqrt(ROL::ROL_EPSILON<RealT>())) {
+      errorFlag = -1;
+    }
+    // gradients
+    scalar_view dofCoords("dofCoords", meshmgr.getNumCells(), basis->getCardinality(), cubature->getDimension());
+    fe.computeDofCoords(dofCoords, cellNodes);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      for (int j=0; j<basis->getCardinality(); ++j) {
+        RealT x = dofCoords(i,j,0);
+        RealT y = dofCoords(i,j,1);
+        RealT z = dofCoords(i,j,2);
+        feCoeffs1(i,j) = x + y + z;
+        feCoeffs2(i,j) = static_cast<RealT>(2)*x + static_cast<RealT>(3)*y + static_cast<RealT>(4)*z;
+      }
+    }
+    fe.evaluateGradient(feGrads1, feCoeffs1);
+    fe.evaluateGradient(feGrads2, feCoeffs2);
+    fe.computeIntegral(feIntegral, feGrads1, feGrads2);
+    RealT gradgrad(0);
+    for (int i=0; i<meshmgr.getNumCells(); ++i) {
+      gradgrad += feIntegral(i);
+    }
+    if (std::abs(gradgrad - static_cast<RealT>(27)) > std::sqrt(ROL::ROL_EPSILON<RealT>())) {
+      errorFlag = -1;
+    }
+    *outStream << std::setprecision(14) << "\nvalvalIntegral=" << valval << "  " << "gradgradIntegral=" << gradgrad << std::endl;
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TEST/test_09K.cpp b/example/PDE-OPT/TEST/test_09K.cpp
new file mode 100644
index 000000000..745c2f9f1
--- /dev/null
+++ b/example/PDE-OPT/TEST/test_09K.cpp
@@ -0,0 +1,139 @@
+// @HEADER
+// *****************************************************************************
+//               Rapid Optimization Library (ROL) Package
+//
+// Copyright 2014 NTESS and the ROL contributors.
+// SPDX-License-Identifier: BSD-3-Clause
+// *****************************************************************************
+// @HEADER
+
+/*! \file  test_09.cpp
+    \brief Unit test for the mesh manager and the degree-of-freedom manager.
+           Mesh type: INTERVAL with LINE CELLS and HGRAD SPACE.
+*/
+
+#include "ROL_Algorithm.hpp"
+#include "ROL_BoundConstraint_SimOpt.hpp"
+#include "ROL_Vector_SimOpt.hpp"
+
+#include "ROL_Stream.hpp"
+#include "ROL_GlobalMPISession.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+
+#include "Intrepid2_HGRAD_LINE_Cn_FEM.hpp"
+
+#include <iostream>
+#include <algorithm>
+
+#include "../TOOLS/meshmanagerK.hpp"
+#include "../TOOLS/dofmanagerK.hpp"
+
+using RealT = double;
+using DeviceT = Kokkos::HostSpace;
+using basis_ptr = Intrepid2::BasisPtr<DeviceT, RealT, RealT>;
+using scalar_view = typename MeshManager<RealT,DeviceT>::scalar_view;
+using int_view = typename MeshManager<RealT,DeviceT>::int_view;
+
+int main(int argc, char *argv[]) {
+
+  ROL::GlobalMPISession mpiSession(&argc, &argv);
+  Kokkos::ScopeGuard kokkosScope(argc, argv);
+  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
+  int iprint     = argc - 1;
+  ROL::Ptr<std::ostream> outStream;
+  ROL::nullstream bhs; // outputs nothing
+  if (iprint > 0)
+    outStream = ROL::makePtrFromRef(std::cout);
+  else
+    outStream = ROL::makePtrFromRef(bhs);
+
+  int errorFlag  = 0;
+
+  // *** Example body.
+  try {
+
+    /*** Read in XML input ***/
+    std::string filename = "input_09.xml";
+    Teuchos::RCP<Teuchos::ParameterList> parlist
+      = Teuchos::rcp( new Teuchos::ParameterList() );
+    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
+
+    /*** Initialize mesh / degree-of-freedom manager. ***/
+    MeshManager_Interval<RealT, DeviceT> meshmgr(*parlist);
+    scalar_view nodes = meshmgr.getNodes();
+    int_view cellToNodeMap = meshmgr.getCellToNodeMap();
+    ROL::Ptr<std::vector<std::vector<std::vector<int> > > > sideSetsPtr = meshmgr.getSideSets();
+
+    std::vector<std::vector<std::vector<int> > >  &sideSets = *sideSetsPtr;
+    // *outStream << "Number of nodes = " << meshmgr.getNumNodes() << std::endl << nodes;
+    // *outStream << "Number of cells = " << meshmgr.getNumCells() << std::endl << cellToNodeMap;
+    // Print mesh info to file.
+    std::ofstream meshfile;
+    meshfile.open("sideset.txt");
+    for (int i=0; i<static_cast<int>(sideSets.size()); ++i) {
+      for (int j=0; j<static_cast<int>(sideSets[i].size()); ++j) {
+        if (sideSets[i][j].size() > 0) {
+          for (int k=0; k<static_cast<int>(sideSets[i][j].size()); ++k) {
+            meshfile << sideSets[i][j][k] << std::endl;
+          }
+        }
+        meshfile << std::endl << std::endl;
+      }
+    }
+    meshfile.close();
+
+    basis_ptr basisPtrL1 = ROL::makePtr<Intrepid2::Basis_HGRAD_LINE_Cn_FEM<DeviceT, RealT, RealT >>(1, Intrepid2::POINTTYPE_EQUISPACED);
+
+    basis_ptr basisPtrL2 = ROL::makePtr<Intrepid2::Basis_HGRAD_LINE_Cn_FEM<DeviceT, RealT, RealT >>(2, Intrepid2::POINTTYPE_EQUISPACED);
+
+    std::vector<basis_ptr> basisPtrs(3, ROL::nullPtr);
+    basisPtrs[0] = basisPtrL2;
+    basisPtrs[1] = basisPtrL1;
+    basisPtrs[2] = basisPtrL2;
+
+    ROL::Ptr<MeshManager<RealT,DeviceT> > meshmgrPtr = ROL::makePtrFromRef(meshmgr);
+
+    DofManager<RealT,DeviceT> dofmgr(meshmgrPtr, basisPtrs);
+
+    *outStream << "Number of node dofs = " << dofmgr.getNumNodeDofs() << std::endl; // << *(dofmgr.getNodeDofs());
+    *outStream << "Number of edge dofs = " << dofmgr.getNumEdgeDofs() << std::endl; // << *(dofmgr.getEdgeDofs());
+    *outStream << "Number of face dofs = " << dofmgr.getNumFaceDofs() << std::endl; // << *(dofmgr.getFaceDofs());
+    *outStream << "Number of void dofs = " << dofmgr.getNumVoidDofs() << std::endl; // << *(dofmgr.getVoidDofs());
+    *outStream << "Total number of dofs = " << dofmgr.getNumDofs() << std::endl;    // << *(dofmgr.getCellDofs());
+
+    std::vector<std::vector<int> > fieldPattern = dofmgr.getFieldPattern();
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField " << i << " pattern:   ";
+      for (int j=0; j<dofmgr.getLocalFieldSize(i); ++j) {
+        *outStream << fieldPattern[i][j] << " ";
+      }
+      *outStream << std::endl;
+    }
+
+    for (int i=0; i<dofmgr.getNumFields(); ++i) {
+      *outStream << "\nField  " << i << std::endl;
+      // *outStream << *(dofmgr.getFieldDofs(i));
+    }
+
+    /*bool correct = true;
+    static const int checkDofs[] = {20, 23, 35, 32, 55, 63, 69, 61, 81};
+    for (int i=0; i<dofmgr.getLocalFieldSize(2); ++i) {
+      correct = correct && ( (*(dofmgr.getFieldDofs(2)))(5,i) == checkDofs[i] );
+    }
+    if (!correct) {
+      errorFlag = -1;
+    }*/
+
+  }
+  catch (std::logic_error& err) {
+    *outStream << err.what() << "\n";
+    errorFlag = -1000;
+  }; // end try
+
+  if (errorFlag != 0)
+    std::cout << "End Result: TEST FAILED\n";
+  else
+    std::cout << "End Result: TEST PASSED\n";
+
+  return 0;
+}
diff --git a/example/PDE-OPT/TOOLS/feK_curl.hpp b/example/PDE-OPT/TOOLS/feK_curl.hpp
index 7c41f28a3..945c1f436 100644
--- a/example/PDE-OPT/TOOLS/feK_curl.hpp
+++ b/example/PDE-OPT/TOOLS/feK_curl.hpp
@@ -81,7 +81,7 @@ class FE_CURL {
     cellTopo_ = ROL::makePtr<shards::CellTopology>(basis_->getBaseCellTopology());
 
     // Compute dimensions of multidimensional array members.
-    c_  = cellNodes_->extent_int(0);
+    c_  = cellNodes_.extent_int(0);
     f_  = basis_->getCardinality();
     p_  = cubature_->getNumPoints();
     d_  = cellTopo_->getDimension();
@@ -126,7 +126,7 @@ class FE_CURL {
     fst::computeCellMeasure(cellWeightedMeasure_, cellJacDet_, cubWeights_);
 
     // Transform reference values into physical space.
-    fst::HCURLtransformVALUE(valPhysical_, valReference_);
+    fst::HCURLtransformVALUE(valPhysical_, cellJacInv_, valReference_);
 
     // Multiply with weighted measure to get weighted values in physical space.
     fst::multiplyMeasure(valPhysicalWeighted_, cellWeightedMeasure_, valPhysical_);
@@ -135,7 +135,7 @@ class FE_CURL {
     fst::HCURLtransformCURL(curlPhysical_, cellJac_, cellJacDet_, curlReference_);
 
     // Multiply with weighted measure to get weighted curls in physical space.
-    fst::multiplyMeasure::multiplyMeasure(curlPhysicalWeighted_, cellWeightedMeasure_, curlPhysical_);
+    fst::multiplyMeasure(curlPhysicalWeighted_, cellWeightedMeasure_, curlPhysical_);
 
     // Compute curl-curl matrices.
     fst::integrate(curlcurlMats_, curlPhysical_, curlPhysicalWeighted_);
@@ -307,7 +307,7 @@ class FE_CURL {
       f2Weighted = scalar_view("f2Weighted", f2.extent(0), f2.extent(1));
     }
     fst::scalarMultiplyDataData(f2Weighted, cellWeightedMeasure_, f2);  // multiply with weighted measure
-                                                              
+
     fst::integrate(integral, f1, f2Weighted, sumInto);                  // compute integral of f1*f2
   }
 
diff --git a/example/PDE-OPT/TOOLS/meshreaderK.hpp b/example/PDE-OPT/TOOLS/meshreaderK.hpp
index f97cd4d63..56c35fa01 100644
--- a/example/PDE-OPT/TOOLS/meshreaderK.hpp
+++ b/example/PDE-OPT/TOOLS/meshreaderK.hpp
@@ -499,7 +499,7 @@ class MeshReader : public MeshManager<Real, DeviceType> {
 
   ROL::Ptr<std::vector<std::vector<std::vector<int> > > > getSideSets (
       const bool verbose = false,
-      std::ostream & outStream = std::cout) const {
+      std::ostream & outStream = std::cout) const override {
     if (verbose) {
       for (int i=0; i<numSideSets_; ++i) {
         outStream << "\nSideset " << i << std::endl;
diff --git a/example/PDE-OPT/helmholtz/CMakeLists.txt b/example/PDE-OPT/helmholtz/CMakeLists.txt
index c84f46fb5..d0d94030c 100644
--- a/example/PDE-OPT/helmholtz/CMakeLists.txt
+++ b/example/PDE-OPT/helmholtz/CMakeLists.txt
@@ -65,7 +65,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid2 AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -76,7 +75,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid2 AND
     ARGS PrintItAll
     NUM_MPI_PROCS 1
     NUM_TOTAL_CORES_USED 1
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
diff --git a/example/PDE-OPT/navier-stokes/CMakeLists.txt b/example/PDE-OPT/navier-stokes/CMakeLists.txt
index c81bef6a4..cf696baf7 100644
--- a/example/PDE-OPT/navier-stokes/CMakeLists.txt
+++ b/example/PDE-OPT/navier-stokes/CMakeLists.txt
@@ -114,10 +114,11 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid2 AND
   TRIBITS_ADD_TEST(
     example_01_Kokkos
     ARGS PrintItAll
-    NUM_MPI_PROCS 4
-    NUM_TOTAL_CORES_USED 4
+    NUM_MPI_PROCS 2
+    NUM_TOTAL_CORES_USED 2
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
+    RUN_SERIAL
   )
 
   #TRIBITS_ADD_TEST(
diff --git a/example/PDE-OPT/navier-stokes/input.xml b/example/PDE-OPT/navier-stokes/input.xml
index 7e54d3410..300af77e3 100644
--- a/example/PDE-OPT/navier-stokes/input.xml
+++ b/example/PDE-OPT/navier-stokes/input.xml
@@ -28,7 +28,7 @@
     <Parameter name="Step width"         type="double"  value="1.0" />
     <Parameter name="Step height"        type="double"  value="0.5" />
     <Parameter name="Observation width"  type="double"  value="3.0" />
-    <Parameter name="Refinement level"   type="int"     value="2"   />
+    <Parameter name="Refinement level"   type="int"     value="1"   />
     <Parameter name="Partition type"     type="int"     value="1"    />
   </ParameterList>
 
@@ -152,7 +152,7 @@
   <!-- ===========  STEP SUBLIST  =========== -->
   <ParameterList name="Step">
     <Parameter name="Type" type="string" value="Fletcher"/>
-    
+
     <!-- ===========  LINE SEARCH  =========== -->
     <ParameterList name="Line Search">
       <Parameter name="Function Evaluation Limit"      type="int"    value="20"    />
@@ -164,7 +164,7 @@
 
       <!-- ===========  DESCENT ALGORITHM SPECIFICATION  =========== -->
       <ParameterList name="Descent Method">
-        <Parameter name="Type"              type="string" value="Newton-Krylov"    /> 
+        <Parameter name="Type"              type="string" value="Newton-Krylov"    />
         <Parameter name="Nonlinear CG Type" type="string" value="Hestenes-Stiefel" />
       </ParameterList>
 
diff --git a/example/PDE-OPT/stefan-boltzmann/CMakeLists.txt b/example/PDE-OPT/stefan-boltzmann/CMakeLists.txt
index a11df5df9..a4df0d42a 100644
--- a/example/PDE-OPT/stefan-boltzmann/CMakeLists.txt
+++ b/example/PDE-OPT/stefan-boltzmann/CMakeLists.txt
@@ -17,7 +17,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -28,7 +27,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -39,7 +37,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -84,7 +81,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid2 AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -95,7 +91,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid2 AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -106,7 +101,6 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid2 AND
     ARGS PrintItAll
     NUM_MPI_PROCS 4
     NUM_TOTAL_CORES_USED 4
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
diff --git a/example/topology-optimization/CMakeLists.txt b/example/topology-optimization/CMakeLists.txt
index 5ad1e4782..b8fb999bc 100644
--- a/example/topology-optimization/CMakeLists.txt
+++ b/example/topology-optimization/CMakeLists.txt
@@ -1,7 +1,7 @@
 
 ROL_ADD_EXECUTABLE(
   example_02
-  SOURCES example_02.cpp 
+  SOURCES example_02.cpp
   ADD_DIR_TO_NAME
 )
 
@@ -24,7 +24,6 @@ if( ${lower_build_type} STREQUAL "release" )
     SOURCES example_01.cpp
     ARGS PrintItAll
     NUM_MPI_PROCS 1
-    CATEGORIES CONTINUOUS
     PASS_REGULAR_EXPRESSION "TEST PASSED"
     ADD_DIR_TO_NAME
   )
@@ -33,7 +32,7 @@ else()
 
   ROL_ADD_EXECUTABLE(
     example_01
-    SOURCES example_01.cpp 
+    SOURCES example_01.cpp
     ADD_DIR_TO_NAME
   )
 
diff --git a/src/step/krylov/ROL_GMRES.hpp b/src/step/krylov/ROL_GMRES.hpp
index 74ba97ceb..7e924359a 100644
--- a/src/step/krylov/ROL_GMRES.hpp
+++ b/src/step/krylov/ROL_GMRES.hpp
@@ -22,6 +22,48 @@
 
 namespace ROL {
 
+template<class Real>
+class VectorArray {
+private:
+  std::vector<Ptr<Vector<Real>>> data_;
+  unsigned size_;
+
+public:
+  VectorArray() : size_(0u) {}
+
+  void allocate(const Vector<Real> &x, unsigned numVec) {
+    data_.clear();
+    for (unsigned i = 0u; i < numVec; ++i)
+      data_.push_back(x.clone());
+    size_ = numVec;
+  }
+
+  const Ptr<Vector<Real>> get(const Vector<Real> &x, unsigned ind) {
+    if (ind < size_) {
+      return data_[ind];
+    }
+    else if (ind == size_) {
+      data_.push_back(x.clone());
+      size_++;
+      return data_[ind];
+    }
+    else {
+      throw Exception::NotImplemented(">>> GMRES : Index out of range!");
+      return nullPtr;
+    }
+  }
+
+  const Ptr<Vector<Real>> get(unsigned ind) {
+    if (size_ > 0u) {
+      return get(*data_[0],ind);
+    }
+    else {
+      throw Exception::NotImplemented(">>> GMRES : No vectors allocated!");
+      return nullPtr;
+    }
+  }
+};
+
 template<class Real>
 class GMRES : public Krylov<Real> {
 
@@ -30,30 +72,31 @@ class GMRES : public Krylov<Real> {
 
 private:
 
-  ROL::Ptr<Vector<Real> > r_;
-  ROL::Ptr<Vector<Real> > z_;
-  ROL::Ptr<Vector<Real> > w_;
+  Ptr<Vector<Real>> r_;
+  Ptr<Vector<Real>> z_;
+  Ptr<Vector<Real>> w_;
 
-  ROL::Ptr<SDMatrix> H_;      // quasi-Hessenberg matrix
-  ROL::Ptr<SDVector> cs_;     // Givens Rotations cosine components
-  ROL::Ptr<SDVector> sn_;     // Givens Rotations sine components
-  ROL::Ptr<SDVector> s_;
-  ROL::Ptr<SDVector> y_;
-  ROL::Ptr<SDVector> cnorm_;
+  Ptr<SDMatrix> H_;      // quasi-Hessenberg matrix
+  Ptr<SDVector> cs_;     // Givens Rotations cosine components
+  Ptr<SDVector> sn_;     // Givens Rotations sine components
+  Ptr<SDVector> s_;
+  Ptr<SDVector> y_;
+  Ptr<SDVector> cnorm_;
 
-  ROL::Ptr<std::vector<Real> > res_;
+  Ptr<std::vector<Real>> res_;
+  Ptr<VectorArray<Real>> V_, Z_;
 
   bool isInitialized_;
   bool useInexact_;
   bool useInitialGuess_;    // If false, inital x will be ignored and zero vec used
   bool printIters_;
-  ROL::Ptr<std::ostream> outStream_;
+  Ptr<std::ostream> outStream_;
 
-  ROL::LAPACK<int,Real> lapack_;
+  LAPACK<int,Real> lapack_;
 
 public:
 
-  GMRES( ROL::ParameterList &parlist ) : Krylov<Real>(parlist), isInitialized_(false), printIters_(false) {
+  GMRES( ParameterList &parlist ) : Krylov<Real>(parlist), isInitialized_(false), printIters_(false) {
 
     using std::vector;
 
@@ -66,14 +109,16 @@ class GMRES : public Krylov<Real> {
     useInitialGuess_ = kList.get("Use Initial Guess",false);
     int maxit = Krylov<Real>::getMaximumIteration();
 
-    H_     = ROL::makePtr<SDMatrix>( maxit+1, maxit );
-    cs_    = ROL::makePtr<SDVector>( maxit );
-    sn_    = ROL::makePtr<SDVector>( maxit );
-    s_     = ROL::makePtr<SDVector>( maxit+1 );
-    y_     = ROL::makePtr<SDVector>( maxit+1 );
-    cnorm_ = ROL::makePtr<SDVector>( maxit );
-    res_   = ROL::makePtr<std::vector<Real>>(maxit+1,zero);
+    H_     = makePtr<SDMatrix>( maxit+1, maxit );
+    cs_    = makePtr<SDVector>( maxit );
+    sn_    = makePtr<SDVector>( maxit );
+    s_     = makePtr<SDVector>( maxit+1 );
+    y_     = makePtr<SDVector>( maxit+1 );
+    cnorm_ = makePtr<SDVector>( maxit );
+    res_   = makePtr<std::vector<Real>>(maxit+1,zero);
 
+    V_ = makePtr<VectorArray<Real>>();
+    Z_ = makePtr<VectorArray<Real>>();
   }
 
   Real run( Vector<Real> &x, LinearOperator<Real> &A, const Vector<Real> &b,
@@ -92,6 +137,9 @@ class GMRES : public Krylov<Real> {
       w_  = b.clone();
       z_  = x.clone();
 
+      V_->allocate(b,std::min(20u,static_cast<unsigned>(maxit))+1u);
+      Z_->allocate(x,std::min(20u,static_cast<unsigned>(maxit)));
+
       isInitialized_ = true;
     }
 
@@ -112,14 +160,10 @@ class GMRES : public Krylov<Real> {
 
     Real temp  = 0;
 
-    std::vector<ROL::Ptr<Vector<Real > > > V;
-    std::vector<ROL::Ptr<Vector<Real > > > Z;
-
     (*res_)[0] = r_->norm();
 
-    if (printIters_) {
+    if (printIters_)
       *outStream_ << "GMRES Iteration " << 0 << ", Residual = " << (*res_)[0] << "\n";
-    }
 
     // This should be a tolerance check
     Real rtol = std::min(absTol,relTol*(*res_)[0]);
@@ -129,42 +173,36 @@ class GMRES : public Krylov<Real> {
       return (*res_)[0];
     }
 
-    V.push_back(b.clone());
-    (V[0])->set(*r_);
-    (V[0])->scale(one/(*res_)[0]);
+    V_->get(0)->set(*r_);
+    V_->get(0)->scale(one/(*res_)[0]);
 
     (*s_)(0) = (*res_)[0];
 
     for( iter=0; iter<maxit; ++iter ) {
 
-      if( useInexact_ ) {
-        itol = rtol/(maxit*(*res_)[iter]);
-      }
-
-      Z.push_back(x.clone());
+      if( useInexact_ ) itol = rtol/(maxit*(*res_)[iter]);
 
       // Apply right preconditioner
-      M.applyInverse(*(Z[iter]),*(V[iter]),itol);
+      M.applyInverse(*(Z_->get(iter)),*(V_->get(iter)),itol);
 
       // Apply operator
-      A.apply(*w_,*(Z[iter]),itol);
+      A.apply(*w_,*(Z_->get(iter)),itol);
 
       // Evaluate coefficients and orthogonalize using Gram-Schmidt
       for( int k=0; k<=iter; ++k ) {
-        (*H_)(k,iter) = w_->dot(*(V[k]));
-        w_->axpy( -(*H_)(k,iter), *(V[k]) );
+        (*H_)(k,iter) = w_->dot(*(V_->get(k)));
+        w_->axpy( -(*H_)(k,iter), *(V_->get(k)) );
       }
 
       (*H_)(iter+1,iter) = w_->norm();
 
-      V.push_back( b.clone() );
-      (V[iter+1])->set(*w_);
-      (V[iter+1])->scale(one/((*H_)(iter+1,iter)));
+      (V_->get(iter+1))->set(*w_);
+      (V_->get(iter+1))->scale(one/((*H_)(iter+1,iter)));
 
       // Apply Givens rotations
       for( int k=0; k<=iter-1; ++k ) {
         temp            =  (*cs_)(k)*(*H_)(k,iter) + (*sn_)(k)*(*H_)(k+1,iter);
-        (*H_)(k+1,iter) = -(*sn_)(k)*(*H_)(k,iter) + (*cs_)(k)*(*H_)(k+1,iter); 
+        (*H_)(k+1,iter) = -(*sn_)(k)*(*H_)(k,iter) + (*cs_)(k)*(*H_)(k+1,iter);
         (*H_)(k,iter)   = temp;
       }
 
@@ -209,7 +247,7 @@ class GMRES : public Krylov<Real> {
       z_->zero();
 
       for( int k=0; k<=iter;++k ) {
-        z_->axpy((*y_)(k),*(Z[k]));
+        z_->axpy((*y_)(k),*(Z_->get(k)));
       }
 
       if( (*res_)[iter+1] <= rtol ) {