Add edge and wedge quadrature with separate triangle and zeta orders

CMakeLists.txt

@@ -220,6 +220,10 @@ 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(TriangleIntegrationTest test/testTriangleIntegration.cpp)
+meshfields_add_exe(WedgeIntegrationTest test/testWedgeIntegration.cpp)
 meshfields_add_exe(Interpolate test/testInterpolate.cpp)
 meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp)
 
@@ -246,6 +250,10 @@ 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(TriangleIntegrationTest ./TriangleIntegrationTest)
+test_func(WedgeIntegrationTest ./WedgeIntegrationTest)
 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_Integrate.hpp

@@ -48,6 +48,102 @@ template <size_t pointSize> class EntityIntegration {
   virtual Integration<pointSize> const *getIntegration(int i) const = 0;
 };
 
+// Gauss-Legendre rules on [-1, 1]. N1 to N6 are copied from SCOREC/core
+// @ 7cd76473 apf/apfIntegrate.cc (EdgeIntegration, lines 27-129). The n-point
+// rule is exact for degree 2n - 1, and the weights sum to 2, the length of
+// the reference edge.
+//! [EdgeIntegration]
+class EdgeIntegration : public EntityIntegration<1> {
+public:
+  class N1 : public Integration<1> {
+  public:
+    virtual int countPoints() const { return 1; }
+    virtual std::vector<IntegrationPoint<1>> getPoints() const {
+      return {IntegrationPoint(Vector1{0}, 2, 1, 0)};
+    }
+    virtual int getAccuracy() const { return 1; }
+  }; // end N1
+  class N2 : public Integration<1> {
+  public:
+    virtual int countPoints() const { return 2; }
+    virtual std::vector<IntegrationPoint<1>> getPoints() const {
+      return {IntegrationPoint(Vector1{0.577350269189626}, 1, 1, 0),
+              IntegrationPoint(Vector1{-0.577350269189626}, 1, 1, 0)};
+    }
+    virtual int getAccuracy() const { return 3; }
+  }; // end N2
+  class N3 : public Integration<1> {
+  public:
+    virtual int countPoints() const { return 3; }
+    virtual std::vector<IntegrationPoint<1>> getPoints() const {
+      return {
+          IntegrationPoint(Vector1{0.000000000000000}, 0.888888888888889, 1, 0),
+          IntegrationPoint(Vector1{0.774596669241483}, 0.555555555555556, 1, 0),
+          IntegrationPoint(Vector1{-0.774596669241483}, 0.555555555555556, 1,
+                           0)};
+    }
+    virtual int getAccuracy() const { return 5; }
+  }; // end N3
+  class N4 : public Integration<1> {
+  public:
+    virtual int countPoints() const { return 4; }
+    virtual std::vector<IntegrationPoint<1>> getPoints() const {
+      return {
+          IntegrationPoint(Vector1{0.339981043584856}, 0.652145154862546, 1, 0),
+          IntegrationPoint(Vector1{-0.339981043584856}, 0.652145154862546, 1,
+                           0),
+          IntegrationPoint(Vector1{0.861136311594053}, 0.347854845137454, 1, 0),
+          IntegrationPoint(Vector1{-0.861136311594053}, 0.347854845137454, 1,
+                           0)};
+    }
+    virtual int getAccuracy() const { return 7; }
+  }; // end N4
+  class N5 : public Integration<1> {
+  public:
+    virtual int countPoints() const { return 5; }
+    virtual std::vector<IntegrationPoint<1>> getPoints() const {
+      return {
+          IntegrationPoint(Vector1{0.000000000000000}, 0.568888888888889, 1, 0),
+          IntegrationPoint(Vector1{0.538469310105683}, 0.478628670499366, 1, 0),
+          IntegrationPoint(Vector1{-0.538469310105683}, 0.478628670499366, 1,
+                           0),
+          IntegrationPoint(Vector1{0.906179845938664}, 0.236926885056189, 1, 0),
+          IntegrationPoint(Vector1{-0.906179845938664}, 0.236926885056189, 1,
+                           0)};
+    }
+    virtual int getAccuracy() const { return 9; }
+  }; // end N5
+  class N6 : public Integration<1> {
+  public:
+    virtual int countPoints() const { return 6; }
+    virtual std::vector<IntegrationPoint<1>> getPoints() const {
+      return {
+          IntegrationPoint(Vector1{0.238619186083197}, 0.467913934572691, 1, 0),
+          IntegrationPoint(Vector1{-0.238619186083197}, 0.467913934572691, 1,
+                           0),
+          IntegrationPoint(Vector1{0.661209386466265}, 0.360761573048139, 1, 0),
+          IntegrationPoint(Vector1{-0.661209386466265}, 0.360761573048139, 1,
+                           0),
+          IntegrationPoint(Vector1{0.932469514203152}, 0.171324492379170, 1, 0),
+          IntegrationPoint(Vector1{-0.932469514203152}, 0.171324492379170, 1,
+                           0)};
+    }
+    virtual int getAccuracy() const { return 11; }
+  }; // end N6
+  virtual int countIntegrations() const { return 6; }
+  virtual Integration<1> const *getIntegration(int i) const {
+    static N1 i1;
+    static N2 i2;
+    static N3 i3;
+    static N4 i4;
+    static N5 i5;
+    static N6 i6;
+    static Integration<1> *integrations[6] = {&i1, &i2, &i3, &i4, &i5, &i6};
+    return integrations[i];
+  }
+};
+//! [EdgeIntegration]
+
 // points are given in the reduced parametric coordinates (xi0=L1, xi1=L2,
 // ...); the redundant barycentric coordinate L0 = 1-sum(xi) is omitted
 //! [TriangleIntegration]
