Reference documentation for deal.II version 9.3.3
polynomials_barycentric.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_simplex_barycentric_polynomials_h
18#define dealii_simplex_barycentric_polynomials_h
19
20#include <deal.II/base/config.h>
21
24#include <deal.II/base/table.h>
25
27
80template <int dim, typename Number = double>
82{
83public:
88
93 const Number coefficient);
94
99 monomial(const unsigned int d);
100
107 void
108 print(std::ostream &out) const;
109
114 degrees() const;
115
120 operator-() const;
121
125 template <typename Number2>
127 operator+(const Number2 &a) const;
128
132 template <typename Number2>
134 operator-(const Number2 &a) const;
135
139 template <typename Number2>
140 BarycentricPolynomial<dim, Number> operator*(const Number2 &a) const;
141
145 template <typename Number2>
147 operator/(const Number2 &a) const;
148
154
160
165 operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
166
171 barycentric_derivative(const unsigned int coordinate) const;
172
177 derivative(const unsigned int coordinate) const;
178
182 Number
183 value(const Point<dim> &point) const;
184
188 std::size_t
189 memory_consumption() const;
190
191protected:
196
206 index_to_indices(const std::size_t & index,
207 const TableIndices<dim + 1> &extent);
208};
209
213template <int dim>
215{
216public:
220 static const unsigned int dimension = dim;
221
226 get_fe_p_basis(const unsigned int degree);
227
232 const std::vector<BarycentricPolynomial<dim>> &polynomials);
233
237 const BarycentricPolynomial<dim> &operator[](const std::size_t i) const;
238
242 void
243 evaluate(const Point<dim> & unit_point,
244 std::vector<double> & values,
247 std::vector<Tensor<3, dim>> &third_derivatives,
248 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
249
253 double
254 compute_value(const unsigned int i, const Point<dim> &p) const override;
255
260 compute_1st_derivative(const unsigned int i,
261 const Point<dim> & p) const override;
262
267 compute_2nd_derivative(const unsigned int i,
268 const Point<dim> & p) const override;
269
274 compute_3rd_derivative(const unsigned int i,
275 const Point<dim> & p) const override;
276
281 compute_4th_derivative(const unsigned int i,
282 const Point<dim> & p) const override;
283
288 compute_grad(const unsigned int i, const Point<dim> &p) const override;
289
295
299 virtual std::size_t
300 memory_consumption() const override;
301
305 std::string
306 name() const override;
307
311 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
312 clone() const override;
313
314protected:
315 std::vector<BarycentricPolynomial<dim>> polys;
316
318
320
322
324};
325
326// non-member template functions for algebra
327
331template <int dim, typename Number1, typename Number2>
334{
335 return bp * Number1(a);
336}
337
341template <int dim, typename Number1, typename Number2>
344{
345 return bp + Number1(a);
346}
347
351template <int dim, typename Number1, typename Number2>
354{
355 return bp - Number1(a);
356}
357
361template <int dim, typename Number>
362std::ostream &
363operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
364{
365 bp.print(out);
366 return out;
367}
368
369// Template function definitions
370
371// BarycentricPolynomial:
372template <int dim, typename Number>
374{
375 TableIndices<dim + 1> extents;
376 for (unsigned int d = 0; d < dim + 1; ++d)
377 extents[d] = 1;
378 coefficients.reinit(extents);
379
380 coefficients(TableIndices<dim + 1>{}) = Number();
381}
382
383
384
385template <int dim, typename Number>
387 const TableIndices<dim + 1> &powers,
388 const Number coefficient)
389{
390 TableIndices<dim + 1> extents;
391 for (unsigned int d = 0; d < dim + 1; ++d)
392 extents[d] = powers[d] + 1;
393 coefficients.reinit(extents);
394
395 coefficients(powers) = coefficient;
396}
397
398
399
400template <int dim, typename Number>
403{
404 AssertIndexRange(d, dim + 1);
405 TableIndices<dim + 1> indices;
406 indices[d] = 1;
407 return BarycentricPolynomial<dim, Number>(indices, Number(1));
408}
409
410
411
412template <int dim, typename Number>
413void
415{
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())
420 {
421 out << coeffs(first);
422 print_plus = true;
423 }
424 for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
425 {
426 const auto indices = index_to_indices(i, coeffs.size());
427 if (coeffs(indices) == Number())
428 continue;
429 if (print_plus)
430 out << " + ";
431 out << coeffs(indices);
432 for (unsigned int d = 0; d < dim + 1; ++d)
433 {
434 if (indices[d] != 0)
435 out << " * t" << d << '^' << indices[d];
436 }
437 print_plus = true;
438 }
439
440 if (!print_plus)
441 out << Number();
442}
443
444
445
446template <int dim, typename Number>
449{
450 auto deg = coefficients.size();
451 for (unsigned int d = 0; d < dim + 1; ++d)
452 deg[d] -= 1;
453 return deg;
454}
455
456
457
458template <int dim, typename Number>
461{
462 return *this * Number(-1);
463}
464
465
466
467template <int dim, typename Number>
468template <typename Number2>
471{
473 result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
474
475 return result;
476}
477
478
479
480template <int dim, typename Number>
481template <typename Number2>
484{
485 return *this + (-a);
486}
487
488
489
490template <int dim, typename Number>
491template <typename Number2>
493 operator*(const Number2 &a) const
494{
495 if (a == Number2())
496 {
498 }
499
501 for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
502 {
503 const auto index = index_to_indices(i, result.coefficients.size());
504 result.coefficients(index) *= a;
505 }
506
507 return result;
508}
509
510
511
512template <int dim, typename Number>
513template <typename Number2>
516{
517 Assert(a != Number2(), ExcDivideByZero());
518 return *this * (Number(1) / Number(a));
519}
520
521
522
523template <int dim, typename Number>
527{
529 for (unsigned int d = 0; d < dim + 1; ++d)
530 {
531 deg[d] = std::max(degrees()[d], augend.degrees()[d]);
532 }
533
534 BarycentricPolynomial<dim, Number> result(deg, Number());
535
536 auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
537 for (std::size_t i = 0; i < in.n_elements(); ++i)
538 {
539 const auto index = index_to_indices(i, in.size());
540 result.coefficients(index) += in(index);
541 }
542 };
543
546 return result;
547}
548
549
550
551template <int dim, typename Number>
555{
556 return *this + (-augend);
557}
558
559
560
561template <int dim, typename Number>
563 operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const
564{
566 for (unsigned int d = 0; d < dim + 1; ++d)
567 {
568 deg[d] = multiplicand.degrees()[d] + degrees()[d];
569 }
570
571 BarycentricPolynomial<dim, Number> result(deg, Number());
572
573 const auto &coef_1 = this->coefficients;
574 const auto &coef_2 = multiplicand.coefficients;
575 auto & coef_out = result.coefficients;
576
577 for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
578 {
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)
581 {
582 const auto index_2 = index_to_indices(i2, coef_2.size());
583
584 TableIndices<dim + 1> index_out;
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);
588 }
589 }
590
591 return result;
592}
593
594
595
596template <int dim, typename Number>
599 const unsigned int coordinate) const
600{
601 AssertIndexRange(coordinate, dim + 1);
602
603 if (degrees()[coordinate] == 0)
605
606 auto deg = degrees();
607 deg[coordinate] -= 1;
610 const auto & coeffs_in = coefficients;
611 auto & coeffs_out = result.coefficients;
612 for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
613 {
614 const auto out_index = index_to_indices(i, coeffs_out.size());
615 auto input_index = out_index;
616 input_index[coordinate] += 1;
617
618 coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
619 }
620
621 return result;
622}
623
624
625
626template <int dim, typename Number>
629 const unsigned int coordinate) const
630{
631 AssertIndexRange(coordinate, dim);
632 return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
633}
634
635
636
637template <int dim, typename Number>
638Number
640{
641 // TODO: this is probably not numerically stable for higher order.
642 // We really need some version of Horner's method.
643 Number result = {};
644
645 // Begin by converting point (which is in Cartesian coordinates) to
646 // barycentric coordinates:
647 std::array<Number, dim + 1> b_point;
648 b_point[0] = 1.0;
649 for (unsigned int d = 0; d < dim; ++d)
650 {
651 b_point[0] -= point[d];
652 b_point[d + 1] = point[d];
653 }
654
655 // Now evaluate the polynomial at the computed barycentric point:
656 for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
657 {
658 const auto indices = index_to_indices(i, coefficients.size());
659 const auto coef = coefficients(indices);
660 if (coef == Number())
661 continue;
662
663 auto temp = Number(1);
664 for (unsigned int d = 0; d < dim + 1; ++d)
665 temp *= std::pow(b_point[d], indices[d]);
666 result += coef * temp;
667 }
668
669 return result;
670}
671
672template <int dim, typename Number>
673std::size_t
675{
676 return coefficients.memory_consumption();
677}
678
679template <int dim, typename Number>
682 const std::size_t & index,
683 const TableIndices<dim + 1> &extent)
684{
686 auto temp = index;
687
688 for (unsigned int n = 0; n < dim + 1; ++n)
689 {
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;
694 temp %= slice_size;
695 }
696 return result;
697}
698
699template <int dim>
701 operator[](const std::size_t i) const
702{
703 AssertIndexRange(i, polys.size());
704 return polys[i];
705}
706
708
709#endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
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)
Definition: point.h:111
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
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 2 > first
Definition: grid_out.cc:4587
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcDivideByZero()
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
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)