Reduced Quintic Implicit

CMakeLists.txt

@@ -172,6 +172,8 @@ meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
+meshfields_add_exe(ReducedQuinticImplicitIntegratorTest test/testReducedQuinticImplicitIntegrator.cpp)
+meshfields_add_exe(ReducedQuinticImplicitEvaluateTest test/testReducedQuinticImplicitEvaluate.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -190,6 +192,8 @@ test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
 test_func(OmegahTetTest, ./OmegahTetTest)
+test_func(ReducedQuinticImplicitIntegratorTest ./ReducedQuinticImplicitIntegratorTest)
+test_func(ReducedQuinticImplicitEvalauteTest ./ReducedQuinticImplicitEvaluateTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)

src/MeshField.hpp

@@ -9,6 +9,7 @@
 #include "Omega_h_file.hpp"    //move
 #include "Omega_h_mesh.hpp"    //move
 #include "Omega_h_simplex.hpp" //move
+#include "MeshField_Shape.hpp"
 
 namespace {
 
@@ -238,8 +239,77 @@ struct QuadraticTetrahedronToField {
   }
 };
 
+struct ReducedQuinticTriangleToField {
+  Omega_h::LOs triVerts;
+
+  ReducedQuinticTriangleToField(Omega_h::Mesh &mesh)
+      : triVerts(mesh.ask_elem_verts()) {
+    if (mesh.dim() != 2 || mesh.family() != OMEGA_H_SIMPLEX) {
+      MeshField::fail(
+          "ReducedQuinticTriangleToField requires 2D simplex mesh\n");
+    }
+  }
+
+  static constexpr KOKKOS_FUNCTION
+  Kokkos::Array<MeshField::Mesh_Topology, 1> getTopology() {
+    return {MeshField::Triangle};
+  }
+
+  KOKKOS_FUNCTION
+  MeshField::ElementToDofHolderMap
+  operator()(MeshField::LO triNodeIdx,
+             MeshField::LO triCompIdx,
+             MeshField::LO tri,
+             MeshField::Mesh_Topology topo) const
+  {
+    assert(topo == MeshField::Triangle);
+
+    MeshField::LO osh_ent = -1;
+    MeshField::LO nodeInHolder = -1;
+
+    if (triNodeIdx == 0) {
+      osh_ent = triVerts[3*tri + 0];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx == 1) {
+      osh_ent = triVerts[3*tri + 1];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx == 2) {
+      osh_ent = triVerts[3*tri + 2];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx >= 3 && triNodeIdx <= 6) {
+      const MeshField::LO v0 = triVerts[3*tri + 0];
+      const MeshField::LO v1 = triVerts[3*tri + 1];
+      osh_ent = v0;  
+      nodeInHolder = triNodeIdx - 3 + 1;
+    }
+    else if (triNodeIdx >= 7 && triNodeIdx <= 10) {
+      const MeshField::LO v1 = triVerts[3*tri + 1];
+      osh_ent = v1;
+      nodeInHolder = triNodeIdx - 7 + 1;
+    }
+    else if (triNodeIdx >= 11 && triNodeIdx <= 14) {
+      const MeshField::LO v2 = triVerts[3*tri + 2];
+      osh_ent = v2;
+      nodeInHolder = triNodeIdx - 11 + 1;
+    }
+    else if (triNodeIdx >= 15 && triNodeIdx <= 20) {
+      const MeshField::LO v0 = triVerts[3*tri + 0];
+      osh_ent = v0;
+      nodeInHolder = triNodeIdx - 15 + 5;
+    }
+    else {
+      MeshField::fail("Invalid node index in ReducedQuinticTriangleToField\n");
+    }
+    return { nodeInHolder, triCompIdx, osh_ent, MeshField::Vertex };
+  }
+};
+
+
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
-  static_assert(ShapeOrder == 1 || ShapeOrder == 2);
+  static_assert(ShapeOrder == 1 || ShapeOrder == 2 || ShapeOrder == 5);
   if constexpr (ShapeOrder == 1) {
     struct result {
       MeshField::LinearTriangleShape shp;
@@ -254,6 +324,13 @@ template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
     };
     return result{MeshField::QuadraticTriangleShape(),
                   QuadraticTriangleToField(mesh)};
+  } else if constexpr (ShapeOrder == 5) {
+    struct result {
+      MeshField::ReducedQuinticImplicitShape shp;
+      ReducedQuinticTriangleToField map;
+    };
+    return result{MeshField::ReducedQuinticImplicitShape(),
+                  ReducedQuinticTriangleToField(mesh)};
   }
 }
 template <int ShapeOrder> auto getTetrahedronElement(Omega_h::Mesh &mesh) {
@@ -275,6 +352,15 @@ template <int ShapeOrder> auto getTetrahedronElement(Omega_h::Mesh &mesh) {
   }
 }
 
