diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index 7572931bfee1c3476e610a6cb262a5e59e4d754c..26f86c7b03010dd408cc3884b7e70eeacf21eea5 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -7,6 +7,7 @@
#include <analysis/QuadratureManager.hpp>
#include <analysis/SquareGaussQuadrature.hpp>
#include <analysis/TriangleGaussQuadrature.hpp>
+#include <geometry/LineTransformation.hpp>
#include <geometry/SquareTransformation.hpp>
#include <geometry/TriangleTransformation.hpp>
#include <utils/Exceptions.hpp>
@@ -596,6 +597,14 @@ integrate(const PolynomialP<N, Dimension> Q, const std::array<TinyVector<Dimensi
} 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) {
+ const QuadratureFormula<1>& lN = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(N));
+ const LineTransformation<2> t{positions[0], positions[1]};
+ double value = 0.;
+ for (size_t i = 0; i < lN.numberOfPoints(); ++i) {
+ value += lN.weight(i) * t.velocityNorm() * Q(t(lN.point(i)));
+ }
+ integral = value;
+
} else if constexpr (P == 3) {
const QuadratureFormula<2>& lN = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(N));
auto point_list = lN.pointList();
@@ -611,9 +620,9 @@ integrate(const PolynomialP<N, Dimension> Q, const std::array<TinyVector<Dimensi
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]));
+ auto value = weight_list[0] * s.jacobianDeterminant(point_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]));
+ value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * Q(s(point_list[i]));
}
integral = value;
}
diff --git a/src/analysis/TaylorPolynomial.hpp b/src/analysis/TaylorPolynomial.hpp
index 57021dad6e440940811147c1954263faa6cdbfc5..3ff9323279dd802f8faee2a62dcc36251419bfd2 100644
--- a/src/analysis/TaylorPolynomial.hpp
+++ b/src/analysis/TaylorPolynomial.hpp
@@ -7,6 +7,7 @@
#include <analysis/QuadratureManager.hpp>
#include <analysis/SquareGaussQuadrature.hpp>
#include <analysis/TriangleGaussQuadrature.hpp>
+#include <geometry/LineTransformation.hpp>
#include <geometry/SquareTransformation.hpp>
#include <geometry/TriangleTransformation.hpp>
#include <utils/Exceptions.hpp>
@@ -71,7 +72,7 @@ class TaylorPolynomial
PUGS_INLINE constexpr bool
operator==(const TaylorPolynomial& q) const
{
- return (m_coefficients == q.m_coefficients && m_x0 = q.m_x0);
+ return (m_coefficients == q.m_coefficients && m_x0 == q.m_x0);
}
PUGS_INLINE constexpr TaylorPolynomial(const TinyVector<size_coef, double>& coefficients,
@@ -96,7 +97,7 @@ class TaylorPolynomial
operator+(const TaylorPolynomial Q) const
{
Assert(m_x0 == Q.m_x0, "You cannot add Taylor polynomials with different origins");
- TaylorPolynomial<N, Dimension> P(m_coefficients);
+ TaylorPolynomial<N, Dimension> P(m_coefficients, m_x0);
for (size_t i = 0; i < size_coef; ++i) {
P.coefficients()[i] += Q.coefficients()[i];
}
@@ -108,6 +109,7 @@ class TaylorPolynomial
{
TaylorPolynomial<N, Dimension> P;
P.coefficients() = -coefficients();
+ P.m_x0 = m_x0;
return P;
}
@@ -115,7 +117,7 @@ class TaylorPolynomial
operator-(const TaylorPolynomial Q) const
{
Assert(m_x0 == Q.m_x0, "You cannot subtract Taylor polynomials with different origins");
- TaylorPolynomial<N, Dimension> P(m_coefficients);
+ TaylorPolynomial<N, Dimension> P(m_coefficients, m_x0);
P = P + (-Q);
return P;
}
@@ -154,7 +156,7 @@ class TaylorPolynomial
operator*(const double& lambda) const
{
TinyVector<size_coef> mult_coefs = lambda * m_coefficients;
- TaylorPolynomial<N, Dimension> M(mult_coefs);
+ TaylorPolynomial<N, Dimension> M(mult_coefs, m_x0);
return M;
}
@@ -307,7 +309,7 @@ class TaylorPolynomial
{
const auto P = *this;
TinyVector<size_coef> coefs(zero);
- TaylorPolynomial<N, Dimension> Q(coefs);
+ TaylorPolynomial<N, Dimension> Q(coefs, m_x0);
if constexpr (N != 0) {
Assert(var < Dimension,
"You can not derive a polynomial with respect to a variable of rank greater than the dimension");
@@ -623,6 +625,13 @@ integrate(const TaylorPolynomial<N, Dimension> Q, const std::array<TinyVector<Di
} 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) {
+ const QuadratureFormula<1>& lN = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(N));
+ const LineTransformation<2> t{positions[0], positions[1]};
+ double value = 0.;
+ for (size_t i = 0; i < lN.numberOfPoints(); ++i) {
+ value += lN.weight(i) * t.velocityNorm() * Q(t(lN.point(i)));
+ }
+ integral = value;
} else if constexpr (P == 3) {
const QuadratureFormula<2>& lN = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(N));
auto point_list = lN.pointList();
@@ -638,9 +647,9 @@ integrate(const TaylorPolynomial<N, Dimension> Q, const std::array<TinyVector<Di
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]));
+ auto value = weight_list[0] * s.jacobianDeterminant(point_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]));
+ value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * Q(s(point_list[i]));
}
integral = value;
}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a23fd187b9d29483f70f9be0dfc307b78f844445..af460c84b2b5ef021907c3ac5da7c09e23ac4469 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -118,6 +118,7 @@ add_executable (unit_tests
test_PETScUtils.cpp
test_Polynomial.cpp
test_PolynomialP.cpp
+ test_TaylorPolynomial.cpp
test_PolynomialBasis.cpp
test_PrimalToDiamondDualConnectivityDataMapper.cpp
test_PrimalToDual1DConnectivityDataMapper.cpp
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index 1533f9be9410de1ea4a31f784960d92b0c5093db..3ab52b3da50334c75404936dc456cdbd20954e9b 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -4,6 +4,7 @@
#include <analysis/PolynomialP.hpp>
#include <analysis/QuadratureManager.hpp>
#include <analysis/SquareGaussQuadrature.hpp>
+#include <analysis/TaylorPolynomial.hpp>
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <utils/PugsAssert.hpp>
@@ -155,6 +156,7 @@ TEST_CASE("PolynomialP", "[analysis]")
TinyVector<6> coef(1, -2, 10, 7, 2, 9);
std::array<TinyVector<2>, 4> positions;
std::array<TinyVector<2>, 3> positions2;
+ std::array<TinyVector<2>, 2> positions4;
positions[0] = TinyVector<2>{0, 0};
positions[1] = TinyVector<2>{0, 0.5};
positions[2] = TinyVector<2>{0.3, 0.7};
@@ -162,6 +164,9 @@ TEST_CASE("PolynomialP", "[analysis]")
positions2[0] = TinyVector<2>{0, 0};
positions2[1] = TinyVector<2>{0, 0.5};
positions2[2] = TinyVector<2>{0.3, 0.7};
+ positions4[0] = TinyVector<2>{0, 0.5};
+ positions4[1] = TinyVector<2>{0.3, 0.7};
+
PolynomialP<2, 2> P(coef);
auto p1 = [](const TinyVector<2>& X) {
const double x = X[0];
@@ -172,9 +177,9 @@ TEST_CASE("PolynomialP", "[analysis]")
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]));
+ auto value = weight_list[0] * s.jacobianDeterminant(point_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]));
+ value += weight_list[i] * s.jacobianDeterminant(point_list[i]) * p1(s(point_list[i]));
}
const QuadratureFormula<2>& t3 = QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(3));
auto point_list2 = t3.pointList();
@@ -184,8 +189,16 @@ TEST_CASE("PolynomialP", "[analysis]")
for (size_t i = 1; i < weight_list2.size(); ++i) {
value2 += weight_list2[i] * p1(t(point_list2[i]));
}
- REQUIRE(value == integrate(P, positions));
+ const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineFormula(GaussQuadratureDescriptor(2));
+ const LineTransformation<2> u{positions4[0], positions4[1]};
+ double value4 = 0.;
+ for (size_t i = 0; i < l1.numberOfPoints(); ++i) {
+ value4 += l1.weight(i) * u.velocityNorm() * p1(u(l1.point(i)));
+ }
+
+ REQUIRE(value == Catch::Approx(integrate(P, positions)));
REQUIRE(value2 == Catch::Approx(integrate(P, positions2)));
+ REQUIRE(value4 == Catch::Approx(integrate(P, positions4)));
}
// // SECTION("product")