diff --git a/README.md b/README.md
index 558cfe0d..166a34db 100644
--- a/README.md
+++ b/README.md
@@ -45,7 +45,7 @@ To build Kokkos with OpenMP and CUDA backend, use:
 cd extern/kokkos
 mkdir build
 cd build
-cmake .. -DCMAKE_INSTALL_PREFIX=../../../installs/kokkos -DKokkos_ENABLE_OPENMP=ON -DKokkos_ENABLE_CUDA=ON -DKokkos_ENABLE_CUDA_LAMBDA=ON -G Ninja
+cmake .. -DCMAKE_INSTALL_PREFIX=../../../installs/kokkos -DCMAKE_CXX_STANDARD=17 -DKokkos_ENABLE_OPENMP=ON -DKokkos_ENABLE_CUDA=ON -DKokkos_ENABLE_CUDA_LAMBDA=ON -DKokkos_ENABLE_CUDA_CONSTEXPR=ON -G Ninja
 ninja install
 ```
 For a complete instruction on installing Kokkos, see [Kokkos
diff --git a/examples/hdiv/hdiv.cpp b/examples/hdiv/hdiv.cpp
index 9fecd23f..5ef853d0 100644
--- a/examples/hdiv/hdiv.cpp
+++ b/examples/hdiv/hdiv.cpp
@@ -207,10 +207,11 @@ void set_geo_spherical(const T alpha, const T R, GeoElemVec& elem_geo) {
   T ulen = 2.0 / nx;
 
   T Xloc[3 * ET::HEX_NVERTS];
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
-    T u = 1.0 * ET::HEX_VERTS_CART[ii][0];
-    T v = 1.0 * ET::HEX_VERTS_CART[ii][1];
-    T w = 1.0 * ET::HEX_VERTS_CART[ii][2];
+    T u = 1.0 * HEX_VERTS_CART[ii][0];
+    T v = 1.0 * HEX_VERTS_CART[ii][1];
+    T w = 1.0 * HEX_VERTS_CART[ii][2];
 
     Xloc[3 * ii] = (2.0 * u - 1.0) * a;
     Xloc[3 * ii + 1] = (2.0 * v - 1.0) * a;
@@ -247,6 +248,7 @@ void set_geo_spherical(const T alpha, const T R, GeoElemVec& elem_geo) {
   }
 
   // We know the 2-direction will be radial...
+  auto HEX_BOUND_VERTS = ET::get_hex_bound_verts();
   for (index_t face = 0; face < ET::HEX_NBOUNDS; face++) {
     for (index_t k = 0; k < nx; k++) {
       for (index_t j = 0; j < nx; j++) {
@@ -274,9 +276,9 @@ void set_geo_spherical(const T alpha, const T R, GeoElemVec& elem_geo) {
             T z0 = 0.0;
 
             for (index_t jj = 0; jj < 4; jj++) {
-              x0 += N[jj] * Xloc[3 * ET::HEX_BOUND_VERTS[face][jj]];
-              y0 += N[jj] * Xloc[3 * ET::HEX_BOUND_VERTS[face][jj] + 1];
-              z0 += N[jj] * Xloc[3 * ET::HEX_BOUND_VERTS[face][jj] + 2];
+              x0 += N[jj] * Xloc[3 * HEX_BOUND_VERTS[face][jj]];
+              y0 += N[jj] * Xloc[3 * HEX_BOUND_VERTS[face][jj] + 1];
+              z0 += N[jj] * Xloc[3 * HEX_BOUND_VERTS[face][jj] + 2];
             }
             T norm = std::sqrt(x0 * x0 + y0 * y0 + z0 * z0);
 
@@ -603,11 +605,12 @@ void find_spherical_error(bool write_sphere = false) {
   }
 
   // For each face add the faces
+  auto HEX_BOUND_VERTS = ET::get_hex_bound_verts();
   for (index_t face = 0; face < ET::HEX_NBOUNDS; face++) {
     for (index_t i = 0; i < 4; i++) {
-      hex[(face + 1) * ET::HEX_NVERTS + i] = hex[ET::HEX_BOUND_VERTS[face][i]];
+      hex[(face + 1) * ET::HEX_NVERTS + i] = hex[HEX_BOUND_VERTS[face][i]];
       hex[(face + 1) * ET::HEX_NVERTS + i + 4] =
-          hex_ext[ET::HEX_BOUND_VERTS[face][i]];
+          hex_ext[HEX_BOUND_VERTS[face][i]];
     }
   }
 
diff --git a/examples/kokkos/scratchpad.cpp b/examples/kokkos/scratchpad.cpp
index 260002b6..400628c0 100644
--- a/examples/kokkos/scratchpad.cpp
+++ b/examples/kokkos/scratchpad.cpp
@@ -6,6 +6,8 @@
 #include "Kokkos_UnorderedMap.hpp"
 #include "a2dobjs.h"
 #include "array.h"
+#include "multiphysics/fequadrature.h"
+#include "multiphysics/fespace.h"
 
 using namespace std;
 using T = double;
@@ -578,22 +580,214 @@ void test_cuda_axpy(int argc, char* argv[]) {
 #endif
 }
 
+void subview() {
+  // constexpr A2D::index_t ndof_per_element = 0;  // won't work
+  constexpr A2D::index_t ndof_per_element = 1;  // ok
+  using ElementDofArray = A2D::MultiArrayNew<A2D::index_t* [ndof_per_element]>;
+  A2D::index_t nelems = 10;
+  ElementDofArray element_dof("element_dof", nelems);
+  A2D::index_t elem = 1;
+  auto elem_dof = Kokkos::subview(element_dof, elem, Kokkos::ALL);
+}
+
+void test_cuda_axpy_with_UVM(int argc, char* argv[]) {
+#ifdef KOKKOS_ENABLE_CUDA
+  using ViewDevice_t =
+      Kokkos::View<T*, Kokkos::LayoutRight, Kokkos::CudaUVMSpace>;
+  if (argc == 1) {
+    std::printf("axpy\nusage: ./scratchpad N, where mat size = 2^N\n");
+    return;
+  }
+
+  // Allocate x and y on host
+  I N = pow(2, atoi(argv[1]));
+  I bytes = N * sizeof(T);
+  double Mbytes = (double)bytes / 1024 / 1024;
+
+  std::printf("Allocating x[%d] on host, size: %.2f MB\n", N, Mbytes);
+  std::printf("Allocating y[%d] on host, size: %.2f MB\n", N, Mbytes);
+
+  // Allocate x and y on device
+  ViewDevice_t x_device("x_device", N);
+  ViewDevice_t y_device("y_device", N);
+
+  for (I i = 0; i < N; i++) {
+    x_device(i) = 2.4;
+    y_device(i) = 0.1;
+  }
+
+  // Perform 100 axpy operations on device and time
+  Kokkos::Timer timer;
+  int repeat = 100;
+  T alpha = 0.01;
+  for (int i = 0; i < repeat; i++) {
+    auto loop_body = KOKKOS_LAMBDA(const I index) {
+      // T alpha = Alpha::get_alpha_array()[0];  // This doesn't work
+      T alpha = Alpha::value_array[0];  // This works
+      auto x_slice = Kokkos::subview(x_device, Kokkos::ALL);
+      y_device(index) += alpha * x_slice(index);
+    };
+    Kokkos::parallel_for("axpy", N, loop_body);
+    Kokkos::fence();
+  }
+
+  double elapse = timer.seconds();
+  printf("averaged time: %.8f ms\n", elapse * 1e3);
+
+  // Compute bandwidth:
+  // x is read once, y is read once and written once
+  double bandwidth = 3 * Mbytes / 1024 * repeat / elapse;
+  printf("averaged bandwidth: %.8f GB/s\n", bandwidth);
+
+  // Print values
+  int len = 10;
+  if (N < len) {
+    len = N;
+  }
+
+  for (I i = 0; i < len; i++) {
+    printf("x[%2d]: %8.2f  y[%2d]: %8.2f\n", i, x_device(i), i, y_device(i));
+  }
+
+  // Check maximum error
+  T max_err = T(0);
+  T val = T(0);
+  for (I i = 0; i < N; i++) {
+    val = fabs(repeat * alpha * 2.4 + 0.1 - y_device(i));
+    if (val > max_err) {
+      max_err = val;
+    }
+  }
+
+  printf("Maximum error: %20.10e\n", max_err);
+#endif
+}
+
+class ParallelVector {
+#ifdef KOKKOS_ENABLE_CUDA
+  using MemSpace = Kokkos::CudaUVMSpace;
+#else
+  using MemSpace = Kokkos::HostSpace;
+#endif
+  using data_t = Kokkos::View<double*, Kokkos::LayoutRight, MemSpace>;
+
+  KOKKOS_FUNCTION void set_one_value(int i, double val) const { data(i) = val; }
+
+ public:
+  ParallelVector(int len) : len(len), data("data", len) {}
+
+  void set_values(double val) {
+    auto loop_body = KOKKOS_CLASS_LAMBDA(int i) { set_one_value(i, val); };
+
+    Kokkos::parallel_for("loop", len, loop_body);
+    Kokkos::fence();
+  }
+
+ private:
+  int len;
+  data_t data;
+};
+
+struct Head {
+  Head(double val) : val(val) {}
+  KOKKOS_FUNCTION double get_val() { return val; };
+
+ private:
+  double val;
+};
+struct Data {
+  Data(Head& head, double val, int id = 0) : val(val), id(id), head(head) {}
+  double val;
+  int id;
+  Head& head;
+};
+
+// This can build, won't work -  ``it is generally not valid to have any pointer
+// or reference members in the functor''
+void test_cuda_functor_pass_by_ref() {
+#ifdef KOKKOS_ENABLE_CUDA
+  Head head(5.6);
+  Data data(head, 4.2);
+  Kokkos::View<double*, Kokkos::LayoutRight, Kokkos::CudaUVMSpace> view("view",
+                                                                        10);
+  auto loop_body = [=] __device__(int i) { view(i) = data.head.get_val(); };
+  Kokkos::parallel_for("loop", 10, loop_body);
+  Kokkos::fence();
+
+  for (int i = 0; i < 10; i++) {
+    std::printf("view[%2d]: %.10f\n", i, view[i]);
+  }
+#endif
+}
+
+void test_KOKKOS_ENABLE_CXX17() {
+#ifdef KOKKOS_ENABLE_CXX17
+  printf("KOKKOS_ENABLE_CXX17 is defined\n");
+#else
+  printf("KOKKOS_ENABLE_CXX17 is not defined\n");
+#endif
+}
+
+void test_fespace() {
+  // Type and literals
+  using T = double;
+  constexpr A2D::index_t dim = 2;
+  constexpr A2D::index_t order = 3;
+
+  using FiniteElementSpace = A2D::FESpace<T, dim, A2D::H1Space<T, dim, dim>>;
+
+  // We loop over 10 elements in parallel
+  int nelems = 10;
+
+#ifdef KOKKOS_ENABLE_CUDA
+  using MemSpace = Kokkos::CudaUVMSpace;
+  printf("using CUDAUVMSpace\n");
+#else
+  using MemSpace = Kokkos::HostSpace;
+#endif
+
+  Kokkos::View<T*, Kokkos::LayoutRight, MemSpace> probe("probe", nelems);
+
+  auto loop_body = KOKKOS_LAMBDA(int i) {
+    FiniteElementSpace cref;
+    A2D::H1Space<T, dim, dim>& s = cref.template get<0>();
+    // A2D::H1Space<T, dim, dim> s;
+    s.get_grad()(0, 0) = 1.2;
+    s.get_grad()(0, 0) *= 1.5 + 1.2;
+    probe(i) = s.get_grad()(0, 0);
+  };
+
+  Kokkos::parallel_for("loop", nelems, loop_body);
+  Kokkos::fence();
+
+  for (int i = 0; i < 10; i++) {
+    std::printf("probe[%2d]: %.10f\n", i, probe(i));
+  }
+}
+
 int main(int argc, char* argv[]) {
   Kokkos::initialize();
   {  // test_axpy(argc, argv);
-     // test_matvec(argc, argv);
-     // test_unordered_set();
-     // test_subview();
-     // test_sort();
-     // test_is_same_layout();
-     // test_complex();
-     // test_is_complex();
-     // test_view_is_allocated();
-     // test_smart_pointer_behavior();
-     // test_copy();
-     // test_parallel_for();
+    // test_matvec(argc, argv);
+    // test_unordered_set();
+    // test_subview();
+    // test_sort();
+    // test_is_same_layout();
+    // test_complex();
+    // test_is_complex();
+    // test_view_is_allocated();
+    // test_smart_pointer_behavior();
+    // test_copy();
+    // test_parallel_for();
     // test_modify_view_from_const_lambda();
-    test_cuda_axpy(argc, argv);
+    // test_cuda_axpy(argc, argv);
+    // subview();
+    // test_cuda_axpy_with_UVM(argc, argv);
+    // ParallelVector pv(10);
+    // pv.set_values(4.2);
+    // test_cuda_functor_pass_by_ref();
+    // test_KOKKOS_ENABLE_CXX17();
+    test_fespace();
   }
   Kokkos::finalize();
-}
\ No newline at end of file
+}
diff --git a/examples/parallel_element/parallel_element.cpp b/examples/parallel_element/parallel_element.cpp
index 8bd22f46..75ee4b7f 100644
--- a/examples/parallel_element/parallel_element.cpp
+++ b/examples/parallel_element/parallel_element.cpp
@@ -32,7 +32,7 @@ class TopoElasticityAnalysis {
   static constexpr int var_dim = spatial_dim;
   static constexpr int block_size = var_dim;
 
-  using Vec_t = SolutionVector<T>;
+  using Vec_t = MultiArrayNew<T *>;
   using BSRMat_t = BSRMat<T, block_size, block_size>;
 
   using Quadrature = QuadGaussQuadrature<order>;
@@ -74,9 +74,9 @@ class TopoElasticityAnalysis {
         mesh(conn),
         geomesh(conn),
         datamesh(conn),
-        sol(mesh.get_num_dof()),
-        geo(geomesh.get_num_dof()),
-        data(datamesh.get_num_dof()),
+        sol("sol", mesh.get_num_dof()),
+        geo("geo", geomesh.get_num_dof()),
+        data("data", datamesh.get_num_dof()),
         elem_sol(mesh, sol),
         elem_geo(geomesh, geo),
         elem_data(datamesh, data) {
@@ -111,7 +111,7 @@ class TopoElasticityAnalysis {
       N = 20;
     }
 
-    Vec_t res(mesh.get_num_dof());
+    Vec_t res("res", mesh.get_num_dof());
     ElemVec elem_res(mesh, res);
     StopWatch watch;
 
@@ -166,7 +166,7 @@ int main(int argc, char *argv[]) {
     MesherRect2D mesher(nx, ny, lx, ly);
     bool randomize = true;
     unsigned int seed = 0;
-    double fraction = 0.2;
+    double fraction = 0.1;
     mesher.set_X_conn<index_t, double>(Xloc.data(), quad.data(), randomize,
                                        seed, fraction);
 
diff --git a/examples/spectral/spectral.cpp b/examples/spectral/spectral.cpp
index 84e3d295..58aaf9f7 100644
--- a/examples/spectral/spectral.cpp
+++ b/examples/spectral/spectral.cpp
@@ -109,7 +109,7 @@ void main_body(std::string type, int ar = 1, double h = 1.0,
   using Quadrature = HexGaussQuadrature<degree + 1>;
   using DataBasis = FEBasis<T>;
   using GeoBasis = FEBasis<T, LagrangeH1HexBasis<T, spatial_dim, degree>>;
-  using DataElemVec = A2D::ElementVector_Serial<T, DataBasis, BasisVecType>;
+  using DataElemVec = A2D::ElementVector_Empty;
   using GeoElemVec = A2D::ElementVector_Serial<T, GeoBasis, BasisVecType>;
   using ElemVec = A2D::ElementVector_Serial<T, Basis, BasisVecType>;
   using FE =
@@ -119,8 +119,7 @@ void main_body(std::string type, int ar = 1, double h = 1.0,
   using LOrderDataBasis = FEBasis<T>;
   using LOrderGeoBasis =
       FEBasis<T, LagrangeH1HexBasis<T, spatial_dim, low_degree>>;
-  using LOrderDataElemVec =
-      A2D::ElementVector_Serial<T, LOrderDataBasis, BasisVecType>;
+  using LOrderDataElemVec = A2D::ElementVector_Empty;
   using LOrderGeoElemVec =
       A2D::ElementVector_Serial<T, LOrderGeoBasis, BasisVecType>;
   using LOrderElemVec = A2D::ElementVector_Serial<T, LOrderBasis, BasisVecType>;
@@ -193,24 +192,21 @@ void main_body(std::string type, int ar = 1, double h = 1.0,
 
   ElementMesh<Basis> mesh(conn);
   ElementMesh<GeoBasis> geomesh(conn);
-  ElementMesh<DataBasis> datamesh(conn);
 
   ElementMesh<LOrderBasis> lorder_mesh(mesh);
   ElementMesh<LOrderGeoBasis> lorder_geomesh(geomesh);
-  ElementMesh<LOrderDataBasis> lorder_datamesh(datamesh);
 
   SolutionVector<T> sol(mesh.get_num_dof());
   SolutionVector<T> geo(geomesh.get_num_dof());
-  SolutionVector<T> data(datamesh.get_num_dof());
 
-  DataElemVec elem_data(datamesh, data);
+  DataElemVec elem_data;
   GeoElemVec elem_geo(geomesh, geo);
   ElemVec elem_sol(mesh, sol);
 
   // Set the geometry from the node locations
   set_geo_from_hex_nodes<GeoBasis>(nhex, hex, Xloc, elem_geo);
 
-  LOrderDataElemVec lorder_elem_data(lorder_datamesh, data);
+  LOrderDataElemVec lorder_elem_data;
   LOrderGeoElemVec lorder_elem_geo(lorder_geomesh, geo);
   LOrderElemVec lorder_elem_sol(lorder_mesh, sol);
 
diff --git a/examples/sphere/sphere.cpp b/examples/sphere/sphere.cpp
index 121b2925..e187ed6e 100644
--- a/examples/sphere/sphere.cpp
+++ b/examples/sphere/sphere.cpp
@@ -24,10 +24,11 @@ void set_geo_spherical(const T alpha, const T R, GeoElemVec& elem_geo) {
   T ulen = 2.0 / nx;
 
   T Xloc[3 * ET::HEX_NVERTS];
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
-    T u = 1.0 * ET::HEX_VERTS_CART[ii][0];
-    T v = 1.0 * ET::HEX_VERTS_CART[ii][1];
-    T w = 1.0 * ET::HEX_VERTS_CART[ii][2];
+    T u = 1.0 * HEX_VERTS_CART[ii][0];
+    T v = 1.0 * HEX_VERTS_CART[ii][1];
+    T w = 1.0 * HEX_VERTS_CART[ii][2];
 
     Xloc[3 * ii] = (2.0 * u - 1.0) * a;
     Xloc[3 * ii + 1] = (2.0 * v - 1.0) * a;
@@ -64,6 +65,7 @@ void set_geo_spherical(const T alpha, const T R, GeoElemVec& elem_geo) {
   }
 
   // We know the 2-direction will be radial...
+  auto HEX_BOUND_VERTS = ET::get_hex_bound_verts();
   for (index_t face = 0; face < ET::HEX_NBOUNDS; face++) {
     for (index_t k = 0; k < nx; k++) {
       for (index_t j = 0; j < nx; j++) {
@@ -91,9 +93,9 @@ void set_geo_spherical(const T alpha, const T R, GeoElemVec& elem_geo) {
             T z0 = 0.0;
 
             for (index_t jj = 0; jj < 4; jj++) {
-              x0 += N[jj] * Xloc[3 * ET::HEX_BOUND_VERTS[face][jj]];
-              y0 += N[jj] * Xloc[3 * ET::HEX_BOUND_VERTS[face][jj] + 1];
-              z0 += N[jj] * Xloc[3 * ET::HEX_BOUND_VERTS[face][jj] + 2];
+              x0 += N[jj] * Xloc[3 * HEX_BOUND_VERTS[face][jj]];
+              y0 += N[jj] * Xloc[3 * HEX_BOUND_VERTS[face][jj] + 1];
+              z0 += N[jj] * Xloc[3 * HEX_BOUND_VERTS[face][jj] + 2];
             }
             T norm = std::sqrt(x0 * x0 + y0 * y0 + z0 * z0);
 
@@ -546,11 +548,12 @@ void find_spherical_error(bool write_sphere = false) {
   }
 
   // For each face add the faces
+  auto HEX_BOUND_VERTS = ET::get_hex_bound_verts();
   for (index_t face = 0; face < ET::HEX_NBOUNDS; face++) {
     for (index_t i = 0; i < 4; i++) {
-      hex[(face + 1) * ET::HEX_NVERTS + i] = hex[ET::HEX_BOUND_VERTS[face][i]];
+      hex[(face + 1) * ET::HEX_NVERTS + i] = hex[HEX_BOUND_VERTS[face][i]];
       hex[(face + 1) * ET::HEX_NVERTS + i + 4] =
-          hex_ext[ET::HEX_BOUND_VERTS[face][i]];
+          hex_ext[HEX_BOUND_VERTS[face][i]];
     }
   }
 
diff --git a/examples/test_verify/verify_L2element.cpp b/examples/test_verify/verify_L2element.cpp
index 4206ffae..04168849 100644
--- a/examples/test_verify/verify_L2element.cpp
+++ b/examples/test_verify/verify_L2element.cpp
@@ -91,13 +91,14 @@ void main_body() {
   index_t *tets = NULL, *wedge = NULL, *pyrmd = NULL;
   index_t hex[8 * nhex];
 
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (index_t k = 0, e = 0; k < nz; k++) {
     for (index_t j = 0; j < ny; j++) {
       for (index_t i = 0; i < nx; i++, e++) {
         for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          hex[8 * e + ii] = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                     j + ET::HEX_VERTS_CART[ii][1],
-                                     k + ET::HEX_VERTS_CART[ii][2]);
+          hex[8 * e + ii] =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
         }
       }
     }
@@ -173,9 +174,9 @@ void main_body() {
 
         // Get the geometry values
         for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          index_t node = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                  j + ET::HEX_VERTS_CART[ii][1],
-                                  k + ET::HEX_VERTS_CART[ii][2]);
+          index_t node =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
 
           // Set the entity DOF
           index_t basis = 0;
diff --git a/examples/test_verify/verify_element.cpp b/examples/test_verify/verify_element.cpp
index eacc2677..4be877a9 100644
--- a/examples/test_verify/verify_element.cpp
+++ b/examples/test_verify/verify_element.cpp
@@ -89,13 +89,14 @@ void main_body() {
   index_t *tets = NULL, *wedge = NULL, *pyrmd = NULL;
   index_t hex[8 * nhex];
 
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (index_t k = 0, e = 0; k < nz; k++) {
     for (index_t j = 0; j < ny; j++) {
       for (index_t i = 0; i < nx; i++, e++) {
         for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          hex[8 * e + ii] = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                     j + ET::HEX_VERTS_CART[ii][1],
-                                     k + ET::HEX_VERTS_CART[ii][2]);
+          hex[8 * e + ii] =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
         }
       }
     }
@@ -167,9 +168,9 @@ void main_body() {
 
         // Get the geometry values
         for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          index_t node = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                  j + ET::HEX_VERTS_CART[ii][1],
-                                  k + ET::HEX_VERTS_CART[ii][2]);
+          index_t node =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
 
           // Set the entity DOF
           index_t basis = 0;
diff --git a/examples/topology_optimization/topology_elasticity.cpp b/examples/topology_optimization/topology_elasticity.cpp
index ef77ce66..636aff1b 100644
--- a/examples/topology_optimization/topology_elasticity.cpp
+++ b/examples/topology_optimization/topology_elasticity.cpp
@@ -176,13 +176,14 @@ create_analysis_box(std::string prefix, double vb_traction_frac, int b_nx,
 
   // Populate hex
   using ET = ElementTypes;
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (int k = 0, e = 0; k < b_nz; k++) {
     for (int j = 0; j < b_ny; j++) {
       for (int i = 0; i < b_nx; i++, e++) {
         for (int ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          hex[8 * e + ii] = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                     j + ET::HEX_VERTS_CART[ii][1],
-                                     k + ET::HEX_VERTS_CART[ii][2]);
+          hex[8 * e + ii] =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
         }
       }
     }
diff --git a/examples/topology_optimization/topology_heat.cpp b/examples/topology_optimization/topology_heat.cpp
index 3146dbc6..4041e0c9 100644
--- a/examples/topology_optimization/topology_heat.cpp
+++ b/examples/topology_optimization/topology_heat.cpp
@@ -50,13 +50,14 @@ void main_body(std::string prefix, double bc_temp, double heat_source,
 
   // Populate hex
   using ET = A2D::ElementTypes;
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (int k = 0, e = 0; k < nz; k++) {
     for (int j = 0; j < ny; j++) {
       for (int i = 0; i < nx; i++, e++) {
         for (int ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          hex[8 * e + ii] = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                     j + ET::HEX_VERTS_CART[ii][1],
-                                     k + ET::HEX_VERTS_CART[ii][2]);
+          hex[8 * e + ii] =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
         }
       }
     }
diff --git a/examples/topology_optimization/toy_grad_check.cpp b/examples/topology_optimization/toy_grad_check.cpp
index 990aebc5..49f570f5 100644
--- a/examples/topology_optimization/toy_grad_check.cpp
+++ b/examples/topology_optimization/toy_grad_check.cpp
@@ -32,14 +32,14 @@ int main(int argc, char *argv[]) {
   int hex[8 * nhex];
 
   using ET = A2D::ElementTypes;
-
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   for (int k = 0, e = 0; k < nz; k++) {
     for (int j = 0; j < ny; j++) {
       for (int i = 0; i < nx; i++, e++) {
         for (int ii = 0; ii < ET::HEX_NVERTS; ii++) {
-          hex[8 * e + ii] = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                     j + ET::HEX_VERTS_CART[ii][1],
-                                     k + ET::HEX_VERTS_CART[ii][2]);
+          hex[8 * e + ii] =
+              node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                       k + HEX_VERTS_CART[ii][2]);
         }
       }
     }
diff --git a/include/multiphysics/febasis.h b/include/multiphysics/febasis.h
index a0998a3c..8a051d46 100644
--- a/include/multiphysics/febasis.h
+++ b/include/multiphysics/febasis.h
@@ -443,14 +443,16 @@ class FEBasis {
    * @param orient Orientation flag indicating the relative orientation
    * @param elem_signs Degrees of freedom for this element
    */
+  template <class ElemSign>
   KOKKOS_FUNCTION static void set_entity_signs(index_t basis,
                                                ET::ElementEntity entity,
                                                index_t index, index_t orient,
-                                               int elem_signs[]) {
+                                               ElemSign& elem_signs) {
     if constexpr (sizeof...(Basis) == 0) {
       return;
     } else {
-      set_entity_signs<0, Basis...>(basis, entity, index, orient, elem_signs);
+      set_entity_signs<0, ElemSign, Basis...>(basis, entity, index, orient,
+                                              elem_signs);
     }
   }
 
@@ -492,11 +494,13 @@ class FEBasis {
    * @param horder_signs The signs from the high-order element
    * @param lorder_signs The output signs for the low-order element
    */
+  template <class HOrderSign, class LOrderSign>
   KOKKOS_FUNCTION static void get_lorder_signs(const index_t n,
-                                               const int horder_signs[],
-                                               int lorder_signs[]) {
+                                               const HOrderSign& horder_signs,
+                                               LOrderSign& lorder_signs) {
     if constexpr (sizeof...(Basis) > 0) {
-      get_lorder_signs_<0, 0, Basis...>(n, horder_signs, lorder_signs);
+      get_lorder_signs_<0, 0, HOrderSign, LOrderSign, Basis...>(n, horder_signs,
+                                                                lorder_signs);
     }
   }
 
@@ -783,21 +787,21 @@ class FEBasis {
     }
   }
 
-  template <index_t r, class First, class... Remain>
+  template <index_t r, class ElemSign, class First, class... Remain>
   KOKKOS_FUNCTION static void set_entity_signs(index_t basis,
                                                ET::ElementEntity entity,
                                                index_t index, index_t orient,
-                                               int elem_signs[]) {
+                                               ElemSign& elem_signs) {
     if (basis == r) {
-      First::template set_entity_signs<get_dof_offset<r>()>(entity, index,
-                                                            orient, elem_signs);
+      First::template set_entity_signs<get_dof_offset<r>(), ElemSign>(
+          entity, index, orient, elem_signs);
     }
     if constexpr (sizeof...(Remain) == 0) {
-      First::template set_entity_signs<get_dof_offset<r>()>(entity, index,
-                                                            orient, elem_signs);
+      First::template set_entity_signs<get_dof_offset<r>(), ElemSign>(
+          entity, index, orient, elem_signs);
     } else {
-      set_entity_signs<r + 1, Remain...>(basis, entity, index, orient,
-                                         elem_signs);
+      set_entity_signs<r + 1, ElemSign, Remain...>(basis, entity, index, orient,
+                                                   elem_signs);
     }
   }
 
@@ -874,16 +878,17 @@ class FEBasis {
     }
   }
 
-  template <index_t hoffset, index_t loffset, class First, class... Remain>
+  template <index_t hoffset, index_t loffset, class HOrderSign,
+            class LOrderSign, class First, class... Remain>
   KOKKOS_FUNCTION static void get_lorder_signs_(const index_t n,
-                                                const int horder_signs[],
-                                                int lorder_signs[]) {
-    First::template get_lorder_signs<hoffset, loffset>(n, horder_signs,
-                                                       lorder_signs);
+                                                const HOrderSign& horder_signs,
+                                                LOrderSign& lorder_signs) {
+    First::template get_lorder_signs<hoffset, loffset, HOrderSign, LOrderSign>(
+        n, horder_signs, lorder_signs);
     if constexpr (sizeof...(Remain) > 0) {
       get_lorder_signs_<hoffset + First::ndof,
-                        loffset + First::LOrderBasis::ndof, Remain...>(
-          n, horder_signs, lorder_signs);
+                        loffset + First::LOrderBasis::ndof, HOrderSign,
+                        LOrderSign, Remain...>(n, horder_signs, lorder_signs);
     }
   }
 };
diff --git a/include/multiphysics/feelementtypes-inl.h b/include/multiphysics/feelementtypes-inl.h
index b5953d82..6a9e9cef 100644
--- a/include/multiphysics/feelementtypes-inl.h
+++ b/include/multiphysics/feelementtypes-inl.h
@@ -19,8 +19,9 @@ namespace A2D {
  */
 template <index_t offset, index_t ndof, index_t nx, class ElemDof,
           class EntityDof>
-void ElementTypes::get_line_bound_dof(index_t b, const ElemDof& element,
-                                      EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_line_bound_dof(index_t b,
+                                                      const ElemDof& element,
+                                                      EntityDof& entity) {
   const index_t start = offset + ndof * (nx - 1) * b;
   for (index_t i = 0; i < ndof; i++) {
     entity[i] = element[start + i];
@@ -41,9 +42,10 @@ void ElementTypes::get_line_bound_dof(index_t b, const ElemDof& element,
  */
 template <index_t offset, index_t ndof, index_t nx, class EntityDof,
           class ElemDof>
-void ElementTypes::set_line_bound_dof(index_t b, const index_t orient,
-                                      const EntityDof& entity,
-                                      ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_line_bound_dof(index_t b,
+                                                      const index_t orient,
+                                                      const EntityDof& entity,
+                                                      ElemDof& element) {
   const index_t start = offset + ndof * (nx - 1) * ((b + orient) % 2);
   for (index_t i = 0; i < ndof; i++) {
     element[start + i] = entity[i];
@@ -64,8 +66,8 @@ void ElementTypes::set_line_bound_dof(index_t b, const index_t orient,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, class ElemDof,
           class EntityDof>
-void ElementTypes::get_line_domain_dof(const ElemDof& element,
-                                       EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_line_domain_dof(const ElemDof& element,
+                                                       EntityDof& entity) {
   if constexpr (ends) {
     for (index_t j = 0; j < nx; j++) {
       for (index_t i = 0; i < ndof; i++) {
@@ -95,9 +97,9 @@ void ElementTypes::get_line_domain_dof(const ElemDof& element,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, class EntityDof,
           class ElemDof>
-void ElementTypes::set_line_domain_dof(const index_t orient,
-                                       const EntityDof& entity,
-                                       ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_line_domain_dof(const index_t orient,
+                                                       const EntityDof& entity,
+                                                       ElemDof& element) {
   index_t index;
   if constexpr (ends) {
     for (index_t j = 0; j < nx; j++) {
@@ -129,8 +131,12 @@ void ElementTypes::set_line_domain_dof(const index_t orient,
  * @param domain_verts
  * @return index_t
  */
-inline index_t ElementTypes::get_quad_domain_orientation(
+KOKKOS_FUNCTION inline index_t ElementTypes::get_quad_domain_orientation(
     const index_t ref_domain_verts[], const index_t domain_verts[]) {
+  static const index_t QUAD_DOMAIN_ORIENTATIONS[][4] = {
+      {0, 1, 2, 3}, {3, 0, 1, 2}, {2, 3, 0, 1}, {1, 2, 3, 0},
+      {0, 3, 2, 1}, {3, 2, 1, 0}, {2, 1, 0, 3}, {1, 0, 3, 2}};
+
   index_t orient = 0;
   for (; orient < NUM_QUAD_DOMAIN_ORIENTATIONS; orient++) {
     if (ref_domain_verts[0] ==
@@ -158,7 +164,7 @@ inline index_t ElementTypes::get_quad_domain_orientation(
  * @param u The 1st coordinate on the reference face
  * @param v The 2nd coordinate on the reference face
  */
-inline void ElementTypes::get_coords_on_quad_ref_element(
+KOKKOS_FUNCTION inline void ElementTypes::get_coords_on_quad_ref_element(
     const index_t orient, const index_t hx, const index_t hy, const index_t x,
     const index_t y, index_t* u, index_t* v) {
   if (orient == 0) {
@@ -191,11 +197,9 @@ inline void ElementTypes::get_coords_on_quad_ref_element(
   }
 }
 
-inline index_t ElementTypes::get_index_on_quad_ref_element(const index_t orient,
-                                                           const index_t hx,
-                                                           const index_t hy,
-                                                           const index_t x,
-                                                           const index_t y) {
+KOKKOS_FUNCTION inline index_t ElementTypes::get_index_on_quad_ref_element(
+    const index_t orient, const index_t hx, const index_t hy, const index_t x,
+    const index_t y) {
   if (orient == 0) {
     // *u = x;
     // *v = y;
@@ -243,7 +247,7 @@ inline index_t ElementTypes::get_index_on_quad_ref_element(const index_t orient,
  * @param j Index along the y-direction
  */
 template <index_t nx, index_t ny>
-int ElementTypes::get_quad_node(const int i, const int j) {
+KOKKOS_FUNCTION int ElementTypes::get_quad_node(const int i, const int j) {
   return i + nx * j;
 }
 
@@ -259,8 +263,9 @@ int ElementTypes::get_quad_node(const int i, const int j) {
  * @return The number of nodes along the edge
  */
 template <index_t nx, index_t ny>
-index_t ElementTypes::get_quad_bound_length(const index_t v0,
-                                            const index_t v1) {
+KOKKOS_FUNCTION index_t ElementTypes::get_quad_bound_length(const index_t v0,
+                                                            const index_t v1) {
+  auto QUAD_VERTS_CART = get_quad_verts_cart();
   if (QUAD_VERTS_CART[v0][0] != QUAD_VERTS_CART[v1][0]) {
     return nx;
   } else {
@@ -283,8 +288,10 @@ index_t ElementTypes::get_quad_bound_length(const index_t v0,
  */
 template <index_t offset, index_t ndof, index_t nx, index_t ny, class ElemDof,
           class EntityDof>
-void ElementTypes::get_quad_vert_dof(index_t v, const ElemDof& element,
-                                     EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_quad_vert_dof(index_t v,
+                                                     const ElemDof& element,
+                                                     EntityDof& entity) {
+  auto QUAD_VERTS_CART = get_quad_verts_cart();
   const index_t node = get_quad_node<nx, ny>((nx - 1) * QUAD_VERTS_CART[v][0],
                                              (ny - 1) * QUAD_VERTS_CART[v][1]);
 
@@ -309,8 +316,10 @@ void ElementTypes::get_quad_vert_dof(index_t v, const ElemDof& element,
  */
 template <index_t offset, index_t ndof, index_t nx, index_t ny, class EntityDof,
           class ElemDof>
-void ElementTypes::set_quad_vert_dof(index_t v, const EntityDof& entity,
-                                     ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_quad_vert_dof(index_t v,
+                                                     const EntityDof& entity,
+                                                     ElemDof& element) {
+  auto QUAD_VERTS_CART = get_quad_verts_cart();
   const index_t node = get_quad_node<nx, ny>((nx - 1) * QUAD_VERTS_CART[v][0],
                                              (ny - 1) * QUAD_VERTS_CART[v][1]);
 
@@ -336,11 +345,14 @@ void ElementTypes::set_quad_vert_dof(index_t v, const EntityDof& entity,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           class ElemDof, class EntityDof>
-void ElementTypes::get_quad_bound_dof(const index_t b, const ElemDof& element,
-                                      EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_quad_bound_dof(const index_t b,
+                                                      const ElemDof& element,
+                                                      EntityDof& entity) {
+  auto QUAD_VERTS_CART = get_quad_verts_cart();
+
   // Get the first and last vertices on the edge
-  const index_t v0 = QUAD_BOUND_VERTS[b][0];
-  const index_t v1 = QUAD_BOUND_VERTS[b][1];
+  const index_t v0 = get_quad_bound_verts()[b][0];
+  const index_t v1 = get_quad_bound_verts()[b][1];
 
   // Get the starting index on the element
   const index_t start = get_quad_node<nx, ny>(
@@ -397,12 +409,15 @@ void ElementTypes::get_quad_bound_dof(const index_t b, const ElemDof& element,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           class EntityDof, class ElemDof>
-void ElementTypes::set_quad_bound_dof(const index_t b, const index_t orient,
-                                      const EntityDof& entity,
-                                      ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_quad_bound_dof(const index_t b,
+                                                      const index_t orient,
+                                                      const EntityDof& entity,
+                                                      ElemDof& element) {
+  auto QUAD_VERTS_CART = get_quad_verts_cart();
+
   // Get the first and last vertices on the edge
-  const index_t v0 = QUAD_BOUND_VERTS[b][orient];
-  const index_t v1 = QUAD_BOUND_VERTS[b][(orient + 1) % 2];
+  const index_t v0 = get_quad_bound_verts()[b][orient];
+  const index_t v1 = get_quad_bound_verts()[b][(orient + 1) % 2];
 
   // Get the starting index on the element
   const index_t start = get_quad_node<nx, ny>(
@@ -454,8 +469,8 @@ void ElementTypes::set_quad_bound_dof(const index_t b, const index_t orient,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           class ElemDof, class EntityDof>
-void ElementTypes::get_quad_domain_dof(const ElemDof& element,
-                                       EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_quad_domain_dof(const ElemDof& element,
+                                                       EntityDof& entity) {
   if constexpr (ends) {
     for (index_t v = 0; v < ny; v++) {
       for (index_t u = 0; u < nx; u++) {
@@ -498,9 +513,9 @@ void ElementTypes::get_quad_domain_dof(const ElemDof& element,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           class EntityDof, class ElemDof>
-void ElementTypes::set_quad_domain_dof(const index_t orient,
-                                       const EntityDof& entity,
-                                       ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_quad_domain_dof(const index_t orient,
+                                                       const EntityDof& entity,
+                                                       ElemDof& element) {
   if constexpr (ends) {
     for (index_t v = 0; v < ny; v++) {
       for (index_t u = 0; u < nx; u++) {
@@ -544,7 +559,8 @@ void ElementTypes::set_quad_domain_dof(const index_t orient,
  * @param k Index along the z-direction
  */
 template <index_t nx, index_t ny, index_t nz>
-int ElementTypes::get_hex_node(const int i, const int j, const int k) {
+KOKKOS_FUNCTION int ElementTypes::get_hex_node(const int i, const int j,
+                                               const int k) {
   return i + nx * (j + ny * k);
 }
 
@@ -561,7 +577,9 @@ int ElementTypes::get_hex_node(const int i, const int j, const int k) {
  * @return The number of nodes along the edge
  */
 template <index_t nx, index_t ny, index_t nz>
-index_t ElementTypes::get_hex_edge_length(const index_t v0, const index_t v1) {
+KOKKOS_FUNCTION index_t ElementTypes::get_hex_edge_length(const index_t v0,
+                                                          const index_t v1) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
   if (HEX_VERTS_CART[v0][0] != HEX_VERTS_CART[v1][0]) {
     return nx;
   } else if (HEX_VERTS_CART[v0][1] != HEX_VERTS_CART[v1][1]) {
@@ -587,8 +605,10 @@ index_t ElementTypes::get_hex_edge_length(const index_t v0, const index_t v1) {
  */
 template <index_t offset, index_t ndof, index_t nx, index_t ny, index_t nz,
           class ElemDof, class EntityDof>
-void ElementTypes::get_hex_vert_dof(index_t v, const ElemDof& element,
-                                    EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_hex_vert_dof(index_t v,
+                                                    const ElemDof& element,
+                                                    EntityDof& entity) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
   const index_t node = get_hex_node<nx, ny, nz>(
       (nx - 1) * HEX_VERTS_CART[v][0], (ny - 1) * HEX_VERTS_CART[v][1],
       (nz - 1) * HEX_VERTS_CART[v][2]);
@@ -615,8 +635,10 @@ void ElementTypes::get_hex_vert_dof(index_t v, const ElemDof& element,
  */
 template <index_t offset, index_t ndof, index_t nx, index_t ny, index_t nz,
           class EntityDof, class ElemDof>
-void ElementTypes::set_hex_vert_dof(index_t v, const EntityDof& entity,
-                                    ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_hex_vert_dof(index_t v,
+                                                    const EntityDof& entity,
+                                                    ElemDof& element) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
   const index_t node = get_hex_node<nx, ny, nz>(
       (nx - 1) * HEX_VERTS_CART[v][0], (ny - 1) * HEX_VERTS_CART[v][1],
       (nz - 1) * HEX_VERTS_CART[v][2]);
@@ -644,11 +666,14 @@ void ElementTypes::set_hex_vert_dof(index_t v, const EntityDof& entity,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           index_t nz, class ElemDof, class EntityDof>
-void ElementTypes::get_hex_edge_dof(const index_t e, const ElemDof& element,
-                                    EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_hex_edge_dof(const index_t e,
+                                                    const ElemDof& element,
+                                                    EntityDof& entity) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
+
   // Get the first and last vertices on the edge
-  const index_t v0 = HEX_EDGE_VERTS[e][0];
-  const index_t v1 = HEX_EDGE_VERTS[e][1];
+  const index_t v0 = get_hex_edge_verts()[e][0];
+  const index_t v1 = get_hex_edge_verts()[e][1];
 
   // Get the starting index on the element
   const index_t start = get_hex_node<nx, ny, nz>(
@@ -708,11 +733,15 @@ void ElementTypes::get_hex_edge_dof(const index_t e, const ElemDof& element,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           index_t nz, class EntityDof, class ElemDof>
-void ElementTypes::set_hex_edge_dof(const index_t e, const index_t orient,
-                                    const EntityDof& entity, ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_hex_edge_dof(const index_t e,
+                                                    const index_t orient,
+                                                    const EntityDof& entity,
+                                                    ElemDof& element) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
+
   // Get the first and last vertices on the edge
-  const index_t v0 = HEX_EDGE_VERTS[e][orient];
-  const index_t v1 = HEX_EDGE_VERTS[e][(orient + 1) % 2];
+  const index_t v0 = get_hex_edge_verts()[e][orient];
+  const index_t v1 = get_hex_edge_verts()[e][(orient + 1) % 2];
 
   // Get the starting index on the element
   const index_t start = get_hex_node<nx, ny, nz>(
@@ -768,12 +797,17 @@ void ElementTypes::set_hex_edge_dof(const index_t e, const index_t orient,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           index_t nz, class ElemDof, class EntityDof>
-void ElementTypes::get_hex_bound_dof(const index_t b, const ElemDof& element,
-                                     EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_hex_bound_dof(const index_t b,
+                                                     const ElemDof& element,
+                                                     EntityDof& entity) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
+
   // Get the origin and bound vert directions
-  const index_t v0 = HEX_BOUND_VERTS[b][0];  // Root vertex
-  const index_t v1 = HEX_BOUND_VERTS[b][1];  // Vertex along the bound u dir
-  const index_t v3 = HEX_BOUND_VERTS[b][3];  // Vertex along the bound v dir
+  const index_t v0 = get_hex_bound_verts()[b][0];  // Root vertex
+  const index_t v1 =
+      get_hex_bound_verts()[b][1];  // Vertex along the bound u dir
+  const index_t v3 =
+      get_hex_bound_verts()[b][3];  // Vertex along the bound v dir
 
   // Get the root node location
   const index_t start = get_hex_node<nx, ny, nz>(
@@ -842,13 +876,18 @@ void ElementTypes::get_hex_bound_dof(const index_t b, const ElemDof& element,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           index_t nz, class EntityDof, class ElemDof>
-void ElementTypes::set_hex_bound_dof(const index_t b, const index_t orient,
-                                     const EntityDof& entity,
-                                     ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_hex_bound_dof(const index_t b,
+                                                     const index_t orient,
+                                                     const EntityDof& entity,
+                                                     ElemDof& element) {
+  auto HEX_VERTS_CART = get_hex_verts_cart();
+
   // Get the origin and bound vert directions
-  const index_t v0 = HEX_BOUND_VERTS[b][0];  // Root vertex
-  const index_t v1 = HEX_BOUND_VERTS[b][1];  // Vertex along the bound u dir
-  const index_t v3 = HEX_BOUND_VERTS[b][3];  // Vertex along the bound v dir
+  const index_t v0 = get_hex_bound_verts()[b][0];  // Root vertex
+  const index_t v1 =
+      get_hex_bound_verts()[b][1];  // Vertex along the bound u dir
+  const index_t v3 =
+      get_hex_bound_verts()[b][3];  // Vertex along the bound v dir
 
   // Get the root node location
   const index_t start = get_hex_node<nx, ny, nz>(
@@ -918,8 +957,8 @@ void ElementTypes::set_hex_bound_dof(const index_t b, const index_t orient,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           index_t nz, class ElemDof, class EntityDof>
-void ElementTypes::get_hex_domain_dof(const ElemDof& element,
-                                      EntityDof& entity) {
+KOKKOS_FUNCTION void ElementTypes::get_hex_domain_dof(const ElemDof& element,
+                                                      EntityDof& entity) {
   if constexpr (ends) {
     for (index_t w = 0; w < nz; w++) {
       for (index_t v = 0; v < ny; v++) {
@@ -967,8 +1006,8 @@ void ElementTypes::get_hex_domain_dof(const ElemDof& element,
  */
 template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
           index_t nz, class EntityDof, class ElemDof>
-void ElementTypes::set_hex_domain_dof(const EntityDof& entity,
-                                      ElemDof& element) {
+KOKKOS_FUNCTION void ElementTypes::set_hex_domain_dof(const EntityDof& entity,
+                                                      ElemDof& element) {
   if constexpr (ends) {
     for (index_t w = 0; w < nz; w++) {
       for (index_t v = 0; v < ny; v++) {
diff --git a/include/multiphysics/feelementtypes.h b/include/multiphysics/feelementtypes.h
index 852ac2ea..cbaef695 100644
--- a/include/multiphysics/feelementtypes.h
+++ b/include/multiphysics/feelementtypes.h
@@ -72,25 +72,30 @@ class ElementTypes {
   // Get the degrees of freedom associated with the bound
   template <index_t offset, index_t ndof, index_t nx, class ElemDof,
             class EntityDof>
-  static void get_line_bound_dof(index_t b, const ElemDof& element,
-                                 EntityDof& entity);
+  KOKKOS_FUNCTION static void get_line_bound_dof(index_t b,
+                                                 const ElemDof& element,
+                                                 EntityDof& entity);
 
   // Get the degrees of freedom associated with a bound
   template <index_t offset, index_t ndof, index_t nx, class EntityDof,
             class ElemDof>
-  static void set_line_bound_dof(index_t b, const index_t orient,
-                                 const EntityDof& entity, ElemDof& element);
+  KOKKOS_FUNCTION static void set_line_bound_dof(index_t b,
+                                                 const index_t orient,
+                                                 const EntityDof& entity,
+                                                 ElemDof& element);
 
   // Get the degrees of freedom from the domain
   template <index_t offset, bool ends, index_t ndof, index_t nx, class ElemDof,
             class EntityDof>
-  static void get_line_domain_dof(const ElemDof& element, EntityDof& entity);
+  KOKKOS_FUNCTION static void get_line_domain_dof(const ElemDof& element,
+                                                  EntityDof& entity);
 
   // Set the degrees of freedom from the domain
   template <index_t offset, bool ends, index_t ndof, index_t nx,
             class EntityDof, class ElemDof>
-  static void set_line_domain_dof(const index_t orient, const EntityDof& entity,
-                                  ElemDof& element);
+  KOKKOS_FUNCTION static void set_line_domain_dof(const index_t orient,
+                                                  const EntityDof& entity,
+                                                  ElemDof& element);
 
   /**
    * @brief Triangle element
@@ -121,11 +126,18 @@ class ElementTypes {
   static const index_t TRI_NVERTS = 3;
 
   // Given bound index, return number of vertices/vertex indices
-  static constexpr index_t TRI_BOUND_NVERTS[] = {2, 2, 2};
-  static constexpr index_t TRI_BOUND_VERTS[][MAX_BOUND_VERTS] = {
-      {1, 2, NO_INDEX, NO_INDEX},
-      {2, 0, NO_INDEX, NO_INDEX},
-      {0, 1, NO_INDEX, NO_INDEX}};
+  KOKKOS_FUNCTION static const index_t* get_tri_bound_nverts() {
+    static constexpr index_t TRI_BOUND_NVERTS[] = {2, 2, 2};
+    return TRI_BOUND_NVERTS;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_tri_bound_verts())[MAX_BOUND_VERTS] {
+    static constexpr index_t TRI_BOUND_VERTS[][MAX_BOUND_VERTS] = {
+        {1, 2, NO_INDEX, NO_INDEX},
+        {2, 0, NO_INDEX, NO_INDEX},
+        {0, 1, NO_INDEX, NO_INDEX}};
+    return TRI_BOUND_VERTS;
+  }
 
   /**
    * @brief Quadrilateral element
@@ -157,82 +169,93 @@ class ElementTypes {
   static const index_t QUAD_NVERTS = 4;
 
   // Given bound index, return number of vertices/vertex indices
-  static constexpr index_t QUAD_BOUND_NVERTS[] = {2, 2, 2, 2};
-  static constexpr index_t QUAD_BOUND_VERTS[][MAX_BOUND_VERTS] = {
-      {0, 1, NO_INDEX, NO_INDEX},
-      {3, 2, NO_INDEX, NO_INDEX},
-      {0, 3, NO_INDEX, NO_INDEX},
-      {1, 2, NO_INDEX, NO_INDEX}};
+  KOKKOS_FUNCTION static const index_t* get_quad_bound_nverts() {
+    static constexpr index_t QUAD_BOUND_NVERTS[] = {2, 2, 2, 2};
+    return QUAD_BOUND_NVERTS;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_quad_bound_verts())[MAX_BOUND_VERTS] {
+    static constexpr index_t QUAD_BOUND_VERTS[][MAX_BOUND_VERTS] = {
+        {0, 1, NO_INDEX, NO_INDEX},
+        {3, 2, NO_INDEX, NO_INDEX},
+        {0, 3, NO_INDEX, NO_INDEX},
+        {1, 2, NO_INDEX, NO_INDEX}};
+    return QUAD_BOUND_VERTS;
+  }
 
   // Cartesian coordinates of the vertices in the reference element
-  static constexpr index_t QUAD_VERTS_CART[][2] = {
-      {0, 0}, {1, 0}, {1, 1}, {0, 1}};
+  KOKKOS_FUNCTION static const index_t (*get_quad_verts_cart())[2] {
+    static constexpr index_t QUAD_VERTS_CART[][2] = {
+        {0, 0}, {1, 0}, {1, 1}, {0, 1}};
+    return QUAD_VERTS_CART;
+  }
 
   static const index_t NUM_QUAD_DOMAIN_ORIENTATIONS = 8;
-  static constexpr index_t QUAD_DOMAIN_ORIENTATIONS[][4] = {
-      {0, 1, 2, 3}, {3, 0, 1, 2}, {2, 3, 0, 1}, {1, 2, 3, 0},
-      {0, 3, 2, 1}, {3, 2, 1, 0}, {2, 1, 0, 3}, {1, 0, 3, 2}};
 
   // Given a reference domain, find the orientation
-  static index_t get_quad_domain_orientation(const index_t ref_domain_verts[],
-                                             const index_t domain_verts[]);
+  KOKKOS_FUNCTION static index_t get_quad_domain_orientation(
+      const index_t ref_domain_verts[], const index_t domain_verts[]);
 
   // Get the coords on quad ref element object
-  static void get_coords_on_quad_ref_element(const index_t orient,
-                                             const index_t hx, const index_t hy,
-                                             const index_t x, const index_t y,
-                                             index_t* u, index_t* v);
+  KOKKOS_FUNCTION static void get_coords_on_quad_ref_element(
+      const index_t orient, const index_t hx, const index_t hy, const index_t x,
+      const index_t y, index_t* u, index_t* v);
 
-  static index_t get_index_on_quad_ref_element(const index_t orient,
-                                               const index_t hx,
-                                               const index_t hy,
-                                               const index_t x,
-                                               const index_t y);
+  KOKKOS_FUNCTION static index_t get_index_on_quad_ref_element(
+      const index_t orient, const index_t hx, const index_t hy, const index_t x,
+      const index_t y);
 
   // Get the local node index (without offset) for a node on an element
   // with nx and ny nodes along the local directions
   template <index_t nx, index_t ny>
-  static inline int get_quad_node(const int i, const int j);
+  KOKKOS_FUNCTION static inline int get_quad_node(const int i, const int j);
 
   // Get the bound length between the vertices v1 and v2
   template <index_t nx, index_t ny>
-  static inline index_t get_quad_bound_length(const index_t v0,
-                                              const index_t v1);
+  KOKKOS_FUNCTION static inline index_t get_quad_bound_length(const index_t v0,
+                                                              const index_t v1);
 
   // Get the degrees of freedom associated with the vertex
   template <index_t offset, index_t ndof, index_t nx, index_t ny, class ElemDof,
             class EntityDof>
-  static void get_quad_vert_dof(index_t v, const ElemDof& element,
-                                EntityDof& entity);
+  KOKKOS_FUNCTION static void get_quad_vert_dof(index_t v,
+                                                const ElemDof& element,
+                                                EntityDof& entity);
 
   // Get the degrees of freedom associated with the vertex
   template <index_t offset, index_t ndof, index_t nx, index_t ny,
             class EntityDof, class ElemDof>
-  static void set_quad_vert_dof(index_t v, const EntityDof& entity,
-                                ElemDof& element);
+  KOKKOS_FUNCTION static void set_quad_vert_dof(index_t v,
+                                                const EntityDof& entity,
+                                                ElemDof& element);
 
   // Get the degrees of freedom from the bound
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             class ElemDof, class EntityDof>
-  static void get_quad_bound_dof(const index_t b, const ElemDof& element,
-                                 EntityDof& entity);
+  KOKKOS_FUNCTION static void get_quad_bound_dof(const index_t b,
+                                                 const ElemDof& element,
+                                                 EntityDof& entity);
 
   // Set the degrees of freedom from the bound
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             class EntityDof, class ElemDof>
-  static void set_quad_bound_dof(const index_t b, const index_t orient,
-                                 const EntityDof& entity, ElemDof& element);
+  KOKKOS_FUNCTION static void set_quad_bound_dof(const index_t b,
+                                                 const index_t orient,
+                                                 const EntityDof& entity,
+                                                 ElemDof& element);
 
   // Get the degrees of freedom from the quad domain
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             class ElemDof, class EntityDof>
-  static void get_quad_domain_dof(const ElemDof& element, EntityDof& entity);
+  KOKKOS_FUNCTION static void get_quad_domain_dof(const ElemDof& element,
+                                                  EntityDof& entity);
 
   // Set the degrees of freedom from the domain into the element
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             class EntityDof, class ElemDof>
-  static void set_quad_domain_dof(const index_t orient, const EntityDof& entity,
-                                  ElemDof& element);
+  KOKKOS_FUNCTION static void set_quad_domain_dof(const index_t orient,
+                                                  const EntityDof& entity,
+                                                  ElemDof& element);
 
   /**
    * @brief Tetrahedral properties
@@ -279,24 +302,41 @@ class ElementTypes {
   static const index_t TET_NVERTS = 4;
 
   // Given edge index, return edge vertex indices
-  static constexpr index_t TET_EDGE_VERTS[][2] = {{0, 1}, {1, 2}, {2, 0},
-                                                  {0, 3}, {1, 3}, {2, 3}};
+  KOKKOS_FUNCTION static const index_t (*get_tet_edge_verts())[2] {
+    static constexpr index_t TET_EDGE_VERTS[][2] = {{0, 1}, {1, 2}, {2, 0},
+                                                    {0, 3}, {1, 3}, {2, 3}};
+    return TET_EDGE_VERTS;
+  }
 
   // Given bounds index, return edge indices
-  static constexpr index_t TET_BOUND_NEDGES[] = {3, 3, 3, 3};
-  static constexpr index_t TET_BOUND_EDGES[][MAX_BOUND_EDGES] = {
-      {1, 5, 4, NO_INDEX},
-      {3, 5, 2, NO_INDEX},
-      {0, 4, 3, NO_INDEX},
-      {2, 1, 0, NO_INDEX}};
+  KOKKOS_FUNCTION static const index_t* get_tet_bound_nedges() {
+    static constexpr index_t TET_BOUND_NEDGES[] = {3, 3, 3, 3};
+    return TET_BOUND_NEDGES;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_tet_bound_edges())[MAX_BOUND_EDGES] {
+    static constexpr index_t TET_BOUND_EDGES[][MAX_BOUND_EDGES] = {
+        {1, 5, 4, NO_INDEX},
+        {3, 5, 2, NO_INDEX},
+        {0, 4, 3, NO_INDEX},
+        {2, 1, 0, NO_INDEX}};
+    return TET_BOUND_EDGES;
+  }
 
   // Given bounds index, return number of vertices/vertex indices
-  static constexpr index_t TET_BOUND_NVERTS[] = {3, 3, 3, 3};
-  static constexpr index_t TET_BOUND_VERTS[][MAX_BOUND_VERTS] = {
-      {1, 2, 3, NO_INDEX},
-      {0, 3, 2, NO_INDEX},
-      {0, 1, 3, NO_INDEX},
-      {0, 2, 1, NO_INDEX}};
+  KOKKOS_FUNCTION static const index_t* get_tet_bound_nverts() {
+    static constexpr index_t TET_BOUND_NVERTS[] = {3, 3, 3, 3};
+    return TET_BOUND_NVERTS;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_tet_bound_verts())[MAX_BOUND_VERTS] {
+    static constexpr index_t TET_BOUND_VERTS[][MAX_BOUND_VERTS] = {
+        {1, 2, 3, NO_INDEX},
+        {0, 3, 2, NO_INDEX},
+        {0, 1, 3, NO_INDEX},
+        {0, 2, 1, NO_INDEX}};
+    return TET_BOUND_VERTS;
+  }
 
   /**
    * @brief Hexahedral properties
@@ -349,81 +389,113 @@ class ElementTypes {
   static const index_t HEX_NVERTS = 8;
 
   // Given edge index, return edge vertex indices
-  static constexpr index_t HEX_EDGE_VERTS[][2] = {
-      {0, 1}, {3, 2}, {4, 5}, {7, 6}, {0, 3}, {1, 2},
-      {4, 7}, {5, 6}, {0, 4}, {1, 5}, {3, 7}, {2, 6}};
+  KOKKOS_FUNCTION static const index_t (*get_hex_edge_verts())[2] {
+    static constexpr index_t HEX_EDGE_VERTS[][2] = {
+        {0, 1}, {3, 2}, {4, 5}, {7, 6}, {0, 3}, {1, 2},
+        {4, 7}, {5, 6}, {0, 4}, {1, 5}, {3, 7}, {2, 6}};
+    return HEX_EDGE_VERTS;
+  }
 
   // Given bound index, return edge indices
-  static constexpr index_t HEX_BOUND_NEDGES[] = {4, 4, 4, 4, 4, 4};
-  static constexpr index_t HEX_BOUND_EDGES[][MAX_BOUND_EDGES] = {
-      {8, 10, 4, 6},  {5, 7, 9, 11}, {0, 2, 8, 9},
-      {1, 3, 11, 10}, {4, 5, 0, 1},  {2, 3, 6, 7}};
+  KOKKOS_FUNCTION static const index_t* get_hex_bound_nedges() {
+    static constexpr index_t HEX_BOUND_NEDGES[] = {4, 4, 4, 4, 4, 4};
+    return HEX_BOUND_NEDGES;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_hex_bound_edges())[MAX_BOUND_EDGES] {
+    static constexpr index_t HEX_BOUND_EDGES[][MAX_BOUND_EDGES] = {
+        {8, 10, 4, 6},  {5, 7, 9, 11}, {0, 2, 8, 9},
+        {1, 3, 11, 10}, {4, 5, 0, 1},  {2, 3, 6, 7}};
+    return HEX_BOUND_EDGES;
+  }
 
   // Given bound index, return number of vertices/vertex indices
-  static constexpr index_t HEX_BOUND_NVERTS[] = {4, 4, 4, 4, 4, 4};
-  static constexpr index_t HEX_BOUND_VERTS[][MAX_BOUND_VERTS] = {
-      {0, 4, 7, 3}, {1, 2, 6, 5}, {0, 1, 5, 4},
-      {3, 7, 6, 2}, {0, 3, 2, 1}, {4, 5, 6, 7}};
+  KOKKOS_FUNCTION static const index_t* get_hex_bound_nverts() {
+    static constexpr index_t HEX_BOUND_NVERTS[] = {4, 4, 4, 4, 4, 4};
+    return HEX_BOUND_NVERTS;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_hex_bound_verts())[MAX_BOUND_VERTS] {
+    static constexpr index_t HEX_BOUND_VERTS[][MAX_BOUND_VERTS] = {
+        {0, 4, 7, 3}, {1, 2, 6, 5}, {0, 1, 5, 4},
+        {3, 7, 6, 2}, {0, 3, 2, 1}, {4, 5, 6, 7}};
+    return HEX_BOUND_VERTS;
+  }
 
   // Cartesian coordinates of the vertices in the reference element
-  static constexpr index_t HEX_VERTS_CART[][3] = {
-      {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
-      {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
+  KOKKOS_FUNCTION static const index_t (*get_hex_verts_cart())[3] {
+    static constexpr index_t HEX_VERTS_CART[][3] = {
+        {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
+        {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
+    return HEX_VERTS_CART;
+  }
 
   // Get the local node index (without offset) for a node on an element
   // with nx, ny and nz nodes along the local directions
   template <index_t nx, index_t ny, index_t nz>
-  static int get_hex_node(const int i, const int j, const int k);
+  KOKKOS_FUNCTION static int get_hex_node(const int i, const int j,
+                                          const int k);
 
   // Get the edge length between the verties v1 and v2
   template <index_t nx, index_t ny, index_t nz>
-  static index_t get_hex_edge_length(const index_t v0, const index_t v1);
+  KOKKOS_FUNCTION static index_t get_hex_edge_length(const index_t v0,
+                                                     const index_t v1);
 
   // Get the degrees of freedom associated with the vertex
   template <index_t offset, index_t ndof, index_t nx, index_t ny, index_t nz,
             class ElemDof, class EntityDof>
-  static void get_hex_vert_dof(index_t v, const ElemDof& element,
-                               EntityDof& entity);
+  KOKKOS_FUNCTION static void get_hex_vert_dof(index_t v,
+                                               const ElemDof& element,
+                                               EntityDof& entity);
 
   // Get the degrees of freedom associated with the vertex
   template <index_t offset, index_t ndof, index_t nx, index_t ny, index_t nz,
             class EntityDof, class ElemDof>
-  static void set_hex_vert_dof(index_t v, const EntityDof& entity,
-                               ElemDof& element);
+  KOKKOS_FUNCTION static void set_hex_vert_dof(index_t v,
+                                               const EntityDof& entity,
+                                               ElemDof& element);
 
   // Get the degrees of freedom from the edge
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             index_t nz, class ElemDof, class EntityDof>
-  static void get_hex_edge_dof(const index_t e, const ElemDof& element,
-                               EntityDof& entity);
+  KOKKOS_FUNCTION static void get_hex_edge_dof(const index_t e,
+                                               const ElemDof& element,
+                                               EntityDof& entity);
 
   // Set the degrees of freedom from the edge
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             index_t nz, class EntityDof, class ElemDof>
-  static void set_hex_edge_dof(const index_t e, const index_t orient,
-                               const EntityDof& entity, ElemDof& element);
+  KOKKOS_FUNCTION static void set_hex_edge_dof(const index_t e,
+                                               const index_t orient,
+                                               const EntityDof& entity,
+                                               ElemDof& element);
 
   // Get the degrees of freedom from the hex bound
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             index_t nz, class ElemDof, class EntityDof>
-  static void get_hex_bound_dof(const index_t f, const ElemDof& element,
-                                EntityDof& entity);
+  KOKKOS_FUNCTION static void get_hex_bound_dof(const index_t f,
+                                                const ElemDof& element,
+                                                EntityDof& entity);
 
   // Set the degrees of freedom from the bound into the element
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             index_t nz, class EntityDof, class ElemDof>
-  static void set_hex_bound_dof(const index_t f, const index_t orient,
-                                const EntityDof& entity, ElemDof& element);
+  KOKKOS_FUNCTION static void set_hex_bound_dof(const index_t f,
+                                                const index_t orient,
+                                                const EntityDof& entity,
+                                                ElemDof& element);
 
   // Get the degrees of freedom from the domain
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             index_t nz, class ElemDof, class EntityDof>
-  static void get_hex_domain_dof(const ElemDof& element, EntityDof& entity);
+  KOKKOS_FUNCTION static void get_hex_domain_dof(const ElemDof& element,
+                                                 EntityDof& entity);
 
   // Set the degrees of freedom into the element array
   template <index_t offset, bool ends, index_t ndof, index_t nx, index_t ny,
             index_t nz, class EntityDof, class ElemDof>
-  static void set_hex_domain_dof(const EntityDof& entity, ElemDof& element);
+  KOKKOS_FUNCTION static void set_hex_domain_dof(const EntityDof& entity,
+                                                 ElemDof& element);
 
   /**
    * @brief Wedge properties
@@ -469,26 +541,43 @@ class ElementTypes {
   static const index_t WEDGE_NVERTS = 6;
 
   // Given edge index, return edge vertex indices
-  static constexpr index_t WEDGE_EDGE_VERTS[][2] = {
-      {0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
+  KOKKOS_FUNCTION static const index_t (*get_wedge_edge_verts())[2] {
+    static constexpr index_t WEDGE_EDGE_VERTS[][2] = {
+        {0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
+    return WEDGE_EDGE_VERTS;
+  }
 
   // Given bound index, return edge indices
-  static constexpr index_t WEDGE_BOUND_NEDGES[] = {3, 3, 4, 4, 4};
-  static constexpr index_t WEDGE_BOUND_EDGES[][MAX_BOUND_EDGES] = {
-      {0, 1, 2, NO_INDEX},
-      {3, 4, 5, NO_INDEX},
-      {6, 3, 7, 0},
-      {7, 4, 8, 1},
-      {2, 8, 5, 6}};
+  KOKKOS_FUNCTION static const index_t* get_wedge_bound_nedges() {
+    static constexpr index_t WEDGE_BOUND_NEDGES[] = {3, 3, 4, 4, 4};
+    return WEDGE_BOUND_NEDGES;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_wedge_bound_edges())[MAX_BOUND_EDGES] {
+    static constexpr index_t WEDGE_BOUND_EDGES[][MAX_BOUND_EDGES] = {
+        {0, 1, 2, NO_INDEX},
+        {3, 4, 5, NO_INDEX},
+        {6, 3, 7, 0},
+        {7, 4, 8, 1},
+        {2, 8, 5, 6}};
+    return WEDGE_BOUND_EDGES;
+  }
 
   // Given bound index, return number of vertices/vertex indices
-  static constexpr index_t WEDGE_BOUND_NVERTS[] = {3, 3, 4, 4, 4};
-  static constexpr index_t WEDGE_BOUND_VERTS[][MAX_BOUND_VERTS] = {
-      {0, 1, 2, NO_INDEX},
-      {3, 4, 5, NO_INDEX},
-      {0, 3, 4, 1},
-      {1, 4, 5, 2},
-      {0, 2, 5, 3}};
+  KOKKOS_FUNCTION static const index_t* get_wedge_bound_nverts() {
+    static constexpr index_t WEDGE_BOUND_NVERTS[] = {3, 3, 4, 4, 4};
+    return WEDGE_BOUND_NVERTS;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_wedge_bound_verts())[MAX_BOUND_VERTS] {
+    static constexpr index_t WEDGE_BOUND_VERTS[][MAX_BOUND_VERTS] = {
+        {0, 1, 2, NO_INDEX},
+        {3, 4, 5, NO_INDEX},
+        {0, 3, 4, 1},
+        {1, 4, 5, 2},
+        {0, 2, 5, 3}};
+    return WEDGE_BOUND_VERTS;
+  }
 
   /**
    * @brief Pyramid properties
@@ -533,27 +622,44 @@ class ElementTypes {
   static const index_t PYRMD_NVERTS = 5;
 
   // Given edge index, return edge vertex indices
-  static constexpr index_t PYRMD_EDGE_VERTS[][2] = {
-      {0, 1}, {1, 2}, {2, 3}, {3, 0}, {0, 4}, {1, 4}, {2, 4}, {3, 4}};
+  KOKKOS_FUNCTION static const index_t (*get_pyrmd_edge_verts())[2] {
+    static constexpr index_t PYRMD_EDGE_VERTS[][2] = {
+        {0, 1}, {1, 2}, {2, 3}, {3, 0}, {0, 4}, {1, 4}, {2, 4}, {3, 4}};
+    return PYRMD_EDGE_VERTS;
+  }
 
   // Given bound index, return edge indices
-  static constexpr index_t PYRMD_BOUND_NEDGES[] = {3, 3, 3, 3, 4};
-  static constexpr index_t PYRMD_BOUND_EDGES[][MAX_BOUND_EDGES] = {
-      {0, 5, 4, NO_INDEX},
-      {1, 6, 5, NO_INDEX},
-      {2, 7, 6, NO_INDEX},
-      {4, 7, 3, NO_INDEX},
-      {3, 2, 1, 0}};
+  KOKKOS_FUNCTION static const index_t* get_pyrmd_bound_nedges() {
+    static constexpr index_t PYRMD_BOUND_NEDGES[] = {3, 3, 3, 3, 4};
+    return PYRMD_BOUND_NEDGES;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_pyrmd_bound_edges())[MAX_BOUND_EDGES] {
+    static constexpr index_t PYRMD_BOUND_EDGES[][MAX_BOUND_EDGES] = {
+        {0, 5, 4, NO_INDEX},
+        {1, 6, 5, NO_INDEX},
+        {2, 7, 6, NO_INDEX},
+        {4, 7, 3, NO_INDEX},
+        {3, 2, 1, 0}};
+    return PYRMD_BOUND_EDGES;
+  }
 
   // Given bound index, return number of vertices/vertex indices
-  static constexpr index_t PYRMD_BOUND_NVERTS[] = {3, 3, 3, 3, 4};
-  static constexpr index_t PYRMD_BOUND_VERTS[][MAX_BOUND_VERTS] = {
-      {0, 1, 4, NO_INDEX},
-      {1, 2, 4, NO_INDEX},
-      {2, 3, 4, NO_INDEX},
-      {0, 4, 3, NO_INDEX},
-      {0, 3, 2, 1}};
-};
+  KOKKOS_FUNCTION static const index_t* get_pyrmd_bound_nverts() {
+    static constexpr index_t PYRMD_BOUND_NVERTS[] = {3, 3, 3, 3, 4};
+    return PYRMD_BOUND_NVERTS;
+  }
+  KOKKOS_FUNCTION static const index_t (
+      *get_pyrmd_bound_verts())[MAX_BOUND_VERTS] {
+    static constexpr index_t PYRMD_BOUND_VERTS[][MAX_BOUND_VERTS] = {
+        {0, 1, 4, NO_INDEX},
+        {1, 2, 4, NO_INDEX},
+        {2, 3, 4, NO_INDEX},
+        {0, 4, 3, NO_INDEX},
+        {0, 3, 2, 1}};
+    return PYRMD_BOUND_VERTS;
+  }
+};  // namespace A2D
 
 }  // namespace A2D
 
diff --git a/include/multiphysics/feelementvector.h b/include/multiphysics/feelementvector.h
index c70e51d3..a7418eb4 100644
--- a/include/multiphysics/feelementvector.h
+++ b/include/multiphysics/feelementvector.h
@@ -137,7 +137,7 @@ class ElementVector_Empty {
  * @tparam VecType type of the solution vector, e.g. SolutionVector<...>
  */
 template <typename T, class Basis, class VecType>
-class ElementVector_Serial : public ElementVector_Empty {
+class ElementVector_Serial final : public ElementVector_Empty {
  public:
   static constexpr ElemVecType evtype = ElemVecType::Serial;
   ElementVector_Serial(ElementMesh<Basis>& mesh, VecType& vec) : mesh(mesh), vec(vec) {}
@@ -247,13 +247,19 @@ class ElementVector_Serial : public ElementVector_Empty {
  * @tparam VecType type of the solution vector, e.g. SolutionVector<...>
  */
 template <typename T, class Basis, class VecType>
-class ElementVector_Parallel : public ElementVector_Empty {
+class ElementVector_Parallel final : public ElementVector_Empty {
   using ElemVecArray_t = MultiArrayNew<T * [Basis::ndof]>;
 
  public:
   static constexpr ElemVecType evtype = ElemVecType::Parallel;
+
+  // Default and copy constructor
+  // Note: copy constructor is needed because we need to capture *this (i.e.
+  //  pass *this by value in the functor) in parallel dispatch
   ElementVector_Parallel(ElementMesh<Basis>& mesh, VecType& vec)
       : mesh(mesh), vec(vec), elem_vec_array("elem_vec_array", mesh.get_num_elements()) {}
+  KOKKOS_FUNCTION ElementVector_Parallel(const ElementVector_Parallel& other)
+      : mesh(other.mesh), vec(other.vec), elem_vec_array(other.elem_vec_array) {}
 
   class FEDof {
    public:
@@ -285,7 +291,7 @@ class ElementVector_Parallel : public ElementVector_Empty {
    * @brief Populate element-view data from global data for all elements
    */
   void get_values() const {
-    auto loop_body = KOKKOS_LAMBDA(const index_t elem) {
+    auto loop_body = KOKKOS_CLASS_LAMBDA(const index_t elem) {
       operate_element_values<ELEM_VALS_OP::GET, Basis::nbasis>(elem);
     };
 
@@ -302,7 +308,7 @@ class ElementVector_Parallel : public ElementVector_Empty {
    * Note: Data consistency is assumed
    */
   void set_values() const {
-    auto loop_body = KOKKOS_LAMBDA(const index_t elem) {
+    auto loop_body = KOKKOS_CLASS_LAMBDA(const index_t elem) {
       operate_element_values<ELEM_VALS_OP::SET, Basis::nbasis>(elem);
     };
 
@@ -318,7 +324,7 @@ class ElementVector_Parallel : public ElementVector_Empty {
    * @brief Add global data from element-view data for all elements
    */
   void add_values() const {
-    auto loop_body = KOKKOS_LAMBDA(const index_t elem) {
+    auto loop_body = KOKKOS_CLASS_LAMBDA(const index_t elem) {
       operate_element_values<ELEM_VALS_OP::ADD, Basis::nbasis>(elem);
     };
 
@@ -341,8 +347,8 @@ class ElementVector_Parallel : public ElementVector_Empty {
   template <ELEM_VALS_OP op, index_t nbasis>
   KOKKOS_FUNCTION void operate_element_values(const index_t& elem_idx) const {
     for (index_t i = 0; i < Basis::template get_ndof<nbasis - 1>(); i++) {
-      const int& sign = mesh.template get_global_dof_sign<nbasis - 1>(elem_idx, i);
-      const index_t& dof_index = mesh.template get_global_dof<nbasis - 1>(elem_idx, i);
+      const int sign = mesh.template get_global_dof_sign<nbasis - 1>(elem_idx, i);
+      const index_t dof_index = mesh.template get_global_dof<nbasis - 1>(elem_idx, i);
       const index_t dof_idx = i + Basis::template get_dof_offset<nbasis - 1>();
       if constexpr (op == ELEM_VALS_OP::GET) {
         elem_vec_array(elem_idx, dof_idx) = sign * vec[dof_index];
@@ -361,8 +367,8 @@ class ElementVector_Parallel : public ElementVector_Empty {
     return;
   }
 
-  ElementMesh<Basis>& mesh;
-  VecType& vec;
+  ElementMesh<Basis> mesh;
+  VecType vec;
   ElemVecArray_t elem_vec_array;  // The heavy-weight storage
 };
 
diff --git a/include/multiphysics/femesh-inl.h b/include/multiphysics/femesh-inl.h
index 0bd10af8..14c8218a 100644
--- a/include/multiphysics/femesh-inl.h
+++ b/include/multiphysics/femesh-inl.h
@@ -1211,9 +1211,10 @@ template <class Basis>
 ElementMesh<Basis>::ElementMesh(MeshConnectivityBase& conn)
     : nelems(conn.get_num_elements()), num_dof(0) {
   // Count up the number of degrees of freedom
-  element_dof = new index_t[nelems * ndof_per_element];
-  element_sign = new int[nelems * ndof_per_element];
-  std::fill(element_dof, element_dof + nelems * ndof_per_element, NO_INDEX);
+  num_dof_offset = DofOffsetArray("num_dof_offset");
+  element_dof = ElementDofArray("element_dof", nelems);
+  element_sign = ElementSignArray("element_sign", nelems);
+  BLAS::fill(element_dof, NO_INDEX);
 
   // Perform a sweep of the elements
   std::vector<index_t> ids(nelems, NO_INDEX), stack(nelems);
@@ -1264,8 +1265,8 @@ ElementMesh<Basis>::ElementMesh(MeshConnectivityBase& conn)
     for (index_t counter = nelems; counter > 0; counter--) {
       index_t elem = stack[counter - 1];
 
-      index_t* elem_dof = &element_dof[elem * ndof_per_element];
-      int* elem_sign = &element_sign[elem * ndof_per_element];
+      auto elem_dof = Kokkos::subview(element_dof, elem, Kokkos::ALL);
+      auto elem_sign = Kokkos::subview(element_sign, elem, Kokkos::ALL);
 
       // The domain DOF are always owned by the element - no need to check
       // for the element that owns them
@@ -1296,15 +1297,16 @@ ElementMesh<Basis>::ElementMesh(MeshConnectivityBase& conn)
           const index_t* owner_bounds;
           index_t nf_owner = conn.get_element_bounds(owner_elem, &owner_bounds);
 
+          auto owner_elem_dof =
+              Kokkos::subview(element_dof, owner_elem, Kokkos::ALL);
+
           for (index_t i = 0; i < nf_owner; i++) {
             if (owner_bounds[i] == bound) {
               index_t ref[ET::MAX_BOUND_VERTS], verts[ET::MAX_BOUND_VERTS];
               index_t nverts = conn.get_element_bound_verts(owner_elem, i, ref);
               conn.get_element_bound_verts(elem, index, verts);
 
-              Basis::get_entity_dof(basis, ET::Bound, i,
-                                    &element_dof[owner_elem * ndof_per_element],
-                                    dof);
+              Basis::get_entity_dof(basis, ET::Bound, i, owner_elem_dof, dof);
 
               if (nverts == 4) {
                 orient = ET::get_quad_domain_orientation(ref, verts);
@@ -1340,15 +1342,16 @@ ElementMesh<Basis>::ElementMesh(MeshConnectivityBase& conn)
           const index_t* owner_edges;
           index_t ne_owner = conn.get_element_edges(owner_elem, &owner_edges);
 
+          auto owner_elem_dof =
+              Kokkos::subview(element_dof, owner_elem, Kokkos::ALL);
+
           for (index_t i = 0; i < ne_owner; i++) {
             if (owner_edges[i] == edge) {
               index_t ref[2], verts[2];
               conn.get_element_edge_verts(owner_elem, i, ref);
               conn.get_element_edge_verts(elem, index, verts);
 
-              Basis::get_entity_dof(basis, ET::Edge, i,
-                                    &element_dof[owner_elem * ndof_per_element],
-                                    dof);
+              Basis::get_entity_dof(basis, ET::Edge, i, owner_elem_dof, dof);
 
               if (ref[0] == verts[1] && ref[1] == verts[0]) {
                 orient = 1;
@@ -1380,11 +1383,12 @@ ElementMesh<Basis>::ElementMesh(MeshConnectivityBase& conn)
           const index_t* owner_verts;
           index_t nv_owner = conn.get_element_verts(owner_elem, &owner_verts);
 
+          auto owner_elem_dof =
+              Kokkos::subview(element_dof, owner_elem, Kokkos::ALL);
+
           for (index_t i = 0; i < nv_owner; i++) {
             if (owner_verts[i] == vert) {
-              Basis::get_entity_dof(basis, ET::Vertex, i,
-                                    &element_dof[owner_elem * ndof_per_element],
-                                    dof);
+              Basis::get_entity_dof(basis, ET::Vertex, i, owner_elem_dof, dof);
               break;
             }
           }
@@ -1421,8 +1425,9 @@ template <class InteriorBasis>
 ElementMesh<Basis>::ElementMesh(const index_t label, MeshConnectivityBase& conn,
                                 ElementMesh<InteriorBasis>& mesh)
     : nelems(conn.get_num_boundary_bounds_with_label(label)) {
-  element_dof = new index_t[nelems * ndof_per_element];
-  element_sign = new int[nelems * ndof_per_element];
+  num_dof_offset = DofOffsetArray("num_dof_offset");
+  element_dof = ElementDofArray("element_dof", nelems);
+  element_sign = ElementSignArray("element_sign", nelems);
 
   // Get the number of boundary bounds
   const index_t* boundary_bounds;
@@ -1449,14 +1454,11 @@ ElementMesh<Basis>::ElementMesh(const index_t label, MeshConnectivityBase& conn,
       }
 
       // Get the signs associated wtih the original mesh
-      const index_t* dof;
-      const int* signs;
-      mesh.get_element_dof(elem, &dof);
-      mesh.get_element_signs(elem, &signs);
+      auto dof = mesh.get_element_dof(elem);
 
       // Set pointers for the entity dof
-      index_t* surf_dof = &element_dof[elem_count * ndof_per_element];
-      int* surf_signs = &element_sign[elem_count * ndof_per_element];
+      auto surf_dof = Kokkos::subview(element_dof, elem_count, Kokkos::ALL);
+      auto surf_signs = Kokkos::subview(element_sign, elem_count, Kokkos::ALL);
 
       // The degree of freedom indices for each entity
       index_t entity_dof[InteriorBasis::ndof];
@@ -1549,11 +1551,12 @@ ElementMesh<Basis>::ElementMesh(const index_t label, MeshConnectivityBase& conn,
  */
 template <class Basis>
 template <class HOrderBasis>
-ElementMesh<Basis>::ElementMesh(ElementMesh<HOrderBasis>& mesh)
+ElementMesh<Basis>::ElementMesh(const ElementMesh<HOrderBasis>& mesh)
     : nelems(HOrderBasis::get_num_lorder_elements() * mesh.get_num_elements()),
       num_dof(mesh.get_num_dof()) {
-  element_dof = new index_t[nelems * ndof_per_element];
-  element_sign = new int[nelems * ndof_per_element];
+  num_dof_offset = DofOffsetArray("num_dof_offset");
+  element_dof = ElementDofArray("element_dof", nelems);
+  element_sign = ElementSignArray("element_sign", nelems);
 
   for (index_t i = 0; i < Basis::nbasis; i++) {
     num_dof_offset[i] = mesh.get_num_cumulative_dof(i);
@@ -1563,20 +1566,18 @@ ElementMesh<Basis>::ElementMesh(ElementMesh<HOrderBasis>& mesh)
   // representation
   for (index_t j = 0; j < mesh.get_num_elements(); j++) {
     // Get the signs associated wtih the high-order mesh
-    const index_t* horder_dof;
-    const int* horder_signs;
-    mesh.get_element_dof(j, &horder_dof);
-    mesh.get_element_signs(j, &horder_signs);
+    auto horder_dof = mesh.get_element_dof(j);
+    auto horder_signs = mesh.get_element_signs(j);
 
     for (index_t i = 0; i < HOrderBasis::get_num_lorder_elements(); i++) {
       index_t elem = i + j * HOrderBasis::get_num_lorder_elements();
 
       // Index into the high-order element
-      index_t* lorder_dof = &element_dof[ndof_per_element * elem];
+      auto lorder_dof = Kokkos::subview(element_dof, elem, Kokkos::ALL);
+      auto lorder_signs = Kokkos::subview(element_sign, elem, Kokkos::ALL);
       HOrderBasis::get_lorder_dof(i, horder_dof, lorder_dof);
 
       // Signs indicating any orientation flip
-      int* lorder_signs = &element_sign[ndof_per_element * elem];
       HOrderBasis::get_lorder_signs(i, horder_signs, lorder_signs);
     }
   }
@@ -1584,16 +1585,14 @@ ElementMesh<Basis>::ElementMesh(ElementMesh<HOrderBasis>& mesh)
 
 template <class Basis>
 template <index_t basis>
-int ElementMesh<Basis>::get_global_dof_sign(index_t elem, index_t index) {
-  return element_sign[ndof_per_element * elem +
-                      Basis::template get_dof_offset<basis>() + index];
+int ElementMesh<Basis>::get_global_dof_sign(index_t elem, index_t index) const {
+  return element_sign(elem, Basis::template get_dof_offset<basis>() + index);
 }
 
 template <class Basis>
 template <index_t basis>
-index_t ElementMesh<Basis>::get_global_dof(index_t elem, index_t index) {
-  return element_dof[ndof_per_element * elem +
-                     Basis::template get_dof_offset<basis>() + index];
+index_t ElementMesh<Basis>::get_global_dof(index_t elem, index_t index) const {
+  return element_dof(elem, Basis::template get_dof_offset<basis>() + index);
 }
 
 template <class Basis>
@@ -1610,9 +1609,8 @@ void ElementMesh<Basis>::create_block_csr(index_t& nrows,
     index_t dof_reduced[Basis::ndof];
     index_t n = 0;
     for (index_t j1 = 0; j1 < ndof_per_elem; j1++, n++) {
-      index_t row = element_dof[i * Basis::ndof + j1] / M;
-      while (j1 + 1 < ndof_per_elem &&
-             row == (element_dof[i * Basis::ndof + j1 + 1] / M)) {
+      index_t row = element_dof(i, j1) / M;
+      while (j1 + 1 < ndof_per_elem && row == (element_dof(i, j1 + 1) / M)) {
         j1++;
       }
       dof_reduced[n] = row;
@@ -1701,8 +1699,7 @@ DirichletBCs<Basis>::DirichletBCs(MeshConnectivityBase& conn,
       conn.get_bound_elements(bound, &elem, &e2);
 
       // Get the degrees of freedom for this element
-      const index_t* elem_dof;
-      mesh.get_element_dof(elem, &elem_dof);
+      auto elem_dof = mesh.get_element_dof(elem);
 
       // Array of dof that may be extracted
       index_t dof[Basis::ndof];
diff --git a/include/multiphysics/femesh.h b/include/multiphysics/femesh.h
index 59c9f56a..ce2a63ae 100644
--- a/include/multiphysics/femesh.h
+++ b/include/multiphysics/femesh.h
@@ -112,7 +112,7 @@ class MetaDataFactory {
         char msg[256];
         std::snprintf(msg, 256,
                       "element enumerator %d does not represent a 2d element",
-                      element);
+                      (int)element);
         throw std::runtime_error(msg);
     }
   }
@@ -132,7 +132,7 @@ class MetaDataFactory {
         char msg[256];
         std::snprintf(msg, 256,
                       "element enumerator %d does not represent a 3d element",
-                      element);
+                      (int)element);
         throw std::runtime_error(msg);
     }
   }
@@ -140,98 +140,98 @@ class MetaDataFactory {
  private:
   static ElemConnMetaData create_tri(index_t* const* tri_bounds,
                                      const index_t* const* tri_verts) {
-    return ElemConnMetaData(ET::TRI_NBOUNDS,       // NBOUNDS
-                            0,                     // NEDGES
-                            ET::TRI_NVERTS,        // NVERTS
-                            ET::TRI_BOUND_NVERTS,  // BOUND_NVERTS
-                            ET::TRI_BOUND_VERTS,   // BOUND_VERTS
-                            nullptr,               // BOUND_NEDGES
-                            nullptr,               // BOUND_EDGES
-                            nullptr,               // EDGE_VERTS
-                            tri_bounds,            // bounds
-                            nullptr,               // edges
-                            tri_verts              // verts
+    return ElemConnMetaData(ET::TRI_NBOUNDS,             // NBOUNDS
+                            0,                           // NEDGES
+                            ET::TRI_NVERTS,              // NVERTS
+                            ET::get_tri_bound_nverts(),  // BOUND_NVERTS
+                            ET::get_tri_bound_verts(),   // BOUND_VERTS
+                            nullptr,                     // BOUND_NEDGES
+                            nullptr,                     // BOUND_EDGES
+                            nullptr,                     // EDGE_VERTS
+                            tri_bounds,                  // bounds
+                            nullptr,                     // edges
+                            tri_verts                    // verts
     );
   }
 
   static ElemConnMetaData create_quad(index_t* const* quad_bounds,
                                       const index_t* const* quad_verts) {
-    return ElemConnMetaData(ET::QUAD_NBOUNDS,       // NBOUNDS
-                            0,                      // NEDGES
-                            ET::QUAD_NVERTS,        // NVERTS
-                            ET::QUAD_BOUND_NVERTS,  // BOUND_NVERTS
-                            ET::QUAD_BOUND_VERTS,   // BOUND_VERTS
-                            nullptr,                // BOUND_NEDGES
-                            nullptr,                // BOUND_EDGES
-                            nullptr,                // EDGE_VERTS
-                            quad_bounds,            // bounds
-                            nullptr,                // edges
-                            quad_verts              // verts
+    return ElemConnMetaData(ET::QUAD_NBOUNDS,             // NBOUNDS
+                            0,                            // NEDGES
+                            ET::QUAD_NVERTS,              // NVERTS
+                            ET::get_quad_bound_nverts(),  // BOUND_NVERTS
+                            ET::get_quad_bound_verts(),   // BOUND_VERTS
+                            nullptr,                      // BOUND_NEDGES
+                            nullptr,                      // BOUND_EDGES
+                            nullptr,                      // EDGE_VERTS
+                            quad_bounds,                  // bounds
+                            nullptr,                      // edges
+                            quad_verts                    // verts
     );
   }
 
   static ElemConnMetaData create_tet(index_t* const* tet_bounds,
                                      index_t* const* tet_edges,
                                      const index_t* const* tet_verts) {
-    return ElemConnMetaData(ET::TET_NBOUNDS,       // NBOUNDS
-                            ET::TET_NEDGES,        // NEDGES
-                            ET::TET_NVERTS,        // NVERTS
-                            ET::TET_BOUND_NVERTS,  // BOUND_NVERTS
-                            ET::TET_BOUND_VERTS,   // BOUND_VERTS
-                            ET::TET_BOUND_NEDGES,  // BOUND_NEDGES
-                            ET::TET_BOUND_EDGES,   // BOUND_EDGES
-                            ET::TET_EDGE_VERTS,    // EDGE_VERTS
-                            tet_bounds,            // bounds
-                            tet_edges,             // edges
-                            tet_verts              // verts
+    return ElemConnMetaData(ET::TET_NBOUNDS,             // NBOUNDS
+                            ET::TET_NEDGES,              // NEDGES
+                            ET::TET_NVERTS,              // NVERTS
+                            ET::get_tet_bound_nverts(),  // BOUND_NVERTS
+                            ET::get_tet_bound_verts(),   // BOUND_VERTS
+                            ET::get_tet_bound_nedges(),  // BOUND_NEDGES
+                            ET::get_tet_bound_edges(),   // BOUND_EDGES
+                            ET::get_tet_edge_verts(),    // EDGE_VERTS
+                            tet_bounds,                  // bounds
+                            tet_edges,                   // edges
+                            tet_verts                    // verts
     );
   }
   static ElemConnMetaData create_hex(index_t* const* hex_bounds,
                                      index_t* const* hex_edges,
                                      const index_t* const* hex_verts) {
-    return ElemConnMetaData(ET::HEX_NBOUNDS,       // NBOUNDS
-                            ET::HEX_NEDGES,        // NEDGES
-                            ET::HEX_NVERTS,        // NVERTS
-                            ET::HEX_BOUND_NVERTS,  // BOUND_NVERTS
-                            ET::HEX_BOUND_VERTS,   // BOUND_VERTS
-                            ET::HEX_BOUND_NEDGES,  // BOUND_NEDGES
-                            ET::HEX_BOUND_EDGES,   // BOUND_EDGES
-                            ET::HEX_EDGE_VERTS,    // EDGE_VERTS
-                            hex_bounds,            // bounds
-                            hex_edges,             // edges
-                            hex_verts              // verts
+    return ElemConnMetaData(ET::HEX_NBOUNDS,             // NBOUNDS
+                            ET::HEX_NEDGES,              // NEDGES
+                            ET::HEX_NVERTS,              // NVERTS
+                            ET::get_hex_bound_nverts(),  // BOUND_NVERTS
+                            ET::get_hex_bound_verts(),   // BOUND_VERTS
+                            ET::get_hex_bound_nedges(),  // BOUND_NEDGES
+                            ET::get_hex_bound_edges(),   // BOUND_EDGES
+                            ET::get_hex_edge_verts(),    // EDGE_VERTS
+                            hex_bounds,                  // bounds
+                            hex_edges,                   // edges
+                            hex_verts                    // verts
     );
   }
   static ElemConnMetaData create_wedge(index_t* const* wedge_bounds,
                                        index_t* const* wedge_edges,
                                        const index_t* const* wedge_verts) {
-    return ElemConnMetaData(ET::WEDGE_NBOUNDS,       // NBOUNDS
-                            ET::WEDGE_NEDGES,        // NEDGES
-                            ET::WEDGE_NVERTS,        // NVERTS
-                            ET::WEDGE_BOUND_NVERTS,  // BOUND_NVERTS
-                            ET::WEDGE_BOUND_VERTS,   // BOUND_VERTS
-                            ET::WEDGE_BOUND_NEDGES,  // BOUND_NEDGES
-                            ET::WEDGE_BOUND_EDGES,   // BOUND_EDGES
-                            ET::WEDGE_EDGE_VERTS,    // EDGE_VERTS
-                            wedge_bounds,            // bounds
-                            wedge_edges,             // edges
-                            wedge_verts              // verts
+    return ElemConnMetaData(ET::WEDGE_NBOUNDS,             // NBOUNDS
+                            ET::WEDGE_NEDGES,              // NEDGES
+                            ET::WEDGE_NVERTS,              // NVERTS
+                            ET::get_wedge_bound_nverts(),  // BOUND_NVERTS
+                            ET::get_wedge_bound_verts(),   // BOUND_VERTS
+                            ET::get_wedge_bound_nedges(),  // BOUND_NEDGES
+                            ET::get_wedge_bound_edges(),   // BOUND_EDGES
+                            ET::get_wedge_edge_verts(),    // EDGE_VERTS
+                            wedge_bounds,                  // bounds
+                            wedge_edges,                   // edges
+                            wedge_verts                    // verts
     );
   }
   static ElemConnMetaData create_pyrmd(index_t* const* pyrmd_bounds,
                                        index_t* const* pyrmd_edges,
                                        const index_t* const* pyrmd_verts) {
-    return ElemConnMetaData(ET::PYRMD_NBOUNDS,       // NBOUNDS
-                            ET::PYRMD_NEDGES,        // NEDGES
-                            ET::PYRMD_NVERTS,        // NVERTS
-                            ET::PYRMD_BOUND_NVERTS,  // BOUND_NVERTS
-                            ET::PYRMD_BOUND_VERTS,   // BOUND_VERTS
-                            ET::PYRMD_BOUND_NEDGES,  // BOUND_NEDGES
-                            ET::PYRMD_BOUND_EDGES,   // BOUND_EDGES
-                            ET::PYRMD_EDGE_VERTS,    // EDGE_VERTS
-                            pyrmd_bounds,            // bounds
-                            pyrmd_edges,             // edges
-                            pyrmd_verts              // verts
+    return ElemConnMetaData(ET::PYRMD_NBOUNDS,             // NBOUNDS
+                            ET::PYRMD_NEDGES,              // NEDGES
+                            ET::PYRMD_NVERTS,              // NVERTS
+                            ET::get_pyrmd_bound_nverts(),  // BOUND_NVERTS
+                            ET::get_pyrmd_bound_verts(),   // BOUND_VERTS
+                            ET::get_pyrmd_bound_nedges(),  // BOUND_NEDGES
+                            ET::get_pyrmd_bound_edges(),   // BOUND_EDGES
+                            ET::get_pyrmd_edge_verts(),    // EDGE_VERTS
+                            pyrmd_bounds,                  // bounds
+                            pyrmd_edges,                   // edges
+                            pyrmd_verts                    // verts
     );
   }
 };
@@ -439,34 +439,51 @@ class ElementMesh {
   // Number of degrees of freedom for each element
   static const index_t ndof_per_element = Basis::ndof;
 
+  // Mult-dimensional array type
+  using DofOffsetArray = MultiArrayNew<index_t[Basis::nbasis]>;
+  using ElementDofArray = MultiArrayNew<index_t* [ndof_per_element]>;
+  using ElementSignArray = MultiArrayNew<int* [ndof_per_element]>;
+
   // Constructors
   ElementMesh(MeshConnectivityBase& conn);
   template <class InteriorBasis>
   ElementMesh(const index_t label, MeshConnectivityBase& conn,
               ElementMesh<InteriorBasis>& mesh);
   template <class HOrderBasis>
-  ElementMesh(ElementMesh<HOrderBasis>& mesh);
-
-  index_t get_num_elements() { return nelems; }
-  index_t get_num_dof() { return num_dof; }
-  index_t get_num_cumulative_dof(index_t basis) {
+  ElementMesh(const ElementMesh<HOrderBasis>& mesh);
+
+  // Copy constructor, this is needed for parallel dispatch
+  // Note, this is also a specialization of the third constructor above
+  template <>
+  KOKKOS_FUNCTION ElementMesh(const ElementMesh<Basis>& other)
+      : nelems(other.nelems),
+        num_dof(other.num_dof),
+        num_dof_offset(other.num_dof_offset),
+        element_dof(other.element_dof),
+        element_sign(other.element_sign) {}
+
+  index_t get_num_elements() const { return nelems; }
+  index_t get_num_dof() const { return num_dof; }
+  index_t get_num_cumulative_dof(index_t basis) const {
     return num_dof_offset[basis];
   }
 
+  // Needed for parallel element execution
   template <index_t basis>
-  int get_global_dof_sign(index_t elem, index_t index);
+  KOKKOS_FUNCTION int get_global_dof_sign(const index_t elem,
+                                          const index_t index) const;
   template <index_t basis>
-  index_t get_global_dof(index_t elem, index_t index);
+  KOKKOS_FUNCTION index_t get_global_dof(const index_t elem,
+                                         const index_t index) const;
 
   // Get the degrees of freedom associated with this element
-  KOKKOS_FUNCTION void get_element_dof(const index_t elem,
-                                       const index_t* dof[]) {
-    *dof = &element_dof[ndof_per_element * elem];
+  KOKKOS_FUNCTION auto get_element_dof(const index_t elem) {
+    return Kokkos::subview(element_dof, elem, Kokkos::ALL);
   }
 
   // Get the signs associated with the degrees of freedom
-  void get_element_signs(const index_t elem, const int* signs[]) {
-    *signs = &element_sign[ndof_per_element * elem];
+  auto get_element_signs(const index_t elem) {
+    return Kokkos::subview(element_sign, elem, Kokkos::ALL);
   }
 
   template <index_t M, index_t basis_offset = Basis::nbasis>
@@ -474,14 +491,14 @@ class ElementMesh {
                         std::vector<index_t>& cols);
 
  private:
-  index_t nelems;                         // Total number of elements
-  index_t num_dof;                        // Total number of degrees of freedom
-  index_t num_dof_offset[Basis::nbasis];  // Cumulative number of degrees of
-                                          // freedom each basis
+  index_t nelems;                 // Total number of elements
+  index_t num_dof;                // Total number of degrees of freedom
+  DofOffsetArray num_dof_offset;  // Cumulative number of degrees of
+                                  // freedom each basis
 
   // Store the degrees of freedom for each element and the element sign
-  index_t* element_dof;
-  int* element_sign;
+  ElementDofArray element_dof;
+  ElementSignArray element_sign;
 };
 
 template <class Basis>
diff --git a/include/multiphysics/hex_tools.h b/include/multiphysics/hex_tools.h
index d9b71f50..e0edf521 100644
--- a/include/multiphysics/hex_tools.h
+++ b/include/multiphysics/hex_tools.h
@@ -12,8 +12,6 @@ namespace A2D {
 template <class GeoBasis, typename I, typename T, class GeoElemVec>
 void set_geo_from_hex_nodes(const index_t nhex, const I hex[], const T Xloc[],
                             GeoElemVec &elem_geo) {
-  constexpr ElemVecType evtype = GeoElemVec::evtype;
-
   for (int e = 0; e < nhex; e++) {
     // Get the geometry values
     typename GeoElemVec::FEDof geo_dof(e, elem_geo);
@@ -54,20 +52,14 @@ void set_geo_from_hex_nodes(const index_t nhex, const I hex[], const T Xloc[],
       }
     }
 
-    if constexpr (evtype == ElemVecType::Serial) {
-      elem_geo.set_element_values(e, geo_dof);
-    }
-  }
-  if constexpr (evtype == ElemVecType::Parallel) {
-    elem_geo.set_values();
+    elem_geo.set_element_values(e, geo_dof);
   }
+  elem_geo.set_values();
 }
 
 template <class GeoBasis, typename I, typename T, class GeoElemVec>
 void set_geo_from_quad_nodes(const index_t nquad, const I quad[],
                              const T Xloc[], GeoElemVec &elem_geo) {
-  constexpr ElemVecType evtype = GeoElemVec::evtype;
-
   for (int e = 0; e < nquad; e++) {
     // Get the geometry values
     typename GeoElemVec::FEDof geo_dof(e, elem_geo);
@@ -98,13 +90,9 @@ void set_geo_from_quad_nodes(const index_t nquad, const I quad[],
       }
     }
 
-    // if constexpr (evtype == ElemVecType::Serial) {
     elem_geo.set_element_values(e, geo_dof);
-    // }
   }
-  // if constexpr (evtype == ElemVecType::Parallel) {
   elem_geo.set_values();
-  // }
 }
 
 template <index_t outputs, index_t degree, typename T, class DataBasis,
@@ -115,6 +103,7 @@ void write_hex_to_vtk(DataElemVec &elem_data, GeoElemVec &elem_geo,
                       const FunctorType &func) {
   Timer timer("write_hex_to_vtk()");
   using ET = ElementTypes;
+  auto HEX_VERTS_CART = ET::get_hex_verts_cart();
   const index_t nex = degree;
   using QuadPts = HexGaussLobattoQuadrature<nex + 1>;
 
@@ -187,9 +176,9 @@ void write_hex_to_vtk(DataElemVec &elem_data, GeoElemVec &elem_geo,
 
           for (index_t ii = 0; ii < ET::HEX_NVERTS; ii++) {
             vtk_conn(counter, ii) =
-                off + vtk_node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                   j + ET::HEX_VERTS_CART[ii][1],
-                                   k + ET::HEX_VERTS_CART[ii][2]);
+                off + vtk_node_num(i + HEX_VERTS_CART[ii][0],
+                                   j + HEX_VERTS_CART[ii][1],
+                                   k + HEX_VERTS_CART[ii][2]);
           }
         }
       }
@@ -278,14 +267,15 @@ void write_quad_to_vtk(DataElemVec &elem_data, GeoElemVec &elem_geo,
       }
     }
 
+    auto QUAD_VERTS_CART = ET::get_quad_verts_cart();
     for (index_t j = 0; j < nex; j++) {
       for (index_t i = 0; i < nex; i++, counter++) {
         index_t off = n * (nex + 1) * (nex + 1);
 
         for (index_t ii = 0; ii < ET::QUAD_NVERTS; ii++) {
           vtk_conn(counter, ii) =
-              off + vtk_node_num(i + ET::QUAD_VERTS_CART[ii][0],
-                                 j + ET::QUAD_VERTS_CART[ii][1]);
+              off + vtk_node_num(i + QUAD_VERTS_CART[ii][0],
+                                 j + QUAD_VERTS_CART[ii][1]);
         }
       }
     }
diff --git a/include/multiphysics/lagrange_hypercube_basis.h b/include/multiphysics/lagrange_hypercube_basis.h
index a9291153..30bbcca7 100644
--- a/include/multiphysics/lagrange_hypercube_basis.h
+++ b/include/multiphysics/lagrange_hypercube_basis.h
@@ -224,10 +224,10 @@ class LagrangeH1HypercubeBasis {
    * @param orient Orientation of the entity relative to the reference
    * @param signs Array of sign values
    */
-  template <index_t offset>
+  template <index_t offset, class ElemSign>
   KOKKOS_FUNCTION static void set_entity_signs(ET::ElementEntity entity,
                                                index_t index, index_t orient,
-                                               int signs[]) {
+                                               ElemSign& signs) {
     int sgns[ndof];
     const index_t entity_ndof = get_entity_ndof(entity, index);
     for (index_t i = 0; i < entity_ndof; i++) {
@@ -324,10 +324,11 @@ class LagrangeH1HypercubeBasis {
    * @param horder_signs The signs for the high-order degrees of freedom
    * @param signs The signs relative to the high-order element
    */
-  template <index_t horder_offset, index_t lorder_offset>
+  template <index_t horder_offset, index_t lorder_offset, class HOrderSign,
+            class LOrderSign>
   KOKKOS_FUNCTION static void get_lorder_signs(const index_t n,
-                                               const int horder_signs[],
-                                               int signs[]) {
+                                               const HOrderSign& horder_signs,
+                                               LOrderSign& signs) {
     for (index_t p = 0; p < indexpow<2, dim>::value * C; p++) {
       signs[lorder_offset + p] = 1;
     }
@@ -685,14 +686,18 @@ class LagrangeH1HypercubeBasis {
         T u1[u1size];
         T u1x[u1size];
         T u1y[u1size];  // not used for dim == 2
-        std::fill(u1, u1 + u1size, T(0.0));
-        std::fill(u1x, u1x + u1size, T(0.0));
-        std::fill(u1y, u1y + u1size, T(0.0));
+        for (index_t k = 0; k < u1size; k++) {
+          u1[k] = T(0.0);
+          u1x[k] = T(0.0);
+          u1y[k] = T(0.0);
+        }
 
         T u0[u0size];
         T u0x[u0size];
-        std::fill(u0, u0 + u0size, T(0.0));
-        std::fill(u0x, u0x + u0size, T(0.0));
+        for (index_t k = 0; k < u0size; k++) {
+          u0[k] = T(0.0);
+          u0x[k] = T(0.0);
+        }
 
         const index_t qouter = dim == 2 ? 1 : q2dim;
         for (index_t q2 = 0; q2 < qouter; q2++) {
@@ -1051,10 +1056,10 @@ class LagrangeL2HypercubeBasis {
    * @param orient Orientation of the entity relative to the reference
    * @param signs Array of sign values
    */
-  template <index_t offset>
+  template <index_t offset, class ElemSign>
   KOKKOS_FUNCTION static void set_entity_signs(ET::ElementEntity entity,
                                                index_t index, index_t orient,
-                                               int signs[]) {
+                                               ElemSign& signs) {
     int sgns[ndof];
     const index_t entity_ndof = get_entity_ndof(entity, index);
     for (index_t i = 0; i < entity_ndof; i++) {
@@ -1124,10 +1129,11 @@ class LagrangeL2HypercubeBasis {
    * @param horder_signs The signs for the high-order degrees of freedom
    * @param signs The signs relative to the high-order element
    */
-  template <index_t horder_offset, index_t lorder_offset>
+  template <index_t horder_offset, index_t lorder_offset, class HOrderSign,
+            class LOrderSign>
   KOKKOS_FUNCTION static void get_lorder_signs(const index_t n,
-                                               const int horder_signs[],
-                                               int signs[]) {
+                                               const HOrderSign& horder_signs,
+                                               LOrderSign& signs) {
     for (index_t p = 0; p < C; p++) {
       signs[lorder_offset + p] = 1;
     }
@@ -1418,9 +1424,14 @@ class LagrangeL2HypercubeBasis {
         }
 
         T u0[u0size];
-        std::fill(u0, u0 + u0size, T(0.0));
+        for (index_t k = 0; k < u0size; k++) {
+          u0[k] = T(0.0);
+        }
+
         T u1[u1size];
-        std::fill(u1, u1 + u1size, T(0.0));
+        for (index_t k = 0; k < u1size; k++) {
+          u1[k] = T(0.0);
+        }
 
         const index_t qouter = dim == 2 ? 1 : q2dim;
         for (index_t q2 = 0; q2 < qouter; q2++) {
diff --git a/include/multiphysics/qhdiv_hex_basis.h b/include/multiphysics/qhdiv_hex_basis.h
index 80b96b02..dcb5867a 100644
--- a/include/multiphysics/qhdiv_hex_basis.h
+++ b/include/multiphysics/qhdiv_hex_basis.h
@@ -194,9 +194,9 @@ class QHdivHexBasis {
    * @param orient Orientation of the entity relative to the reference
    * @param signs Array of sign values
    */
-  template <index_t offset>
+  template <index_t offset, class ElemSign>
   static void set_entity_signs(ET::ElementEntity entity, index_t index,
-                               index_t orient, int signs[]) {
+                               index_t orient, ElemSign& signs) {
     int sgns[ndof];
     const index_t entity_ndof = get_entity_ndof(entity, index);
     for (index_t k = 0; k < entity_ndof; k++) {
@@ -273,9 +273,10 @@ class QHdivHexBasis {
    * @param horder_signs The signs for the high-order degrees of freedom
    * @param signs The signs relative to the high-order element
    */
-  template <index_t horder_offset, index_t lorder_offset>
-  static void get_lorder_signs(const index_t n, const int horder_signs[],
-                               int signs[]) {
+  template <index_t horder_offset, index_t lorder_offset, class HOrderSign,
+            class LOrderSign>
+  static void get_lorder_signs(const index_t n, const HOrderSign& horder_signs,
+                               LOrderSign& signs) {
     const index_t i = n % degree;
     const index_t j = (n % (degree * degree)) / degree;
     const index_t k = n / (degree * degree);
diff --git a/include/utils/a2dmesh.h b/include/utils/a2dmesh.h
index 9ccd887d..4762de27 100644
--- a/include/utils/a2dmesh.h
+++ b/include/utils/a2dmesh.h
@@ -10,6 +10,7 @@
 #define A2D_MESH_H
 
 #include <cstdlib>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -94,9 +95,9 @@ class MesherRect2D {
     T dx = lx / nx;
     T dy = ly / ny;
 
-    if (randomize) {
-      std::srand(random_seed);
-    }
+    // This is a portable rng engine
+    std::mt19937 engine(random_seed);
+    std::uniform_real_distribution<T> rng(-1.0, 1.0);
 
     // Set X
     for (I j = 0; j < ny + 1; j++) {
@@ -106,8 +107,8 @@ class MesherRect2D {
 
         if (randomize) {
           if (i and j and (i - nx) and (j - ny)) {
-            px = (T(std::rand()) / RAND_MAX - 0.5) * fraction * dx;
-            py = (T(std::rand()) / RAND_MAX - 0.5) * fraction * dy;
+            px = rng(engine) * fraction * dx;
+            py = rng(engine) * fraction * dy;
           }
         }
 
@@ -118,11 +119,12 @@ class MesherRect2D {
 
     // Set connectivity
     using ET = ElementTypes;
+    auto QUAD_VERTS_CART = ET::get_quad_verts_cart();
     for (I j = 0, e = 0; j < ny; j++) {
       for (I i = 0; i < nx; i++, e++) {
         for (I ii = 0; ii < ET::QUAD_NVERTS; ii++) {
-          quad[4 * e + ii] = node_num(i + ET::QUAD_VERTS_CART[ii][0],
-                                      j + ET::QUAD_VERTS_CART[ii][1]);
+          quad[4 * e + ii] =
+              node_num(i + QUAD_VERTS_CART[ii][0], j + QUAD_VERTS_CART[ii][1]);
         }
       }
     }
@@ -299,13 +301,14 @@ class MesherBrick3D {
 
     // Set connectivity
     using ET = ElementTypes;
+    auto HEX_VERTS_CART = ET::get_hex_verts_cart();
     for (I k = 0, e = 0; k < nz; k++) {
       for (I j = 0; j < ny; j++) {
         for (I i = 0; i < nx; i++, e++) {
           for (I ii = 0; ii < ET::HEX_NVERTS; ii++) {
-            hex[8 * e + ii] = node_num(i + ET::HEX_VERTS_CART[ii][0],
-                                       j + ET::HEX_VERTS_CART[ii][1],
-                                       k + ET::HEX_VERTS_CART[ii][2]);
+            hex[8 * e + ii] =
+                node_num(i + HEX_VERTS_CART[ii][0], j + HEX_VERTS_CART[ii][1],
+                         k + HEX_VERTS_CART[ii][2]);
           }
         }
       }