diff --git a/src/analysis/PolynomialP.hpp b/src/analysis/PolynomialP.hpp
index 253e6369b9f0d4d087e7cc8eab14c70b3b04230e..e2830df47c1aef4919c8dc160c887b038b99c7fb 100644
--- a/src/analysis/PolynomialP.hpp
+++ b/src/analysis/PolynomialP.hpp
@@ -161,11 +161,12 @@ 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 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];
- // return (N + 1) * (N + 2) * (N + 3) / 6;
+ // throw NotImplementedError("Not yet Available in 3D");
+ static_assert(Dimension == 3);
+ size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
+ size_t total_sub_degree = relative_pos[1] + relative_pos[2];
+ return total_degree * (total_degree + 1) * (total_degree + 2) / 6 +
+ total_sub_degree * (total_sub_degree + 1) / 2 + relative_pos[2];
}
return absolute_coefs[absolute_position];
@@ -191,11 +192,12 @@ 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 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];
- // return (N + 1) * (N + 2) * (N + 3) / 6;
+ // throw NotImplementedError("Not yet Available in 3D");
+ static_assert(Dimension == 3);
+ size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
+ size_t total_sub_degree = relative_pos[1] + relative_pos[2];
+ absolute_position = total_degree * (total_degree + 1) * (total_degree + 2) / 6 +
+ total_sub_degree * (total_sub_degree + 1) / 2 + relative_pos[2];
}
return absolute_coefs[absolute_position];
@@ -219,11 +221,12 @@ class PolynomialP
} else if constexpr (Dimension == 2) {
abs_pos = total_degree * (total_degree + 1) / 2 + relative_pos[1];
} else {
- 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];
- // return (N + 1) * (N + 2) * (N + 3) / 6;
+ static_assert(Dimension == 3);
+ size_t total_degree = relative_pos[0] + relative_pos[1] + relative_pos[2];
+ size_t total_sub_degree = relative_pos[1] + relative_pos[2];
+ abs_pos = total_degree * (total_degree + 1) * (total_degree + 2) / 6 +
+ total_sub_degree * (total_sub_degree + 1) / 2 + relative_pos[2];
+ // throw NotImplementedError("Not yet Available in 3D");
}
return abs_pos;
@@ -275,22 +278,19 @@ class PolynomialP
}
}
- PUGS_INLINE constexpr auto
+ PUGS_INLINE constexpr PolynomialP<N, Dimension>
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 {
+ if constexpr (N != 0) {
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) {
+ for (size_t i = 0; i < size_coef - 1; ++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) {
@@ -313,11 +313,49 @@ class PolynomialP
}
}
}
- return Q;
} else {
- throw NotImplementedError("Not yet Available in 3D");
+ static_assert(Dimension == 3);
+ if (var == 0) {
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N - i; ++j) {
+ for (size_t k = 0; k < N - i - j; ++k) {
+ TinyVector<Dimension, size_t> relative_pos(i, j, k);
+ TinyVector<Dimension, size_t> relative_posp(i + 1, j, k);
+ 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 if (var == 1) {
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N - i; ++j) {
+ for (size_t k = 0; k < N - i - j; ++k) {
+ TinyVector<Dimension, size_t> relative_pos(i, j, k);
+ TinyVector<Dimension, size_t> relative_posp(i, j + 1, k);
+ 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];
+ }
+ }
+ }
+ } else {
+ for (size_t i = 0; i < N; ++i) {
+ for (size_t j = 0; j < N - i; ++j) {
+ for (size_t k = 0; k < N - i - j; ++k) {
+ TinyVector<Dimension, size_t> relative_pos(i, j, k);
+ TinyVector<Dimension, size_t> relative_posp(i, j, k + 1);
+ size_t absolute_position = Q.absolute_position(relative_pos);
+ size_t absolute_positionp = P.absolute_position(relative_posp);
+ Q.coefficients()[absolute_position] = double(k + 1) * m_coefficients[absolute_positionp];
+ }
+ }
+ }
+ }
+ // throw NotImplementedError("Not yet Available in 3D");
}
}
+ return Q;
}
PUGS_INLINE
@@ -445,7 +483,92 @@ class PolynomialP
return os;
} else {
- throw NotImplementedError("Not yet Available in 3D");
+ // size_t i = 0;
+ // size_t j = 0;
+ // size_t k = N;
+ // TinyVector<Dimension, size_t> rel_pos(i, j, k);
+ // double coef = P[rel_pos];
+ // if (coef != 0.) {
+ // if (coef < 0.) {
+ // os << " - ";
+ // }
+ // if (coef != 1 && coef != -1) {
+ // os << std::abs(coef);
+ // }
+ // os << "z^" << N;
+ // }
+ // size_t degree = N;
+ // for (size_t l = size_coef - 1; l > 0; --l) {
+ // if (k > 0) {
+ // k--;
+ // if (j < k) {
+ // j++;
+ // } else {
+ // j--;
+ // i++;
+ // }
+ // } else {
+ // degree--;
+ // k = degree;
+ // i = 0;
+ // j = 0;
+ // }
+
+ // rel_pos = TinyVector<Dimension, size_t>(i, j, k);
+ // double coef = P[rel_pos];
+ // if (coef != 0.) {
+ // if (coef > 0.) {
+ // os << " + ";
+ // } else if (coef < 0.) {
+ // os << " - ";
+ // }
+ // if ((coef != 1 && coef != -1) || (i == 0 && j == 0 && k == 0)) {
+ // os << std::abs(coef);
+ // }
+ // if (i == 0 && j == 0 && k == 0)
+ // continue;
+ // if (i == 0 && j == 0) {
+ // if (k != 1) {
+ // os << "z^" << j;
+ // } else {
+ // os << "z";
+ // }
+ // } else if (i == 0 && k == 0) {
+ // if (j == 1) {
+ // os << "y";
+ // } else {
+ // os << "y^" << i;
+ // }
+ // } else if (j == 0 && k == 0) {
+ // if (i == 1) {
+ // os << "x";
+ // } else {
+ // os << "x^" << i;
+ // }
+ // } else {
+ // if (i == 1 && j == 1 && k == 1) {
+ // os << "xyz";
+ // } else if (i == 1) {
+ // os << "x"
+ // << "y^" << j << "z^" << k;
+ // } else if (j == 1) {
+ // os << "x^" << i << "y"
+ // << "z^" << k;
+ // } else if (k == 1) {
+ // os << "x^" << i << "y^" << j << "z";
+ // } else {
+ // os << "x^" << i << "y^" << j << "z^" << k;
+ // }
+ // }
+ // all_coef_zero = false;
+ // }
+ //
+ for (size_t l = 0; l < size_coef; ++l) {
+ double coef = P.coefficients()[l];
+ os << coef << " ";
+ }
+ return os;
+ // throw NotImplementedError("Not yet Available in 3D");
}
}
@@ -453,173 +576,4 @@ class PolynomialP
~PolynomialP() = default;
};
-// template <size_t M>
-// PUGS_INLINE constexpr friend void
-// divide(const PolynomialP<N>& P1, const PolynomialP<M>& P2, PolynomialP<N>& Q, PolynomialP<N>& R)
-// {
-// const size_t Nr = P1.realDegree();
-// const size_t Mr = P2.realDegree();
-// R.coefficients() = P1.coefficients();
-// Q.coefficients() = zero;
-// for (ssize_t k = Nr - Mr; k >= 0; --k) {
-// Q.coefficients()[k] = R.coefficients()[Mr + k] / P2.coefficients()[Mr];
-// for (ssize_t j = Mr + k; j >= k; --j) {
-// R.coefficients()[j] -= Q.coefficients()[k] * P2.coefficients()[j - k];
-// }
-// }
-// for (size_t j = Mr; j <= Nr; ++j) {
-// R.coefficients()[j] = 0;
-// }
-// }
-
-// PUGS_INLINE
-// constexpr friend PolynomialP<N + 1>
-// primitive(const PolynomialP<N>& P)
-// {
-// TinyVector<N + 2> coefs;
-// for (size_t i = 0; i < N + 1; ++i) {
-// coefs[i + 1] = P.coefficients()[i] / double(i + 1);
-// }
-// coefs[0] = 0;
-// return PolynomialP<N + 1>{coefs};
-// }
-
-// PUGS_INLINE
-// constexpr friend std::ostream&
-// operator<<(std::ostream& os, const PolynomialP<N>& P)
-// {
-// // os << "P(x) = ";
-// bool all_coef_zero = true;
-// if (N == 0) {
-// os << P.coefficients()[0];
-// return os;
-// }
-// if (N != 1) {
-// if (P.coefficients()[N] != 0.) {
-// if (P.coefficients()[N] < 0.) {
-// os << "- ";
-// }
-// if (P.coefficients()[N] != 1 && P.coefficients()[N] != -1) {
-// os << std::abs(P.coefficients()[N]);
-// }
-// os << "x^" << N;
-// all_coef_zero = false;
-// }
-// }
-// for (size_t i = N - 1; i > 1; --i) {
-// if (P.coefficients()[i] != 0.) {
-// if (P.coefficients()[i] > 0.) {
-// os << " + ";
-// } else if (P.coefficients()[i] < 0.) {
-// os << " - ";
-// }
-// if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) {
-// os << std::abs(P.coefficients()[i]);
-// }
-// os << "x^" << i;
-// all_coef_zero = false;
-// }
-// }
-// if (P.coefficients()[1] != 0.) {
-// if (P.coefficients()[1] > 0. && N != 1) {
-// os << " + ";
-// } else if (P.coefficients()[1] < 0.) {
-// os << " - ";
-// }
-// if (P.coefficients()[1] != 1 && P.coefficients()[1] != -1) {
-// os << std::abs(P.coefficients()[1]);
-// }
-// os << "x";
-// all_coef_zero = false;
-// }
-// if (P.coefficients()[0] != 0. || all_coef_zero) {
-// if (P.coefficients()[0] > 0.) {
-// os << " + ";
-// } else if (P.coefficients()[0] < 0.) {
-// os << " - ";
-// }
-// os << std::abs(P.coefficients()[0]);
-// }
-// return os;
-// }
-
-// PUGS_INLINE
-// constexpr friend void
-// lagrangeBasis(const TinyVector<N + 1> zeros, TinyVector<N + 1, PolynomialP<N>>& basis)
-// {
-// PolynomialP<N> lj;
-// for (size_t j = 0; j < N + 1; ++j) {
-// basis[j] = lagrangePolynomialP(zeros, j);
-// }
-// }
-
-// PUGS_INLINE
-// constexpr friend PolynomialP<N>
-// lagrangeToCanonical(const TinyVector<N + 1> lagrange_coefs, const TinyVector<N + 1, PolynomialP<N>>& basis)
-// {
-// PolynomialP<N> P(zero);
-// // lagrangeBasis({0, 0, 0}, basis);
-// for (size_t j = 0; j < N + 1; ++j) {
-// P += basis[j] * lagrange_coefs[j];
-// }
-// return P;
-// }
-
-// template <size_t N>
-// PUGS_INLINE constexpr PolynomialP<N> lagrangePolynomialP(const TinyVector<N + 1> zeros, const size_t k);
-
-// template <size_t N>
-// PUGS_INLINE constexpr TinyVector<N, PolynomialP<N - 1>>
-// lagrangeBasis(const TinyVector<N>& zeros)
-// {
-// static_assert(N >= 1, "invalid degree");
-// TinyVector<N, PolynomialP<N - 1>> basis;
-// PolynomialP<N - 1> lj;
-// for (size_t j = 0; j < N; ++j) {
-// basis[j] = lagrangePolynomialP<N - 1>(zeros, j);
-// }
-// return basis;
-// }
-
-// template <size_t N>
-// PUGS_INLINE constexpr double
-// integrate(const PolynomialP<N>& P, const double& xinf, const double& xsup)
-// {
-// PolynomialP<N + 1> Q = primitive(P);
-// return (Q(xsup) - Q(xinf));
-// }
-
-// template <size_t N>
-// PUGS_INLINE constexpr double
-// symmetricIntegrate(const PolynomialP<N>& P, const double& delta)
-// {
-// Assert(delta > 0, "A positive delta is needed for symmetricIntegrate");
-// double integral = 0.;
-// for (size_t i = 0; i <= N; ++i) {
-// if (i % 2 == 0)
-// integral += 2. * P.coefficients()[i] * std::pow(delta, i + 1) / (i + 1);
-// }
-// return integral;
-// }
-
-// template <size_t N>
-// PUGS_INLINE constexpr PolynomialP<N>
-// lagrangePolynomialP(const TinyVector<N + 1> zeros, const size_t k)
-// {
-// for (size_t i = 0; i < N; ++i) {
-// Assert(zeros[i] < zeros[i + 1], "Interpolation values must be strictly increasing in Lagrange polynomialPs");
-// }
-// PolynomialP<N> lk;
-// lk.coefficients() = zero;
-// lk.coefficients()[0] = 1;
-// for (size_t i = 0; i < N + 1; ++i) {
-// if (k == i)
-// continue;
-// double factor = 1. / (zeros[k] - zeros[i]);
-// PolynomialP<1> P({-zeros[i] * factor, factor});
-// lk *= P;
-// }
-// return lk;
-// }
-
#endif // POLYNOMIALP_HPP