Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tensor_product_polynomials.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2018 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 #include <deal.II/base/polynomial.h>
25 #include <deal.II/base/tensor.h>
26 #include <deal.II/base/utilities.h>
27 
28 #include <vector>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
64 template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
66 {
67 public:
72  static const unsigned int dimension = dim;
73 
80  template <class Pol>
81  TensorProductPolynomials(const std::vector<Pol> &pols);
82 
86  void
87  output_indices(std::ostream &out) const;
88 
93  void
94  set_numbering(const std::vector<unsigned int> &renumber);
95 
99  const std::vector<unsigned int> &
100  get_numbering() const;
101 
105  const std::vector<unsigned int> &
106  get_numbering_inverse() const;
107 
120  void
121  compute(const Point<dim> & unit_point,
122  std::vector<double> & values,
123  std::vector<Tensor<1, dim>> &grads,
124  std::vector<Tensor<2, dim>> &grad_grads,
125  std::vector<Tensor<3, dim>> &third_derivatives,
126  std::vector<Tensor<4, dim>> &fourth_derivatives) const;
127 
140  double
141  compute_value(const unsigned int i, const Point<dim> &p) const;
142 
157  template <int order>
159  compute_derivative(const unsigned int i, const Point<dim> &p) const;
160 
174  compute_grad(const unsigned int i, const Point<dim> &p) const;
175 
189  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
190 
195  unsigned int
196  n() const;
197 
198 
199 protected:
203  std::vector<PolynomialType> polynomials;
204 
208  unsigned int n_tensor_pols;
209 
213  std::vector<unsigned int> index_map;
214 
218  std::vector<unsigned int> index_map_inverse;
219 
226  // fix to avoid compiler warnings about zero length arrays
227  void
228  compute_index(const unsigned int i,
229  unsigned int (&indices)[(dim > 0 ? dim : 1)]) const;
230 };
231 
232 
233 
261 template <int dim>
263 {
264 public:
278  const std::vector<std::vector<Polynomials::Polynomial<double>>>
279  &base_polynomials);
280 
294  void
295  compute(const Point<dim> & unit_point,
296  std::vector<double> & values,
297  std::vector<Tensor<1, dim>> &grads,
298  std::vector<Tensor<2, dim>> &grad_grads,
299  std::vector<Tensor<3, dim>> &third_derivatives,
300  std::vector<Tensor<4, dim>> &fourth_derivatives) const;
301 
314  double
315  compute_value(const unsigned int i, const Point<dim> &p) const;
316 
331  template <int order>
333  compute_derivative(const unsigned int i, const Point<dim> &p) const;
334 
348  compute_grad(const unsigned int i, const Point<dim> &p) const;
349 
363  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
364 
369  unsigned int
370  n() const;
371 
372 private:
376  const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
377 
382  const unsigned int n_tensor_pols;
383 
390  void
391  compute_index(const unsigned int i, unsigned int (&indices)[dim]) const;
392 
396  static unsigned int
398  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
399 };
400 
403 #ifndef DOXYGEN
404 
405 
406 /* ---------------- template and inline functions ---------- */
407 
408 
409 template <int dim, typename PolynomialType>
410 template <class Pol>
412  const std::vector<Pol> &pols)
413  : polynomials(pols.begin(), pols.end())
414  , n_tensor_pols(Utilities::fixed_power<dim>(pols.size()))
415  , index_map(n_tensor_pols)
416  , index_map_inverse(n_tensor_pols)
417 {
418  // per default set this index map to identity. This map can be changed by
419  // the user through the set_numbering() function
420  for (unsigned int i = 0; i < n_tensor_pols; ++i)
421  {
422  index_map[i] = i;
423  index_map_inverse[i] = i;
424  }
425 }
426 
427 
428 
429 template <int dim, typename PolynomialType>
430 inline unsigned int
432 {
433  if (dim == 0)
435  else
436  return n_tensor_pols;
437 }
438 
439 
440 
441 template <int dim, typename PolynomialType>
442 inline const std::vector<unsigned int> &
444 {
445  return index_map;
446 }
447 
448 
449 template <int dim, typename PolynomialType>
450 inline const std::vector<unsigned int> &
452 {
453  return index_map_inverse;
454 }
455 
456 template <int dim, typename PolynomialType>
457 template <int order>
460  const unsigned int i,
461  const Point<dim> & p) const
462 {
463  unsigned int indices[dim];
464  compute_index(i, indices);
465 
466  double v[dim][5];
467  {
468  std::vector<double> tmp(5);
469  for (unsigned int d = 0; d < dim; ++d)
470  {
471  polynomials[indices[d]].value(p(d), tmp);
472  v[d][0] = tmp[0];
473  v[d][1] = tmp[1];
474  v[d][2] = tmp[2];
475  v[d][3] = tmp[3];
476  v[d][4] = tmp[4];
477  }
478  }
479 
480  Tensor<order, dim> derivative;
481  switch (order)
482  {
483  case 1:
484  {
485  Tensor<1, dim> &derivative_1 =
486  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
487  for (unsigned int d = 0; d < dim; ++d)
488  {
489  derivative_1[d] = 1.;
490  for (unsigned int x = 0; x < dim; ++x)
491  {
492  unsigned int x_order = 0;
493  if (d == x)
494  ++x_order;
495 
496  derivative_1[d] *= v[x][x_order];
497  }
498  }
499 
500  return derivative;
501  }
502  case 2:
503  {
504  Tensor<2, dim> &derivative_2 =
505  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
506  for (unsigned int d1 = 0; d1 < dim; ++d1)
507  for (unsigned int d2 = 0; d2 < dim; ++d2)
508  {
509  derivative_2[d1][d2] = 1.;
510  for (unsigned int x = 0; x < dim; ++x)
511  {
512  unsigned int x_order = 0;
513  if (d1 == x)
514  ++x_order;
515  if (d2 == x)
516  ++x_order;
517 
518  derivative_2[d1][d2] *= v[x][x_order];
519  }
520  }
521 
522  return derivative;
523  }
524  case 3:
525  {
526  Tensor<3, dim> &derivative_3 =
527  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
528  for (unsigned int d1 = 0; d1 < dim; ++d1)
529  for (unsigned int d2 = 0; d2 < dim; ++d2)
530  for (unsigned int d3 = 0; d3 < dim; ++d3)
531  {
532  derivative_3[d1][d2][d3] = 1.;
533  for (unsigned int x = 0; x < dim; ++x)
534  {
535  unsigned int x_order = 0;
536  if (d1 == x)
537  ++x_order;
538  if (d2 == x)
539  ++x_order;
540  if (d3 == x)
541  ++x_order;
542 
543  derivative_3[d1][d2][d3] *= v[x][x_order];
544  }
545  }
546 
547  return derivative;
548  }
549  case 4:
550  {
551  Tensor<4, dim> &derivative_4 =
552  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
553  for (unsigned int d1 = 0; d1 < dim; ++d1)
554  for (unsigned int d2 = 0; d2 < dim; ++d2)
555  for (unsigned int d3 = 0; d3 < dim; ++d3)
556  for (unsigned int d4 = 0; d4 < dim; ++d4)
557  {
558  derivative_4[d1][d2][d3][d4] = 1.;
559  for (unsigned int x = 0; x < dim; ++x)
560  {
561  unsigned int x_order = 0;
562  if (d1 == x)
563  ++x_order;
564  if (d2 == x)
565  ++x_order;
566  if (d3 == x)
567  ++x_order;
568  if (d4 == x)
569  ++x_order;
570 
571  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
572  }
573  }
574 
575  return derivative;
576  }
577  default:
578  {
579  Assert(false, ExcNotImplemented());
580  return derivative;
581  }
582  }
583 }
584 
585 template <int dim>
586 template <int order>
589  const Point<dim> & p) const
590 {
591  unsigned int indices[dim];
592  compute_index(i, indices);
593 
594  std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
595  for (unsigned int d = 0; d < dim; ++d)
596  polynomials[d][indices[d]].value(p(d), v[d]);
597 
598  Tensor<order, dim> derivative;
599  switch (order)
600  {
601  case 1:
602  {
603  Tensor<1, dim> &derivative_1 =
604  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
605  for (unsigned int d = 0; d < dim; ++d)
606  {
607  derivative_1[d] = 1.;
608  for (unsigned int x = 0; x < dim; ++x)
609  {
610  unsigned int x_order = 0;
611  if (d == x)
612  ++x_order;
613 
614  derivative_1[d] *= v[x][x_order];
615  }
616  }
617 
618  return derivative;
619  }
620  case 2:
621  {
622  Tensor<2, dim> &derivative_2 =
623  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
624  for (unsigned int d1 = 0; d1 < dim; ++d1)
625  for (unsigned int d2 = 0; d2 < dim; ++d2)
626  {
627  derivative_2[d1][d2] = 1.;
628  for (unsigned int x = 0; x < dim; ++x)
629  {
630  unsigned int x_order = 0;
631  if (d1 == x)
632  ++x_order;
633  if (d2 == x)
634  ++x_order;
635 
636  derivative_2[d1][d2] *= v[x][x_order];
637  }
638  }
639 
640  return derivative;
641  }
642  case 3:
643  {
644  Tensor<3, dim> &derivative_3 =
645  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
646  for (unsigned int d1 = 0; d1 < dim; ++d1)
647  for (unsigned int d2 = 0; d2 < dim; ++d2)
648  for (unsigned int d3 = 0; d3 < dim; ++d3)
649  {
650  derivative_3[d1][d2][d3] = 1.;
651  for (unsigned int x = 0; x < dim; ++x)
652  {
653  unsigned int x_order = 0;
654  if (d1 == x)
655  ++x_order;
656  if (d2 == x)
657  ++x_order;
658  if (d3 == x)
659  ++x_order;
660 
661  derivative_3[d1][d2][d3] *= v[x][x_order];
662  }
663  }
664 
665  return derivative;
666  }
667  case 4:
668  {
669  Tensor<4, dim> &derivative_4 =
670  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
671  for (unsigned int d1 = 0; d1 < dim; ++d1)
672  for (unsigned int d2 = 0; d2 < dim; ++d2)
673  for (unsigned int d3 = 0; d3 < dim; ++d3)
674  for (unsigned int d4 = 0; d4 < dim; ++d4)
675  {
676  derivative_4[d1][d2][d3][d4] = 1.;
677  for (unsigned int x = 0; x < dim; ++x)
678  {
679  unsigned int x_order = 0;
680  if (d1 == x)
681  ++x_order;
682  if (d2 == x)
683  ++x_order;
684  if (d3 == x)
685  ++x_order;
686  if (d4 == x)
687  ++x_order;
688 
689  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
690  }
691  }
692 
693  return derivative;
694  }
695  default:
696  {
697  Assert(false, ExcNotImplemented());
698  return derivative;
699  }
700  }
701 }
702 
703 
704 
705 #endif // DOXYGEN
706 DEAL_II_NAMESPACE_CLOSE
707 
708 #endif
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
void compute(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
unsigned int n() const
static const unsigned int invalid_unsigned_int
Definition: types.h:173
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)
void compute(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
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
TensorProductPolynomials(const std::vector< Pol > &pols)
void compute_index(const unsigned int i, unsigned int(&indices)[dim]) const
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
#define Assert(cond, exc)
Definition: exceptions.h:1407
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)
const std::vector< unsigned int > & get_numbering() const
Definition: cuda.h:31
static const unsigned int dimension
std::vector< unsigned int > index_map
double compute_value(const unsigned int i, const Point< dim > &p) const
const std::vector< unsigned int > & get_numbering_inverse() const
static ::ExceptionBase & ExcNotImplemented()
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
void compute_index(const unsigned int i, unsigned int(&indices)[(dim > 0 ? dim :1)]) const
double compute_value(const unsigned int i, const Point< dim > &p) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const