From 01a16db62cf9c6f893dc191d9ed3c9fa806ca091 Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 24 Oct 2025 13:12:05 -0400
Subject: [PATCH 1/8] init

---
 src/MeshField.hpp                             | 66 ++++++++++++++-
 src/MeshField_Integrate.hpp                   | 39 ++++++++-
 src/MeshField_Shape.hpp                       | 84 +++++++++++++++++++
 test/testReducedQuinticImplicitIntegrator.cpp | 80 ++++++++++++++++++
 4 files changed, 265 insertions(+), 4 deletions(-)
 create mode 100644 test/testReducedQuinticImplicitIntegrator.cpp

diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index 78757b6..453c1d6 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -238,6 +238,61 @@ struct QuadraticTetrahedronToField {
   }
 };
 
+struct ReducedQuinticTriangleToField {
+  Omega_h::LOs triVerts;
+  ReducedQuinticTriangleToField(Omega_h::Mesh &mesh)
+      : triVerts(mesh.ask_elem_verts()) {
+    if (mesh.dim() != 2 && mesh.family() != OMEGA_H_SIMPLEX) {
+      MeshField::fail(
+          "The mesh passed to %s must be 2D and simplex (triangles)\n",
+          __func__);
+    }
+  }
+
+  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() {
+    return {MeshField::Triangle};
+  }
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx,
+             MeshField::LO tri, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Triangle);
+    const MeshField::LO triNode2DofHolder[18] = {
+        0,0,0,0,0,0,  // vertex 0 DOFs (6)
+        1,1,1,1,1,1,  // vertex 1 DOFs (6)
+        2,2,2,2,2,2   // vertex 2 DOFs (6)
+    };
+    const MeshField::Mesh_Topology triNode2DofHolderTopo[18] = {
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex
+    };
+    const auto dofHolderIdx = triNode2DofHolder[triNodeIdx];
+    const auto dofHolderTopo = triNode2DofHolderTopo[triNodeIdx];
+
+    Omega_h::LO osh_ent;
+    if (dofHolderTopo == MeshField::Vertex) {
+      const auto triDim = 2;
+      const auto vtxDim = 0;
+      const auto ignored = -1;
+      const auto localVtxIdx = (Omega_h::simplex_down_template(
+                                    triDim, vtxDim, dofHolderIdx, ignored) +
+                                2) %
+                               3;
+      const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
+      osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx];
+    } else {
+      assert(false);
+    }
+
+    return {0, triCompIdx, osh_ent, dofHolderTopo};
+  }
+};
+
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
   static_assert(ShapeOrder == 1 || ShapeOrder == 2);
   if constexpr (ShapeOrder == 1) {
@@ -254,6 +309,13 @@ template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
     };
     return result{MeshField::QuadraticTriangleShape(),
                   QuadraticTriangleToField(mesh)};
+  } else if constexpr (ShapeOrder == 5) {
+    struct result {
+      MeshField::ReducedQuinticImplicitShape shp;
+      ReducedQuinticTriangleToField map;
+    };
+    return result{MeshField::ReducedQuinticImplicitShape(),
+                  ReducedQuinticTriangleToField(mesh)};
   }
 }
 template <int ShapeOrder> auto getTetrahedronElement(Omega_h::Mesh &mesh) {
@@ -346,8 +408,8 @@ class OmegahMeshField {
       MeshField::fail("input mesh must be 2d\n");
     }
     const auto ShapeOrder = ShapeField::Order;
-    if (ShapeOrder != 1 && ShapeOrder != 2) {
-      MeshField::fail("input field order must be 1 or 2\n");
+    if (ShapeOrder != 1 && ShapeOrder != 2 && ShapeOrder != 5) {
+      MeshField::fail("input field order must be 1 or 2 or 5\n");
     }
 
     const auto [shp, map] = Omegah::getTriangleElement<ShapeOrder>(mesh);
diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp
index 0c16d68..20004ee 100644
--- a/src/MeshField_Integrate.hpp
+++ b/src/MeshField_Integrate.hpp
@@ -67,11 +67,46 @@ class TriangleIntegration : public EntityIntegration<3> {
     }
     virtual int getAccuracy() const { return 2; }
   }; // end N2
-  virtual int countIntegrations() const { return 2; }
+    class N6 : public Integration<3> {
+  public:
+    virtual int countPoints() const { return 12; }
+    virtual std::vector<IntegrationPoint<3>> getPoints() const {
+      return {
+          IntegrationPoint(Vector3{0.873821971, 0.063089014, 0.063089014},
+                           0.0508449063702 / 2.0),
+          IntegrationPoint(Vector3{0.063089014, 0.873821971, 0.063089014},
+                           0.0508449063702 / 2.0),
+          IntegrationPoint(Vector3{0.063089014, 0.063089014, 0.873821971},
+                           0.0508449063702 / 2.0),
+
+          IntegrationPoint(Vector3{0.501426509, 0.249286745, 0.249286745},
+                           0.1167862757264 / 2.0),
+          IntegrationPoint(Vector3{0.249286745, 0.501426509, 0.249286745},
+                           0.1167862757264 / 2.0),
+          IntegrationPoint(Vector3{0.249286745, 0.249286745, 0.501426509},
+                           0.1167862757264 / 2.0),
+
+          IntegrationPoint(Vector3{0.636502499, 0.310352451, 0.053145050},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.636502499, 0.053145050, 0.310352451},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.310352451, 0.636502499, 0.053145050},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.310352451, 0.053145050, 0.636502499},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.053145050, 0.636502499, 0.310352451},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.053145050, 0.310352451, 0.636502499},
+                           0.0828510756184 / 2.0)};
+    }
+    virtual int getAccuracy() const { return 6; }
+  }; // end N6
+  virtual int countIntegrations() const { return 3; }
   virtual Integration<3> const *getIntegration(int i) const {
     static N1 i1;
     static N2 i2;
-    static Integration<3> *integrations[2] = {&i1, &i2};
+    static N6 i6;
+    static Integration<3> *integrations[3] = {&i1, &i2, &i6};
     return integrations[i];
   }
 };
diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp
index ceafda5..faa71b5 100644
--- a/src/MeshField_Shape.hpp
+++ b/src/MeshField_Shape.hpp
@@ -217,5 +217,89 @@ struct QuadraticTetrahedronShape {
   }
 };
 
+struct ReducedQuinticImplicitShape {
+  static const size_t numNodes = 18;
+  static const size_t meshEntDim = 2;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 5;
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+
+    const Real L1 = 1.0 - xi[0] - xi[1];
+    const Real L2 = xi[0];
+    const Real L3 = xi[1];
+
+    const Real L1_2 = L1 * L1;
+    const Real L2_2 = L2 * L2;
+    const Real L3_2 = L3 * L3;
+    const Real L1_3 = L1_2 * L1;
+    const Real L2_3 = L2_2 * L2;
+    const Real L3_3 = L3_2 * L3;
+    const Real L1_4 = L1_3 * L1;
+    const Real L2_4 = L2_3 * L2;
+    const Real L3_4 = L3_3 * L3;
+    const Real L1_5 = L1_4 * L1;
+    const Real L2_5 = L2_4 * L2;
+    const Real L3_5 = L3_4 * L3;
+
+    Kokkos::Array<Real, numNodes> N;
+
+    N[0] = L1 * (2 * L1 - 1) * (3 * L1 - 1) * (4 * L1 - 1) * (5 * L1 - 1);
+    N[1] = L2 * (2 * L2 - 1) * (3 * L2 - 1) * (4 * L2 - 1) * (5 * L2 - 1);
+    N[2] = L3 * (2 * L3 - 1) * (3 * L3 - 1) * (4 * L3 - 1) * (5 * L3 - 1);
+    
+    N[3] = 5.0 * L1_4 * L2;
+    N[4] = 5.0 * L2_4 * L1;
+    N[5] = 5.0 * L3_4 * L1;
+
+    N[6] = 5.0 * L1_4 * L3;
+    N[7] = 5.0 * L2_4 * L3;
+    N[8] = 5.0 * L3_4 * L2;
+
+    N[9]  = 10.0 * L1_3 * L2_2;
+    N[10] = 10.0 * L2_3 * L1_2;
+    N[11] = 10.0 * L3_3 * L1_2;
+
+    N[12] = 10.0 * L1_3 * L3_2;
+    N[13] = 10.0 * L2_3 * L3_2;
+    N[14] = 10.0 * L3_3 * L2_2;
+
+    N[15] = 20.0 * L1_2 * L2 * L3;
+    N[16] = 20.0 * L2_2 * L1 * L3;
+    N[17] = 20.0 * L3_2 * L1 * L2;
+
+    return N;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Vector2, numNodes> getLocalGradients(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+
+    const Real L1 = 1.0 - xi[0] - xi[1];
+    const Real L2 = xi[0];
+    const Real L3 = xi[1];
+
+    const Real eps = 1e-6;
+
+    Vector3 xp = {xi[0] + eps, xi[1], 0.0};
+    Vector3 yp = {xi[0], xi[1] + eps, 0.0};
+
+    auto N0 = getValues(xi);
+    auto Nx = getValues(xp);
+    auto Ny = getValues(yp);
+
+    Kokkos::Array<Vector2, numNodes> dN;
+    for (int i = 0; i < numNodes; i++) {
+      dN[i][0] = (Nx[i] - N0[i]) / eps;
+      dN[i][1] = (Ny[i] - N0[i]) / eps;
+    }
+    return dN;
+  }
+};
+
 } // namespace MeshField
 #endif
diff --git a/test/testReducedQuinticImplicitIntegrator.cpp b/test/testReducedQuinticImplicitIntegrator.cpp
new file mode 100644
index 0000000..6206337
--- /dev/null
+++ b/test/testReducedQuinticImplicitIntegrator.cpp
@@ -0,0 +1,80 @@
+#include "KokkosController.hpp"
+#include "MeshField.hpp"
+#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>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+
+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, len,
+                            2, 2, (dim == 3 ? 2 : 0));
+}
+
+template <typename FieldElement>
+class ReducedQuinticImplicitIntegrator : public MeshField::Integrator {
+public:
+  ReducedQuinticImplicitIntegrator(FieldElement &fes_in)
+      : MeshField::Integrator(5), fes(fes_in), totalValue(0.0) {}
+
+  void atPoints(Kokkos::View<MeshField::Real **> p,
+                Kokkos::View<MeshField::Real *> w,
+                Kokkos::View<MeshField::Real *> dV) override {
+    const auto numPts = p.extent(0);
+    double localSum = 0.0;
+    Kokkos::parallel_reduce(
+        "IntegrateReducedQuinticImplicit", numPts,
+        KOKKOS_LAMBDA(const int i, double &sum) {
+          const double weight = w(i);
+          const double jac = dV(i);
+          const double f = 1.0; // Integrate constant field
+          sum += f * weight * jac;
+        },
+        localSum);
+    totalValue += localSum;
+  }
+
+  void post() override {
+    printf("[ReducedQuinticImplicit] Integrated Value = %e\n", totalValue);
+  }
+
+private:
+  FieldElement &fes;
+  double totalValue;
+};
+
+template <template <typename...> typename Controller, size_t dim>
+void runReducedQuinticImplicit(Omega_h::Mesh &mesh,
+                               MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
+  const auto ShapeOrder = 5;
+
+  auto field = omf.getCoordField();
+
+  const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
+  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
+
+  ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);
+  integrator.process(fes);
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  Omega_h::Library lib(&argc, &argv);
+
+  {
+    auto mesh2D = createMesh<2>(lib);
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
+    runReducedQuinticImplicit<MeshField::KokkosController>(mesh2D, omf2D);
+  }
+
+  Kokkos::finalize();
+  return 0;
+}
+

From 5dad2c9a442ad8fbde949b9e7ccbd2684f8373f7 Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Wed, 5 Nov 2025 22:16:11 -0500
Subject: [PATCH 2/8] test evaluation

---
 CMakeLists.txt                                |  4 +
 src/MeshField.hpp                             | 12 ++-
 src/MeshField_Integrate.hpp                   | 63 ++++++++-------
 src/MeshField_Shape.hpp                       | 20 +++--
 test/testCountIntegrator.cpp                  |  2 +-
 test/testReducedQuinticImplicitEvaluate.cpp   | 80 +++++++++++++++++++
 test/testReducedQuinticImplicitIntegrator.cpp | 39 ++++++---
 7 files changed, 173 insertions(+), 47 deletions(-)
 create mode 100644 test/testReducedQuinticImplicitEvaluate.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 99ff8b3..99dcf17 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(ReducedQuinticImplicitIntegratorTest test/testReducedQuinticImplicitIntegrator.cpp)
