Reference documentation for deal.II version GIT a3f17f8a20 2023-06-02 10: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 
26 #include <limits>
27 
29 
82 template <int dim, typename Number = double>
84 {
85 public:
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 
156  operator+(const BarycentricPolynomial<dim, Number> &augend) const;
157 
162  operator-(const BarycentricPolynomial<dim, Number> &augend) const;
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 
194 protected:
199 
208  static TableIndices<dim + 1>
209  index_to_indices(const std::size_t & index,
210  const TableIndices<dim + 1> &extent);
211 };
212 
216 template <int dim>
218 {
219 public:
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,
274  std::vector<Tensor<1, dim>> &grads,
275  std::vector<Tensor<2, dim>> &grad_grads,
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 
323  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
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 
343 protected:
344  std::vector<PolyType> polys;
345  std::vector<GradType> poly_grads;
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 
356 template <int dim, typename Number1, typename Number2>
358 operator*(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
359 {
360  return bp * Number1(a);
361 }
362 
366 template <int dim, typename Number1, typename Number2>
368 operator+(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
369 {
370  return bp + Number1(a);
371 }
372 
376 template <int dim, typename Number1, typename Number2>
378 operator-(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
379 {
380  return bp - Number1(a);
381 }
382 
386 template <int dim, typename Number>
387 std::ostream &
388 operator<<(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:
397 template <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 
410 template <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 
425 template <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 
437 template <int dim, typename Number>
438 void
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 
471 template <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 
483 template <int dim, typename Number>
486 {
487  return *this * Number(-1);
488 }
489 
490 
491 
492 template <int dim, typename Number>
493 template <typename Number2>
496 {
498  result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
499 
500  return result;
501 }
502 
503 
504 
505 template <int dim, typename Number>
506 template <typename Number2>
509 {
510  return *this + (-a);
511 }
512 
513 
514 
515 template <int dim, typename Number>
516 template <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 
537 template <int dim, typename Number>
538 template <typename Number2>
541 {
542  Assert(a != Number2(), ExcDivideByZero());
543  return *this * (Number(1) / Number(a));
544 }
545 
546 
547 
548 template <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 
569  add_coefficients(this->coefficients);
570  add_coefficients(augend.coefficients);
571  return result;
572 }
573 
574 
575 
576 template <int dim, typename Number>
579  const BarycentricPolynomial<dim, Number> &augend) const
580 {
581  return *this + (-augend);
582 }
583 
584 
585 
586 template <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 
622 template <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;
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 
652 template <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 
663 template <int dim, typename Number>
664 Number
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 
698 template <int dim, typename Number>
699 std::size_t
701 {
702  return coefficients.memory_consumption();
703 }
704 
705 template <int dim, typename Number>
708  const std::size_t & index,
709  const TableIndices<dim + 1> &extent)
710 {
711  TableIndices<dim + 1> result;
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 
725 template <int dim>
727 BarycentricPolynomials<dim>::operator[](const std::size_t i) const
728 {
729  AssertIndexRange(i, polys.size());
730  return polys[i];
731 }
732 
734 
735 #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:112
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
virtual unsigned int degree() const
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr DEAL_II_HOST 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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 2 > first
Definition: grid_out.cc:4615
static ::ExceptionBase & ExcDivideByZero()
#define Assert(cond, exc)
Definition: exceptions.h:1614
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:189
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)