diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index bcb89042b4b124a99e4e0ccbb4e0758690f26653..76d6bff82f6539b917fd5e2082b17a9b5a1469d3 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -161,7 +161,7 @@ class PolynomialP
size_t total_degree = relative_pos[0] + relative_pos[1];
absolute_position = total_degree * (total_degree + 1) / 2 + relative_pos[1];
} else {
- throw UnexpectedError("Not yet Available in 3D");
+ throw NotImplementedError("Not yet Available in 3D");
// static_assert(Dimension == 3);
// size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
// return (total_degree + 1) * (total_degree + 2) * (total_degree + 3) / 6 + relative_pos[1];
@@ -191,7 +191,7 @@ class PolynomialP
size_t total_degree = relative_pos[0] + relative_pos[1];
absolute_position = total_degree * (total_degree + 1) / 2 + relative_pos[1];
} else {
- throw UnexpectedError("Not yet Available in 3D");
+ throw NotImplementedError("Not yet Available in 3D");
// static_assert(Dimension == 3);
// size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
// return (total_degree + 1) * (total_degree + 2) * (total_degree + 3) / 6 + relative_pos[1];
@@ -203,7 +203,7 @@ class PolynomialP
PUGS_INLINE
constexpr double
- absolute_position(const TinyVector<Dimension, size_t> relative_pos)
+ absolute_position(const TinyVector<Dimension, size_t> relative_pos) const
{
size_t total_degree = 0;
for (size_t i = 0; i < Dimension; ++i) {
@@ -220,7 +220,7 @@ class PolynomialP
size_t total_degree = relative_pos[0] + relative_pos[1];
abs_pos = total_degree * (total_degree + 1) / 2 + relative_pos[1];
} else {
- throw UnexpectedError("Not yet Available in 3D");
+ throw NotImplementedError("Not yet Available in 3D");
// static_assert(Dimension == 3);
// size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
// return (total_degree + 1) * (total_degree + 2) * (total_degree + 3) / 6 + relative_pos[1];
@@ -257,7 +257,7 @@ class PolynomialP
value += valuex;
}
} else {
- throw UnexpectedError("Not yet Available in 3D");
+ throw NotImplementedError("Not yet Available in 3D");
}
return value;
@@ -276,45 +276,50 @@ class PolynomialP
}
}
- // PUGS_INLINE constexpr auto
- // derivative(const PolynomialP& P, const size_t var)
- // {
- // TinyVector<size_coef> coefs(zero);
- // PolynomialP<N, Dimension> Q(coefs);
- // if constexpr (N == 0) {
- // return Q;
- // } else {
- // Assert(var < Dimension,
- // "You can not derive a polynomial with respect to a variable of rank greater than the dimension");
- // if constexpr (Dimension == 1) {
- // for (size_t i = 0; i < size_coef; ++i) {
- // coefs[i] = double(i + 1) * P.coefficients()[i + 1];
- // }
- // return Q;
- // } else if constexpr (Dimension == 2) {
- // if (var == 0) {
- // for (size_t i = 0; i < N; ++i) {
- // for (size_t j = 0; j < N; ++i) {
- // TinyVector<Dimension, size_t> relative_pos(i, j);
- // TinyVector<Dimension, size_t> relative_posp(i + 1, j);
- // Q[relative_pos] = double(i + 1) * P[relative_posp];
- // }
- // }
- // } else {
- // for (size_t i = 0; i < N; ++i) {
- // for (size_t j = 0; j < N; ++i) {
- // TinyVector<Dimension, size_t> relative_pos(i, j);
- // TinyVector<Dimension, size_t> relative_posp(i, j + 1);
- // Q[relative_pos] = double(j + 1) * P[relative_posp];
- // }
- // }
- // }
- // return Q;
- // } else {
- // throw UnexpectedError("Not yet Available in 3D");
- // }
- // }
- // }
+ PUGS_INLINE constexpr auto
+ derivative(const size_t var) const
+ {
+ const auto P = *this;
+ TinyVector<size_coef> coefs(zero);
+ PolynomialP<N, Dimension> Q(coefs);
+ if constexpr (N == 0) {
+ return Q;
+ } else {
+ Assert(var < Dimension,
+ "You can not derive a polynomial with respect to a variable of rank greater than the dimension");
+ if constexpr (Dimension == 1) {
+ for (size_t i = 0; i < size_coef; ++i) {
+ coefs[i] = double(i + 1) * P.coefficients()[i + 1];
+ }
+ return Q;
+ } else if constexpr (Dimension == 2) {
+ if (var == 0) {
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N - i; ++j) {
+ TinyVector<Dimension, size_t> relative_pos(i, j);
+ TinyVector<Dimension, size_t> relative_posp(i + 1, j);
+ size_t absolute_position = Q.absolute_position(relative_pos);
+ size_t absolute_positionp = P.absolute_position(relative_posp);
+ Q.coefficients()[absolute_position] = double(i + 1) * m_coefficients[absolute_positionp];
+ }
+ }
+ } else {
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N - i; ++j) {
+ TinyVector<Dimension, size_t> relative_pos(i, j);
+ TinyVector<Dimension, size_t> relative_posp(i, j + 1);
+ size_t absolute_position = Q.absolute_position(relative_pos);
+ size_t absolute_positionp = P.absolute_position(relative_posp);
+ Q.coefficients()[absolute_position] = double(j + 1) * m_coefficients[absolute_positionp];
+ }
+ }
+ }
+ return Q;
+ } else {
+ throw UnexpectedError("Not yet Available in 3D");
+ }
+ }
+ }
PUGS_INLINE constexpr PolynomialP() noexcept = default;
~PolynomialP() = default;
diff --git a/tests/test_PolynomialP.cpp b/tests/test_PolynomialP.cpp
index 1855c7b671d4201f55e471f2655bd86e6a4f9c5e..b568179e3597bd4e80ff79442629f011500f9074 100644
--- a/tests/test_PolynomialP.cpp
+++ b/tests/test_PolynomialP.cpp
@@ -119,6 +119,19 @@ TEST_CASE("PolynomialP", "[analysis]")
REQUIRE(Q(pos) == -24);
}
+ SECTION("derivation")
+ {
+ 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);
+ PolynomialP<2, 2> P(coef);
+ PolynomialP<2, 2> Q(coef2);
+ PolynomialP<2, 2> R(coef3);
+
+ REQUIRE(Q == P.derivative(0));
+ REQUIRE(R == P.derivative(1));
+ }
+
// // SECTION("product")
// {
// Polynomial<2> P(2, 3, 4);