Add Wedge topology and LinearWedgeShape

CMakeLists.txt

@@ -220,6 +220,8 @@ meshfields_add_exe(IntegratorTest test/testIntegrator.cpp)
 meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp)
 
 meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp)
+meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp)
+meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp)
 meshfields_add_exe(Interpolate test/testInterpolate.cpp)
 meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp)
 
@@ -246,6 +248,8 @@ test_func(BoxIntegratorTest ./BoxIntegratorTest)
 smoke_test_func(OmegahTriTests ./OmegahTriTests)
 test_func(OmegahTetTest ./OmegahTetTest)
 test_func(IntegratorMeasure ./IntegratorMeasure)
+test_func(LinearWedgeShapeTest ./LinearWedgeShapeTest)
+test_func(LinearWedgeElementTest ./LinearWedgeElementTest)
 test_func(Interpolate ./Interpolate)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error

src/MeshField_Defines.hpp

@@ -13,7 +13,8 @@ enum Mesh_Topology {
   Tetrahedron,
   Hexahedron,
   Prism,
-  Pyramid
+  Pyramid,
+  Wedge
 };
 static bool Debug = false;
 const Real MachinePrecision = 1e-15;

src/MeshField_Shape.hpp

@@ -557,5 +557,113 @@ struct QuadraticTetrahedronShape {
   }
 };
 
+/**
+ * @brief Linear shape functions for 3D wedge (triangular prism) elements
+ *
+ * Defines linear basis functions for a wedge element with 6 nodes: the tensor product 
+ * of the linear (P1) triangle in (xi0, xi1) and the linear edge in
+ * zeta = xi2. 
+ * Parametric coordinate range: xi0, xi1 \f$\in\f$ [0, 1] with
+ * xi0 + xi1 \f$\le\f$ 1, and zeta \f$\in\f$ [-1, 1]
+ *
+ * Node ordering (node a+3 is joined to node a by an edge along \zeta):
+ * - Node 0 at (0, 0 , -1)
+ * - Node 1 at (1, 0 , -1)
+ * - Node 2 at (0, 1 , -1)
+ * - Node 3 at (0, 0 , 1)
+ * - Node 4 at (1, 0 , 1)
+ * - Node 5 at (0, 1 , 1)
+ *
+ * Shape functions use barycentric coordinates L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 of the triangle and, for a = 0, 1, 2:
+ * - \f$N_a     = L_a (1 - \zeta) / 2\f$
+ * - \f$N_{a+3} = L_a (1 + \zeta) / 2\f$
+ *
+ * Shape functions and ordering from:
+ * Zienkiewicz, Taylor, and Zhu
+ * 'The Finite Element Method: Its Basis and Fundamentals', 2013, Section 6.4,
+ * Fig. 6.16a: \f$N_a = \frac{1}{2} L_a (1 + \zeta_a \zeta)\f$, with 
+ * \f$\zeta_a = -1\f$ for nodes 0-2 and \f$+1\f$ for nodes 3-5 
+ */
+struct LinearWedgeShape {
+  // clang-format off
+  static const size_t numNodes = 6;                         ///< Number of nodes (6 for linear wedge)
+  static const size_t meshEntDim = 3;                       ///< Mesh entity dimension (3D)
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};  ///< DOFs located at vertices
+  constexpr static size_t Order = 1;                        ///< Polynomial order (linear)
+  // clang-format on
+
+  /**
+   * @brief Get parametric coordinates of element nodes
+   * @return Array of node coordinates [node0_xi0, node0_xi1, node0_zeta, ...]
+   */
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes * meshEntDim> getNodeParametricCoords() const {
+    // clang-format off
+    return {0,0,-1,   //node 0
+            1,0,-1,   //node 1
+            0,1,-1,   //node 2
+            0,0, 1,   //node 3
+            1,0, 1,   //node 4
+            0,1, 1};  //node 5
+    // clang-format on
+  }
+
+  /**
+   * @brief Evaluate shape functions at parametric coordinate
+   * @param xi Parametric coordinates (xi0, xi1, zeta) where xi0, xi1
+   *        \f$\in\f$ [0,1], xi0 + xi1 \f$\le\f$ 1 and zeta \f$\in\f$ [-1,1]
+   * @return Array of shape function values [N0, ..., N5] where 
+   * N_a = L_a ( 1 - \zeta)/2 and N_{a+3} = L_a (1+\zeta)/2 for a = 0,1,2
+   * with L0 = 1 - x0 - xi1, L1 = xi0 and L2 = xi1
+   */
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
+    assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol));
+    assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol));
+    assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol));
+    assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol));
+    const Real L0 = 1 - xi[0] - xi[1];
+    assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol));
+    const Real down = (1 - xi[2]) / 2; // linear edge function of the -1 face
+    const Real up = (1 + xi[2]) / 2;   // linear edge function of the +1 face
+    // clang-format off
+    return {L0*down,     //node 0
+            xi[0]*down,  //node 1
+            xi[1]*down,  //node 2
+            L0*up,       //node 3
+            xi[0]*up,    //node 4
+            xi[1]*up};   //node 5
+    // clang-format on
+  }
+
+  /**
+   * @brief Get gradients of shape functions in parametric coordinates
+   * @param xi Parametric coordinates (xi0, xi1, zeta), see getValues
+   * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN0/dzeta, dN1/dxi0, ...]
+   */
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, meshEntDim * numNodes>
+  getLocalGradients(Vector3 const &xi) const {
+    assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol));
+    assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol));
+    assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol));
+    assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol));
+    const Real L0 = 1 - xi[0] - xi[1];
+    assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol));
+    const Real L1 = xi[0];
+    const Real L2 = xi[1];
+    const Real down = (1 - xi[2]) / 2;
+    const Real up = (1 + xi[2]) / 2;
+    // clang-format off
+    return {-down, -down, -L0/2,  //node 0
+             down,     0, -L1/2,  //node 1
+                0,  down, -L2/2,  //node 2
+              -up,   -up,  L0/2,  //node 3
+               up,     0,  L1/2,  //node 4
+                0,    up,  L2/2}; //node 5
+    // clang-format on
+  }
+};
+
 } // namespace MeshField
 #endif

