Jk/linear tet

CMakeLists.txt

@@ -170,6 +170,7 @@ meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
 meshfields_add_exe(OmegahElementTests test/testOmegahElement.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
+meshfields_add_exe(TetTest test/testTet.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -186,6 +187,7 @@ test_func(ElementJacobian2d ./ElementJacobian2d)
 test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahElementTests ./OmegahElementTests)
 test_func(PointMapping ./PointMapping)
+test_func(TetTest, ./TetTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)

src/MeshField.hpp

@@ -89,6 +89,37 @@
     return {0, triCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo
   }
 };
+struct LinearTetrahedronToVertexField {
+  Omega_h::LOs tetVerts;
+  LinearTetrahedronToVertexField(Omega_h::Mesh &mesh)
+      : tetVerts(mesh.ask_elem_verts()) {
+    if (mesh.dim() != 3 && mesh.family() != OMEGA_H_SIMPLEX) {
+      MeshField::fail(
+          "The mesh passed to %s must be 3D and simplex (tetrahedron)\n",
+          __func__);
+    }
+  }
+  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() const {
+    return {MeshField::Tetrahedron};
+  }
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx,
+             MeshField::LO tet, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Tetrahedron);
+    const auto tetDim = 3;
+    const auto vtxDim = 0;
+    const auto ignored = -1;
+    const auto localVtxIdx =
+        (Omega_h::simplex_down_template(tetDim, vtxDim, tetNodeIdx, ignored) +
+         3) %
+        4;
+    const auto tetToVtxDegree = Omega_h::simplex_degree(tetDim, vtxDim);
+    const MeshField::LO vtx = tetVerts[(tet * tetToVtxDegree) + localVtxIdx];
+    return {0, tetCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo
+  }
+};
 struct QuadraticTriangleToField {
   Omega_h::LOs triVerts;
   Omega_h::LOs triEdges;
@@ -151,6 +182,59 @@
   }
 };
 
+struct QuadraticTetrahedronToField {
+  Omega_h::LOs tetVerts;
+  Omega_h::LOs tetEdges;
+  QuadraticTetrahedronToField(Omega_h::Mesh &mesh)
+      : tetVerts(mesh.ask_elem_verts()),
+        tetEdges(mesh.ask_down(mesh.dim(), 1).ab2b) {
+    if (mesh.dim() != 3 && mesh.family() != OMEGA_H_SIMPLEX) {
+      MeshField::fail(
+          "The mesh passed to %s must be 3D and simplex (tetrahedron)",
+          __func__);
+    }
+  }
+
+  KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() const {
+    return {MeshField::Tetrahedron};
+  }
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx,
+             MeshField::LO tet, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Tetrahedron);
+    const MeshField::LO tetNode2DofHolder[10] = {0, 1, 2, 3, 3, 0, 4, 5, 2, 1};
+    const MeshField::Mesh_Topology tetNode2DofHolderTopo[10] = {
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Edge,   MeshField::Edge,
+        MeshField::Edge,   MeshField::Edge,   MeshField::Edge,
+        MeshField::Edge};
+    const auto dofHolderIdx = tetNode2DofHolder[tetNodeIdx];
+    const auto dofHolderTopo = tetNode2DofHolderTopo[tetNodeIdx];
+    Omega_h::LO osh_ent;
+    if (dofHolderTopo == MeshField::Vertex) {
+      const auto tetDim = 3;
+      const auto vtxDim = 0;
+      const auto ignored = -1;
+      const auto localVtxIdx = (Omega_h::simplex_down_template(
+                                    tetDim, vtxDim, dofHolderIdx, ignored) +
+                                3) %
+                               4;
+      const auto tetToVtxDegree = Omega_h::simplex_degree(tetDim, vtxDim);
+      osh_ent = tetVerts[(tet * tetToVtxDegree) + localVtxIdx];
+    } else if (dofHolderTopo == MeshField::Edge) {
+      const auto tetDim = 3;
+      const auto edgeDim = 1;
+      const auto tetToEdgeDegree = Omega_h::simplex_degree(tetDim, edgeDim);
+      osh_ent = tetEdges[(tet * tetToEdgeDegree) + dofHolderIdx];
+    } else {
+      assert(false);
+    }
+    return {0, tetCompIdx, osh_ent, dofHolderTopo};
+  }
+};
+
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
   static_assert(ShapeOrder == 1 || ShapeOrder == 2);
   if constexpr (ShapeOrder == 1) {
@@ -169,6 +253,24 @@
                   QuadraticTriangleToField(mesh)};
   }
 }
