Fixes to Integration

CMakeLists.txt

@@ -4,6 +4,7 @@ project(meshfields VERSION 1.0.1 LANGUAGES CXX)
 
 include(GNUInstallDirs)
 include(CMakePackageConfigHelpers)
+include(CheckGit.cmake)
 
 option(MeshFields_USE_Cabana "Build with the Cabana storage backend" OFF)
 
@@ -76,6 +77,9 @@ configure_file(
   "${CMAKE_CURRENT_SOURCE_DIR}/src/MeshField_Config.hpp.in"
   "${CMAKE_CURRENT_BINARY_DIR}/MeshField_Config.hpp")
 
+
+CheckGitSetup()
+target_sources(meshfields PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/MeshField_GetHash.cpp)
 if(MeshFields_USE_Cabana)
   target_link_libraries(meshfields INTERFACE Cabana::Core)
   target_compile_definitions(meshfields INTERFACE MESHFIELDS_ENABLE_CABANA)
@@ -212,6 +216,9 @@ meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
 meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
+meshfields_add_exe(IntegratorTest test/testIntegrator.cpp)
+meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp)
+
 meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp)
 meshfields_add_exe(Interpolate test/testInterpolate.cpp)
 meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp)
@@ -220,6 +227,7 @@ if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
   meshfields_add_smoke_exe(CabanaTests test/testCabana.cpp)
   meshfields_add_exe(CabanaControllerLifetimeTest test/testCabanaControllerLifetime.cpp)
+  meshfields_add_exe(IntegratorPerformance test/testIntegratorPerformance.cpp)
   smoke_test_func(CabanaTests ./CabanaTests)
   test_func(ControllerPerformance ./ControllerPerformance)
   test_func(CabanaControllerLifetimeTest ./CabanaControllerLifetimeTest)
@@ -233,6 +241,8 @@ test_func(ElementJacobian1d ./ElementJacobian1d)
 test_func(ElementJacobian2d ./ElementJacobian2d)
 test_func(ElementJacobian3d ./ElementJacobian3d)
 test_func(CountIntegrator ./CountIntegrator)
+test_func(IntegratorTest ./IntegratorTest)
+test_func(BoxIntegratorTest ./BoxIntegratorTest)
 smoke_test_func(OmegahTriTests ./OmegahTriTests)
 test_func(OmegahTetTest ./OmegahTetTest)
 test_func(IntegratorMeasure ./IntegratorMeasure)
@@ -244,7 +254,7 @@ else()
   will_fail_test_func(ExceptionTest ./ExceptionTest)
 endif()
 
