add quadratic tet

src/MeshField_Shape.hpp

@@ -13,7 +13,10 @@ namespace MeshField {
 struct LinearEdgeShape {
   KOKKOS_INLINE_FUNCTION
   Kokkos::Array<Real, 2> getValues(Kokkos::Array<Real, 2> const &xi) const {
-    return {(1.0 - xi[0]) / 2.0, (1.0 + xi[0]) / 2.0};
+    // clang-format off
+    return {(1.0 - xi[0]) / 2.0,
+            (1.0 + xi[0]) / 2.0};
+    // clang-format on
   }
   static const size_t numNodes = 2;
   static const size_t numComponentsPerDof = 1;
@@ -24,7 +27,11 @@ struct LinearEdgeShape {
 struct LinearTriangleShape {
   KOKKOS_INLINE_FUNCTION
   Kokkos::Array<Real, 3> getValues(Kokkos::Array<Real, 3> const &xi) const {
-    return {1 - xi[0] - xi[1], xi[0], xi[1]};
+    // clang-format off
+    return {1 - xi[0] - xi[1],
+            xi[0],
+            xi[1]};
+    // clang-format on
   }
   static const size_t numNodes = 3;
   static const size_t numComponentsPerDof = 1;
@@ -36,9 +43,14 @@ struct QuadraticTriangleShape {
   KOKKOS_INLINE_FUNCTION
   Kokkos::Array<Real, 6> getValues(Kokkos::Array<Real, 3> const &xi) const {
     const Real xi2 = 1 - xi[0] - xi[1];
-    return {xi2 * (2 * xi2 - 1),     xi[0] * (2 * xi[0] - 1),
-            xi[1] * (2 * xi[1] - 1), 4 * xi[0] * xi2,
-            4 * xi[0] * xi[1],       4 * xi[1] * xi2};
+    // clang-format off
+    return {xi2 * (2 * xi2 - 1),
+            xi[0] * (2 * xi[0] - 1),
+            xi[1] * (2 * xi[1] - 1),
+            4 * xi[0] * xi2,
+            4 * xi[0] * xi[1],
+            4 * xi[1] * xi2};
+    // clang-format on
   }
 
   static const size_t numNodes = 6;
@@ -48,5 +60,32 @@ struct QuadraticTriangleShape {
   constexpr static size_t NumDofHolders[2] = {3, 3};
   constexpr static size_t DofsPerHolder[2] = {1, 1};
 };
+
+struct QuadraticTetrahedronShape {
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, 10> getValues(Kokkos::Array<Real, 4> const &xi) const {
+    const Real xi3 = 1 - xi[0] - xi[1] - xi[2];
+    // clang-format off
+    return {xi3*(2*xi3-1),
+            xi[0]*(2*xi[0]-1),
+            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[2]*xi[0],
+            4*xi[1]*xi[2]};
+    // clang-format on
+  }
+
+  static const size_t numNodes = 10;
+  static const size_t numComponentsPerDof = 1;
+  static const size_t meshEntDim = 3;
+  constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
+  constexpr static size_t NumDofHolders[2] = {4, 6};
+  constexpr static size_t DofsPerHolder[2] = {1, 1};
+};
+
 } // namespace MeshField
 #endif

src/MeshField_ShapeField.hpp

@@ -74,8 +74,8 @@ auto CreateLagrangeField(MeshInfo &meshInfo) {
   static_assert(
       (order == 1 || order == 2),
       "CreateLagrangeField only supports linear and quadratic fields\n");
-  static_assert((dim == 1 || dim == 2),
-                "CreateLagrangeField only supports 1d and 2d meshes\n");
+  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)) {
     assert(meshInfo.numVtx > 0);
@@ -90,7 +90,7 @@ auto CreateLagrangeField(MeshInfo &meshInfo) {
         ShapeField<MeshField<Ctrlr>, LinearTriangleShape, LA>;
     LinearLagrangeShapeField llsf(kokkosMeshField, meshInfo, {vtxField});
     return llsf;
-  } else if constexpr (order == 2 && dim == 2) {
+  } else if constexpr (order == 2 && (dim == 2 || dim == 3)) {
     assert(meshInfo.numVtx > 0);
     assert(meshInfo.numEdge > 0);
     using Ctrlr = Controller::KokkosController<MemorySpace, ExecutionSpace,

test/testElement.cpp

@@ -105,14 +105,16 @@ struct QuadraticTriangleToField {
              MeshField::LO ent, MeshField::Mesh_Topology topo) const {
     assert(topo == MeshField::Triangle);
     assert(ent == 0);
+    // clang-format off
     // hardcoded using the following numbering in MeshFields
     //         2
     //       /   \
     //      5     4
     //    /        \
     //   0 --  3 -- 1
     MeshField::LO triNode2DofHolder[6] = {/*vertices*/ 0, 1, 2,
-                                          /*edges*/ 0,    1, 2};
+                                          /*edges*/ 0, 1, 2};
+    // clang-format on
     MeshField::Mesh_Topology triNode2DofHolderTopo[6] = {
         /*vertices*/ MeshField::Vertex,
         MeshField::Vertex,
@@ -126,7 +128,7 @@ struct QuadraticTriangleToField {
   }
 };
 
-// evaluate a field at the specified local coordinate for each triangle using
+// evaluate a field at the specified local coordinate for one triangle using
 // quadratic shape functions
 void quadraticTriangleLocalPointEval() {
   MeshField::MeshInfo meshInfo;
@@ -147,6 +149,58 @@ void quadraticTriangleLocalPointEval() {
   auto x = MeshField::evaluate(f, lc);
 }
 
+struct QuadraticTetrahedronToField {
+  constexpr static MeshField::Mesh_Topology Topology[1] = {
+      MeshField::Tetrahedron};
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx,
+             MeshField::LO ent, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Tetrahedron);
+    assert(ent == 0);
+    // clang-format off
+    MeshField::LO tetNode2DofHolder[10] = {
+        /*vertices*/ 0, 1, 2, 3,
+        /*edges*/ 0, 1, 2, 3, 4, 5};
+    // clang-format on
+    MeshField::Mesh_Topology tetNode2DofHolderTopo[10] = {
+        /*vertices*/ MeshField::Vertex,
+        MeshField::Vertex,
+        MeshField::Vertex,
+        /*edges*/ MeshField::Edge,
+        MeshField::Edge,
+        MeshField::Edge,
+        MeshField::Edge,
+        MeshField::Edge,
+        MeshField::Edge};
+    const auto dofHolder = tetNode2DofHolder[tetNodeIdx];
+    const auto dofHolderTopo = tetNode2DofHolderTopo[tetNodeIdx];
+    return {0, 0, dofHolder, dofHolderTopo};
+  }
+};
+
+// evaluate a field at the specified local coordinate for one tet using
+// linear shape functions
+void quadraticTetrahedronLocalPointEval() {
+  MeshField::MeshInfo meshInfo;
+  meshInfo.numVtx = 4;
+  meshInfo.numEdge = 6;
+  meshInfo.numTri = 4;
+  meshInfo.numTet = 1;
+  auto field =
+      MeshField::CreateLagrangeField<ExecutionSpace, MeshField::Real, 2, 3>(
+          meshInfo);
+
+  MeshField::Element elm{MeshField::QuadraticTetrahedronShape(),
+                         QuadraticTetrahedronToField()};
+
+  MeshField::FieldElement f(meshInfo.numTet, field, elm);
+
+  Kokkos::View<MeshField::Real[1][4]> lc("localCoords");
+  Kokkos::deep_copy(lc, 0.5);
+  auto x = MeshField::evaluate(f, lc);
+}
+
 int main(int argc, char **argv) {
   Kokkos::initialize(argc, argv);
   triangleLocalPointEval();