+template <int ShapeOrder> auto getTetrahedronElement(Omega_h::Mesh &mesh) {
+  static_assert(ShapeOrder == 1 || ShapeOrder == 2);
+  if constexpr (ShapeOrder == 1) {
+    struct result {
+      MeshField::LinearTetrahedronShape shp;
+      LinearTetrahedronToVertexField map;
+    };
+    return result{MeshField::LinearTetrahedronShape(),
+                  LinearTetrahedronToVertexField(mesh)};
+  } else if constexpr (ShapeOrder == 2) {
+    struct result {
+      MeshField::QuadraticTetrahedronShape shp;
+      QuadraticTetrahedronToField map;
+    };
+    return result{MeshField::QuadraticTetrahedronShape(),
+                  QuadraticTetrahedronToField(mesh)};
+  }
+}
 
 } // namespace Omegah
 
@@ -252,6 +354,32 @@
     auto eval = MeshField::evaluate(f, localCoords, offsets);
     return eval;
   }
+
+  template <typename ViewType, typename ShapeField>
+  auto tetrahedronLocalPointEval(ViewType localCoords, size_t NumPtsPerElem,
+                                 ShapeField field) {
+    auto offsets = createOffsets(meshInfo.numTet, NumPtsPerElem);
+    auto eval = tetrahedronLocalPointEval(localCoords, offsets, field);
+    return eval;
+  }
+
+  template <typename ViewType, typename ShapeField>
+  auto tetrahedronLocalPointEval(ViewType localCoords,
+                                 Kokkos::View<LO *> offsets, ShapeField field) {
+    const auto MeshDim = 3;
+    if (mesh.dim() != MeshDim) {
+      MeshField::fail("input mesh must be 3d\n");
+    }
+    const auto ShapeOrder = ShapeField::Order;
+    if (ShapeOrder != 1 && ShapeOrder != 2) {
+      MeshField::fail("input field order must be 1 or 2\n");
+    }
+    const auto [shp, map] = Omegah::getTetrahedronElement<ShapeOrder>(mesh);
+    MeshField::FieldElement<ShapeField, decltype(shp), decltype(map)> f(
+        meshInfo.numTet, field, shp, map);
+    auto eval = MeshField::evaluate(f, localCoords, offsets);
+    return eval;
+  }
 };
 
 } // namespace MeshField

src/MeshField_Shape.hpp

@@ -105,6 +105,25 @@ struct LinearTriangleCoordinateShape {
   }
 };
 
