Merge pull request #65 from SCOREC/jk/linearTet

cwsmith · web-flow · commit fc4ef7d328aa · 2025-08-05T09:36:18.000-04:00
add quadratic and linear tet shape functions
diff --git a/CMakeLists.txt b/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(OmegahTetTest test/testOmegahTet.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(OmegahTetTest, ./OmegahTetTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)
diff --git a/src/MeshField.hpp b/src/MeshField.hpp
@@ -89,6 +89,37 @@ struct LinearTriangleToVertexField {
     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 QuadraticTriangleToField {
   }
 };
 
+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, 4, 5, 0, 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 @@ template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
                   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 @@ class OmegahMeshField {
     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
diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp
@@ -105,6 +105,34 @@ 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));
+    // clang-format off
+    return {1 - xi[0] - xi[1] - xi[2], 
+            xi[0], xi[1], 
+            xi[2]};
+    // clang-format on
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, meshEntDim * numNodes> getLocalGradients() const {
+    // clang-format off
+    return {-1, -1, -1, 
+             1,  0,  0, 
+             0,  1,  0, 
+             0,  0,  1};
+    // clang-format on
+  }
+};
+
 struct QuadraticTriangleShape {
   static const size_t numNodes = 6;
   static const size_t meshEntDim = 2;
@@ -151,7 +179,8 @@ struct QuadraticTetrahedronShape {
   constexpr static size_t NumDofHolders[2] = {4, 6};
   constexpr static size_t DofsPerHolder[2] = {1, 1};
   constexpr static size_t Order = 2;
-
+  // ordering taken from mfem
+  // see mfem/mfem fem/fe/fe_fixed_order.cpp @597cba8
   KOKKOS_INLINE_FUNCTION
   Kokkos::Array<Real, numNodes> getValues(Vector4 const &xi) const {
     assert(greaterThanOrEqualZero(xi));
@@ -163,9 +192,9 @@ struct QuadraticTetrahedronShape {
             xi[1]*(2*xi[1]-1),
             xi[2]*(2*xi[2]-1),
             4*xi[0]*xi3,
-            4*xi[0]*xi[1],
             4*xi[1]*xi3,
             4*xi[2]*xi3,
+            4*xi[0]*xi[1],
             4*xi[2]*xi[0],
             4*xi[1]*xi[2]};
     // clang-format on
@@ -183,9 +212,9 @@ struct QuadraticTetrahedronShape {
             0,4*xi[1]-1,0,
             0,0,4*xi[2]-1,
             4*xi3-4*xi[0],-4*xi[0],-4*xi[0],
-            4*xi[1],4*xi[0],0,
             -4*xi[1],4*xi3-4*xi[1],-4*xi[1],
             -4*xi[2],-4*xi[2],4*xi3-4*xi[2],
+            4*xi[1],4*xi[0],0,
             4*xi[2],0,4*xi[0],
             0,4*xi[2],4*xi[1]};
     // clang-format on
diff --git a/src/MeshField_ShapeField.hpp b/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,12 @@ 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>;
+    // clang-format off
+    using LinearLagrangeShapeField = std::conditional_t<
+        dim == 3, 
+        ShapeField<numComp, Ctrlr, LinearTetrahedronShape, LA>,
+        ShapeField<numComp, Ctrlr, LinearTriangleShape, LA>>;
+    // clang-format on
     LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField});
     return llsf;
   } else if constexpr (order == 2 && (dim == 2 || dim == 3)) {
@@ -242,8 +246,12 @@ 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>;
+    // clang-format off
+    using QuadraticLagrangeShapeField = std::conditional_t<
+        dim == 3, 
+        ShapeField<numComp, Ctrlr, QuadraticTetrahedronShape, QA>,
+        ShapeField<numComp, Ctrlr, QuadraticTriangleShape, QA>>;
+    // clang-format on
     QuadraticLagrangeShapeField qlsf(kk_ctrl, meshInfo, {vtxField, edgeField});
     return qlsf;
   } else {
diff --git a/test/testOmegahTet.cpp b/test/testOmegahTet.cpp
