17#ifndef dealii_simplex_barycentric_polynomials_h
18#define dealii_simplex_barycentric_polynomials_h
80template <
int dim,
typename Number =
double>
93 const Number coefficient);
108 print(std::ostream &out)
const;
125 template <
typename Number2>
132 template <
typename Number2>
139 template <
typename Number2>
145 template <
typename Number2>
177 derivative(
const unsigned int coordinate)
const;
244 std::vector<double> &
values,
306 name()
const override;
311 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
312 clone()
const override;
315 std::vector<BarycentricPolynomial<dim>>
polys;
331template <
int dim,
typename Number1,
typename Number2>
335 return bp * Number1(a);
341template <
int dim,
typename Number1,
typename Number2>
345 return bp + Number1(a);
351template <
int dim,
typename Number1,
typename Number2>
355 return bp - Number1(a);
361template <
int dim,
typename Number>
372template <
int dim,
typename Number>
376 for (
unsigned int d = 0;
d < dim + 1; ++
d)
378 coefficients.reinit(extents);
385template <
int dim,
typename Number>
388 const Number coefficient)
391 for (
unsigned int d = 0;
d < dim + 1; ++
d)
392 extents[
d] = powers[
d] + 1;
393 coefficients.reinit(extents);
395 coefficients(powers) = coefficient;
400template <
int dim,
typename Number>
412template <
int dim,
typename Number>
416 const auto &coeffs = this->coefficients;
417 auto first = index_to_indices(0, coeffs.size());
418 bool print_plus =
false;
419 if (coeffs(
first) != Number())
421 out << coeffs(
first);
424 for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
426 const auto indices = index_to_indices(i, coeffs.size());
427 if (coeffs(indices) == Number())
431 out << coeffs(indices);
432 for (
unsigned int d = 0;
d < dim + 1; ++
d)
435 out <<
" * t" <<
d <<
'^' << indices[
d];
446template <
int dim,
typename Number>
450 auto deg = coefficients.size();
451 for (
unsigned int d = 0;
d < dim + 1; ++
d)
458template <
int dim,
typename Number>
462 return *
this * Number(-1);
467template <
int dim,
typename Number>
468template <
typename Number2>
480template <
int dim,
typename Number>
481template <
typename Number2>
490template <
int dim,
typename Number>
491template <
typename Number2>
501 for (std::size_t i = 0; i < result.
coefficients.n_elements(); ++i)
503 const auto index = index_to_indices(i, result.
coefficients.size());
512template <
int dim,
typename Number>
513template <
typename Number2>
518 return *
this * (Number(1) / Number(a));
523template <
int dim,
typename Number>
529 for (
unsigned int d = 0;
d < dim + 1; ++
d)
537 for (std::size_t i = 0; i < in.n_elements(); ++i)
539 const auto index = index_to_indices(i, in.size());
544 add_coefficients(this->coefficients);
551template <
int dim,
typename Number>
556 return *
this + (-augend);
561template <
int dim,
typename Number>
566 for (
unsigned int d = 0;
d < dim + 1; ++
d)
568 deg[
d] = multiplicand.
degrees()[
d] + degrees()[
d];
573 const auto &coef_1 = this->coefficients;
577 for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
579 const auto index_1 = index_to_indices(i1, coef_1.size());
580 for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
582 const auto index_2 = index_to_indices(i2, coef_2.size());
585 for (
unsigned int d = 0;
d < dim + 1; ++
d)
586 index_out[
d] = index_1[
d] + index_2[
d];
587 coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
596template <
int dim,
typename Number>
599 const unsigned int coordinate)
const
603 if (degrees()[coordinate] == 0)
606 auto deg = degrees();
607 deg[coordinate] -= 1;
610 const auto & coeffs_in = coefficients;
612 for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
614 const auto out_index = index_to_indices(i, coeffs_out.size());
615 auto input_index = out_index;
616 input_index[coordinate] += 1;
618 coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
626template <
int dim,
typename Number>
629 const unsigned int coordinate)
const
632 return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
637template <
int dim,
typename Number>
647 std::array<Number, dim + 1> b_point;
649 for (
unsigned int d = 0;
d < dim; ++
d)
656 for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
658 const auto indices = index_to_indices(i, coefficients.size());
659 const auto coef = coefficients(indices);
660 if (coef == Number())
663 auto temp = Number(1);
664 for (
unsigned int d = 0;
d < dim + 1; ++
d)
666 result += coef * temp;
672template <
int dim,
typename Number>
676 return coefficients.memory_consumption();
679template <
int dim,
typename Number>
682 const std::size_t & index,
688 for (
unsigned int n = 0; n < dim + 1; ++n)
690 std::size_t slice_size = 1;
691 for (
unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
692 slice_size *= extent[n2];
693 result[n] = temp / slice_size;
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
BarycentricPolynomial< dim, Number > operator*(const Number2 &a) const
std::size_t memory_consumption() const
BarycentricPolynomial< dim, Number > operator+(const Number2 &a) const
TableIndices< dim+1 > degrees() const
Table< dim+1, Number > coefficients
BarycentricPolynomial< dim, Number > barycentric_derivative(const unsigned int coordinate) const
Number value(const Point< dim > &point) const
BarycentricPolynomial< dim, Number > operator/(const Number2 &a) const
static TableIndices< dim+1 > index_to_indices(const std::size_t &index, const TableIndices< dim+1 > &extent)
BarycentricPolynomial< dim, Number > operator-() const
void print(std::ostream &out) const
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
BarycentricPolynomial< dim, Number > derivative(const unsigned int coordinate) const
Table< 2, BarycentricPolynomial< dim > > poly_grads
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
BarycentricPolynomials(const std::vector< BarycentricPolynomial< dim > > &polynomials)
virtual std::size_t memory_consumption() const override
std::vector< BarycentricPolynomial< dim > > polys
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
std::string name() const override
static const unsigned int dimension
Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
Table< 3, BarycentricPolynomial< dim > > poly_hessians
Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim > > &grads, std::vector< Tensor< 2, dim > > &grad_grads, std::vector< Tensor< 3, dim > > &third_derivatives, std::vector< Tensor< 4, dim > > &fourth_derivatives) const override
double compute_value(const unsigned int i, const Point< dim > &p) const override
const BarycentricPolynomial< dim > & operator[](const std::size_t i) const
Table< 4, BarycentricPolynomial< dim > > poly_third_derivatives
Table< 5, BarycentricPolynomial< dim > > poly_fourth_derivatives
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
virtual unsigned int degree() const
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)