+struct LinearTetrahedronShape {
+  static const size_t numNodes = 4;
+  static const size_t meshEntDim = 3;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 1;
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector4 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+    return {1 - xi[0] - xi[1] - xi[2], xi[0], xi[1], xi[2]};
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, meshEntDim * numNodes> getLocalGradients() const {
+    return {-1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1};
+  }
+};
+
 struct QuadraticTriangleShape {
   static const size_t numNodes = 6;
   static const size_t meshEntDim = 2;

src/MeshField_ShapeField.hpp

@@ -172,7 +172,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
   static_assert((dim == 1 || dim == 2 || dim == 3),
                 "CreateLagrangeField only supports 1d, 2d, and 3d meshes\n");
   using MemorySpace = typename ExecutionSpace::memory_space;
-  if constexpr (order == 1 && (dim == 1 || dim == 2)) {
+  if constexpr (order == 1 && (dim == 1 || dim == 2 || dim == 3)) {
     if (meshInfo.numVtx <= 0) {
       fail("mesh has no vertices\n");
     }
@@ -201,8 +201,9 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
 #endif
     auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
     using LA = LinearAccessor<decltype(vtxField)>;
-    using LinearLagrangeShapeField =
-        ShapeField<numComp, Ctrlr, LinearTriangleShape, LA>;
+    using LinearLagrangeShapeField = std::conditional_t<
+        dim == 3, ShapeField<numComp, Ctrlr, LinearTetrahedronShape, LA>,
+        ShapeField<numComp, Ctrlr, LinearTriangleShape, LA>>;
     LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField});
     return llsf;
   } else if constexpr (order == 2 && (dim == 2 || dim == 3)) {
@@ -242,8 +243,9 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
     auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
     auto edgeField = MeshField::makeField<Ctrlr, 1>(kk_ctrl);
     using QA = QuadraticAccessor<decltype(vtxField), decltype(edgeField)>;
-    using QuadraticLagrangeShapeField =
-        ShapeField<numComp, Ctrlr, QuadraticTriangleShape, QA>;
+    using QuadraticLagrangeShapeField = std::conditional_t<
+        dim == 3, ShapeField<numComp, Ctrlr, QuadraticTetrahedronShape, QA>,
+        ShapeField<numComp, Ctrlr, QuadraticTriangleShape, QA>>;
     QuadraticLagrangeShapeField qlsf(kk_ctrl, meshInfo, {vtxField, edgeField});
     return qlsf;
   } else {

test/testTet.cpp

@@ -0,0 +1,139 @@
+#include <Kokkos_Core.hpp>
+#include <MeshField.hpp>
+#include <Omega_h_build.hpp>
+#include <Omega_h_mesh.hpp>
+template <size_t order> void runTest(Omega_h::Mesh mesh) {
+  MeshField::OmegahMeshField<Kokkos::DefaultExecutionSpace, 3> mesh_field(mesh);
+  auto shape_field =
+      mesh_field.template CreateLagrangeField<double, order, 1>();
+  auto dim = mesh.dim();
+  auto coords = mesh.coords();
+  auto f = KOKKOS_LAMBDA(double x, double y, double z)->double {
+    return 2 * x + 3 * y + 4 * z;
+  };
+  auto edge2vtx = mesh.ask_down(1, 0).ab2b;
+  auto edgeMap = mesh.ask_down(dim, 1).ab2b;
+  Kokkos::parallel_for(
+      mesh.nverts(), KOKKOS_LAMBDA(int vtx) {
+        auto x = coords[vtx * dim];
+        auto y = coords[vtx * dim + 1];
+        auto z = coords[vtx * dim + 2];
+        shape_field(vtx, 0, 0, MeshField::Mesh_Topology::Vertex) = f(x, y, z);
+      });
+  if (order == 2) {
+    Kokkos::parallel_for(
+        mesh.nedges(), KOKKOS_LAMBDA(int edge) {
+          const auto left = edge2vtx[edge * 2];
+          const auto right = edge2vtx[edge * 2 + 1];
+          const auto x = (coords[left * dim] + coords[right * dim]) / 2.0;
+          const auto y =
+              (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0;
+          const auto z =
+              (coords[left * dim + 2] + coords[right * dim + 2]) / 2.0;
+          shape_field(edge, 0, 0, MeshField::Mesh_Topology::Edge) = f(x, y, z);
+        });
+  }
+  const auto numNodesPerElem = order == 2 ? 10 : 4;
+  Kokkos::View<double **> local_coords("", mesh.nelems() * numNodesPerElem, 4);
+  Kokkos::parallel_for(
+      "set", mesh.nelems() * numNodesPerElem, KOKKOS_LAMBDA(const int &i) {
+        const auto val = i % numNodesPerElem;
+        local_coords(i, 0) = (val == 0);
+        local_coords(i, 1) = (val == 1);
+        local_coords(i, 2) = (val == 2);
+        local_coords(i, 3) = (val == 3);
+        if constexpr (order == 2) {
+          if (val == 4) {
+            local_coords(i, 0) = 1 / 2.0;
+            local_coords(i, 1) = 1 / 2.0;
+          } else if (val == 5) {
+            local_coords(i, 1) = 1 / 2.0;
+            local_coords(i, 2) = 1 / 2.0;
+          } else if (val == 6) {
+            local_coords(i, 0) = 1 / 2.0;
+            local_coords(i, 2) = 1 / 2.0;
+          } else if (val == 7) {
+            local_coords(i, 0) = 1 / 2.0;
+            local_coords(i, 3) = 1 / 2.0;
+          } else if (val == 8) {
+            local_coords(i, 1) = 1 / 2.0;
+            local_coords(i, 3) = 1 / 2.0;
+          } else if (val == 9) {
+            local_coords(i, 2) = 1 / 2.0;
+            local_coords(i, 3) = 1 / 2.0;
+          }
+        }
+      });
+  auto eval_results = mesh_field.tetrahedronLocalPointEval(
+      local_coords, numNodesPerElem, shape_field);
+
+  int errors = 0;
+  const auto tetVerts = mesh.ask_elem_verts();
+  const auto coordField = mesh_field.getCoordField();
+  Kokkos::parallel_reduce(
+      "test", mesh.nelems(),
+      KOKKOS_LAMBDA(const int &tet, int &errors) {
+        for (int node = 0; node < 4; ++node) {
+          const auto tetDim = 3;
+          const auto vtxDim = 0;
+          const auto ignored = -1;
+          const auto localVtxIdx =
+              Omega_h::simplex_down_template(tetDim, vtxDim, node, ignored);
+          const auto tetToVtxDegree = Omega_h::simplex_degree(tetDim, vtxDim);
+          int vtx = tetVerts[(tet * tetToVtxDegree) + localVtxIdx];
+          auto x = coords[vtx * dim];
+          auto y = coords[vtx * dim + 1];
+          auto z = coords[vtx * dim + 2];
+          auto expected = f(x, y, z);
+          auto actual = eval_results(tet * numNodesPerElem + node, 0);
+          if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) {
+            ++errors;
+
+            Kokkos::printf(
+                "expected: %lf, actual: %lf, element: %d, node(vtx): %d\n",
+                expected, actual, tet, node);
+          }
+        }
+        for (int node = 4; node < numNodesPerElem; ++node) {
+          const auto tetDim = 3;
+          const auto edgeDim = 1;
+          const auto tetToEdgeDegree = Omega_h::simplex_degree(tetDim, edgeDim);
+          const MeshField::LO tetNode2DofHolder[6] = {0, 1, 2, 3, 4, 5};
+          int edge =
+              edgeMap[(tet * tetToEdgeDegree) + tetNode2DofHolder[node - 4]];
+          auto left = edge2vtx[edge * 2];
+          auto right = edge2vtx[edge * 2 + 1];
+          const auto x = (coords[left * dim] + coords[right * dim]) / 2.0;
+          const auto y =
+              (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0;
+          const auto z =
+              (coords[left * dim + 2] + coords[right * dim + 2]) / 2.0;
+          auto expected = f(x, y, z);
+          auto actual = eval_results(tet * numNodesPerElem + node, 0);
+          if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) {
+            ++errors;
+            Kokkos::printf(
+                "expected: %lf, actual: %lf, element: %d, node(edge): %d\n",
+                expected, actual, tet, node);
+          }
+        }
+      },
+      errors);
+  if (errors > 0) {
+    MeshField::fail("One or more mappings did not match\n");
+  }
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  auto lib = Omega_h::Library(&argc, &argv);
+  {
+    auto world = lib.world();
+    auto mesh =
+        Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 10, false);
+    runTest<1>(mesh);
+    runTest<2>(mesh);
+  }
+  Kokkos::finalize();
+  return 0;
+}