Reference documentation for deal.II version Git d3aed38b93 2021-10-28 13:33:27 +0200
\(\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 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:
221  static const unsigned int dimension = dim;
222 
227  get_fe_p_basis(const unsigned int degree);
228 
233  const std::vector<BarycentricPolynomial<dim>> &polynomials);
234 
239  operator[](const std::size_t i) const;
240 
244  void
245  evaluate(const Point<dim> & unit_point,
246  std::vector<double> & values,
247  std::vector<Tensor<1, dim>> &grads,
248  std::vector<Tensor<2, dim>> &grad_grads,
249  std::vector<Tensor<3, dim>> &third_derivatives,
250  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
251 
255  double
256  compute_value(const unsigned int i, const Point<dim> &p) const override;
257 
262  compute_1st_derivative(const unsigned int i,
263  const Point<dim> & p) const override;
264 
269  compute_2nd_derivative(const unsigned int i,
270  const Point<dim> & p) const override;
271 
276  compute_3rd_derivative(const unsigned int i,
277  const Point<dim> & p) const override;
278 
283  compute_4th_derivative(const unsigned int i,
284  const Point<dim> & p) const override;
285 
290  compute_grad(const unsigned int i, const Point<dim> &p) const override;
291 
296  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
297 
301  virtual std::size_t
302  memory_consumption() const override;
303 
307  std::string
308  name() const override;
309 
313  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
314  clone() const override;
315 
316 protected:
317  std::vector<BarycentricPolynomial<dim>> polys;
318 
320 
322 
324 
326 };
327 
328 // non-member template functions for algebra
329 
333 template <int dim, typename Number1, typename Number2>
335 operator*(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
336 {
337  return bp * Number1(a);
338 }
339 
343 template <int dim, typename Number1, typename Number2>
345 operator+(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
346 {
347  return bp + Number1(a);
348 }
349 
353 template <int dim, typename Number1, typename Number2>
355 operator-(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
356 {
357  return bp - Number1(a);
358 }
359 
363 template <int dim, typename Number>
364 std::ostream &
365 operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
366 {
367  bp.print(out);
368  return out;
369 }
370 
371 // Template function definitions
372 
373 // BarycentricPolynomial:
374 template <int dim, typename Number>
376 {
377  TableIndices<dim + 1> extents;
378  for (unsigned int d = 0; d < dim + 1; ++d)
379  extents[d] = 1;
380  coefficients.reinit(extents);
381 
382  coefficients(TableIndices<dim + 1>{}) = Number();
383 }
384 
385 
386 
387 template <int dim, typename Number>
389  const TableIndices<dim + 1> &powers,
390  const Number coefficient)
391 {
392  TableIndices<dim + 1> extents;
393  for (unsigned int d = 0; d < dim + 1; ++d)
394  extents[d] = powers[d] + 1;
395  coefficients.reinit(extents);
396 
397  coefficients(powers) = coefficient;
398 }
399 
400 
401 
402 template <int dim, typename Number>
405 {
406  AssertIndexRange(d, dim + 1);
407  TableIndices<dim + 1> indices;
408  indices[d] = 1;
409  return BarycentricPolynomial<dim, Number>(indices, Number(1));
410 }
411 
412 
413 
414 template <int dim, typename Number>
415 void
417 {
418  const auto &coeffs = this->coefficients;
419  auto first = index_to_indices(0, coeffs.size());
420  bool print_plus = false;
421  if (coeffs(first) != Number())
422  {
423  out << coeffs(first);
424  print_plus = true;
425  }
426  for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
427  {
428  const auto indices = index_to_indices(i, coeffs.size());
429  if (coeffs(indices) == Number())
430  continue;
431  if (print_plus)
432  out << " + ";
433  out << coeffs(indices);
434  for (unsigned int d = 0; d < dim + 1; ++d)
435  {
436  if (indices[d] != 0)
437  out << " * t" << d << '^' << indices[d];
438  }
439  print_plus = true;
440  }
441 
442  if (!print_plus)
443  out << Number();
444 }
445 
446 
447 
448 template <int dim, typename Number>
451 {
452  auto deg = coefficients.size();
453  for (unsigned int d = 0; d < dim + 1; ++d)
454  deg[d] -= 1;
455  return deg;
456 }
457 
458 
459 
460 template <int dim, typename Number>
463 {
464  return *this * Number(-1);
465 }
466 
467 
468 
469 template <int dim, typename Number>
470 template <typename Number2>
473 {
475  result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
476 
477  return result;
478 }
479 
480 
481 
482 template <int dim, typename Number>
483 template <typename Number2>
486 {
487  return *this + (-a);
488 }
489 
490 
491 
492 template <int dim, typename Number>
493 template <typename Number2>
496 {
497  if (a == Number2())
498  {
500  }
501 
503  for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
504  {
505  const auto index = index_to_indices(i, result.coefficients.size());
506  result.coefficients(index) *= a;
507  }
508 
509  return result;
510 }
511 
512 
513 
514 template <int dim, typename Number>
515 template <typename Number2>
518 {
519  Assert(a != Number2(), ExcDivideByZero());
520  return *this * (Number(1) / Number(a));
521 }
522 
523 
524 
525 template <int dim, typename Number>
528  const BarycentricPolynomial<dim, Number> &augend) const
529 {
531  for (unsigned int d = 0; d < dim + 1; ++d)
532  {
533  deg[d] = std::max(degrees()[d], augend.degrees()[d]);
534  }
535 
536  BarycentricPolynomial<dim, Number> result(deg, Number());
537 
538  auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
539  for (std::size_t i = 0; i < in.n_elements(); ++i)
540  {
541  const auto index = index_to_indices(i, in.size());
542  result.coefficients(index) += in(index);
543  }
544  };
545 
546  add_coefficients(this->coefficients);
547  add_coefficients(augend.coefficients);
548  return result;
549 }
550 
551 
552 
553 template <int dim, typename Number>
556  const BarycentricPolynomial<dim, Number> &augend) const
557 {
558  return *this + (-augend);
559 }
560 
561 
562 
563 template <int dim, typename Number>
566  const BarycentricPolynomial<dim, Number> &multiplicand) const
567 {
569  for (unsigned int d = 0; d < dim + 1; ++d)
570  {
571  deg[d] = multiplicand.degrees()[d] + degrees()[d];
572  }
573 
574  BarycentricPolynomial<dim, Number> result(deg, Number());
575 
576  const auto &coef_1 = this->coefficients;
577  const auto &coef_2 = multiplicand.coefficients;
578  auto & coef_out = result.coefficients;
579 
580  for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
581  {
582  const auto index_1 = index_to_indices(i1, coef_1.size());
583  for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
584  {
585  const auto index_2 = index_to_indices(i2, coef_2.size());
586 
587  TableIndices<dim + 1> index_out;
588  for (unsigned int d = 0; d < dim + 1; ++d)
589  index_out[d] = index_1[d] + index_2[d];
590  coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
591  }
592  }
593 
594  return result;
595 }
596 
597 
598 
599 template <int dim, typename Number>
602  const unsigned int coordinate) const
603 {
604  AssertIndexRange(coordinate, dim + 1);
605 
606  if (degrees()[coordinate] == 0)
608 
609  auto deg = degrees();
610  deg[coordinate] -= 1;
613  const auto & coeffs_in = coefficients;
614  auto & coeffs_out = result.coefficients;
615  for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
616  {
617  const auto out_index = index_to_indices(i, coeffs_out.size());
618  auto input_index = out_index;
619  input_index[coordinate] += 1;
620 
621  coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
622  }
623 
624  return result;
625 }
626 
627 
628 
629 template <int dim, typename Number>
632  const unsigned int coordinate) const
633 {
634  AssertIndexRange(coordinate, dim);
635  return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
636 }
637 
638 
639 
640 template <int dim, typename Number>
641 Number
643 {
644  // TODO: this is probably not numerically stable for higher order.
645  // We really need some version of Horner's method.
646  Number result = {};
647 
648  // Begin by converting point (which is in Cartesian coordinates) to
649  // barycentric coordinates:
650  std::array<Number, dim + 1> b_point;
651  b_point[0] = 1.0;
652  for (unsigned int d = 0; d < dim; ++d)
653  {
654  b_point[0] -= point[d];
655  b_point[d + 1] = point[d];
656  }
657 
658  // Now evaluate the polynomial at the computed barycentric point:
659  for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
660  {
661  const auto indices = index_to_indices(i, coefficients.size());
662  const auto coef = coefficients(indices);
663  if (coef == Number())
664  continue;
665 
666  auto temp = Number(1);
667  for (unsigned int d = 0; d < dim + 1; ++d)
668  temp *= std::pow(b_point[d], indices[d]);
669  result += coef * temp;
670  }
671 
672  return result;
673 }
674 
675 template <int dim, typename Number>
676 std::size_t
678 {
679  return coefficients.memory_consumption();
680 }
681 
682 template <int dim, typename Number>
685  const std::size_t & index,
686  const TableIndices<dim + 1> &extent)
687 {
688  TableIndices<dim + 1> result;
689  auto temp = index;
690 
691  for (unsigned int n = 0; n < dim + 1; ++n)
692  {
693  std::size_t slice_size = 1;
694  for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
695  slice_size *= extent[n2];
696  result[n] = temp / slice_size;
697  temp %= slice_size;
698  }
699  return result;
700 }
701 
702 template <int dim>
704 BarycentricPolynomials<dim>::operator[](const std::size_t i) const
705 {
706  AssertIndexRange(i, polys.size());
707  return polys[i];
708 }
709 
711 
712 #endif
static TableIndices< dim+1 > index_to_indices(const std::size_t &index, const TableIndices< dim+1 > &extent)
Table< 2, BarycentricPolynomial< dim > > poly_grads
Table< 5, BarycentricPolynomial< dim > > poly_fourth_derivatives
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
BarycentricPolynomial< dim, Number > operator*(const Number2 &a) const
Table< 4, BarycentricPolynomial< dim > > poly_third_derivatives
std::size_t memory_consumption() const
BarycentricPolynomial< dim, Number > operator/(const Number2 &a) const
static ::ExceptionBase & ExcDivideByZero()
BarycentricPolynomial< dim, Number > barycentric_derivative(const unsigned int coordinate) const
#define Assert(cond, exc)
Definition: exceptions.h:1461
BarycentricPolynomial< dim, Number > operator+(const Number2 &a) const
BarycentricPolynomial< dim, Number > derivative(const unsigned int coordinate) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
const BarycentricPolynomial< dim > & operator[](const std::size_t i) const
void print(std::ostream &out) const
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:185
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number value(const Point< dim > &point) const
Point< 2 > first
Definition: grid_out.cc:4586
TableIndices< dim+1 > degrees() const
Table< dim+1, Number > coefficients
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
Table< 3, BarycentricPolynomial< dim > > poly_hessians
BarycentricPolynomial< dim, Number > operator-() const
std::vector< BarycentricPolynomial< dim > > polys