Skip to content
Snippets Groups Projects
Commit e0d7cfed authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

add features for 3D polynomials except offstream

parent 2cfa61f3
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/analysis/PolynomialP.hpp
+ 147
− 193
View file @ e0d7cfed
......@@ -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,12 +313,50 @@ 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
constexpr friend std::ostream&
......@@ -445,181 +483,97 @@ class PolynomialP
return os;
} else {
throw NotImplementedError("Not yet Available in 3D");
}
}
PUGS_INLINE constexpr PolynomialP() noexcept = default;
~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];
// }
// 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 << " - ";
// }
// for (size_t j = Mr; j <= Nr; ++j) {
// R.coefficients()[j] = 0;
// if (coef != 1 && coef != -1) {
// os << std::abs(coef);
// }
// os << "z^" << N;
// }
// 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);
// 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++;
// }
// coefs[0] = 0;
// return PolynomialP<N + 1>{coefs};
// } else {
// degree--;
// k = degree;
// i = 0;
// j = 0;
// }
// 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.) {
// rel_pos = TinyVector<Dimension, size_t>(i, j, k);
// double coef = P[rel_pos];
// if (coef != 0.) {
// if (coef > 0.) {
// os << " + ";
// } else if (P.coefficients()[i] < 0.) {
// } else if (coef < 0.) {
// os << " - ";
// }
// if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) {
// os << std::abs(P.coefficients()[i]);
// 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;
// all_coef_zero = false;
// }
// }
// if (P.coefficients()[1] != 0.) {
// if (P.coefficients()[1] > 0. && N != 1) {
// os << " + ";
// } else if (P.coefficients()[1] < 0.) {
// os << " - ";
// } 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;
// }
// 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;
// }
//
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");
}
}
// 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;
// }
PUGS_INLINE constexpr PolynomialP() noexcept = default;
~PolynomialP() = default;
};
#endif // POLYNOMIALP_HPP
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment