Reduced Quintic Implicit

src/MeshField.hpp

@@ -238,6 +238,61 @@ 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(
+          "The mesh passed to %s must be 2D and simplex (triangles)\n",
+          __func__);
+    }
+  }
+
+  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);
+    const MeshField::LO triNode2DofHolder[18] = {
+        0,0,0,0,0,0,  // vertex 0 DOFs (6)
+        1,1,1,1,1,1,  // vertex 1 DOFs (6)
+        2,2,2,2,2,2   // vertex 2 DOFs (6)
+    };
+    const MeshField::Mesh_Topology triNode2DofHolderTopo[18] = {
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
+        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex
+    };
+    const auto dofHolderIdx = triNode2DofHolder[triNodeIdx];
+    const auto dofHolderTopo = triNode2DofHolderTopo[triNodeIdx];
+
+    Omega_h::LO osh_ent;
+    if (dofHolderTopo == MeshField::Vertex) {
+      const auto triDim = 2;
+      const auto vtxDim = 0;
+      const auto ignored = -1;
+      const auto localVtxIdx = (Omega_h::simplex_down_template(
+                                    triDim, vtxDim, dofHolderIdx, ignored) +
+                                2) %
+                               3;
+      const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
+      osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx];
+    } else {
+      assert(false);
+    }
+
+    return {0, triCompIdx, osh_ent, dofHolderTopo};
+  }
+};
+
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
   static_assert(ShapeOrder == 1 || ShapeOrder == 2);
   if constexpr (ShapeOrder == 1) {
@@ -254,6 +309,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) {
@@ -346,8 +408,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,46 @@ 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.873821971, 0.063089014, 0.063089014},
+                           0.0508449063702 / 2.0),
+          IntegrationPoint(Vector3{0.063089014, 0.873821971, 0.063089014},
+                           0.0508449063702 / 2.0),
+          IntegrationPoint(Vector3{0.063089014, 0.063089014, 0.873821971},
+                           0.0508449063702 / 2.0),
+
+          IntegrationPoint(Vector3{0.501426509, 0.249286745, 0.249286745},
+                           0.1167862757264 / 2.0),
+          IntegrationPoint(Vector3{0.249286745, 0.501426509, 0.249286745},
+                           0.1167862757264 / 2.0),
+          IntegrationPoint(Vector3{0.249286745, 0.249286745, 0.501426509},
+                           0.1167862757264 / 2.0),
+
+          IntegrationPoint(Vector3{0.636502499, 0.310352451, 0.053145050},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.636502499, 0.053145050, 0.310352451},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.310352451, 0.636502499, 0.053145050},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.310352451, 0.053145050, 0.636502499},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.053145050, 0.636502499, 0.310352451},
+                           0.0828510756184 / 2.0),
+          IntegrationPoint(Vector3{0.053145050, 0.310352451, 0.636502499},
+                           0.0828510756184 / 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

@@ -217,5 +217,89 @@ struct QuadraticTetrahedronShape {
   }
 };
 
