From e401f7ec30f90b01f859c188c415694d7878811e Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Wed, 18 Jun 2025 15:57:31 -0400
Subject: [PATCH 01/14] Added num comp parameter to template

---
 src/MeshField.hpp            |  4 ++--
 src/MeshField_ShapeField.hpp | 19 ++++++++++---------
 test/testElement.cpp         |  8 ++++----
 test/testOmegahElement.cpp   |  6 +++---
 4 files changed, 19 insertions(+), 18 deletions(-)

diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index 69844ae..6bcaaed 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -189,11 +189,11 @@ class OmegahMeshField {
     static_assert(dim == 1 || dim == 2 || dim == 3);
   }
 
-  template <typename DataType, size_t order>
+  template <typename DataType, size_t order, size_t numComp>
   // Ordering of field indexing changed to 'entity, node, component'
   auto CreateLagrangeField() {
     return MeshField::CreateLagrangeField<ExecutionSpace, Controller, DataType,
-                                          order, dim>(meshInfo);
+                                          order, dim, numComp>(meshInfo);
   }
 
   auto getCoordField() { return coordField; }
diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp
index 4bd0254..b7c1bff 100644
--- a/src/MeshField_ShapeField.hpp
+++ b/src/MeshField_ShapeField.hpp
@@ -158,7 +158,7 @@ template <typename VtxAccessor> struct LinearAccessor {
 template <typename ExecutionSpace,
           template <typename...>
           typename Controller = MeshField::KokkosController,
-          typename DataType, size_t order, size_t dim>
+          typename DataType, size_t order, size_t dim, size_t numComp>
 auto CreateLagrangeField(const MeshInfo &meshInfo) {
   static_assert((std::is_same_v<Real4, DataType> == true ||
                  std::is_same_v<Real8, DataType> == true),
@@ -179,10 +179,10 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
         std::is_same_v<
             Controller<ExecutionSpace, MemorySpace, DataType>,
             MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
-        Controller<ExecutionSpace, MemorySpace, DataType[1][1]>,
+        Controller<ExecutionSpace, MemorySpace, DataType[1][numComp]>,
         Controller<MemorySpace, ExecutionSpace, DataType ***>>;
     // 1 dof with 1 component per vtx
-    auto createController = [](const int numComp, auto numVtx) {
+    auto createController = [](auto numVtx) {
       if constexpr (std::is_same_v<
                         Controller<ExecutionSpace, MemorySpace, DataType>,
                         MeshField::CabanaController<ExecutionSpace, MemorySpace,
@@ -192,7 +192,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
         return Ctrlr({/*field 0*/ numVtx, 1, numComp});
       }
     };
-    Ctrlr kk_ctrl = createController(1, meshInfo.numVtx);
+    Ctrlr kk_ctrl = createController(meshInfo.numVtx);
 #else
     using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
     Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1});
@@ -214,10 +214,11 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
         std::is_same_v<
             Controller<ExecutionSpace, MemorySpace, DataType>,
             MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
-        Controller<ExecutionSpace, MemorySpace, DataType[1][1], DataType[1][1]>,
+        Controller<ExecutionSpace, MemorySpace, DataType[1][numComp],
+                   DataType[1][numComp]>,
         Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>>;
     // 1 dof with 1 comp per vtx/edge
-    auto createController = [](const int numComp, auto numVtx, auto numEdge) {
+    auto createController = [](auto numVtx, auto numEdge) {
       if constexpr (std::is_same_v<
                         Controller<ExecutionSpace, MemorySpace, DataType>,
                         MeshField::CabanaController<ExecutionSpace, MemorySpace,
@@ -228,12 +229,12 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
                       /*field 1*/ numEdge, 1, numComp});
       }
     };
-    Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge);
+    Ctrlr kk_ctrl = createController(meshInfo.numVtx, meshInfo.numEdge);
 #else
     using Ctrlr =
         Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>;
-    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1,
-                   /*field 1*/ meshInfo.numEdge, 1, 1});
+    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp,
+                   /*field 1*/ meshInfo.numEdge, 1, numComp});
 #endif
     auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
     auto edgeField = MeshField::makeField<Ctrlr, 1>(kk_ctrl);
