diff --git a/CMakeLists.txt b/CMakeLists.txt
index b287e40..99ff8b3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -166,6 +166,7 @@ meshfields_add_exe(SerializationTests test/testSerialize.cpp)
 meshfields_add_exe(ElementTests test/testElement.cpp)
 meshfields_add_exe(ElementJacobian1d test/testElementJacobian1d.cpp)
 meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
+meshfields_add_exe(ElementJacobian3d test/testElementJacobian3d.cpp)
 meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
 meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
@@ -184,6 +185,7 @@ test_func(SerializationTests ./SerializationTests)
 test_func(ElementTests ./ElementTests)
 test_func(ElementJacobian1d ./ElementJacobian1d)
 test_func(ElementJacobian2d ./ElementJacobian2d)
+test_func(ElementJacobian3d ./ElementJacobian3d)
 test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index 63b51da..78757b6 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -71,8 +71,8 @@ struct LinearTriangleToVertexField {
     }
   }
 
-  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
-  getTopology() const {
+  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() {
     return {MeshField::Triangle};
   }
 
@@ -102,8 +102,8 @@ struct LinearTetrahedronToVertexField {
           __func__);
     }
   }
-  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
-  getTopology() const {
+  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() {
     return {MeshField::Tetrahedron};
   }
 
@@ -136,8 +136,8 @@ struct QuadraticTriangleToField {
     }
   }
 
-  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
-  getTopology() const {
+  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() {
     return {MeshField::Triangle};
   }
 
@@ -198,8 +198,8 @@ struct QuadraticTetrahedronToField {
     }
   }
 