+struct ReducedQuinticImplicitShape {
+  static const size_t numNodes = 18;
+  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];
+
+    const Real L1_2 = L1 * L1;
+    const Real L2_2 = L2 * L2;
+    const Real L3_2 = L3 * L3;
+    const Real L1_3 = L1_2 * L1;
+    const Real L2_3 = L2_2 * L2;
+    const Real L3_3 = L3_2 * L3;
+    const Real L1_4 = L1_3 * L1;
+    const Real L2_4 = L2_3 * L2;
+    const Real L3_4 = L3_3 * L3;
+    const Real L1_5 = L1_4 * L1;
+    const Real L2_5 = L2_4 * L2;
+    const Real L3_5 = L3_4 * L3;
+
+    Kokkos::Array<Real, numNodes> N;
+
+    N[0] = L1 * (2 * L1 - 1) * (3 * L1 - 1) * (4 * L1 - 1) * (5 * L1 - 1);
+    N[1] = L2 * (2 * L2 - 1) * (3 * L2 - 1) * (4 * L2 - 1) * (5 * L2 - 1);
+    N[2] = L3 * (2 * L3 - 1) * (3 * L3 - 1) * (4 * L3 - 1) * (5 * L3 - 1);
+
+    N[3] = 5.0 * L1_4 * L2;
+    N[4] = 5.0 * L2_4 * L1;
+    N[5] = 5.0 * L3_4 * L1;
+
+    N[6] = 5.0 * L1_4 * L3;
+    N[7] = 5.0 * L2_4 * L3;
+    N[8] = 5.0 * L3_4 * L2;
+
+    N[9]  = 10.0 * L1_3 * L2_2;
+    N[10] = 10.0 * L2_3 * L1_2;
+    N[11] = 10.0 * L3_3 * L1_2;
+
+    N[12] = 10.0 * L1_3 * L3_2;
+    N[13] = 10.0 * L2_3 * L3_2;
+    N[14] = 10.0 * L3_3 * L2_2;
+
+    N[15] = 20.0 * L1_2 * L2 * L3;
+    N[16] = 20.0 * L2_2 * L1 * L3;
+    N[17] = 20.0 * L3_2 * L1 * L2;
+
+    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];
+
+    const Real eps = 1e-6;
+
+    Vector3 xp = {xi[0] + eps, xi[1], 0.0};
+    Vector3 yp = {xi[0], xi[1] + eps, 0.0};
+
+    auto N0 = getValues(xi);
+    auto Nx = getValues(xp);
+    auto Ny = getValues(yp);
+
+    Kokkos::Array<Vector2, numNodes> dN;
+    for (int i = 0; i < numNodes; i++) {
+      dN[i][0] = (Nx[i] - N0[i]) / eps;
+      dN[i][1] = (Ny[i] - N0[i]) / eps;
+    }
+    return dN;
+  }
+};
+
 } // namespace MeshField
 #endif

test/testReducedQuinticImplicitIntegrator.cpp

@@ -0,0 +1,80 @@
+#include "KokkosController.hpp"
+#include "MeshField.hpp"
+#include "Omega_h_build.hpp"
+#include "Omega_h_file.hpp"
+#include "Omega_h_simplex.hpp"
+#include <Kokkos_Core.hpp>
+#include <MeshField_Integrate.hpp>
+#include <iostream>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+
+template <size_t dim>
+Omega_h::Mesh createMesh(Omega_h::Library &lib) {
+  auto world = lib.world();
+  const auto family = OMEGA_H_SIMPLEX;
+  auto len = 1.0;
+  return Omega_h::build_box(world, family, len, len, len,
+                            2, 2, (dim == 3 ? 2 : 0));
+}
+
+template <typename FieldElement>
+class ReducedQuinticImplicitIntegrator : public MeshField::Integrator {
+public:
+  ReducedQuinticImplicitIntegrator(FieldElement &fes_in)
+      : MeshField::Integrator(5), fes(fes_in), totalValue(0.0) {}
+
+  void atPoints(Kokkos::View<MeshField::Real **> p,
+                Kokkos::View<MeshField::Real *> w,
+                Kokkos::View<MeshField::Real *> dV) override {
+    const auto numPts = p.extent(0);
+    double localSum = 0.0;
+    Kokkos::parallel_reduce(
+        "IntegrateReducedQuinticImplicit", numPts,
+        KOKKOS_LAMBDA(const int i, double &sum) {
+          const double weight = w(i);
+          const double jac = dV(i);
+          const double f = 1.0; // Integrate constant field
+          sum += f * weight * jac;
+        },
+        localSum);
+    totalValue += localSum;
+  }
+
+  void post() override {
+    printf("[ReducedQuinticImplicit] Integrated Value = %e\n", totalValue);
+  }
+
+private:
+  FieldElement &fes;
+  double totalValue;
+};
+
+template <template <typename...> typename Controller, size_t dim>
+void runReducedQuinticImplicit(Omega_h::Mesh &mesh,
+                               MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
+  const auto ShapeOrder = 5;
+
+  auto field = omf.getCoordField();
+
+  const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
+  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
+
+  ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);
+  integrator.process(fes);
+}
+
+int main(int argc, char **argv) {
+  Kokkos::initialize(argc, argv);
+  Omega_h::Library lib(&argc, &argv);
+
+  {
+    auto mesh2D = createMesh<2>(lib);
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
+    runReducedQuinticImplicit<MeshField::KokkosController>(mesh2D, omf2D);
+  }
+
+  Kokkos::finalize();
+  return 0;
+}
+