Fixes to Integration

CMakeLists.txt

@@ -4,6 +4,7 @@ project(meshfields VERSION 0.1.0 LANGUAGES CXX)
 
 include(GNUInstallDirs)
 include(CMakePackageConfigHelpers)
+include(CheckGit.cmake)
 
 option(MeshFields_USE_Cabana "Build with the Cabana storage backend" OFF)
 
@@ -67,6 +68,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)
@@ -172,12 +176,16 @@ meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
+meshfields_add_exe(IntegratorTest test/testIntegrator.cpp)
+meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
   meshfields_add_exe(CabanaTests test/testCabana.cpp)
+  meshfields_add_exe(IntegratorPerformance test/testIntegratorPerformance.cpp)
   test_func(CabanaTests ./CabanaTests)
   test_func(ControllerPerformance ./ControllerPerformance)
+  test_func(IntegratorPerformance ./IntegratorPerformance)
 endif()
 
 test_func(KokkosTests ./KokkosTests)
@@ -190,6 +198,9 @@ test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
 test_func(OmegahTetTest, ./OmegahTetTest)
+test_func(IntegratorTest ./IntegratorTest)
+test_func(BoxIntegratorTest ./BoxIntegratorTest)
+
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)

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,76 @@
+\documentclass{article}
+\usepackage{graphicx}
+\usepackage[letterpaper,textwidth=5in,right=2.5in,textheight=9in]{geometry}
+
+\usepackage{amsmath,amssymb, changepage}
+\usepackage{enumerate}
+\usepackage{upgreek}
+\pagestyle{empty}
+\usepackage{listings}
+\usepackage{tikz}
+\begin{document}
+\noindent
+Joshua Kloepfer\\
+9/30/2025\\
+\textbf{Meshfields DOF Holder and Integration Point Ordering}\\\\
+Linear Triangle:\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0) -- (2.5, 2.5);
+    \draw[gray, thick] (0, 0) -- (5, 0);
+    \draw[gray, thick] (2.5, 2.5) -- (5, 0);
+    \filldraw[black] (0, 0) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (2.5, 2.5) circle (2pt) node[anchor=south]{V2};
+\end{tikzpicture}
+\\
+Quadratic Triangle\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0) -- (2.5, 2.5);
+    \draw[gray, thick] (0, 0) -- (5, 0);
+    \draw[gray, thick] (2.5, 2.5) -- (5, 0);
+    \filldraw[black] (0, 0) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (2.5, 2.5) circle (2pt) node[anchor=south]{V2};
+    \filldraw[black] (2.5, 0) circle (2pt) node[anchor=north]{E0};
+    \filldraw[black] (3.75, 1.25) circle (2pt) node[anchor=west]{E1};
+    \filldraw[black] (1.25, 1.25) circle (2pt) node[anchor=east]{E2};
+\end{tikzpicture}
+\\
+Linear Tetrahedron:\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0, 0) -- (5, 0, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (5, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (5, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (0, 5, 0) -- (0, 0, 5);
+    \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=east]{V2};
+    \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V3};
+\end{tikzpicture}
+\\
+Vertex ordering must follow right hand rule.\\
+Quadratic Tetrahedron:\\
+\begin{tikzpicture}
+    \draw[gray, thick] (0, 0, 0) -- (5, 0, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (5, 0, 0) -- (0, 5, 0);
+    \draw[gray, thick] (0, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (5, 0, 0) -- (0, 0, 5);
+    \draw[gray, thick] (0, 5, 0) -- (0, 0, 5);
+    \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V0};
+    \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1};
+    \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=east]{V2};
+    \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V3};
+
+    \filldraw[black] (2.5, 0, 2.5) circle (2pt) node[anchor=north]{E0};
+    \filldraw[black] (0, 0, 2.5) circle (2pt) node[anchor=east]{E1};
+    \filldraw[black] (0, 2.5, 2.5) circle (2pt) node[anchor=east]{E2};
+    \filldraw[black] (2.5, 0, 0) circle (2pt) node[anchor=south]{E3};
+    \filldraw[black] (2.5, 2.5, 0) circle (2pt) node[anchor=east]{E4};
+    \filldraw[black] (0, 2.5, 0) circle (2pt) node[anchor=west]{E5};
+\end{tikzpicture}
+\\
+Vertex ordering must follow right hand rule and edge ordering must follow diagram.
+\end{document}