diff --git a/test/testElement.cpp b/test/testElement.cpp
index 57984ae..4b0340d 100644
--- a/test/testElement.cpp
+++ b/test/testElement.cpp
@@ -41,7 +41,7 @@ void triangleLocalPointEval() {
   const auto numElms = 3; // provided by the mesh
   const MeshField::MeshInfo meshInfo{.numVtx = 5, .numTri = 3};
   auto field = MeshField::CreateLagrangeField<
-      ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 2>(
+      ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 2, 1>(
       meshInfo);
 
   MeshField::FieldElement f(numElms, field, MeshField::LinearTriangleShape(),
@@ -85,7 +85,7 @@ struct LinearEdgeToVertexField {
 void edgeLocalPointEval() {
   const MeshField::MeshInfo meshInfo{.numVtx = 5, .numEdge = 7, .dim = 1};
   auto field = MeshField::CreateLagrangeField<
-      ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 1>(
+      ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 1, 1>(
       meshInfo);
 
   MeshField::FieldElement f(meshInfo.numEdge, field,
@@ -137,7 +137,7 @@ void quadraticTriangleLocalPointEval() {
   const MeshField::MeshInfo meshInfo{
       .numVtx = 3, .numEdge = 3, .numTri = 1, .dim = 2};
   auto field = MeshField::CreateLagrangeField<
-      ExecutionSpace, MeshField::KokkosController, MeshField::Real, 2, 2>(
+      ExecutionSpace, MeshField::KokkosController, MeshField::Real, 2, 2, 1>(
       meshInfo);
 
   MeshField::FieldElement f(meshInfo.numTri, field,
@@ -192,7 +192,7 @@ void quadraticTetrahedronLocalPointEval() {
   auto field =
       MeshField::CreateLagrangeField<ExecutionSpace,
                                      MeshField::KokkosController,
-                                     MeshField::Real, ShapeOrder, MeshDim>(
+                                     MeshField::Real, ShapeOrder, MeshDim, 1>(
           meshInfo);
 
   MeshField::FieldElement f(meshInfo.numTet, field,
diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 09d8ab7..104b0a2 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -172,7 +172,7 @@ void doRun(Omega_h::Mesh &mesh,
       const auto ShapeOrder = 1;
 
       auto field =
-          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder>();
+          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder, 1>();
       LinearFunction func = LinearFunction();
       setVertices(mesh, func, field);
       using FieldType = decltype(field);
@@ -187,7 +187,7 @@ void doRun(Omega_h::Mesh &mesh,
     {
       const auto ShapeOrder = 2;
       auto field =
-          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder>();
+          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder, 1>();
       auto func = QuadraticFunction();
       setVertices(mesh, func, field);
       setEdges(mesh, func, field);
@@ -203,7 +203,7 @@ void doRun(Omega_h::Mesh &mesh,
     {
       const auto ShapeOrder = 2;
       auto field =
-          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder>();
+          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder, 1>();
       auto func = LinearFunction();
       setVertices(mesh, func, field);
       setEdges(mesh, func, field);

From 3e05f711230e91af7fd193ebada7af41bb98687a Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Tue, 24 Jun 2025 12:26:19 -0400
Subject: [PATCH 02/14] Adding templating for number of components

---
 src/MeshField.hpp          |  5 +++--
 src/MeshField_Element.hpp  | 11 ++++++-----
 test/testOmegahElement.cpp | 22 +++++++++++++++++++---
 3 files changed, 28 insertions(+), 10 deletions(-)

diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index 6bcaaed..8e54b72 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -220,7 +220,7 @@ class OmegahMeshField {
   }
 
   // evaluate a field at the specified local coordinate for each triangle
-  template <typename ViewType, typename ShapeField>
+  template <typename ViewType, typename ShapeField, size_t numComp>
   auto triangleLocalPointEval(ViewType localCoords, size_t NumPtsPerElem,
                               ShapeField field) {
     const auto MeshDim = 2;
@@ -234,7 +234,8 @@ class OmegahMeshField {
 
     const auto [shp, map] = Omegah::getTriangleElement<ShapeOrder>(mesh);
 
-    MeshField::FieldElement f(meshInfo.numTri, field, shp, map);
+    MeshField::FieldElement<ShapeField, decltype(shp), decltype(map), numComp>
+        f(meshInfo.numTri, field, shp, map);
     auto offsets = createOffsets(meshInfo.numTri, NumPtsPerElem);
     auto eval = MeshField::evaluate(f, localCoords, offsets);
     return eval;
diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp
index 06d79a6..0acf938 100644
--- a/src/MeshField_Element.hpp
+++ b/src/MeshField_Element.hpp
@@ -124,7 +124,8 @@ struct ElementToDofHolderMap {
  * @param elmIn see ElementType
  */
 template <typename FieldAccessor, typename ShapeType,
-          typename ElementDofHolderAccessor>
+          typename ElementDofHolderAccessor,
+          size_t numComponentsPerDof = ShapeType::numComponentsPerDof>
 struct FieldElement {
   const size_t numMeshEnts;
   const FieldAccessor field;
@@ -150,8 +151,8 @@ struct FieldElement {
   };
   using ValArray =
       Kokkos::Array<typename baseType<typename FieldAccessor::BaseType>::type,
-                    ShapeType::numComponentsPerDof>;
-  static const size_t NumComponents = ShapeType::numComponentsPerDof;
+                    numComponentsPerDof>;
+  static const size_t NumComponents = numComponentsPerDof;
 
   /**
    * @brief
@@ -173,11 +174,11 @@ struct FieldElement {
     assert(ent < numMeshEnts);
     ValArray c;
     const auto shapeValues = shapeFn.getValues(localCoord);
-    for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci)
+    for (int ci = 0; ci < numComponentsPerDof; ++ci)
       c[ci] = 0;
     for (auto topo : elm2dof.getTopology()) { // element topology
       for (int ni = 0; ni < shapeFn.numNodes; ++ni) {
-        for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci) {
+        for (int ci = 0; ci < numComponentsPerDof; ++ci) {
           auto map = elm2dof(ni, ci, ent, topo);
           const auto fval =
               field(map.entity, map.node, map.component, map.topo);
diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 104b0a2..2612072 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -176,7 +176,7 @@ void doRun(Omega_h::Mesh &mesh,
       LinearFunction func = LinearFunction();
       setVertices(mesh, func, field);
       using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
+      auto result = omf.template triangleLocalPointEval<ViewType, FieldType, 1>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 LinearFunction{});
@@ -192,7 +192,7 @@ void doRun(Omega_h::Mesh &mesh,
       setVertices(mesh, func, field);
       setEdges(mesh, func, field);
       using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
+      auto result = omf.template triangleLocalPointEval<ViewType, FieldType, 1>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 QuadraticFunction{});
@@ -208,13 +208,29 @@ void doRun(Omega_h::Mesh &mesh,
       setVertices(mesh, func, field);
       setEdges(mesh, func, field);
       using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
+      auto result = omf.template triangleLocalPointEval<ViewType, FieldType, 1>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 LinearFunction{});
       if (failed)
         doFail("quadratic", "linear", testCase.name);
     }
+
+    {
+      const auto ShapeOrder = 1;
+
+      auto field =
+          omf.template CreateLagrangeField<MeshField::Real, ShapeOrder, 2>();
+      LinearFunction func = LinearFunction();
+      setVertices(mesh, func, field);
+      using FieldType = decltype(field);
+      auto result = omf.template triangleLocalPointEval<ViewType, FieldType, 2>(
+          testCase.coords, testCase.NumPtsPerElem, field);
+      auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
+                                LinearFunction{});
+      if (failed)
+        doFail("linear", "linear", testCase.name);
+    }
   }
 }
 

From 390c5608425b9e1faf0aa95b5dfd4eb1f26a0725 Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Wed, 25 Jun 2025 14:04:20 -0400
Subject: [PATCH 03/14] Tests properly check for multiple components

---
 test/testOmegahElement.cpp | 42 +++++++++++++++++++++++---------------
 1 file changed, 25 insertions(+), 17 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 6ceec9d..16dbbcb 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -46,7 +46,7 @@ struct TestCoords {
 
 template <typename Result, typename CoordField, typename AnalyticFunction>
 bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
-                 TestCoords testCase, AnalyticFunction func) {
+                 TestCoords testCase, AnalyticFunction func, size_t numComp) {
   const auto numPtsPerElem = testCase.NumPtsPerElem;
   MeshField::FieldElement fcoords(
       mesh.nfaces(), coordField, MeshField::LinearTriangleCoordinateShape(),
@@ -63,16 +63,20 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
         for (auto pt = first; pt < last; pt++) {
           const auto x = globalCoords(pt, 0);
           const auto y = globalCoords(pt, 1);
-          const auto expected = func(x, y);
-          const auto computed = result(pt, 0);
-          MeshField::LO isError = 0;
-          if (Kokkos::fabs(computed - expected) > MeshField::MachinePrecision) {
-            isError = 1;
-            Kokkos::printf("result for elm %d, pt %d, does not match: expected "
-                           "%f computed %f\n",
-                           ent, pt, expected, computed);
+          for (int i = 0; i < numComp; ++i) {
+            const auto expected = func(x, y);
+            const auto computed = result(pt, i);
+            MeshField::LO isError = 0;
+            if (Kokkos::fabs(computed - expected) >
+                MeshField::MachinePrecision) {
+              isError = 1;
+              Kokkos::printf(
+                  "result for elm %d, pt %d, does not match: expected "
+                  "%f computed %f\n",
+                  ent, pt, expected, computed);
+            }
+            lerrors += isError;
           }
-          lerrors += isError;
         }
       },
       numErrors);
@@ -88,7 +92,9 @@ void setVertices(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) {
     // - TODO should be encoded in the field?
     const auto x = coords[vtx * MeshDim];
     const auto y = coords[vtx * MeshDim + 1];
-    field(vtx, 0, 0, MeshField::Vertex) = func(x, y);
+    for (int i = 0; i < field.numComp; ++i) {
+      field(vtx, 0, i, MeshField::Vertex) = func(x, y);
+    }
   };
   MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()},
                           setFieldAtVertices, "setFieldAtVertices");
@@ -109,7 +115,9 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) {
     const auto x = (coords[left * MeshDim] + coords[right * MeshDim]) / 2.0;
     const auto y =
         (coords[left * MeshDim + 1] + coords[right * MeshDim + 1]) / 2.0;
-    field(edge, 0, 0, MeshField::Edge) = func(x, y);
+    for (int i = 0; i < field.numComp; ++i) {
+      field(edge, 0, 0, MeshField::Edge) = func(x, y);
+    }
   };
   MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()},
                           setFieldAtEdges, "setFieldAtEdges");
@@ -179,7 +187,7 @@ void doRun(Omega_h::Mesh &mesh,
       auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{});
+                                LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name);
     }
@@ -196,7 +204,7 @@ void doRun(Omega_h::Mesh &mesh,
       auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                QuadraticFunction{});
+                                QuadraticFunction{}, numComponents);
       if (failed)
         doFail("quadratic", "quadratic", testCase.name);
     }
@@ -213,7 +221,7 @@ void doRun(Omega_h::Mesh &mesh,
       auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{});
+                                LinearFunction{}, numComponents);
       if (failed)
         doFail("quadratic", "linear", testCase.name);
     }
@@ -229,7 +237,7 @@ void doRun(Omega_h::Mesh &mesh,
       auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{});
+                                LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name);
     }
@@ -244,7 +252,7 @@ void doRun(Omega_h::Mesh &mesh,
       auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
           testCase.coords, testCase.NumPtsPerElem, field);
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{});
+                                LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name);
     }

From 385993b21a5f1cb9b55a42d4441e02f09b655c6c Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Wed, 25 Jun 2025 15:47:34 -0400
Subject: [PATCH 04/14] adjusted fail message

---
 test/testOmegahElement.cpp | 21 +++++++++++++--------
 1 file changed, 13 insertions(+), 8 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 16dbbcb..4376f45 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -141,10 +141,10 @@ createElmAreaCoords(size_t numElements,
 }
 
 void doFail(std::string_view order, std::string_view function,
-            std::string_view location) {
+            std::string_view location, std::string_view numComp) {
   std::stringstream ss;
-  ss << order << " field evaluation with " << function
-     << " analytic function at " << location << " points failed\n";
+  ss << order << " field evaluation with " << numComp << "components and "
+     << function << " analytic function at " << location << " points failed\n";
   std::string msg = ss.str();
   MeshField::fail(msg);
 }
@@ -189,7 +189,8 @@ void doRun(Omega_h::Mesh &mesh,
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 LinearFunction{}, numComponents);
       if (failed)
-        doFail("linear", "linear", testCase.name);
+        doFail("linear", "linear", testCase.name,
+               std::to_string(numComponents));
     }
 
     {
@@ -206,7 +207,8 @@ void doRun(Omega_h::Mesh &mesh,
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 QuadraticFunction{}, numComponents);
       if (failed)
-        doFail("quadratic", "quadratic", testCase.name);
+        doFail("quadratic", "quadratic", testCase.name,
+               std::to_string(numComponents));
     }
 
     {
@@ -223,7 +225,8 @@ void doRun(Omega_h::Mesh &mesh,
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 LinearFunction{}, numComponents);
       if (failed)
-        doFail("quadratic", "linear", testCase.name);
+        doFail("quadratic", "linear", testCase.name,
+               std::to_string(numComponents));
     }
 
     {
@@ -239,7 +242,8 @@ void doRun(Omega_h::Mesh &mesh,
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 LinearFunction{}, numComponents);
       if (failed)
-        doFail("linear", "linear", testCase.name);
+        doFail("linear", "linear", testCase.name,
+               std::to_string(numComponents));
     }
     {
       const auto ShapeOrder = 1;
@@ -254,7 +258,8 @@ void doRun(Omega_h::Mesh &mesh,
       auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
                                 LinearFunction{}, numComponents);
       if (failed)
-        doFail("linear", "linear", testCase.name);
+        doFail("linear", "linear", testCase.name,
+               std::to_string(numComponents));
     }
   }
 }

From 85750023497085270465567771d64c3b3c46142a Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Wed, 25 Jun 2025 17:11:49 -0400
Subject: [PATCH 05/14] adding debug print

---
 test/testOmegahElement.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 4376f45..7709597 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -80,6 +80,7 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
         }
       },
       numErrors);
+  printf("numErrors: %d\n", numErrors);
   return (numErrors > 0);
 }
 
@@ -143,7 +144,7 @@ createElmAreaCoords(size_t numElements,
 void doFail(std::string_view order, std::string_view function,
             std::string_view location, std::string_view numComp) {
   std::stringstream ss;
-  ss << order << " field evaluation with " << numComp << "components and "
+  ss << order << " field evaluation with " << numComp << " components and "
      << function << " analytic function at " << location << " points failed\n";
   std::string msg = ss.str();
   MeshField::fail(msg);

From c9ed213a246813262fa37ea2d212c491e039a826 Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 10:53:52 -0400
Subject: [PATCH 06/14] Testing

---
 test/testOmegahElement.cpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 7709597..6abc4a8 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -80,8 +80,7 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
         }
       },
       numErrors);
-  printf("numErrors: %d\n", numErrors);
-  return (numErrors > 0);
+  return true;
 }
 
 template <typename AnalyticFunction, typename ShapeField>

From f9b9d5e9c807ee84d131edcf1a8455f876e9daac Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 11:11:36 -0400
Subject: [PATCH 07/14] Testing

---
 test/testOmegahElement.cpp | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 6abc4a8..bd241d2 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -45,7 +45,7 @@ struct TestCoords {
 };
 
 template <typename Result, typename CoordField, typename AnalyticFunction>
-bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
+MeshField::LO checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
                  TestCoords testCase, AnalyticFunction func, size_t numComp) {
   const auto numPtsPerElem = testCase.NumPtsPerElem;
   MeshField::FieldElement fcoords(
@@ -80,7 +80,7 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
         }
       },
       numErrors);
-  return true;
+  return numErrors;
 }
 
 template <typename AnalyticFunction, typename ShapeField>
@@ -141,10 +141,10 @@ createElmAreaCoords(size_t numElements,
 }
 
 void doFail(std::string_view order, std::string_view function,
-            std::string_view location, std::string_view numComp) {
+            std::string_view location, std::string_view numComp, std::string_view numErrors) {
   std::stringstream ss;
   ss << order << " field evaluation with " << numComp << " components and "
-     << function << " analytic function at " << location << " points failed\n";
+     << function << " analytic function at " << location << " points failed with " << numErrors << std::endl;
   std::string msg = ss.str();
   MeshField::fail(msg);
 }
@@ -190,7 +190,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents));
+               std::to_string(numComponents), std::to_string(failed));
     }
 
     {
@@ -208,7 +208,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 QuadraticFunction{}, numComponents);
       if (failed)
         doFail("quadratic", "quadratic", testCase.name,
-               std::to_string(numComponents));
+               std::to_string(numComponents), std::to_string(failed));
     }
 
     {
@@ -226,7 +226,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("quadratic", "linear", testCase.name,
-               std::to_string(numComponents));
+               std::to_string(numComponents), std::to_string(failed));
     }
 
     {
@@ -243,7 +243,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents));
+               std::to_string(numComponents), std::to_string(failed));
     }
     {
       const auto ShapeOrder = 1;
@@ -259,7 +259,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents));
+               std::to_string(numComponents), std::to_string(failed));
     }
   }
 }

From fef38e8889cacc27095c127380517f993a0f23b9 Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 11:17:18 -0400
Subject: [PATCH 08/14] Testing

---
 test/testOmegahElement.cpp | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index bd241d2..5ba9bb9 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -45,7 +45,7 @@ struct TestCoords {
 };
 
 template <typename Result, typename CoordField, typename AnalyticFunction>
-MeshField::LO checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
+bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
                  TestCoords testCase, AnalyticFunction func, size_t numComp) {
   const auto numPtsPerElem = testCase.NumPtsPerElem;
   MeshField::FieldElement fcoords(
@@ -65,7 +65,7 @@ MeshField::LO checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordFi
           const auto y = globalCoords(pt, 1);
           for (int i = 0; i < numComp; ++i) {
             const auto expected = func(x, y);
-            const auto computed = result(pt, i);
+            const auto computed = result(pt, 0);
             MeshField::LO isError = 0;
             if (Kokkos::fabs(computed - expected) >
                 MeshField::MachinePrecision) {
@@ -80,7 +80,7 @@ MeshField::LO checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordFi
         }
       },
       numErrors);
-  return numErrors;
+  return (numErrors > 0);
 }
 
 template <typename AnalyticFunction, typename ShapeField>
@@ -141,10 +141,10 @@ createElmAreaCoords(size_t numElements,
 }
 
 void doFail(std::string_view order, std::string_view function,
-            std::string_view location, std::string_view numComp, std::string_view numErrors) {
+            std::string_view location, std::string_view numComp) {
   std::stringstream ss;
   ss << order << " field evaluation with " << numComp << " components and "
-     << function << " analytic function at " << location << " points failed with " << numErrors << std::endl;
+     << function << " analytic function at " << location << " points failed\n";
   std::string msg = ss.str();
   MeshField::fail(msg);
 }
@@ -190,7 +190,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents), std::to_string(failed));
+               std::to_string(numComponents));
     }
 
     {
@@ -208,7 +208,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 QuadraticFunction{}, numComponents);
       if (failed)
         doFail("quadratic", "quadratic", testCase.name,
-               std::to_string(numComponents), std::to_string(failed));
+               std::to_string(numComponents));
     }
 
     {
@@ -226,7 +226,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("quadratic", "linear", testCase.name,
-               std::to_string(numComponents), std::to_string(failed));
+               std::to_string(numComponents));
     }
 
     {
@@ -243,7 +243,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents), std::to_string(failed));
+               std::to_string(numComponents));
     }
     {
       const auto ShapeOrder = 1;
@@ -259,7 +259,7 @@ void doRun(Omega_h::Mesh &mesh,
                                 LinearFunction{}, numComponents);
       if (failed)
         doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents), std::to_string(failed));
+               std::to_string(numComponents));
     }
   }
 }

From feeca54ab3be172caea0133a4679cab72cb66373 Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 11:21:48 -0400
Subject: [PATCH 09/14] Testing

---
 test/testOmegahElement.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 5ba9bb9..e8e4725 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -63,9 +63,9 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
         for (auto pt = first; pt < last; pt++) {
           const auto x = globalCoords(pt, 0);
           const auto y = globalCoords(pt, 1);
+          const auto expected = func(x, y);
           for (int i = 0; i < numComp; ++i) {
-            const auto expected = func(x, y);
-            const auto computed = result(pt, 0);
+            const auto computed = result(pt, i);
             MeshField::LO isError = 0;
             if (Kokkos::fabs(computed - expected) >
                 MeshField::MachinePrecision) {

From 35dcc3326efb89446dcd1b1c11d16e26b887bb98 Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 11:26:43 -0400
Subject: [PATCH 10/14] Testing

---
 test/testOmegahElement.cpp | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index e8e4725..e6a32f8 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -64,8 +64,7 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
           const auto x = globalCoords(pt, 0);
           const auto y = globalCoords(pt, 1);
           const auto expected = func(x, y);
-          for (int i = 0; i < numComp; ++i) {
-            const auto computed = result(pt, i);
+            const auto computed = result(pt, 0);
             MeshField::LO isError = 0;
             if (Kokkos::fabs(computed - expected) >
                 MeshField::MachinePrecision) {
@@ -77,7 +76,6 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
             }
             lerrors += isError;
           }
-        }
       },
       numErrors);
   return (numErrors > 0);

From 55ee87ecc06bc31a40a9d72ea51ac93e4a3a955f Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 11:31:23 -0400
Subject: [PATCH 11/14] Testing

---
 test/.testOmegahElement.cpp.swp | Bin 0 -> 36864 bytes
 test/testOmegahElement.cpp      |   4 ++--
 2 files changed, 2 insertions(+), 2 deletions(-)
 create mode 100644 test/.testOmegahElement.cpp.swp

diff --git a/test/.testOmegahElement.cpp.swp b/test/.testOmegahElement.cpp.swp
new file mode 100644
index 0000000000000000000000000000000000000000..456d0f1eda2e3ff5bb479f0a0926aac374da7a1b
GIT binary patch
literal 36864
zcmeI53y>Ved4LCS>=1-;;wl^<gysawp0wIKolJsH(ov_=354!+LJ~?ULR{u{XYVH3
z-C1UKPdYBbBz74mPDODN6G)sGo1_A9Yy%ZSB{nuLQ?`pL5<BIBI2Ee!BaY(-sA6Mm
zqS*QV?wQ%y$4Ms;225|=w~y(bp6<V2^Z)(#+|1zKZIkNq!jQx34UY56<^6N_{fFax
za{YSeKqc_%zo|ncbjeq7&MO}9!l)1xgU~A$!n$AZO3kkLH^m7KRJ~|!tM64x(H3uB
z|0?TW0}U@~=(nk=H|x%gSG=lMYZQw0dRJohfR%ui!15)~Y?RktezvoIXx*Tkw6gDF
z^~RsOcKPSqm{tN-0#*W60#*W60#*W60#*W60<U-pH0rN&-byJ>k4rKZf1cLy`P%q*
zU;KGxNBB?0zw_~@-x2=0`1k7g^F1BmuZe#rpP%Rm|LOSmug0IB>Igq0{+(Rs{I~!~
zIpSdZVI^QCU?pHBU?pHBU?pHBU?pHBU?pHBU?uRXlz>}uoMR~Omx$ub`hTMTKY5Ph
zd<JfZ1F#9sfu~>ZIR65lg%83JsKd1|3I$jT&z|i#{{Z*HU%=bpEpRbh2<O9j@WU03
z^Pli%a3jn^1R-pIi{SfbInLMM18@ht58egW!&VrA0Vu#%e#UYB6t03d!qGDw=T0cY
zPFMkF!M89dd>%do0bBvEhbNH3N8q>NS`fK>597y!a0E(l8JrIyx4(x!fa_r=48kRF
zKKv~F2!qDs@F08?{t0e|0ayXw#&B^jEW#wLgl90ud<Krd@50+*GyEUMjQ@q_;Gx*C
z(eSGEireti#>PV3tGQKAjn>@CLc=d^ZPtnnKd9wZ8nk=Pt$W(@V#IL{J8H>S<TZ4(
z(Z*gcBrPw>tG@1oef=9qWHNE*7#^;>p<AhVmFaR2=Ee_s#irz7cik;|x&FLb96ZdM
zxPKIC^CTSQ`VZ$x;jkpbhkjQ|dxPZW8;A00Ft7Z4_6m|8?@+(mq{@^pN0bGf3C;>a
zmGf&2<r6r#fsc(Qr9!P)9Sf@U4QkCAzh5ots}zG;)KG4-5vT>?7lR-yMf<p>TCFd)
z#jmP0YTbU5vrOzmOT><24h*QZYR}Y`DHY8H%}Pnlc*?64gOXQLeoZyzJd>#_bS{WJ
z>ors<D9;6zk{7Ca5cwi9<u-KOs$Z%HQmCM;ywa@4Svt?%bwReaPi@~mwR?JEbaL~S
z(VWS%oV$u*nR4rMO<9SrE5iKBX(48cDU~WRh+~!fszjL&{F36@Nvp_dN1f@p;bDod
zR!NHat`?}PDz93t&xk;>^a9&&=rx*QP35TJ@h}X+NR6mLU6guA$;-LEG-Oy^TDq>*
zmsjbLxZS88M|7)H`js?F`H}suRK+uwT(e2}k-ofXEatq<*ybqiOnFrY?@-+<HoQ`v
zrdO>u`IcrUIbS~eE-j}K)TL!9fRhb3-}^!@YE~M$3wV_(C6$Zlk5b@8suna<)om2#
zhSTe+zK(RZ2L}!W5nb2_t|K%qAkAT2dm$wQQulhepc%NuIlt!Z#1G7mXvax+lW6kt
z%PN;#sqD@~v{Dn${hmZiYBG@%6<O|%?M4wJL(T{yClaQsZ`!0X2f9y^BbhdzT!B(%
zPV=+poiMg(eaR`A=t+f>c<0DvVrSyYh%{RXX57k{9+z}^*0nW|Ol%pEx?*SZp+`pD
z^=TAb-VdY325BH0RmC-5*R09bR3>qWketJSNaxwvh<19R_~x}*LgoaWI%m+*8A^M`
zUviw%y(80_CZ^P?xR<rE*Qfmrcj*>Sca|ZiP%AUTFq=r{NZU;^7lk+@GXD_No^(18
zklVox4&n03#PR0cN|VVoS2lwrGbs*=V{YUXk^^n!(-oxUoUIPE$<YaCQ#SsM!^0Ew
z6*t@y`fhEu;_V5{?56utP72na)_Hp|&3vWfwat$0>FLJf*C>G=dX0E8<<sZKPL?Vz
zqlhU}n{DMoPN@KP^lk{Iq`dDDN$JL$p0nh{5cxNH(+z1WGEGZO&jdk5W$L;2JY(8G
z+*k8TUrv`PE1q6cj5#rLW7Ik|bFlM9O75AY$w!US@Nm>%oSap1$EA;%B6LY5X^eY)
z@8HIa=y!Q;MXlbqe*cKcYqY^3#b&XQ72Kk^NsA}BbgG?c`jyi3OmHZ7FbFFQ^<}r}
zR~GV2NZjn{x4{B8PV4ob83=T%B#0_pSD>-CRkAsspsC66D@UiVn%=!_a_7YOuSqn0
zqP`B{6@R9nzstQqnja~aDuy2IOuo{Z*InMKunB)Abi)NYvOgnLYh-9<(q#IUn=bB9
z)~MW2VUX|UmabD=LEQ`81~-=dxy~SN13M2b^vCUZky6AX&bA#B+jfjkZ@qTM*q&`u
zI}(v!({xLu+`47%olTRsLu7M7B-44p7cG0HQzy6clDB3gZEw=62H}Ffxu$NCmbS$!
zyW9)3-+T`bS9Rp+NN&Y9#0j*=IZ+ynkSqOKvC=G|ZA7$I&#jJR4*5a%F_Wq3G5o>{
zbM<;AdSBdprXwE>#ycdzeXZ2e;Id!Qr}hfd9U(r^tCjq6=KL{t#;v(yL9G!6Xe42h
ze!1k8)#Ui@tF~?%pV+c{dVI&|=85s?vC++=J4Q3-r?orXx^2~IsDXhHS<ZV`^o+eV
zFqcyC>jbe~zhY7~(T^WH@rM^tN&Mr+P6`~?Sf8niPJxL2e-?WE7>NFlx~pIJ^7}7A
zbpKhn2F5|+Z3g?Y60j1m60j1m60j1m60j1m60j1m60j1m5_pwKKulS&Axj(95<^TG
z%WDO*pjlco2Cq`Et69t10Xwzw5pDL(tDss=O}|Qh8FR4Mo!Tr<#*myED>k;9bn<Gb
z%`nz=HAnRSi_qt{i%yUJ&yuu$eVpHa2RFuH&$6aJ#P>Xiz2Lw2eJkvQemDc3WIg`@
z@Sy--V(tFD@O!WeR>Fnw3vezx!+QPa;bs_ybKxj!_K(6}!Ctrk{+adq--jDv9L|L|
zz&UUxe4TasharNs@FR>EUxa&L27U?to%Q<L;2p3BE`^KWJorBA^k0Ct!6rBd{u@C&
z1R?wie1moQyWmzR!DTQ2YoH%;@LKpP>+}!6d*Qd>5bT08;3#YM55hmfeehv;C)^HF
zt=iOD6=N#v(H-VgF`i20FFGrYr7Kxp>}pY+hOZa5MRi^6WJ=;6-x$}{(6@Lq%~=G4
zkz5<xjX~R=;ZSObj!n(ERc<(2%q#hr*HTNmq+erk-E>M5C5|DB<W<+NY4-sWQk-)`
zwR*LVpv}P)%=BUNPJJaWy&ljArp=OotVPQ*Mo?olD{3v^C|3EmV(!cJRrI0{OINLL
z>9UHfCnPJpvQ}E=43>tS<XTI8w%dpoakFU8h8Z~GHf9o)RwhsGFpDAaxD)DSH9h_m
z@Jwi}i-`M1lFKB6l9}XKx+0byXj=<mkwPBvk&bD)`pm`12jg2I^#ovHMX!cU1Vy*9
zQ=A{hS%=#wD=2&Lamg#a#3NrvWU)$K#vYy?C&iLitccr7uKluE*9GZmWD75n&djw7
zQ%51m_M>K0H<v$|_qJ*oX{DkCqBGkmk1;(IOt`ZlVlm!3*%@RPD@n0$#Z-t3Nm4oO
zi+HUTN%HHH){u##t=u{RBE)O_i;k4cw)aG&B8{|p7xm(1t5Xzgdl&h1zFrRPeQ4S3
z-I46S;0Av~vYu$i_;_)?b4+`Y$M|tP<)U~->oveHZyGQ|GV=qu|H&ErQ~#}{2LHIX
zu*`hwJ43c-q$BC`z31s#wzs!EOW*dK`YzEW{?^TmG@RwXUWzHAJ$J{wgU=g>pY#U5
z{nS?AQpdHU(92aT5dHt7=*PF9H;ewCkM;li(D9Ez6RrcX5jgN^^!z`B8MqwIgpZ-?
z-wo5S4xU4||1LZYAA%ZG-~eog95`?f`u+l}hYR2_bp7{(3m3ze(Cr_GN8us(I6Me9
z!3J0fkCSi7`<r1W?0`vl6Ff_~{ub_r8(<s!B77a4|9-d=ejToc)8L!v{dd4Sp$>at
z1w4%Ye=`)|JotI|0lNMd;cmDSeh03Hvq60KZ-)r3hGXdWpM(a;%QO@-yJ+=}`Ne8t
zuS!kqgp+ZbCT7&!L=mRdq{DqWml~}^tJy8D-y~EGt0K;;k@2O~4ba3VDba>w9bU$9
z8Yz1diWOAOwIqb|xC&}<%Qfvc`@Sc2_FQE9N3uB}ulCG^p4ajh#^R9)?mHZB`^cu0
ze9xpfuXlXp)f8^!jo9;;l+5A2jHiBbZt|H)10#YL`a#&0N)nKylesC<qCzr@nRJTS
zbVBw(upOo&J84LglFo(p!_;qhT^>%YAA}_b+2Vq)rjz7(k$C1%>$t9{jNoLr&^Ceb
zN_CVk+DH>^>=3rG`?YySchzXQ9ACIsdf5VYCdN*lG@oh7bUsbBq(0BjC*gG>fH%wA
zWkzH=;iRV0q;-o;9lzD$IyIoyORhS<WD-ehX$cv&ml@$i4e<+&JxOCf<qkoUls2gJ
zikeN2c)B-J3+Nxsy0M6<gMMR9iK99DG|FL6)!Tk>UDSKH^q?sFfs$E6x06Zwsy22F
z4==<1KbHDHj;K|#u7EjH@8i*?<CdeYq#Ezw$@)xpZ7?h7JZxlhRlLDW@x4!xWgXBJ
z+I1j_P6=CQ@U2&p?NX3kQGF7xf^AT5lhSGF9a2b7Po`{<5p_R`Bp#c(*i@w>8QWUy
z(l0x#wBU7W@eQbnS4?fU6se@Q$z&Y9jf*R(Cng0`f=rC*c?^|l2IHGYV;VL=@zpZb
zpv}TvN$9-yI6o14zRf$lQcQA`9cHqoj7l*ldxKW19a-0XyqRUQ%#iG;Gn2tSHioQL
z>j@Z&_q3TZ9{2opz0OZ|&n@Hp!CvPNX0ESiGsXbfwxYN9>Fq0RMUgEn*?lasVNn~A
z%>Ws1{W8y#W&mG$7fPjSX2@jOplw$pZEJZi>N$4UjXH&{CxRsP<7^E;<bBfA<844a
z_kOiDyo`pSyU=YKUf+53WXpjRvd8$GYi-9!@}GB3oY~Bn5hu`AIm_N~suk6?-KP4;
z3)yPwX>-tVwwg)=<{SN&8C>PU(k`XPJ^!7i(puQ&Mt6?(RlKrvF1bds!IXX$2ZsLa
zTxV!MlN$RDOyADjoUlX8k%TjcWE)uT>+E0h8vB<@-{gf|T!>(Z?W6nT3ppPmsia$H
zuC2?}BX-<@neu+Tn_{KrBIL1jAJ<j8hzyT=8{8(~L>8o5mAhp(Ixt-d4%TuKKy~I7
zg5tBk;LZ%qNR?$IWP&)gbMsAnG9mizAahF6N_ua>)kW+%d!NkyThafIplgc`F8Y6;
z)=<%H(e?it=3yVm{{Lst^=0q>`{3Ph19)&9Tn&Rz0NMZl9qa%ffVaYLz$Ne$I{!m(
zE6l+#oDRph_Ooz5)S(7dkUaInX<%Pg0#*W60#*W60#*W60#*W60;fa*soU8J9-x$o
z*JTmwy5Y_|M4+{Mc-b0o%YJSjxumr-`^crN!%FJzD*EBe>5)rfj}W7mekxMB*3q^=
zleI_M2Plca^h;_oHt0W(2Pm1xM`2+}9x>FfMw~e1&g3Vil4n&JQ-tXMt63|53Pk^3
z9S^EsM&Ew~J`3_3z<0w9a0#r0Z=&OW3T^}O1sH$};S1>b55ppOAkX<b9ln8{|2cRE
zjKFi~^moA@!yiEeyJ00<2*<EV{2SZ?qae@n`!?(S55qjj`oGTeDVYRpB31%c0#*W6
z0#*W60zbJDFgHb7r?7jV<bjE@j`ix?17$t$vKm6>d&b+%dfv(7p{?g#>^x~b@3PN)
zd?`Hd3Wbr@!`#POetpIhEMrBG*)3O%X0;V2<5{LShT<PZ1w7Kk&sX&S$C>b_@KYE4
zKl%T9&!OZ00KNnN3Ll33@FqBlp8p8^4cr8>?*9xr{#{UnYe1guFFOBM;VbY2d<<@b
zE8#SF0zLl$_!Bq)qW6CpJ^wwh8$|D40Z*Xo%RYd|;B)X1coz)Acd=W18t#CD&<{`H
z<Ntn`hDo>pzJ-7P*WgL`1pGOO-~VgiF?{^*fm`4$@I(Cj?}Y&N!x}h-U;p>uf8Z!Q
z4UfVH;U*};7+eA8!6V31o&zAy?!Orppb1*eFWn869_RJ_kIekE1mpiT=}Evz96!nE
zT+2Ka);`9-KE^<|!6n|>w^`c#eCY-c?Y*4Z)U(S0_7*<tW@jH`U?1CI{_TzS+GU1o
zA5vaPdsJCVYkO?R@>NzxhS&DBHV)|fz;qkyiFc%CH*4ER)#!h)^5cC}jdAf^dVpWv
zG+@kK_Av(1Ct9x;)ETlpqmlee5c?Pd{s9Ef=Xz0}-gN4_#4GwS1}8m5bUxhU{{ctZ
BgmVA@

literal 0
HcmV?d00001

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index e6a32f8..471da90 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -229,7 +229,7 @@ void doRun(Omega_h::Mesh &mesh,
 
     {
       const auto ShapeOrder = 1;
-      const auto numComponents = 2;
+      const auto numComponents = 1;
       auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
                                                     numComponents>();
       LinearFunction func = LinearFunction();
@@ -245,7 +245,7 @@ void doRun(Omega_h::Mesh &mesh,
     }
     {
       const auto ShapeOrder = 1;
-      const auto numComponents = 3;
+      const auto numComponents = 1;
       auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
                                                     numComponents>();
       LinearFunction func = LinearFunction();

From a156cf0c158757110a5abefb1c8df52644646876 Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Thu, 26 Jun 2025 11:43:00 -0400
Subject: [PATCH 12/14] Testing

---
 test/.testOmegahElement.cpp.swp | Bin 36864 -> 0 bytes
 test/testOmegahElement.cpp      |  10 ++++++----
 2 files changed, 6 insertions(+), 4 deletions(-)
 delete mode 100644 test/.testOmegahElement.cpp.swp

diff --git a/test/.testOmegahElement.cpp.swp b/test/.testOmegahElement.cpp.swp
deleted file mode 100644
index 456d0f1eda2e3ff5bb479f0a0926aac374da7a1b..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 36864
zcmeI53y>Ved4LCS>=1-;;wl^<gysawp0wIKolJsH(ov_=354!+LJ~?ULR{u{XYVH3
z-C1UKPdYBbBz74mPDODN6G)sGo1_A9Yy%ZSB{nuLQ?`pL5<BIBI2Ee!BaY(-sA6Mm
zqS*QV?wQ%y$4Ms;225|=w~y(bp6<V2^Z)(#+|1zKZIkNq!jQx34UY56<^6N_{fFax
za{YSeKqc_%zo|ncbjeq7&MO}9!l)1xgU~A$!n$AZO3kkLH^m7KRJ~|!tM64x(H3uB
z|0?TW0}U@~=(nk=H|x%gSG=lMYZQw0dRJohfR%ui!15)~Y?RktezvoIXx*Tkw6gDF
z^~RsOcKPSqm{tN-0#*W60#*W60#*W60#*W60<U-pH0rN&-byJ>k4rKZf1cLy`P%q*
zU;KGxNBB?0zw_~@-x2=0`1k7g^F1BmuZe#rpP%Rm|LOSmug0IB>Igq0{+(Rs{I~!~
zIpSdZVI^QCU?pHBU?pHBU?pHBU?pHBU?pHBU?uRXlz>}uoMR~Omx$ub`hTMTKY5Ph
zd<JfZ1F#9sfu~>ZIR65lg%83JsKd1|3I$jT&z|i#{{Z*HU%=bpEpRbh2<O9j@WU03
z^Pli%a3jn^1R-pIi{SfbInLMM18@ht58egW!&VrA0Vu#%e#UYB6t03d!qGDw=T0cY
zPFMkF!M89dd>%do0bBvEhbNH3N8q>NS`fK>597y!a0E(l8JrIyx4(x!fa_r=48kRF
zKKv~F2!qDs@F08?{t0e|0ayXw#&B^jEW#wLgl90ud<Krd@50+*GyEUMjQ@q_;Gx*C
z(eSGEireti#>PV3tGQKAjn>@CLc=d^ZPtnnKd9wZ8nk=Pt$W(@V#IL{J8H>S<TZ4(
z(Z*gcBrPw>tG@1oef=9qWHNE*7#^;>p<AhVmFaR2=Ee_s#irz7cik;|x&FLb96ZdM
zxPKIC^CTSQ`VZ$x;jkpbhkjQ|dxPZW8;A00Ft7Z4_6m|8?@+(mq{@^pN0bGf3C;>a
zmGf&2<r6r#fsc(Qr9!P)9Sf@U4QkCAzh5ots}zG;)KG4-5vT>?7lR-yMf<p>TCFd)
z#jmP0YTbU5vrOzmOT><24h*QZYR}Y`DHY8H%}Pnlc*?64gOXQLeoZyzJd>#_bS{WJ
z>ors<D9;6zk{7Ca5cwi9<u-KOs$Z%HQmCM;ywa@4Svt?%bwReaPi@~mwR?JEbaL~S
z(VWS%oV$u*nR4rMO<9SrE5iKBX(48cDU~WRh+~!fszjL&{F36@Nvp_dN1f@p;bDod
zR!NHat`?}PDz93t&xk;>^a9&&=rx*QP35TJ@h}X+NR6mLU6guA$;-LEG-Oy^TDq>*
zmsjbLxZS88M|7)H`js?F`H}suRK+uwT(e2}k-ofXEatq<*ybqiOnFrY?@-+<HoQ`v
zrdO>u`IcrUIbS~eE-j}K)TL!9fRhb3-}^!@YE~M$3wV_(C6$Zlk5b@8suna<)om2#
zhSTe+zK(RZ2L}!W5nb2_t|K%qAkAT2dm$wQQulhepc%NuIlt!Z#1G7mXvax+lW6kt
z%PN;#sqD@~v{Dn${hmZiYBG@%6<O|%?M4wJL(T{yClaQsZ`!0X2f9y^BbhdzT!B(%
zPV=+poiMg(eaR`A=t+f>c<0DvVrSyYh%{RXX57k{9+z}^*0nW|Ol%pEx?*SZp+`pD
z^=TAb-VdY325BH0RmC-5*R09bR3>qWketJSNaxwvh<19R_~x}*LgoaWI%m+*8A^M`
zUviw%y(80_CZ^P?xR<rE*Qfmrcj*>Sca|ZiP%AUTFq=r{NZU;^7lk+@GXD_No^(18
zklVox4&n03#PR0cN|VVoS2lwrGbs*=V{YUXk^^n!(-oxUoUIPE$<YaCQ#SsM!^0Ew
z6*t@y`fhEu;_V5{?56utP72na)_Hp|&3vWfwat$0>FLJf*C>G=dX0E8<<sZKPL?Vz
zqlhU}n{DMoPN@KP^lk{Iq`dDDN$JL$p0nh{5cxNH(+z1WGEGZO&jdk5W$L;2JY(8G
z+*k8TUrv`PE1q6cj5#rLW7Ik|bFlM9O75AY$w!US@Nm>%oSap1$EA;%B6LY5X^eY)
z@8HIa=y!Q;MXlbqe*cKcYqY^3#b&XQ72Kk^NsA}BbgG?c`jyi3OmHZ7FbFFQ^<}r}
zR~GV2NZjn{x4{B8PV4ob83=T%B#0_pSD>-CRkAsspsC66D@UiVn%=!_a_7YOuSqn0
zqP`B{6@R9nzstQqnja~aDuy2IOuo{Z*InMKunB)Abi)NYvOgnLYh-9<(q#IUn=bB9
z)~MW2VUX|UmabD=LEQ`81~-=dxy~SN13M2b^vCUZky6AX&bA#B+jfjkZ@qTM*q&`u
zI}(v!({xLu+`47%olTRsLu7M7B-44p7cG0HQzy6clDB3gZEw=62H}Ffxu$NCmbS$!
zyW9)3-+T`bS9Rp+NN&Y9#0j*=IZ+ynkSqOKvC=G|ZA7$I&#jJR4*5a%F_Wq3G5o>{
zbM<;AdSBdprXwE>#ycdzeXZ2e;Id!Qr}hfd9U(r^tCjq6=KL{t#;v(yL9G!6Xe42h
ze!1k8)#Ui@tF~?%pV+c{dVI&|=85s?vC++=J4Q3-r?orXx^2~IsDXhHS<ZV`^o+eV
zFqcyC>jbe~zhY7~(T^WH@rM^tN&Mr+P6`~?Sf8niPJxL2e-?WE7>NFlx~pIJ^7}7A
zbpKhn2F5|+Z3g?Y60j1m60j1m60j1m60j1m60j1m60j1m5_pwKKulS&Axj(95<^TG
z%WDO*pjlco2Cq`Et69t10Xwzw5pDL(tDss=O}|Qh8FR4Mo!Tr<#*myED>k;9bn<Gb
z%`nz=HAnRSi_qt{i%yUJ&yuu$eVpHa2RFuH&$6aJ#P>Xiz2Lw2eJkvQemDc3WIg`@
z@Sy--V(tFD@O!WeR>Fnw3vezx!+QPa;bs_ybKxj!_K(6}!Ctrk{+adq--jDv9L|L|
zz&UUxe4TasharNs@FR>EUxa&L27U?to%Q<L;2p3BE`^KWJorBA^k0Ct!6rBd{u@C&
z1R?wie1moQyWmzR!DTQ2YoH%;@LKpP>+}!6d*Qd>5bT08;3#YM55hmfeehv;C)^HF
zt=iOD6=N#v(H-VgF`i20FFGrYr7Kxp>}pY+hOZa5MRi^6WJ=;6-x$}{(6@Lq%~=G4
zkz5<xjX~R=;ZSObj!n(ERc<(2%q#hr*HTNmq+erk-E>M5C5|DB<W<+NY4-sWQk-)`
zwR*LVpv}P)%=BUNPJJaWy&ljArp=OotVPQ*Mo?olD{3v^C|3EmV(!cJRrI0{OINLL
z>9UHfCnPJpvQ}E=43>tS<XTI8w%dpoakFU8h8Z~GHf9o)RwhsGFpDAaxD)DSH9h_m
z@Jwi}i-`M1lFKB6l9}XKx+0byXj=<mkwPBvk&bD)`pm`12jg2I^#ovHMX!cU1Vy*9
zQ=A{hS%=#wD=2&Lamg#a#3NrvWU)$K#vYy?C&iLitccr7uKluE*9GZmWD75n&djw7
zQ%51m_M>K0H<v$|_qJ*oX{DkCqBGkmk1;(IOt`ZlVlm!3*%@RPD@n0$#Z-t3Nm4oO
zi+HUTN%HHH){u##t=u{RBE)O_i;k4cw)aG&B8{|p7xm(1t5Xzgdl&h1zFrRPeQ4S3
z-I46S;0Av~vYu$i_;_)?b4+`Y$M|tP<)U~->oveHZyGQ|GV=qu|H&ErQ~#}{2LHIX
zu*`hwJ43c-q$BC`z31s#wzs!EOW*dK`YzEW{?^TmG@RwXUWzHAJ$J{wgU=g>pY#U5
z{nS?AQpdHU(92aT5dHt7=*PF9H;ewCkM;li(D9Ez6RrcX5jgN^^!z`B8MqwIgpZ-?
z-wo5S4xU4||1LZYAA%ZG-~eog95`?f`u+l}hYR2_bp7{(3m3ze(Cr_GN8us(I6Me9
z!3J0fkCSi7`<r1W?0`vl6Ff_~{ub_r8(<s!B77a4|9-d=ejToc)8L!v{dd4Sp$>at
z1w4%Ye=`)|JotI|0lNMd;cmDSeh03Hvq60KZ-)r3hGXdWpM(a;%QO@-yJ+=}`Ne8t
zuS!kqgp+ZbCT7&!L=mRdq{DqWml~}^tJy8D-y~EGt0K;;k@2O~4ba3VDba>w9bU$9
z8Yz1diWOAOwIqb|xC&}<%Qfvc`@Sc2_FQE9N3uB}ulCG^p4ajh#^R9)?mHZB`^cu0
ze9xpfuXlXp)f8^!jo9;;l+5A2jHiBbZt|H)10#YL`a#&0N)nKylesC<qCzr@nRJTS
zbVBw(upOo&J84LglFo(p!_;qhT^>%YAA}_b+2Vq)rjz7(k$C1%>$t9{jNoLr&^Ceb
zN_CVk+DH>^>=3rG`?YySchzXQ9ACIsdf5VYCdN*lG@oh7bUsbBq(0BjC*gG>fH%wA
zWkzH=;iRV0q;-o;9lzD$IyIoyORhS<WD-ehX$cv&ml@$i4e<+&JxOCf<qkoUls2gJ
zikeN2c)B-J3+Nxsy0M6<gMMR9iK99DG|FL6)!Tk>UDSKH^q?sFfs$E6x06Zwsy22F
z4==<1KbHDHj;K|#u7EjH@8i*?<CdeYq#Ezw$@)xpZ7?h7JZxlhRlLDW@x4!xWgXBJ
z+I1j_P6=CQ@U2&p?NX3kQGF7xf^AT5lhSGF9a2b7Po`{<5p_R`Bp#c(*i@w>8QWUy
z(l0x#wBU7W@eQbnS4?fU6se@Q$z&Y9jf*R(Cng0`f=rC*c?^|l2IHGYV;VL=@zpZb
zpv}TvN$9-yI6o14zRf$lQcQA`9cHqoj7l*ldxKW19a-0XyqRUQ%#iG;Gn2tSHioQL
z>j@Z&_q3TZ9{2opz0OZ|&n@Hp!CvPNX0ESiGsXbfwxYN9>Fq0RMUgEn*?lasVNn~A
z%>Ws1{W8y#W&mG$7fPjSX2@jOplw$pZEJZi>N$4UjXH&{CxRsP<7^E;<bBfA<844a
z_kOiDyo`pSyU=YKUf+53WXpjRvd8$GYi-9!@}GB3oY~Bn5hu`AIm_N~suk6?-KP4;
z3)yPwX>-tVwwg)=<{SN&8C>PU(k`XPJ^!7i(puQ&Mt6?(RlKrvF1bds!IXX$2ZsLa
zTxV!MlN$RDOyADjoUlX8k%TjcWE)uT>+E0h8vB<@-{gf|T!>(Z?W6nT3ppPmsia$H
zuC2?}BX-<@neu+Tn_{KrBIL1jAJ<j8hzyT=8{8(~L>8o5mAhp(Ixt-d4%TuKKy~I7
zg5tBk;LZ%qNR?$IWP&)gbMsAnG9mizAahF6N_ua>)kW+%d!NkyThafIplgc`F8Y6;
z)=<%H(e?it=3yVm{{Lst^=0q>`{3Ph19)&9Tn&Rz0NMZl9qa%ffVaYLz$Ne$I{!m(
zE6l+#oDRph_Ooz5)S(7dkUaInX<%Pg0#*W60#*W60#*W60#*W60;fa*soU8J9-x$o
z*JTmwy5Y_|M4+{Mc-b0o%YJSjxumr-`^crN!%FJzD*EBe>5)rfj}W7mekxMB*3q^=
zleI_M2Plca^h;_oHt0W(2Pm1xM`2+}9x>FfMw~e1&g3Vil4n&JQ-tXMt63|53Pk^3
z9S^EsM&Ew~J`3_3z<0w9a0#r0Z=&OW3T^}O1sH$};S1>b55ppOAkX<b9ln8{|2cRE
zjKFi~^moA@!yiEeyJ00<2*<EV{2SZ?qae@n`!?(S55qjj`oGTeDVYRpB31%c0#*W6
z0#*W60zbJDFgHb7r?7jV<bjE@j`ix?17$t$vKm6>d&b+%dfv(7p{?g#>^x~b@3PN)
zd?`Hd3Wbr@!`#POetpIhEMrBG*)3O%X0;V2<5{LShT<PZ1w7Kk&sX&S$C>b_@KYE4
zKl%T9&!OZ00KNnN3Ll33@FqBlp8p8^4cr8>?*9xr{#{UnYe1guFFOBM;VbY2d<<@b
zE8#SF0zLl$_!Bq)qW6CpJ^wwh8$|D40Z*Xo%RYd|;B)X1coz)Acd=W18t#CD&<{`H
z<Ntn`hDo>pzJ-7P*WgL`1pGOO-~VgiF?{^*fm`4$@I(Cj?}Y&N!x}h-U;p>uf8Z!Q
z4UfVH;U*};7+eA8!6V31o&zAy?!Orppb1*eFWn869_RJ_kIekE1mpiT=}Evz96!nE
zT+2Ka);`9-KE^<|!6n|>w^`c#eCY-c?Y*4Z)U(S0_7*<tW@jH`U?1CI{_TzS+GU1o
zA5vaPdsJCVYkO?R@>NzxhS&DBHV)|fz;qkyiFc%CH*4ER)#!h)^5cC}jdAf^dVpWv
zG+@kK_Av(1Ct9x;)ETlpqmlee5c?Pd{s9Ef=Xz0}-gN4_#4GwS1}8m5bUxhU{{ctZ
BgmVA@

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 471da90..45c0c00 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -45,7 +45,7 @@ struct TestCoords {
 };
 
 template <typename Result, typename CoordField, typename AnalyticFunction>
-bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
+bool checkResult(Omega_h::Mesh &mesh, Result &result, CoordField coordField,
                  TestCoords testCase, AnalyticFunction func, size_t numComp) {
   const auto numPtsPerElem = testCase.NumPtsPerElem;
   MeshField::FieldElement fcoords(
@@ -64,7 +64,8 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
           const auto x = globalCoords(pt, 0);
           const auto y = globalCoords(pt, 1);
           const auto expected = func(x, y);
-            const auto computed = result(pt, 0);
+          for (int i = 0; i < numComp; ++i) {
+            const auto computed = result(pt, i);
             MeshField::LO isError = 0;
             if (Kokkos::fabs(computed - expected) >
                 MeshField::MachinePrecision) {
@@ -76,6 +77,7 @@ bool checkResult(Omega_h::Mesh &mesh, Result result, CoordField coordField,
             }
             lerrors += isError;
           }
+        }
       },
       numErrors);
   return (numErrors > 0);
@@ -229,7 +231,7 @@ void doRun(Omega_h::Mesh &mesh,
 
     {
       const auto ShapeOrder = 1;
-      const auto numComponents = 1;
+      const auto numComponents = 2;
       auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
                                                     numComponents>();
       LinearFunction func = LinearFunction();
@@ -245,7 +247,7 @@ void doRun(Omega_h::Mesh &mesh,
     }
     {
       const auto ShapeOrder = 1;
-      const auto numComponents = 1;
+      const auto numComponents = 3;
       auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
                                                     numComponents>();
       LinearFunction func = LinearFunction();

From cbc524985c3619f5b2f1747b8a3abceae6752d7e Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Fri, 27 Jun 2025 14:17:18 -0400
Subject: [PATCH 13/14] Fixed Shape Field

---
 src/MeshField_ShapeField.hpp | 2 +-
 test/testOmegahElement.cpp   | 4 +---
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp
index bfc3fb1..523aad1 100644
--- a/src/MeshField_ShapeField.hpp
+++ b/src/MeshField_ShapeField.hpp
@@ -197,7 +197,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
     Ctrlr kk_ctrl = createController(meshInfo.numVtx);
 #else
     using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
-    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1});
+    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp});
 #endif
     auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
     using LA = LinearAccessor<decltype(vtxField)>;
diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 45c0c00..0f39aa8 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -115,9 +115,7 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) {
     const auto x = (coords[left * MeshDim] + coords[right * MeshDim]) / 2.0;
     const auto y =
         (coords[left * MeshDim + 1] + coords[right * MeshDim + 1]) / 2.0;
-    for (int i = 0; i < field.numComp; ++i) {
-      field(edge, 0, 0, MeshField::Edge) = func(x, y);
-    }
+    field(edge, 0, 0, MeshField::Edge) = func(x, y);
   };
   MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()},
                           setFieldAtEdges, "setFieldAtEdges");

From 0c5d459bd763eeab7f0b17892c3e053d1c80bd4f Mon Sep 17 00:00:00 2001
From: Kloepfer <kloepj2@checkers.scorec.rpi.edu>
Date: Mon, 30 Jun 2025 12:19:28 -0400
Subject: [PATCH 14/14] added tests and created helper function

---
 test/testOmegahElement.cpp | 113 ++++++++++++++++---------------------
 1 file changed, 50 insertions(+), 63 deletions(-)

diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp
index 0f39aa8..8e055ed 100644
--- a/test/testOmegahElement.cpp
+++ b/test/testOmegahElement.cpp
@@ -115,7 +115,9 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) {
     const auto x = (coords[left * MeshDim] + coords[right * MeshDim]) / 2.0;
     const auto y =
         (coords[left * MeshDim + 1] + coords[right * MeshDim + 1]) / 2.0;
-    field(edge, 0, 0, MeshField::Edge) = func(x, y);
+    for (int i = 0; i < field.numComp; ++i) {
+      field(edge, 0, i, MeshField::Edge) = func(x, y);
+    }
   };
   MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()},
                           setFieldAtEdges, "setFieldAtEdges");
@@ -146,6 +148,35 @@ void doFail(std::string_view order, std::string_view function,
   std::string msg = ss.str();
   MeshField::fail(msg);
 }
+template <size_t ShapeOrder, size_t numComponents,
+          template <typename...> typename Controller>
+void runTest(Omega_h::Mesh &mesh,
+             MeshField::OmegahMeshField<ExecutionSpace, 2, Controller> &omf,
+             auto testCase, auto function) {
+  using functionType = decltype(function);
+  using ViewType = decltype(testCase.coords);
+  auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
+                                                numComponents>();
+  using FieldType = decltype(field);
+  setVertices(mesh, function, field);
+  if constexpr (ShapeOrder == 2) {
+    setEdges(mesh, function, field);
+  }
+  auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
+      testCase.coords, testCase.NumPtsPerElem, field);
+  auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
+                            decltype(function){}, numComponents);
+  if (failed) {
+    std::string fieldErr = ShapeOrder == 1 ? "linear" : "quadratic";
+    std::string functionErr;
+    if constexpr (std::is_same_v<functionType, LinearFunction>) {
+      functionErr = "linear";
+    } else {
+      functionErr = "quadratic";
+    }
+    doFail(fieldErr, functionErr, testCase.name, std::to_string(numComponents));
+  }
+}
 
 template <template <typename...> typename Controller>
 void doRun(Omega_h::Mesh &mesh,
@@ -177,87 +208,43 @@ void doRun(Omega_h::Mesh &mesh,
     {
       const auto ShapeOrder = 1;
       const auto numComponents = 1;
-      auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
-                                                    numComponents>();
-      LinearFunction func = LinearFunction();
-      setVertices(mesh, func, field);
-      using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
-          testCase.coords, testCase.NumPtsPerElem, field);
-      auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{}, numComponents);
-      if (failed)
-        doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents));
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase, LinearFunction());
     }
 
     {
       const auto ShapeOrder = 2;
       const auto numComponents = 1;
-      auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
-                                                    numComponents>();
-      auto func = QuadraticFunction();
-      setVertices(mesh, func, field);
-      setEdges(mesh, func, field);
-      using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
-          testCase.coords, testCase.NumPtsPerElem, field);
-      auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                QuadraticFunction{}, numComponents);
-      if (failed)
-        doFail("quadratic", "quadratic", testCase.name,
-               std::to_string(numComponents));
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase,
+                                         QuadraticFunction());
     }
 
     {
       const auto ShapeOrder = 2;
       const auto numComponents = 1;
-      auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
-                                                    numComponents>();
-      auto func = LinearFunction();
-      setVertices(mesh, func, field);
-      setEdges(mesh, func, field);
-      using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
-          testCase.coords, testCase.NumPtsPerElem, field);
-      auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{}, numComponents);
-      if (failed)
-        doFail("quadratic", "linear", testCase.name,
-               std::to_string(numComponents));
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase, LinearFunction());
     }
 
     {
       const auto ShapeOrder = 1;
       const auto numComponents = 2;
-      auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
-                                                    numComponents>();
-      LinearFunction func = LinearFunction();
-      setVertices(mesh, func, field);
-      using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
-          testCase.coords, testCase.NumPtsPerElem, field);
-      auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{}, numComponents);
-      if (failed)
-        doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents));
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase, LinearFunction());
     }
     {
       const auto ShapeOrder = 1;
       const auto numComponents = 3;
-      auto field = omf.template CreateLagrangeField<MeshField::Real, ShapeOrder,
-                                                    numComponents>();
-      LinearFunction func = LinearFunction();
-      setVertices(mesh, func, field);
-      using FieldType = decltype(field);
-      auto result = omf.template triangleLocalPointEval<ViewType, FieldType>(
-          testCase.coords, testCase.NumPtsPerElem, field);
-      auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
-                                LinearFunction{}, numComponents);
-      if (failed)
-        doFail("linear", "linear", testCase.name,
-               std::to_string(numComponents));
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase, LinearFunction());
+    }
+    {
+      const auto ShapeOrder = 2;
+      const auto numComponents = 2;
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase,
+                                         QuadraticFunction());
+    }
+
+    {
+      const auto ShapeOrder = 2;
+      const auto numComponents = 3;
+      runTest<ShapeOrder, numComponents>(mesh, omf, testCase, LinearFunction());
     }
   }
 }