-meshfields_add_exe(OmegahSprAdapt test/testSprThwaitesAdapt.cpp)
+meshfields_add_exe(OmegahSprAdapt test/testSprThwaitesAdapt.cpp) 
 test_func(OmegahSprAdapt ./OmegahSprAdapt
   ${CMAKE_SOURCE_DIR}/meshes/thwaites_basal_effectiveStrain.osh/
   thwaitesAdapted #output prefix

CheckGit.cmake

@@ -0,0 +1,61 @@
+set(pre_configure_file ${CMAKE_CURRENT_SOURCE_DIR}/src/MeshField_GetHash.cpp.in)
+set(post_configure_file ${CMAKE_CURRENT_BINARY_DIR}/MeshField_GetHash.cpp)
+
+function(CheckGitWrite git_hash)
+    file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/git-state.txt "${git_hash}\n")
+endfunction()
+
+function(CheckGitRead git_hash)
+    if (EXISTS ${CMAKE_CURRENT_BINARY_DIR}/git-state.txt)
+        file(STRINGS ${CMAKE_CURRENT_BINARY_DIR}/git-state.txt CONTENT)
+        LIST(GET CONTENT 0 var)
+
+        set(${git_hash} ${var} PARENT_SCOPE)
+    endif ()
+endfunction()
+
+function(CheckGitVersion)
+    execute_process(COMMAND git rev-parse HEAD
+      WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}
+      OUTPUT_VARIABLE GIT_HASH
+      OUTPUT_STRIP_TRAILING_WHITESPACE
+    )
+    CheckGitRead(GIT_HASH_CACHE)
+    if (NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR})
+        file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+    endif ()
+
+    if (NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/MeshField_GetHash.hpp)
+        file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/MeshField_GetHash.hpp DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+    endif()
+
+    if (NOT DEFINED GIT_HASH_CACHE)
+        set(GIT_HASH_CACHE "INVALID")
+    endif ()
+
+    if (NOT ${GIT_HASH} STREQUAL ${GIT_HASH_CACHE} OR NOT EXISTS ${post_configure_file})
+        CheckGitWrite(${GIT_HASH})
+
+        configure_file(${pre_configure_file} ${post_configure_file} @ONLY)
+    endif ()
+
+endfunction()
+
+function(CheckGitSetup)
+
+    add_custom_target(AlwaysCheckGit COMMAND ${CMAKE_COMMAND}
+        -DRUN_CHECK_GIT_VERSION=1
+        -Dpre_configure_dir=${CMAKE_CURRENT_SOURCE_DIR}
+        -Dpost_configure_file=${CMAKE_CURRENT_BINARY_DIR}
+        -DGIT_HASH_CACHE=${GIT_HASH_CACHE}
+        -P ${CMAKE_CURRENT_SOURCE_DIR}/CheckGit.cmake
+        BYPRODUCTS ${post_configure_file}
+        )
+
+
+    CheckGitVersion()
+endfunction()
+
+if (RUN_CHECK_GIT_VERSION)
+  CheckGitVersion()
+endif ()

docs/Ordering.tex

@@ -0,0 +1,87 @@
+\documentclass{article}
+\usepackage[letterpaper,textwidth=5in,right=2.5in,textheight=9in]{geometry}
+\usepackage{tikz}
+\pagestyle{empty}
+\begin{document}
+\noindent
+Joshua Kloepfer\\
+9/30/2025\\
+\textbf{Meshfields DOF Holder and Integration Point Ordering}\\\\
+The node ordering and parametric coordinates below follow
+Zienkiewicz, Taylor, and Zhu, ``The Finite Element Method: Its Basis and
+Fundamentals'', 2013, and match the shape functions defined in
+\texttt{src/MeshField\_Shape.hpp}.\\
+\\
+Linear Triangle:\\
+\begin{tikzpicture}[scale=0.8]
+    \draw[gray, thick] (0, 0) -- (0, 5);
+    \draw[gray, thick] (0, 0) -- (5, 0);
+    \draw[gray, thick] (0, 5) -- (5, 0);
+    \draw[->] (-1, 0) -- (-1, 2) node[anchor=east]{$\xi_1$};
+    \draw[->] (0, -1) -- (2, -1) node[anchor=north]{$\xi_0$};
+    \filldraw[black] (0, 0) circle (2pt) node[anchor=north east]{V0 $(0,0)$};
+    \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1 $(1,0)$};
+    \filldraw[black] (0, 5) circle (2pt) node[anchor=south]{V2 $(0,1)$};
+\end{tikzpicture}
+\\
+Quadratic Triangle\\
+\begin{tikzpicture}[scale=0.8]
+    \draw[gray, thick] (0, 0) -- (0, 5);
+    \draw[gray, thick] (0, 0) -- (5, 0);
+    \draw[gray, thick] (0, 5) -- (5, 0);
+    \draw[->] (-1, 0) -- (-1, 2) node[anchor=east]{$\xi_1$};
+    \draw[->] (0, -1) -- (2, -1) node[anchor=north]{$\xi_0$};
+    \filldraw[black] (0, 0) circle (2pt) node[anchor=north east]{V0 $(0,0)$};
+    \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1 $(1,0)$};
+    \filldraw[black] (0, 5) circle (2pt) node[anchor=south]{V2 $(0,1)$};
+    \filldraw[black] (2.5, 0) circle (2pt) node[anchor=north]{E0};
+    \filldraw[black] (2.5, 2.5) circle (2pt) node[anchor=south west]{E1};
+    \filldraw[black] (0, 2.5) circle (2pt) node[anchor=east]{E2};
+\end{tikzpicture}
+\\
+Edge nodes: E0 on V0--V1, E1 on V1--V2, E2 on V0--V2;
+these are nodes 3, 4 and 5 of \texttt{QuadraticTriangleShape}.\\
+Linear Tetrahedron:\\
+\begin{tikzpicture}[scale=0.8]
+    \draw[gray, thick] (0, 0, 0) -- (5, 0, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (5, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (5, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (0, 5, 0) -- (0, 0, 5);
+    \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=north east]{V0 $(0,0,0)$};
+    \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1 $(1,0,0)$};
+    \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V2 $(0,1,0)$};
+    \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V3 $(0,0,1)$};
+\end{tikzpicture}
+\\
+Vertex node $i$ is at $L_i = 1$: V0 at the origin, V1 on $\xi_0$,
+V2 on $\xi_1$, V3 on $\xi_2$.  Vertex ordering must follow the right hand
+rule, so that the Jacobian determinant is positive.\\
+Quadratic Tetrahedron:\\
+\begin{tikzpicture}[scale=0.8]
+    \draw[gray, thick] (0, 0, 0) -- (5, 0, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (5, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (5, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (0, 5, 0) -- (0, 0, 5);
+    \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=north east]{V0};
+    \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V2};
+    \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V3};
+
+    \filldraw[black] (2.5, 0, 0) circle (2pt) node[anchor=north]{E0};
+    \filldraw[black] (0, 2.5, 0) circle (2pt) node[anchor=west]{E1};
+    \filldraw[black] (0, 0, 2.5) circle (2pt) node[anchor=east]{E2};
+    \filldraw[black] (2.5, 2.5, 0) circle (2pt) node[anchor=west]{E3};
+    \filldraw[black] (0, 2.5, 2.5) circle (2pt) node[anchor=west]{E4};
+    \filldraw[black] (2.5, 0, 2.5) circle (2pt) node[anchor=north]{E5};
+\end{tikzpicture}
+\\
+Edge nodes, corresponding to nodes 4--9 of
+\texttt{QuadraticTetrahedronShape}: E0 on V0--V1, E1 on V0--V2,
+E2 on V0--V3, E3 on V1--V2, E4 on V2--V3, E5 on V1--V3.
+Vertex ordering must follow the right hand rule and edge ordering must
+follow the diagram.
+\end{document}

