Reference documentation for deal.II version 9.3.0
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>
140  BarycentricPolynomial<dim, Number> operator*(const Number2 &a) const;
141
145  template <typename Number2>
147  operator/(const Number2 &a) const;
148
153  operator+(const BarycentricPolynomial<dim, Number> &augend) const;
154
159  operator-(const BarycentricPolynomial<dim, Number> &augend) const;
160
165  operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
166
171  barycentric_derivative(const unsigned int coordinate) const;
172
177  derivative(const unsigned int coordinate) const;
178
182  Number
183  value(const Point<dim> &point) const;
184
188  std::size_t
189  memory_consumption() const;
190
191 protected:
196
205  static TableIndices<dim + 1>
206  index_to_indices(const std::size_t & index,
207  const TableIndices<dim + 1> &extent);
208 };
209
213 template <int dim>
215 {
216 public:
220  static const unsigned int dimension = dim;
221
226  get_fe_p_basis(const unsigned int degree);
227
232  const std::vector<BarycentricPolynomial<dim>> &polynomials);
233
237  const BarycentricPolynomial<dim> &operator[](const std::size_t i) const;
238
242  void
243  evaluate(const Point<dim> & unit_point,
244  std::vector<double> & values,
245  std::vector<Tensor<1, dim>> &grads,
246  std::vector<Tensor<2, dim>> &grad_grads,
247  std::vector<Tensor<3, dim>> &third_derivatives,
248  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
249
253  double
254  compute_value(const unsigned int i, const Point<dim> &p) const override;
255
260  compute_1st_derivative(const unsigned int i,
261  const Point<dim> & p) const override;
262
267  compute_2nd_derivative(const unsigned int i,
268  const Point<dim> & p) const override;
269
274  compute_3rd_derivative(const unsigned int i,
275  const Point<dim> & p) const override;
276
281  compute_4th_derivative(const unsigned int i,
282  const Point<dim> & p) const override;
283
288  compute_grad(const unsigned int i, const Point<dim> &p) const override;
289
294  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
295
299  virtual std::size_t
300  memory_consumption() const override;
301
305  std::string
306  name() const override;
307
311  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
312  clone() const override;
313
314 protected:
315  std::vector<BarycentricPolynomial<dim>> polys;
316
318
320
322
324 };
325
326 // non-member template functions for algebra
327
331 template <int dim, typename Number1, typename Number2>
333 operator*(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
334 {
335  return bp * Number1(a);
336 }
337
341 template <int dim, typename Number1, typename Number2>
343 operator+(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
344 {
345  return bp + Number1(a);
346 }
347
351 template <int dim, typename Number1, typename Number2>
353 operator-(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
354 {
355  return bp - Number1(a);
356 }
357
361 template <int dim, typename Number>
362 std::ostream &
363 operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
364 {
365  bp.print(out);
366  return out;
367 }
368
369 // Template function definitions
370
371 // BarycentricPolynomial:
372 template <int dim, typename Number>
374 {
375  TableIndices<dim + 1> extents;
376  for (unsigned int d = 0; d < dim + 1; ++d)
377  extents[d] = 1;
378  coefficients.reinit(extents);
379
380  coefficients(TableIndices<dim + 1>{}) = Number();
381 }
382
383
384
385 template <int dim, typename Number>
387  const TableIndices<dim + 1> &powers,
388  const Number coefficient)
389 {
390  TableIndices<dim + 1> extents;
391  for (unsigned int d = 0; d < dim + 1; ++d)
392  extents[d] = powers[d] + 1;
393  coefficients.reinit(extents);
394
395  coefficients(powers) = coefficient;
396 }
397
398
399
400 template <int dim, typename Number>
403 {
404  AssertIndexRange(d, dim + 1);
405  TableIndices<dim + 1> indices;
406  indices[d] = 1;
407  return BarycentricPolynomial<dim, Number>(indices, Number(1));
408 }
409
410
411
412 template <int dim, typename Number>
413 void
415 {
416  const auto &coeffs = this->coefficients;
417  auto first = index_to_indices(0, coeffs.size());
418  bool print_plus = false;
419  if (coeffs(first) != Number())
420  {
421  out << coeffs(first);
422  print_plus = true;
423  }
424  for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
425  {
426  const auto indices = index_to_indices(i, coeffs.size());
427  if (coeffs(indices) == Number())
428  continue;
429  if (print_plus)
430  out << " + ";
431  out << coeffs(indices);
432  for (unsigned int d = 0; d < dim + 1; ++d)
433  {
434  if (indices[d] != 0)
435  out << " * t" << d << '^' << indices[d];
436  }
437  print_plus = true;
438  }
439
440  if (!print_plus)
441  out << Number();
442 }
443
444
445
446 template <int dim, typename Number>
449 {
450  auto deg = coefficients.size();
451  for (unsigned int d = 0; d < dim + 1; ++d)
452  deg[d] -= 1;
453  return deg;
454 }
455
456
457
458 template <int dim, typename Number>
461 {
462  return *this * Number(-1);
463 }
464
465
466
467 template <int dim, typename Number>
468 template <typename Number2>
471 {
473  result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
474
475  return result;
476 }
477
478
479
480 template <int dim, typename Number>
481 template <typename Number2>
484 {
485  return *this + (-a);
486 }
487
488
489
490 template <int dim, typename Number>
491 template <typename Number2>
493  operator*(const Number2 &a) const
494 {
495  if (a == Number2())
496  {
498  }
499
501  for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
502  {
503  const auto index = index_to_indices(i, result.coefficients.size());
504  result.coefficients(index) *= a;
505  }
506
507  return result;
508 }
509
510
511
512 template <int dim, typename Number>
513 template <typename Number2>
516 {
517  Assert(a != Number2(), ExcDivideByZero());
518  return *this * (Number(1) / Number(a));
519 }
520
521
522
523 template <int dim, typename Number>
527 {
529  for (unsigned int d = 0; d < dim + 1; ++d)
530  {
531  deg[d] = std::max(degrees()[d], augend.degrees()[d]);
532  }
533
534  BarycentricPolynomial<dim, Number> result(deg, Number());
535
536  auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
537  for (std::size_t i = 0; i < in.n_elements(); ++i)
538  {
539  const auto index = index_to_indices(i, in.size());
540  result.coefficients(index) += in(index);
541  }
542  };
543
544  add_coefficients(this->coefficients);
545  add_coefficients(augend.coefficients);
546  return result;
547 }
548
549
550
551 template <int dim, typename Number>
555 {
556  return *this + (-augend);
557 }
558
559
560
561 template <int dim, typename Number>
564 {
566  for (unsigned int d = 0; d < dim + 1; ++d)
567  {
568  deg[d] = multiplicand.degrees()[d] + degrees()[d];
569  }
570
571  BarycentricPolynomial<dim, Number> result(deg, Number());
572
573  const auto &coef_1 = this->coefficients;
574  const auto &coef_2 = multiplicand.coefficients;
575  auto & coef_out = result.coefficients;
576
577  for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
578  {
579  const auto index_1 = index_to_indices(i1, coef_1.size());
580  for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
581  {
582  const auto index_2 = index_to_indices(i2, coef_2.size());
583
584  TableIndices<dim + 1> index_out;
585  for (unsigned int d = 0; d < dim + 1; ++d)
586  index_out[d] = index_1[d] + index_2[d];
587  coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
588  }
589  }
590
591  return result;
592 }
593
594
595
596 template <int dim, typename Number>
599  const unsigned int coordinate) const
600 {
601  AssertIndexRange(coordinate, dim + 1);
602
603  if (degrees()[coordinate] == 0)
605
606  auto deg = degrees();
607  deg[coordinate] -= 1;
610  const auto & coeffs_in = coefficients;
611  auto & coeffs_out = result.coefficients;
612  for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
613  {
614  const auto out_index = index_to_indices(i, coeffs_out.size());
615  auto input_index = out_index;
616  input_index[coordinate] += 1;
617
618  coeffs_out(out_index) = coeffs_in(input_index) * input_index[coordinate];
619  }
620
621  return result;
622 }
623
624
625
626 template <int dim, typename Number>
629  const unsigned int coordinate) const
630 {
631  AssertIndexRange(coordinate, dim);
632  return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
633 }
634
635
636
637 template <int dim, typename Number>
638 Number
640 {
641  // TODO: this is probably not numerically stable for higher order.
642  // We really need some version of Horner's method.
643  Number result = {};
644
645  // Begin by converting point (which is in Cartesian coordinates) to
646  // barycentric coordinates:
647  std::array<Number, dim + 1> b_point;
648  b_point[0] = 1.0;
649  for (unsigned int d = 0; d < dim; ++d)
650  {
651  b_point[0] -= point[d];
652  b_point[d + 1] = point[d];
653  }
654
655  // Now evaluate the polynomial at the computed barycentric point:
656  for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
657  {
658  const auto indices = index_to_indices(i, coefficients.size());
659  const auto coef = coefficients(indices);
660  if (coef == Number())
661  continue;
662
663  auto temp = Number(1);
664  for (unsigned int d = 0; d < dim + 1; ++d)
665  temp *= std::pow(b_point[d], indices[d]);
666  result += coef * temp;
667  }
668
669  return result;
670 }
671
672 template <int dim, typename Number>
673 std::size_t
675 {
676  return coefficients.memory_consumption();
677 }
678
679 template <int dim, typename Number>
682  const std::size_t & index,
683  const TableIndices<dim + 1> &extent)
684 {
685  TableIndices<dim + 1> result;
686  auto temp = index;
687
688  for (unsigned int n = 0; n < dim + 1; ++n)
689  {
690  std::size_t slice_size = 1;
691  for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
692  slice_size *= extent[n2];
693  result[n] = temp / slice_size;
694  temp %= slice_size;
695  }
696  return result;
697 }
698
699 template <int dim>
701  operator[](const std::size_t i) const
702 {
703  AssertIndexRange(i, polys.size());
704  return polys[i];
705 }
706
708
709 #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:1690
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:1465
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:395
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:188
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:4587
TableIndices< dim+1 > degrees() const
Table< dim+1, Number > coefficients
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
Table< 3, BarycentricPolynomial< dim > > poly_hessians
BarycentricPolynomial< dim, Number > operator-() const
std::vector< BarycentricPolynomial< dim > > polys