Reference documentation for deal.II version 9.2.0
\(\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\}}\)
tensor_product_polynomials.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2020 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 #ifndef dealii_tensor_product_polynomials_h
17 #define dealii_tensor_product_polynomials_h
18 
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/point.h>
26 #include <deal.II/base/tensor.h>
27 #include <deal.II/base/utilities.h>
28 
29 #include <vector>
30 
32 
33 // Forward declarations for friends
34 // TODO: We may be able to modify these classes so they aren't
35 // required to be friends
36 template <int dim>
38 template <int dim>
40 
73 template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
75 {
76 public:
81  static const unsigned int dimension = dim;
82 
89  template <class Pol>
90  TensorProductPolynomials(const std::vector<Pol> &pols);
91 
95  void
96  output_indices(std::ostream &out) const;
97 
102  void
103  set_numbering(const std::vector<unsigned int> &renumber);
104 
108  const std::vector<unsigned int> &
109  get_numbering() const;
110 
114  const std::vector<unsigned int> &
115  get_numbering_inverse() const;
116 
129  void
130  evaluate(const Point<dim> & unit_point,
131  std::vector<double> & values,
132  std::vector<Tensor<1, dim>> &grads,
133  std::vector<Tensor<2, dim>> &grad_grads,
134  std::vector<Tensor<3, dim>> &third_derivatives,
135  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
136 
149  double
150  compute_value(const unsigned int i, const Point<dim> &p) const;
151 
166  template <int order>
168  compute_derivative(const unsigned int i, const Point<dim> &p) const;
169 
183  compute_grad(const unsigned int i, const Point<dim> &p) const;
184 
198  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
199 
203  std::string
204  name() const override;
205 
209  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
210  clone() const override;
211 
215  virtual std::size_t
216  memory_consumption() const override;
217 
218 protected:
222  std::vector<PolynomialType> polynomials;
223 
227  std::vector<unsigned int> index_map;
228 
232  std::vector<unsigned int> index_map_inverse;
233 
240  void
241  compute_index(const unsigned int i,
242  std::array<unsigned int, dim> &indices) const;
243 
249 
254  friend class TensorProductPolynomialsConst<dim>;
255 };
256 
257 
258 
286 template <int dim>
288 {
289 public:
306  const std::vector<std::vector<Polynomials::Polynomial<double>>>
307  &base_polynomials);
308 
322  void
323  evaluate(const Point<dim> & unit_point,
324  std::vector<double> & values,
325  std::vector<Tensor<1, dim>> &grads,
326  std::vector<Tensor<2, dim>> &grad_grads,
327  std::vector<Tensor<3, dim>> &third_derivatives,
328  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
329 
342  double
343  compute_value(const unsigned int i, const Point<dim> &p) const;
344 
359  template <int order>
361  compute_derivative(const unsigned int i, const Point<dim> &p) const;
362 
376  compute_grad(const unsigned int i, const Point<dim> &p) const;
377 
391  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
392 
396  std::string
397  name() const override;
398 
402  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
403  clone() const override;
404 
405 private:
409  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
410 
417  void
418  compute_index(const unsigned int i,
419  std::array<unsigned int, dim> &indices) const;
420 
424  static unsigned int
425  get_n_tensor_pols(
426  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
427 };
428 
431 #ifndef DOXYGEN
432 
433 
434 /* ---------------- template and inline functions ---------- */
435 
436 
437 template <int dim, typename PolynomialType>
438 template <class Pol>
440  const std::vector<Pol> &pols)
441  : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
442  , polynomials(pols.begin(), pols.end())
443  , index_map(this->n())
444  , index_map_inverse(this->n())
445 {
446  // per default set this index map to identity. This map can be changed by
447  // the user through the set_numbering() function
448  for (unsigned int i = 0; i < this->n(); ++i)
449  {
450  index_map[i] = i;
451  index_map_inverse[i] = i;
452  }
453 }
454 
455 
456 template <int dim, typename PolynomialType>
457 inline const std::vector<unsigned int> &
459 {
460  return index_map;
461 }
462 
463 
464 template <int dim, typename PolynomialType>
465 inline const std::vector<unsigned int> &
467 {
468  return index_map_inverse;
469 }
470 
471 
472 template <int dim, typename PolynomialType>
473 inline std::string
475 {
476  return "TensorProductPolynomials";
477 }
478 
479 
480 template <int dim, typename PolynomialType>
481 template <int order>
484  const unsigned int i,
485  const Point<dim> & p) const
486 {
487  std::array<unsigned int, dim> indices;
488  compute_index(i, indices);
489 
490  std::array<std::array<double, 5>, dim> v;
491  {
492  std::vector<double> tmp(5);
493  for (unsigned int d = 0; d < dim; ++d)
494  {
495  polynomials[indices[d]].value(p(d), tmp);
496  v[d][0] = tmp[0];
497  v[d][1] = tmp[1];
498  v[d][2] = tmp[2];
499  v[d][3] = tmp[3];
500  v[d][4] = tmp[4];
501  }
502  }
503 
504  Tensor<order, dim> derivative;
505  switch (order)
506  {
507  case 1:
508  {
509  Tensor<1, dim> &derivative_1 =
510  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
511  for (unsigned int d = 0; d < dim; ++d)
512  {
513  derivative_1[d] = 1.;
514  for (unsigned int x = 0; x < dim; ++x)
515  {
516  unsigned int x_order = 0;
517  if (d == x)
518  ++x_order;
519 
520  derivative_1[d] *= v[x][x_order];
521  }
522  }
523 
524  return derivative;
525  }
526  case 2:
527  {
528  Tensor<2, dim> &derivative_2 =
529  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
530  for (unsigned int d1 = 0; d1 < dim; ++d1)
531  for (unsigned int d2 = 0; d2 < dim; ++d2)
532  {
533  derivative_2[d1][d2] = 1.;
534  for (unsigned int x = 0; x < dim; ++x)
535  {
536  unsigned int x_order = 0;
537  if (d1 == x)
538  ++x_order;
539  if (d2 == x)
540  ++x_order;
541 
542  derivative_2[d1][d2] *= v[x][x_order];
543  }
544  }
545 
546  return derivative;
547  }
548  case 3:
549  {
550  Tensor<3, dim> &derivative_3 =
551  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
552  for (unsigned int d1 = 0; d1 < dim; ++d1)
553  for (unsigned int d2 = 0; d2 < dim; ++d2)
554  for (unsigned int d3 = 0; d3 < dim; ++d3)
555  {
556  derivative_3[d1][d2][d3] = 1.;
557  for (unsigned int x = 0; x < dim; ++x)
558  {
559  unsigned int x_order = 0;
560  if (d1 == x)
561  ++x_order;
562  if (d2 == x)
563  ++x_order;
564  if (d3 == x)
565  ++x_order;
566 
567  derivative_3[d1][d2][d3] *= v[x][x_order];
568  }
569  }
570 
571  return derivative;
572  }
573  case 4:
574  {
575  Tensor<4, dim> &derivative_4 =
576  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
577  for (unsigned int d1 = 0; d1 < dim; ++d1)
578  for (unsigned int d2 = 0; d2 < dim; ++d2)
579  for (unsigned int d3 = 0; d3 < dim; ++d3)
580  for (unsigned int d4 = 0; d4 < dim; ++d4)
581  {
582  derivative_4[d1][d2][d3][d4] = 1.;
583  for (unsigned int x = 0; x < dim; ++x)
584  {
585  unsigned int x_order = 0;
586  if (d1 == x)
587  ++x_order;
588  if (d2 == x)
589  ++x_order;
590  if (d3 == x)
591  ++x_order;
592  if (d4 == x)
593  ++x_order;
594 
595  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
596  }
597  }
598 
599  return derivative;
600  }
601  default:
602  {
603  Assert(false, ExcNotImplemented());
604  return derivative;
605  }
606  }
607 }
608 
609 template <int dim>
610 template <int order>
613  const Point<dim> & p) const
614 {
615  std::array<unsigned int, dim> indices;
616  compute_index(i, indices);
617 
618  std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
619  for (unsigned int d = 0; d < dim; ++d)
620  polynomials[d][indices[d]].value(p(d), v[d]);
621 
622  Tensor<order, dim> derivative;
623  switch (order)
624  {
625  case 1:
626  {
627  Tensor<1, dim> &derivative_1 =
628  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
629  for (unsigned int d = 0; d < dim; ++d)
630  {
631  derivative_1[d] = 1.;
632  for (unsigned int x = 0; x < dim; ++x)
633  {
634  unsigned int x_order = 0;
635  if (d == x)
636  ++x_order;
637 
638  derivative_1[d] *= v[x][x_order];
639  }
640  }
641 
642  return derivative;
643  }
644  case 2:
645  {
646  Tensor<2, dim> &derivative_2 =
647  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
648  for (unsigned int d1 = 0; d1 < dim; ++d1)
649  for (unsigned int d2 = 0; d2 < dim; ++d2)
650  {
651  derivative_2[d1][d2] = 1.;
652  for (unsigned int x = 0; x < dim; ++x)
653  {
654  unsigned int x_order = 0;
655  if (d1 == x)
656  ++x_order;
657  if (d2 == x)
658  ++x_order;
659 
660  derivative_2[d1][d2] *= v[x][x_order];
661  }
662  }
663 
664  return derivative;
665  }
666  case 3:
667  {
668  Tensor<3, dim> &derivative_3 =
669  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
670  for (unsigned int d1 = 0; d1 < dim; ++d1)
671  for (unsigned int d2 = 0; d2 < dim; ++d2)
672  for (unsigned int d3 = 0; d3 < dim; ++d3)
673  {
674  derivative_3[d1][d2][d3] = 1.;
675  for (unsigned int x = 0; x < dim; ++x)
676  {
677  unsigned int x_order = 0;
678  if (d1 == x)
679  ++x_order;
680  if (d2 == x)
681  ++x_order;
682  if (d3 == x)
683  ++x_order;
684 
685  derivative_3[d1][d2][d3] *= v[x][x_order];
686  }
687  }
688 
689  return derivative;
690  }
691  case 4:
692  {
693  Tensor<4, dim> &derivative_4 =
694  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
695  for (unsigned int d1 = 0; d1 < dim; ++d1)
696  for (unsigned int d2 = 0; d2 < dim; ++d2)
697  for (unsigned int d3 = 0; d3 < dim; ++d3)
698  for (unsigned int d4 = 0; d4 < dim; ++d4)
699  {
700  derivative_4[d1][d2][d3][d4] = 1.;
701  for (unsigned int x = 0; x < dim; ++x)
702  {
703  unsigned int x_order = 0;
704  if (d1 == x)
705  ++x_order;
706  if (d2 == x)
707  ++x_order;
708  if (d3 == x)
709  ++x_order;
710  if (d4 == x)
711  ++x_order;
712 
713  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
714  }
715  }
716 
717  return derivative;
718  }
719  default:
720  {
721  Assert(false, ExcNotImplemented());
722  return derivative;
723  }
724  }
725 }
726 
727 
728 template <int dim>
729 inline std::string
731 {
732  return "AnisotropicPolynomials";
733 }
734 
735 
736 
737 #endif // DOXYGEN
739 
740 #endif
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
std::vector< PolynomialType > polynomials
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
void output_indices(std::ostream &out) const
void set_numbering(const std::vector< unsigned int > &renumber)
TensorProductPolynomials(const std::vector< Pol > &pols)
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
T fixed_power(const T t)
Definition: utilities.h:1072
#define Assert(cond, exc)
Definition: exceptions.h:1419
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
virtual std::size_t memory_consumption() const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
const std::vector< std::vector< Polynomials::Polynomial< double > > > polynomials
std::vector< unsigned int > index_map_inverse
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
const std::vector< unsigned int > & get_numbering() const
std::string name() const override
Definition: cuda.h:31
std::string name() const override
static const unsigned int dimension
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
std::vector< unsigned int > index_map
const std::vector< unsigned int > & get_numbering_inverse() const
static ::ExceptionBase & ExcNotImplemented()
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
static const bool value