docs/addingShapeFunctionsForOmegah.md

@@ -52,6 +52,9 @@ omitted:
 | Triangle (2D) | \f$(\xi_0,\xi_1)\f$  | \f$[0,1]^2\f$, \f$\xi_0+\xi_1 \le 1\f$ |
 | Tet (3D)      | \f$(\xi_0,\xi_1,\xi_2)\f$ | \f$[0,1]^3\f$, \f$\sum \xi_i \le 1\f$ |
 
+Diagrams of the vertex and edge node ordering, with parametric coordinates, for
+the linear and quadratic triangle and tetrahedron are in `docs/Ordering.tex`.
+
 ---
 
 ## Step 1 - Define the Shape Function Struct {#step1}
@@ -107,10 +110,8 @@ The mapping struct lives in `src/MeshField.hpp` inside
 (as numbered by the shape function) into the on-process mesh entity indices
 (as numbered by Omega\_h) that hold the corresponding DOFs.
 
-### Vertex-Ordering Correction
+### Mapping Omega\_h Vertices and Edges to Nodes
 
-As previously stated, Omega\_h numbers vertices and edges within an element 
-differently from the MeshFields canonical ordering used by the shape functions.
 The Omega\_h canonical orderings for simplices are depicted below:
 
 <a href="omegah_tri.png"><img src="omegah_tri.png" width="20%" alt="Omega_h triangle ordering"/></a>
@@ -124,27 +125,50 @@ Omega_h tetrahedron ordering.
 Meshes composed of hypercubes (quadrilaterals and hexahedrons), pyramids, and prisms are supported by Omega_h, but they cannot be adapted.
 Version 1.1.0 of MeshFields only supports Omega\_h simplices.
 
