Reference documentation for deal.II version GIT c415589cf0 2022-08-14 18:50:02+00:00
\(\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\}}\)
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 
80 template <int dim, typename Number = double>
82 {
83 public:
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 
154  operator+(const BarycentricPolynomial<dim, Number> &augend) const;
155 
160  operator-(const BarycentricPolynomial<dim, Number> &augend) const;
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 
192 protected:
197 
206  static TableIndices<dim + 1>
207  index_to_indices(const std::size_t & index,
208  const TableIndices<dim + 1> &extent);
209 };
210 
214 template <int dim>
216 {
217 public:
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 
341 protected:
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 
354 template <int dim, typename Number1, typename Number2>
356 operator*(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
357 {
358  return bp * Number1(a);
359 }
360 
364 template <int dim, typename Number1, typename Number2>
366 operator+(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
367 {
368  return bp + Number1(a);
369 }
370 
374 template <int dim, typename Number1, typename Number2>
376 operator-(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
377 {
378  return bp - Number1(a);
379 }
380 
384 template <int dim, typename Number>
385 std::ostream &
386 operator<<(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:
395 template <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 
408 template <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 
423 template <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 
435 template <int dim, typename Number>
436 void
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 
469 template <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 
481 template <int dim, typename Number>
484 {
485  return *this * Number(-1);
486 }
487 
488 
489 
490 template <int dim, typename Number>
491 template <typename Number2>
494 {
496  result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
497 
498  return result;
499 }
500 
501 
502 
503 template <int dim, typename Number>
504 template <typename Number2>
507 {
508  return *this + (-a);
509 }
510 
511 
512 
513 template <int dim, typename Number>
514 template <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 
535 template <int dim, typename Number>
536 template <typename Number2>
539 {
540  Assert(a != Number2(), ExcDivideByZero());
541  return *this * (Number(1) / Number(a));
542 }
543 
544 
545 
546 template <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 
574 template <int dim, typename Number>
577  const BarycentricPolynomial<dim, Number> &augend) const
578 {
579  return *this + (-augend);
580 }
581 
582 
583 
584 template <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 
620 template <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;
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 
650 template <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 
661 template <int dim, typename Number>
662 Number
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 
696 template <int dim, typename Number>
697 std::size_t
699 {
700  return coefficients.memory_consumption();
701 }
702 
703 template <int dim, typename Number>
706  const std::size_t & index,
707  const TableIndices<dim + 1> &extent)
708 {
709  TableIndices<dim + 1> result;
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 
723 template <int dim>
725 BarycentricPolynomials<dim>::operator[](const std::size_t i) const
726 {
727  AssertIndexRange(i, polys.size());
728  return polys[i];
729 }
730 
732 
733 #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
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
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
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
BarycentricPolynomials(const std::vector< BarycentricPolynomial< dim >> &polynomials)
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
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 2 > first
Definition: grid_out.cc:4604
static ::ExceptionBase & ExcDivideByZero()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)