mfemTest/IntegratorMfemTest.cpp

@@ -0,0 +1,58 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <chrono>
+
+using namespace std;
+using namespace mfem;
+
+double intFunc(const Vector &x) {
+  return x[0];
+}
+
+int main(int argc, char *argv[]) {
+
+  OptionsParser args(argc, argv);
+  int order = 1;
+  int size = 1;
+  int numberOfElements = 1;
+  int type = 1;
+  args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+  args.AddOption(&size, "-s", "--size", "Size of mesh");
+  args.AddOption(&numberOfElements, "-n", "--numElem", "number of elements in each direction");
+  args.AddOption(&type, "-t", "--type", "type of element(1 for tet, 0 for triangle");
+  args.Parse();
+
+  if (type == 1) {
+    Mesh mesh = Mesh::MakeCartesian3D(numberOfElements, numberOfElements, numberOfElements, Element::TETRAHEDRON, size, size, size);
+    H1_FECollection fec(order, mesh.Dimension());
+    FiniteElementSpace fespace(&mesh, &fec);
+
+    FunctionCoefficient custom(intFunc);
+
+    LinearForm b(&fespace);
+
+    b.AddDomainIntegrator(new DomainLFIntegrator(custom)); 
+    auto start = std::chrono::high_resolution_clock::now();
+    b.Assemble(); 
+    auto end = std::chrono::high_resolution_clock::now();
+    cout << "RESULT: " << mfem::GetGitStr() << "," << mesh.GetNE() << "," << order << "," << b.Sum() << "," << std::chrono::duration<double>(end - start).count() << endl;
+
+  }
+  else { 
+    Mesh mesh = Mesh::MakeCartesian2D(numberOfElements, numberOfElements, Element::TRIANGLE, size, size);
+    H1_FECollection fec(order, mesh.Dimension());
+    FiniteElementSpace fespace(&mesh, &fec);
+
+    FunctionCoefficient custom(intFunc);
+
+    LinearForm b(&fespace);
+
+    b.AddDomainIntegrator(new DomainLFIntegrator(custom)); 
+    auto start = std::chrono::high_resolution_clock::now();
+    b.Assemble(); 
+    auto end = std::chrono::high_resolution_clock::now();
+    cout << "RESULT: " << mfem::GetGitStr() << "," << mesh.GetNE() << "," << order << "," << b.Sum() << "," << std::chrono::duration<double>(end - start).count() << endl;
+  }
+  return 0;
+}

mfemTest/README.md

@@ -0,0 +1,7 @@
+Run using g++ IntegratorMfemTest.cpp  -I<path-to-mfem-build> -L<path-to-mfem-build>  -l mfem
+in order to run, use -m to specify the mesh file, -o to specify order, -s to specify size, -n to specify number of elements in each direction, and -t to specify type of element with 1 for tetrahedron and 0 for triangle.
+mfem version: https://github.com/mfem/mfem/tree/2caa75e35a54df93d19f23655170254556dfc081
+For the build:
+mkdir <mfem-build-dir> ; cd <mfem-build-dir>
+cmake <mfem-source-dir>
+make -j 4

src/MeshField_Element.hpp

@@ -303,6 +303,30 @@ struct FieldElement {
     return Kokkos::View<Real *>("foo", J.extent(0));
   }
 