-The existing linear-triangle and linear-tetrahedron mappings correct for
-the difference in Omega\_h vs MeshFields ordering with a cyclic rotation.
+**Vertices.** The MeshFields vertex node ordering for triangles and tetrahedra
+matches the Omega\_h element-to-vertex ordering, so vertex node \f$i\f$ maps
+directly to local vertex \f$i\f$ with no rotation:
 
 ```cpp
-// For triangles (triDim=2, vtxDim=0):
+// For triangles (triDim=2) and tetrahedra (tetDim=3), vtxDim=0:
 const auto localVtxIdx =
-    (Omega_h::simplex_down_template(triDim, vtxDim, nodeIdx, /*ignored=*/-1) + 2) % 3;
+    Omega_h::simplex_down_template(elemDim, vtxDim, nodeIdx, /*ignored=*/-1);
+```
 
-// For tetrahedra (tetDim=3, vtxDim=0):
-const auto localVtxIdx =
-    (Omega_h::simplex_down_template(tetDim, vtxDim, nodeIdx, /*ignored=*/-1) + 3) % 4;
+**Edges.** Edge nodes do not necessarily follow the Omega\_h element-to-edge
+ordering, so the mapping struct uses a lookup table from node index to the
+Omega\_h local edge index.
+For the quadratic triangle the MeshFields edges (E0: V0-V1, E1: V1-V2,
+E2: V0-V2) coincide with the Omega\_h edges, so the table is the identity.
+For the quadratic tetrahedron the MeshFields edges (nodes 4-9: V0-V1, V0-V2,
+V0-V3, V1-V2, V2-V3, V1-V3) are reordered to the Omega\_h edges
+(0: V0-V1, 1: V1-V2, 2: V2-V0, 3: V0-V3, 4: V1-V3, 5: V2-V3):
+
+```cpp
+// from QuadraticTetrahedronToField
+const MeshField::LO tetNode2DofHolder[10] = {0, 1, 2, 3, 0, 2, 3, 1, 5, 4};
 ```
 