@@ -74,11 +170,40 @@ class TriangleIntegration : public EntityIntegration<2> {
     }
     virtual int getAccuracy() const { return 2; }
   }; // end N2
-  virtual int countIntegrations() const { return 2; }
+
+  // Symmetric Gaussian quadrature rules for the triangle from
+  //   D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature
+  //   rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985).
+  // N4 is the p = 4 rule of Appendix II (p. 1140), copied from
+  // SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). The weights
+  // sum to one, so they are divided by 2, the area of the reference triangle.
+  // The p = 3 rule has a negative weight, so order 3 uses N4.
+  class N4 : public Integration<2> {
+  public:
+    virtual int countPoints() const { return 6; }
+    virtual std::vector<IntegrationPoint<2>> getPoints() const {
+      return {IntegrationPoint(Vector2{0.816847572980459, 0.091576213509771},
+                               0.109951743655322 / 2.0, 2, 0),
+              IntegrationPoint(Vector2{0.091576213509771, 0.816847572980459},
+                               0.109951743655322 / 2.0, 2, 0),
+              IntegrationPoint(Vector2{0.091576213509771, 0.091576213509771},
+                               0.109951743655322 / 2.0, 2, 0),
+              IntegrationPoint(Vector2{0.108103018168070, 0.445948490915965},
+                               0.223381589678011 / 2.0, 2, 0),
+              IntegrationPoint(Vector2{0.445948490915965, 0.108103018168070},
+                               0.223381589678011 / 2.0, 2, 0),
+              IntegrationPoint(Vector2{0.445948490915965, 0.445948490915965},
+                               0.223381589678011 / 2.0, 2, 0)};
+    }
+    virtual int getAccuracy() const { return 4; }
+  }; // end N4
+
+  virtual int countIntegrations() const { return 3; }
   virtual Integration<2> const *getIntegration(int i) const {
     static N1 i1;
     static N2 i2;
-    static Integration<2> *integrations[2] = {&i1, &i2};
+    static N4 i4;
+    static Integration<2> *integrations[3] = {&i1, &i2, &i4};
     return integrations[i];
   }
 };
