Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+00:00
Searching...
No Matches
polynomials_barycentric.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_simplex_barycentric_polynomials_h
17#define dealii_simplex_barycentric_polynomials_h
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/table.h>
24
25#include <limits>
26
28
81template <int dim, typename Number = double>
83{
84public:
89
94 const Number coefficient);
95
100 monomial(const unsigned int d);
101
108 void
109 print(std::ostream &out) const;
110
115 degrees() const;
116
121 operator-() const;
122
126 template <typename Number2>
128 operator+(const Number2 &a) const;
129
133 template <typename Number2>
135 operator-(const Number2 &a) const;
136
140 template <typename Number2>
142 operator*(const Number2 &a) const;
143
147 template <typename Number2>
149 operator/(const Number2 &a) const;
150
156
162
167 operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
168
173 barycentric_derivative(const unsigned int coordinate) const;
174
179 derivative(const unsigned int coordinate) const;
180
184 Number
185 value(const Point<dim> &point) const;
186
190 std::size_t
191 memory_consumption() const;
192
193protected:
198
208 index_to_indices(const std::size_t &index,
209 const TableIndices<dim + 1> &extents);
210};
211
215template <int dim>
217{
218public:
223
227 using GradType = std::array<PolyType, dim>;
228
232 using HessianType = std::array<GradType, dim>;
233
237 using ThirdDerivativesType = std::array<HessianType, dim>;
238
242 using FourthDerivativesType = std::array<ThirdDerivativesType, dim>;
243
247 static constexpr unsigned int dimension = dim;
248
253 get_fe_p_basis(const unsigned int degree);
254
259 const std::vector<BarycentricPolynomial<dim>> &polynomials);
260
265 operator[](const std::size_t i) const;
266
270 void
271 evaluate(const Point<dim> &unit_point,
272 std::vector<double> &values,
275 std::vector<Tensor<3, dim>> &third_derivatives,
276 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
277
281 double
282 compute_value(const unsigned int i, const Point<dim> &p) const override;
283
288 compute_1st_derivative(const unsigned int i,
289 const Point<dim> &p) const override;
290
295 compute_2nd_derivative(const unsigned int i,
296 const Point<dim> &p) const override;
297
302 compute_3rd_derivative(const unsigned int i,
303 const Point<dim> &p) const override;
304
309 compute_4th_derivative(const unsigned int i,
310 const Point<dim> &p) const override;
311
316 compute_grad(const unsigned int i, const Point<dim> &p) const override;
317
323
327 virtual std::size_t
328 memory_consumption() const override;
329
333 std::string
334 name() const override;
335
339 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
340 clone() const override;
341
342protected:
343 std::vector<PolyType> polys;
345 std::vector<HessianType> poly_hessians;
346 std::vector<ThirdDerivativesType> poly_third_derivatives;
347 std::vector<FourthDerivativesType> poly_fourth_derivatives;
348};
349
350// non-member template functions for algebra
351
355template <int dim, typename Number1, typename Number2>
358{
359 return bp * Number1(a);
360}
361
365template <int dim, typename Number1, typename Number2>
368{
369 return bp + Number1(a);
370}
371
375template <int dim, typename Number1, typename Number2>
378{
379 return bp - Number1(a);
380}
381
385template <int dim, typename Number>
386std::ostream &
387operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
388{
389 bp.print(out);
390 return out;
391}
392
393// Template function definitions
394
395// BarycentricPolynomial:
396template <int dim, typename Number>
398{
399 TableIndices<dim + 1> extents;
400 for (unsigned int d = 0; d < dim + 1; ++d)
401 extents[d] = 1;
402 coefficients.reinit(extents);
403
404 coefficients(TableIndices<dim + 1>{}) = Number();
405}
406
407
408
409template <int dim, typename Number>
411 const TableIndices<dim + 1> &powers,
412 const Number coefficient)
413{
414 TableIndices<dim + 1> extents;
415 for (unsigned int d = 0; d < dim + 1; ++d)
416 extents[d] = powers[d] + 1;
417 coefficients.reinit(extents);
418
419 coefficients(powers) = coefficient;
420}
421
422
423
424template <int dim, typename Number>
427{
428 AssertIndexRange(d, dim + 1);
429 TableIndices<dim + 1> indices;
430 indices[d] = 1;
431 return BarycentricPolynomial<dim, Number>(indices, Number(1));
432}
433
434
435
436template <int dim, typename Number>
437void
439{
440 const auto &coeffs = this->coefficients;
441 auto first = index_to_indices(0, coeffs.size());
442 bool print_plus = false;
443 if (coeffs(first) != Number())
444 {
445 out << coeffs(first);
446 print_plus = true;
447 }
448 for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
449 {
450 const auto indices = index_to_indices(i, coeffs.size());
451 if (coeffs(indices) == Number())
452 continue;
453 if (print_plus)
454 out << " + ";
455 out << coeffs(indices);
456 for (unsigned int d = 0; d < dim + 1; ++d)
457 {
458 if (indices[d] != 0)
459 out << " * t" << d << '^' << indices[d];
460 }
461 print_plus = true;
462 }
463
464 if (!print_plus)
465 out << Number();
466}
467
468
469
470template <int dim, typename Number>
473{
474 auto deg = coefficients.size();
475 for (unsigned int d = 0; d < dim + 1; ++d)
476 deg[d] -= 1;
477 return deg;
478}
479
480
481
482template <int dim, typename Number>
485{
486 return *this * Number(-1);
487}
488
489
490
491template <int dim, typename Number>
492template <typename Number2>
495{
497 result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
498
499 return result;
500}
501
502
503
504template <int dim, typename Number>
505template <typename Number2>
508{
509 return *this + (-a);
510}
511
512
513
514template <int dim, typename Number>
515template <typename Number2>
518{
519 if (a == Number2())
520 {
522 }
523
525 for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
526 {
527 const auto index = index_to_indices(i, result.coefficients.size());
528 result.coefficients(index) *= a;
529 }
530
531 return result;
532}
533
534
535
536template <int dim, typename Number>
537template <typename Number2>
540{
541 Assert(a != Number2(), ExcDivideByZero());
542 return *this * (Number(1) / Number(a));
543}
544
545
546
547template <int dim, typename Number>
550 const BarycentricPolynomial<dim, Number> &augend) const
551{
553 for (unsigned int d = 0; d < dim + 1; ++d)
554 {
555 deg[d] = std::max(degrees()[d], augend.degrees()[d]);
556 }
557
558 BarycentricPolynomial<dim, Number> result(deg, Number());
559
560 auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
561 for (std::size_t i = 0; i < in.n_elements(); ++i)
562 {
563 const auto index = index_to_indices(i, in.size());
564 result.coefficients(index) += in(index);
565 }
566 };
567
570 return result;
571}
572
573
574
575template <int dim, typename Number>
578 const BarycentricPolynomial<dim, Number> &augend) const
579{
580 return *this + (-augend);
581}
582
583
584
585template <int dim, typename Number>
588 const BarycentricPolynomial<dim, Number> &multiplicand) const
589{
591 for (unsigned int d = 0; d < dim + 1; ++d)
592 {
593 deg[d] = multiplicand.degrees()[d] + degrees()[d];
594 }
595
596 BarycentricPolynomial<dim, Number> result(deg, Number());
597
598 const auto &coef_1 = this->coefficients;
599 const auto &coef_2 = multiplicand.coefficients;
600 auto &coef_out = result.coefficients;
601
602 for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
603 {
604 const auto index_1 = index_to_indices(i1, coef_1.size());
605 for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
606 {
607 const auto index_2 = index_to_indices(i2, coef_2.size());
608
609 TableIndices<dim + 1> index_out;
610 for (unsigned int d = 0; d < dim + 1; ++d)
611 index_out[d] = index_1[d] + index_2[d];
612 coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
613 }
614 }
615
616 return result;
617}
618
619
620
621template <int dim, typename Number>
624 const unsigned int coordinate) const
625{
626 AssertIndexRange(coordinate, dim + 1);
627
628 if (degrees()[coordinate] == 0)
630
631 auto deg = degrees();
632 deg[coordinate] -= 1;
634 std::numeric_limits<Number>::max());
635 const auto &coeffs_in = coefficients;
636 auto &coeffs_out = result.coefficients;
637 for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
638 {
639 const auto out_index = index_to_indices(i, coeffs_out.size());
640 auto input_index = out_index;
641 input_index[coordinate] += 1;
642
643 coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
644 }
645
646 return result;
647}
648
649
650
651template <int dim, typename Number>
654 const unsigned int coordinate) const
655{
656 AssertIndexRange(coordinate, dim);
657 return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
658}
659
660
661
662template <int dim, typename Number>
663Number
665{
666 // TODO: this is probably not numerically stable for higher order.
667 // We really need some version of Horner's method.
668 Number result = {};
669
670 // Begin by converting point (which is in Cartesian coordinates) to
671 // barycentric coordinates:
672 std::array<Number, dim + 1> b_point;
673 b_point[0] = 1.0;
674 for (unsigned int d = 0; d < dim; ++d)
675 {
676 b_point[0] -= point[d];
677 b_point[d + 1] = point[d];
678 }
679
680 // Now evaluate the polynomial at the computed barycentric point:
681 for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
682 {
683 const auto indices = index_to_indices(i, coefficients.size());
684 const auto coef = coefficients(indices);
685 if (coef == Number())
686 continue;
687
688 auto temp = Number(1);
689 for (unsigned int d = 0; d < dim + 1; ++d)
690 temp *= Utilities::pow(b_point[d], indices[d]);
691 result += coef * temp;
692 }
693
694 return result;
695}
696
697
698
699template <int dim, typename Number>
700std::size_t
702{
703 return coefficients.memory_consumption();
704}
705
706
707
708template <int dim, typename Number>
711 const std::size_t &index,
712 const TableIndices<dim + 1> &extents)
713{
715 auto temp = index;
716
717 for (unsigned int n = 0; n < dim + 1; ++n)
718 {
719 std::size_t slice_size = 1;
720 for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
721 slice_size *= extents[n2];
722 result[n] = temp / slice_size;
723 temp %= slice_size;
724 }
725 return result;
726}
727
728
729
730template <int dim>
733{
734 AssertIndexRange(i, polys.size());
735 return polys[i];
736}
737
739
740#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
static TableIndices< dim+1 > index_to_indices(const std::size_t &index, const TableIndices< dim+1 > &extents)
BarycentricPolynomial< dim, Number > operator/(const Number2 &a) const
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:111
virtual unsigned int degree() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > first
Definition grid_out.cc:4623
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
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)