16 #ifndef dealii_polynomial_h
17 #define dealii_polynomial_h
30 #include <shared_mutex>
65 template <
typename number>
92 const unsigned int evaluation_point);
122 value(
const number x, std::vector<number> &
values)
const;
142 template <
typename Number2>
145 const unsigned int n_derivatives,
160 template <std::
size_t n_entries,
typename Number2>
163 const unsigned int n_derivatives,
164 std::array<Number2, n_entries> *
values)
const;
182 scale(
const number factor);
199 template <
typename number2>
201 shift(
const number2 offset);
250 print(std::ostream &out)
const;
257 template <
class Archive>
277 template <
typename number2>
329 template <
typename number>
337 Monomial(
const unsigned int n,
const double coefficient = 1.);
345 static std::vector<Polynomial<number>>
352 static std::vector<number>
353 make_vector(
unsigned int n,
const double coefficient);
390 static std::vector<Polynomial<double>>
400 const unsigned int support_point,
401 std::vector<double> &a);
412 std::vector<Polynomial<double>>
444 static std::vector<Polynomial<double>>
474 Lobatto(
const unsigned int p = 0);
480 static std::vector<Polynomial<double>>
549 static std::vector<Polynomial<double>>
563 static const std::vector<double> &
574 static std::vector<std::unique_ptr<const std::vector<double>>>
627 static std::vector<Polynomial<double>>
742 const unsigned int index);
748 static std::vector<Polynomial<double>>
763 template <
typename Number>
783 template <
typename Number>
797 template <
typename number>
799 : in_lagrange_product_form(false)
800 , lagrange_weight(1.)
805 template <
typename number>
809 if (in_lagrange_product_form ==
true)
811 return lagrange_support_points.size();
816 return coefficients.size() - 1;
822 template <
typename number>
826 if (in_lagrange_product_form ==
false)
831 const unsigned int m = coefficients.size();
832 number value = coefficients.back();
833 for (
int k = m - 2; k >= 0; --k)
834 value = value * x + coefficients[k];
840 const unsigned int m = lagrange_support_points.size();
842 for (
unsigned int j = 0; j < m; ++j)
843 value *= x - lagrange_support_points[j];
844 value *= lagrange_weight;
851 template <
typename number>
852 template <
typename Number2>
855 const unsigned int n_derivatives,
858 values_of_array(std::array<Number2, 1ul>{{x}},
860 reinterpret_cast<std::array<Number2, 1ul> *
>(
values));
865 template <
typename number>
866 template <std::
size_t n_entries,
typename Number2>
869 const std::array<Number2, n_entries> &x,
870 const unsigned int n_derivatives,
871 std::array<Number2, n_entries> *
values)
const
874 if (in_lagrange_product_form ==
true)
879 const unsigned int n_supp = lagrange_support_points.size();
880 const number weight = lagrange_weight;
881 switch (n_derivatives)
884 for (
unsigned int e = 0;
e < n_entries; ++
e)
886 for (
unsigned int k = 1; k <= n_derivatives; ++k)
887 for (
unsigned int e = 0;
e < n_entries; ++
e)
889 for (
unsigned int i = 0; i < n_supp; ++i)
891 std::array<Number2, n_entries> v = x;
892 for (
unsigned int e = 0;
e < n_entries; ++
e)
893 v[
e] -= lagrange_support_points[i];
901 for (
unsigned int k = n_derivatives; k > 0; --k)
902 for (
unsigned int e = 0;
e < n_entries; ++
e)
904 for (
unsigned int e = 0;
e < n_entries; ++
e)
910 number k_factorial = 2;
911 for (
unsigned int k = 2; k <= n_derivatives; ++k)
913 for (
unsigned int e = 0;
e < n_entries; ++
e)
915 k_factorial *=
static_cast<number
>(k + 1);
927 std::array<Number2, n_entries> value;
928 for (
unsigned int e = 0;
e < n_entries; ++
e)
930 for (
unsigned int i = 0; i < n_supp; ++i)
931 for (
unsigned int e = 0;
e < n_entries; ++
e)
932 value[
e] *= (x[
e] - lagrange_support_points[i]);
934 for (
unsigned int e = 0;
e < n_entries; ++
e)
941 std::array<Number2, n_entries> value, derivative = {};
942 for (
unsigned int e = 0;
e < n_entries; ++
e)
944 for (
unsigned int i = 0; i < n_supp; ++i)
945 for (
unsigned int e = 0;
e < n_entries; ++
e)
947 const Number2 v = x[
e] - lagrange_support_points[i];
948 derivative[
e] = derivative[
e] * v + value[
e];
952 for (
unsigned int e = 0;
e < n_entries; ++
e)
962 std::array<Number2, n_entries> value, derivative = {},
964 for (
unsigned int e = 0;
e < n_entries; ++
e)
966 for (
unsigned int i = 0; i < n_supp; ++i)
967 for (
unsigned int e = 0;
e < n_entries; ++
e)
969 const Number2 v = x[
e] - lagrange_support_points[i];
971 derivative[
e] = derivative[
e] * v + value[
e];
975 for (
unsigned int e = 0;
e < n_entries; ++
e)
991 const unsigned int m = coefficients.size();
992 std::vector<std::array<Number2, n_entries>> a(coefficients.size());
993 for (
unsigned int i = 0; i < coefficients.size(); ++i)
994 for (
unsigned int e = 0;
e < n_entries; ++
e)
995 a[i][
e] = coefficients[i];
997 unsigned int j_factorial = 1;
1002 const unsigned int min_valuessize_m =
std::min(n_derivatives + 1, m);
1003 for (
unsigned int j = 0; j < min_valuessize_m; ++j)
1005 for (
int k = m - 2; k >=
static_cast<int>(j); --k)
1006 for (
unsigned int e = 0;
e < n_entries; ++
e)
1007 a[k][
e] += x[
e] * a[k + 1][
e];
1008 for (
unsigned int e = 0;
e < n_entries; ++
e)
1009 values[j][
e] =
static_cast<number
>(j_factorial) * a[j][
e];
1011 j_factorial *= j + 1;
1015 for (
unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
1016 for (
unsigned int e = 0;
e < n_entries; ++
e)
1022 template <
typename number>
1023 template <
class Archive>
1030 ar &in_lagrange_product_form;
1031 ar &lagrange_support_points;
1032 ar &lagrange_weight;
1037 template <
typename Number>
1044 Assert(alpha >= 0 && beta >= 0,
1051 const Number xeval = Number(-1) + 2. * x;
1057 p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
1061 for (
unsigned int i = 1; i < degree; ++i)
1063 const Number v = 2 * i + (alpha + beta);
1064 const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
1065 const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
1066 const Number a3 = v * (v + 1) * (v + 2);
1067 const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
1069 const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
1078 template <
typename Number>
1084 std::vector<Number> x(degree, 0.5);
1095 const Number tolerance =
1105 const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
1106 for (
unsigned int k = 0; k < n_points; ++k)
1110 Number r = 0.5 - 0.5 * std::cos(
static_cast<Number
>(2 * k + 1) /
1113 r = (r + x[k - 1]) / 2;
1116 for (
unsigned int it = 1; it < 1000; ++it)
1119 for (
unsigned int i = 0; i < k; ++i)
1120 s += 1. / (r - x[i]);
1124 (alpha + beta + degree + 1) *
1129 const Number delta = f / (f * s - J_x);
1132 std::abs(delta) < tolerance)
1137 if (it == converged + 1)
1142 ExcMessage(
"Newton iteration for zero of Jacobi polynomial "
1143 "did not converge."));
1149 for (
unsigned int k = n_points; k < degree; ++k)
1150 x[k] = 1.0 - x[degree - k - 1];
HermiteInterpolation(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
HermiteLikeInterpolation(const unsigned int degree, const unsigned int index)
static std::vector< std::unique_ptr< const std::vector< double > > > recursive_coefficients
Hierarchical(const unsigned int p)
static void compute_coefficients(const unsigned int p)
static const std::vector< double > & get_coefficients(const unsigned int p)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::shared_mutex coefficients_lock
static void compute_coefficients(const unsigned int n, const unsigned int support_point, std::vector< double > &a)
LagrangeEquidistant(const unsigned int n, const unsigned int support_point)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Legendre(const unsigned int p)
std::vector< double > compute_coefficients(const unsigned int p)
Lobatto(const unsigned int p=0)
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int p)
static std::vector< Polynomial< number > > generate_complete_basis(const unsigned int degree)
Monomial(const unsigned int n, const double coefficient=1.)
static std::vector< number > make_vector(unsigned int n, const double coefficient)
number value(const number x) const
bool operator==(const Polynomial< number > &p) const
std::vector< number > coefficients
Polynomial< number > primitive() const
Polynomial< number > & operator+=(const Polynomial< number > &p)
void values_of_array(const std::array< Number2, n_entries > &points, const unsigned int n_derivatives, std::array< Number2, n_entries > *values) const
Polynomial< number > derivative() const
void transform_into_standard_form()
void scale(const number factor)
Polynomial< number > & operator-=(const Polynomial< number > &p)
std::vector< number > lagrange_support_points
void shift(const number2 offset)
void print(std::ostream &out) const
bool in_lagrange_product_form
void serialize(Archive &ar, const unsigned int version)
static void multiply(std::vector< number > &coefficients, const number factor)
void value(const Number2 x, const unsigned int n_derivatives, Number2 *values) const
Polynomial< number > & operator*=(const double s)
virtual std::size_t memory_consumption() const
unsigned int degree() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
std::vector< Number > jacobi_polynomial_roots(const unsigned int degree, const int alpha, const int beta)
static constexpr double PI
static const unsigned int invalid_unsigned_int