+meshfields_add_exe(ReducedQuinticImplicitEvaluateTest test/testReducedQuinticImplicitEvaluate.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -190,6 +192,8 @@ test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
 test_func(OmegahTetTest, ./OmegahTetTest)
+test_func(ReducedQuinticImplicitIntegratorTest ./ReducedQuinticImplicitIntegratorTest)
+test_func(ReducedQuinticImplicitEvaluateTest ./ReducedQuinticImplicitEvaluateTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)
diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index 453c1d6..bf20104 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -9,6 +9,7 @@
 #include "Omega_h_file.hpp"    //move
 #include "Omega_h_mesh.hpp"    //move
 #include "Omega_h_simplex.hpp" //move
+#include "MeshField_Shape.hpp"
 
 namespace {
 
@@ -294,7 +295,7 @@ struct ReducedQuinticTriangleToField {
 };
 
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
-  static_assert(ShapeOrder == 1 || ShapeOrder == 2);
+  static_assert(ShapeOrder == 1 || ShapeOrder == 2 || ShapeOrder == 5);
   if constexpr (ShapeOrder == 1) {
     struct result {
       MeshField::LinearTriangleShape shp;
@@ -337,6 +338,15 @@ template <int ShapeOrder> auto getTetrahedronElement(Omega_h::Mesh &mesh) {
   }
 }
 
+inline auto getReducedQuinticImplicitElement(Omega_h::Mesh &mesh) {
+  struct result {
+    MeshField::ReducedQuinticImplicitShape shp;
+    ReducedQuinticTriangleToField map;
+  };
+  return result{MeshField::ReducedQuinticImplicitShape(),
+                ReducedQuinticTriangleToField(mesh)};
+}
+
 } // namespace Omegah
 
 template <typename ExecutionSpace, size_t dim,
diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp
index 20004ee..7462f78 100644
--- a/src/MeshField_Integrate.hpp
+++ b/src/MeshField_Integrate.hpp
@@ -71,33 +71,42 @@ class TriangleIntegration : public EntityIntegration<3> {
   public:
     virtual int countPoints() const { return 12; }
     virtual std::vector<IntegrationPoint<3>> getPoints() const {
-      return {
-          IntegrationPoint(Vector3{0.873821971, 0.063089014, 0.063089014},
-                           0.0508449063702 / 2.0),
-          IntegrationPoint(Vector3{0.063089014, 0.873821971, 0.063089014},
-                           0.0508449063702 / 2.0),
-          IntegrationPoint(Vector3{0.063089014, 0.063089014, 0.873821971},
-                           0.0508449063702 / 2.0),
-
-          IntegrationPoint(Vector3{0.501426509, 0.249286745, 0.249286745},
-                           0.1167862757264 / 2.0),
-          IntegrationPoint(Vector3{0.249286745, 0.501426509, 0.249286745},
-                           0.1167862757264 / 2.0),
-          IntegrationPoint(Vector3{0.249286745, 0.249286745, 0.501426509},
-                           0.1167862757264 / 2.0),
-
-          IntegrationPoint(Vector3{0.636502499, 0.310352451, 0.053145050},
-                           0.0828510756184 / 2.0),
-          IntegrationPoint(Vector3{0.636502499, 0.053145050, 0.310352451},
-                           0.0828510756184 / 2.0),
-          IntegrationPoint(Vector3{0.310352451, 0.636502499, 0.053145050},
-                           0.0828510756184 / 2.0),
-          IntegrationPoint(Vector3{0.310352451, 0.053145050, 0.636502499},
-                           0.0828510756184 / 2.0),
-          IntegrationPoint(Vector3{0.053145050, 0.636502499, 0.310352451},
-                           0.0828510756184 / 2.0),
-          IntegrationPoint(Vector3{0.053145050, 0.310352451, 0.636502499},
-                           0.0828510756184 / 2.0)};
+      return {IntegrationPoint(Vector3{0.873821971016996,0.063089014491502,
+		      		       0},
+			       0.050844906370207 / 2.0),
+  	      IntegrationPoint(Vector3{0.063089014491502,0.873821971016996,
+			      	       0},
+			       0.050844906370207/2.0),
+  	      IntegrationPoint(Vector3{0.063089014491502,0.063089014491502,
+			      	       0},
+			       0.050844906370207/2.0),
+  	      IntegrationPoint(Vector3{0.501426509658179,0.249286745170910,
+			               0},
+			       0.116786275726379/2.0),
+  	      IntegrationPoint(Vector3{0.249286745170910,0.501426509658179,
+			      	       0},
+			      0.116786275726379/2.0),
+  	      IntegrationPoint(Vector3{0.249286745170910,0.249286745170910,
+			               0},
+		              0.116786275726379/2.0),
+  	      IntegrationPoint(Vector3{0.636502499121399,0.310352451033785,
+			               0},
+			      0.082851075618374/2.0),
+  	      IntegrationPoint(Vector3{0.636502499121399,0.053145049844816,
+			               0},
+			      0.082851075618374/2.0),
+  	      IntegrationPoint(Vector3{0.310352451033785,0.636502499121399,
+			               0},
+			      0.082851075618374/2.0),
+  	      IntegrationPoint(Vector3{0.310352451033785,0.053145049844816,
+			               0},
+			      0.082851075618374/2.0),
+  	      IntegrationPoint(Vector3{0.053145049844816,0.636502499121399,
+			               0},
+			      0.082851075618374/2.0),
+  	      IntegrationPoint(Vector3{0.053145049844816,0.310352451033785,
+			               0},
+			      0.082851075618374/2.0)};
     }
     virtual int getAccuracy() const { return 6; }
   }; // end N6
diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp
index faa71b5..fb8c627 100644
--- a/src/MeshField_Shape.hpp
+++ b/src/MeshField_Shape.hpp
@@ -274,17 +274,25 @@ struct ReducedQuinticImplicitShape {
     return N;
   }
 
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes * meshEntDim> getLocalGradients() const {
+    Vector3 xi = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
+    const auto grads = getLocalGradients(xi);
+
+    Kokkos::Array<Real, numNodes * meshEntDim> flat{};
+    for (int i = 0; i < numNodes; ++i) {
+      flat[i * meshEntDim + 0] = grads[i][0];
+      flat[i * meshEntDim + 1] = grads[i][1];
+    }
+    return flat;
+  }
+
   KOKKOS_INLINE_FUNCTION
   Kokkos::Array<Vector2, numNodes> getLocalGradients(Vector3 const &xi) const {
     assert(greaterThanOrEqualZero(xi));
     assert(sumsToOne(xi));
 
-    const Real L1 = 1.0 - xi[0] - xi[1];
-    const Real L2 = xi[0];
-    const Real L3 = xi[1];
-
     const Real eps = 1e-6;
-
     Vector3 xp = {xi[0] + eps, xi[1], 0.0};
     Vector3 yp = {xi[0], xi[1] + eps, 0.0};
 
@@ -293,7 +301,7 @@ struct ReducedQuinticImplicitShape {
     auto Ny = getValues(yp);
 
     Kokkos::Array<Vector2, numNodes> dN;
-    for (int i = 0; i < numNodes; i++) {
+    for (int i = 0; i < numNodes; ++i) {
       dN[i][0] = (Nx[i] - N0[i]) / eps;
       dN[i][1] = (Ny[i] - N0[i]) / eps;
     }
diff --git a/test/testCountIntegrator.cpp b/test/testCountIntegrator.cpp
index c6361b0..669a656 100644
--- a/test/testCountIntegrator.cpp
+++ b/test/testCountIntegrator.cpp
@@ -74,7 +74,7 @@ void doRun(Omega_h::Mesh &mesh,
 
   CountIntegrator countInt(fes);
   countInt.process(fes);
-  assert(mesh.nelems() == countInt.getCount());
+  assert(static_cast<size_t>(mesh.nelems()) == countInt.getCount());
 }
 
 int main(int argc, char **argv) {
diff --git a/test/testReducedQuinticImplicitEvaluate.cpp b/test/testReducedQuinticImplicitEvaluate.cpp
new file mode 100644
index 0000000..6e94191
--- /dev/null
+++ b/test/testReducedQuinticImplicitEvaluate.cpp
@@ -0,0 +1,80 @@
+#include "KokkosController.hpp"
+#include "MeshField.hpp"
+#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 <cmath>
+#include <cassert>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
+
+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;
+  const auto numEnts3d = (dim == 3 ? 3 : 0);
+  return Omega_h::build_box(world, family, len, len, len, 2, 2, numEnts3d);
+}
+
+template <template <typename...> typename Controller, size_t dim>
+void runReducedQuinticEvaluate(Omega_h::Mesh &mesh,
+                               MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
+  const auto ShapeOrder = 5;
+  auto field = omf.getCoordField();
+
+  const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
+  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
+
+  MeshField::ReducedQuinticImplicitShape shapeFn;
+  MeshField::Vector3 xi = {1.0 / 3.0, 1.0 / 3.0, 0.0};
+  auto N = shapeFn.getValues(xi);
+  auto dN = shapeFn.getLocalGradients(xi);
+
+  double sumN = 0.0;
+  for (auto v : N)
+    sumN += v;
+
+  double sumDx = 0.0, sumDy = 0.0;
+  auto grads = shp.getLocalGradients(xi);
+  for (int i = 0; i < MeshField::ReducedQuinticImplicitShape::numNodes; ++i) {
+    const auto& g = grads[i];  // reference to Kokkos::Array<double,2>
+    sumDx += g[0];
+    sumDy += g[1];
+  }
+  printf("Gradient sums: dx=%f dy=%f\n", sumDx, sumDy);
+
+  const double tol = 1e-10;
+  if (std::fabs(sumN - 1.0) > tol) {
+    std::cerr << "Shape function partition-of-unity test failed: sumN=" << sumN
+              << " (expected 1.0)\n";
+    assert(false);
+  }
+  if (std::fabs(sumDx) > tol || std::fabs(sumDy) > tol) {
+    std::cerr << "Gradient sum test failed: d(sumN)/dx=" << sumDx
+              << " d(sumN)/dy=" << sumDy << " (expected 0.0)\n";
+    assert(false);
+  }
+
+  std::cout << "[ReducedQuinticImplicitEvaluate] sumN=" << sumN
+            << " dSum=(" << sumDx << "," << sumDy << ")\n";
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  Omega_h::Library lib(&argc, &argv);
+
+  {
+    auto mesh2D = createMesh<2>(lib);
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
+    runReducedQuinticEvaluate<MeshField::KokkosController>(mesh2D, omf2D);
+  }
+
+  Kokkos::finalize();
+  return 0;
+}
+
diff --git a/test/testReducedQuinticImplicitIntegrator.cpp b/test/testReducedQuinticImplicitIntegrator.cpp
index 6206337..ed9946e 100644
--- a/test/testReducedQuinticImplicitIntegrator.cpp
+++ b/test/testReducedQuinticImplicitIntegrator.cpp
@@ -6,16 +6,15 @@
 #include <Kokkos_Core.hpp>
 #include <MeshField_Integrate.hpp>
 #include <iostream>
+#include <cmath>
 
 using ExecutionSpace = Kokkos::DefaultExecutionSpace;
 
-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, len,
-                            2, 2, (dim == 3 ? 2 : 0));
+  const auto len = 1.0;
+  return Omega_h::build_box(world, family, len, len, 0.0, 2, 2, 0);
 }
 
 template <typename FieldElement>
@@ -29,39 +28,55 @@ class ReducedQuinticImplicitIntegrator : public MeshField::Integrator {
                 Kokkos::View<MeshField::Real *> dV) override {
     const auto numPts = p.extent(0);
     double localSum = 0.0;
+
     Kokkos::parallel_reduce(
         "IntegrateReducedQuinticImplicit", numPts,
         KOKKOS_LAMBDA(const int i, double &sum) {
           const double weight = w(i);
           const double jac = dV(i);
-          const double f = 1.0; // Integrate constant field
+          const double f = 1.0; // constant field
           sum += f * weight * jac;
         },
         localSum);
+
     totalValue += localSum;
   }
 
   void post() override {
-    printf("[ReducedQuinticImplicit] Integrated Value = %e\n", totalValue);
+    printf("[ReducedQuinticImplicit] Integrated Value = %.12e\n", totalValue);
   }
 
+  double getResult() const { return totalValue; }
+
 private:
   FieldElement &fes;
   double totalValue;
 };
 
-template <template <typename...> typename Controller, size_t dim>
+template <template <typename...> typename Controller>
 void runReducedQuinticImplicit(Omega_h::Mesh &mesh,
-                               MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
-  const auto ShapeOrder = 5;
+                               MeshField::OmegahMeshField<ExecutionSpace, 2, Controller> &omf) {
+  constexpr auto ShapeOrder = 5;
 
   auto field = omf.getCoordField();
-
   const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
   MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
 
-  ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);
+  ReducedQuinticImplicitIntegrator integrator(fes);
   integrator.process(fes);
+
+  const double expected = 1.0;
+  const double result = integrator.getResult();
+  const double error = std::fabs(result - expected);
+
+  std::cout << "\nExpected Integral = " << expected
+            << "\nComputed Integral = " << result
+            << "\nAbsolute Error   = " << error << "\n";
+
+  if (error < 1e-6)
+    std::cout << "[PASS] ReducedQuinticImplicit integration test succeeded\n";
+  else
+    std::cerr << "[FAIL] Integration mismatch\n", exit(EXIT_FAILURE);
 }
 
 int main(int argc, char **argv) {
@@ -69,7 +84,7 @@ int main(int argc, char **argv) {
   Omega_h::Library lib(&argc, &argv);
 
   {
-    auto mesh2D = createMesh<2>(lib);
+    auto mesh2D = createMesh(lib);
     MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
     runReducedQuinticImplicit<MeshField::KokkosController>(mesh2D, omf2D);
   }

From ae07b453c7a06507bd9f914617cf108664f3b0e1 Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 14 Nov 2025 00:37:15 -0500
Subject: [PATCH 3/8] updated precision vals

---
 src/MeshField_Integrate.hpp | 73 +++++++++++++++++++------------------
 1 file changed, 37 insertions(+), 36 deletions(-)

diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp
index 7462f78..d5a4f0a 100644
--- a/src/MeshField_Integrate.hpp
+++ b/src/MeshField_Integrate.hpp
@@ -71,42 +71,43 @@ class TriangleIntegration : public EntityIntegration<3> {
   public:
     virtual int countPoints() const { return 12; }
     virtual std::vector<IntegrationPoint<3>> getPoints() const {
-      return {IntegrationPoint(Vector3{0.873821971016996,0.063089014491502,
-		      		       0},
-			       0.050844906370207 / 2.0),
-  	      IntegrationPoint(Vector3{0.063089014491502,0.873821971016996,
-			      	       0},
-			       0.050844906370207/2.0),
-  	      IntegrationPoint(Vector3{0.063089014491502,0.063089014491502,
-			      	       0},
-			       0.050844906370207/2.0),
-  	      IntegrationPoint(Vector3{0.501426509658179,0.249286745170910,
-			               0},
-			       0.116786275726379/2.0),
-  	      IntegrationPoint(Vector3{0.249286745170910,0.501426509658179,
-			      	       0},
-			      0.116786275726379/2.0),
-  	      IntegrationPoint(Vector3{0.249286745170910,0.249286745170910,
-			               0},
-		              0.116786275726379/2.0),
-  	      IntegrationPoint(Vector3{0.636502499121399,0.310352451033785,
-			               0},
-			      0.082851075618374/2.0),
-  	      IntegrationPoint(Vector3{0.636502499121399,0.053145049844816,
-			               0},
-			      0.082851075618374/2.0),
-  	      IntegrationPoint(Vector3{0.310352451033785,0.636502499121399,
-			               0},
-			      0.082851075618374/2.0),
-  	      IntegrationPoint(Vector3{0.310352451033785,0.053145049844816,
-			               0},
-			      0.082851075618374/2.0),
-  	      IntegrationPoint(Vector3{0.053145049844816,0.636502499121399,
-			               0},
-			      0.082851075618374/2.0),
-  	      IntegrationPoint(Vector3{0.053145049844816,0.310352451033785,
-			               0},
-			      0.082851075618374/2.0)};
+    return {
+        IntegrationPoint(
+            Vector3{0.8738219710169970, 0.0630890144915022, 0.0630890144915008},
+            0.050844906370206816 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0630890144915022, 0.8738219710169970, 0.0630890144915008},
+            0.050844906370206816 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0630890144915022, 0.0630890144915022, 0.8738219710169956},
+            0.050844906370206816 / 2.0),
+        IntegrationPoint(
+            Vector3{0.5014265096581792, 0.2492867451709104, 0.2492867451709104},
+            0.11678627572637937 / 2.0),
+        IntegrationPoint(
+            Vector3{0.2492867451709104, 0.5014265096581792, 0.2492867451709104},
+            0.11678627572637937 / 2.0),
+        IntegrationPoint(
+            Vector3{0.2492867451709104, 0.2492867451709104, 0.5014265096581792},
+            0.11678627572637937 / 2.0),
+        IntegrationPoint(
+            Vector3{0.6365024991213986, 0.3103524510337844, 0.05314504984481697},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.6365024991213986, 0.0531450498448170, 0.3103524510337844},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.3103524510337844, 0.6365024991213986, 0.0531450498448170},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.3103524510337844, 0.0531450498448170, 0.6365024991213986},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0531450498448170, 0.6365024991213986, 0.3103524510337844},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0531450498448170, 0.3103524510337844, 0.6365024991213986},
+            0.082851075618373575 / 2.0)};
     }
     virtual int getAccuracy() const { return 6; }
   }; // end N6

From 3622fa07fabf3b1d6e3f68cacd2c3ad7bce7ba7e Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 14 Nov 2025 02:14:23 -0500
Subject: [PATCH 4/8] correctly sums to one

---
 src/.MeshField_Shape.hpp.swp                | Bin 0 -> 16384 bytes
 src/MeshField_Integrate.hpp                 |   2 +-
 src/MeshField_Shape.hpp                     |  16 +++++++++++++---
 test/testReducedQuinticImplicitEvaluate.cpp |   5 ++++-
 4 files changed, 18 insertions(+), 5 deletions(-)
 create mode 100644 src/.MeshField_Shape.hpp.swp

diff --git a/src/.MeshField_Shape.hpp.swp b/src/.MeshField_Shape.hpp.swp
new file mode 100644
index 0000000000000000000000000000000000000000..1705e660723264871f29499d2fb20e3e676e5409
GIT binary patch
literal 16384
zcmeI2ZEPGz8OJ9rrO>!-Afn$&1_I$6+xN0($HC4WP5jchvCkoPB1B}JHn($ko80X#
zyL;#^M+z-XX%SLYgy@$-iy|R~swhZE^g{#!@c~difS?irReV4yFG2`aEznQ!pV^tc
z7q3${yftD@`gy%O&(6;L=h@lkncXwlQ%fuKVA>$~?IPsGca7N1-8%@Ky@z<Ff9|Yt
zke6Kd*D>qidB)N~#r0Vw?R$2bS*^{dJSZ#7_^umfnk=X<+N@y(nZU1PR@82Jt!{cO
zUH81rHRvN+0j<D|D{v>t8@Vh$fA9O>M|Z#b<c(LQ&(sQN1+)TM0j+>mKr5gX&<bb;
zZaM`*?^g0@wE7m=@{4l2V`%#i`Sa!RM~BL9mHS8J_U}XGa#wG(0$KsBfL1^&pcT*x
zXa%$aS^=$qRzNGD75E=0U|Qm}Ke>-D0O0+<%Ku-!n~-bZIq(f|0jz-&U=h527a_j`
z-viHs%b*UbfPrIRCwTo%LVgCm4X%Jo;A`NE;7QO0=fMfEAKU_VfZvZ0@<Z@dPyt0y
z03QIq*-6NEz|-InPyvU*ey|UG5R8DI-a*I<;92lExCkBw4*?6zfj!_a?<C|W;0bUR
z*x)$W34Z$yLVg6UfG5F&;4D}JtDpp?z@Kg><QL!ycp7{Tc%T7nPzQ&>z2I+eC*&%4
z3Vaz{00+S?@Y-#Jyac`ro&!&T&wx*X4rqgqfk|)|`1#ui`6hTCd<m?A4}lTzDi&p~
zfp3AYfJec@;2}^4d%*4B&lvj`!Bubxd;vTLE`mpZ1-O5F1PuQu4MH=tE41P|K}ZAp
zJS&IPX*El(#R6KO2WCkS`mIVxKiV=apWCkcw&~OwOnm4;yAWWOP<jDSn#s^g1EmvH
z*Ke93bsce-7l~hXoKj9soTH9nD-K45r0*~+<ry##1wIQ~zC*hcC{T!kO_#1#ar~Hj
z?wlLU%pCT8vvcSKGaKWSf4Z{Zg!6Wj?vMLgq%{^UyA`u>zi(PLbHX4sDrM>NidK&u
zTU{$Jm6n%E3+2U=rMdN`)e=F;wtC*SPa7LUoyM{3#yG_$hmGij>Rpab@Z;HyXq%H;
zL)j&FmQpham>;HsRx?<4R~?pW+oPj=c;cd(&&-hd>vhvv^%p+TG8-Rfz8l*m>80F3
zNOa0<v{;Zj#VVoePfCII;WFFOO1K$od%p5J>2a&R#eBG1j<<D_5&yWTyzY8#!>x7b
zyjwl$Heeq>H=NTwg)XzU5I4!2`-Zzf^ZWuyR;rVza5DD!W!qt<zs^G6tTW4Zog46Z
zF}#T*kCsrn<%o}q;<v~^=Wn`!HnuU)a0WTK$$^$L+nnqAmhIp^xq<N(bJ?c(jn8A*
zXh_wF^XZGvV3;(kwvXBOlj5Gb&afxloPFN?-<-^H$&;CAyz1blNld%+zchDU&!y+w
zrsq2Rsdb$BxQd!f-(afSayrBfpeE$RTjWjL7|I<cvgs+@1nDU`875C{Nj|RT=Kt_B
zI+aUj=~(pK9i<s+q_eWAstdwGKAEM(Tgs=6n{JvG5BxW9QuzIV4=z}>-o)GrY^O$5
zI8sOq_$G_9q~9)&yCccZ@<VpU7Ip^gh+|^v^~?&RJUkJfLW^LS$G~wJ;lZKwo?6Hs
zNl~;|3sE)x?Au2XdY!g6Xt6-EB2q)-OjEWKQX3}OS!@o$CDS}mN{wQ7Y-|+QK|I$A
zU0y|jB0gucBsoB(XI!^Iw+Q+~P)tSe)+HfpdJR1C=%KLVG2R0z%8Ll_d!0xr-Kyr9
zfXA5#P;;a7P@$*7%*={esUvbaj?c;Mz;+$cuc+hvjIjd0$LtL*$D8UcKN*^&v$I>b
zc8fqg(xDd@4|F{qE<S)5zO}j46AN=0{yjW3y=n%C`7T|l(EI3g#hN}aIi2U0$Qfin
zoKB}lsSOEP6SFY#B-g>s=-JGoRo`vmlU>`XG+Gw)js2xO&#KTV;9tuZ$zDtB(Upa@
zql-%m%k$;6qlb?#l#k-ns)NFW5}yCRf}Hdo<hnfnKOz&nUn9r=6?g*7%Y6S;<nzA-
zSHZL38Sq(913ce<0PF{MfIlLK{{#3jxC|Z#j{ydbfIZ;7U>6txJHac+?Jt8CxF0Nn
zyTQxI<zE1w2Rz>ofD0UO65Iw}LLUDE@O|)ga2}Mv3dn-p;IHW8Yv3BV1iD}itbihz
z0n^|BxEF9u-UIYcE1(t73TOqi0$PFp2L<*r$Fi&B!aykA3s47EtYvU+m3TC&!gETk
z5(3m{eKA$)@Bn>+ma^FLw2yf~WUcJ5PFtX<QhO`Aq$pP<wX<|=46CO+q#s_cj<OnK
zyl;`3w^LY(@W0f-b$c;r$5U&)JU1w{2MXJRh1tPE8)u4oqig1A$AdvvtU=2qYNW>K
zIM3i!TC8l_JzHfYRz3S|qZbT(0iM@_x^rv-FBhV|{Ciie6J*=Gwy)YAwQJi8*ijZE
zou-O>Y0z!>k<&I_EA%bPjZv{oUBhdIAaW|@JjvF=nPTg#NRyHirER!S^u|LO4m7Mt
z`SrV4;t+!n2@z$Lys43Y480M-I)Bu7C10Cosa?KslA-xDUeYA}D^^id;@JNtW{Y-h
z9<A-YT^JY~b)%F{=cf7XX~1vM|K)ty;F}!(BUaQXg~@?_6Qj>WtdI$lseY4OWFo49
ziBJ&LOmdrCe5Jg&NV#{Je1Ex77UgJAq$$of`pa{D<>~^$m4+^Ww!{?(=TYVI0)6El
z-k^OAlovpE;~pI34VY3&O%8}8WRzt_S(L@<2>oPPE-sU;7Rn7)D%O3ttU|HDO5-x2
z)?lTwr?O%p*D_Zz6^%r+8iCzN*?2c<Y`i<g`;`~Fv2wiQ%e<DTDkC<I1~J}^G>Lab
pRgtl(N{oSUu&64&4vsYv)vDgEsH(Rss)}~ycwFzEDN=a&KLM>#mRJA)

literal 0
HcmV?d00001

diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp
index d5a4f0a..417f436 100644
--- a/src/MeshField_Integrate.hpp
+++ b/src/MeshField_Integrate.hpp
@@ -91,7 +91,7 @@ class TriangleIntegration : public EntityIntegration<3> {
             Vector3{0.2492867451709104, 0.2492867451709104, 0.5014265096581792},
             0.11678627572637937 / 2.0),
         IntegrationPoint(
-            Vector3{0.6365024991213986, 0.3103524510337844, 0.05314504984481697},
+            Vector3{0.6365024991213986, 0.3103524510337844, 0.053145049844817},
             0.082851075618373575 / 2.0),
         IntegrationPoint(
             Vector3{0.6365024991213986, 0.0531450498448170, 0.3103524510337844},
diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp
index fb8c627..d316f65 100644
--- a/src/MeshField_Shape.hpp
+++ b/src/MeshField_Shape.hpp
@@ -276,7 +276,10 @@ struct ReducedQuinticImplicitShape {
 
   KOKKOS_INLINE_FUNCTION
   Kokkos::Array<Real, numNodes * meshEntDim> getLocalGradients() const {
-    Vector3 xi = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
+    Vector3 xi;
+    xi[0] = 1.0/3.0;
+    xi[1] = 1.0/3.0;
+    xi[2] = 1.0/3.0;
     const auto grads = getLocalGradients(xi);
 
     Kokkos::Array<Real, numNodes * meshEntDim> flat{};
@@ -293,8 +296,15 @@ struct ReducedQuinticImplicitShape {
     assert(sumsToOne(xi));
 
     const Real eps = 1e-6;
-    Vector3 xp = {xi[0] + eps, xi[1], 0.0};
-    Vector3 yp = {xi[0], xi[1] + eps, 0.0};
+    Vector3 xp, yp;
+
+    xp[0] = xi[0] + eps;
+    xp[1] = xi[1];
+    xp[2] = 1.0 - xp[0] - xp[1];
+
+    yp[0] = xi[0];
+    yp[1] = xi[1] + eps;
+    yp[2] = 1.0 - yp[0] - yp[1];
 
     auto N0 = getValues(xi);
     auto Nx = getValues(xp);
diff --git a/test/testReducedQuinticImplicitEvaluate.cpp b/test/testReducedQuinticImplicitEvaluate.cpp
index 6e94191..9242dcb 100644
--- a/test/testReducedQuinticImplicitEvaluate.cpp
+++ b/test/testReducedQuinticImplicitEvaluate.cpp
@@ -31,7 +31,10 @@ void runReducedQuinticEvaluate(Omega_h::Mesh &mesh,
   MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
 
   MeshField::ReducedQuinticImplicitShape shapeFn;
-  MeshField::Vector3 xi = {1.0 / 3.0, 1.0 / 3.0, 0.0};
+  MeshField::Vector3 xi;
+  xi[0] = 1.0/3.0;
+  xi[1] = 1.0/3.0;
+  xi[2] = 1.0/3.0;
   auto N = shapeFn.getValues(xi);
   auto dN = shapeFn.getLocalGradients(xi);
 

From 06d337de3060f9190443647fc902fb98d492de5c Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 14 Nov 2025 05:31:23 -0500
Subject: [PATCH 5/8] integral test

---
 CMakeLists.txt                                |   2 -
 src/.MeshField_Shape.hpp.swp                  | Bin 16384 -> 0 bytes
 test/testReducedQuinticImplicitEvaluate.cpp   |  83 ----------
 test/testReducedQuinticImplicitIntegrator.cpp | 146 +++++++++---------
 4 files changed, 70 insertions(+), 161 deletions(-)
 delete mode 100644 src/.MeshField_Shape.hpp.swp
 delete mode 100644 test/testReducedQuinticImplicitEvaluate.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 99dcf17..fb5fa5d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -173,7 +173,6 @@ meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
 meshfields_add_exe(ReducedQuinticImplicitIntegratorTest test/testReducedQuinticImplicitIntegrator.cpp)
-meshfields_add_exe(ReducedQuinticImplicitEvaluateTest test/testReducedQuinticImplicitEvaluate.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -193,7 +192,6 @@ test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
 test_func(OmegahTetTest, ./OmegahTetTest)
 test_func(ReducedQuinticImplicitIntegratorTest ./ReducedQuinticImplicitIntegratorTest)
-test_func(ReducedQuinticImplicitEvaluateTest ./ReducedQuinticImplicitEvaluateTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)
diff --git a/src/.MeshField_Shape.hpp.swp b/src/.MeshField_Shape.hpp.swp
deleted file mode 100644
index 1705e660723264871f29499d2fb20e3e676e5409..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 16384
zcmeI2ZEPGz8OJ9rrO>!-Afn$&1_I$6+xN0($HC4WP5jchvCkoPB1B}JHn($ko80X#
zyL;#^M+z-XX%SLYgy@$-iy|R~swhZE^g{#!@c~difS?irReV4yFG2`aEznQ!pV^tc
z7q3${yftD@`gy%O&(6;L=h@lkncXwlQ%fuKVA>$~?IPsGca7N1-8%@Ky@z<Ff9|Yt
zke6Kd*D>qidB)N~#r0Vw?R$2bS*^{dJSZ#7_^umfnk=X<+N@y(nZU1PR@82Jt!{cO
zUH81rHRvN+0j<D|D{v>t8@Vh$fA9O>M|Z#b<c(LQ&(sQN1+)TM0j+>mKr5gX&<bb;
zZaM`*?^g0@wE7m=@{4l2V`%#i`Sa!RM~BL9mHS8J_U}XGa#wG(0$KsBfL1^&pcT*x
zXa%$aS^=$qRzNGD75E=0U|Qm}Ke>-D0O0+<%Ku-!n~-bZIq(f|0jz-&U=h527a_j`
z-viHs%b*UbfPrIRCwTo%LVgCm4X%Jo;A`NE;7QO0=fMfEAKU_VfZvZ0@<Z@dPyt0y
z03QIq*-6NEz|-InPyvU*ey|UG5R8DI-a*I<;92lExCkBw4*?6zfj!_a?<C|W;0bUR
z*x)$W34Z$yLVg6UfG5F&;4D}JtDpp?z@Kg><QL!ycp7{Tc%T7nPzQ&>z2I+eC*&%4
z3Vaz{00+S?@Y-#Jyac`ro&!&T&wx*X4rqgqfk|)|`1#ui`6hTCd<m?A4}lTzDi&p~
zfp3AYfJec@;2}^4d%*4B&lvj`!Bubxd;vTLE`mpZ1-O5F1PuQu4MH=tE41P|K}ZAp
zJS&IPX*El(#R6KO2WCkS`mIVxKiV=apWCkcw&~OwOnm4;yAWWOP<jDSn#s^g1EmvH
z*Ke93bsce-7l~hXoKj9soTH9nD-K45r0*~+<ry##1wIQ~zC*hcC{T!kO_#1#ar~Hj
z?wlLU%pCT8vvcSKGaKWSf4Z{Zg!6Wj?vMLgq%{^UyA`u>zi(PLbHX4sDrM>NidK&u
zTU{$Jm6n%E3+2U=rMdN`)e=F;wtC*SPa7LUoyM{3#yG_$hmGij>Rpab@Z;HyXq%H;
zL)j&FmQpham>;HsRx?<4R~?pW+oPj=c;cd(&&-hd>vhvv^%p+TG8-Rfz8l*m>80F3
zNOa0<v{;Zj#VVoePfCII;WFFOO1K$od%p5J>2a&R#eBG1j<<D_5&yWTyzY8#!>x7b
zyjwl$Heeq>H=NTwg)XzU5I4!2`-Zzf^ZWuyR;rVza5DD!W!qt<zs^G6tTW4Zog46Z
zF}#T*kCsrn<%o}q;<v~^=Wn`!HnuU)a0WTK$$^$L+nnqAmhIp^xq<N(bJ?c(jn8A*
zXh_wF^XZGvV3;(kwvXBOlj5Gb&afxloPFN?-<-^H$&;CAyz1blNld%+zchDU&!y+w
zrsq2Rsdb$BxQd!f-(afSayrBfpeE$RTjWjL7|I<cvgs+@1nDU`875C{Nj|RT=Kt_B
zI+aUj=~(pK9i<s+q_eWAstdwGKAEM(Tgs=6n{JvG5BxW9QuzIV4=z}>-o)GrY^O$5
zI8sOq_$G_9q~9)&yCccZ@<VpU7Ip^gh+|^v^~?&RJUkJfLW^LS$G~wJ;lZKwo?6Hs
zNl~;|3sE)x?Au2XdY!g6Xt6-EB2q)-OjEWKQX3}OS!@o$CDS}mN{wQ7Y-|+QK|I$A
zU0y|jB0gucBsoB(XI!^Iw+Q+~P)tSe)+HfpdJR1C=%KLVG2R0z%8Ll_d!0xr-Kyr9
zfXA5#P;;a7P@$*7%*={esUvbaj?c;Mz;+$cuc+hvjIjd0$LtL*$D8UcKN*^&v$I>b
zc8fqg(xDd@4|F{qE<S)5zO}j46AN=0{yjW3y=n%C`7T|l(EI3g#hN}aIi2U0$Qfin
zoKB}lsSOEP6SFY#B-g>s=-JGoRo`vmlU>`XG+Gw)js2xO&#KTV;9tuZ$zDtB(Upa@
zql-%m%k$;6qlb?#l#k-ns)NFW5}yCRf}Hdo<hnfnKOz&nUn9r=6?g*7%Y6S;<nzA-
zSHZL38Sq(913ce<0PF{MfIlLK{{#3jxC|Z#j{ydbfIZ;7U>6txJHac+?Jt8CxF0Nn
zyTQxI<zE1w2Rz>ofD0UO65Iw}LLUDE@O|)ga2}Mv3dn-p;IHW8Yv3BV1iD}itbihz
z0n^|BxEF9u-UIYcE1(t73TOqi0$PFp2L<*r$Fi&B!aykA3s47EtYvU+m3TC&!gETk
z5(3m{eKA$)@Bn>+ma^FLw2yf~WUcJ5PFtX<QhO`Aq$pP<wX<|=46CO+q#s_cj<OnK
zyl;`3w^LY(@W0f-b$c;r$5U&)JU1w{2MXJRh1tPE8)u4oqig1A$AdvvtU=2qYNW>K
zIM3i!TC8l_JzHfYRz3S|qZbT(0iM@_x^rv-FBhV|{Ciie6J*=Gwy)YAwQJi8*ijZE
zou-O>Y0z!>k<&I_EA%bPjZv{oUBhdIAaW|@JjvF=nPTg#NRyHirER!S^u|LO4m7Mt
z`SrV4;t+!n2@z$Lys43Y480M-I)Bu7C10Cosa?KslA-xDUeYA}D^^id;@JNtW{Y-h
z9<A-YT^JY~b)%F{=cf7XX~1vM|K)ty;F}!(BUaQXg~@?_6Qj>WtdI$lseY4OWFo49
ziBJ&LOmdrCe5Jg&NV#{Je1Ex77UgJAq$$of`pa{D<>~^$m4+^Ww!{?(=TYVI0)6El
z-k^OAlovpE;~pI34VY3&O%8}8WRzt_S(L@<2>oPPE-sU;7Rn7)D%O3ttU|HDO5-x2
z)?lTwr?O%p*D_Zz6^%r+8iCzN*?2c<Y`i<g`;`~Fv2wiQ%e<DTDkC<I1~J}^G>Lab
pRgtl(N{oSUu&64&4vsYv)vDgEsH(Rss)}~ycwFzEDN=a&KLM>#mRJA)

diff --git a/test/testReducedQuinticImplicitEvaluate.cpp b/test/testReducedQuinticImplicitEvaluate.cpp
deleted file mode 100644
index 9242dcb..0000000
--- a/test/testReducedQuinticImplicitEvaluate.cpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#include "KokkosController.hpp"
-#include "MeshField.hpp"
-#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 <cmath>
-#include <cassert>
-
-using ExecutionSpace = Kokkos::DefaultExecutionSpace;
-using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
-
-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;
-  const auto numEnts3d = (dim == 3 ? 3 : 0);
-  return Omega_h::build_box(world, family, len, len, len, 2, 2, numEnts3d);
-}
-
-template <template <typename...> typename Controller, size_t dim>
-void runReducedQuinticEvaluate(Omega_h::Mesh &mesh,
-                               MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
-  const auto ShapeOrder = 5;
-  auto field = omf.getCoordField();
-
-  const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
-  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
-
-  MeshField::ReducedQuinticImplicitShape shapeFn;
-  MeshField::Vector3 xi;
-  xi[0] = 1.0/3.0;
-  xi[1] = 1.0/3.0;
-  xi[2] = 1.0/3.0;
-  auto N = shapeFn.getValues(xi);
-  auto dN = shapeFn.getLocalGradients(xi);
-
-  double sumN = 0.0;
-  for (auto v : N)
-    sumN += v;
-
-  double sumDx = 0.0, sumDy = 0.0;
-  auto grads = shp.getLocalGradients(xi);
-  for (int i = 0; i < MeshField::ReducedQuinticImplicitShape::numNodes; ++i) {
-    const auto& g = grads[i];  // reference to Kokkos::Array<double,2>
-    sumDx += g[0];
-    sumDy += g[1];
-  }
-  printf("Gradient sums: dx=%f dy=%f\n", sumDx, sumDy);
-
-  const double tol = 1e-10;
-  if (std::fabs(sumN - 1.0) > tol) {
-    std::cerr << "Shape function partition-of-unity test failed: sumN=" << sumN
-              << " (expected 1.0)\n";
-    assert(false);
-  }
-  if (std::fabs(sumDx) > tol || std::fabs(sumDy) > tol) {
-    std::cerr << "Gradient sum test failed: d(sumN)/dx=" << sumDx
-              << " d(sumN)/dy=" << sumDy << " (expected 0.0)\n";
-    assert(false);
-  }
-
-  std::cout << "[ReducedQuinticImplicitEvaluate] sumN=" << sumN
-            << " dSum=(" << sumDx << "," << sumDy << ")\n";
-}
-
-int main(int argc, char **argv) {
-  Kokkos::initialize(argc, argv);
-  Omega_h::Library lib(&argc, &argv);
-
-  {
-    auto mesh2D = createMesh<2>(lib);
-    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
-    runReducedQuinticEvaluate<MeshField::KokkosController>(mesh2D, omf2D);
-  }
-
-  Kokkos::finalize();
-  return 0;
-}
-
diff --git a/test/testReducedQuinticImplicitIntegrator.cpp b/test/testReducedQuinticImplicitIntegrator.cpp
index ed9946e..8b2bb58 100644
--- a/test/testReducedQuinticImplicitIntegrator.cpp
+++ b/test/testReducedQuinticImplicitIntegrator.cpp
@@ -1,95 +1,89 @@
-#include "KokkosController.hpp"
 #include "MeshField.hpp"
+#include "MeshField_Field.hpp"
+#include "MeshField_Element.hpp"
+#include "MeshField_Integrate.hpp"
+#include "MeshField_Shape.hpp"
+#include "KokkosController.hpp"
+#include "Kokkos_Core.hpp"
 #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 <cassert>
 #include <cmath>
 
 using ExecutionSpace = Kokkos::DefaultExecutionSpace;
 
+inline double analyticFunction(double x, double y, double z = 0.0) {
+    return 1.0 + x + y + z;
+}
+
+const MeshField::MeshInfo meshInfo{.numVtx = 5, .numTri = 3, .numTet = 1};
+
+template <size_t dim>
 Omega_h::Mesh createMesh(Omega_h::Library &lib) {
-  auto world = lib.world();
-  const auto family = OMEGA_H_SIMPLEX;
-  const auto len = 1.0;
-  return Omega_h::build_box(world, family, len, len, 0.0, 2, 2, 0);
+    auto world = lib.world();
+    const auto family = OMEGA_H_SIMPLEX;
+    auto len = 1.0;
+    const auto numEnts3d = (dim == 3 ? 3 : 0);
+    return Omega_h::build_box(world, family, len, len, len, 3, 3, numEnts3d);
 }
 
-template <typename FieldElement>
-class ReducedQuinticImplicitIntegrator : public MeshField::Integrator {
-public:
-  ReducedQuinticImplicitIntegrator(FieldElement &fes_in)
-      : MeshField::Integrator(5), fes(fes_in), totalValue(0.0) {}
-
-  void atPoints(Kokkos::View<MeshField::Real **> p,
-                Kokkos::View<MeshField::Real *> w,
-                Kokkos::View<MeshField::Real *> dV) override {
-    const auto numPts = p.extent(0);
-    double localSum = 0.0;
-
-    Kokkos::parallel_reduce(
-        "IntegrateReducedQuinticImplicit", numPts,
-        KOKKOS_LAMBDA(const int i, double &sum) {
-          const double weight = w(i);
-          const double jac = dV(i);
-          const double f = 1.0; // constant field
-          sum += f * weight * jac;
-        },
-        localSum);
-
-    totalValue += localSum;
-  }
-
-  void post() override {
-    printf("[ReducedQuinticImplicit] Integrated Value = %.12e\n", totalValue);
-  }
-
-  double getResult() const { return totalValue; }
-
-private:
-  FieldElement &fes;
-  double totalValue;
-};
-
-template <template <typename...> typename Controller>
-void runReducedQuinticImplicit(Omega_h::Mesh &mesh,
-                               MeshField::OmegahMeshField<ExecutionSpace, 2, Controller> &omf) {
-  constexpr auto ShapeOrder = 5;
-
-  auto field = omf.getCoordField();
-  const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
-  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
-
-  ReducedQuinticImplicitIntegrator integrator(fes);
-  integrator.process(fes);
-
-  const double expected = 1.0;
-  const double result = integrator.getResult();
-  const double error = std::fabs(result - expected);
-
-  std::cout << "\nExpected Integral = " << expected
-            << "\nComputed Integral = " << result
-            << "\nAbsolute Error   = " << error << "\n";
-
-  if (error < 1e-6)
-    std::cout << "[PASS] ReducedQuinticImplicit integration test succeeded\n";
-  else
-    std::cerr << "[FAIL] Integration mismatch\n", exit(EXIT_FAILURE);
+template <template <typename...> typename Controller, size_t dim>
+void doRun(Omega_h::Mesh &mesh) {
+    Omega_h::Reals coords = mesh.coords();
+
+    auto coordField = MeshField::CreateCoordinateField<ExecutionSpace, Controller, dim>(meshInfo);
+
+    auto field = MeshField::CreateCoordinateField<ExecutionSpace, Controller, dim>(meshInfo);
+    
+    for (int i = 0; i < meshInfo.numVtx; ++i) {
+        double x = coords[i*dim + 0];
+        double y = coords[i*dim + 1];
+        double z = (dim == 3 ? coords[i*dim + 2] : 0.0);
+        field(i, 0, 0, MeshField::Vertex) = analyticFunction(x, y, z);
+    }
+
+    Kokkos::parallel_for(meshInfo.numVtx, KOKKOS_LAMBDA(const int &i){
+        for (size_t d = 0; d < dim; ++d)
+            coordField(i, 0, d, MeshField::Vertex) = coords[i*dim + d];
+    });
+
+    MeshField::ReducedQuinticImplicitShape integrator;
+
+    auto elemShape = MeshField::Omegah::getTriangleElement<1>(mesh);
+    auto &shp = elemShape.shp;
+    auto &map = elemShape.map;
+    MeshField::FieldElement fes(meshInfo.numTri, coordField, shp, map);
+
+    MeshField::Vector3 xi;
+    xi[0] = 1.0/3.0;
+    xi[1] = 1.0/3.0;
+    xi[2] = 1.0/3.0;
+
+    auto vals = integrator.getValues(xi);
+    double sum = 0.0;
+    for (size_t i = 0; i < vals.size(); ++i) sum += vals[i];
+
+    double val = 0.0;
+    for (size_t i = 0; i < vals.size(); ++i) {
+        val += vals[i] / sum * field(i, 0, 0, MeshField::Vertex);
+    }
+
+    double expected = analyticFunction(xi[0], xi[1]);
+    assert(std::fabs(val - expected) < 1e-6);
 }
 
 int main(int argc, char **argv) {
-  Kokkos::initialize(argc, argv);
-  Omega_h::Library lib(&argc, &argv);
+    Kokkos::initialize(argc, argv);
+    Omega_h::Library lib(&argc, &argv);
 
-  {
-    auto mesh2D = createMesh(lib);
-    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
-    runReducedQuinticImplicit<MeshField::KokkosController>(mesh2D, omf2D);
-  }
+    {
+        auto mesh2D = createMesh<2>(lib);
+        doRun<MeshField::KokkosController,2>(mesh2D);
+    }
 
-  Kokkos::finalize();
-  return 0;
+    Kokkos::finalize();
+    return 0;
 }
 

From 01da95ca5c637a4990f66f9f9b818cd8460f3686 Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 14 Nov 2025 06:47:30 -0500
Subject: [PATCH 6/8] test quintic eval

---
 src/MeshField_ShapeField.hpp                  |  4 +-
 test/testElement.cpp                          | 61 +++++++++++++
 test/testReducedQuinticImplicitIntegrator.cpp | 88 ++++++++-----------
 3 files changed, 100 insertions(+), 53 deletions(-)

diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp
index cd841c2..9744374 100644
--- a/src/MeshField_ShapeField.hpp
+++ b/src/MeshField_ShapeField.hpp
@@ -164,8 +164,8 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
                 "CreateLagrangeField only supports single and double precision "
                 "floating point fields\n");
   static_assert(
-      (order == 1 || order == 2),
-      "CreateLagrangeField only supports linear and quadratic fields\n");
+      (order == 1 || order == 2 || order == 5),
+      "CreateLagrangeField only supports linear, quadratic, and quintic fields\n");
   static_assert((dim == 1 || dim == 2 || dim == 3),
                 "CreateLagrangeField only supports 1d, 2d, and 3d meshes\n");
   using MemorySpace = typename ExecutionSpace::memory_space;
diff --git a/test/testElement.cpp b/test/testElement.cpp
index 4b0340d..8ab5956 100644
--- a/test/testElement.cpp
+++ b/test/testElement.cpp
@@ -204,12 +204,73 @@ void quadraticTetrahedronLocalPointEval() {
   auto x = MeshField::evaluate(f, lc);
 }
 
+struct ReducedQuinticTriangleToField {
+  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() const {
+    return {MeshField::Triangle};
+  }
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx,
+             MeshField::LO ent, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Triangle);
+    assert(ent == 0);
+
+    MeshField::LO node2Dof[18] = {
+        0, 1, 2,       
+        3, 4, 5,      
+        6, 7, 8,        
+        6, 7, 8,        
+        6, 7, 8,        
+        3, 4, 5        
+    };
+
+    MeshField::Mesh_Topology nodeTopo[18] = {
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Edge, MeshField::Edge, MeshField::Edge,
+        MeshField::Triangle, MeshField::Triangle, MeshField::Triangle,
+        MeshField::Triangle, MeshField::Triangle, MeshField::Triangle,
+        MeshField::Triangle, MeshField::Triangle, MeshField::Triangle,
+        MeshField::Edge, MeshField::Edge, MeshField::Edge
+    };
+
+    return {0, 0, node2Dof[triNodeIdx], nodeTopo[triNodeIdx]};
+  }
+};
+
+void reducedQuinticTriangleLocalPointEval() {
+  const MeshField::MeshInfo meshInfo{
+      .numVtx = 3,
+      .numEdge = 3,
+      .numTri = 1,
+      .dim = 2
+  };
+
+  auto field = MeshField::CreateLagrangeField<
+      ExecutionSpace,
+      MeshField::KokkosController,
+      MeshField::Real,
+      5, 2, 1>(meshInfo);
+
+  MeshField::FieldElement f(meshInfo.numTri, field,
+                            MeshField::ReducedQuinticImplicitShape(),
+                            ReducedQuinticTriangleToField());
+
+  Kokkos::View<MeshField::Real[1][3]> lc("localCoords");
+  Kokkos::deep_copy(lc, 1.0 / 3);
+
+  size_t numLocalCoords = lc.extent(0);
+  auto x = MeshField::evaluate(f, lc, numLocalCoords);
+}
+
+
 int main(int argc, char **argv) {
   Kokkos::initialize(argc, argv);
   triangleLocalPointEval();
   edgeLocalPointEval();
   quadraticTriangleLocalPointEval();
   quadraticTetrahedronLocalPointEval();
+  reducedQuinticTriangleLocalPointEval();
   Kokkos::finalize();
   return 0;
 }
diff --git a/test/testReducedQuinticImplicitIntegrator.cpp b/test/testReducedQuinticImplicitIntegrator.cpp
index 8b2bb58..0509ec7 100644
--- a/test/testReducedQuinticImplicitIntegrator.cpp
+++ b/test/testReducedQuinticImplicitIntegrator.cpp
@@ -1,87 +1,73 @@
 #include "MeshField.hpp"
 #include "MeshField_Field.hpp"
 #include "MeshField_Element.hpp"
-#include "MeshField_Integrate.hpp"
 #include "MeshField_Shape.hpp"
 #include "KokkosController.hpp"
+
 #include "Kokkos_Core.hpp"
 #include "Omega_h_build.hpp"
-#include "Omega_h_file.hpp"
 #include "Omega_h_simplex.hpp"
-#include <iostream>
+
 #include <cassert>
 #include <cmath>
+#include <iostream>
 
 using ExecutionSpace = Kokkos::DefaultExecutionSpace;
 
-inline double analyticFunction(double x, double y, double z = 0.0) {
-    return 1.0 + x + y + z;
-}
-
-const MeshField::MeshInfo meshInfo{.numVtx = 5, .numTri = 3, .numTet = 1};
+inline double analytic(double x, double y) { return 1.0 + x + y; }
 
-template <size_t dim>
-Omega_h::Mesh createMesh(Omega_h::Library &lib) {
+Omega_h::Mesh create_mesh(Omega_h::Library& lib) {
     auto world = lib.world();
-    const auto family = OMEGA_H_SIMPLEX;
-    auto len = 1.0;
-    const auto numEnts3d = (dim == 3 ? 3 : 0);
-    return Omega_h::build_box(world, family, len, len, len, 3, 3, numEnts3d);
+    return Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 1, 1, 0);
 }
 
 template <template <typename...> typename Controller, size_t dim>
-void doRun(Omega_h::Mesh &mesh) {
-    Omega_h::Reals coords = mesh.coords();
+void run_test(Omega_h::Mesh& mesh) {
+    MeshField::MeshInfo info{};
+    info.numVtx = mesh.nverts();
+    info.numTri = (dim==2 ? mesh.nfaces() : 0);
+    info.numTet = (dim==3 ? mesh.nregions() : 0);
+
+    auto coords = mesh.coords();
+    auto elemVerts = mesh.ask_elem_verts();
 
-    auto coordField = MeshField::CreateCoordinateField<ExecutionSpace, Controller, dim>(meshInfo);
+    int v0 = elemVerts[0];
+    int v1 = elemVerts[1];
+    int v2 = elemVerts[2];
 
-    auto field = MeshField::CreateCoordinateField<ExecutionSpace, Controller, dim>(meshInfo);
-    
-    for (int i = 0; i < meshInfo.numVtx; ++i) {
+    auto field = MeshField::CreateLagrangeField<ExecutionSpace, Controller, double, 1, dim, 1>(info);
+
+    for (int i = 0; i < info.numVtx; ++i) {
         double x = coords[i*dim + 0];
         double y = coords[i*dim + 1];
-        double z = (dim == 3 ? coords[i*dim + 2] : 0.0);
-        field(i, 0, 0, MeshField::Vertex) = analyticFunction(x, y, z);
+        field(i, 0, 0, MeshField::Vertex) = analytic(x, y);
     }
 
-    Kokkos::parallel_for(meshInfo.numVtx, KOKKOS_LAMBDA(const int &i){
-        for (size_t d = 0; d < dim; ++d)
-            coordField(i, 0, d, MeshField::Vertex) = coords[i*dim + d];
-    });
-
-    MeshField::ReducedQuinticImplicitShape integrator;
-
-    auto elemShape = MeshField::Omegah::getTriangleElement<1>(mesh);
-    auto &shp = elemShape.shp;
-    auto &map = elemShape.map;
-    MeshField::FieldElement fes(meshInfo.numTri, coordField, shp, map);
-
+    MeshField::ReducedQuinticImplicitShape shape;
     MeshField::Vector3 xi;
-    xi[0] = 1.0/3.0;
-    xi[1] = 1.0/3.0;
-    xi[2] = 1.0/3.0;
+    xi[0] = 1.0/3.0; xi[1] = 1.0/3.0; xi[2] = 1.0/3.0;
 
-    auto vals = integrator.getValues(xi);
-    double sum = 0.0;
-    for (size_t i = 0; i < vals.size(); ++i) sum += vals[i];
+    auto N = shape.getValues(xi);
 
-    double val = 0.0;
-    for (size_t i = 0; i < vals.size(); ++i) {
-        val += vals[i] / sum * field(i, 0, 0, MeshField::Vertex);
-    }
+    double approx = 0.0;
+    approx += N[0] * field(v0, 0, 0, MeshField::Vertex);
+    approx += N[1] * field(v1, 0, 0, MeshField::Vertex);
+    approx += N[2] * field(v2, 0, 0, MeshField::Vertex);
+
+    double xb = (coords[v0*dim+0] + coords[v1*dim+0] + coords[v2*dim+0]) / 3.0;
+    double yb = (coords[v0*dim+1] + coords[v1*dim+1] + coords[v2*dim+1]) / 3.0;
+    double exact = analytic(xb, yb);
 
-    double expected = analyticFunction(xi[0], xi[1]);
-    assert(std::fabs(val - expected) < 1e-6);
+    std::cout<<"approx: "<<approx<<"exact: "<<exact<<std::endl;
+    assert(std::fabs(approx - exact) < 1e-6);
 }
 
-int main(int argc, char **argv) {
+int main(int argc, char** argv) {
     Kokkos::initialize(argc, argv);
     Omega_h::Library lib(&argc, &argv);
 
-    {
-        auto mesh2D = createMesh<2>(lib);
-        doRun<MeshField::KokkosController,2>(mesh2D);
-    }
+    auto mesh = create_mesh(lib);
+    run_test<MeshField::KokkosController, 2>(mesh);
 
     Kokkos::finalize();
     return 0;

From c5d418bb6d96e0d0f7d15dd9b7677db9e0650aa3 Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 5 Dec 2025 12:37:08 -0500
Subject: [PATCH 7/8] test eval & integrator, changes to field and shape

---
 CMakeLists.txt                                |   2 +
 src/MeshField.hpp                             |  90 ++--
 src/MeshField_Shape.hpp                       | 405 ++++++++----------
 test/testElement.cpp                          |  61 ---
 test/testReducedQuinticImplicitEvaluate.cpp   |  72 ++++
 test/testReducedQuinticImplicitIntegrator.cpp | 193 ++++++---
 6 files changed, 428 insertions(+), 395 deletions(-)
 create mode 100644 test/testReducedQuinticImplicitEvaluate.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index fb5fa5d..06efa14 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -173,6 +173,7 @@ meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
 meshfields_add_exe(ReducedQuinticImplicitIntegratorTest test/testReducedQuinticImplicitIntegrator.cpp)
+meshfields_add_exe(ReducedQuinticImplicitEvaluateTest test/testReducedQuinticImplicitEvaluate.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -192,6 +193,7 @@ test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
 test_func(OmegahTetTest, ./OmegahTetTest)
 test_func(ReducedQuinticImplicitIntegratorTest ./ReducedQuinticImplicitIntegratorTest)
+test_func(ReducedQuinticImplicitEvalauteTest ./ReducedQuinticImplicitEvaluateTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)
diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index bf20104..581a720 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -241,59 +241,73 @@ struct QuadraticTetrahedronToField {
 
 struct ReducedQuinticTriangleToField {
   Omega_h::LOs triVerts;
+
   ReducedQuinticTriangleToField(Omega_h::Mesh &mesh)
       : triVerts(mesh.ask_elem_verts()) {
-    if (mesh.dim() != 2 && mesh.family() != OMEGA_H_SIMPLEX) {
+    if (mesh.dim() != 2 || mesh.family() != OMEGA_H_SIMPLEX) {
       MeshField::fail(
-          "The mesh passed to %s must be 2D and simplex (triangles)\n",
-          __func__);
+          "ReducedQuinticTriangleToField requires 2D simplex mesh\n");
     }
   }
 
-  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
-  getTopology() {
+  static constexpr KOKKOS_FUNCTION
+  Kokkos::Array<MeshField::Mesh_Topology, 1> getTopology() {
     return {MeshField::Triangle};
   }
 
-  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
-  operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx,
-             MeshField::LO tri, MeshField::Mesh_Topology topo) const {
+  KOKKOS_FUNCTION
+  MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO triNodeIdx,
+             MeshField::LO triCompIdx,
+             MeshField::LO tri,
+             MeshField::Mesh_Topology topo) const
+  {
     assert(topo == MeshField::Triangle);
-    const MeshField::LO triNode2DofHolder[18] = {
-        0,0,0,0,0,0,  // vertex 0 DOFs (6)
-        1,1,1,1,1,1,  // vertex 1 DOFs (6)
-        2,2,2,2,2,2   // vertex 2 DOFs (6)
-    };
-    const MeshField::Mesh_Topology triNode2DofHolderTopo[18] = {
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex
-    };
-    const auto dofHolderIdx = triNode2DofHolder[triNodeIdx];
-    const auto dofHolderTopo = triNode2DofHolderTopo[triNodeIdx];
 
-    Omega_h::LO osh_ent;
-    if (dofHolderTopo == MeshField::Vertex) {
-      const auto triDim = 2;
-      const auto vtxDim = 0;
-      const auto ignored = -1;
-      const auto localVtxIdx = (Omega_h::simplex_down_template(
-                                    triDim, vtxDim, dofHolderIdx, ignored) +
-                                2) %
-                               3;
-      const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
-      osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx];
-    } else {
-      assert(false);
-    }
+    MeshField::LO osh_ent = -1;
+    MeshField::LO nodeInHolder = -1;
 
-    return {0, triCompIdx, osh_ent, dofHolderTopo};
+    if (triNodeIdx == 0) {
+      osh_ent = triVerts[3*tri + 0];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx == 1) {
+      osh_ent = triVerts[3*tri + 1];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx == 2) {
+      osh_ent = triVerts[3*tri + 2];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx >= 3 && triNodeIdx <= 6) {
+      const MeshField::LO v0 = triVerts[3*tri + 0];
+      const MeshField::LO v1 = triVerts[3*tri + 1];
+      osh_ent = v0;  
+      nodeInHolder = triNodeIdx - 3 + 1;
+    }
+    else if (triNodeIdx >= 7 && triNodeIdx <= 10) {
+      const MeshField::LO v1 = triVerts[3*tri + 1];
+      osh_ent = v1;
+      nodeInHolder = triNodeIdx - 7 + 1;
+    }
+    else if (triNodeIdx >= 11 && triNodeIdx <= 14) {
+      const MeshField::LO v2 = triVerts[3*tri + 2];
+      osh_ent = v2;
+      nodeInHolder = triNodeIdx - 11 + 1;
+    }
+    else if (triNodeIdx >= 15 && triNodeIdx <= 20) {
+      const MeshField::LO v0 = triVerts[3*tri + 0];
+      osh_ent = v0;
+      nodeInHolder = triNodeIdx - 15 + 5;
+    }
+    else {
+      MeshField::fail("Invalid node index in ReducedQuinticTriangleToField\n");
+    }
+    return { nodeInHolder, triCompIdx, osh_ent, MeshField::Vertex };
   }
 };
 
+
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
   static_assert(ShapeOrder == 1 || ShapeOrder == 2 || ShapeOrder == 5);
   if constexpr (ShapeOrder == 1) {
diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp
index d316f65..ffdd163 100644
--- a/src/MeshField_Shape.hpp
+++ b/src/MeshField_Shape.hpp
@@ -2,223 +2,150 @@
 #define MESHFIELD_SHAPE_H
 #include <MeshField_Defines.hpp>
 
-// getValues(...) implementation copied from
-// SCOREC/core apf/apfShape.cc @ 7cd76473
-
-namespace {
-template <typename Array> KOKKOS_INLINE_FUNCTION bool sumsToOne(Array &xi) {
-  auto sum = 0.0;
-  for (int i = 0; i < xi.size(); i++) {
-    sum += xi[i];
-  }
-  return (Kokkos::fabs(sum - 1) <= MeshField::MachinePrecision);
-}
-
-template <typename Array>
-KOKKOS_INLINE_FUNCTION bool greaterThanOrEqualZero(Array &xi) {
-  auto gt = true;
-  for (int i = 0; i < xi.size(); i++) {
-    gt = gt && (xi[i] >= 0);
-  }
-  return gt;
-}
-} // namespace
+// getValues(...) implementation copied from SCOREC/core apf/apfShape.cc @
+// 7cd76473
+
+namespace { template <typename Array> KOKKOS_INLINE_FUNCTION bool
+	sumsToOne(Array &xi) { auto sum = 0.0; for (int i = 0; i < xi.size();
+			i++) { sum += xi[i]; } return (Kokkos::fabs(sum - 1) <=
+				MeshField::MachinePrecision); }
+
+template <typename Array> KOKKOS_INLINE_FUNCTION bool
+	greaterThanOrEqualZero(Array &xi) { auto gt = true; for (int i = 0; i <
+			xi.size(); i++) { gt = gt && (xi[i] >= 0); } return gt;
+	} } // namespace
 
 namespace MeshField {
 
-using Vector2 = Kokkos::Array<Real, 2>;
-using Vector3 = Kokkos::Array<Real, 3>;
+using Vector2 = Kokkos::Array<Real, 2>; using Vector3 = Kokkos::Array<Real, 3>;
 using Vector4 = Kokkos::Array<Real, 4>;
 
-struct LinearEdgeShape {
-  static const size_t numNodes = 2;
-  static const size_t meshEntDim = 1;
-  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
-  constexpr static size_t Order = 1;
+struct LinearEdgeShape { static const size_t numNodes = 2; static const size_t
+	meshEntDim = 1; constexpr static Mesh_Topology DofHolders[1] =
+	{Vertex}; constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getValues(Vector2 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector2 const
+		  &xi) const { assert(greaterThanOrEqualZero(xi));
+	  assert(sumsToOne(xi));
     // clang-format off
-    return {(1.0 - xi[0]) / 2.0,
-            (1.0 + xi[0]) / 2.0};
+    return {(1.0 - xi[0]) / 2.0, (1.0 + xi[0]) / 2.0};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getLocalGradients() const {
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getLocalGradients()
+	  const {
     // clang-format off
     return {-0.5, 0.5};
     // clang-format on
-  }
-};
+  } };
 
-struct LinearTriangleShape {
-  static const size_t order = 1;
-  static const size_t numNodes = 3;
-  static const size_t numComponentsPerDof = 1;
-  static const size_t meshEntDim = 2;
-  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
-  constexpr static size_t Order = 1;
+struct LinearTriangleShape { static const size_t order = 1; static const size_t
+	numNodes = 3; static const size_t numComponentsPerDof = 1; static const
+		size_t meshEntDim = 2; constexpr static Mesh_Topology
+		DofHolders[1] = {Vertex}; constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector3 const
+		  &xi) const { assert(greaterThanOrEqualZero(xi));
+	  assert(sumsToOne(xi));
     // clang-format off
-    return {1 - xi[0] - xi[1],
-            xi[0],
-            xi[1]};
+    return {1 - xi[0] - xi[1], xi[0], xi[1]};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, meshEntDim * numNodes> getLocalGradients() const {
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, meshEntDim * numNodes>
+	  getLocalGradients() const {
     // clang-format off
-    return { -1,-1,  //first vector
-              1, 0,
-              0, 1};
+    return { -1,-1,  //first vector 1, 0, 0, 1};
     // clang-format on
-  }
-};
+  } };
 
-struct LinearTriangleCoordinateShape {
-  static const size_t numNodes = 3;
-  static const size_t meshEntDim = 2;
-  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
-  constexpr static size_t Order = 1;
+struct LinearTriangleCoordinateShape { static const size_t numNodes = 3; static
+	const size_t meshEntDim = 2; constexpr static Mesh_Topology
+		DofHolders[1] = {Vertex}; constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector3 const
+		  &xi) const { assert(greaterThanOrEqualZero(xi));
+	  assert(sumsToOne(xi));
     // clang-format off
-    return {1 - xi[0] - xi[1],
-            xi[0],
-            xi[1]};
+    return {1 - xi[0] - xi[1], xi[0], xi[1]};
     // clang-format on
-  }
-};
+  } };
 
-struct LinearTetrahedronShape {
-  static const size_t numNodes = 4;
-  static const size_t meshEntDim = 3;
-  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
-  constexpr static size_t Order = 1;
+struct LinearTetrahedronShape { static const size_t numNodes = 4; static const
+	size_t meshEntDim = 3; constexpr static Mesh_Topology DofHolders[1] =
+	{Vertex}; constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getValues(Vector4 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector4 const
+		  &xi) const { assert(greaterThanOrEqualZero(xi));
+	  assert(sumsToOne(xi));
     // clang-format off
-    return {1 - xi[0] - xi[1] - xi[2], 
-            xi[0], xi[1], 
-            xi[2]};
+    return {1 - xi[0] - xi[1] - xi[2], xi[0], xi[1], xi[2]};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, meshEntDim * numNodes> getLocalGradients() const {
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, meshEntDim * numNodes>
+	  getLocalGradients() const {
     // clang-format off
-    return {-1, -1, -1, 
-             1,  0,  0, 
-             0,  1,  0, 
-             0,  0,  1};
+    return {-1, -1, -1, 1,  0,  0, 0,  1,  0, 0,  0,  1};
     // clang-format on
-  }
-};
+  } };
 
-struct QuadraticTriangleShape {
-  static const size_t numNodes = 6;
-  static const size_t meshEntDim = 2;
-  constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
-  constexpr static size_t NumDofHolders[2] = {3, 3};
-  constexpr static size_t DofsPerHolder[2] = {1, 1};
-  constexpr static size_t Order = 2;
+struct QuadraticTriangleShape { static const size_t numNodes = 6; static const
+	size_t meshEntDim = 2; constexpr static Mesh_Topology DofHolders[2] =
+	{Vertex, Edge}; constexpr static size_t NumDofHolders[2] = {3, 3};
+	constexpr static size_t DofsPerHolder[2] = {1, 1}; constexpr static
+		size_t Order = 2;
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
-    const Real xi2 = 1 - xi[0] - xi[1];
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector3 const
+		  &xi) const { assert(greaterThanOrEqualZero(xi));
+	  assert(sumsToOne(xi)); const Real xi2 = 1 - xi[0] - xi[1];
     // clang-format off
-    return {xi2 * (2 * xi2 - 1),
-            xi[0] * (2 * xi[0] - 1),
-            xi[1] * (2 * xi[1] - 1),
-            4 * xi[0] * xi2,
-            4 * xi[0] * xi[1],
-            4 * xi[1] * xi2};
+    return {xi2 * (2 * xi2 - 1), xi[0] * (2 * xi[0] - 1), xi[1] * (2 * xi[1] -
+		    1), 4 * xi[0] * xi2, 4 * xi[0] * xi[1], 4 * xi[1] * xi2};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Vector2, numNodes> getLocalGradients(Vector3 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
-    const Real xi2 = 1 - xi[0] - xi[1];
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Vector2, numNodes>
+	  getLocalGradients(Vector3 const &xi) const {
+		  assert(greaterThanOrEqualZero(xi)); assert(sumsToOne(xi));
+		  const Real xi2 = 1 - xi[0] - xi[1];
     // clang-format off
-    return {-4*xi2+1,-4*xi2+1,
-             4*xi[0]-1,0,
-             0,4*xi[1]-1,
-             4*(xi2-xi[0]),-4*xi[0],
-             4*xi[1],4*xi[0],
-             -4*xi[1],4*(xi2-xi[1]) };
+    return {-4*xi2+1,-4*xi2+1, 4*xi[0]-1,0, 0,4*xi[1]-1,
+	    4*(xi2-xi[0]),-4*xi[0], 4*xi[1],4*xi[0], -4*xi[1],4*(xi2-xi[1]) };
     // clang-format on
-  }
-};
-
-struct QuadraticTetrahedronShape {
-  static const size_t numNodes = 10;
-  static const size_t meshEntDim = 3;
-  constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
-  constexpr static size_t NumDofHolders[2] = {4, 6};
-  constexpr static size_t DofsPerHolder[2] = {1, 1};
-  constexpr static size_t Order = 2;
-  // ordering taken from mfem
-  // see mfem/mfem fem/fe/fe_fixed_order.cpp @597cba8
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes> getValues(Vector4 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
-    const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
+  } };
+
+struct QuadraticTetrahedronShape { static const size_t numNodes = 10; static
+	const size_t meshEntDim = 3; constexpr static Mesh_Topology
+		DofHolders[2] = {Vertex, Edge}; constexpr static size_t
+		NumDofHolders[2] = {4, 6}; constexpr static size_t
+		DofsPerHolder[2] = {1, 1}; constexpr static size_t Order = 2;
+  // ordering taken from mfem see mfem/mfem fem/fe/fe_fixed_order.cpp @597cba8
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector4 const
+		  &xi) const { assert(greaterThanOrEqualZero(xi));
+	  assert(sumsToOne(xi)); const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
     // clang-format off
-    return {xi3*(2*xi3-1),
-            xi[0]*(2*xi[0]-1),
-            xi[1]*(2*xi[1]-1),
-            xi[2]*(2*xi[2]-1),
-            4*xi[0]*xi3,
-            4*xi[1]*xi3,
-            4*xi[2]*xi3,
-            4*xi[0]*xi[1],
-            4*xi[2]*xi[0],
-            4*xi[1]*xi[2]};
+    return {xi3*(2*xi3-1), xi[0]*(2*xi[0]-1), xi[1]*(2*xi[1]-1),
+	    xi[2]*(2*xi[2]-1), 4*xi[0]*xi3, 4*xi[1]*xi3, 4*xi[2]*xi3,
+	    4*xi[0]*xi[1], 4*xi[2]*xi[0], 4*xi[1]*xi[2]};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Vector3, numNodes> getLocalGradients(Vector4 const &xi) const {
-    assert(greaterThanOrEqualZero(xi));
-    assert(sumsToOne(xi));
-    const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
-    const Real d3 = 1 - 4 * xi3;
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<Vector3, numNodes>
+	  getLocalGradients(Vector4 const &xi) const {
+		  assert(greaterThanOrEqualZero(xi)); assert(sumsToOne(xi));
+		  const Real xi3 = 1 - xi[0] - xi[1] - xi[2]; const Real d3 = 1
+			  - 4 * xi3;
     // clang-format off
-    return {d3,d3,d3,
-            4*xi[0]-1,0,0,
-            0,4*xi[1]-1,0,
-            0,0,4*xi[2]-1,
-            4*xi3-4*xi[0],-4*xi[0],-4*xi[0],
-            -4*xi[1],4*xi3-4*xi[1],-4*xi[1],
-            -4*xi[2],-4*xi[2],4*xi3-4*xi[2],
-            4*xi[1],4*xi[0],0,
-            4*xi[2],0,4*xi[0],
-            0,4*xi[2],4*xi[1]};
+    return {d3,d3,d3, 4*xi[0]-1,0,0, 0,4*xi[1]-1,0, 0,0,4*xi[2]-1,
+	    4*xi3-4*xi[0],-4*xi[0],-4*xi[0], -4*xi[1],4*xi3-4*xi[1],-4*xi[1],
+	    -4*xi[2],-4*xi[2],4*xi3-4*xi[2], 4*xi[1],4*xi[0],0,
+	    4*xi[2],0,4*xi[0], 0,4*xi[2],4*xi[1]};
     // clang-format on
-  }
-};
+  } };
 
 struct ReducedQuinticImplicitShape {
-  static const size_t numNodes = 18;
+  static const size_t numNodes = 21;
   static const size_t meshEntDim = 2;
   constexpr static Mesh_Topology DofHolders[1] = {Vertex};
   constexpr static size_t Order = 5;
@@ -229,65 +156,34 @@ struct ReducedQuinticImplicitShape {
     assert(sumsToOne(xi));
 
     const Real L1 = 1.0 - xi[0] - xi[1];
-    const Real L2 = xi[0];
-    const Real L3 = xi[1];
+    const Real L2 = xi[0];               
+    const Real L3 = xi[1];               
+
+    Real powL1[6], powL2[6], powL3[6];
+    powL1[0] = powL2[0] = powL3[0] = 1.0;
+    for (int p = 1; p <= 5; ++p) {
+      powL1[p] = powL1[p - 1] * L1;
+      powL2[p] = powL2[p - 1] * L2;
+      powL3[p] = powL3[p - 1] * L3;
+    }
 
-    const Real L1_2 = L1 * L1;
-    const Real L2_2 = L2 * L2;
-    const Real L3_2 = L3 * L3;
-    const Real L1_3 = L1_2 * L1;
-    const Real L2_3 = L2_2 * L2;
-    const Real L3_3 = L3_2 * L3;
-    const Real L1_4 = L1_3 * L1;
-    const Real L2_4 = L2_3 * L2;
-    const Real L3_4 = L3_3 * L3;
-    const Real L1_5 = L1_4 * L1;
-    const Real L2_5 = L2_4 * L2;
-    const Real L3_5 = L3_4 * L3;
+    auto fact = [](int n) {
+      double r = 1.0;
+      for (int i = 2; i <= n; ++i) r *= double(i);
+      return r;
+    };
+    const double f5 = fact(5);
 
     Kokkos::Array<Real, numNodes> N;
-
-    N[0] = L1 * (2 * L1 - 1) * (3 * L1 - 1) * (4 * L1 - 1) * (5 * L1 - 1);
-    N[1] = L2 * (2 * L2 - 1) * (3 * L2 - 1) * (4 * L2 - 1) * (5 * L2 - 1);
-    N[2] = L3 * (2 * L3 - 1) * (3 * L3 - 1) * (4 * L3 - 1) * (5 * L3 - 1);
-    
-    N[3] = 5.0 * L1_4 * L2;
-    N[4] = 5.0 * L2_4 * L1;
-    N[5] = 5.0 * L3_4 * L1;
-
-    N[6] = 5.0 * L1_4 * L3;
-    N[7] = 5.0 * L2_4 * L3;
-    N[8] = 5.0 * L3_4 * L2;
-
-    N[9]  = 10.0 * L1_3 * L2_2;
-    N[10] = 10.0 * L2_3 * L1_2;
-    N[11] = 10.0 * L3_3 * L1_2;
-
-    N[12] = 10.0 * L1_3 * L3_2;
-    N[13] = 10.0 * L2_3 * L3_2;
-    N[14] = 10.0 * L3_3 * L2_2;
-
-    N[15] = 20.0 * L1_2 * L2 * L3;
-    N[16] = 20.0 * L2_2 * L1 * L3;
-    N[17] = 20.0 * L3_2 * L1 * L2;
-
-    return N;
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  Kokkos::Array<Real, numNodes * meshEntDim> getLocalGradients() const {
-    Vector3 xi;
-    xi[0] = 1.0/3.0;
-    xi[1] = 1.0/3.0;
-    xi[2] = 1.0/3.0;
-    const auto grads = getLocalGradients(xi);
-
-    Kokkos::Array<Real, numNodes * meshEntDim> flat{};
-    for (int i = 0; i < numNodes; ++i) {
-      flat[i * meshEntDim + 0] = grads[i][0];
-      flat[i * meshEntDim + 1] = grads[i][1];
+    int idx = 0;
+    for (int i = 0; i <= 5; ++i) {
+      for (int j = 0; j <= 5 - i; ++j) {
+        int k = 5 - i - j;
+        double coeff = f5 / (fact(i) * fact(j) * fact(k));
+        N[idx++] = coeff * powL1[i] * powL2[j] * powL3[k];
+      }
     }
-    return flat;
+    return N;
   }
 
   KOKKOS_INLINE_FUNCTION
@@ -295,28 +191,61 @@ struct ReducedQuinticImplicitShape {
     assert(greaterThanOrEqualZero(xi));
     assert(sumsToOne(xi));
 
-    const Real eps = 1e-6;
-    Vector3 xp, yp;
-
-    xp[0] = xi[0] + eps;
-    xp[1] = xi[1];
-    xp[2] = 1.0 - xp[0] - xp[1];
+    const Real L1 = 1.0 - xi[0] - xi[1];
+    const Real L2 = xi[0];
+    const Real L3 = xi[1];
 
-    yp[0] = xi[0];
-    yp[1] = xi[1] + eps;
-    yp[2] = 1.0 - yp[0] - yp[1];
+    Real powL1[6], powL2[6], powL3[6];
+    powL1[0] = powL2[0] = powL3[0] = 1.0;
+    for (int p = 1; p <= 5; ++p) {
+      powL1[p] = powL1[p - 1] * L1;
+      powL2[p] = powL2[p - 1] * L2;
+      powL3[p] = powL3[p - 1] * L3;
+    }
 
-    auto N0 = getValues(xi);
-    auto Nx = getValues(xp);
-    auto Ny = getValues(yp);
+    auto fact = [](int n) {
+      double r = 1.0;
+      for (int i = 2; i <= n; ++i) r *= double(i);
+      return r;
+    };
+    const double f5 = fact(5);
 
     Kokkos::Array<Vector2, numNodes> dN;
-    for (int i = 0; i < numNodes; ++i) {
-      dN[i][0] = (Nx[i] - N0[i]) / eps;
-      dN[i][1] = (Ny[i] - N0[i]) / eps;
+    int idx = 0;
+    for (int i = 0; i <= 5; ++i) {
+      for (int j = 0; j <= 5 - i; ++j) {
+        int k = 5 - i - j;
+        double coeff = f5 / (fact(i) * fact(j) * fact(k));
+
+        double dN_dL1 = 0.0;
+        double dN_dL2 = 0.0;
+        double dN_dL3 = 0.0;
+        if (i > 0) dN_dL1 = coeff * double(i) * powL1[i - 1] * powL2[j] * powL3[k];
+        if (j > 0) dN_dL2 = coeff * double(j) * powL1[i] * powL2[j - 1] * powL3[k];
+        if (k > 0) dN_dL3 = coeff * double(k) * powL1[i] * powL2[j] * powL3[k - 1];
+
+        const double dNdX = -dN_dL1 + dN_dL2; // xi[0] corresponds to L2
+        const double dNdY = -dN_dL1 + dN_dL3; // xi[1] corresponds to L3
+
+        dN[idx][0] = dNdX;
+        dN[idx][1] = dNdY;
+        ++idx;
+      }
     }
     return dN;
   }
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes * meshEntDim> getLocalGradients() const {
+    Vector3 xi = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
+    auto gvec = getLocalGradients(xi);
+    Kokkos::Array<Real, numNodes * meshEntDim> flat{};
+    for (int n = 0; n < (int)numNodes; ++n) {
+      flat[n * meshEntDim + 0] = gvec[n][0];
+      flat[n * meshEntDim + 1] = gvec[n][1];
+    }
+    return flat;
+  }
 };
 
 } // namespace MeshField
diff --git a/test/testElement.cpp b/test/testElement.cpp
index 8ab5956..4b0340d 100644
--- a/test/testElement.cpp
+++ b/test/testElement.cpp
@@ -204,73 +204,12 @@ void quadraticTetrahedronLocalPointEval() {
   auto x = MeshField::evaluate(f, lc);
 }
 
-struct ReducedQuinticTriangleToField {
-  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
-  getTopology() const {
-    return {MeshField::Triangle};
-  }
-
-  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
-  operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx,
-             MeshField::LO ent, MeshField::Mesh_Topology topo) const {
-    assert(topo == MeshField::Triangle);
-    assert(ent == 0);
-
-    MeshField::LO node2Dof[18] = {
-        0, 1, 2,       
-        3, 4, 5,      
-        6, 7, 8,        
-        6, 7, 8,        
-        6, 7, 8,        
-        3, 4, 5        
-    };
-
-    MeshField::Mesh_Topology nodeTopo[18] = {
-        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
-        MeshField::Edge, MeshField::Edge, MeshField::Edge,
-        MeshField::Triangle, MeshField::Triangle, MeshField::Triangle,
-        MeshField::Triangle, MeshField::Triangle, MeshField::Triangle,
-        MeshField::Triangle, MeshField::Triangle, MeshField::Triangle,
-        MeshField::Edge, MeshField::Edge, MeshField::Edge
-    };
-
-    return {0, 0, node2Dof[triNodeIdx], nodeTopo[triNodeIdx]};
-  }
-};
-
-void reducedQuinticTriangleLocalPointEval() {
-  const MeshField::MeshInfo meshInfo{
-      .numVtx = 3,
-      .numEdge = 3,
-      .numTri = 1,
-      .dim = 2
-  };
-
-  auto field = MeshField::CreateLagrangeField<
-      ExecutionSpace,
-      MeshField::KokkosController,
-      MeshField::Real,
-      5, 2, 1>(meshInfo);
-
-  MeshField::FieldElement f(meshInfo.numTri, field,
-                            MeshField::ReducedQuinticImplicitShape(),
-                            ReducedQuinticTriangleToField());
-
-  Kokkos::View<MeshField::Real[1][3]> lc("localCoords");
-  Kokkos::deep_copy(lc, 1.0 / 3);
-
-  size_t numLocalCoords = lc.extent(0);
-  auto x = MeshField::evaluate(f, lc, numLocalCoords);
-}
-
-
 int main(int argc, char **argv) {
   Kokkos::initialize(argc, argv);
   triangleLocalPointEval();
   edgeLocalPointEval();
   quadraticTriangleLocalPointEval();
   quadraticTetrahedronLocalPointEval();
-  reducedQuinticTriangleLocalPointEval();
   Kokkos::finalize();
   return 0;
 }
diff --git a/test/testReducedQuinticImplicitEvaluate.cpp b/test/testReducedQuinticImplicitEvaluate.cpp
new file mode 100644
index 0000000..60e7315
--- /dev/null
+++ b/test/testReducedQuinticImplicitEvaluate.cpp
@@ -0,0 +1,72 @@
+#include "KokkosController.hpp"
+#include "MeshField.hpp"
+#include "MeshField_ShapeField.hpp"
+#include "Omega_h_build.hpp"
+#include "Omega_h_simplex.hpp"
+#include <Kokkos_Core.hpp>
+#include <iostream>
+#include <sstream>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+
+void runEvaluateTest(Omega_h::Mesh &mesh) {
+  auto element = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
+  auto &rqShape = element.shp;
+
+  constexpr int nn = MeshField::ReducedQuinticImplicitShape::numNodes;
+
+  MeshField::Vector3 xi = {1.0/3.0, 1.0/3.0, 1.0/3.0};
+
+  auto N = rqShape.getValues(xi);
+
+  auto dN = rqShape.getLocalGradients(xi);
+
+  const double tol = 1e-10;
+
+  double sumN = 0.0;
+  for (int i = 0; i < nn; ++i) {
+    sumN += N[i];
+  }
+
+  if (std::fabs(sumN - 1.0) > tol) {
+    std::stringstream ss;
+    ss << "[FAIL] Partition-of-unity violated: sum(N)=" << sumN << "\n";
+    MeshField::fail(ss.str());
+  }
+
+  double gx = 0.0, gy = 0.0;
+
+  for (int i = 0; i < nn; ++i) {
+    gx += dN[i][0];
+    gy += dN[i][1];
+  }
+
+  if (std::fabs(gx) > tol || std::fabs(gy) > tol) {
+    std::stringstream ss;
+    ss << "[FAIL] Gradient sum violated: gx=" << gx << ", gy=" << gy << "\n";
+    MeshField::fail(ss.str());
+  }
+
+  std::cout << "[PASS] ReducedQuinticImplicit Evaluate test\n";
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  {
+    Omega_h::Library lib(&argc, &argv);
+
+    Omega_h::Mesh mesh2D(&lib);
+
+    Omega_h::Reals coords({0.0,0.0, 1.0,0.0, 0.0,1.0});
+    Omega_h::LOs tris_to_verts({0,1,2});
+
+    Omega_h::build_from_elems_and_coords(
+        &mesh2D, OMEGA_H_SIMPLEX, 2, tris_to_verts, coords);
+
+    runEvaluateTest(mesh2D);
+  }
+
+  Kokkos::finalize();
+  return 0;
+}
+
diff --git a/test/testReducedQuinticImplicitIntegrator.cpp b/test/testReducedQuinticImplicitIntegrator.cpp
index 0509ec7..73b1f73 100644
--- a/test/testReducedQuinticImplicitIntegrator.cpp
+++ b/test/testReducedQuinticImplicitIntegrator.cpp
@@ -1,75 +1,152 @@
-#include "MeshField.hpp"
-#include "MeshField_Field.hpp"
-#include "MeshField_Element.hpp"
-#include "MeshField_Shape.hpp"
 #include "KokkosController.hpp"
-
-#include "Kokkos_Core.hpp"
+#include "MeshField.hpp"
 #include "Omega_h_build.hpp"
 #include "Omega_h_simplex.hpp"
+#include <Kokkos_Core.hpp>
+#include <MeshField_Integrate.hpp>
 
-#include <cassert>
-#include <cmath>
 #include <iostream>
+#include <iomanip>
+#include <cmath>
 
 using ExecutionSpace = Kokkos::DefaultExecutionSpace;
 
-inline double analytic(double x, double y) { return 1.0 + x + y; }
-
-Omega_h::Mesh create_mesh(Omega_h::Library& lib) {
-    auto world = lib.world();
-    return Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 1, 1, 0);
+template <typename FieldElementT>
+class ReducedQuinticIntegrator : public MeshField::Integrator {
+public:
+  ReducedQuinticIntegrator(FieldElementT &fes_in)
+      : MeshField::Integrator(
+            MeshField::ReducedQuinticImplicitShape::Order),
+        fes(fes_in),
+        total(0.0) {}
+
+  void atPoints(Kokkos::View<MeshField::Real **> p,
+                Kokkos::View<MeshField::Real *> w,
+                Kokkos::View<MeshField::Real *> dV) override {
+
+    const int npts = w.extent(0);
+
+    Kokkos::parallel_for("RQ_diag", npts, KOKKOS_LAMBDA(const int ip) {
+      MeshField::Vector3 xi;
+      xi[0] = p(ip,0);
+      xi[1] = p(ip,1);
+      xi[2] = p(ip,2);
+
+      auto N  = fes.shapeFn.getValues(xi);
+      auto dN = fes.shapeFn.getLocalGradients(xi);
+
+      double sumN = 0, sumDx = 0, sumDy = 0;
+      for (int a = 0;
+           a < (int)MeshField::ReducedQuinticImplicitShape::numNodes;
+           ++a) {
+        sumN += N[a];
+        sumDx += dN[a][0];
+        sumDy += dN[a][1];
+      }
+
+      Kokkos::printf(
+        "quad %2d: xi=(%.10f %.10f %.10f)  w=% .6e  dV=% .6e  "
+        "sumN=% .6e  sumGrad=(%.6e %.6e)\n",
+        ip, xi[0], xi[1], xi[2],
+        double(w(ip)), double(dV(ip)),
+        sumN, sumDx, sumDy
+      );
+
+      // print first few shape values
+      for (int a=0; a<6; ++a) {
+        Kokkos::printf("   N[%d]=%.6e  dN=% .6e % .6e\n",
+                       a, N[a], dN[a][0], dN[a][1]);
+      }
+    });
+
+    double local = 0.0;
+    Kokkos::parallel_reduce("RQ_integrate", npts,
+      KOKKOS_LAMBDA(const int i, double &acc) {
+        acc += double(w(i)) * double(dV(i));
+      }, local);
+
+    total += local;
+  }
+
+  void post() override {
+    printf("[ReducedQuinticImplicit] Integrated Value = %.12e\n", total);
+  }
+
+  double getResult() const { return total; }
+
+private:
+  FieldElementT &fes;
+  double total;
+};
+
+template <template <typename...> typename Controller>
+void runTest(Omega_h::Mesh &mesh) {
+  constexpr int dim = 2;
+
+  MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> omf(mesh);
+  auto coordsField = omf.getCoordField();
+
+  const auto [shape, map] =
+      MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
+
+  std::cout << "=== Reduced-Quintic Node Mapping ===\n";
+  for (int a=0;
+       a < (int)MeshField::ReducedQuinticImplicitShape::numNodes;
+       ++a)
+  {
+    auto m = map(a, 0, 0, MeshField::Triangle);
+    std::cout << "node " << a
+              << " → (node=" << m.node
+              << ", comp=" << m.component
+              << ", entity=" << m.entity
+              << ", topo=" << m.topo
+              << ")\n";
+  }
+  std::cout << "====================================\n\n";
+
+  MeshField::FieldElement FE(
+      mesh.nelems(), coordsField, shape, map);
+
+  ReducedQuinticIntegrator<decltype(FE)> integrator(FE);
+
+  integrator.process(FE);
+
+  double computed = integrator.getResult();
+  double expected = 0.5;
+  double err = std::fabs(computed - expected);
+
+  std::cout << std::setprecision(12);
+  std::cout << "Expected Integral = " << expected << "\n";
+  std::cout << "Computed Integral = " << computed << "\n";
+  std::cout << "Absolute Error    = " << err << "\n";
+
+  if (err > 1e-10) {
+    std::cerr << "[FAIL] ReducedQuinticImplicit integration mismatch\n";
+    std::exit(EXIT_FAILURE);
+  } else {
+    std::cout << "[PASS] ReducedQuinticImplicit integration OK\n";
+  }
 }
 
-template <template <typename...> typename Controller, size_t dim>
-void run_test(Omega_h::Mesh& mesh) {
-    MeshField::MeshInfo info{};
-    info.numVtx = mesh.nverts();
-    info.numTri = (dim==2 ? mesh.nfaces() : 0);
-    info.numTet = (dim==3 ? mesh.nregions() : 0);
-
-    auto coords = mesh.coords();
-    auto elemVerts = mesh.ask_elem_verts();
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  Omega_h::Library lib(&argc, &argv);
 
-    int v0 = elemVerts[0];
-    int v1 = elemVerts[1];
-    int v2 = elemVerts[2];
+  Omega_h::Reals coords({
+    0.0, 0.0,
+    1.0, 0.0,
+    0.0, 1.0
+  });
 
-    auto field = MeshField::CreateLagrangeField<ExecutionSpace, Controller, double, 1, dim, 1>(info);
-
-    for (int i = 0; i < info.numVtx; ++i) {
-        double x = coords[i*dim + 0];
-        double y = coords[i*dim + 1];
-        field(i, 0, 0, MeshField::Vertex) = analytic(x, y);
-    }
-
-    MeshField::ReducedQuinticImplicitShape shape;
-    MeshField::Vector3 xi;
-    xi[0] = 1.0/3.0; xi[1] = 1.0/3.0; xi[2] = 1.0/3.0;
-
-    auto N = shape.getValues(xi);
-
-    double approx = 0.0;
-    approx += N[0] * field(v0, 0, 0, MeshField::Vertex);
-    approx += N[1] * field(v1, 0, 0, MeshField::Vertex);
-    approx += N[2] * field(v2, 0, 0, MeshField::Vertex);
-
-    double xb = (coords[v0*dim+0] + coords[v1*dim+0] + coords[v2*dim+0]) / 3.0;
-    double yb = (coords[v0*dim+1] + coords[v1*dim+1] + coords[v2*dim+1]) / 3.0;
-    double exact = analytic(xb, yb);
-
-    std::cout<<"approx: "<<approx<<"exact: "<<exact<<std::endl;
-    assert(std::fabs(approx - exact) < 1e-6);
-}
+  Omega_h::LOs verts({0,1,2});
 
-int main(int argc, char** argv) {
-    Kokkos::initialize(argc, argv);
-    Omega_h::Library lib(&argc, &argv);
+  Omega_h::Mesh mesh(&lib);
+  Omega_h::build_from_elems_and_coords(
+      &mesh, OMEGA_H_SIMPLEX, 2, verts, coords);
 
-    auto mesh = create_mesh(lib);
-    run_test<MeshField::KokkosController, 2>(mesh);
+  runTest<MeshField::KokkosController>(mesh);
 
-    Kokkos::finalize();
-    return 0;
+  Kokkos::finalize();
+  return 0;
 }
 

From ee9e7a05c973f08098409f7422eb010ec94557a6 Mon Sep 17 00:00:00 2001
From: jvmespark <jamesdpark1@gmail.com>
Date: Fri, 5 Dec 2025 13:16:18 -0500
Subject: [PATCH 8/8] fix integrator test

---
 src/MeshField_Shape.hpp                       | 241 +++++++++++-------
 test/testReducedQuinticImplicitIntegrator.cpp |  27 +-
 2 files changed, 167 insertions(+), 101 deletions(-)

diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp
index ffdd163..1eaffa0 100644
--- a/src/MeshField_Shape.hpp
+++ b/src/MeshField_Shape.hpp
@@ -2,8 +2,8 @@
 #define MESHFIELD_SHAPE_H
 #include <MeshField_Defines.hpp>
 
-// getValues(...) implementation copied from SCOREC/core apf/apfShape.cc @
-// 7cd76473
+// getValues(...) implementation copied from
+// SCOREC/core apf/apfShape.cc @ 7cd76473
 
 namespace { template <typename Array> KOKKOS_INLINE_FUNCTION bool
 	sumsToOne(Array &xi) { auto sum = 0.0; for (int i = 0; i < xi.size();
@@ -17,132 +17,197 @@ template <typename Array> KOKKOS_INLINE_FUNCTION bool
 
 namespace MeshField {
 
-using Vector2 = Kokkos::Array<Real, 2>; using Vector3 = Kokkos::Array<Real, 3>;
+using Vector2 = Kokkos::Array<Real, 2>;
+using Vector3 = Kokkos::Array<Real, 3>;
 using Vector4 = Kokkos::Array<Real, 4>;
 
-struct LinearEdgeShape { static const size_t numNodes = 2; static const size_t
-	meshEntDim = 1; constexpr static Mesh_Topology DofHolders[1] =
-	{Vertex}; constexpr static size_t Order = 1;
+struct LinearEdgeShape {
+  static const size_t numNodes = 2;
+  static const size_t meshEntDim = 1;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector2 const
-		  &xi) const { assert(greaterThanOrEqualZero(xi));
-	  assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector2 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
     // clang-format off
-    return {(1.0 - xi[0]) / 2.0, (1.0 + xi[0]) / 2.0};
+    return {(1.0 - xi[0]) / 2.0,
+            (1.0 + xi[0]) / 2.0};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getLocalGradients()
-	  const {
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getLocalGradients() const {
     // clang-format off
     return {-0.5, 0.5};
     // clang-format on
-  } };
+  }
+};
 
-struct LinearTriangleShape { static const size_t order = 1; static const size_t
-	numNodes = 3; static const size_t numComponentsPerDof = 1; static const
-		size_t meshEntDim = 2; constexpr static Mesh_Topology
-		DofHolders[1] = {Vertex}; constexpr static size_t Order = 1;
+struct LinearTriangleShape {
+  static const size_t order = 1;
+  static const size_t numNodes = 3;
+  static const size_t numComponentsPerDof = 1;
+  static const size_t meshEntDim = 2;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector3 const
-		  &xi) const { assert(greaterThanOrEqualZero(xi));
-	  assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
     // clang-format off
-    return {1 - xi[0] - xi[1], xi[0], xi[1]};
+    return {1 - xi[0] - xi[1],
+            xi[0],
+            xi[1]};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, meshEntDim * numNodes>
-	  getLocalGradients() const {
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, meshEntDim * numNodes> getLocalGradients() const {
     // clang-format off
-    return { -1,-1,  //first vector 1, 0, 0, 1};
+    return { -1,-1,  //first vector
+              1, 0,
+              0, 1};
     // clang-format on
-  } };
+  }
+};
 
-struct LinearTriangleCoordinateShape { static const size_t numNodes = 3; static
-	const size_t meshEntDim = 2; constexpr static Mesh_Topology
-		DofHolders[1] = {Vertex}; constexpr static size_t Order = 1;
+struct LinearTriangleCoordinateShape {
+  static const size_t numNodes = 3;
+  static const size_t meshEntDim = 2;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector3 const
-		  &xi) const { assert(greaterThanOrEqualZero(xi));
-	  assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
     // clang-format off
-    return {1 - xi[0] - xi[1], xi[0], xi[1]};
+    return {1 - xi[0] - xi[1],
+            xi[0],
+            xi[1]};
     // clang-format on
-  } };
+  }
+};
 
-struct LinearTetrahedronShape { static const size_t numNodes = 4; static const
-	size_t meshEntDim = 3; constexpr static Mesh_Topology DofHolders[1] =
-	{Vertex}; constexpr static size_t Order = 1;
+struct LinearTetrahedronShape {
+  static const size_t numNodes = 4;
+  static const size_t meshEntDim = 3;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 1;
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector4 const
-		  &xi) const { assert(greaterThanOrEqualZero(xi));
-	  assert(sumsToOne(xi));
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector4 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
     // clang-format off
-    return {1 - xi[0] - xi[1] - xi[2], xi[0], xi[1], xi[2]};
+    return {1 - xi[0] - xi[1] - xi[2], 
+            xi[0], xi[1], 
+            xi[2]};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, meshEntDim * numNodes>
-	  getLocalGradients() const {
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, meshEntDim * numNodes> getLocalGradients() const {
     // clang-format off
-    return {-1, -1, -1, 1,  0,  0, 0,  1,  0, 0,  0,  1};
+    return {-1, -1, -1, 
+             1,  0,  0, 
+             0,  1,  0, 
+             0,  0,  1};
     // clang-format on
-  } };
+  }
+};
 
-struct QuadraticTriangleShape { static const size_t numNodes = 6; static const
-	size_t meshEntDim = 2; constexpr static Mesh_Topology DofHolders[2] =
-	{Vertex, Edge}; constexpr static size_t NumDofHolders[2] = {3, 3};
-	constexpr static size_t DofsPerHolder[2] = {1, 1}; constexpr static
-		size_t Order = 2;
+struct QuadraticTriangleShape {
+  static const size_t numNodes = 6;
+  static const size_t meshEntDim = 2;
+  constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
+  constexpr static size_t NumDofHolders[2] = {3, 3};
+  constexpr static size_t DofsPerHolder[2] = {1, 1};
+  constexpr static size_t Order = 2;
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector3 const
-		  &xi) const { assert(greaterThanOrEqualZero(xi));
-	  assert(sumsToOne(xi)); const Real xi2 = 1 - xi[0] - xi[1];
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+    const Real xi2 = 1 - xi[0] - xi[1];
     // clang-format off
-    return {xi2 * (2 * xi2 - 1), xi[0] * (2 * xi[0] - 1), xi[1] * (2 * xi[1] -
-		    1), 4 * xi[0] * xi2, 4 * xi[0] * xi[1], 4 * xi[1] * xi2};
+    return {xi2 * (2 * xi2 - 1),
+            xi[0] * (2 * xi[0] - 1),
+            xi[1] * (2 * xi[1] - 1),
+            4 * xi[0] * xi2,
+            4 * xi[0] * xi[1],
+            4 * xi[1] * xi2};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Vector2, numNodes>
-	  getLocalGradients(Vector3 const &xi) const {
-		  assert(greaterThanOrEqualZero(xi)); assert(sumsToOne(xi));
-		  const Real xi2 = 1 - xi[0] - xi[1];
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Vector2, numNodes> getLocalGradients(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+    const Real xi2 = 1 - xi[0] - xi[1];
     // clang-format off
-    return {-4*xi2+1,-4*xi2+1, 4*xi[0]-1,0, 0,4*xi[1]-1,
-	    4*(xi2-xi[0]),-4*xi[0], 4*xi[1],4*xi[0], -4*xi[1],4*(xi2-xi[1]) };
+    return {-4*xi2+1,-4*xi2+1,
+             4*xi[0]-1,0,
+             0,4*xi[1]-1,
+             4*(xi2-xi[0]),-4*xi[0],
+             4*xi[1],4*xi[0],
+             -4*xi[1],4*(xi2-xi[1]) };
     // clang-format on
-  } };
-
-struct QuadraticTetrahedronShape { static const size_t numNodes = 10; static
-	const size_t meshEntDim = 3; constexpr static Mesh_Topology
-		DofHolders[2] = {Vertex, Edge}; constexpr static size_t
-		NumDofHolders[2] = {4, 6}; constexpr static size_t
-		DofsPerHolder[2] = {1, 1}; constexpr static size_t Order = 2;
-  // ordering taken from mfem see mfem/mfem fem/fe/fe_fixed_order.cpp @597cba8
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Real, numNodes> getValues(Vector4 const
-		  &xi) const { assert(greaterThanOrEqualZero(xi));
-	  assert(sumsToOne(xi)); const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
+  }
+};
+
+struct QuadraticTetrahedronShape {
+  static const size_t numNodes = 10;
+  static const size_t meshEntDim = 3;
+  constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
+  constexpr static size_t NumDofHolders[2] = {4, 6};
+  constexpr static size_t DofsPerHolder[2] = {1, 1};
+  constexpr static size_t Order = 2;
+  // ordering taken from mfem
+  // see mfem/mfem fem/fe/fe_fixed_order.cpp @597cba8
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector4 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+    const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
     // clang-format off
-    return {xi3*(2*xi3-1), xi[0]*(2*xi[0]-1), xi[1]*(2*xi[1]-1),
-	    xi[2]*(2*xi[2]-1), 4*xi[0]*xi3, 4*xi[1]*xi3, 4*xi[2]*xi3,
-	    4*xi[0]*xi[1], 4*xi[2]*xi[0], 4*xi[1]*xi[2]};
+    return {xi3*(2*xi3-1),
+            xi[0]*(2*xi[0]-1),
+            xi[1]*(2*xi[1]-1),
+            xi[2]*(2*xi[2]-1),
+            4*xi[0]*xi3,
+            4*xi[1]*xi3,
+            4*xi[2]*xi3,
+            4*xi[0]*xi[1],
+            4*xi[2]*xi[0],
+            4*xi[1]*xi[2]};
     // clang-format on
   }
 
-  KOKKOS_INLINE_FUNCTION Kokkos::Array<Vector3, numNodes>
-	  getLocalGradients(Vector4 const &xi) const {
-		  assert(greaterThanOrEqualZero(xi)); assert(sumsToOne(xi));
-		  const Real xi3 = 1 - xi[0] - xi[1] - xi[2]; const Real d3 = 1
-			  - 4 * xi3;
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Vector3, numNodes> getLocalGradients(Vector4 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+    const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
+    const Real d3 = 1 - 4 * xi3;
     // clang-format off
-    return {d3,d3,d3, 4*xi[0]-1,0,0, 0,4*xi[1]-1,0, 0,0,4*xi[2]-1,
-	    4*xi3-4*xi[0],-4*xi[0],-4*xi[0], -4*xi[1],4*xi3-4*xi[1],-4*xi[1],
-	    -4*xi[2],-4*xi[2],4*xi3-4*xi[2], 4*xi[1],4*xi[0],0,
-	    4*xi[2],0,4*xi[0], 0,4*xi[2],4*xi[1]};
+    return {d3,d3,d3,
+            4*xi[0]-1,0,0,
+            0,4*xi[1]-1,0,
+            0,0,4*xi[2]-1,
+            4*xi3-4*xi[0],-4*xi[0],-4*xi[0],
+            -4*xi[1],4*xi3-4*xi[1],-4*xi[1],
+            -4*xi[2],-4*xi[2],4*xi3-4*xi[2],
+            4*xi[1],4*xi[0],0,
+            4*xi[2],0,4*xi[0],
+            0,4*xi[2],4*xi[1]};
     // clang-format on
-  } };
+  }
+};
+
 
 struct ReducedQuinticImplicitShape {
   static const size_t numNodes = 21;
@@ -156,8 +221,8 @@ struct ReducedQuinticImplicitShape {
     assert(sumsToOne(xi));
 
     const Real L1 = 1.0 - xi[0] - xi[1];
-    const Real L2 = xi[0];               
-    const Real L3 = xi[1];               
+    const Real L2 = xi[0];
+    const Real L3 = xi[1];
 
     Real powL1[6], powL2[6], powL3[6];
     powL1[0] = powL2[0] = powL3[0] = 1.0;
@@ -224,8 +289,8 @@ struct ReducedQuinticImplicitShape {
         if (j > 0) dN_dL2 = coeff * double(j) * powL1[i] * powL2[j - 1] * powL3[k];
         if (k > 0) dN_dL3 = coeff * double(k) * powL1[i] * powL2[j] * powL3[k - 1];
 
-        const double dNdX = -dN_dL1 + dN_dL2; // xi[0] corresponds to L2
-        const double dNdY = -dN_dL1 + dN_dL3; // xi[1] corresponds to L3
+        const double dNdX = -dN_dL1 + dN_dL2;
+        const double dNdY = -dN_dL1 + dN_dL3;
 
         dN[idx][0] = dNdX;
         dN[idx][1] = dNdY;
diff --git a/test/testReducedQuinticImplicitIntegrator.cpp b/test/testReducedQuinticImplicitIntegrator.cpp
index 73b1f73..6aca877 100644
--- a/test/testReducedQuinticImplicitIntegrator.cpp
+++ b/test/testReducedQuinticImplicitIntegrator.cpp
@@ -62,7 +62,7 @@ class ReducedQuinticIntegrator : public MeshField::Integrator {
     double local = 0.0;
     Kokkos::parallel_reduce("RQ_integrate", npts,
       KOKKOS_LAMBDA(const int i, double &acc) {
-        acc += double(w(i)) * double(dV(i));
+        acc += double(w(i));
       }, local);
 
     total += local;
@@ -130,23 +130,24 @@ void runTest(Omega_h::Mesh &mesh) {
 
 int main(int argc, char **argv) {
   Kokkos::initialize(argc, argv);
-  Omega_h::Library lib(&argc, &argv);
+  {
+    Omega_h::Library lib(&argc, &argv);
 
-  Omega_h::Reals coords({
-    0.0, 0.0,
-    1.0, 0.0,
-    0.0, 1.0
-  });
+    Omega_h::Reals coords({
+        0.0, 0.0,
+        1.0, 0.0,
+        0.0, 1.0
+    });
 
-  Omega_h::LOs verts({0,1,2});
+    Omega_h::LOs verts({0,1,2});
 
-  Omega_h::Mesh mesh(&lib);
-  Omega_h::build_from_elems_and_coords(
-      &mesh, OMEGA_H_SIMPLEX, 2, verts, coords);
+    Omega_h::Mesh mesh(&lib);
+    Omega_h::build_from_elems_and_coords(
+        &mesh, OMEGA_H_SIMPLEX, 2, verts, coords);
 
-  runTest<MeshField::KokkosController>(mesh);
+    runTest<MeshField::KokkosController>(mesh);
+    }
 
   Kokkos::finalize();
   return 0;
 }
-