Reference documentation for deal.II version GIT relicensing-837-gaf73e64026 2024-06-15 22:40:02+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\}}\)
Loading...
Searching...
No Matches
tensor_product_polynomials.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tensor_product_polynomials_h
16#define dealii_tensor_product_polynomials_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/point.h>
26#include <deal.II/base/tensor.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
36template <int dim>
38template <int dim>
40
74template <int dim, typename PolynomialType = Polynomials::Polynomial<double>>
76{
77public:
82 static constexpr unsigned int dimension = dim;
83
90 template <class Pol>
91 TensorProductPolynomials(const std::vector<Pol> &pols);
92
96 void
97 output_indices(std::ostream &out) const;
98
103 void
104 set_numbering(const std::vector<unsigned int> &renumber);
105
109 const std::vector<unsigned int> &
111
115 const std::vector<unsigned int> &
117
130 void
131 evaluate(const Point<dim> &unit_point,
132 std::vector<double> &values,
133 std::vector<Tensor<1, dim>> &grads,
134 std::vector<Tensor<2, dim>> &grad_grads,
135 std::vector<Tensor<3, dim>> &third_derivatives,
136 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
137
150 double
151 compute_value(const unsigned int i, const Point<dim> &p) const override;
152
167 template <int order>
169 compute_derivative(const unsigned int i, const Point<dim> &p) const;
170
174 virtual Tensor<1, dim>
175 compute_1st_derivative(const unsigned int i,
176 const Point<dim> &p) const override;
177
181 virtual Tensor<2, dim>
182 compute_2nd_derivative(const unsigned int i,
183 const Point<dim> &p) const override;
184
188 virtual Tensor<3, dim>
189 compute_3rd_derivative(const unsigned int i,
190 const Point<dim> &p) const override;
191
195 virtual Tensor<4, dim>
196 compute_4th_derivative(const unsigned int i,
197 const Point<dim> &p) const override;
198
212 compute_grad(const unsigned int i, const Point<dim> &p) const override;
213
227 compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
228
232 std::string
233 name() const override;
234
238 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
239 clone() const override;
240
244 virtual std::size_t
245 memory_consumption() const override;
246
251 std::vector<PolynomialType>
253
254protected:
258 std::vector<PolynomialType> polynomials;
259
263 std::vector<unsigned int> index_map;
264
268 std::vector<unsigned int> index_map_inverse;
269
276 void
277 compute_index(const unsigned int i,
278 std::array<unsigned int, dim> &indices) const;
279
284 friend class TensorProductPolynomialsBubbles<dim>;
285
290 friend class TensorProductPolynomialsConst<dim>;
291};
292
293
294
320template <int dim>
322{
323public:
340 const std::vector<std::vector<Polynomials::Polynomial<double>>>
341 &base_polynomials);
342
347 void
348 set_numbering(const std::vector<unsigned int> &renumber);
349
353 const std::vector<unsigned int> &
354 get_numbering() const;
355
359 const std::vector<unsigned int> &
360 get_numbering_inverse() const;
361
375 void
376 evaluate(const Point<dim> &unit_point,
377 std::vector<double> &values,
378 std::vector<Tensor<1, dim>> &grads,
379 std::vector<Tensor<2, dim>> &grad_grads,
380 std::vector<Tensor<3, dim>> &third_derivatives,
381 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
382
395 double
396 compute_value(const unsigned int i, const Point<dim> &p) const override;
397
412 template <int order>
414 compute_derivative(const unsigned int i, const Point<dim> &p) const;
415
419 virtual Tensor<1, dim>
420 compute_1st_derivative(const unsigned int i,
421 const Point<dim> &p) const override;
422
426 virtual Tensor<2, dim>
427 compute_2nd_derivative(const unsigned int i,
428 const Point<dim> &p) const override;
429
433 virtual Tensor<3, dim>
434 compute_3rd_derivative(const unsigned int i,
435 const Point<dim> &p) const override;
436
440 virtual Tensor<4, dim>
441 compute_4th_derivative(const unsigned int i,
442 const Point<dim> &p) const override;
443
457 compute_grad(const unsigned int i, const Point<dim> &p) const override;
458
472 compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
473
477 std::string
478 name() const override;
479
483 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
484 clone() const override;
485
486private:
490 const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
491
495 std::vector<unsigned int> index_map;
496
500 std::vector<unsigned int> index_map_inverse;
501
508 void
509 compute_index(const unsigned int i,
510 std::array<unsigned int, dim> &indices) const;
511
515 static unsigned int
517 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
518};
519
522#ifndef DOXYGEN
523
524
525/* ---------------- template and inline functions ---------- */
526
527
528template <int dim, typename PolynomialType>
529template <class Pol>
531 const std::vector<Pol> &pols)
532 : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
533 , polynomials(pols.begin(), pols.end())
534 , index_map(this->n())
535 , index_map_inverse(this->n())
536{
537 // per default set this index map to identity. This map can be changed by
538 // the user through the set_numbering() function
539 for (unsigned int i = 0; i < this->n(); ++i)
540 {
541 index_map[i] = i;
542 index_map_inverse[i] = i;
543 }
544}
545
546
547template <int dim, typename PolynomialType>
548inline const std::vector<unsigned int> &
550{
551 return index_map;
552}
553
554
555template <int dim, typename PolynomialType>
556inline const std::vector<unsigned int> &
558{
559 return index_map_inverse;
560}
561
562
563template <int dim, typename PolynomialType>
564inline std::string
566{
567 return "TensorProductPolynomials";
568}
569
570
571template <int dim, typename PolynomialType>
572template <int order>
575 const unsigned int i,
576 const Point<dim> &p) const
577{
578 std::array<unsigned int, dim> indices;
579 compute_index(i, indices);
580
582 {
583 std::vector<double> tmp(5);
584 for (unsigned int d = 0; d < dim; ++d)
585 {
586 polynomials[indices[d]].value(p[d], tmp);
587 v[d][0] = tmp[0];
588 v[d][1] = tmp[1];
589 v[d][2] = tmp[2];
590 v[d][3] = tmp[3];
591 v[d][4] = tmp[4];
592 }
593 }
594
595 Tensor<order, dim> derivative;
596 switch (order)
597 {
598 case 1:
599 {
600 Tensor<1, dim> &derivative_1 =
601 *reinterpret_cast<Tensor<1, dim> *>(&derivative);
602 for (unsigned int d = 0; d < dim; ++d)
603 {
604 derivative_1[d] = 1.;
605 for (unsigned int x = 0; x < dim; ++x)
606 {
607 unsigned int x_order = 0;
608 if (d == x)
609 ++x_order;
610
611 derivative_1[d] *= v[x][x_order];
612 }
613 }
614
615 return derivative;
616 }
617 case 2:
618 {
619 Tensor<2, dim> &derivative_2 =
620 *reinterpret_cast<Tensor<2, dim> *>(&derivative);
621 for (unsigned int d1 = 0; d1 < dim; ++d1)
622 for (unsigned int d2 = 0; d2 < dim; ++d2)
623 {
624 derivative_2[d1][d2] = 1.;
625 for (unsigned int x = 0; x < dim; ++x)
626 {
627 unsigned int x_order = 0;
628 if (d1 == x)
629 ++x_order;
630 if (d2 == x)
631 ++x_order;
632
633 derivative_2[d1][d2] *= v[x][x_order];
634 }
635 }
636
637 return derivative;
638 }
639 case 3:
640 {
641 Tensor<3, dim> &derivative_3 =
642 *reinterpret_cast<Tensor<3, dim> *>(&derivative);
643 for (unsigned int d1 = 0; d1 < dim; ++d1)
644 for (unsigned int d2 = 0; d2 < dim; ++d2)
645 for (unsigned int d3 = 0; d3 < dim; ++d3)
646 {
647 derivative_3[d1][d2][d3] = 1.;
648 for (unsigned int x = 0; x < dim; ++x)
649 {
650 unsigned int x_order = 0;
651 if (d1 == x)
652 ++x_order;
653 if (d2 == x)
654 ++x_order;
655 if (d3 == x)
656 ++x_order;
657
658 derivative_3[d1][d2][d3] *= v[x][x_order];
659 }
660 }
661
662 return derivative;
663 }
664 case 4:
665 {
666 Tensor<4, dim> &derivative_4 =
667 *reinterpret_cast<Tensor<4, dim> *>(&derivative);
668 for (unsigned int d1 = 0; d1 < dim; ++d1)
669 for (unsigned int d2 = 0; d2 < dim; ++d2)
670 for (unsigned int d3 = 0; d3 < dim; ++d3)
671 for (unsigned int d4 = 0; d4 < dim; ++d4)
672 {
673 derivative_4[d1][d2][d3][d4] = 1.;
674 for (unsigned int x = 0; x < dim; ++x)
675 {
676 unsigned int x_order = 0;
677 if (d1 == x)
678 ++x_order;
679 if (d2 == x)
680 ++x_order;
681 if (d3 == x)
682 ++x_order;
683 if (d4 == x)
684 ++x_order;
685
686 derivative_4[d1][d2][d3][d4] *= v[x][x_order];
687 }
688 }
689
690 return derivative;
691 }
692 default:
693 {
695 return derivative;
696 }
697 }
698}
699
700
701
702template <>
703template <int order>
706 compute_derivative(const unsigned int, const Point<0> &) const
707{
709
710 return {};
711}
712
713
714
715template <int dim, typename PolynomialType>
716inline Tensor<1, dim>
718 const unsigned int i,
719 const Point<dim> &p) const
720{
721 return compute_derivative<1>(i, p);
722}
723
724
725
726template <int dim, typename PolynomialType>
727inline Tensor<2, dim>
729 const unsigned int i,
730 const Point<dim> &p) const
731{
732 return compute_derivative<2>(i, p);
733}
734
735
736
737template <int dim, typename PolynomialType>
738inline Tensor<3, dim>
740 const unsigned int i,
741 const Point<dim> &p) const
742{
743 return compute_derivative<3>(i, p);
744}
745
746
747
748template <int dim, typename PolynomialType>
749inline Tensor<4, dim>
751 const unsigned int i,
752 const Point<dim> &p) const
753{
754 return compute_derivative<4>(i, p);
755}
756
757
758
759template <int dim>
760template <int order>
763 const Point<dim> &p) const
764{
765 std::array<unsigned int, dim> indices;
766 compute_index(i, indices);
767
768 std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
769 for (unsigned int d = 0; d < dim; ++d)
770 polynomials[d][indices[d]].value(p[d], v[d]);
771
772 Tensor<order, dim> derivative;
773 switch (order)
774 {
775 case 1:
776 {
777 Tensor<1, dim> &derivative_1 =
778 *reinterpret_cast<Tensor<1, dim> *>(&derivative);
779 for (unsigned int d = 0; d < dim; ++d)
780 {
781 derivative_1[d] = 1.;
782 for (unsigned int x = 0; x < dim; ++x)
783 {
784 unsigned int x_order = 0;
785 if (d == x)
786 ++x_order;
787
788 derivative_1[d] *= v[x][x_order];
789 }
790 }
791
792 return derivative;
793 }
794 case 2:
795 {
796 Tensor<2, dim> &derivative_2 =
797 *reinterpret_cast<Tensor<2, dim> *>(&derivative);
798 for (unsigned int d1 = 0; d1 < dim; ++d1)
799 for (unsigned int d2 = 0; d2 < dim; ++d2)
800 {
801 derivative_2[d1][d2] = 1.;
802 for (unsigned int x = 0; x < dim; ++x)
803 {
804 unsigned int x_order = 0;
805 if (d1 == x)
806 ++x_order;
807 if (d2 == x)
808 ++x_order;
809
810 derivative_2[d1][d2] *= v[x][x_order];
811 }
812 }
813
814 return derivative;
815 }
816 case 3:
817 {
818 Tensor<3, dim> &derivative_3 =
819 *reinterpret_cast<Tensor<3, dim> *>(&derivative);
820 for (unsigned int d1 = 0; d1 < dim; ++d1)
821 for (unsigned int d2 = 0; d2 < dim; ++d2)
822 for (unsigned int d3 = 0; d3 < dim; ++d3)
823 {
824 derivative_3[d1][d2][d3] = 1.;
825 for (unsigned int x = 0; x < dim; ++x)
826 {
827 unsigned int x_order = 0;
828 if (d1 == x)
829 ++x_order;
830 if (d2 == x)
831 ++x_order;
832 if (d3 == x)
833 ++x_order;
834
835 derivative_3[d1][d2][d3] *= v[x][x_order];
836 }
837 }
838
839 return derivative;
840 }
841 case 4:
842 {
843 Tensor<4, dim> &derivative_4 =
844 *reinterpret_cast<Tensor<4, dim> *>(&derivative);
845 for (unsigned int d1 = 0; d1 < dim; ++d1)
846 for (unsigned int d2 = 0; d2 < dim; ++d2)
847 for (unsigned int d3 = 0; d3 < dim; ++d3)
848 for (unsigned int d4 = 0; d4 < dim; ++d4)
849 {
850 derivative_4[d1][d2][d3][d4] = 1.;
851 for (unsigned int x = 0; x < dim; ++x)
852 {
853 unsigned int x_order = 0;
854 if (d1 == x)
855 ++x_order;
856 if (d2 == x)
857 ++x_order;
858 if (d3 == x)
859 ++x_order;
860 if (d4 == x)
861 ++x_order;
862
863 derivative_4[d1][d2][d3][d4] *= v[x][x_order];
864 }
865 }
866
867 return derivative;
868 }
869 default:
870 {
872 return derivative;
873 }
874 }
875}
876
877
878
879template <>
880template <int order>
883 const Point<0> &) const
884{
886
887 return {};
888}
889
890
891
892template <int dim>
893inline Tensor<1, dim>
895 const Point<dim> &p) const
896{
897 return compute_derivative<1>(i, p);
898}
899
900
901
902template <int dim>
903inline Tensor<2, dim>
905 const Point<dim> &p) const
906{
907 return compute_derivative<2>(i, p);
908}
909
910
911
912template <int dim>
913inline Tensor<3, dim>
915 const Point<dim> &p) const
916{
917 return compute_derivative<3>(i, p);
918}
919
920
921
922template <int dim>
923inline Tensor<4, dim>
925 const Point<dim> &p) const
926{
927 return compute_derivative<4>(i, p);
928}
929
930
931
932template <int dim>
933inline std::string
935{
936 return "AnisotropicPolynomials";
937}
938
939
940
941#endif // DOXYGEN
943
944#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
void set_numbering(const std::vector< unsigned int > &renumber)
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)
std::vector< unsigned int > index_map
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
const std::vector< unsigned int > & get_numbering() const
std::string name() const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map_inverse
virtual Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
const std::vector< unsigned int > & get_numbering_inverse() const
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
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
static constexpr unsigned int dimension
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107