src/MeshField_ShapeField.hpp

@@ -25,6 +25,7 @@ struct MeshInfo {
   int numHex = 0;     // entDim = 3
   int numPrism = 0;   // entDim = 3
   int numPyramid = 0; // entDim = 3
+  int numWedge = 0;   // entDim = 3
   int dim = 0;
 
   /** @brief runtime map of number of entities based on topology */
@@ -46,6 +47,8 @@ struct MeshInfo {
       return numPrism;
     case Pyramid:
       return numPyramid;
+    case Wedge:
+      return numWedge;
     }
     return 0;
   }

test/testLinearWedgeElement.cpp

@@ -0,0 +1,193 @@
+#include "KokkosController.hpp"
+#include "MeshField_Element.hpp"
+#include "MeshField_Field.hpp"
+#include "MeshField_For.hpp"
+#include "MeshField_ShapeField.hpp"
+#include <Kokkos_Core.hpp>
+#include <Omega_h_vector.hpp>
+#include <array>
+#include <cassert>
+#include <iostream>
+
+// Checks MeshField::LinearWedgeShape inside a MeshField::FieldElement: field
+// evaluation (MeshField::evaluate) and the Jacobian (MeshField::getJacobians,
+// MeshField::getJacobianDeterminants) for two wedges, compared against the
+// geometric map x(xi) differentiated by hand. Failures are counted explicitly
+// so the test also fails in builds where assert() is disabled.
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MeshField::Real;
+using Vec3 = Omega_h::Vector<3>;
+
+constexpr int numWedges = 2;
+constexpr int numNodes = MeshField::LinearWedgeShape::numNodes;
+
+bool isClose(Real a, Real b) {
+  return Omega_h::are_close(a, b, 1e-12, 1e-12);
+}
+
+// Each wedge has its own six vertices: node n of wedge w is vertex 6*w + n.
+// Nodes 0-2 are the zeta = -1 triangle and nodes 3-5 the zeta = +1 triangle,
+// matching the LinearWedgeShape node ordering.
+struct WedgeToVertexField {
+  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
+  getTopology() {
+    return {MeshField::Wedge};
+  }
+
+  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO wedgeNodeIdx, MeshField::LO wedgeCompIdx,
+             MeshField::LO wedge, MeshField::Mesh_Topology topo) const {
+    assert(topo == MeshField::Wedge);
+    const MeshField::LO vtx = numNodes * wedge + wedgeNodeIdx;
+    return {0, wedgeCompIdx, vtx, MeshField::Vertex};
+  }
+};
+
+// clang-format off
+// wedge 0: right prism, triangle (0,0),(2,0),(0,3) swept by height 4
+//   x(xi) = (2 xi0, 3 xi1, 2 (1 + zeta)), J = diag(2, 3, 2), det J = 12
+// wedge 1: the top triangle is shifted, scaled and tilted, so J varies with
+//   xi0, xi1 and zeta
+const Vec3 wedgeCoords[numWedges][numNodes] = {
+  {{0, 0, 0},     {2, 0, 0},     {0, 3, 0},
+   {0, 0, 4},     {2, 0, 4},     {0, 3, 4}},
+  {{0, 0, 0},     {1, 0, 0},     {0, 1, 0},
+   {0.2, 0.1, 1}, {1.7, 0.3, 1.2}, {0.1, 1.4, 0.9}}};
+const char *wedgeNames[numWedges] = {"right prism", "skewed wedge"};
+
+// the six nodes, the centroid, and interior points
+const Vec3 localPts[] = {
+  {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1},
+  {1. / 3, 1. / 3, 0},
+  {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}};
+// clang-format on
+constexpr int numPts = static_cast<int>(sizeof(localPts) / sizeof(localPts[0]));
+constexpr int centroidPt = 6;
+
+// x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the
+// bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates
+Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) {
+  const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]};
+  Vec3 x = Omega_h::zero_vector<3>();
+  for (int a = 0; a < 3; ++a) {
+    x += L[a] * (1 - xi[2]) / 2 * nodes[a];
+    x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3];
+  }
+  return x;
+}
+
+// rows of the MeshField Jacobian, J(i,j) = dx_j/dxi_i, from differentiating
+// x(xi) directly:
+//   dx/dxi0  = (1-zeta)/2 (B1 - B0) + (1+zeta)/2 (T1 - T0)
+//   dx/dxi1  = (1-zeta)/2 (B2 - B0) + (1+zeta)/2 (T2 - T0)
+//   dx/dzeta = 1/2 sum_a L_a (T_a - B_a)
+std::array<Vec3, 3> expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) {
+  const Real down = (1 - xi[2]) / 2;
+  const Real up = (1 + xi[2]) / 2;
+  const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]};
+  const Vec3 dXi0 = down * (nodes[1] - nodes[0]) + up * (nodes[4] - nodes[3]);
+  const Vec3 dXi1 = down * (nodes[2] - nodes[0]) + up * (nodes[5] - nodes[3]);
+  Vec3 dZeta = Omega_h::zero_vector<3>();
+  for (int a = 0; a < 3; ++a)
+    dZeta += L[a] / 2 * (nodes[a + 3] - nodes[a]);
+  return {dXi0, dXi1, dZeta};
+}
+
+// Writes the wedge vertex coordinates into the coordinate field.
+template <typename CoordField>
+void setVtxCoords(CoordField field) {
+  const int numVtx = numWedges * numNodes;
+  Kokkos::View<Real *[3], Kokkos::HostSpace> coords_h("coords_h", numVtx);
+  for (int w = 0; w < numWedges; ++w)
+    for (int n = 0; n < numNodes; ++n)
+      for (int d = 0; d < 3; ++d)
+        coords_h(numNodes * w + n, d) = wedgeCoords[w][n][d];
+  auto coords = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), coords_h);
+  auto setCoordField = KOKKOS_LAMBDA(const int &vtx) {
+    for (int d = 0; d < 3; ++d)
+      field(vtx, 0, d, MeshField::Vertex) = coords(vtx, d);
+  };
+  MeshField::parallel_for<typename CoordField::Ctrlr>(
+      {0}, {numVtx}, setCoordField, "setCoordField");
+}
+
+// Compares the evaluated positions, Jacobians and determinants of wedge w with
+// the hand-differentiated map at every local point. Returns the failure count.
+template <class XView, class JView, class DetView>
+int checkWedge(int w, const XView &x_h, const JView &J_h,
+               const DetView &detJ_h) {
+  const auto &nodes = wedgeCoords[w];
+  int evalFailed = 0, jacFailed = 0, detFailed = 0;
+  for (int p = 0; p < numPts; ++p) {
+    const int row = w * numPts + p;
+    // at a node, x is that node; at the centroid, x is the mean of the six
+    // vertices; everywhere else x follows the map
+    Vec3 xExpected = expectedX(nodes, localPts[p]);
+    if (p < numNodes)
+      xExpected = nodes[p];
+    if (p == centroidPt) {
+      xExpected = Omega_h::zero_vector<3>();
+      for (int n = 0; n < numNodes; ++n)
+        xExpected += nodes[n] / Real(numNodes);
+    }
+    for (int d = 0; d < 3; ++d)
+      evalFailed += !isClose(x_h(row, d), xExpected[d]);
+
+    const auto jExpected = expectedJ(nodes, localPts[p]);
+    for (int i = 0; i < 3; ++i)
+      for (int j = 0; j < 3; ++j)
+        jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]);
+    const Real detExpected =
+        jExpected[0] * Omega_h::cross(jExpected[1], jExpected[2]);
+    detFailed += !isClose(detJ_h(row), detExpected);
+  }
+  std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w]
+            << ": evaluate failures " << evalFailed << ", Jacobian failures "
+            << jacFailed << ", det J failures " << detFailed << "\n";
+  return evalFailed + jacFailed + detFailed;
+}
+
+// Builds the coordinate field and the FieldElement, runs the evaluation,
+// Jacobian and determinant, and checks every wedge. Must be called between
+// Kokkos::initialize and Kokkos::finalize.
+int runLinearWedgeElementChecks() {
+  const MeshField::MeshInfo meshInfo{
+      .numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3};
+  auto coordField = MeshField::CreateCoordinateField<
+      ExecutionSpace, MeshField::KokkosController, 3>(meshInfo);
+  setVtxCoords(coordField.field);
+
+  MeshField::FieldElement fe(numWedges, coordField.field,
+                             MeshField::LinearWedgeShape(),
+                             WedgeToVertexField());
+
+  // the same local points in every wedge, stored element-major
+  Kokkos::View<Real **> lc("localCoords", numWedges * numPts, 3);
+  auto lc_h = Kokkos::create_mirror_view(lc);
+  for (int w = 0; w < numWedges; ++w)
+    for (int p = 0; p < numPts; ++p)
+      for (int d = 0; d < 3; ++d)
+        lc_h(w * numPts + p, d) = localPts[p][d];
+  Kokkos::deep_copy(lc, lc_h);
+
+  const auto x = MeshField::evaluate(fe, lc, numPts);
+  const auto J = MeshField::getJacobians(fe, lc, numPts);
+  const auto detJ = MeshField::getJacobianDeterminants(fe, J);
+  const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x);
+  const auto J_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J);
+  const auto detJ_h =
+      Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ);
+
+  int numFailed = 0;
+  for (int w = 0; w < numWedges; ++w)
+    numFailed += checkWedge(w, x_h, J_h, detJ_h);
+  return numFailed;
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  const int numFailed = runLinearWedgeElementChecks();
+  Kokkos::finalize();
+  return numFailed ? 1 : 0;
+}