Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
polynomials_barycentric.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2021 - 2022 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
9// Public License as published by the Free Software Foundation; either
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>
141 operator*(const Number2 &a) const;
142
146 template <typename Number2>
148 operator/(const Number2 &a) const;
149
155
161
166 operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
167
172 barycentric_derivative(const unsigned int coordinate) const;
173
178 derivative(const unsigned int coordinate) const;
179
183 Number
184 value(const Point<dim> &point) const;
185
189 std::size_t
190 memory_consumption() const;
191
192protected:
197
207 index_to_indices(const std::size_t & index,
208 const TableIndices<dim + 1> &extent);
209};
210
214template <int dim>
216{
217public:
222
226 using GradType = std::array<PolyType, dim>;
227
231 using HessianType = std::array<GradType, dim>;
232
236 using ThirdDerivativesType = std::array<HessianType, dim>;
237
241 using FourthDerivativesType = std::array<ThirdDerivativesType, dim>;
242
246 static constexpr unsigned int dimension = dim;
247
252 get_fe_p_basis(const unsigned int degree);
253
258 const std::vector<BarycentricPolynomial<dim>> &polynomials);
259
264 operator[](const std::size_t i) const;
265
269 void
270 evaluate(const Point<dim> & unit_point,
271 std::vector<double> & values,
272 std::vector<Tensor<1, dim>> &grads,
273 std::vector<Tensor<2, dim>> &grad_grads,
274 std::vector<Tensor<3, dim>> &third_derivatives,
275 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
276
280 double
281 compute_value(const unsigned int i, const Point<dim> &p) const override;
282
287 compute_1st_derivative(const unsigned int i,
288 const Point<dim> & p) const override;
289
294 compute_2nd_derivative(const unsigned int i,
295 const Point<dim> & p) const override;
296
301 compute_3rd_derivative(const unsigned int i,
302 const Point<dim> & p) const override;
303
308 compute_4th_derivative(const unsigned int i,
309 const Point<dim> & p) const override;
310
315 compute_grad(const unsigned int i, const Point<dim> &p) const override;
316
321 compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
322
326 virtual std::size_t
327 memory_consumption() const override;
328
332 std::string
333 name() const override;
334
338 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
339 clone() const override;
340
341protected:
342 std::vector<PolyType> polys;
343 std::vector<GradType> poly_grads;
344 std::vector<HessianType> poly_hessians;
345 std::vector<ThirdDerivativesType> poly_third_derivatives;
346 std::vector<FourthDerivativesType> poly_fourth_derivatives;
347};
348
349// non-member template functions for algebra
350
354template <int dim, typename Number1, typename Number2>
357{
358 return bp * Number1(a);
359}
360
364template <int dim, typename Number1, typename Number2>
367{
368 return bp + Number1(a);
369}
370
374template <int dim, typename Number1, typename Number2>
377{
378 return bp - Number1(a);
379}
380
384template <int dim, typename Number>
385std::ostream &
386operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
387{
388 bp.print(out);
389 return out;
390}
391
392// Template function definitions
393
394// BarycentricPolynomial:
395template <int dim, typename Number>
397{
398 TableIndices<dim + 1> extents;
399 for (unsigned int d = 0; d < dim + 1; ++d)
400 extents[d] = 1;
401 coefficients.reinit(extents);
402
403 coefficients(TableIndices<dim + 1>{}) = Number();
404}
405
406
407
408template <int dim, typename Number>
410 const TableIndices<dim + 1> &powers,
411 const Number coefficient)
412{
413 TableIndices<dim + 1> extents;
414 for (unsigned int d = 0; d < dim + 1; ++d)
415 extents[d] = powers[d] + 1;
416 coefficients.reinit(extents);
417
418 coefficients(powers) = coefficient;
419}
420
421
422
423template <int dim, typename Number>
426{
427 AssertIndexRange(d, dim + 1);
428 TableIndices<dim + 1> indices;
429 indices[d] = 1;
430 return BarycentricPolynomial<dim, Number>(indices, Number(1));
431}
432
433
434
435template <int dim, typename Number>
436void
438{
439 const auto &coeffs = this->coefficients;
440 auto first = index_to_indices(0, coeffs.size());
441 bool print_plus = false;
442 if (coeffs(first) != Number())
443 {
444 out << coeffs(first);
445 print_plus = true;
446 }
447 for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
448 {
449 const auto indices = index_to_indices(i, coeffs.size());
450 if (coeffs(indices) == Number())
451 continue;
452 if (print_plus)
453 out << " + ";
454 out << coeffs(indices);
455 for (unsigned int d = 0; d < dim + 1; ++d)
456 {
457 if (indices[d] != 0)
458 out << " * t" << d << '^' << indices[d];
459 }
460 print_plus = true;
461 }
462
463 if (!print_plus)
464 out << Number();
465}
466
467
468
469template <int dim, typename Number>
472{
473 auto deg = coefficients.size();
474 for (unsigned int d = 0; d < dim + 1; ++d)
475 deg[d] -= 1;
476 return deg;
477}
478
479
480
481template <int dim, typename Number>
484{
485 return *this * Number(-1);
486}
487
488
489
490template <int dim, typename Number>
491template <typename Number2>
494{
496 result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
497
498 return result;
499}
500
501
502
503template <int dim, typename Number>
504template <typename Number2>
507{
508 return *this + (-a);
509}
510
511
512
513template <int dim, typename Number>
514template <typename Number2>
517{
518 if (a == Number2())
519 {
521 }
522
524 for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
525 {
526 const auto index = index_to_indices(i, result.coefficients.size());
527 result.coefficients(index) *= a;
528 }
529
530 return result;
531}
532
533
534
535template <int dim, typename Number>
536template <typename Number2>
539{
540 Assert(a != Number2(), ExcDivideByZero());
541 return *this * (Number(1) / Number(a));
542}
543
544
545
546template <int dim, typename Number>
549 const BarycentricPolynomial<dim, Number> &augend) const
550{
552 for (unsigned int d = 0; d < dim + 1; ++d)
553 {
554 deg[d] = std::max(degrees()[d], augend.degrees()[d]);
555 }
556
557 BarycentricPolynomial<dim, Number> result(deg, Number());
558
559 auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
560 for (std::size_t i = 0; i < in.n_elements(); ++i)
561 {
562 const auto index = index_to_indices(i, in.size());
563 result.coefficients(index) += in(index);
564 }
565 };
566
567 add_coefficients(this->coefficients);
568 add_coefficients(augend.coefficients);
569 return result;
570}
571
572
573
574template <int dim, typename Number>
577 const BarycentricPolynomial<dim, Number> &augend) const
578{
579 return *this + (-augend);
580}
581
582
583
584template <int dim, typename Number>
587 const BarycentricPolynomial<dim, Number> &multiplicand) const
588{
590 for (unsigned int d = 0; d < dim + 1; ++d)
591 {
592 deg[d] = multiplicand.degrees()[d] + degrees()[d];
593 }
594
595 BarycentricPolynomial<dim, Number> result(deg, Number());
596
597 const auto &coef_1 = this->coefficients;
598 const auto &coef_2 = multiplicand.coefficients;
599 auto & coef_out = result.coefficients;
600
601 for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
602 {
603 const auto index_1 = index_to_indices(i1, coef_1.size());
604 for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
605 {
606 const auto index_2 = index_to_indices(i2, coef_2.size());
607
608 TableIndices<dim + 1> index_out;
609 for (unsigned int d = 0; d < dim + 1; ++d)
610 index_out[d] = index_1[d] + index_2[d];
611 coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
612 }
613 }
614
615 return result;
616}
617
618
619
620template <int dim, typename Number>
623 const unsigned int coordinate) const
624{
625 AssertIndexRange(coordinate, dim + 1);
626
627 if (degrees()[coordinate] == 0)
629
630 auto deg = degrees();
631 deg[coordinate] -= 1;
633 std::numeric_limits<Number>::max());
634 const auto & coeffs_in = coefficients;
635 auto & coeffs_out = result.coefficients;
636 for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
637 {
638 const auto out_index = index_to_indices(i, coeffs_out.size());
639 auto input_index = out_index;
640 input_index[coordinate] += 1;
641
642 coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
643 }
644
645 return result;
646}
647
648
649
650template <int dim, typename Number>
653 const unsigned int coordinate) const
654{
655 AssertIndexRange(coordinate, dim);
656 return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
657}
658
659
660
661template <int dim, typename Number>
662Number
664{
665 // TODO: this is probably not numerically stable for higher order.
666 // We really need some version of Horner's method.
667 Number result = {};
668
669 // Begin by converting point (which is in Cartesian coordinates) to
670 // barycentric coordinates:
671 std::array<Number, dim + 1> b_point;
672 b_point[0] = 1.0;
673 for (unsigned int d = 0; d < dim; ++d)
674 {
675 b_point[0] -= point[d];
676 b_point[d + 1] = point[d];
677 }
678
679 // Now evaluate the polynomial at the computed barycentric point:
680 for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
681 {
682 const auto indices = index_to_indices(i, coefficients.size());
683 const auto coef = coefficients(indices);
684 if (coef == Number())
685 continue;
686
687 auto temp = Number(1);
688 for (unsigned int d = 0; d < dim + 1; ++d)
689 temp *= std::pow(b_point[d], indices[d]);
690 result += coef * temp;
691 }
692
693 return result;
694}
695
696template <int dim, typename Number>
697std::size_t
699{
700 return coefficients.memory_consumption();
701}
702
703template <int dim, typename Number>
706 const std::size_t & index,
707 const TableIndices<dim + 1> &extent)
708{
710 auto temp = index;
711
712 for (unsigned int n = 0; n < dim + 1; ++n)
713 {
714 std::size_t slice_size = 1;
715 for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
716 slice_size *= extent[n2];
717 result[n] = temp / slice_size;
718 temp %= slice_size;
719 }
720 return result;
721}
722
723template <int dim>
726{
727 AssertIndexRange(i, polys.size());
728 return polys[i];
729}
730
732
733#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
std::array< PolyType, dim > GradType
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::vector< GradType > poly_grads
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
std::array< GradType, dim > HessianType
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:111
virtual unsigned int degree() const
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 2 > first
Definition: grid_out.cc:4603
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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)