Reference documentation for deal.II version Git 689de043d4 2020-08-10 16:46:15 +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\}}\)
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 
70 template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
72 {
73 public:
78  static const unsigned int dimension = dim;
79 
86  template <class Pol>
87  TensorProductPolynomials(const std::vector<Pol> &pols);
88 
92  void
93  output_indices(std::ostream &out) const;
94 
99  void
100  set_numbering(const std::vector<unsigned int> &renumber);
101 
105  const std::vector<unsigned int> &
106  get_numbering() const;
107 
111  const std::vector<unsigned int> &
112  get_numbering_inverse() const;
113 
126  void
127  evaluate(const Point<dim> & unit_point,
128  std::vector<double> & values,
129  std::vector<Tensor<1, dim>> &grads,
130  std::vector<Tensor<2, dim>> &grad_grads,
131  std::vector<Tensor<3, dim>> &third_derivatives,
132  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
133 
146  double
147  compute_value(const unsigned int i, const Point<dim> &p) const override;
148 
163  template <int order>
165  compute_derivative(const unsigned int i, const Point<dim> &p) const;
166 
170  virtual Tensor<1, dim>
171  compute_1st_derivative(const unsigned int i,
172  const Point<dim> & p) const override;
173 
177  virtual Tensor<2, dim>
178  compute_2nd_derivative(const unsigned int i,
179  const Point<dim> & p) const override;
180 
184  virtual Tensor<3, dim>
185  compute_3rd_derivative(const unsigned int i,
186  const Point<dim> & p) const override;
187 
191  virtual Tensor<4, dim>
192  compute_4th_derivative(const unsigned int i,
193  const Point<dim> & p) const override;
194 
208  compute_grad(const unsigned int i, const Point<dim> &p) const override;
209 
223  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
224 
228  std::string
229  name() const override;
230 
234  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
235  clone() const override;
236 
240  virtual std::size_t
241  memory_consumption() const override;
242 
243 protected:
247  std::vector<PolynomialType> polynomials;
248 
252  std::vector<unsigned int> index_map;
253 
257  std::vector<unsigned int> index_map_inverse;
258 
265  void
266  compute_index(const unsigned int i,
267  std::array<unsigned int, dim> &indices) const;
268 
274 
279  friend class TensorProductPolynomialsConst<dim>;
280 };
281 
282 
283 
309 template <int dim>
311 {
312 public:
329  const std::vector<std::vector<Polynomials::Polynomial<double>>>
330  &base_polynomials);
331 
345  void
346  evaluate(const Point<dim> & unit_point,
347  std::vector<double> & values,
348  std::vector<Tensor<1, dim>> &grads,
349  std::vector<Tensor<2, dim>> &grad_grads,
350  std::vector<Tensor<3, dim>> &third_derivatives,
351  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
352 
365  double
366  compute_value(const unsigned int i, const Point<dim> &p) const override;
367 
382  template <int order>
384  compute_derivative(const unsigned int i, const Point<dim> &p) const;
385 
389  virtual Tensor<1, dim>
390  compute_1st_derivative(const unsigned int i,
391  const Point<dim> & p) const override;
392 
396  virtual Tensor<2, dim>
397  compute_2nd_derivative(const unsigned int i,
398  const Point<dim> & p) const override;
399 
403  virtual Tensor<3, dim>
404  compute_3rd_derivative(const unsigned int i,
405  const Point<dim> & p) const override;
406 
410  virtual Tensor<4, dim>
411  compute_4th_derivative(const unsigned int i,
412  const Point<dim> & p) const override;
413 
427  compute_grad(const unsigned int i, const Point<dim> &p) const override;
428 
442  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
443 
447  std::string
448  name() const override;
449 
453  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
454  clone() const override;
455 
456 private:
460  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
461 
468  void
469  compute_index(const unsigned int i,
470  std::array<unsigned int, dim> &indices) const;
471 
475  static unsigned int
476  get_n_tensor_pols(
477  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
478 };
479 
482 #ifndef DOXYGEN
483 
484 
485 /* ---------------- template and inline functions ---------- */
486 
487 
488 template <int dim, typename PolynomialType>
489 template <class Pol>
491  const std::vector<Pol> &pols)
492  : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
493  , polynomials(pols.begin(), pols.end())
494  , index_map(this->n())
495  , index_map_inverse(this->n())
496 {
497  // per default set this index map to identity. This map can be changed by
498  // the user through the set_numbering() function
499  for (unsigned int i = 0; i < this->n(); ++i)
500  {
501  index_map[i] = i;
502  index_map_inverse[i] = i;
503  }
504 }
505 
506 
507 template <int dim, typename PolynomialType>
508 inline const std::vector<unsigned int> &
510 {
511  return index_map;
512 }
513 
514 
515 template <int dim, typename PolynomialType>
516 inline const std::vector<unsigned int> &
518 {
519  return index_map_inverse;
520 }
521 
522 
523 template <int dim, typename PolynomialType>
524 inline std::string
526 {
527  return "TensorProductPolynomials";
528 }
529 
530 
531 template <int dim, typename PolynomialType>
532 template <int order>
535  const unsigned int i,
536  const Point<dim> & p) const
537 {
538  std::array<unsigned int, dim> indices;
539  compute_index(i, indices);
540 
541  std::array<std::array<double, 5>, dim> v;
542  {
543  std::vector<double> tmp(5);
544  for (unsigned int d = 0; d < dim; ++d)
545  {
546  polynomials[indices[d]].value(p(d), tmp);
547  v[d][0] = tmp[0];
548  v[d][1] = tmp[1];
549  v[d][2] = tmp[2];
550  v[d][3] = tmp[3];
551  v[d][4] = tmp[4];
552  }
553  }
554 
555  Tensor<order, dim> derivative;
556  switch (order)
557  {
558  case 1:
559  {
560  Tensor<1, dim> &derivative_1 =
561  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
562  for (unsigned int d = 0; d < dim; ++d)
563  {
564  derivative_1[d] = 1.;
565  for (unsigned int x = 0; x < dim; ++x)
566  {
567  unsigned int x_order = 0;
568  if (d == x)
569  ++x_order;
570 
571  derivative_1[d] *= v[x][x_order];
572  }
573  }
574 
575  return derivative;
576  }
577  case 2:
578  {
579  Tensor<2, dim> &derivative_2 =
580  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
581  for (unsigned int d1 = 0; d1 < dim; ++d1)
582  for (unsigned int d2 = 0; d2 < dim; ++d2)
583  {
584  derivative_2[d1][d2] = 1.;
585  for (unsigned int x = 0; x < dim; ++x)
586  {
587  unsigned int x_order = 0;
588  if (d1 == x)
589  ++x_order;
590  if (d2 == x)
591  ++x_order;
592 
593  derivative_2[d1][d2] *= v[x][x_order];
594  }
595  }
596 
597  return derivative;
598  }
599  case 3:
600  {
601  Tensor<3, dim> &derivative_3 =
602  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
603  for (unsigned int d1 = 0; d1 < dim; ++d1)
604  for (unsigned int d2 = 0; d2 < dim; ++d2)
605  for (unsigned int d3 = 0; d3 < dim; ++d3)
606  {
607  derivative_3[d1][d2][d3] = 1.;
608  for (unsigned int x = 0; x < dim; ++x)
609  {
610  unsigned int x_order = 0;
611  if (d1 == x)
612  ++x_order;
613  if (d2 == x)
614  ++x_order;
615  if (d3 == x)
616  ++x_order;
617 
618  derivative_3[d1][d2][d3] *= v[x][x_order];
619  }
620  }
621 
622  return derivative;
623  }
624  case 4:
625  {
626  Tensor<4, dim> &derivative_4 =
627  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
628  for (unsigned int d1 = 0; d1 < dim; ++d1)
629  for (unsigned int d2 = 0; d2 < dim; ++d2)
630  for (unsigned int d3 = 0; d3 < dim; ++d3)
631  for (unsigned int d4 = 0; d4 < dim; ++d4)
632  {
633  derivative_4[d1][d2][d3][d4] = 1.;
634  for (unsigned int x = 0; x < dim; ++x)
635  {
636  unsigned int x_order = 0;
637  if (d1 == x)
638  ++x_order;
639  if (d2 == x)
640  ++x_order;
641  if (d3 == x)
642  ++x_order;
643  if (d4 == x)
644  ++x_order;
645 
646  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
647  }
648  }
649 
650  return derivative;
651  }
652  default:
653  {
654  Assert(false, ExcNotImplemented());
655  return derivative;
656  }
657  }
658 }
659 
660 
661 
662 template <int dim, typename PolynomialType>
663 inline Tensor<1, dim>
665  const unsigned int i,
666  const Point<dim> & p) const
667 {
668  return compute_derivative<1>(i, p);
669 }
670 
671 
672 
673 template <int dim, typename PolynomialType>
674 inline Tensor<2, dim>
676  const unsigned int i,
677  const Point<dim> & p) const
678 {
679  return compute_derivative<2>(i, p);
680 }
681 
682 
683 
684 template <int dim, typename PolynomialType>
685 inline Tensor<3, dim>
687  const unsigned int i,
688  const Point<dim> & p) const
689 {
690  return compute_derivative<3>(i, p);
691 }
692 
693 
694 
695 template <int dim, typename PolynomialType>
696 inline Tensor<4, dim>
698  const unsigned int i,
699  const Point<dim> & p) const
700 {
701  return compute_derivative<4>(i, p);
702 }
703 
704 
705 
706 template <int dim>
707 template <int order>
710  const Point<dim> & p) const
711 {
712  std::array<unsigned int, dim> indices;
713  compute_index(i, indices);
714 
715  std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
716  for (unsigned int d = 0; d < dim; ++d)
717  polynomials[d][indices[d]].value(p(d), v[d]);
718 
719  Tensor<order, dim> derivative;
720  switch (order)
721  {
722  case 1:
723  {
724  Tensor<1, dim> &derivative_1 =
725  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
726  for (unsigned int d = 0; d < dim; ++d)
727  {
728  derivative_1[d] = 1.;
729  for (unsigned int x = 0; x < dim; ++x)
730  {
731  unsigned int x_order = 0;
732  if (d == x)
733  ++x_order;
734 
735  derivative_1[d] *= v[x][x_order];
736  }
737  }
738 
739  return derivative;
740  }
741  case 2:
742  {
743  Tensor<2, dim> &derivative_2 =
744  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
745  for (unsigned int d1 = 0; d1 < dim; ++d1)
746  for (unsigned int d2 = 0; d2 < dim; ++d2)
747  {
748  derivative_2[d1][d2] = 1.;
749  for (unsigned int x = 0; x < dim; ++x)
750  {
751  unsigned int x_order = 0;
752  if (d1 == x)
753  ++x_order;
754  if (d2 == x)
755  ++x_order;
756 
757  derivative_2[d1][d2] *= v[x][x_order];
758  }
759  }
760 
761  return derivative;
762  }
763  case 3:
764  {
765  Tensor<3, dim> &derivative_3 =
766  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
767  for (unsigned int d1 = 0; d1 < dim; ++d1)
768  for (unsigned int d2 = 0; d2 < dim; ++d2)
769  for (unsigned int d3 = 0; d3 < dim; ++d3)
770  {
771  derivative_3[d1][d2][d3] = 1.;
772  for (unsigned int x = 0; x < dim; ++x)
773  {
774  unsigned int x_order = 0;
775  if (d1 == x)
776  ++x_order;
777  if (d2 == x)
778  ++x_order;
779  if (d3 == x)
780  ++x_order;
781 
782  derivative_3[d1][d2][d3] *= v[x][x_order];
783  }
784  }
785 
786  return derivative;
787  }
788  case 4:
789  {
790  Tensor<4, dim> &derivative_4 =
791  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
792  for (unsigned int d1 = 0; d1 < dim; ++d1)
793  for (unsigned int d2 = 0; d2 < dim; ++d2)
794  for (unsigned int d3 = 0; d3 < dim; ++d3)
795  for (unsigned int d4 = 0; d4 < dim; ++d4)
796  {
797  derivative_4[d1][d2][d3][d4] = 1.;
798  for (unsigned int x = 0; x < dim; ++x)
799  {
800  unsigned int x_order = 0;
801  if (d1 == x)
802  ++x_order;
803  if (d2 == x)
804  ++x_order;
805  if (d3 == x)
806  ++x_order;
807  if (d4 == x)
808  ++x_order;
809 
810  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
811  }
812  }
813 
814  return derivative;
815  }
816  default:
817  {
818  Assert(false, ExcNotImplemented());
819  return derivative;
820  }
821  }
822 }
823 
824 
825 
826 template <int dim>
827 inline Tensor<1, dim>
829  const Point<dim> & p) const
830 {
831  return compute_derivative<1>(i, p);
832 }
833 
834 
835 
836 template <int dim>
837 inline Tensor<2, dim>
839  const Point<dim> & p) const
840 {
841  return compute_derivative<2>(i, p);
842 }
843 
844 
845 
846 template <int dim>
847 inline Tensor<3, dim>
849  const Point<dim> & p) const
850 {
851  return compute_derivative<3>(i, p);
852 }
853 
854 
855 
856 template <int dim>
857 inline Tensor<4, dim>
859  const Point<dim> & p) const
860 {
861  return compute_derivative<4>(i, p);
862 }
863 
864 
865 
866 template <int dim>
867 inline std::string
869 {
870  return "AnisotropicPolynomials";
871 }
872 
873 
874 
875 #endif // DOXYGEN
877 
878 #endif
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:1045
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
#define Assert(cond, exc)
Definition: exceptions.h:1411
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:363
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
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
std::string name() const override
double compute_value(const unsigned int i, const Point< dim > &p) const override
Definition: cuda.h:31
std::string name() const override
static const unsigned int dimension
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
const std::vector< unsigned int > & get_numbering_inverse() const
virtual Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
static ::ExceptionBase & ExcNotImplemented()
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
virtual Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override