Add integration in 2D on quads and triangles for PolynomialP

src/analysis/PolynomialP.hpp

@@ -2,6 +2,13 @@
 #define POLYNOMIALP_HPP
 
 #include <algebra/TinyVector.hpp>
+#include <analysis/CubeGaussQuadrature.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <analysis/SquareGaussQuadrature.hpp>
+#include <analysis/TriangleGaussQuadrature.hpp>
+#include <geometry/SquareTransformation.hpp>
+#include <geometry/TriangleTransformation.hpp>
 #include <utils/Exceptions.hpp>
 
 template <size_t N, size_t Dimension>
@@ -576,4 +583,45 @@ class PolynomialP
   ~PolynomialP()                               = default;
 };
 
+template <size_t N, size_t Dimension, size_t P>
+PUGS_INLINE double
+integrate(const PolynomialP<N, Dimension> Q, const std::array<TinyVector<Dimension>, P>& positions)
+{
+  double integral = 0.;
+  static_assert(P > 1, "For the integration, number of positions should be greater or equal to 2");
+  static_assert(N >= 0, "invalid degree");
+  if constexpr (Dimension == 1) {
+    static_assert(P == 2, "In 1D number of positions should be 2");
+    throw NotImplementedError("Not yet Available in 1D");
+  } else if constexpr (Dimension == 2) {
+    static_assert(P <= 4, "In 2D number of positions should be lesser or equal to 4");
+    if constexpr (P == 2) {
+    } else if constexpr (P == 3) {
+      const QuadratureFormula<2>& lN = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(N));
+      auto point_list                = lN.pointList();
+      auto weight_list               = lN.weightList();
+      TriangleTransformation<2> t{positions[0], positions[1], positions[2]};
+      auto value = weight_list[0] * Q(t(point_list[0]));
+      for (size_t i = 1; i < weight_list.size(); ++i) {
+        value += weight_list[i] * Q(t(point_list[i]));
+      }
+      integral = value;
+    } else {
+      const QuadratureFormula<2>& lN = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(N));
+      auto point_list                = lN.pointList();
+      auto weight_list               = lN.weightList();
+      SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
+      auto value = weight_list[0] * Q(s(point_list[0]));
+      for (size_t i = 1; i < weight_list.size(); ++i) {
+        value += weight_list[i] * Q(s(point_list[i]));
+      }
+      integral = value;
+    }
+  } else {
+    static_assert(Dimension == 3, "Dimension should be <=3");
+    throw NotImplementedError("Not yet Available in 3D");
+  }
+  return integral;
+}
+
 #endif   // POLYNOMIALP_HPP

tests/test_PolynomialP.cpp

-#include <catch2/catch_test_macros.hpp>
-
 #include <Kokkos_Core.hpp>
-
-#include <utils/PugsAssert.hpp>
-#include <utils/Types.hpp>
-
 #include <algebra/TinyMatrix.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
 #include <analysis/PolynomialP.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <analysis/SquareGaussQuadrature.hpp>
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
 
 // Instantiate to ensure full coverage is performed
 template class PolynomialP<0, 2>;
