Reference documentation for deal.II version 9.5.0
Searching...
No Matches
polynomials_barycentric.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2021 - 2023 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
26#include <limits>
27
29
82template <int dim, typename Number = double>
84{
85public:
90
95 const Number coefficient);
96
101 monomial(const unsigned int d);
102
109 void
110 print(std::ostream &out) const;
111
116 degrees() const;
117
122 operator-() const;
123
127 template <typename Number2>
129 operator+(const Number2 &a) const;
130
134 template <typename Number2>
136 operator-(const Number2 &a) const;
137
141 template <typename Number2>
143 operator*(const Number2 &a) const;
144
148 template <typename Number2>
150 operator/(const Number2 &a) const;
151
157
163
168 operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
169
174 barycentric_derivative(const unsigned int coordinate) const;
175
180 derivative(const unsigned int coordinate) const;
181
185 Number
186 value(const Point<dim> &point) const;
187
191 std::size_t
192 memory_consumption() const;
193
194protected:
199
209 index_to_indices(const std::size_t & index,
210 const TableIndices<dim + 1> &extent);
211};
212
216template <int dim>
218{
219public:
224
228 using GradType = std::array<PolyType, dim>;
229
233 using HessianType = std::array<GradType, dim>;
234
238 using ThirdDerivativesType = std::array<HessianType, dim>;
239
243 using FourthDerivativesType = std::array<ThirdDerivativesType, dim>;
244
248 static constexpr unsigned int dimension = dim;
249
254 get_fe_p_basis(const unsigned int degree);
255
260 const std::vector<BarycentricPolynomial<dim>> &polynomials);
261
266 operator[](const std::size_t i) const;
267
271 void
272 evaluate(const Point<dim> & unit_point,
273 std::vector<double> & values,
276 std::vector<Tensor<3, dim>> &third_derivatives,
277 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
278
282 double
283 compute_value(const unsigned int i, const Point<dim> &p) const override;
284
289 compute_1st_derivative(const unsigned int i,
290 const Point<dim> & p) const override;
291
296 compute_2nd_derivative(const unsigned int i,
297 const Point<dim> & p) const override;
298
303 compute_3rd_derivative(const unsigned int i,
304 const Point<dim> & p) const override;
305
310 compute_4th_derivative(const unsigned int i,
311 const Point<dim> & p) const override;
312
317 compute_grad(const unsigned int i, const Point<dim> &p) const override;
318
324
328 virtual std::size_t
329 memory_consumption() const override;
330
334 std::string
335 name() const override;
336
340 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
341 clone() const override;
342
343protected:
344 std::vector<PolyType> polys;
346 std::vector<HessianType> poly_hessians;
347 std::vector<ThirdDerivativesType> poly_third_derivatives;
348 std::vector<FourthDerivativesType> poly_fourth_derivatives;
349};
350
351// non-member template functions for algebra
352
356template <int dim, typename Number1, typename Number2>
359{
360 return bp * Number1(a);
361}
362
366template <int dim, typename Number1, typename Number2>
369{
370 return bp + Number1(a);
371}
372
376template <int dim, typename Number1, typename Number2>
379{
380 return bp - Number1(a);
381}
382
386template <int dim, typename Number>
387std::ostream &
388operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
389{
390 bp.print(out);
391 return out;
392}
393
394// Template function definitions
395
396// BarycentricPolynomial:
397template <int dim, typename Number>
399{
400 TableIndices<dim + 1> extents;
401 for (unsigned int d = 0; d < dim + 1; ++d)
402 extents[d] = 1;
403 coefficients.reinit(extents);
404
405 coefficients(TableIndices<dim + 1>{}) = Number();
406}
407
408
409
410template <int dim, typename Number>
412 const TableIndices<dim + 1> &powers,
413 const Number coefficient)
414{
415 TableIndices<dim + 1> extents;
416 for (unsigned int d = 0; d < dim + 1; ++d)
417 extents[d] = powers[d] + 1;
418 coefficients.reinit(extents);
419
420 coefficients(powers) = coefficient;
421}
422
423
424
425template <int dim, typename Number>
428{
429 AssertIndexRange(d, dim + 1);
430 TableIndices<dim + 1> indices;
431 indices[d] = 1;
432 return BarycentricPolynomial<dim, Number>(indices, Number(1));
433}
434
435
436
437template <int dim, typename Number>
438void
440{
441 const auto &coeffs = this->coefficients;
442 auto first = index_to_indices(0, coeffs.size());
443 bool print_plus = false;
444 if (coeffs(first) != Number())
445 {
446 out << coeffs(first);
447 print_plus = true;
448 }
449 for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
450 {
451 const auto indices = index_to_indices(i, coeffs.size());
452 if (coeffs(indices) == Number())
453 continue;
454 if (print_plus)
455 out << " + ";
456 out << coeffs(indices);
457 for (unsigned int d = 0; d < dim + 1; ++d)
458 {
459 if (indices[d] != 0)
460 out << " * t" << d << '^' << indices[d];
461 }
462 print_plus = true;
463 }
464
465 if (!print_plus)
466 out << Number();
467}
468
469
470
471template <int dim, typename Number>
474{
475 auto deg = coefficients.size();
476 for (unsigned int d = 0; d < dim + 1; ++d)
477 deg[d] -= 1;
478 return deg;
479}
480
481
482
483template <int dim, typename Number>
486{
487 return *this * Number(-1);
488}
489
490
491
492template <int dim, typename Number>
493template <typename Number2>
496{
498 result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
499
500 return result;
501}
502
503
504
505template <int dim, typename Number>
506template <typename Number2>
509{
510 return *this + (-a);
511}
512
513
514
515template <int dim, typename Number>
516template <typename Number2>
519{
520 if (a == Number2())
521 {
523 }
524
526 for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
527 {
528 const auto index = index_to_indices(i, result.coefficients.size());
529 result.coefficients(index) *= a;
530 }
531
532 return result;
533}
534
535
536
537template <int dim, typename Number>
538template <typename Number2>
541{
542 Assert(a != Number2(), ExcDivideByZero());
543 return *this * (Number(1) / Number(a));
544}
545
546
547
548template <int dim, typename Number>
551 const BarycentricPolynomial<dim, Number> &augend) const
552{
554 for (unsigned int d = 0; d < dim + 1; ++d)
555 {
556 deg[d] = std::max(degrees()[d], augend.degrees()[d]);
557 }
558
559 BarycentricPolynomial<dim, Number> result(deg, Number());
560
561 auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
562 for (std::size_t i = 0; i < in.n_elements(); ++i)
563 {
564 const auto index = index_to_indices(i, in.size());
565 result.coefficients(index) += in(index);
566 }
567 };
568
571 return result;
572}
573
574
575
576template <int dim, typename Number>
579 const BarycentricPolynomial<dim, Number> &augend) const
580{
581 return *this + (-augend);
582}
583
584
585
586template <int dim, typename Number>
589 const BarycentricPolynomial<dim, Number> &multiplicand) const
590{
592 for (unsigned int d = 0; d < dim + 1; ++d)
593 {
594 deg[d] = multiplicand.degrees()[d] + degrees()[d];
595 }
596
597 BarycentricPolynomial<dim, Number> result(deg, Number());
598
599 const auto &coef_1 = this->coefficients;
600 const auto &coef_2 = multiplicand.coefficients;
601 auto & coef_out = result.coefficients;
602
603 for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
604 {
605 const auto index_1 = index_to_indices(i1, coef_1.size());
606 for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
607 {
608 const auto index_2 = index_to_indices(i2, coef_2.size());
609
610 TableIndices<dim + 1> index_out;
611 for (unsigned int d = 0; d < dim + 1; ++d)
612 index_out[d] = index_1[d] + index_2[d];
613 coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
614 }
615 }
616
617 return result;
618}
619
620
621
622template <int dim, typename Number>
625 const unsigned int coordinate) const
626{
627 AssertIndexRange(coordinate, dim + 1);
628
629 if (degrees()[coordinate] == 0)
631
632 auto deg = degrees();
633 deg[coordinate] -= 1;
635 std::numeric_limits<Number>::max());
636 const auto & coeffs_in = coefficients;
637 auto & coeffs_out = result.coefficients;
638 for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
639 {
640 const auto out_index = index_to_indices(i, coeffs_out.size());
641 auto input_index = out_index;
642 input_index[coordinate] += 1;
643
644 coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
645 }
646
647 return result;
648}
649
650
651
652template <int dim, typename Number>
655 const unsigned int coordinate) const
656{
657 AssertIndexRange(coordinate, dim);
658 return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
659}
660
661
662
663template <int dim, typename Number>
664Number
666{
667 // TODO: this is probably not numerically stable for higher order.
668 // We really need some version of Horner's method.
669 Number result = {};
670
671 // Begin by converting point (which is in Cartesian coordinates) to
672 // barycentric coordinates:
673 std::array<Number, dim + 1> b_point;
674 b_point[0] = 1.0;
675 for (unsigned int d = 0; d < dim; ++d)
676 {
677 b_point[0] -= point[d];
678 b_point[d + 1] = point[d];
679 }
680
681 // Now evaluate the polynomial at the computed barycentric point:
682 for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
683 {
684 const auto indices = index_to_indices(i, coefficients.size());
685 const auto coef = coefficients(indices);
686 if (coef == Number())
687 continue;
688
689 auto temp = Number(1);
690 for (unsigned int d = 0; d < dim + 1; ++d)
691 temp *= std::pow(b_point[d], indices[d]);
692 result += coef * temp;
693 }
694
695 return result;
696}
697
698template <int dim, typename Number>
699std::size_t
701{
702 return coefficients.memory_consumption();
703}
704
705template <int dim, typename Number>
708 const std::size_t & index,
709 const TableIndices<dim + 1> &extent)
710{
712 auto temp = index;
713
714 for (unsigned int n = 0; n < dim + 1; ++n)
715 {
716 std::size_t slice_size = 1;
717 for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
718 slice_size *= extent[n2];
719 result[n] = temp / slice_size;
720 temp %= slice_size;
721 }
722 return result;
723}
724
725template <int dim>
728{
729 AssertIndexRange(i, polys.size());
730 return polys[i];
731}
732
734
735#endif
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
std::array< HessianType, dim > ThirdDerivativesType
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::size_t memory_consumption() const override
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
std::array< ThirdDerivativesType, dim > FourthDerivativesType
std::string name() const override
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
Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
static constexpr unsigned int dimension
std::vector< PolyType > polys
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
std::vector< ThirdDerivativesType > poly_third_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
std::vector< HessianType > poly_hessians
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
std::vector< FourthDerivativesType > poly_fourth_derivatives
Definition point.h:112
virtual unsigned int degree() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > first
Definition grid_out.cc:4615
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
std::ostream & operator<<(std::ostream &out, const BarycentricPolynomial< dim, Number > &bp)
BarycentricPolynomial< dim, Number1 > operator*(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)