Reference documentation for deal.II version Git d2841aca35 2021-05-14 14:24:54 -0600
\(\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 - 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 #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/ndarray.h>
24 #include <deal.II/base/point.h>
27 #include <deal.II/base/tensor.h>
28 #include <deal.II/base/utilities.h>
29 
30 #include <vector>
31 
33 
34 // Forward declarations for friends
35 // TODO: We may be able to modify these classes so they aren't
36 // required to be friends
37 template <int dim>
39 template <int dim>
41 
75 template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
77 {
78 public:
83  static const unsigned int dimension = dim;
84 
91  template <class Pol>
92  TensorProductPolynomials(const std::vector<Pol> &pols);
93 
97  void
98  output_indices(std::ostream &out) const;
99 
104  void
105  set_numbering(const std::vector<unsigned int> &renumber);
106 
110  const std::vector<unsigned int> &
111  get_numbering() const;
112 
116  const std::vector<unsigned int> &
117  get_numbering_inverse() const;
118 
131  void
132  evaluate(const Point<dim> & unit_point,
133  std::vector<double> & values,
134  std::vector<Tensor<1, dim>> &grads,
135  std::vector<Tensor<2, dim>> &grad_grads,
136  std::vector<Tensor<3, dim>> &third_derivatives,
137  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
138 
151  double
152  compute_value(const unsigned int i, const Point<dim> &p) const override;
153 
168  template <int order>
170  compute_derivative(const unsigned int i, const Point<dim> &p) const;
171 
175  virtual Tensor<1, dim>
176  compute_1st_derivative(const unsigned int i,
177  const Point<dim> & p) const override;
178 
182  virtual Tensor<2, dim>
183  compute_2nd_derivative(const unsigned int i,
184  const Point<dim> & p) const override;
185 
189  virtual Tensor<3, dim>
190  compute_3rd_derivative(const unsigned int i,
191  const Point<dim> & p) const override;
192 
196  virtual Tensor<4, dim>
197  compute_4th_derivative(const unsigned int i,
198  const Point<dim> & p) const override;
199 
213  compute_grad(const unsigned int i, const Point<dim> &p) const override;
214 
228  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
229 
233  std::string
234  name() const override;
235 
239  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
240  clone() const override;
241 
245  virtual std::size_t
246  memory_consumption() const override;
247 
248 protected:
252  std::vector<PolynomialType> polynomials;
253 
257  std::vector<unsigned int> index_map;
258 
262  std::vector<unsigned int> index_map_inverse;
263 
270  void
271  compute_index(const unsigned int i,
272  std::array<unsigned int, dim> &indices) const;
273 
279 
284  friend class TensorProductPolynomialsConst<dim>;
285 };
286 
287 
288 
314 template <int dim>
316 {
317 public:
334  const std::vector<std::vector<Polynomials::Polynomial<double>>>
335  &base_polynomials);
336 
350  void
351  evaluate(const Point<dim> & unit_point,
352  std::vector<double> & values,
353  std::vector<Tensor<1, dim>> &grads,
354  std::vector<Tensor<2, dim>> &grad_grads,
355  std::vector<Tensor<3, dim>> &third_derivatives,
356  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
357 
370  double
371  compute_value(const unsigned int i, const Point<dim> &p) const override;
372 
387  template <int order>
389  compute_derivative(const unsigned int i, const Point<dim> &p) const;
390 
394  virtual Tensor<1, dim>
395  compute_1st_derivative(const unsigned int i,
396  const Point<dim> & p) const override;
397 
401  virtual Tensor<2, dim>
402  compute_2nd_derivative(const unsigned int i,
403  const Point<dim> & p) const override;
404 
408  virtual Tensor<3, dim>
409  compute_3rd_derivative(const unsigned int i,
410  const Point<dim> & p) const override;
411 
415  virtual Tensor<4, dim>
416  compute_4th_derivative(const unsigned int i,
417  const Point<dim> & p) const override;
418 
432  compute_grad(const unsigned int i, const Point<dim> &p) const override;
433 
447  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
448 
452  std::string
453  name() const override;
454 
458  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
459  clone() const override;
460 
461 private:
465  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
466 
473  void
474  compute_index(const unsigned int i,
475  std::array<unsigned int, dim> &indices) const;
476 
480  static unsigned int
481  get_n_tensor_pols(
482  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
483 };
484 
487 #ifndef DOXYGEN
488 
489 
490 /* ---------------- template and inline functions ---------- */
491 
492 
493 template <int dim, typename PolynomialType>
494 template <class Pol>
496  const std::vector<Pol> &pols)
497  : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
498  , polynomials(pols.begin(), pols.end())
499  , index_map(this->n())
500  , index_map_inverse(this->n())
501 {
502  // per default set this index map to identity. This map can be changed by
503  // the user through the set_numbering() function
504  for (unsigned int i = 0; i < this->n(); ++i)
505  {
506  index_map[i] = i;
507  index_map_inverse[i] = i;
508  }
509 }
510 
511 
512 template <int dim, typename PolynomialType>
513 inline const std::vector<unsigned int> &
515 {
516  return index_map;
517 }
518 
519 
520 template <int dim, typename PolynomialType>
521 inline const std::vector<unsigned int> &
523 {
524  return index_map_inverse;
525 }
526 
527 
528 template <int dim, typename PolynomialType>
529 inline std::string
531 {
532  return "TensorProductPolynomials";
533 }
534 
535 
536 template <int dim, typename PolynomialType>
537 template <int order>
540  const unsigned int i,
541  const Point<dim> & p) const
542 {
543  std::array<unsigned int, dim> indices;
544  compute_index(i, indices);
545 
547  {
548  std::vector<double> tmp(5);
549  for (unsigned int d = 0; d < dim; ++d)
550  {
551  polynomials[indices[d]].value(p(d), tmp);
552  v[d][0] = tmp[0];
553  v[d][1] = tmp[1];
554  v[d][2] = tmp[2];
555  v[d][3] = tmp[3];
556  v[d][4] = tmp[4];
557  }
558  }
559 
560  Tensor<order, dim> derivative;
561  switch (order)
562  {
563  case 1:
564  {
565  Tensor<1, dim> &derivative_1 =
566  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
567  for (unsigned int d = 0; d < dim; ++d)
568  {
569  derivative_1[d] = 1.;
570  for (unsigned int x = 0; x < dim; ++x)
571  {
572  unsigned int x_order = 0;
573  if (d == x)
574  ++x_order;
575 
576  derivative_1[d] *= v[x][x_order];
577  }
578  }
579 
580  return derivative;
581  }
582  case 2:
583  {
584  Tensor<2, dim> &derivative_2 =
585  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
586  for (unsigned int d1 = 0; d1 < dim; ++d1)
587  for (unsigned int d2 = 0; d2 < dim; ++d2)
588  {
589  derivative_2[d1][d2] = 1.;
590  for (unsigned int x = 0; x < dim; ++x)
591  {
592  unsigned int x_order = 0;
593  if (d1 == x)
594  ++x_order;
595  if (d2 == x)
596  ++x_order;
597 
598  derivative_2[d1][d2] *= v[x][x_order];
599  }
600  }
601 
602  return derivative;
603  }
604  case 3:
605  {
606  Tensor<3, dim> &derivative_3 =
607  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
608  for (unsigned int d1 = 0; d1 < dim; ++d1)
609  for (unsigned int d2 = 0; d2 < dim; ++d2)
610  for (unsigned int d3 = 0; d3 < dim; ++d3)
611  {
612  derivative_3[d1][d2][d3] = 1.;
613  for (unsigned int x = 0; x < dim; ++x)
614  {
615  unsigned int x_order = 0;
616  if (d1 == x)
617  ++x_order;
618  if (d2 == x)
619  ++x_order;
620  if (d3 == x)
621  ++x_order;
622 
623  derivative_3[d1][d2][d3] *= v[x][x_order];
624  }
625  }
626 
627  return derivative;
628  }
629  case 4:
630  {
631  Tensor<4, dim> &derivative_4 =
632  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
633  for (unsigned int d1 = 0; d1 < dim; ++d1)
634  for (unsigned int d2 = 0; d2 < dim; ++d2)
635  for (unsigned int d3 = 0; d3 < dim; ++d3)
636  for (unsigned int d4 = 0; d4 < dim; ++d4)
637  {
638  derivative_4[d1][d2][d3][d4] = 1.;
639  for (unsigned int x = 0; x < dim; ++x)
640  {
641  unsigned int x_order = 0;
642  if (d1 == x)
643  ++x_order;
644  if (d2 == x)
645  ++x_order;
646  if (d3 == x)
647  ++x_order;
648  if (d4 == x)
649  ++x_order;
650 
651  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
652  }
653  }
654 
655  return derivative;
656  }
657  default:
658  {
659  Assert(false, ExcNotImplemented());
660  return derivative;
661  }
662  }
663 }
664 
665 
666 
667 template <>
668 template <int order>
671  compute_derivative(const unsigned int, const Point<0> &) const
672 {
673  AssertThrow(false, ExcNotImplemented());
674 
675  return {};
676 }
677 
678 
679 
680 template <int dim, typename PolynomialType>
681 inline Tensor<1, dim>
683  const unsigned int i,
684  const Point<dim> & p) const
685 {
686  return compute_derivative<1>(i, p);
687 }
688 
689 
690 
691 template <int dim, typename PolynomialType>
692 inline Tensor<2, dim>
694  const unsigned int i,
695  const Point<dim> & p) const
696 {
697  return compute_derivative<2>(i, p);
698 }
699 
700 
701 
702 template <int dim, typename PolynomialType>
703 inline Tensor<3, dim>
705  const unsigned int i,
706  const Point<dim> & p) const
707 {
708  return compute_derivative<3>(i, p);
709 }
710 
711 
712 
713 template <int dim, typename PolynomialType>
714 inline Tensor<4, dim>
716  const unsigned int i,
717  const Point<dim> & p) const
718 {
719  return compute_derivative<4>(i, p);
720 }
721 
722 
723 
724 template <int dim>
725 template <int order>
728  const Point<dim> & p) const
729 {
730  std::array<unsigned int, dim> indices;
731  compute_index(i, indices);
732 
733  std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
734  for (unsigned int d = 0; d < dim; ++d)
735  polynomials[d][indices[d]].value(p(d), v[d]);
736 
737  Tensor<order, dim> derivative;
738  switch (order)
739  {
740  case 1:
741  {
742  Tensor<1, dim> &derivative_1 =
743  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
744  for (unsigned int d = 0; d < dim; ++d)
745  {
746  derivative_1[d] = 1.;
747  for (unsigned int x = 0; x < dim; ++x)
748  {
749  unsigned int x_order = 0;
750  if (d == x)
751  ++x_order;
752 
753  derivative_1[d] *= v[x][x_order];
754  }
755  }
756 
757  return derivative;
758  }
759  case 2:
760  {
761  Tensor<2, dim> &derivative_2 =
762  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
763  for (unsigned int d1 = 0; d1 < dim; ++d1)
764  for (unsigned int d2 = 0; d2 < dim; ++d2)
765  {
766  derivative_2[d1][d2] = 1.;
767  for (unsigned int x = 0; x < dim; ++x)
768  {
769  unsigned int x_order = 0;
770  if (d1 == x)
771  ++x_order;
772  if (d2 == x)
773  ++x_order;
774 
775  derivative_2[d1][d2] *= v[x][x_order];
776  }
777  }
778 
779  return derivative;
780  }
781  case 3:
782  {
783  Tensor<3, dim> &derivative_3 =
784  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
785  for (unsigned int d1 = 0; d1 < dim; ++d1)
786  for (unsigned int d2 = 0; d2 < dim; ++d2)
787  for (unsigned int d3 = 0; d3 < dim; ++d3)
788  {
789  derivative_3[d1][d2][d3] = 1.;
790  for (unsigned int x = 0; x < dim; ++x)
791  {
792  unsigned int x_order = 0;
793  if (d1 == x)
794  ++x_order;
795  if (d2 == x)
796  ++x_order;
797  if (d3 == x)
798  ++x_order;
799 
800  derivative_3[d1][d2][d3] *= v[x][x_order];
801  }
802  }
803 
804  return derivative;
805  }
806  case 4:
807  {
808  Tensor<4, dim> &derivative_4 =
809  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
810  for (unsigned int d1 = 0; d1 < dim; ++d1)
811  for (unsigned int d2 = 0; d2 < dim; ++d2)
812  for (unsigned int d3 = 0; d3 < dim; ++d3)
813  for (unsigned int d4 = 0; d4 < dim; ++d4)
814  {
815  derivative_4[d1][d2][d3][d4] = 1.;
816  for (unsigned int x = 0; x < dim; ++x)
817  {
818  unsigned int x_order = 0;
819  if (d1 == x)
820  ++x_order;
821  if (d2 == x)
822  ++x_order;
823  if (d3 == x)
824  ++x_order;
825  if (d4 == x)
826  ++x_order;
827 
828  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
829  }
830  }
831 
832  return derivative;
833  }
834  default:
835  {
836  Assert(false, ExcNotImplemented());
837  return derivative;
838  }
839  }
840 }
841 
842 
843 
844 template <>
845 template <int order>
848  const Point<0> &) const
849 {
850  AssertThrow(false, ExcNotImplemented());
851 
852  return {};
853 }
854 
855 
856 
857 template <int dim>
858 inline Tensor<1, dim>
860  const Point<dim> & p) const
861 {
862  return compute_derivative<1>(i, p);
863 }
864 
865 
866 
867 template <int dim>
868 inline Tensor<2, dim>
870  const Point<dim> & p) const
871 {
872  return compute_derivative<2>(i, p);
873 }
874 
875 
876 
877 template <int dim>
878 inline Tensor<3, dim>
880  const Point<dim> & p) const
881 {
882  return compute_derivative<3>(i, p);
883 }
884 
885 
886 
887 template <int dim>
888 inline Tensor<4, dim>
890  const Point<dim> & p) const
891 {
892  return compute_derivative<4>(i, p);
893 }
894 
895 
896 
897 template <int dim>
898 inline std::string
900 {
901  return "AnisotropicPolynomials";
902 }
903 
904 
905 
906 #endif // DOXYGEN
908 
909 #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:1575
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:1081
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
#define Assert(cond, exc)
Definition: exceptions.h:1465
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:396
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
std::string name() const override
static const unsigned int dimension
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
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
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override