diff --git a/CMakeLists.txt b/CMakeLists.txt
index 99ff8b35..1d65cd14 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -172,6 +172,8 @@ meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
+meshfields_add_exe(IntegratorTest test/testIntegrator.cpp)
+meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -190,6 +192,9 @@ test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
 test_func(OmegahTetTest, ./OmegahTetTest)
+test_func(IntegratorTest ./IntegratorTest)
+test_func(BoxIntegratorTest ./BoxIntegratorTest)
+
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)
diff --git a/docs/Ordering.pdf b/docs/Ordering.pdf
new file mode 100644
index 00000000..941cf6a7
Binary files /dev/null and b/docs/Ordering.pdf differ
diff --git a/docs/Ordering.tex b/docs/Ordering.tex
new file mode 100644
index 00000000..bf6002b9
--- /dev/null
+++ b/docs/Ordering.tex
@@ -0,0 +1,76 @@
+\documentclass{article}
+\usepackage{graphicx}
+\usepackage[letterpaper,textwidth=5in,right=2.5in,textheight=9in]{geometry}
+
+\usepackage{amsmath,amssymb, changepage}
+\usepackage{enumerate}
+\usepackage{upgreek}
+\pagestyle{empty}
+\usepackage{listings}
+\usepackage{tikz}
+\begin{document}
+\noindent
+Joshua Kloepfer\\
+9/30/2025\\
+\textbf{Meshfields DOF Holder and Integration Point Ordering}\\\\
+Linear Triangle:\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0) -- (2.5, 2.5);
+    \draw[gray, thick] (0, 0) -- (5, 0);
+    \draw[gray, thick] (2.5, 2.5) -- (5, 0);
+    \filldraw[black] (0, 0) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (2.5, 2.5) circle (2pt) node[anchor=south]{V2};
+\end{tikzpicture}
+\\
+Quadratic Triangle\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0) -- (2.5, 2.5);
+    \draw[gray, thick] (0, 0) -- (5, 0);
+    \draw[gray, thick] (2.5, 2.5) -- (5, 0);
+    \filldraw[black] (0, 0) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (2.5, 2.5) circle (2pt) node[anchor=south]{V2};
+    \filldraw[black] (2.5, 0) circle (2pt) node[anchor=north]{E0};
+    \filldraw[black] (3.75, 1.25) circle (2pt) node[anchor=west]{E1};
+    \filldraw[black] (1.25, 1.25) circle (2pt) node[anchor=east]{E2};
+\end{tikzpicture}
+\\
+Linear Tetrahedron:\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0, 0) -- (5, 0, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (5, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (5, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (0, 5, 0) -- (0, 0, 5);
+    \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=east]{V2};
+    \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V3};
+\end{tikzpicture}
+\\
+Vertex ordering must follow right hand rule.\\
+Quadratic Tetrahedron:\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0, 0) -- (5, 0, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (5, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (5, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (0, 5, 0) -- (0, 0, 5);
+    \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=east]{V2};
+    \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V3};
+    
+    \filldraw[black] (2.5, 0, 2.5) circle (2pt) node[anchor=north]{E0};
+    \filldraw[black] (0, 0, 2.5) circle (2pt) node[anchor=east]{E1};
+    \filldraw[black] (0, 2.5, 2.5) circle (2pt) node[anchor=east]{E2};
+    \filldraw[black] (2.5, 0, 0) circle (2pt) node[anchor=south]{E3};
+    \filldraw[black] (2.5, 2.5, 0) circle (2pt) node[anchor=east]{E4};
+    \filldraw[black] (0, 2.5, 0) circle (2pt) node[anchor=west]{E5};
+\end{tikzpicture}
+\\
+Vertex ordering must follow right hand rule and edge ordering must follow diagram.
+\end{document}
diff --git a/mfemTest/IntegratorMfemLargeTest.cpp b/mfemTest/IntegratorMfemLargeTest.cpp
new file mode 100644
index 00000000..a55a6922
--- /dev/null
+++ b/mfemTest/IntegratorMfemLargeTest.cpp
@@ -0,0 +1,27 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double intFunc(const Vector &x) {
+  return x[0];
+}
+
+int main(int argc, char *argv[]) {
+  Mesh mesh = Mesh::MakeCartesian3D(100, 100, 100, Element::TETRAHEDRON, 100, 100, 100, false);
+  //6 million element mesh
+  H1_FECollection fec(1, mesh.Dimension());
+  FiniteElementSpace fespace(&mesh, &fec);
+  
+  FunctionCoefficient custom(intFunc);
+
+  LinearForm b(&fespace);
+
+  b.AddDomainIntegrator(new DomainLFIntegrator(custom)); 
+
+  b.Assemble(); 
+  cout << b.Sum() << endl;
+  return 0;
+}
diff --git a/mfemTest/IntegratorMfemTest.cpp b/mfemTest/IntegratorMfemTest.cpp
new file mode 100644
index 00000000..efb7c855
--- /dev/null
+++ b/mfemTest/IntegratorMfemTest.cpp
@@ -0,0 +1,37 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double intFunc(const Vector &x) {
+  return x[0]; 
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../data/star.mesh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+   Mesh mesh(mesh_file);
+
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+  
+   FunctionCoefficient custom(intFunc);
+
+   LinearForm b(&fespace);
+
+   b.AddDomainIntegrator(new DomainLFIntegrator(custom)); 
+
+   b.Assemble(); 
+
+   cout << b.Sum() << endl;
+   return 0;
+}
diff --git a/mfemTest/README.md b/mfemTest/README.md
new file mode 100644
index 00000000..07d2e9d7
--- /dev/null
+++ b/mfemTest/README.md
@@ -0,0 +1,7 @@
+Run using g++ IntegratorMfemTest.cpp  -I<path-to-mfem-build> -L<path-to-mfem-build>  -l mfem
+in order to run, use -m to specify the mesh file and -o to specify order
+mfem version: https://github.com/mfem/mfem/tree/2caa75e35a54df93d19f23655170254556dfc081
+For the build:
+mkdir <mfem-build-dir> ; cd <mfem-build-dir>
+cmake <mfem-source-dir>
+make -j 4
diff --git a/mfemTest/box.mesh b/mfemTest/box.mesh
new file mode 100644
index 00000000..5fc3abcb
--- /dev/null
+++ b/mfemTest/box.mesh
@@ -0,0 +1,53 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+6
+1 4 7 0 3 1
+1 4 7 0 1 5
+1 4 7 0 5 4
+1 4 7 0 2 3
+1 4 7 0 6 2
+1 4 7 0 4 6
+
+boundary
+12
+1 2 3 0 2
+1 2 0 3 1
+6 2 7 4 5
+6 2 4 7 6
+5 2 6 0 4
+5 2 0 6 2
+3 2 7 1 3
+3 2 1 7 5
+2 2 5 0 1
+2 2 0 5 4
+4 2 7 2 6
+4 2 2 7 3
+
+vertices
+8
+3
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+0 0 1
+1 0 1
+0 1 1
+1 1 1
diff --git a/mfemTest/square.mesh b/mfemTest/square.mesh
new file mode 100644
index 00000000..c59e1b80
--- /dev/null
+++ b/mfemTest/square.mesh
@@ -0,0 +1,36 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+2
+1 2 0 1 3
+1 2 2 3 1
+
+boundary
+4
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 0
+
+vertices
+4
+2
+0 0
+1 0
+1 1
+0 1
diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp
index 2197a0a6..34d4756c 100644
--- a/src/MeshField_Element.hpp
+++ b/src/MeshField_Element.hpp
@@ -303,6 +303,30 @@ struct FieldElement {
     return Kokkos::View<Real *>("foo", J.extent(0));
   }
 
+  template <size_t MeshEntDim>
+  KOKKOS_INLINE_FUNCTION auto getGradients(Kokkos::View<Real **> lc,
+                                           size_t pt) const {
+    if constexpr (ShapeType::Order == 1) {
+      return shapeFn.getLocalGradients();
+    } else {
+      Kokkos::Array<Real, MeshEntDim + 1> coord;
+      for (int i = 0; i < MeshEntDim + 1; ++i) {
+        coord[i] = lc(pt, i);
+      }
+      return shapeFn.getLocalGradients(coord);
+    }
+  }
+
+  KOKKOS_INLINE_FUNCTION auto setNodalGradients(
+      Kokkos::View<Real * [ShapeType::numNodes][MeshEntDim]> nodalGradients,
+      auto grad, size_t pt, size_t node, size_t d) const {
+    if constexpr (ShapeType::Order == 1) {
+      nodalGradients(pt, node, d) = grad[node * MeshEntDim + d];
+    } else {
+      nodalGradients(pt, node, d) = grad[node][d];
+    }
+  }
+
   /**
    * @brief
    * Given an array of parametric coordinates 'localCoords', one per mesh
@@ -380,22 +404,22 @@ struct FieldElement {
       // one matrix per point
       Kokkos::View<Real ***> res("result", numPts, MeshEntDim, MeshEntDim);
       Kokkos::deep_copy(res, 0.0); // initialize all entries to zero
-
       // fill the views of node coordinates and node gradients
       Kokkos::View<Real * [ShapeType::numNodes][MeshEntDim]> nodeCoords(
           "nodeCoords", numPts);
       Kokkos::View<Real * [ShapeType::numNodes][MeshEntDim]> nodalGradients(
           "nodalGradients", numPts);
-      const auto grad = shapeFn.getLocalGradients();
       Kokkos::parallel_for(
           numMeshEnts, KOKKOS_CLASS_LAMBDA(const int ent) {
             const auto vals = getNodeValues(ent);
             assert(vals.size() == MeshEntDim * ShapeType::numNodes);
             for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) {
+              auto ptCoords = Kokkos::subview(localCoords, pt, Kokkos::ALL());
+              const auto grad = getGradients<MeshEntDim>(localCoords, pt);
               for (size_t node = 0; node < ShapeType::numNodes; node++) {
                 for (size_t d = 0; d < MeshEntDim; d++) {
                   nodeCoords(pt, node, d) = vals[node * MeshEntDim + d];
-                  nodalGradients(pt, node, d) = grad[node * MeshEntDim + d];
+                  setNodalGradients(nodalGradients, grad, pt, node, d);
                 }
               }
             }
diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp
index 0c16d68c..849638b7 100644
--- a/src/MeshField_Integrate.hpp
+++ b/src/MeshField_Integrate.hpp
@@ -92,16 +92,16 @@ class TetrahedronIntegration : public EntityIntegration<4> {
     virtual std::vector<IntegrationPoint<4>> getPoints() const {
 
       return {IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
-                                       0.138196601125011, 0.585410196624969},
+                                       0.138196601125011, 0.585410196624967},
                                0.25 / 6.0),
-              IntegrationPoint(Vector4{0.585410196624969, 0.138196601125011,
+              IntegrationPoint(Vector4{0.585410196624967, 0.138196601125011,
                                        0.138196601125011, 0.138196601125011},
                                0.25 / 6.0),
-              IntegrationPoint(Vector4{0.138196601125011, 0.585410196624969,
+              IntegrationPoint(Vector4{0.138196601125011, 0.585410196624967,
                                        0.138196601125011, 0.138196601125011},
                                0.25 / 6.0),
               IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
-                                       0.585410196624969, 0.138196601125011},
+                                       0.585410196624967, 0.138196601125011},
                                0.25 / 6.0)};
     }
     virtual int getAccuracy() const { return 2; }
diff --git a/test/testBoxIntegrator.cpp b/test/testBoxIntegrator.cpp
new file mode 100644
index 00000000..fa05b4ac
--- /dev/null
+++ b/test/testBoxIntegrator.cpp
@@ -0,0 +1,119 @@
+#include "KokkosController.hpp"
+#include "MeshField.hpp"
+#include "MeshField_Element.hpp"    //remove?
+#include "MeshField_Fail.hpp"       //remove?
+#include "MeshField_For.hpp"        //remove?
+#include "MeshField_ShapeField.hpp" //remove?
+#include "Omega_h_build.hpp"
+#include "Omega_h_file.hpp"
+#include "Omega_h_simplex.hpp"
+#include <Kokkos_Core.hpp>
+#include <MeshField_Integrate.hpp>
+#include <iostream>
+#include <sstream>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
+
+template <typename FieldElement>
+class testIntegrator : public MeshField::Integrator {
+private:
+  testIntegrator(){};
+
+protected:
+  MeshField::Real integral;
+  FieldElement &fes;
+
+public:
+  MeshField::Real getIntegral() { return integral; }
+  testIntegrator(FieldElement &fes_in, int order)
+      : Integrator(order), integral(0), fes(fes_in){};
+  void atPoints(Kokkos::View<MeshField::Real **> p,
+                Kokkos::View<MeshField::Real *> w,
+                Kokkos::View<MeshField::Real *> dV) {
+    const auto globalCoords = MeshField::evaluate(fes, p, p.extent(0) / fes.numMeshEnts);
+    Kokkos::parallel_reduce(
+        "integrate", p.extent(0),
+        KOKKOS_LAMBDA(const int &ent, MeshField::Real &integ) {
+          const auto x = globalCoords(ent, 0);
+          const auto y = globalCoords(ent, 1);
+          const auto z = globalCoords.extent(1) == 3 ? globalCoords(ent, 2) : 1;
+          integ += w(ent) * x * dV(ent);
+        },
+        integral);
+  }
+};
+
+template <template <typename...> typename Controller, size_t ShapeOrder,
+          size_t dim>
+void doRun(Omega_h::Mesh &mesh,
+           MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
+  auto field =
+      omf.template CreateLagrangeField<MeshField::Real, ShapeOrder, dim>();
+  auto coords = mesh.coords();
+  Kokkos::parallel_for(
+      mesh.nverts(), KOKKOS_LAMBDA(int vtx) {
+        field(vtx, 0, 0, MeshField::Vertex) = coords[vtx * dim];
+        field(vtx, 0, 1, MeshField::Vertex) = coords[vtx * dim + 1];
+        if constexpr (dim == 3) {
+          field(vtx, 0, 2, MeshField::Vertex) = coords[vtx * dim + 2];
+        }
+      });
+  if (ShapeOrder == 2) {
+    auto edge2vtx = mesh.ask_down(1, 0).ab2b;
+    auto edgeMap = mesh.ask_down(dim, 1).ab2b;
+    Kokkos::parallel_for(
+        mesh.nedges(), KOKKOS_LAMBDA(int edge) {
+          const auto left = edge2vtx[edge * 2];
+          const auto right = edge2vtx[edge * 2 + 1];
+          const auto x = (coords[left * dim] + coords[right * dim]) / 2.0;
+          const auto y =
+              (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0;
+          field(edge, 0, 0, MeshField::Edge) = x;
+          field(edge, 0, 1, MeshField::Edge) = y;
+          if constexpr (dim == 3) {
+            const auto z =
+                (coords[left * dim + 2] + coords[right * dim + 2]) / 2.0;
+            field(edge, 0, 2, MeshField::Edge) = z;
+          }
+        });
+  }
+
+  auto shapeSet = [&]() {
+    if constexpr (dim == 3) {
+      return MeshField::Omegah::getTetrahedronElement<ShapeOrder>(mesh);
+    } else {
+      return MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
+    }
+  };
+  const auto [shp, map] = shapeSet();
+  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
+  testIntegrator testInt(fes, ShapeOrder);
+  testInt.process(fes);
+  std::cout << testInt.getIntegral() << std::endl;
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  auto lib = Omega_h::Library(&argc, &argv);
+  auto world = lib.world();
+  const auto family = OMEGA_H_SIMPLEX;
+  auto len = 1.0;
+#ifdef MESHFIELDS_ENABLE_CABANA
+  {
+    Omega_h::Mesh mesh3D = Omega_h::build_box(world, family, 1.0, 1.0, 1.0, len, len, len);
+    MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::CabanaController>
+        omf3D(mesh3D);
+    Omega_h::Mesh mesh2D = Omega_h::build_box(world, family, 1.0, 1.0, 0.0, len, len, 0.0);
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::CabanaController> omf2D(mesh2D);
+    doRun<MeshField::CabanaController, 1>(mesh2D, omf2D);
+    doRun<MeshField::CabanaController, 2>(mesh2D, omf2D);
+    doRun<MeshField::CabanaController, 1>(mesh3D, omf3D);
+    doRun<MeshField::CabanaController, 2>(mesh3D, omf3D);
+  }
+#endif
+  {
+  }
+  Kokkos::finalize();
+  return 0;
+}
diff --git a/test/testIntegrator.cpp b/test/testIntegrator.cpp
new file mode 100644
index 00000000..cd7b4916
--- /dev/null
+++ b/test/testIntegrator.cpp
@@ -0,0 +1,198 @@
+#include "KokkosController.hpp"
+#include "MeshField.hpp"
+#include "MeshField_Element.hpp"    //remove?
+#include "MeshField_Fail.hpp"       //remove?
+#include "MeshField_For.hpp"        //remove?
+#include "MeshField_ShapeField.hpp" //remove?
+#include "Omega_h_build.hpp"
+#include "Omega_h_file.hpp"
+#include "Omega_h_simplex.hpp"
+#include <Kokkos_Core.hpp>
+#include <MeshField_Integrate.hpp>
+#include <iostream>
+#include <sstream>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
+
+template <typename FieldElement>
+class testIntegrator : public MeshField::Integrator {
+private:
+  testIntegrator(){};
+
+protected:
+  MeshField::Real integral;
+  FieldElement &fes;
+  int m;
+  int n;
+  int l;
+
+public:
+  MeshField::Real getIntegral() { return integral; }
+  testIntegrator(FieldElement &fes_in, int _m, int _n, int _l, int order)
+      : Integrator(order), integral(0), fes(fes_in), m(_m), n(_n), l(_l){};
+  void atPoints(Kokkos::View<MeshField::Real **> p,
+                Kokkos::View<MeshField::Real *> w,
+                Kokkos::View<MeshField::Real *> dV) {
+    const auto globalCoords = MeshField::evaluate(fes, p, p.extent(0));
+    const int _m = m;
+    const int _n = n;
+    const int _l = l;
+    Kokkos::parallel_reduce(
+        "integrate", p.extent(0),
+        KOKKOS_LAMBDA(const int &ent, MeshField::Real &integ) {
+          const auto x = globalCoords(ent, 0);
+          const auto y = globalCoords(ent, 1);
+          const auto z = globalCoords.extent(1) == 3 ? globalCoords(ent, 2) : 1;
+          integ += w(ent) * pow(x, _m) * pow(y, _n) * pow(z, _l) * dV(ent);
+        },
+        integral);
+  }
+};
+
+template <template <typename...> typename Controller, size_t ShapeOrder,
+          size_t dim>
+void doRun(Omega_h::Mesh &mesh,
+           MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
+  auto field =
+      omf.template CreateLagrangeField<MeshField::Real, ShapeOrder, dim>();
+  auto coords = mesh.coords();
+  Kokkos::parallel_for(
+      mesh.nverts(), KOKKOS_LAMBDA(int vtx) {
+        field(vtx, 0, 0, MeshField::Vertex) = coords[vtx * dim];
+        field(vtx, 0, 1, MeshField::Vertex) = coords[vtx * dim + 1];
+        if constexpr (dim == 3) {
+          field(vtx, 0, 2, MeshField::Vertex) = coords[vtx * dim + 2];
+        }
+      });
+  if (ShapeOrder == 2) {
+    auto edge2vtx = mesh.ask_down(1, 0).ab2b;
+    auto edgeMap = mesh.ask_down(dim, 1).ab2b;
+    Kokkos::parallel_for(
+        mesh.nedges(), KOKKOS_LAMBDA(int edge) {
+          const auto left = edge2vtx[edge * 2];
+          const auto right = edge2vtx[edge * 2 + 1];
+          const auto x = (coords[left * dim] + coords[right * dim]) / 2.0;
+          const auto y =
+              (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0;
+          field(edge, 0, 0, MeshField::Edge) = x;
+          field(edge, 0, 1, MeshField::Edge) = y;
+          if constexpr (dim == 3) {
+            const auto z =
+                (coords[left * dim + 2] + coords[right * dim + 2]) / 2.0;
+            field(edge, 0, 2, MeshField::Edge) = z;
+          }
+        });
+  }
+
+  auto shapeSet = [&]() {
+    if constexpr (dim == 3) {
+      return MeshField::Omegah::getTetrahedronElement<ShapeOrder>(mesh);
+    } else {
+      return MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
+    }
+  };
+  const auto [shp, map] = shapeSet();
+  MeshField::FieldElement fes(1, field, shp, map);
+  const int maxn = 32;
+  int binom[maxn + 1][maxn + 1];
+  for (int n = 0; n <= maxn; n++) {
+    binom[n][0] = binom[n][n] = 1;
+    for (int k = 1; k < n; k++) {
+      binom[n][k] = binom[n - 1][k] + binom[n - 1][k - 1];
+    }
+  }
+  if (dim == 2) {
+    for (int p = 1; p <= ShapeOrder; ++p) {
+      for (int m = p; m >= 0; --m) {
+        int n = p - m;
+        testIntegrator testInt(fes, m, n, 0, p);
+        testInt.process(fes);
+        double exact = 1.0 / binom[p][m] / (p + 1) / (p + 2);
+        if (Kokkos::fabs(testInt.getIntegral() - exact) >
+            MeshField::MachinePrecision) {
+          std::stringstream ss;
+          ss << "Integration over "
+             << ((dim == 3) ? "Tetrahedron " : "Triangle ") << "with order "
+             << ShapeOrder << " failed\n";
+          MeshField::fail(ss.str());
+        }
+      }
+    }
+  } else {
+    for (int p = 1; p <= ShapeOrder; ++p) {
+      for (int l = p; l >= 0; --l) {
+        for (int m = p - l; m >= 0; --m) {
+          int n = p - l - m;
+          testIntegrator testInt(fes, l, m, n, p);
+          testInt.process(fes);
+          double exact = 1.0 / binom[p][l + m] / binom[l + m][l] / (p + 1) /
+                         (p + 2) / (p + 3);
+          if (Kokkos::fabs(testInt.getIntegral() - exact) >
+              MeshField::MachinePrecision) {
+            std::stringstream ss;
+            ss << "Integration over "
+               << ((dim == 3) ? "Tetrahedron " : "Triangle ") << "with order "
+               << ShapeOrder << " failed\n";
+            MeshField::fail(ss.str());
+          }
+        }
+      }
+    }
+  }
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  auto lib = Omega_h::Library(&argc, &argv);
+#ifdef MESHFIELDS_ENABLE_CABANA
+  {
+    Omega_h::Reals coords2D({0.0, 0.0, 1.0, 0.0, 0.0, 1.0});
+    Omega_h::Reals coords3D(
+        {0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0});
+
+    Omega_h::LOs tris_to_verts({0, 1, 2});
+    Omega_h::LOs tets_to_verts({3, 0, 1, 2});
+    Omega_h::Mesh mesh2D(&lib);
+    Omega_h::Mesh mesh3D(&lib);
+    Omega_h::build_from_elems_and_coords(&mesh2D, OMEGA_H_SIMPLEX, 2,
+                                         tris_to_verts, coords2D);
+    Omega_h::build_from_elems_and_coords(&mesh3D, OMEGA_H_SIMPLEX, 3,
+                                         tets_to_verts, coords3D);
+
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::CabanaController>
+        omf2D(mesh2D);
+    doRun<MeshField::CabanaController, 1>(mesh2D, omf2D);
+    doRun<MeshField::CabanaController, 2>(mesh2D, omf2D);
+    MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::CabanaController>
+        omf3D(mesh3D);
+    doRun<MeshField::CabanaController, 1>(mesh3D, omf3D);
+    doRun<MeshField::CabanaController, 2>(mesh3D, omf3D);
+  }
+#endif
+  {
+    Omega_h::Reals coords2D({0.0, 0.0, 1.0, 0.0, 0.0, 1.0});
+    Omega_h::Reals coords3D(
+        {0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0});
+
+    Omega_h::LOs tris_to_verts({0, 1, 2});
+    Omega_h::LOs tets_to_verts({3, 0, 1, 2});
+    Omega_h::Mesh mesh2D(&lib);
+    Omega_h::Mesh mesh3D(&lib);
+    Omega_h::build_from_elems_and_coords(&mesh2D, OMEGA_H_SIMPLEX, 2,
+                                         tris_to_verts, coords2D);
+    Omega_h::build_from_elems_and_coords(&mesh3D, OMEGA_H_SIMPLEX, 3,
+                                         tets_to_verts, coords3D);
+
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController>
+        omf2D(mesh2D);
+    doRun<MeshField::KokkosController, 1>(mesh2D, omf2D);
+    doRun<MeshField::KokkosController, 2>(mesh2D, omf2D);
+    MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::KokkosController>
+        omf3D(mesh3D);
+    doRun<MeshField::KokkosController, 1>(mesh3D, omf3D);
+    doRun<MeshField::KokkosController, 2>(mesh3D, omf3D);
+  }
+  Kokkos::finalize();
+  return 0;
+}