Fixes to Integration

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(IntegratorTest test/testIntegrator.cpp)
+meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -190,6 +192,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)

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/IntegratorMfemLargeTest.cpp

@@ -0,0 +1,27 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double intFunc(const Vector &x) {
+  return x[0];
+}
+
+int main(int argc, char *argv[]) {
+  Mesh mesh = Mesh::MakeCartesian3D(100, 100, 100, Element::TETRAHEDRON, 100, 100, 100, false);
+  //6 million element mesh
+  H1_FECollection fec(1, mesh.Dimension());
+  FiniteElementSpace fespace(&mesh, &fec);
+
+  FunctionCoefficient custom(intFunc);
+
+  LinearForm b(&fespace);
+
+  b.AddDomainIntegrator(new DomainLFIntegrator(custom)); 
+
+  b.Assemble(); 
+  cout << b.Sum() << endl;
+  return 0;
+}

mfemTest/IntegratorMfemTest.cpp

@@ -0,0 +1,37 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double intFunc(const Vector &x) {
+  return x[0]; 
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../data/star.mesh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+   Mesh mesh(mesh_file);
+
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   FunctionCoefficient custom(intFunc);
+
+   LinearForm b(&fespace);
+
+   b.AddDomainIntegrator(new DomainLFIntegrator(custom)); 
+
+   b.Assemble(); 
+
+   cout << b.Sum() << 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 and -o to specify order
+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

mfemTest/box.mesh

@@ -0,0 +1,53 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+6
+1 4 7 0 3 1
+1 4 7 0 1 5
+1 4 7 0 5 4
+1 4 7 0 2 3
+1 4 7 0 6 2
+1 4 7 0 4 6
+
+boundary
+12
+1 2 3 0 2
+1 2 0 3 1
+6 2 7 4 5
+6 2 4 7 6
+5 2 6 0 4
+5 2 0 6 2
+3 2 7 1 3
+3 2 1 7 5
+2 2 5 0 1
+2 2 0 5 4
+4 2 7 2 6
+4 2 2 7 3
+
+vertices
+8
+3
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+0 0 1
+1 0 1
+0 1 1
+1 1 1

mfemTest/square.mesh

@@ -0,0 +1,36 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+2
+1 2 0 1 3
+1 2 2 3 1
+
+boundary
+4
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 0
+
+vertices
+4
+2
+0 0
+1 0
+1 1
+0 1

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_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; }