Reference documentation for deal.II version GIT 58febcd5cf 2023-09-30 20:00:01+00:00
\(\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 - 2022 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 constexpr 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> &
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 
252  std::vector<PolynomialType>
254 
255 protected:
259  std::vector<PolynomialType> polynomials;
260 
264  std::vector<unsigned int> index_map;
265 
269  std::vector<unsigned int> index_map_inverse;
270 
277  void
278  compute_index(const unsigned int i,
279  std::array<unsigned int, dim> &indices) const;
280 
285  friend class TensorProductPolynomialsBubbles<dim>;
286 
291  friend class TensorProductPolynomialsConst<dim>;
292 };
293 
294 
295 
321 template <int dim>
323 {
324 public:
341  const std::vector<std::vector<Polynomials::Polynomial<double>>>
342  &base_polynomials);
343 
357  void
358  evaluate(const Point<dim> &unit_point,
359  std::vector<double> &values,
360  std::vector<Tensor<1, dim>> &grads,
361  std::vector<Tensor<2, dim>> &grad_grads,
362  std::vector<Tensor<3, dim>> &third_derivatives,
363  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
364 
377  double
378  compute_value(const unsigned int i, const Point<dim> &p) const override;
379 
394  template <int order>
396  compute_derivative(const unsigned int i, const Point<dim> &p) const;
397 
401  virtual Tensor<1, dim>
402  compute_1st_derivative(const unsigned int i,
403  const Point<dim> &p) const override;
404 
408  virtual Tensor<2, dim>
409  compute_2nd_derivative(const unsigned int i,
410  const Point<dim> &p) const override;
411 
415  virtual Tensor<3, dim>
416  compute_3rd_derivative(const unsigned int i,
417  const Point<dim> &p) const override;
418 
422  virtual Tensor<4, dim>
423  compute_4th_derivative(const unsigned int i,
424  const Point<dim> &p) const override;
425 
439  compute_grad(const unsigned int i, const Point<dim> &p) const override;
440 
454  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
455 
459  std::string
460  name() const override;
461 
465  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
466  clone() const override;
467 
468 private:
472  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
473 
480  void
481  compute_index(const unsigned int i,
482  std::array<unsigned int, dim> &indices) const;
483 
487  static unsigned int
489  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
490 };
491 
494 #ifndef DOXYGEN
495 
496 
497 /* ---------------- template and inline functions ---------- */
498 
499 
500 template <int dim, typename PolynomialType>
501 template <class Pol>
503  const std::vector<Pol> &pols)
504  : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
505  , polynomials(pols.begin(), pols.end())
506  , index_map(this->n())
507  , index_map_inverse(this->n())
508 {
509  // per default set this index map to identity. This map can be changed by
510  // the user through the set_numbering() function
511  for (unsigned int i = 0; i < this->n(); ++i)
512  {
513  index_map[i] = i;
514  index_map_inverse[i] = i;
515  }
516 }
517 
518 
519 template <int dim, typename PolynomialType>
520 inline const std::vector<unsigned int> &
522 {
523  return index_map;
524 }
525 
526 
527 template <int dim, typename PolynomialType>
528 inline const std::vector<unsigned int> &
530 {
531  return index_map_inverse;
532 }
533 
534 
535 template <int dim, typename PolynomialType>
536 inline std::string
538 {
539  return "TensorProductPolynomials";
540 }
541 
542 
543 template <int dim, typename PolynomialType>
544 template <int order>
547  const unsigned int i,
548  const Point<dim> &p) const
549 {
550  std::array<unsigned int, dim> indices;
551  compute_index(i, indices);
552 
554  {
555  std::vector<double> tmp(5);
556  for (unsigned int d = 0; d < dim; ++d)
557  {
558  polynomials[indices[d]].value(p(d), tmp);
559  v[d][0] = tmp[0];
560  v[d][1] = tmp[1];
561  v[d][2] = tmp[2];
562  v[d][3] = tmp[3];
563  v[d][4] = tmp[4];
564  }
565  }
566 
567  Tensor<order, dim> derivative;
568  switch (order)
569  {
570  case 1:
571  {
572  Tensor<1, dim> &derivative_1 =
573  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
574  for (unsigned int d = 0; d < dim; ++d)
575  {
576  derivative_1[d] = 1.;
577  for (unsigned int x = 0; x < dim; ++x)
578  {
579  unsigned int x_order = 0;
580  if (d == x)
581  ++x_order;
582 
583  derivative_1[d] *= v[x][x_order];
584  }
585  }
586 
587  return derivative;
588  }
589  case 2:
590  {
591  Tensor<2, dim> &derivative_2 =
592  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
593  for (unsigned int d1 = 0; d1 < dim; ++d1)
594  for (unsigned int d2 = 0; d2 < dim; ++d2)
595  {
596  derivative_2[d1][d2] = 1.;
597  for (unsigned int x = 0; x < dim; ++x)
598  {
599  unsigned int x_order = 0;
600  if (d1 == x)
601  ++x_order;
602  if (d2 == x)
603  ++x_order;
604 
605  derivative_2[d1][d2] *= v[x][x_order];
606  }
607  }
608 
609  return derivative;
610  }
611  case 3:
612  {
613  Tensor<3, dim> &derivative_3 =
614  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
615  for (unsigned int d1 = 0; d1 < dim; ++d1)
616  for (unsigned int d2 = 0; d2 < dim; ++d2)
617  for (unsigned int d3 = 0; d3 < dim; ++d3)
618  {
619  derivative_3[d1][d2][d3] = 1.;
620  for (unsigned int x = 0; x < dim; ++x)
621  {
622  unsigned int x_order = 0;
623  if (d1 == x)
624  ++x_order;
625  if (d2 == x)
626  ++x_order;
627  if (d3 == x)
628  ++x_order;
629 
630  derivative_3[d1][d2][d3] *= v[x][x_order];
631  }
632  }
633 
634  return derivative;
635  }
636  case 4:
637  {
638  Tensor<4, dim> &derivative_4 =
639  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
640  for (unsigned int d1 = 0; d1 < dim; ++d1)
641  for (unsigned int d2 = 0; d2 < dim; ++d2)
642  for (unsigned int d3 = 0; d3 < dim; ++d3)
643  for (unsigned int d4 = 0; d4 < dim; ++d4)
644  {
645  derivative_4[d1][d2][d3][d4] = 1.;
646  for (unsigned int x = 0; x < dim; ++x)
647  {
648  unsigned int x_order = 0;
649  if (d1 == x)
650  ++x_order;
651  if (d2 == x)
652  ++x_order;
653  if (d3 == x)
654  ++x_order;
655  if (d4 == x)
656  ++x_order;
657 
658  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
659  }
660  }
661 
662  return derivative;
663  }
664  default:
665  {
666  Assert(false, ExcNotImplemented());
667  return derivative;
668  }
669  }
670 }
671 
672 
673 
674 template <>
675 template <int order>
678  compute_derivative(const unsigned int, const Point<0> &) const
679 {
680  AssertThrow(false, ExcNotImplemented());
681 
682  return {};
683 }
684 
685 
686 
687 template <int dim, typename PolynomialType>
688 inline Tensor<1, dim>
690  const unsigned int i,
691  const Point<dim> &p) const
692 {
693  return compute_derivative<1>(i, p);
694 }
695 
696 
697 
698 template <int dim, typename PolynomialType>
699 inline Tensor<2, dim>
701  const unsigned int i,
702  const Point<dim> &p) const
703 {
704  return compute_derivative<2>(i, p);
705 }
706 
707 
708 
709 template <int dim, typename PolynomialType>
710 inline Tensor<3, dim>
712  const unsigned int i,
713  const Point<dim> &p) const
714 {
715  return compute_derivative<3>(i, p);
716 }
717 
718 
719 
720 template <int dim, typename PolynomialType>
721 inline Tensor<4, dim>
723  const unsigned int i,
724  const Point<dim> &p) const
725 {
726  return compute_derivative<4>(i, p);
727 }
728 
729 
730 
731 template <int dim>
732 template <int order>
735  const Point<dim> &p) const
736 {
737  std::array<unsigned int, dim> indices;
738  compute_index(i, indices);
739 
740  std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
741  for (unsigned int d = 0; d < dim; ++d)
742  polynomials[d][indices[d]].value(p(d), v[d]);
743 
744  Tensor<order, dim> derivative;
745  switch (order)
746  {
747  case 1:
748  {
749  Tensor<1, dim> &derivative_1 =
750  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
751  for (unsigned int d = 0; d < dim; ++d)
752  {
753  derivative_1[d] = 1.;
754  for (unsigned int x = 0; x < dim; ++x)
755  {
756  unsigned int x_order = 0;
757  if (d == x)
758  ++x_order;
759 
760  derivative_1[d] *= v[x][x_order];
761  }
762  }
763 
764  return derivative;
765  }
766  case 2:
767  {
768  Tensor<2, dim> &derivative_2 =
769  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
770  for (unsigned int d1 = 0; d1 < dim; ++d1)
771  for (unsigned int d2 = 0; d2 < dim; ++d2)
772  {
773  derivative_2[d1][d2] = 1.;
774  for (unsigned int x = 0; x < dim; ++x)
775  {
776  unsigned int x_order = 0;
777  if (d1 == x)
778  ++x_order;
779  if (d2 == x)
780  ++x_order;
781 
782  derivative_2[d1][d2] *= v[x][x_order];
783  }
784  }
785 
786  return derivative;
787  }
788  case 3:
789  {
790  Tensor<3, dim> &derivative_3 =
791  *reinterpret_cast<Tensor<3, 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  {
796  derivative_3[d1][d2][d3] = 1.;
797  for (unsigned int x = 0; x < dim; ++x)
798  {
799  unsigned int x_order = 0;
800  if (d1 == x)
801  ++x_order;
802  if (d2 == x)
803  ++x_order;
804  if (d3 == x)
805  ++x_order;
806 
807  derivative_3[d1][d2][d3] *= v[x][x_order];
808  }
809  }
810 
811  return derivative;
812  }
813  case 4:
814  {
815  Tensor<4, dim> &derivative_4 =
816  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
817  for (unsigned int d1 = 0; d1 < dim; ++d1)
818  for (unsigned int d2 = 0; d2 < dim; ++d2)
819  for (unsigned int d3 = 0; d3 < dim; ++d3)
820  for (unsigned int d4 = 0; d4 < dim; ++d4)
821  {
822  derivative_4[d1][d2][d3][d4] = 1.;
823  for (unsigned int x = 0; x < dim; ++x)
824  {
825  unsigned int x_order = 0;
826  if (d1 == x)
827  ++x_order;
828  if (d2 == x)
829  ++x_order;
830  if (d3 == x)
831  ++x_order;
832  if (d4 == x)
833  ++x_order;
834 
835  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
836  }
837  }
838 
839  return derivative;
840  }
841  default:
842  {
843  Assert(false, ExcNotImplemented());
844  return derivative;
845  }
846  }
847 }
848 
849 
850 
851 template <>
852 template <int order>
855  const Point<0> &) const
856 {
857  AssertThrow(false, ExcNotImplemented());
858 
859  return {};
860 }
861 
862 
863 
864 template <int dim>
865 inline Tensor<1, dim>
867  const Point<dim> &p) const
868 {
869  return compute_derivative<1>(i, p);
870 }
871 
872 
873 
874 template <int dim>
875 inline Tensor<2, dim>
877  const Point<dim> &p) const
878 {
879  return compute_derivative<2>(i, p);
880 }
881 
882 
883 
884 template <int dim>
885 inline Tensor<3, dim>
887  const Point<dim> &p) const
888 {
889  return compute_derivative<3>(i, p);
890 }
891 
892 
893 
894 template <int dim>
895 inline Tensor<4, dim>
897  const Point<dim> &p) const
898 {
899  return compute_derivative<4>(i, p);
900 }
901 
902 
903 
904 template <int dim>
905 inline std::string
907 {
908  return "AnisotropicPolynomials";
909 }
910 
911 
912 
913 #endif // DOXYGEN
915 
916 #endif
double compute_value(const unsigned int i, const Point< dim > &p) const override
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
const std::vector< std::vector< Polynomials::Polynomial< double > > > polynomials
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
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 override
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
std::string name() const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
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 Tensor< 2, dim > compute_2nd_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
Definition: point.h:112
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
void output_indices(std::ostream &out) const
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
double compute_value(const unsigned int i, const Point< dim > &p) const override
virtual std::size_t memory_consumption() const override
const std::vector< unsigned int > & get_numbering_inverse() const
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
const std::vector< unsigned int > & get_numbering() const
virtual Tensor< 3, dim > compute_3rd_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
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
std::vector< unsigned int > index_map
virtual Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map_inverse
std::vector< PolynomialType > get_underlying_polynomials() const
std::vector< PolynomialType > polynomials
void set_numbering(const std::vector< unsigned int > &renumber)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
std::string name() const override
static constexpr unsigned int dimension
virtual Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
TensorProductPolynomials(const std::vector< Pol > &pols)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T fixed_power(const T t)
Definition: utilities.h:975
Tensor< order, dim > compute_derivative(const unsigned int, const Point< dim > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108