+  template <size_t MeshEntDim>
+  KOKKOS_INLINE_FUNCTION auto getGradients(Kokkos::View<Real **> lc,
+                                           size_t pt) const {
+    if constexpr (ShapeType::Order == 1) {
+      return shapeFn.getLocalGradients();
+    } else {
+      Kokkos::Array<Real, MeshEntDim + 1> coord;
+      for (int i = 0; i < MeshEntDim + 1; ++i) {
+        coord[i] = lc(pt, i);
+      }
+      return shapeFn.getLocalGradients(coord);
+    }
+  }
+
+  KOKKOS_INLINE_FUNCTION auto setNodalGradients(
+      Kokkos::View<Real * [ShapeType::numNodes][MeshEntDim]> nodalGradients,
+      auto grad, size_t pt, size_t node, size_t d) const {
+    if constexpr (ShapeType::Order == 1) {
+      nodalGradients(pt, node, d) = grad[node * MeshEntDim + d];
+    } else {
+      nodalGradients(pt, node, d) = grad[node][d];
+    }
+  }
+
   /**
    * @brief
    * Given an array of parametric coordinates 'localCoords', one per mesh
@@ -380,22 +404,22 @@ struct FieldElement {
       // one matrix per point
       Kokkos::View<Real ***> res("result", numPts, MeshEntDim, MeshEntDim);
       Kokkos::deep_copy(res, 0.0); // initialize all entries to zero
-
       // fill the views of node coordinates and node gradients
       Kokkos::View<Real * [ShapeType::numNodes][MeshEntDim]> nodeCoords(
           "nodeCoords", numPts);
       Kokkos::View<Real * [ShapeType::numNodes][MeshEntDim]> nodalGradients(
           "nodalGradients", numPts);
-      const auto grad = shapeFn.getLocalGradients();
       Kokkos::parallel_for(
           numMeshEnts, KOKKOS_CLASS_LAMBDA(const int ent) {
             const auto vals = getNodeValues(ent);
             assert(vals.size() == MeshEntDim * ShapeType::numNodes);
             for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) {
+              auto ptCoords = Kokkos::subview(localCoords, pt, Kokkos::ALL());
+              const auto grad = getGradients<MeshEntDim>(localCoords, pt);
               for (size_t node = 0; node < ShapeType::numNodes; node++) {
                 for (size_t d = 0; d < MeshEntDim; d++) {
                   nodeCoords(pt, node, d) = vals[node * MeshEntDim + d];
-                  nodalGradients(pt, node, d) = grad[node * MeshEntDim + d];
+                  setNodalGradients(nodalGradients, grad, pt, node, d);
                 }
               }
             }

src/MeshField_GetHash.cpp.in

@@ -0,0 +1,2 @@
+#include "MeshField_GetHash.hpp"
+const char * GitHash = "@GIT_HASH@";

src/MeshField_GetHash.hpp

@@ -0,0 +1,6 @@
+#ifndef GETHASH
+#define GETHASH
+
+extern const char * GitHash;
+
+#endif

src/MeshField_Integrate.hpp

@@ -92,16 +92,16 @@ class TetrahedronIntegration : public EntityIntegration<4> {
     virtual std::vector<IntegrationPoint<4>> getPoints() const {
 
       return {IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
-                                       0.138196601125011, 0.585410196624969},
+                                       0.138196601125011, 0.585410196624967},
                                0.25 / 6.0),
-              IntegrationPoint(Vector4{0.585410196624969, 0.138196601125011,
+              IntegrationPoint(Vector4{0.585410196624967, 0.138196601125011,
                                        0.138196601125011, 0.138196601125011},
                                0.25 / 6.0),
-              IntegrationPoint(Vector4{0.138196601125011, 0.585410196624969,
+              IntegrationPoint(Vector4{0.138196601125011, 0.585410196624967,
                                        0.138196601125011, 0.138196601125011},
                                0.25 / 6.0),
               IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
-                                       0.585410196624969, 0.138196601125011},
+                                       0.585410196624967, 0.138196601125011},
                                0.25 / 6.0)};
     }
     virtual int getAccuracy() const { return 2; }