+inline auto getReducedQuinticImplicitElement(Omega_h::Mesh &mesh) {
+  struct result {
+    MeshField::ReducedQuinticImplicitShape shp;
+    ReducedQuinticTriangleToField map;
+  };
+  return result{MeshField::ReducedQuinticImplicitShape(),
+                ReducedQuinticTriangleToField(mesh)};
+}
+
 } // namespace Omegah
 
 template <typename ExecutionSpace, size_t dim,
@@ -346,8 +432,8 @@ class OmegahMeshField {
       MeshField::fail("input mesh must be 2d\n");
     }
     const auto ShapeOrder = ShapeField::Order;
-    if (ShapeOrder != 1 && ShapeOrder != 2) {
-      MeshField::fail("input field order must be 1 or 2\n");
+    if (ShapeOrder != 1 && ShapeOrder != 2 && ShapeOrder != 5) {
+      MeshField::fail("input field order must be 1 or 2 or 5\n");
     }
 
     const auto [shp, map] = Omegah::getTriangleElement<ShapeOrder>(mesh);

src/MeshField_Integrate.hpp

@@ -67,11 +67,56 @@ class TriangleIntegration : public EntityIntegration<3> {
     }
     virtual int getAccuracy() const { return 2; }
   }; // end N2
-  virtual int countIntegrations() const { return 2; }
+    class N6 : public Integration<3> {
+  public:
+    virtual int countPoints() const { return 12; }
+    virtual std::vector<IntegrationPoint<3>> getPoints() const {
+    return {
+        IntegrationPoint(
+            Vector3{0.8738219710169970, 0.0630890144915022, 0.0630890144915008},
+            0.050844906370206816 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0630890144915022, 0.8738219710169970, 0.0630890144915008},
+            0.050844906370206816 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0630890144915022, 0.0630890144915022, 0.8738219710169956},
+            0.050844906370206816 / 2.0),
+        IntegrationPoint(
+            Vector3{0.5014265096581792, 0.2492867451709104, 0.2492867451709104},
+            0.11678627572637937 / 2.0),
+        IntegrationPoint(
+            Vector3{0.2492867451709104, 0.5014265096581792, 0.2492867451709104},
+            0.11678627572637937 / 2.0),
+        IntegrationPoint(
+            Vector3{0.2492867451709104, 0.2492867451709104, 0.5014265096581792},
+            0.11678627572637937 / 2.0),
+        IntegrationPoint(
+            Vector3{0.6365024991213986, 0.3103524510337844, 0.053145049844817},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.6365024991213986, 0.0531450498448170, 0.3103524510337844},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.3103524510337844, 0.6365024991213986, 0.0531450498448170},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.3103524510337844, 0.0531450498448170, 0.6365024991213986},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0531450498448170, 0.6365024991213986, 0.3103524510337844},
+            0.082851075618373575 / 2.0),
+        IntegrationPoint(
+            Vector3{0.0531450498448170, 0.3103524510337844, 0.6365024991213986},
+            0.082851075618373575 / 2.0)};
+    }
+    virtual int getAccuracy() const { return 6; }
+  }; // end N6
+  virtual int countIntegrations() const { return 3; }
   virtual Integration<3> const *getIntegration(int i) const {
     static N1 i1;
     static N2 i2;
-    static Integration<3> *integrations[2] = {&i1, &i2};
+    static N6 i6;
+    static Integration<3> *integrations[3] = {&i1, &i2, &i6};
     return integrations[i];
   }
 };

src/MeshField_Shape.hpp

@@ -5,24 +5,15 @@
 // getValues(...) implementation copied from
 // SCOREC/core apf/apfShape.cc @ 7cd76473
 
-namespace {
-template <typename Array> KOKKOS_INLINE_FUNCTION bool sumsToOne(Array &xi) {
-  auto sum = 0.0;
-  for (int i = 0; i < xi.size(); i++) {
-    sum += xi[i];
-  }
-  return (Kokkos::fabs(sum - 1) <= MeshField::MachinePrecision);
-}
-
-template <typename Array>
-KOKKOS_INLINE_FUNCTION bool greaterThanOrEqualZero(Array &xi) {
-  auto gt = true;
-  for (int i = 0; i < xi.size(); i++) {
-    gt = gt && (xi[i] >= 0);
-  }
-  return gt;
-}
-} // namespace
+namespace { template <typename Array> KOKKOS_INLINE_FUNCTION bool
+	sumsToOne(Array &xi) { auto sum = 0.0; for (int i = 0; i < xi.size();
+			i++) { sum += xi[i]; } return (Kokkos::fabs(sum - 1) <=
+				MeshField::MachinePrecision); }
+
+template <typename Array> KOKKOS_INLINE_FUNCTION bool
+	greaterThanOrEqualZero(Array &xi) { auto gt = true; for (int i = 0; i <
+			xi.size(); i++) { gt = gt && (xi[i] >= 0); } return gt;
+	} } // namespace
 
 namespace MeshField {
 
@@ -217,5 +208,110 @@ struct QuadraticTetrahedronShape {
   }
 };
 