-You must determine the correct rotation offset for your topology and DOF holders by
-comparing the Omega\_h canonical vertex/edge ordering (see
+You must determine the correct vertex and edge mapping for your topology and
+DOF holders by comparing the Omega\_h canonical vertex/edge ordering (see
 figures above based on `Omega_h_simplex.hpp` and the `simplex_down_template` function)
 with the node ordering established by `getNodeParametricCoords()` in your shape
 struct.  It is a good idea to add a unit test that evaluates the shape functions at
 the parametric coordinates of each node and checks that the value for
 that node is 1.0 and all others are 0.0.
 
+### Jacobian Orientation
+
+@ref MeshField::FieldElement "FieldElement" returns **signed** Jacobian
+determinants for 2D and 3D elements.  The integrator
+(`src/MeshField_Integrate.hpp`) uses the determinants as-is, and debug builds
+assert that each one is non-negative.
+The vertex ordering produced by your mapping must therefore follow the
+right-hand rule so that the Jacobian determinant of a valid (un-inverted)
+Omega\_h element is positive.
+An incorrect vertex permutation results in a negative determinant and a
+failed assertion in debug builds.
+
 ### Interface Requirements
 
 | Member | Description |
@@ -225,7 +249,7 @@ The steps within that branch are:
 2. **Allocate storage** -- construct a `Controller` sized for each DOF-holding entity type.  Each field in the controller needs entries for `numEntities * dofsPerHolder * numComp`.
 3. **Create fields** -- call `MeshField::makeField<Ctrlr, N>(ctrl)` for each field index `N` to obtain typed slice objects.
 4. **Construct the Accessor** -- pass the slice objects to your Accessor's constructor (e.g., `QuadraticAccessor{vtxField, edgeField}`).
-5. **Return** `FieldWithController<Ctrlr, ShapeField<numComp, YourShape, YourAccessor>>{ctrl, ShapeField(meshInfo, accessor)}`.
+5. **Return** `FieldWithController<Ctrlr, ShapeField<numComp, YourShape, Ctrlr, YourAccessor>>{ctrl, ShapeField(meshInfo, accessor)}`.
 
 See the existing `order == 2 && dim == 2` branch (quadratic triangle) as a
 concrete reference:

pumiTest/CMakeLists.txt

@@ -0,0 +1,17 @@
+cmake_minimum_required(VERSION 3.16)
+project(PumiTest LANGUAGES CXX C)
+
+find_package(MPI REQUIRED)
+
+find_package(SCOREC REQUIRED)
+
+add_executable(PumiTest main.cpp)
+
+target_link_libraries(PumiTest PRIVATE
+    MPI::MPI_CXX
+    SCOREC::apf
+    SCOREC::mds
+    SCOREC::gmi
+    SCOREC::pcu
+    SCOREC::lion
+)

pumiTest/main.cpp

@@ -0,0 +1,56 @@
+#include <apf.h>
+#include <apfMDS.h>
+#include <apfMesh2.h>
+#include <apfBox.h>
+#include <pcu_util.h>
+#include <pumi_version.h>
+#include <PCU.h>
+#include <iostream>
+#include <gmi_mesh.h>
+#include <gmi_null.h>
+#include <chrono>
+
+class CountIntegrator : public apf::Integrator {
+  protected:
+    double numEnt;
+    apf::MeshElement* currentMeshElement;
+  public:
+    double getCount() {return numEnt;}
+    void resetCount() { numEnt = 0; }
+    CountIntegrator(size_t order) : Integrator(order), numEnt(0) {};
+    void inElement(apf::MeshElement* me) override {
+      currentMeshElement = me;
+    }
+    void atPoint(apf::Vector3 const& point, double weight, double jacobian) override {
+      apf::Vector3 globalCoord;
+      apf::mapLocalToGlobal(currentMeshElement, point, globalCoord);
+      numEnt += weight * globalCoord[0] * jacobian;
+    }
+};
+int main(int argc, char ** argv) {
+  pcu::Init(&argc, &argv);
+  {
+  pcu::PCU PCUObj;
+  size_t order = atoi(argv[1]);
+  size_t size = atoi(argv[2]);
+  size_t numberOfElems = atoi(argv[3]);
+  size_t dim = atoi(argv[4]);
+  size_t runs = atoi(argv[5]);
+  gmi_register_mesh();
+  gmi_register_null();
+  apf::Mesh2* mesh = apf::makeMdsBox(numberOfElems, numberOfElems, dim == 3 ? numberOfElems : 0, size, size, dim == 3 ? size : 0, 1, &PCUObj);
+  CountIntegrator * countInt = new CountIntegrator(order);
+  for (int i = 0; i < runs; ++i) {
+    countInt->resetCount();
+    auto start = std::chrono::high_resolution_clock::now();
+    countInt->process(mesh);
+    auto end = std::chrono::high_resolution_clock::now();
+    std::cout << "RESULT: " << pumi_git_sha() << "," << mesh->count(mesh->getDimension()) << "," << order << "," << countInt->getCount() << "," << std::chrono::duration<double>(end - start).count() << std::endl;
+  }
+  delete countInt;
+  mesh->destroyNative();
+  apf::destroyMesh(mesh);
+  }
+  pcu::Finalize();
+  return 0;
+}

src/KokkosController.hpp

@@ -62,7 +62,7 @@ template <class SliceType, class T> struct KokkosSliceWrapper {
   }
 };
 
-template <class MemorySpace, class ExecutionSpace, typename... Ts>
+template <class MemorySpace, class ExeSpace, typename... Ts>
 class KokkosController {
 
   // type definitions
@@ -182,7 +182,8 @@ class KokkosController {
   std::tuple<Kokkos::View<Ts, MemorySpace>...> values_;
 
 public:
-  typedef ExecutionSpace exe;
+  typedef ExeSpace exe;
+  using ExecutionSpace=exe;
   static const int MAX_RANK = 5;
 
   KokkosController() {