@@ -91,13 +92,25 @@ TEST_CASE("PolynomialP", "[analysis]")
   {
     TinyVector<6> coef(1, -2, 10, 7, 2, 9);
     TinyVector<10> coef2(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+    TinyVector<10> coef3(2, -4, -1, -3, 3, -5, -6, 0, 1, 7);
+    TinyVector<20> coef4(2, -4, -1, -3, 3, -5, -6, -2, 1, 7, 3, -2, 1, 2.5, 6, -9, 0.5, 4, -5, -8);
     PolynomialP<2, 2> P(coef);
     PolynomialP<3, 2> Q(coef2);
+    PolynomialP<2, 3> R(coef3);
+    PolynomialP<3, 3> S(coef4);
     TinyVector<2, size_t> relative_coef(1, 1);
     TinyVector<2, size_t> relative_coef2(1, 2);
+    TinyVector<3, size_t> relative_coef3(1, 0, 1);
+    TinyVector<3, size_t> relative_coef3b(0, 0, 2);
+    TinyVector<3, size_t> relative_coef4(1, 1, 1);
+    TinyVector<3, size_t> relative_coef5(0, 2, 1);
     REQUIRE(P[relative_coef] == 2);
     REQUIRE(Q[relative_coef2] == 1);
     REQUIRE(Q[relative_coef] == 3);
+    REQUIRE(R[relative_coef3] == -6);
+    REQUIRE(R[relative_coef3b] == 7);
+    REQUIRE(S[relative_coef4] == 6);
+    REQUIRE(S[relative_coef5] == 4);
   }
 
   SECTION("evaluation")
@@ -124,12 +137,55 @@ TEST_CASE("PolynomialP", "[analysis]")
     TinyVector<6> coef(1, -2, 10, 7, 2, 9);
     TinyVector<6> coef2(-2, 14, 2, 0, 0, 0);
     TinyVector<6> coef3(10, 2, 18, 0, 0, 0);
+    TinyVector<10> coef4(2, -4, -1, -3, 3, -5, -6, 1, 1, 7);
+    TinyVector<10> coef5(-1, -5, 2, 1, 0, 0, 0, 0, 0, 0);
     PolynomialP<2, 2> P(coef);
     PolynomialP<2, 2> Q(coef2);
     PolynomialP<2, 2> R(coef3);
+    PolynomialP<2, 3> S(coef4);
+    PolynomialP<2, 3> T(coef5);
 
     REQUIRE(Q == P.derivative(0));
     REQUIRE(R == P.derivative(1));
+    REQUIRE(T == S.derivative(1));
+  }
+
+  SECTION("integrate")
+  {
+    TinyVector<6> coef(1, -2, 10, 7, 2, 9);
+    std::array<TinyVector<2>, 4> positions;
+    std::array<TinyVector<2>, 3> positions2;
+    positions[0]  = TinyVector<2>{0, 0};
+    positions[1]  = TinyVector<2>{0, 0.5};
+    positions[2]  = TinyVector<2>{0.3, 0.7};
+    positions[3]  = TinyVector<2>{0.4, 0.1};
+    positions2[0] = TinyVector<2>{0, 0};
+    positions2[1] = TinyVector<2>{0, 0.5};
+    positions2[2] = TinyVector<2>{0.3, 0.7};
+    PolynomialP<2, 2> P(coef);
+    auto p1 = [](const TinyVector<2>& X) {
+      const double x = X[0];
+      const double y = X[1];
+      return 1 - 2. * x + 10 * y + 7 * x * x + 2 * x * y + 9 * y * y;
+    };
+    const QuadratureFormula<2>& l3 = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(3));
+    auto point_list                = l3.pointList();
+    auto weight_list               = l3.weightList();
+    SquareTransformation<2> s{positions[0], positions[1], positions[2], positions[3]};
+    auto value = weight_list[0] * p1(s(point_list[0]));
+    for (size_t i = 1; i < weight_list.size(); ++i) {
+      value += weight_list[i] * p1(s(point_list[i]));
+    }
+    const QuadratureFormula<2>& t3 = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(3));
+    auto point_list2               = t3.pointList();
+    auto weight_list2              = t3.weightList();
+    TriangleTransformation<2> t{positions2[0], positions2[1], positions2[2]};
+    auto value2 = weight_list2[0] * p1(t(point_list2[0]));
+    for (size_t i = 1; i < weight_list2.size(); ++i) {
+      value2 += weight_list2[i] * p1(t(point_list2[i]));
+    }
+    REQUIRE(value == integrate(P, positions));
+    REQUIRE(value2 == Catch::Approx(integrate(P, positions2)));
   }
 
   // // SECTION("product")