@@ -122,19 +247,53 @@ class TetrahedronIntegration : public EntityIntegration<3> {
     return integrations[i];
   }
 };
+
+// Wedge rules are the tensor product of a triangle rule in (xi0, xi1) and an
+// edge rule in zeta. The rule for (triOrder, zetaOrder) integrates
+// xi0^a xi1^b zeta^c exactly for a + b <= triOrder and c <= zetaOrder. The
+// weights sum to 1, the volume of the reference wedge (area 1/2 times length
+// 2).
+inline std::vector<IntegrationPoint<3>>
+getWedgeIntegrationPoints(int triOrder, int zetaOrder) {
+  const auto *tri = TriangleIntegration().getAccurate(triOrder);
+  if (tri == nullptr) {
+    fail("getWedgeIntegrationPoints: no triangle rule of order %d\n", triOrder);
+  }
+  const auto *edge = EdgeIntegration().getAccurate(zetaOrder);
+  if (edge == nullptr) {
+    fail("getWedgeIntegrationPoints: no edge rule of order %d\n", zetaOrder);
+  }
+  std::vector<IntegrationPoint<3>> points;
+  for (auto const &t : tri->getPoints())
+    for (auto const &z : edge->getPoints())
+      points.push_back(
+          IntegrationPoint(Vector3{t.param[0], t.param[1], z.param[0]},
+                           t.weight * z.weight, 3, 0));
+  return points;
+}
+
 //! [getIntegration]
 template <Mesh_Topology topo> auto const getIntegration() {
-  if constexpr (topo == Triangle) {
+  if constexpr (topo == Edge) {
+    return std::make_shared<EdgeIntegration>();
+  } else if constexpr (topo == Triangle) {
     return std::make_shared<TriangleIntegration>();
   } else if constexpr (topo == Tetrahedron) {
     return std::make_shared<TetrahedronIntegration>();
   }
   fail("getIntegration does not support given topology\n");
 }
 //! [getIntegration]
-template <Mesh_Topology topo> auto getIntegrationPoints(int order) {
-  auto ip = getIntegration<topo>()->getAccurate(order)->getPoints();
-  return ip;
+
+// zetaOrder is used only for wedges and defaults to order
+template <Mesh_Topology topo>
+auto getIntegrationPoints(int order, int zetaOrder = -1) {
+  if constexpr (topo == Wedge) {
+    return getWedgeIntegrationPoints(order, zetaOrder < 0 ? order : zetaOrder);
+  } else {
+    auto ip = getIntegration<topo>()->getAccurate(order)->getPoints();
+    return ip;
+  }
 }
 
 template <typename FieldElement>
@@ -208,7 +367,12 @@ auto getJacobianDeterminants(FieldElement &fes,
 class Integrator {
 public:
   /** \brief Construct an Integrator given an order of accuracy. */
-  Integrator(int o) : order(o) {}
+  Integrator(int o) : order(o), zetaOrder(o) {}
+  /** \brief Construct an Integrator for wedges, with separate orders of
+   * accuracy for the triangle (xi0, xi1) and zeta directions. */
+  Integrator(int triOrderIn, int zetaOrderIn)
+      : order(triOrderIn), zetaOrder(zetaOrderIn) {}
+
   virtual ~Integrator(){};
   /** \brief User callback: process entry
    *
@@ -248,7 +412,7 @@ class Integrator {
   template <typename FieldElement> void process(FieldElement &fes) {
     constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology();
     pre();
-    auto ip = getIntegrationPoints<topo[0]>(order);
+    auto ip = getIntegrationPoints<topo[0]>(order, zetaOrder);
     auto localCoords = getIntegrationPointLocalCoords(fes, ip);
     auto weights = getIntegrationPointWeights(fes, ip);
     auto dV = getJacobianDeterminants(fes, localCoords);
@@ -265,6 +429,7 @@ class Integrator {
 
 protected:
   int order;
+  int zetaOrder;
 };
 
 } // namespace MeshField

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