-  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
-  getTopology() const {
+  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() {
     return {MeshField::Tetrahedron};
   }
 
diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp
index b43c500..f29c7d8 100644
--- a/src/MeshField_Element.hpp
+++ b/src/MeshField_Element.hpp
@@ -256,9 +256,17 @@ struct FieldElement {
          "tangent vectors" in 3D, the volume of their
          parallelpiped, which is the differential volume
          of the coordinate field */
-      fail("getJacobianDeterminant doesn't yet support 3d.  "
-           "The given matrices have dimension %lu x %lu \n",
-           J.extent(0), J.extent(1));
+      Kokkos::View<Real *> determinants("3dJacobianDeterminants", J.extent(0));
+      Kokkos::parallel_for(
+          J.extent(0), KOKKOS_LAMBDA(const int i) {
+            const auto Ji = Kokkos::subview(J, i, Kokkos::ALL(), Kokkos::ALL());
+            const auto cofactorSum =
+                Ji(0, 0) * (Ji(1, 1) * Ji(2, 2) - Ji(1, 2) * Ji(2, 1)) -
+                Ji(0, 1) * (Ji(1, 0) * Ji(2, 2) - Ji(1, 2) * Ji(2, 0)) +
+                Ji(0, 2) * (Ji(1, 0) * Ji(2, 1) - Ji(1, 1) * Ji(2, 0));
+            determinants(i) = cofactorSum;
+          });
+      return determinants;
     }
     if (dimension == 2) {
       if (J.extent(1) != 2) {
@@ -349,8 +357,9 @@ struct FieldElement {
       fail("The input array of offsets must have size = %zu\n",
            numMeshEnts + 1);
     }
-    if (MeshEntDim != 1 && MeshEntDim != 2) {
-      fail("getJacobians only currently supports 1d and 2d meshes.  Input mesh "
+    if (MeshEntDim != 1 && MeshEntDim != 2 && MeshEntDim != 3) {
+      fail("getJacobians only currently supports 1d, 2d, and 3d meshes.  Input "
+           "mesh "
            "has %zu dimensions.\n",
            numMeshEnts);
     }
@@ -366,7 +375,7 @@ struct FieldElement {
             }
           });
       return res;
-    } else if constexpr (MeshEntDim == 2) {
+    } else if constexpr (MeshEntDim == 2 || MeshEntDim == 3) {
       const auto numPts = MeshFieldUtil::getLastValue(offsets);
       // one matrix per point
       Kokkos::View<Real ***> res("result", numPts, MeshEntDim, MeshEntDim);
diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp
index 15fa8bc..0c16d68 100644
--- a/src/MeshField_Integrate.hpp
+++ b/src/MeshField_Integrate.hpp
@@ -10,50 +10,51 @@
 
 namespace MeshField {
 // directly copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.[h|cc]
-struct IntegrationPoint {
-  IntegrationPoint(Vector3 const &p, double w) : param(p), weight(w) {}
-  Vector3 param;
+template <size_t pointSize> struct IntegrationPoint {
+  // template parameter pointSize specifies the length of the integration point
+  // array for one point
+  IntegrationPoint(Kokkos::Array<Real, pointSize> const &p, double w)
+      : param(p), weight(w) {}
+  Kokkos::Array<Real, pointSize> param;
   double weight;
 };
-
-class Integration {
+template <size_t pointSize> class Integration {
 public:
   virtual ~Integration() {}
   virtual int countPoints() const = 0;
-  virtual std::vector<IntegrationPoint> getPoints() const = 0;
+  virtual std::vector<IntegrationPoint<pointSize>> getPoints() const = 0;
   virtual int getAccuracy() const = 0;
 };
-
-class EntityIntegration {
+template <size_t pointSize> class EntityIntegration {
 public:
   virtual ~EntityIntegration() {}
-  Integration const *getAccurate(int minimumAccuracy) const {
+  Integration<pointSize> const *getAccurate(int minimumAccuracy) const {
     int n = countIntegrations();
     for (int i = 0; i < n; ++i) {
-      Integration const *integration = getIntegration(i);
+      Integration<pointSize> const *integration = getIntegration(i);
       if (integration->getAccuracy() >= minimumAccuracy)
         return integration;
     }
     return NULL;
   }
   virtual int countIntegrations() const = 0;
-  virtual Integration const *getIntegration(int i) const = 0;
+  virtual Integration<pointSize> const *getIntegration(int i) const = 0;
 };
 
-class TriangleIntegration : public EntityIntegration {
+class TriangleIntegration : public EntityIntegration<3> {
 public:
-  class N1 : public Integration {
+  class N1 : public Integration<3> {
   public:
     virtual int countPoints() const { return 1; }
-    virtual std::vector<IntegrationPoint> getPoints() const {
+    virtual std::vector<IntegrationPoint<3>> getPoints() const {
       return {IntegrationPoint(Vector3{1. / 3., 1. / 3., 1. / 3.}, 1.0 / 2.0)};
     }
     virtual int getAccuracy() const { return 1; }
   }; // end N1
-  class N2 : public Integration {
+  class N2 : public Integration<3> {
   public:
     virtual int countPoints() const { return 3; }
-    virtual std::vector<IntegrationPoint> getPoints() const {
+    virtual std::vector<IntegrationPoint<3>> getPoints() const {
       return {IntegrationPoint(Vector3{0.666666666666667, 0.166666666666667,
                                        0.166666666666667},
                                1. / 3. / 2.0),
@@ -67,31 +68,69 @@ class TriangleIntegration : public EntityIntegration {
     virtual int getAccuracy() const { return 2; }
   }; // end N2
   virtual int countIntegrations() const { return 2; }
-  virtual Integration const *getIntegration(int i) const {
+  virtual Integration<3> const *getIntegration(int i) const {
     static N1 i1;
     static N2 i2;
-    static Integration *integrations[2] = {&i1, &i2};
+    static Integration<3> *integrations[2] = {&i1, &i2};
     return integrations[i];
   }
 };
 
-std::shared_ptr<EntityIntegration> const getIntegration(Mesh_Topology topo) {
-  if (topo != Triangle) {
-    fail("getIntegration only supports triangles (Mesh_Topology::Triangle)\n");
+class TetrahedronIntegration : public EntityIntegration<4> {
+public:
+  class N1 : public Integration<4> {
+  public:
+    virtual int countPoints() const { return 1; }
+    virtual std::vector<IntegrationPoint<4>> getPoints() const {
+      return {IntegrationPoint(Vector4{0.25, 0.25, 0.25, 0.25}, 1.0 / 6.0)};
+    }
+    virtual int getAccuracy() const { return 1; }
+  };
+  class N2 : public Integration<4> {
+  public:
+    virtual int countPoints() const { return 4; }
+    virtual std::vector<IntegrationPoint<4>> getPoints() const {
+
+      return {IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
+                                       0.138196601125011, 0.585410196624969},
+                               0.25 / 6.0),
+              IntegrationPoint(Vector4{0.585410196624969, 0.138196601125011,
+                                       0.138196601125011, 0.138196601125011},
+                               0.25 / 6.0),
+              IntegrationPoint(Vector4{0.138196601125011, 0.585410196624969,
+                                       0.138196601125011, 0.138196601125011},
+                               0.25 / 6.0),
+              IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
+                                       0.585410196624969, 0.138196601125011},
+                               0.25 / 6.0)};
+    }
+    virtual int getAccuracy() const { return 2; }
+  };
+  virtual int countIntegrations() const { return 2; }
+  virtual Integration<4> const *getIntegration(int i) const {
+    static N1 i1;
+    static N2 i2;
+    static Integration<4> *integrations[2] = {&i1, &i2};
+    return integrations[i];
   }
-  return std::make_shared<TriangleIntegration>();
+};
+template <Mesh_Topology topo> auto const getIntegration() {
+  if constexpr (topo == Triangle) {
+    return std::make_shared<TriangleIntegration>();
+  } else if constexpr (topo == Tetrahedron) {
+    return std::make_shared<TetrahedronIntegration>();
+  }
+  fail("getIntegration does not support given topology\n");
 }
-
-std::vector<IntegrationPoint> getIntegrationPoints(Mesh_Topology topo,
-                                                   int order) {
-  auto ip = getIntegration(topo)->getAccurate(order)->getPoints();
+template <Mesh_Topology topo> auto getIntegrationPoints(int order) {
+  auto ip = getIntegration<topo>()->getAccurate(order)->getPoints();
   return ip;
 }
 
 template <typename FieldElement>
-Kokkos::View<MeshField::Real **>
-getIntegrationPointLocalCoords(FieldElement &fes,
-                               std::vector<IntegrationPoint> ip) {
+Kokkos::View<MeshField::Real **> getIntegrationPointLocalCoords(
+    FieldElement &fes,
+    std::vector<IntegrationPoint<FieldElement::MeshEntDim + 1>> ip) {
   const auto numPtsPerElm = ip.size();
   const auto numMeshEnts = fes.numMeshEnts;
   const auto meshEntDim = fes.MeshEntDim;
@@ -115,9 +154,9 @@ getIntegrationPointLocalCoords(FieldElement &fes,
 }
 
 template <typename FieldElement>
-Kokkos::View<Real *>
-getIntegrationPointWeights(FieldElement &fes,
-                           std::vector<IntegrationPoint> ip) {
+Kokkos::View<Real *> getIntegrationPointWeights(
+    FieldElement &fes,
+    std::vector<IntegrationPoint<FieldElement::MeshEntDim + 1>> ip) {
   const auto numPtsPerElm = ip.size();
   const auto numMeshEnts = fes.numMeshEnts;
   const auto meshEntDim = fes.MeshEntDim;
@@ -194,13 +233,9 @@ class Integrator {
    * FIXME make the sensible
    * */
   template <typename FieldElement> void process(FieldElement &fes) {
-    const auto topo = fes.elm2dof.getTopology();
-    if (topo[0] != MeshField::Triangle) {
-      throw std::invalid_argument(
-          "Integrator::process only supports triangles.");
-    }
+    constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology();
     pre();
-    auto ip = getIntegrationPoints(topo[0], order);
+    auto ip = getIntegrationPoints<topo[0]>(order);
     auto localCoords = getIntegrationPointLocalCoords(fes, ip);
     auto weights = getIntegrationPointWeights(fes, ip);
     auto dV = getJacobianDeterminants(fes, localCoords, ip.size());
diff --git a/test/testCountIntegrator.cpp b/test/testCountIntegrator.cpp
index 68aab23..c6361b0 100644
--- a/test/testCountIntegrator.cpp
+++ b/test/testCountIntegrator.cpp
@@ -14,12 +14,12 @@
 
 using ExecutionSpace = Kokkos::DefaultExecutionSpace;
 using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
-
-Omega_h::Mesh createMeshTri18(Omega_h::Library &lib) {
+template <size_t dim> Omega_h::Mesh createMesh(Omega_h::Library &lib) {
   auto world = lib.world();
   const auto family = OMEGA_H_SIMPLEX;
   auto len = 1.0;
-  return Omega_h::build_box(world, family, len, len, 0.0, 3, 3, 0);
+  const auto numEnts3d = (dim == 3 ? 3 : 0);
+  return Omega_h::build_box(world, family, len, len, len, 3, 3, numEnts3d);
 }
 
 template <typename AnalyticFunction, typename ShapeField>
@@ -57,13 +57,19 @@ class CountIntegrator : public MeshField::Integrator {
   }
 };
 
-template <template <typename...> typename Controller>
+template <template <typename...> typename Controller, size_t dim>
 void doRun(Omega_h::Mesh &mesh,
-           MeshField::OmegahMeshField<ExecutionSpace, 2, Controller> &omf) {
+           MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
   const auto ShapeOrder = 1;
   auto field = omf.getCoordField();
-  const auto [shp, map] =
-      MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
+  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);
 
   CountIntegrator countInt(fes);
@@ -74,18 +80,28 @@ void doRun(Omega_h::Mesh &mesh,
 int main(int argc, char **argv) {
   Kokkos::initialize(argc, argv);
   auto lib = Omega_h::Library(&argc, &argv);
-  auto mesh = createMeshTri18(lib);
 #ifdef MESHFIELDS_ENABLE_CABANA
   {
+    auto mesh2D = createMesh<2>(lib);
+    auto mesh3D = createMesh<3>(lib);
     MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::CabanaController>
-        omf(mesh);
-    doRun<MeshField::CabanaController>(mesh, omf);
+        omf2D(mesh2D);
+    doRun<MeshField::CabanaController>(mesh2D, omf2D);
+    MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::CabanaController>
+        omf3D(mesh3D);
+    doRun<MeshField::CabanaController>(mesh3D, omf3D);
   }
 #endif
   {
+    auto mesh2D = createMesh<2>(lib);
+    auto mesh3D = createMesh<3>(lib);
     MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController>
-        omf(mesh);
-    doRun<MeshField::KokkosController>(mesh, omf);
+        omf2D(mesh2D);
+    doRun<MeshField::KokkosController>(mesh2D, omf2D);
+    MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::KokkosController>
+        omf3D(mesh3D);
+    doRun<MeshField::KokkosController>(mesh3D, omf3D);
   }
+  Kokkos::finalize();
   return 0;
 }
diff --git a/test/testElementJacobian2d.cpp b/test/testElementJacobian2d.cpp
index 1025f8c..733e56d 100644
--- a/test/testElementJacobian2d.cpp
+++ b/test/testElementJacobian2d.cpp
@@ -84,7 +84,7 @@ void triJacobian() {
                               LinearTriangleToVertexField());
 
     Kokkos::View<MeshField::Real *[3]> lc("localCoords", 1);
-    Kokkos::deep_copy(lc, 1.0 / 2);
+    Kokkos::deep_copy(lc, 1.0 / 3);
     const auto numPtsPerElement = 1;
     const auto J = MeshField::getJacobians(f, lc, numPtsPerElement);
     const auto J_h =
diff --git a/test/testElementJacobian3d.cpp b/test/testElementJacobian3d.cpp
new file mode 100644
index 0000000..6803433
--- /dev/null
+++ b/test/testElementJacobian3d.cpp
@@ -0,0 +1,118 @@
+#include "KokkosController.hpp"
+#include "MeshField_Element.hpp"
+#include "MeshField_Field.hpp"
+#include "MeshField_For.hpp"
+#include "MeshField_ShapeField.hpp"
+#include <Kokkos_Core.hpp>
+#include <iostream>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
+
+struct LinearTetrahedronToVertexField {
+  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() const {
+    return {MeshField::Tetrahedron};
+  }
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx,
+             MeshField::LO tet, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Tetrahedron);
+    MeshField::LO tetNode2Vtx[1][4] = {{0, 1, 2, 3}};
+    const MeshField::LO vtx = tetNode2Vtx[tet][tetNodeIdx];
+    return {0, tetCompIdx, vtx, MeshField::Vertex};
+  }
+};
+
+class TetrahedronTestCase {
+public:
+  TetrahedronTestCase(std::vector<MeshField::Real> coords_in,
+                      std::vector<MeshField::Real> jacobian_in,
+                      MeshField::Real jacobianDeterminant_in)
+      : coords(coords_in), jacobian(jacobian_in),
+        jacobianDeterminant(jacobianDeterminant_in) {
+    assert(coords.size() == 12);
+    assert(jacobian.size() == 9);
+  };
+  const std::vector<MeshField::Real> coords;
+  const std::vector<MeshField::Real> jacobian;
+  const MeshField::Real jacobianDeterminant;
+
+private:
+  TetrahedronTestCase();
+};
+
+template <typename ShapeField>
+void setVtxCoords(size_t numVerts, size_t meshDim, TetrahedronTestCase testTri,
+                  ShapeField field) {
+  Kokkos::View<MeshField::Real *, Kokkos::HostSpace> coords_h(
+      "coords_h", numVerts * meshDim);
+  for (size_t i = 0; i < numVerts * meshDim; i++)
+    coords_h[i] = testTri.coords.at(i);
+  auto coords = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), coords_h);
+  auto setCoordField = KOKKOS_LAMBDA(const int &vtx) {
+    field(vtx, 0, 0, MeshField::Vertex) = coords(vtx * meshDim);
+    field(vtx, 0, 1, MeshField::Vertex) = coords(vtx * meshDim + 1);
+    field(vtx, 0, 2, MeshField::Vertex) = coords(vtx * meshDim + 2);
+  };
+  MeshField::parallel_for(ExecutionSpace(), {0}, {numVerts}, setCoordField,
+                          "setCoordField");
+}
+
+void tetJacobian() {
+  const MeshField::MeshInfo meshInfo{
+      .numVtx = 4, .numEdge = 6, .numTet = 1, .dim = 3};
+  auto coordField = MeshField::CreateCoordinateField<
+      ExecutionSpace, MeshField::KokkosController, 3>(meshInfo);
+  TetrahedronTestCase identityTet({0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1},
+                                  {1, 0, 0, 0, 1, 0, 0, 0, 1}, 1);
+  TetrahedronTestCase skewedTet({0, 0, 0, 5, 1, 2, 3, 4, 6, 10, 8, 9},
+                                {5, 1, 2, 3, 4, 6, 10, 8, 9}, -59);
+  for (auto testCase : {identityTet, skewedTet}) {
+    setVtxCoords(meshInfo.numVtx, meshInfo.dim, testCase, coordField);
+
+    MeshField::FieldElement f(meshInfo.numTet, coordField,
+                              MeshField::LinearTetrahedronShape(),
+                              LinearTetrahedronToVertexField());
+
+    Kokkos::View<MeshField::Real *[4]> lc("localCoords", 1);
+    Kokkos::deep_copy(lc, 1.0 / 4);
+    const auto numPtsPerElement = 1;
+    const auto J = MeshField::getJacobians(f, lc, numPtsPerElement);
+    const auto J_h =
+        Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J);
+    assert(J_h.rank() == 3);
+    assert(J_h.extent(0) == 1);
+    assert(J_h.extent(1) == 3);
+    assert(J_h.extent(2) == 3);
+    std::cout << "tet jacobian\n"
+              << J_h(0, 0, 0) << " " << J_h(0, 0, 1) << " " << J_h(0, 0, 2)
+              << "\n"
+              << J_h(0, 1, 0) << " " << J_h(0, 1, 1) << " " << J_h(0, 1, 2)
+              << "\n"
+              << J_h(0, 2, 0) << " " << J_h(0, 2, 1) << " " << J_h(0, 2, 2)
+              << "\n";
+    const auto determinants = MeshField::getJacobianDeterminants(f, J);
+    const auto determinants_h =
+        Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), determinants);
+    std::cout << "tet jacobian determinant " << determinants_h(0) << "\n";
+    for (int i = 0; i < 3; ++i) {
+      for (int j = 0; j < 3; ++j) {
+        assert(std::fabs(J_h(0, i, j) - testCase.jacobian.at(i * 3 + j)) <=
+               MeshField::MachinePrecision);
+      }
+    }
+    assert(determinants.rank() == 1);
+    assert(determinants.extent(0) == 1);
+    assert(std::fabs(determinants_h(0) - testCase.jacobianDeterminant) <=
+           MeshField::MachinePrecision);
+  }
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  tetJacobian();
+  Kokkos::finalize();
+  return 0;
+}