Reference documentation for deal.II version 9.3.3
\(\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
24#include <deal.II/base/point.h>
27#include <deal.II/base/tensor.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
37template <int dim>
39template <int dim>
41
75template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
77{
78public:
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> &
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
255protected:
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
321template <int dim>
323{
324public:
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
468private:
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
500template <int dim, typename PolynomialType>
501template <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
519template <int dim, typename PolynomialType>
520inline const std::vector<unsigned int> &
522{
523 return index_map;
524}
525
526
527template <int dim, typename PolynomialType>
528inline const std::vector<unsigned int> &
530{
531 return index_map_inverse;
532}
533
534
535template <int dim, typename PolynomialType>
536inline std::string
538{
539 return "TensorProductPolynomials";
540}
541
542
543template <int dim, typename PolynomialType>
544template <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
674template <>
675template <int order>
678 compute_derivative(const unsigned int, const Point<0> &) const
679{
681
682 return {};
683}
684
685
686
687template <int dim, typename PolynomialType>
688inline 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
698template <int dim, typename PolynomialType>
699inline 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
709template <int dim, typename PolynomialType>
710inline 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
720template <int dim, typename PolynomialType>
721inline 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
731template <int dim>
732template <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
851template <>
852template <int order>
855 const Point<0> &) const
856{
858
859 return {};
860}
861
862
863
864template <int dim>
865inline Tensor<1, dim>
867 const Point<dim> & p) const
868{
869 return compute_derivative<1>(i, p);
870}
871
872
873
874template <int dim>
875inline Tensor<2, dim>
877 const Point<dim> & p) const
878{
879 return compute_derivative<2>(i, p);
880}
881
882
883
884template <int dim>
885inline Tensor<3, dim>
887 const Point<dim> & p) const
888{
889 return compute_derivative<3>(i, p);
890}
891
892
893
894template <int dim>
895inline Tensor<4, dim>
897 const Point<dim> & p) const
898{
899 return compute_derivative<4>(i, p);
900}
901
902
903
904template <int dim>
905inline std::string
907{
908 return "AnisotropicPolynomials";
909}
910
911
912
913#endif // DOXYGEN
915
916#endif
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
double compute_value(const unsigned int i, const Point< dim > &p) const override
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double > > > &base_polynomials)
const std::vector< std::vector< Polynomials::Polynomial< double > > > polynomials
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double > > > &pols)
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
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< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
std::string name() const override
Tensor< 1, dim > compute_grad(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
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Definition: point.h:111
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
virtual Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
const std::vector< unsigned int > & get_numbering() const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map
std::vector< unsigned int > index_map_inverse
std::vector< PolynomialType > get_underlying_polynomials() const
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
static const unsigned int dimension
std::vector< PolynomialType > polynomials
void set_numbering(const std::vector< unsigned int > &renumber)
const std::vector< unsigned int > & get_numbering_inverse() const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
std::string name() const override
virtual Tensor< 4, dim > compute_4th_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:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
unsigned int compute_index()
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T fixed_power(const T t)
Definition: utilities.h:1081
Tensor< order, dim > compute_derivative(const unsigned int, const Point< dim > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108