diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index e2830df47c1aef4919c8dc160c887b038b99c7fb..7572931bfee1c3476e610a6cb262a5e59e4d754c 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/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
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index b568179e3597bd4e80ff79442629f011500f9074..1533f9be9410de1ea4a31f784960d92b0c5093db 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -1,12 +1,13 @@
-#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")