Reference documentation for deal.II version GIT relicensing-414-gae6ebd04e8 2024-04-17 07:40:05+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>>>
342
356 void
357 evaluate(const Point<dim> &unit_point,
358 std::vector<double> &values,
359 std::vector<Tensor<1, dim>> &grads,
360 std::vector<Tensor<2, dim>> &grad_grads,
361 std::vector<Tensor<3, dim>> &third_derivatives,
362 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
363
376 double
377 compute_value(const unsigned int i, const Point<dim> &p) const override;
378
393 template <int order>
395 compute_derivative(const unsigned int i, const Point<dim> &p) const;
396
400 virtual Tensor<1, dim>
401 compute_1st_derivative(const unsigned int i,
402 const Point<dim> &p) const override;
403
407 virtual Tensor<2, dim>
408 compute_2nd_derivative(const unsigned int i,
409 const Point<dim> &p) const override;
410
414 virtual Tensor<3, dim>
415 compute_3rd_derivative(const unsigned int i,
416 const Point<dim> &p) const override;
417
421 virtual Tensor<4, dim>
422 compute_4th_derivative(const unsigned int i,
423 const Point<dim> &p) const override;
424
438 compute_grad(const unsigned int i, const Point<dim> &p) const override;
439
453 compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
454
458 std::string
459 name() const override;
460
464 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
465 clone() const override;
466
467private:
471 const std::vector<std::vector<Polynomials::Polynomial<double>>> polynomials;
472
479 void
480 compute_index(const unsigned int i,
481 std::array<unsigned int, dim> &indices) const;
482
486 static unsigned int
488 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols);
489};
490
493#ifndef DOXYGEN
494
495
496/* ---------------- template and inline functions ---------- */
497
498
499template <int dim, typename PolynomialType>
500template <class Pol>
502 const std::vector<Pol> &pols)
503 : ScalarPolynomialsBase<dim>(1, Utilities::fixed_power<dim>(pols.size()))
504 , polynomials(pols.begin(), pols.end())
505 , index_map(this->n())
506 , index_map_inverse(this->n())
507{
508 // per default set this index map to identity. This map can be changed by
509 // the user through the set_numbering() function
510 for (unsigned int i = 0; i < this->n(); ++i)
511 {
512 index_map[i] = i;
513 index_map_inverse[i] = i;
514 }
515}
516
517
518template <int dim, typename PolynomialType>
519inline const std::vector<unsigned int> &
521{
522 return index_map;
523}
524
525
526template <int dim, typename PolynomialType>
527inline const std::vector<unsigned int> &
529{
530 return index_map_inverse;
531}
532
533
534template <int dim, typename PolynomialType>
535inline std::string
537{
538 return "TensorProductPolynomials";
539}
540
541
542template <int dim, typename PolynomialType>
543template <int order>
546 const unsigned int i,
547 const Point<dim> &p) const
548{
549 std::array<unsigned int, dim> indices;
550 compute_index(i, indices);
551
553 {
554 std::vector<double> tmp(5);
555 for (unsigned int d = 0; d < dim; ++d)
556 {
557 polynomials[indices[d]].value(p[d], tmp);
558 v[d][0] = tmp[0];
559 v[d][1] = tmp[1];
560 v[d][2] = tmp[2];
561 v[d][3] = tmp[3];
562 v[d][4] = tmp[4];
563 }
564 }
565
566 Tensor<order, dim> derivative;
567 switch (order)
568 {
569 case 1:
570 {
572 *reinterpret_cast<Tensor<1, dim> *>(&derivative);
573 for (unsigned int d = 0; d < dim; ++d)
574 {
575 derivative_1[d] = 1.;
576 for (unsigned int x = 0; x < dim; ++x)
577 {
578 unsigned int x_order = 0;
579 if (d == x)
580 ++x_order;
581
582 derivative_1[d] *= v[x][x_order];
583 }
584 }
585
586 return derivative;
587 }
588 case 2:
589 {
591 *reinterpret_cast<Tensor<2, dim> *>(&derivative);
592 for (unsigned int d1 = 0; d1 < dim; ++d1)
593 for (unsigned int d2 = 0; d2 < dim; ++d2)
594 {
595 derivative_2[d1][d2] = 1.;
596 for (unsigned int x = 0; x < dim; ++x)
597 {
598 unsigned int x_order = 0;
599 if (d1 == x)
600 ++x_order;
601 if (d2 == x)
602 ++x_order;
603
604 derivative_2[d1][d2] *= v[x][x_order];
605 }
606 }
607
608 return derivative;
609 }
610 case 3:
611 {
613 *reinterpret_cast<Tensor<3, dim> *>(&derivative);
614 for (unsigned int d1 = 0; d1 < dim; ++d1)
615 for (unsigned int d2 = 0; d2 < dim; ++d2)
616 for (unsigned int d3 = 0; d3 < dim; ++d3)
617 {
618 derivative_3[d1][d2][d3] = 1.;
619 for (unsigned int x = 0; x < dim; ++x)
620 {
621 unsigned int x_order = 0;
622 if (d1 == x)
623 ++x_order;
624 if (d2 == x)
625 ++x_order;
626 if (d3 == x)
627 ++x_order;
628
629 derivative_3[d1][d2][d3] *= v[x][x_order];
630 }
631 }
632
633 return derivative;
634 }
635 case 4:
636 {
638 *reinterpret_cast<Tensor<4, dim> *>(&derivative);
639 for (unsigned int d1 = 0; d1 < dim; ++d1)
640 for (unsigned int d2 = 0; d2 < dim; ++d2)
641 for (unsigned int d3 = 0; d3 < dim; ++d3)
642 for (unsigned int d4 = 0; d4 < dim; ++d4)
643 {
644 derivative_4[d1][d2][d3][d4] = 1.;
645 for (unsigned int x = 0; x < dim; ++x)
646 {
647 unsigned int x_order = 0;
648 if (d1 == x)
649 ++x_order;
650 if (d2 == x)
651 ++x_order;
652 if (d3 == x)
653 ++x_order;
654 if (d4 == x)
655 ++x_order;
656
657 derivative_4[d1][d2][d3][d4] *= v[x][x_order];
658 }
659 }
660
661 return derivative;
662 }
663 default:
664 {
666 return derivative;
667 }
668 }
669}
670
671
672
673template <>
674template <int order>
677 compute_derivative(const unsigned int, const Point<0> &) const
678{
680
681 return {};
682}
683
684
685
686template <int dim, typename PolynomialType>
687inline Tensor<1, dim>
689 const unsigned int i,
690 const Point<dim> &p) const
691{
692 return compute_derivative<1>(i, p);
693}
694
695
696
697template <int dim, typename PolynomialType>
698inline Tensor<2, dim>
700 const unsigned int i,
701 const Point<dim> &p) const
702{
703 return compute_derivative<2>(i, p);
704}
705
706
707
708template <int dim, typename PolynomialType>
709inline Tensor<3, dim>
711 const unsigned int i,
712 const Point<dim> &p) const
713{
714 return compute_derivative<3>(i, p);
715}
716
717
718
719template <int dim, typename PolynomialType>
720inline Tensor<4, dim>
722 const unsigned int i,
723 const Point<dim> &p) const
724{
725 return compute_derivative<4>(i, p);
726}
727
728
729
730template <int dim>
731template <int order>
734 const Point<dim> &p) const
735{
736 std::array<unsigned int, dim> indices;
737 compute_index(i, indices);
738
739 std::vector<std::vector<double>> v(dim, std::vector<double>(order + 1));
740 for (unsigned int d = 0; d < dim; ++d)
741 polynomials[d][indices[d]].value(p[d], v[d]);
742
743 Tensor<order, dim> derivative;
744 switch (order)
745 {
746 case 1:
747 {
749 *reinterpret_cast<Tensor<1, dim> *>(&derivative);
750 for (unsigned int d = 0; d < dim; ++d)
751 {
752 derivative_1[d] = 1.;
753 for (unsigned int x = 0; x < dim; ++x)
754 {
755 unsigned int x_order = 0;
756 if (d == x)
757 ++x_order;
758
759 derivative_1[d] *= v[x][x_order];
760 }
761 }
762
763 return derivative;
764 }
765 case 2:
766 {
768 *reinterpret_cast<Tensor<2, dim> *>(&derivative);
769 for (unsigned int d1 = 0; d1 < dim; ++d1)
770 for (unsigned int d2 = 0; d2 < dim; ++d2)
771 {
772 derivative_2[d1][d2] = 1.;
773 for (unsigned int x = 0; x < dim; ++x)
774 {
775 unsigned int x_order = 0;
776 if (d1 == x)
777 ++x_order;
778 if (d2 == x)
779 ++x_order;
780
781 derivative_2[d1][d2] *= v[x][x_order];
782 }
783 }
784
785 return derivative;
786 }
787 case 3:
788 {
790 *reinterpret_cast<Tensor<3, dim> *>(&derivative);
791 for (unsigned int d1 = 0; d1 < dim; ++d1)
792 for (unsigned int d2 = 0; d2 < dim; ++d2)
793 for (unsigned int d3 = 0; d3 < dim; ++d3)
794 {
795 derivative_3[d1][d2][d3] = 1.;
796 for (unsigned int x = 0; x < dim; ++x)
797 {
798 unsigned int x_order = 0;
799 if (d1 == x)
800 ++x_order;
801 if (d2 == x)
802 ++x_order;
803 if (d3 == x)
804 ++x_order;
805
806 derivative_3[d1][d2][d3] *= v[x][x_order];
807 }
808 }
809
810 return derivative;
811 }
812 case 4:
813 {
815 *reinterpret_cast<Tensor<4, dim> *>(&derivative);
816 for (unsigned int d1 = 0; d1 < dim; ++d1)
817 for (unsigned int d2 = 0; d2 < dim; ++d2)
818 for (unsigned int d3 = 0; d3 < dim; ++d3)
819 for (unsigned int d4 = 0; d4 < dim; ++d4)
820 {
821 derivative_4[d1][d2][d3][d4] = 1.;
822 for (unsigned int x = 0; x < dim; ++x)
823 {
824 unsigned int x_order = 0;
825 if (d1 == x)
826 ++x_order;
827 if (d2 == x)
828 ++x_order;
829 if (d3 == x)
830 ++x_order;
831 if (d4 == x)
832 ++x_order;
833
834 derivative_4[d1][d2][d3][d4] *= v[x][x_order];
835 }
836 }
837
838 return derivative;
839 }
840 default:
841 {
843 return derivative;
844 }
845 }
846}
847
848
849
850template <>
851template <int order>
854 const Point<0> &) const
855{
857
858 return {};
859}
860
861
862
863template <int dim>
864inline Tensor<1, dim>
866 const Point<dim> &p) const
867{
868 return compute_derivative<1>(i, p);
869}
870
871
872
873template <int dim>
874inline Tensor<2, dim>
876 const Point<dim> &p) const
877{
878 return compute_derivative<2>(i, p);
879}
880
881
882
883template <int dim>
884inline Tensor<3, dim>
886 const Point<dim> &p) const
887{
888 return compute_derivative<3>(i, p);
889}
890
891
892
893template <int dim>
894inline Tensor<4, dim>
896 const Point<dim> &p) const
897{
898 return compute_derivative<4>(i, p);
899}
900
901
902
903template <int dim>
904inline std::string
906{
907 return "AnisotropicPolynomials";
908}
909
910
911
912#endif // DOXYGEN
914
915#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
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
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)
friend class Tensor
Definition tensor.h:882
#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)