Reference documentation for deal.II version Git 6d63218887 2020-10-30 17:17:53 -0400
\(\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 <>
663 template <int order>
666  compute_derivative(const unsigned int, const Point<0> &) const
667 {
668  AssertThrow(false, ExcNotImplemented());
669 
670  return {};
671 }
672 
673 
674 
675 template <int dim, typename PolynomialType>
676 inline Tensor<1, dim>
678  const unsigned int i,
679  const Point<dim> & p) const
680 {
681  return compute_derivative<1>(i, p);
682 }
683 
684 
685 
686 template <int dim, typename PolynomialType>
687 inline Tensor<2, dim>
689  const unsigned int i,
690  const Point<dim> & p) const
691 {
692  return compute_derivative<2>(i, p);
693 }
694 
695 
696 
697 template <int dim, typename PolynomialType>
698 inline Tensor<3, dim>
700  const unsigned int i,
701  const Point<dim> & p) const
702 {
703  return compute_derivative<3>(i, p);
704 }
705 
706 
707 
708 template <int dim, typename PolynomialType>
709 inline Tensor<4, dim>
711  const unsigned int i,
712  const Point<dim> & p) const
713 {
714  return compute_derivative<4>(i, p);
715 }
716 
717 
718 
719 template <int dim>
720 template <int order>
723  const Point<dim> & p) const
724 {
725  std::array<unsigned int, dim> indices;
726  compute_index(i, indices);
727 
728  std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
729  for (unsigned int d = 0; d < dim; ++d)
730  polynomials[d][indices[d]].value(p(d), v[d]);
731 
732  Tensor<order, dim> derivative;
733  switch (order)
734  {
735  case 1:
736  {
737  Tensor<1, dim> &derivative_1 =
738  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
739  for (unsigned int d = 0; d < dim; ++d)
740  {
741  derivative_1[d] = 1.;
742  for (unsigned int x = 0; x < dim; ++x)
743  {
744  unsigned int x_order = 0;
745  if (d == x)
746  ++x_order;
747 
748  derivative_1[d] *= v[x][x_order];
749  }
750  }
751 
752  return derivative;
753  }
754  case 2:
755  {
756  Tensor<2, dim> &derivative_2 =
757  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
758  for (unsigned int d1 = 0; d1 < dim; ++d1)
759  for (unsigned int d2 = 0; d2 < dim; ++d2)
760  {
761  derivative_2[d1][d2] = 1.;
762  for (unsigned int x = 0; x < dim; ++x)
763  {
764  unsigned int x_order = 0;
765  if (d1 == x)
766  ++x_order;
767  if (d2 == x)
768  ++x_order;
769 
770  derivative_2[d1][d2] *= v[x][x_order];
771  }
772  }
773 
774  return derivative;
775  }
776  case 3:
777  {
778  Tensor<3, dim> &derivative_3 =
779  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
780  for (unsigned int d1 = 0; d1 < dim; ++d1)
781  for (unsigned int d2 = 0; d2 < dim; ++d2)
782  for (unsigned int d3 = 0; d3 < dim; ++d3)
783  {
784  derivative_3[d1][d2][d3] = 1.;
785  for (unsigned int x = 0; x < dim; ++x)
786  {
787  unsigned int x_order = 0;
788  if (d1 == x)
789  ++x_order;
790  if (d2 == x)
791  ++x_order;
792  if (d3 == x)
793  ++x_order;
794 
795  derivative_3[d1][d2][d3] *= v[x][x_order];
796  }
797  }
798 
799  return derivative;
800  }
801  case 4:
802  {
803  Tensor<4, dim> &derivative_4 =
804  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
805  for (unsigned int d1 = 0; d1 < dim; ++d1)
806  for (unsigned int d2 = 0; d2 < dim; ++d2)
807  for (unsigned int d3 = 0; d3 < dim; ++d3)
808  for (unsigned int d4 = 0; d4 < dim; ++d4)
809  {
810  derivative_4[d1][d2][d3][d4] = 1.;
811  for (unsigned int x = 0; x < dim; ++x)
812  {
813  unsigned int x_order = 0;
814  if (d1 == x)
815  ++x_order;
816  if (d2 == x)
817  ++x_order;
818  if (d3 == x)
819  ++x_order;
820  if (d4 == x)
821  ++x_order;
822 
823  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
824  }
825  }
826 
827  return derivative;
828  }
829  default:
830  {
831  Assert(false, ExcNotImplemented());
832  return derivative;
833  }
834  }
835 }
836 
837 
838 
839 template <>
840 template <int order>
843  const Point<0> &) const
844 {
845  AssertThrow(false, ExcNotImplemented());
846 
847  return {};
848 }
849 
850 
851 
852 template <int dim>
853 inline Tensor<1, dim>
855  const Point<dim> & p) const
856 {
857  return compute_derivative<1>(i, p);
858 }
859 
860 
861 
862 template <int dim>
863 inline Tensor<2, dim>
865  const Point<dim> & p) const
866 {
867  return compute_derivative<2>(i, p);
868 }
869 
870 
871 
872 template <int dim>
873 inline Tensor<3, dim>
875  const Point<dim> & p) const
876 {
877  return compute_derivative<3>(i, p);
878 }
879 
880 
881 
882 template <int dim>
883 inline Tensor<4, dim>
885  const Point<dim> & p) const
886 {
887  return compute_derivative<4>(i, p);
888 }
889 
890 
891 
892 template <int dim>
893 inline std::string
895 {
896  return "AnisotropicPolynomials";
897 }
898 
899 
900 
901 #endif // DOXYGEN
903 
904 #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)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
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:1423
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:369
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:32
std::string name() const override
static const unsigned int dimension
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
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