+
+struct ReducedQuinticImplicitShape {
+  static const size_t numNodes = 21;
+  static const size_t meshEntDim = 2;
+  constexpr static Mesh_Topology DofHolders[1] = {Vertex};
+  constexpr static size_t Order = 5;
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes> getValues(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+
+    const Real L1 = 1.0 - xi[0] - xi[1];
+    const Real L2 = xi[0];
+    const Real L3 = xi[1];
+
+    Real powL1[6], powL2[6], powL3[6];
+    powL1[0] = powL2[0] = powL3[0] = 1.0;
+    for (int p = 1; p <= 5; ++p) {
+      powL1[p] = powL1[p - 1] * L1;
+      powL2[p] = powL2[p - 1] * L2;
+      powL3[p] = powL3[p - 1] * L3;
+    }
+
+    auto fact = [](int n) {
+      double r = 1.0;
+      for (int i = 2; i <= n; ++i) r *= double(i);
+      return r;
+    };
+    const double f5 = fact(5);
+
+    Kokkos::Array<Real, numNodes> N;
+    int idx = 0;
+    for (int i = 0; i <= 5; ++i) {
+      for (int j = 0; j <= 5 - i; ++j) {
+        int k = 5 - i - j;
+        double coeff = f5 / (fact(i) * fact(j) * fact(k));
+        N[idx++] = coeff * powL1[i] * powL2[j] * powL3[k];
+      }
+    }
+    return N;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Vector2, numNodes> getLocalGradients(Vector3 const &xi) const {
+    assert(greaterThanOrEqualZero(xi));
+    assert(sumsToOne(xi));
+
+    const Real L1 = 1.0 - xi[0] - xi[1];
+    const Real L2 = xi[0];
+    const Real L3 = xi[1];
+
+    Real powL1[6], powL2[6], powL3[6];
+    powL1[0] = powL2[0] = powL3[0] = 1.0;
+    for (int p = 1; p <= 5; ++p) {
+      powL1[p] = powL1[p - 1] * L1;
+      powL2[p] = powL2[p - 1] * L2;
+      powL3[p] = powL3[p - 1] * L3;
+    }
+
+    auto fact = [](int n) {
+      double r = 1.0;
+      for (int i = 2; i <= n; ++i) r *= double(i);
+      return r;
+    };
+    const double f5 = fact(5);
+
+    Kokkos::Array<Vector2, numNodes> dN;
+    int idx = 0;
+    for (int i = 0; i <= 5; ++i) {
+      for (int j = 0; j <= 5 - i; ++j) {
+        int k = 5 - i - j;
+        double coeff = f5 / (fact(i) * fact(j) * fact(k));
+
+        double dN_dL1 = 0.0;
+        double dN_dL2 = 0.0;
+        double dN_dL3 = 0.0;
+        if (i > 0) dN_dL1 = coeff * double(i) * powL1[i - 1] * powL2[j] * powL3[k];
+        if (j > 0) dN_dL2 = coeff * double(j) * powL1[i] * powL2[j - 1] * powL3[k];
+        if (k > 0) dN_dL3 = coeff * double(k) * powL1[i] * powL2[j] * powL3[k - 1];
+
+        const double dNdX = -dN_dL1 + dN_dL2;
+        const double dNdY = -dN_dL1 + dN_dL3;
+
+        dN[idx][0] = dNdX;
+        dN[idx][1] = dNdY;
+        ++idx;
+      }
+    }
+    return dN;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  Kokkos::Array<Real, numNodes * meshEntDim> getLocalGradients() const {
+    Vector3 xi = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
+    auto gvec = getLocalGradients(xi);
+    Kokkos::Array<Real, numNodes * meshEntDim> flat{};
+    for (int n = 0; n < (int)numNodes; ++n) {
+      flat[n * meshEntDim + 0] = gvec[n][0];
+      flat[n * meshEntDim + 1] = gvec[n][1];
+    }
+    return flat;
+  }
+};
+
 } // namespace MeshField
 #endif

src/MeshField_ShapeField.hpp

@@ -164,8 +164,8 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
                 "CreateLagrangeField only supports single and double precision "
                 "floating point fields\n");
   static_assert(
-      (order == 1 || order == 2),
-      "CreateLagrangeField only supports linear and quadratic fields\n");
+      (order == 1 || order == 2 || order == 5),
+      "CreateLagrangeField only supports linear, quadratic, and quintic fields\n");
   static_assert((dim == 1 || dim == 2 || dim == 3),
                 "CreateLagrangeField only supports 1d, 2d, and 3d meshes\n");
   using MemorySpace = typename ExecutionSpace::memory_space;

test/testCountIntegrator.cpp

@@ -74,7 +74,7 @@ void doRun(Omega_h::Mesh &mesh,
 
   CountIntegrator countInt(fes);
   countInt.process(fes);
-  assert(mesh.nelems() == countInt.getCount());
+  assert(static_cast<size_t>(mesh.nelems()) == countInt.getCount());
 }
 
 int main(int argc, char **argv) {