Reference documentation for deal.II version Git 42bfb73b5a 2019-11-14 07:57:06 +0100
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tensor.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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_h
17 #define dealii_tensor_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/exceptions.h>
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/std_cxx14/utility.h>
24 #include <deal.II/base/table_indices.h>
25 #include <deal.II/base/template_constraints.h>
26 #include <deal.II/base/tensor_accessors.h>
27 #include <deal.II/base/utilities.h>
28 
29 #ifdef DEAL_II_WITH_ADOLC
30 # include <adolc/adouble.h> // Taped double
31 #endif
32 
33 #include <cmath>
34 #include <ostream>
35 #include <vector>
36 
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 // Forward declarations:
41 #ifndef DOXYGEN
42 template <int dim, typename Number>
43 class Point;
44 template <int rank_, int dim, typename Number = double>
45 class Tensor;
46 template <typename Number>
47 class Vector;
48 namespace Differentiation
49 {
50  namespace SD
51  {
52  class Expression;
53  }
54 } // namespace Differentiation
55 #endif
56 
57 #ifndef DOXYGEN
58 // Overload invalid tensor types of negative rank that come up during
59 // overload resolution of operator* and related contraction variants.
60 template <int dim, typename Number>
61 class Tensor<-2, dim, Number>
62 {};
63 
64 template <int dim, typename Number>
65 class Tensor<-1, dim, Number>
66 {};
67 #endif /* DOXYGEN */
68 
69 
100 template <int dim, typename Number>
101 class Tensor<0, dim, Number>
102 {
103 public:
112  static constexpr unsigned int dimension = dim;
113 
117  static constexpr unsigned int rank = 0;
118 
122  static constexpr unsigned int n_independent_components = 1;
123 
133 
138  using value_type = Number;
139 
145  using array_type = Number;
146 
152  constexpr DEAL_II_CUDA_HOST_DEV
153  Tensor();
154 
162  template <typename OtherNumber>
163  constexpr DEAL_II_CUDA_HOST_DEV
164  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
165 
171  template <typename OtherNumber>
172  constexpr DEAL_II_CUDA_HOST_DEV
173  Tensor(const OtherNumber &initializer);
174 
178  Number *
179  begin_raw();
180 
184  const Number *
185  begin_raw() const;
186 
190  Number *
191  end_raw();
192 
197  const Number *
198  end_raw() const;
199 
209  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator Number &();
210 
219  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator const Number &() const;
220 
228  template <typename OtherNumber>
229  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
230  operator=(const Tensor<0, dim, OtherNumber> &rhs);
231 
232 #ifdef __INTEL_COMPILER
233 
241  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
242  operator=(const Tensor<0, dim, Number> &rhs);
243 #endif
244 
251  template <typename OtherNumber>
252  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
253  operator=(const OtherNumber &d);
254 
258  template <typename OtherNumber>
259  DEAL_II_CONSTEXPR bool
260  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
261 
265  template <typename OtherNumber>
266  constexpr bool
267  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
268 
274  template <typename OtherNumber>
275  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
276  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
277 
283  template <typename OtherNumber>
284  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
285  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
286 
292  template <typename OtherNumber>
293  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
294  operator*=(const OtherNumber &factor);
295 
301  template <typename OtherNumber>
302  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
303  operator/=(const OtherNumber &factor);
304 
310  constexpr DEAL_II_CUDA_HOST_DEV Tensor
311  operator-() const;
312 
325  DEAL_II_CONSTEXPR void
326  clear();
327 
333  real_type
334  norm() const;
335 
342  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV real_type
343  norm_square() const;
344 
349  template <class Archive>
350  void
351  serialize(Archive &ar, const unsigned int version);
352 
357  using tensor_type = Number;
358 
359 private:
363  Number value;
364 
368  template <typename OtherNumber>
369  void
370  unroll_recursion(Vector<OtherNumber> &result,
371  unsigned int & start_index) const;
372 
373  // Allow an arbitrary Tensor to access the underlying values.
374  template <int, int, typename>
375  friend class Tensor;
376 };
377 
378 
379 
421 template <int rank_, int dim, typename Number>
422 class Tensor
423 {
424 public:
433  static constexpr unsigned int dimension = dim;
434 
438  static constexpr unsigned int rank = rank_;
439 
444  static constexpr unsigned int n_independent_components =
445  Tensor<rank_ - 1, dim>::n_independent_components * dim;
446 
452  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
453 
458  using array_type =
459  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
460 
466  constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
467  Tensor() = default;
468 
474  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
475  const array_type &initializer);
476 
484  template <typename OtherNumber>
485  constexpr DEAL_II_CUDA_HOST_DEV
486  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
487 
491  template <typename OtherNumber>
492  constexpr Tensor(
493  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
494 
498  template <typename OtherNumber>
499  constexpr
500  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
501 
507  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV value_type &
508  operator[](const unsigned int i);
509 
515  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
516  operator[](const unsigned int i) const;
517 
521  DEAL_II_CONSTEXPR const Number &
522  operator[](const TableIndices<rank_> &indices) const;
523 
527  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
528 
532  Number *
533  begin_raw();
534 
538  const Number *
539  begin_raw() const;
540 
544  Number *
545  end_raw();
546 
550  const Number *
551  end_raw() const;
552 
560  template <typename OtherNumber>
561  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
562  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
563 
570  DEAL_II_CONSTEXPR Tensor &
571  operator=(const Number &d);
572 
576  template <typename OtherNumber>
577  DEAL_II_CONSTEXPR bool
579 
583  template <typename OtherNumber>
584  constexpr bool
586 
592  template <typename OtherNumber>
593  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
594  operator+=(const Tensor<rank_, dim, OtherNumber> &);
595 
601  template <typename OtherNumber>
602  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
603  operator-=(const Tensor<rank_, dim, OtherNumber> &);
604 
611  template <typename OtherNumber>
612  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
613  operator*=(const OtherNumber &factor);
614 
620  template <typename OtherNumber>
621  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
622  operator/=(const OtherNumber &factor);
623 
629  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor
630  operator-() const;
631 
644  DEAL_II_CONSTEXPR void
645  clear();
646 
654  DEAL_II_CUDA_HOST_DEV
656  norm() const;
657 
664  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV
666  norm_square() const;
667 
675  template <typename OtherNumber>
676  void
677  unroll(Vector<OtherNumber> &result) const;
678 
683  static DEAL_II_CONSTEXPR unsigned int
684  component_to_unrolled_index(const TableIndices<rank_> &indices);
685 
690  static DEAL_II_CONSTEXPR TableIndices<rank_>
691  unrolled_to_component_indices(const unsigned int i);
692 
697  static constexpr std::size_t
698  memory_consumption();
699 
704  template <class Archive>
705  void
706  serialize(Archive &ar, const unsigned int version);
707 
713 
714 private:
718  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
719  // ... avoid a compiler warning in case of dim == 0 and ensure that the
720  // array always has positive size.
721 
725  template <typename OtherNumber>
726  void
727  unroll_recursion(Vector<OtherNumber> &result,
728  unsigned int & start_index) const;
729 
736  template <typename ArrayLike, std::size_t... Indices>
737  constexpr DEAL_II_CUDA_HOST_DEV
738  Tensor(const ArrayLike &initializer, std_cxx14::index_sequence<Indices...>);
739 
740  // Allow an arbitrary Tensor to access the underlying values.
741  template <int, int, typename>
742  friend class Tensor;
743 
744  // Point is allowed access to the coordinates. This is supposed to improve
745  // speed.
746  friend class Point<dim, Number>;
747 };
748 
749 
750 namespace internal
751 {
756  template <int rank, int dim, typename T>
757  struct NumberType<Tensor<rank, dim, T>>
758  {
759  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
760  value(const Tensor<rank, dim, T> &t)
761  {
762  return t;
763  }
764 
765  static DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE Tensor<rank, dim, T>
766  value(const T &t)
767  {
769  tmp = t;
770  return tmp;
771  }
772  };
773 } // namespace internal
774 
775 
776 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
777 
778 
779 template <int dim, typename Number>
780 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
782  // Some auto-differentiable numbers need explicit
783  // zero initialization such as adtl::adouble.
784  : Tensor{0.0}
785 {}
786 
787 
788 
789 template <int dim, typename Number>
790 template <typename OtherNumber>
791 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
792  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
793  : value(internal::NumberType<Number>::value(initializer))
794 {}
795 
796 
797 
798 template <int dim, typename Number>
799 template <typename OtherNumber>
800 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
802  : Tensor{p.value}
803 {}
804 
805 
806 
807 template <int dim, typename Number>
808 inline Number *
810 {
811  return std::addressof(value);
812 }
813 
814 
815 
816 template <int dim, typename Number>
817 inline const Number *
819 {
820  return std::addressof(value);
821 }
822 
823 
824 
825 template <int dim, typename Number>
826 inline Number *
828 {
830 }
831 
832 
833 
834 template <int dim, typename Number>
835 const Number *
837 {
839 }
840 
841 
842 
843 template <int dim, typename Number>
844 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
845  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
846 {
847  // We cannot use Assert inside a CUDA kernel
848 #ifndef __CUDA_ARCH__
849  Assert(dim != 0,
850  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
851 #endif
852  return value;
853 }
854 
855 
856 template <int dim, typename Number>
857 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
858  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator const Number &() const
859 {
860  // We cannot use Assert inside a CUDA kernel
861 #ifndef __CUDA_ARCH__
862  Assert(dim != 0,
863  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
864 #endif
865  return value;
866 }
867 
868 
869 template <int dim, typename Number>
870 template <typename OtherNumber>
871 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
872  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
874 {
875  value = internal::NumberType<Number>::value(p);
876  return *this;
877 }
878 
879 
880 #ifdef __INTEL_COMPILER
881 template <int dim, typename Number>
882 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
883  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
885 {
886  value = p.value;
887  return *this;
888 }
889 #endif
890 
891 
892 template <int dim, typename Number>
893 template <typename OtherNumber>
894 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
895  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
896  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
897 {
898  value = internal::NumberType<Number>::value(d);
899  return *this;
900 }
901 
902 
903 template <int dim, typename Number>
904 template <typename OtherNumber>
905 DEAL_II_CONSTEXPR inline bool
907 {
908 #if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
909  Assert(!(std::is_same<Number, adouble>::value ||
910  std::is_same<OtherNumber, adouble>::value),
911  ExcMessage(
912  "The Tensor equality operator for ADOL-C taped numbers has not yet "
913  "been extended to support advanced branching."));
914 #endif
915 
916  return numbers::values_are_equal(value, p.value);
917 }
918 
919 
920 template <int dim, typename Number>
921 template <typename OtherNumber>
922 constexpr bool
924 {
925  return !((*this) == p);
926 }
927 
928 
929 template <int dim, typename Number>
930 template <typename OtherNumber>
931 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
932  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
934 {
935  value += p.value;
936  return *this;
937 }
938 
939 
940 template <int dim, typename Number>
941 template <typename OtherNumber>
942 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
943  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
945 {
946  value -= p.value;
947  return *this;
948 }
949 
950 
951 
952 namespace internal
953 {
954  namespace ComplexWorkaround
955  {
956  template <typename Number, typename OtherNumber>
957  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
958  multiply_assign_scalar(Number &val, const OtherNumber &s)
959  {
960  val *= s;
961  }
962 
963 #ifdef __CUDA_ARCH__
964  template <typename Number, typename OtherNumber>
965  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
966  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
967  {
968  printf("This function is not implemented for std::complex<Number>!\n");
969  assert(false);
970  }
971 #endif
972  } // namespace ComplexWorkaround
973 } // namespace internal
974 
975 
976 template <int dim, typename Number>
977 template <typename OtherNumber>
978 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
979  DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
980  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
981 {
982  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
983  return *this;
984 }
985 
986 
987 
988 template <int dim, typename Number>
989 template <typename OtherNumber>
990 DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
991 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
992 {
993  value /= s;
994  return *this;
995 }
996 
997 
998 template <int dim, typename Number>
999 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>
1001 {
1002  return -value;
1003 }
1004 
1005 
1006 template <int dim, typename Number>
1007 inline DEAL_II_ALWAYS_INLINE typename Tensor<0, dim, Number>::real_type
1009 {
1010  Assert(dim != 0,
1011  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1012  return numbers::NumberTraits<Number>::abs(value);
1013 }
1014 
1015 
1016 template <int dim, typename Number>
1017 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1020 {
1021  // We cannot use Assert inside a CUDA kernel
1022 #ifndef __CUDA_ARCH__
1023  Assert(dim != 0,
1024  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1025 #endif
1027 }
1028 
1029 
1030 template <int dim, typename Number>
1031 template <typename OtherNumber>
1032 inline void
1033 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1034  unsigned int & index) const
1035 {
1036  Assert(dim != 0,
1037  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1038  result[index] = value;
1039  ++index;
1040 }
1041 
1042 
1043 template <int dim, typename Number>
1044 DEAL_II_CONSTEXPR inline void
1046 {
1047  // Some auto-differentiable numbers need explicit
1048  // zero initialization.
1049  value = internal::NumberType<Number>::value(0.0);
1050 }
1051 
1052 
1053 template <int dim, typename Number>
1054 template <class Archive>
1055 inline void
1056 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1057 {
1058  ar &value;
1059 }
1060 
1061 
1062 
1063 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1064 
1065 template <int rank_, int dim, typename Number>
1066 template <typename ArrayLike, std::size_t... indices>
1067 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1068  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1069  std_cxx14::index_sequence<indices...>)
1070  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1071 {
1072  static_assert(sizeof...(indices) == dim,
1073  "dim should match the number of indices");
1074 }
1075 
1076 
1077 template <int rank_, int dim, typename Number>
1078 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1080  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1081 {}
1082 
1083 
1084 template <int rank_, int dim, typename Number>
1085 template <typename OtherNumber>
1086 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1088  const Tensor<rank_, dim, OtherNumber> &initializer)
1089  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1090 {}
1091 
1092 
1093 template <int rank_, int dim, typename Number>
1094 template <typename OtherNumber>
1095 constexpr DEAL_II_ALWAYS_INLINE
1097  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1098  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
1099 {}
1100 
1101 
1102 template <int rank_, int dim, typename Number>
1103 template <typename OtherNumber>
1104 constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
1105  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1106 {
1107  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1108 }
1109 
1110 
1111 
1112 namespace internal
1113 {
1114  namespace TensorSubscriptor
1115  {
1116  template <typename ArrayElementType, int dim>
1117  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1118  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1119  subscript(ArrayElementType * values,
1120  const unsigned int i,
1121  std::integral_constant<int, dim>)
1122  {
1123  // We cannot use Assert in a CUDA kernel
1124 #ifndef __CUDA_ARCH__
1125  Assert(i < dim, ExcIndexRange(i, 0, dim));
1126 #endif
1127  return values[i];
1128  }
1129 
1130  // The variables within this struct will be referenced in the next function.
1131  // It is a workaround that allows returning a reference to a static variable
1132  // while allowing constexpr evaluation of the function.
1133  // It has to be defined outside the function because constexpr functions
1134  // cannot define static variables
1135  template <typename ArrayElementType>
1136  struct Uninitialized
1137  {
1138  static ArrayElementType value;
1139  };
1140 
1141  template <typename Type>
1142  Type Uninitialized<Type>::value;
1143 
1144  template <typename ArrayElementType>
1145  DEAL_II_CONSTEXPR inline ArrayElementType &
1146  subscript(ArrayElementType *,
1147  const unsigned int,
1148  std::integral_constant<int, 0>)
1149  {
1150  Assert(
1151  false,
1152  ExcMessage(
1153  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1154  return Uninitialized<ArrayElementType>::value;
1155  }
1156  } // namespace TensorSubscriptor
1157 } // namespace internal
1158 
1159 
1160 template <int rank_, int dim, typename Number>
1161 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV //
1163  operator[](const unsigned int i)
1164 {
1165  return ::internal::TensorSubscriptor::subscript(
1166  values, i, std::integral_constant<int, dim>());
1167 }
1168 
1169 
1170 template <int rank_, int dim, typename Number>
1171 constexpr DEAL_II_ALWAYS_INLINE
1172  DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
1173  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1174 {
1175  return values[i];
1176 }
1177 
1178 
1179 template <int rank_, int dim, typename Number>
1180 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1182  operator[](const TableIndices<rank_> &indices) const
1183 {
1184  Assert(dim != 0,
1185  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1186 
1187  return TensorAccessors::extract<rank_>(*this, indices);
1188 }
1189 
1190 
1191 
1192 template <int rank_, int dim, typename Number>
1193 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1195 {
1196  Assert(dim != 0,
1197  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1198 
1199  return TensorAccessors::extract<rank_>(*this, indices);
1200 }
1201 
1202 
1203 
1204 template <int rank_, int dim, typename Number>
1205 inline Number *
1207 {
1208  return std::addressof(
1209  this->operator[](this->unrolled_to_component_indices(0)));
1210 }
1211 
1212 
1213 
1214 template <int rank_, int dim, typename Number>
1215 inline const Number *
1217 {
1218  return std::addressof(
1219  this->operator[](this->unrolled_to_component_indices(0)));
1220 }
1221 
1222 
1223 
1224 template <int rank_, int dim, typename Number>
1225 inline Number *
1227 {
1229 }
1230 
1231 
1232 
1233 template <int rank_, int dim, typename Number>
1234 inline const Number *
1236 {
1238 }
1239 
1240 
1241 
1242 template <int rank_, int dim, typename Number>
1243 template <typename OtherNumber>
1244 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1246 {
1247  // The following loop could be written more concisely using std::copy, but
1248  // that function is only constexpr from C++20 on.
1249  for (unsigned int i = 0; i < dim; ++i)
1250  values[i] = t.values[i];
1251  return *this;
1252 }
1253 
1254 
1255 template <int rank_, int dim, typename Number>
1256 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1258 {
1260  ExcMessage("Only assignment with zero is allowed"));
1261  (void)d;
1262 
1263  for (unsigned int i = 0; i < dim; ++i)
1264  values[i] = internal::NumberType<Number>::value(0.0);
1265  return *this;
1266 }
1267 
1268 
1269 template <int rank_, int dim, typename Number>
1270 template <typename OtherNumber>
1271 DEAL_II_CONSTEXPR inline bool
1274 {
1275  for (unsigned int i = 0; i < dim; ++i)
1276  if (values[i] != p.values[i])
1277  return false;
1278  return true;
1279 }
1280 
1281 
1282 // At some places in the library, we have Point<0> for formal reasons
1283 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1284 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1285 // in the above function that the loop end check always fails, we
1286 // implement this function here
1287 template <>
1288 template <>
1289 DEAL_II_CONSTEXPR inline bool
1291 {
1292  return true;
1293 }
1294 
1295 
1296 template <int rank_, int dim, typename Number>
1297 template <typename OtherNumber>
1298 constexpr bool
1301 {
1302  return !((*this) == p);
1303 }
1304 
1305 
1306 template <int rank_, int dim, typename Number>
1307 template <typename OtherNumber>
1308 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1309  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1312 {
1313  for (unsigned int i = 0; i < dim; ++i)
1314  values[i] += p.values[i];
1315  return *this;
1316 }
1317 
1318 
1319 template <int rank_, int dim, typename Number>
1320 template <typename OtherNumber>
1321 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1322  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1325 {
1326  for (unsigned int i = 0; i < dim; ++i)
1327  values[i] -= p.values[i];
1328  return *this;
1329 }
1330 
1331 
1332 template <int rank_, int dim, typename Number>
1333 template <typename OtherNumber>
1334 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1335  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1336  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1337 {
1338  for (unsigned int i = 0; i < dim; ++i)
1339  values[i] *= s;
1340  return *this;
1341 }
1342 
1343 
1344 namespace internal
1345 {
1346  namespace TensorImplementation
1347  {
1348  template <int rank,
1349  int dim,
1350  typename Number,
1351  typename OtherNumber,
1352  typename std::enable_if<
1353  !std::is_integral<
1354  typename ProductType<Number, OtherNumber>::type>::value &&
1355  !std::is_same<Number, Differentiation::SD::Expression>::value,
1356  int>::type = 0>
1357  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1358  division_operator(Tensor<rank, dim, Number> (&t)[dim],
1359  const OtherNumber &factor)
1360  {
1361  const Number inverse_factor = Number(1.) / factor;
1362  // recurse over the base objects
1363  for (unsigned int d = 0; d < dim; ++d)
1364  t[d] *= inverse_factor;
1365  }
1366 
1367 
1368  template <int rank,
1369  int dim,
1370  typename Number,
1371  typename OtherNumber,
1372  typename std::enable_if<
1373  std::is_integral<
1374  typename ProductType<Number, OtherNumber>::type>::value ||
1375  std::is_same<Number, Differentiation::SD::Expression>::value,
1376  int>::type = 0>
1377  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1378  division_operator(::Tensor<rank, dim, Number> (&t)[dim],
1379  const OtherNumber &factor)
1380  {
1381  // recurse over the base objects
1382  for (unsigned int d = 0; d < dim; ++d)
1383  t[d] /= factor;
1384  }
1385  } // namespace TensorImplementation
1386 } // namespace internal
1387 
1388 
1389 template <int rank_, int dim, typename Number>
1390 template <typename OtherNumber>
1391 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1392  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1393  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1394 {
1395  internal::TensorImplementation::division_operator(values, s);
1396  return *this;
1397 }
1398 
1399 
1400 template <int rank_, int dim, typename Number>
1401 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1402  DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number>
1404 {
1406 
1407  for (unsigned int i = 0; i < dim; ++i)
1408  tmp.values[i] = -values[i];
1409 
1410  return tmp;
1411 }
1412 
1413 
1414 template <int rank_, int dim, typename Number>
1417 {
1418  return std::sqrt(norm_square());
1419 }
1420 
1421 
1422 template <int rank_, int dim, typename Number>
1423 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1426 {
1427  typename numbers::NumberTraits<Number>::real_type s = internal::NumberType<
1428  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1429  for (unsigned int i = 0; i < dim; ++i)
1430  s += values[i].norm_square();
1431 
1432  return s;
1433 }
1434 
1435 
1436 template <int rank_, int dim, typename Number>
1437 template <typename OtherNumber>
1438 inline void
1439 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1440 {
1441  AssertDimension(result.size(),
1442  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1443 
1444  unsigned int index = 0;
1445  unroll_recursion(result, index);
1446 }
1447 
1448 
1449 template <int rank_, int dim, typename Number>
1450 template <typename OtherNumber>
1451 inline void
1453  unsigned int & index) const
1454 {
1455  for (unsigned int i = 0; i < dim; ++i)
1456  values[i].unroll_recursion(result, index);
1457 }
1458 
1459 
1460 template <int rank_, int dim, typename Number>
1461 DEAL_II_CONSTEXPR inline unsigned int
1463  const TableIndices<rank_> &indices)
1464 {
1465  unsigned int index = 0;
1466  for (int r = 0; r < rank_; ++r)
1467  index = index * dim + indices[r];
1468 
1469  return index;
1470 }
1471 
1472 
1473 
1474 namespace internal
1475 {
1476  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1477  // and rank=2. Make sure we don't have compiler warnings.
1478 
1479  template <int dim>
1480  inline DEAL_II_CONSTEXPR unsigned int
1481  mod(const unsigned int x)
1482  {
1483  return x % dim;
1484  }
1485 
1486  template <>
1487  inline unsigned int
1488  mod<0>(const unsigned int x)
1489  {
1490  Assert(false, ExcInternalError());
1491  return x;
1492  }
1493 
1494  template <int dim>
1495  inline DEAL_II_CONSTEXPR unsigned int
1496  div(const unsigned int x)
1497  {
1498  return x / dim;
1499  }
1500 
1501  template <>
1502  inline unsigned int
1503  div<0>(const unsigned int x)
1504  {
1505  Assert(false, ExcInternalError());
1506  return x;
1507  }
1508 
1509 } // namespace internal
1510 
1511 
1512 
1513 template <int rank_, int dim, typename Number>
1514 DEAL_II_CONSTEXPR inline TableIndices<rank_>
1516 {
1519 
1520  TableIndices<rank_> indices;
1521 
1522  unsigned int remainder = i;
1523  for (int r = rank_ - 1; r >= 0; --r)
1524  {
1525  indices[r] = internal::mod<dim>(remainder);
1526  remainder = internal::div<dim>(remainder);
1527  }
1528  Assert(remainder == 0, ExcInternalError());
1529 
1530  return indices;
1531 }
1532 
1533 
1534 template <int rank_, int dim, typename Number>
1535 DEAL_II_CONSTEXPR inline void
1537 {
1538  for (unsigned int i = 0; i < dim; ++i)
1539  values[i] = internal::NumberType<Number>::value(0.0);
1540 }
1541 
1542 
1543 template <int rank_, int dim, typename Number>
1544 constexpr std::size_t
1546 {
1547  return sizeof(Tensor<rank_, dim, Number>);
1548 }
1549 
1550 
1551 template <int rank_, int dim, typename Number>
1552 template <class Archive>
1553 inline void
1554 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1555 {
1556  ar &values;
1557 }
1558 
1559 
1560 /* ----------------- Non-member functions operating on tensors. ------------ */
1561 
1566 
1574 template <int rank_, int dim, typename Number>
1575 inline std::ostream &
1576 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1577 {
1578  for (unsigned int i = 0; i < dim; ++i)
1579  {
1580  out << p[i];
1581  if (i != dim - 1)
1582  out << ' ';
1583  }
1584 
1585  return out;
1586 }
1587 
1588 
1595 template <int dim, typename Number>
1596 inline std::ostream &
1597 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1598 {
1599  out << static_cast<const Number &>(p);
1600  return out;
1601 }
1602 
1603 
1605 
1609 
1610 
1621 template <int dim, typename Number, typename Other>
1622 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1623  typename ProductType<Other, Number>::type
1624  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1625 {
1626  return object * static_cast<const Number &>(t);
1627 }
1628 
1629 
1630 
1641 template <int dim, typename Number, typename Other>
1642 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1643  typename ProductType<Number, Other>::type
1644  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1645 {
1646  return static_cast<const Number &>(t) * object;
1647 }
1648 
1649 
1661 template <int dim, typename Number, typename OtherNumber>
1662 DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
1663  typename ProductType<Number, OtherNumber>::type
1665  const Tensor<0, dim, OtherNumber> &src2)
1666 {
1667  return static_cast<const Number &>(src1) *
1668  static_cast<const OtherNumber &>(src2);
1669 }
1670 
1671 
1679 template <int dim, typename Number, typename OtherNumber>
1680 DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
1681  Tensor<0,
1682  dim,
1683  typename ProductType<Number,
1684  typename EnableIfScalar<OtherNumber>::type>::type>
1685  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1686 {
1687  return static_cast<const Number &>(t) / factor;
1688 }
1689 
1690 
1698 template <int dim, typename Number, typename OtherNumber>
1699 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1702  const Tensor<0, dim, OtherNumber> &q)
1703 {
1704  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1705 }
1706 
1707 
1715 template <int dim, typename Number, typename OtherNumber>
1716 constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1719  const Tensor<0, dim, OtherNumber> &q)
1720 {
1721  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1722 }
1723 
1724 
1737 template <int rank, int dim, typename Number, typename OtherNumber>
1738 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1739  Tensor<rank,
1740  dim,
1741  typename ProductType<Number,
1742  typename EnableIfScalar<OtherNumber>::type>::type>
1743  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1744 {
1745  // recurse over the base objects
1747  for (unsigned int d = 0; d < dim; ++d)
1748  tt[d] = t[d] * factor;
1749  return tt;
1750 }
1751 
1752 
1765 template <int rank, int dim, typename Number, typename OtherNumber>
1766 DEAL_II_CUDA_HOST_DEV DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1767  Tensor<rank,
1768  dim,
1770  OtherNumber>::type>
1771  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1772 {
1773  // simply forward to the operator above
1774  return t * factor;
1775 }
1776 
1777 
1778 namespace internal
1779 {
1780  namespace TensorImplementation
1781  {
1782  template <int rank,
1783  int dim,
1784  typename Number,
1785  typename OtherNumber,
1786  typename std::enable_if<
1787  !std::is_integral<
1788  typename ProductType<Number, OtherNumber>::type>::value,
1789  int>::type = 0>
1790  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1792  division_operator(const Tensor<rank, dim, Number> &t,
1793  const OtherNumber & factor)
1794  {
1796  const Number inverse_factor = Number(1.) / factor;
1797  // recurse over the base objects
1798  for (unsigned int d = 0; d < dim; ++d)
1799  tt[d] = t[d] * inverse_factor;
1800  return tt;
1801  }
1802 
1803 
1804  template <int rank,
1805  int dim,
1806  typename Number,
1807  typename OtherNumber,
1808  typename std::enable_if<
1809  std::is_integral<
1810  typename ProductType<Number, OtherNumber>::type>::value,
1811  int>::type = 0>
1812  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1814  division_operator(const Tensor<rank, dim, Number> &t,
1815  const OtherNumber & factor)
1816  {
1818  // recurse over the base objects
1819  for (unsigned int d = 0; d < dim; ++d)
1820  tt[d] = t[d] / factor;
1821  return tt;
1822  }
1823  } // namespace TensorImplementation
1824 } // namespace internal
1825 
1826 
1836 template <int rank, int dim, typename Number, typename OtherNumber>
1837 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1838  Tensor<rank,
1839  dim,
1840  typename ProductType<Number,
1841  typename EnableIfScalar<OtherNumber>::type>::type>
1842  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1843 {
1844  return internal::TensorImplementation::division_operator(t, factor);
1845 }
1846 
1847 
1857 template <int rank, int dim, typename Number, typename OtherNumber>
1858 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1862 {
1864 
1865  for (unsigned int i = 0; i < dim; ++i)
1866  tmp[i] += q[i];
1867 
1868  return tmp;
1869 }
1870 
1871 
1881 template <int rank, int dim, typename Number, typename OtherNumber>
1882 DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
1886 {
1888 
1889  for (unsigned int i = 0; i < dim; ++i)
1890  tmp[i] -= q[i];
1891 
1892  return tmp;
1893 }
1894 
1901 template <int dim, typename Number, typename OtherNumber>
1902 inline DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE
1905  const Tensor<0, dim, OtherNumber> &src2)
1906 {
1908 
1909  tmp *= src2;
1910 
1911  return tmp;
1912 }
1913 
1919 template <int dim, typename Number, typename OtherNumber>
1920 inline DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE
1923  const Tensor<1, dim, OtherNumber> &src2)
1924 {
1926 
1927  for (unsigned int i = 0; i < dim; ++i)
1928  tmp[i] *= src2[i];
1929 
1930  return tmp;
1931 }
1932 
1949 template <int rank, int dim, typename Number, typename OtherNumber>
1950 inline DEAL_II_CONSTEXPR DEAL_II_ALWAYS_INLINE
1953  const Tensor<rank, dim, OtherNumber> &src2)
1954 {
1956 
1957  for (unsigned int i = 0; i < dim; ++i)
1958  tmp[i] = schur_product(src1[i], src2[i]);
1959 
1960  return tmp;
1961 }
1962 
1964 
1968 
1969 
1993 template <int rank_1,
1994  int rank_2,
1995  int dim,
1996  typename Number,
1997  typename OtherNumber>
1998 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1999  typename Tensor<rank_1 + rank_2 - 2,
2000  dim,
2001  typename ProductType<Number, OtherNumber>::type>::tensor_type
2004 {
2005  typename Tensor<rank_1 + rank_2 - 2,
2006  dim,
2007  typename ProductType<Number, OtherNumber>::type>::tensor_type
2008  result{};
2009 
2010  TensorAccessors::internal::
2011  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2012  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2013  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2014 
2015  return result;
2016 }
2017 
2018 
2048 template <int index_1,
2049  int index_2,
2050  int rank_1,
2051  int rank_2,
2052  int dim,
2053  typename Number,
2054  typename OtherNumber>
2055 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2056  typename Tensor<rank_1 + rank_2 - 2,
2057  dim,
2058  typename ProductType<Number, OtherNumber>::type>::tensor_type
2061 {
2062  Assert(0 <= index_1 && index_1 < rank_1,
2063  ExcMessage(
2064  "The specified index_1 must lie within the range [0,rank_1)"));
2065  Assert(0 <= index_2 && index_2 < rank_2,
2066  ExcMessage(
2067  "The specified index_2 must lie within the range [0,rank_2)"));
2068 
2069  using namespace TensorAccessors;
2070  using namespace TensorAccessors::internal;
2071 
2072  // Reorder index_1 to the end of src1:
2073  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
2074  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2075 
2076  // Reorder index_2 to the end of src2:
2077  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2078  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2079 
2080  typename Tensor<rank_1 + rank_2 - 2,
2081  dim,
2082  typename ProductType<Number, OtherNumber>::type>::tensor_type
2083  result{};
2084  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2085  return result;
2086 }
2087 
2088 
2120 template <int index_1,
2121  int index_2,
2122  int index_3,
2123  int index_4,
2124  int rank_1,
2125  int rank_2,
2126  int dim,
2127  typename Number,
2128  typename OtherNumber>
2129 DEAL_II_CONSTEXPR inline
2130  typename Tensor<rank_1 + rank_2 - 4,
2131  dim,
2132  typename ProductType<Number, OtherNumber>::type>::tensor_type
2135 {
2136  Assert(0 <= index_1 && index_1 < rank_1,
2137  ExcMessage(
2138  "The specified index_1 must lie within the range [0,rank_1)"));
2139  Assert(0 <= index_3 && index_3 < rank_1,
2140  ExcMessage(
2141  "The specified index_3 must lie within the range [0,rank_1)"));
2142  Assert(index_1 != index_3,
2143  ExcMessage("index_1 and index_3 must not be the same"));
2144  Assert(0 <= index_2 && index_2 < rank_2,
2145  ExcMessage(
2146  "The specified index_2 must lie within the range [0,rank_2)"));
2147  Assert(0 <= index_4 && index_4 < rank_2,
2148  ExcMessage(
2149  "The specified index_4 must lie within the range [0,rank_2)"));
2150  Assert(index_2 != index_4,
2151  ExcMessage("index_2 and index_4 must not be the same"));
2152 
2153  using namespace TensorAccessors;
2154  using namespace TensorAccessors::internal;
2155 
2156  // Reorder index_1 to the end of src1:
2157  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
2158  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2159 
2160  // Reorder index_2 to the end of src2:
2161  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2162  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2163 
2164  // Now, reorder index_3 to the end of src1. We have to make sure to
2165  // preserve the original ordering: index_1 has been removed. If
2166  // index_3 > index_1, we have to use (index_3 - 1) instead:
2167  ReorderedIndexView<
2168  (index_3 < index_1 ? index_3 : index_3 - 1),
2169  rank_1,
2170  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2171  reord_3 =
2172  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2173  index_3 - 1,
2174  rank_1 > (reord_1);
2175 
2176  // Now, reorder index_4 to the end of src2. We have to make sure to
2177  // preserve the original ordering: index_2 has been removed. If
2178  // index_4 > index_2, we have to use (index_4 - 1) instead:
2179  ReorderedIndexView<
2180  (index_4 < index_2 ? index_4 : index_4 - 1),
2181  rank_2,
2182  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2183  reord_4 =
2184  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2185  index_4 - 1,
2186  rank_2 > (reord_2);
2187 
2188  typename Tensor<rank_1 + rank_2 - 4,
2189  dim,
2190  typename ProductType<Number, OtherNumber>::type>::tensor_type
2191  result{};
2192  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2193  return result;
2194 }
2195 
2196 
2210 template <int rank, int dim, typename Number, typename OtherNumber>
2211 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2212  typename ProductType<Number, OtherNumber>::type
2214  const Tensor<rank, dim, OtherNumber> &right)
2215 {
2216  typename ProductType<Number, OtherNumber>::type result{};
2217  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2218  return result;
2219 }
2220 
2221 
2240 template <template <int, int, typename> class TensorT1,
2241  template <int, int, typename> class TensorT2,
2242  template <int, int, typename> class TensorT3,
2243  int rank_1,
2244  int rank_2,
2245  int dim,
2246  typename T1,
2247  typename T2,
2248  typename T3>
2249 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2251  contract3(const TensorT1<rank_1, dim, T1> & left,
2252  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2253  const TensorT3<rank_2, dim, T3> & right)
2254 {
2255  using return_type =
2257  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2258  middle,
2259  right);
2260 }
2261 
2262 
2274 template <int rank_1,
2275  int rank_2,
2276  int dim,
2277  typename Number,
2278  typename OtherNumber>
2279 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2283 {
2284  typename Tensor<rank_1 + rank_2,
2285  dim,
2286  typename ProductType<Number, OtherNumber>::type>::tensor_type
2287  result{};
2288  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2289  return result;
2290 }
2291 
2292 
2294 
2298 
2299 
2311 template <int dim, typename Number>
2312 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2314 {
2315  Assert(dim == 2, ExcInternalError());
2316 
2317  Tensor<1, dim, Number> result;
2318 
2319  result[0] = src[1];
2320  result[1] = -src[0];
2321 
2322  return result;
2323 }
2324 
2325 
2336 template <int dim, typename Number1, typename Number2>
2337 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2340  const Tensor<1, dim, Number2> &src2)
2341 {
2342  Assert(dim == 3, ExcInternalError());
2343 
2345 
2346  result[0] = src1[1] * src2[2] - src1[2] * src2[1];
2347  result[1] = src1[2] * src2[0] - src1[0] * src2[2];
2348  result[2] = src1[0] * src2[1] - src1[1] * src2[0];
2349 
2350  return result;
2351 }
2352 
2353 
2355 
2359 
2360 
2367 template <int dim, typename Number>
2368 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2370 {
2371  // Compute the determinant using the Laplace expansion of the
2372  // determinant. We expand along the last row.
2373  Number det = internal::NumberType<Number>::value(0.0);
2374 
2375  for (unsigned int k = 0; k < dim; ++k)
2376  {
2377  Tensor<2, dim - 1, Number> minor;
2378  for (unsigned int i = 0; i < dim - 1; ++i)
2379  for (unsigned int j = 0; j < dim - 1; ++j)
2380  minor[i][j] = t[i][j < k ? j : j + 1];
2381 
2382  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2383 
2384  det += t[dim - 1][k] * cofactor;
2385  }
2386 
2387  return ((dim % 2 == 0) ? 1. : -1.) * det;
2388 }
2389 
2395 template <typename Number>
2396 constexpr DEAL_II_ALWAYS_INLINE Number
2398 {
2399  return t[0][0];
2400 }
2401 
2402 
2410 template <int dim, typename Number>
2411 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2413 {
2414  Number t = d[0][0];
2415  for (unsigned int i = 1; i < dim; ++i)
2416  t += d[i][i];
2417  return t;
2418 }
2419 
2420 
2430 template <int dim, typename Number>
2431 DEAL_II_CONSTEXPR inline Tensor<2, dim, Number>
2433 {
2434  Number return_tensor[dim][dim];
2435 
2436  // if desired, take over the
2437  // inversion of a 4x4 tensor
2438  // from the FullMatrix
2439  AssertThrow(false, ExcNotImplemented());
2440 
2441  return Tensor<2, dim, Number>(return_tensor);
2442 }
2443 
2444 
2445 #ifndef DOXYGEN
2446 
2447 template <typename Number>
2448 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, 1, Number>
2449  invert(const Tensor<2, 1, Number> &t)
2450 {
2451  Tensor<2, 1, Number> return_tensor;
2452 
2453  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2454 
2455  return return_tensor;
2456 }
2457 
2458 
2459 template <typename Number>
2460 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, 2, Number>
2461  invert(const Tensor<2, 2, Number> &t)
2462 {
2463  Tensor<2, 2, Number> return_tensor;
2464 
2465  // this is Maple output,
2466  // thus a bit unstructured
2467  const Number inv_det_t = internal::NumberType<Number>::value(
2468  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2469  return_tensor[0][0] = t[1][1];
2470  return_tensor[0][1] = -t[0][1];
2471  return_tensor[1][0] = -t[1][0];
2472  return_tensor[1][1] = t[0][0];
2473  return_tensor *= inv_det_t;
2474 
2475  return return_tensor;
2476 }
2477 
2478 
2479 template <typename Number>
2480 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, 3, Number>
2481  invert(const Tensor<2, 3, Number> &t)
2482 {
2483  Tensor<2, 3, Number> return_tensor;
2484 
2485  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2486  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2487  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2488  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2489  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2490  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2491  inv_det_t = internal::NumberType<Number>::value(
2492  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2493  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2494  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2495  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2496  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2497  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2498  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2499  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2500  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2501  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2502  return_tensor[1][1] =
2503  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2504  return_tensor[1][2] = t00 - t6;
2505  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2506  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2507  return_tensor[2][1] =
2508  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2509  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2510  return_tensor *= inv_det_t;
2511 
2512  return return_tensor;
2513 }
2514 
2515 #endif /* DOXYGEN */
2516 
2517 
2524 template <int dim, typename Number>
2525 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
2527 {
2529  for (unsigned int i = 0; i < dim; ++i)
2530  {
2531  tt[i][i] = t[i][i];
2532  for (unsigned int j = i + 1; j < dim; ++j)
2533  {
2534  tt[i][j] = t[j][i];
2535  tt[j][i] = t[i][j];
2536  };
2537  }
2538  return tt;
2539 }
2540 
2541 
2556 template <int dim, typename Number>
2557 constexpr Tensor<2, dim, Number>
2559 {
2560  return determinant(t) * invert(t);
2561 }
2562 
2563 
2578 template <int dim, typename Number>
2579 constexpr Tensor<2, dim, Number>
2581 {
2582  return transpose(adjugate(t));
2583 }
2584 
2585 
2593 template <int dim, typename Number>
2594 inline Number
2596 {
2597  Number max = internal::NumberType<Number>::value(0.0);
2598  for (unsigned int j = 0; j < dim; ++j)
2599  {
2600  Number sum = internal::NumberType<Number>::value(0.0);
2601  for (unsigned int i = 0; i < dim; ++i)
2602  sum += std::fabs(t[i][j]);
2603 
2604  if (sum > max)
2605  max = sum;
2606  }
2607 
2608  return max;
2609 }
2610 
2611 
2619 template <int dim, typename Number>
2620 inline Number
2622 {
2623  Number max = internal::NumberType<Number>::value(0.0);
2624  for (unsigned int i = 0; i < dim; ++i)
2625  {
2626  Number sum = internal::NumberType<Number>::value(0.0);
2627  for (unsigned int j = 0; j < dim; ++j)
2628  sum += std::fabs(t[i][j]);
2629 
2630  if (sum > max)
2631  max = sum;
2632  }
2633 
2634  return max;
2635 }
2636 
2638 
2639 
2640 #ifndef DOXYGEN
2641 
2642 
2643 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2644 
2645 // Specialization of functions for ADOL-C number types when
2646 // the advanced branching feature is used
2647 template <int dim>
2648 inline adouble
2650 {
2651  adouble max = internal::NumberType<adouble>::value(0.0);
2652  for (unsigned int j = 0; j < dim; ++j)
2653  {
2654  adouble sum = internal::NumberType<adouble>::value(0.0);
2655  for (unsigned int i = 0; i < dim; ++i)
2656  sum += std::fabs(t[i][j]);
2657 
2658  condassign(max, (sum > max), sum, max);
2659  }
2660 
2661  return max;
2662 }
2663 
2664 
2665 template <int dim>
2666 inline adouble
2668 {
2669  adouble max = internal::NumberType<adouble>::value(0.0);
2670  for (unsigned int i = 0; i < dim; ++i)
2671  {
2672  adouble sum = internal::NumberType<adouble>::value(0.0);
2673  for (unsigned int j = 0; j < dim; ++j)
2674  sum += std::fabs(t[i][j]);
2675 
2676  condassign(max, (sum > max), sum, max);
2677  }
2678 
2679  return max;
2680 }
2681 
2682 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2683 
2684 
2685 #endif // DOXYGEN
2686 
2687 DEAL_II_NAMESPACE_CLOSE
2688 
2689 #endif
constexpr Tensor()=default
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< rank, dim, Number > &src1, const Tensor< rank, dim, OtherNumber > &src2)
Definition: tensor.h:1952
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
DEAL_II_CONSTEXPR ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition: tensor.h:1644
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
DEAL_II_CONSTEXPR ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2213
DEAL_II_CONSTEXPR ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1624
DEAL_II_CONSTEXPR Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2526
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2580
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:943
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type contract(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2059
DEAL_II_CONSTEXPR Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2432
DEAL_II_CONSTEXPR Tensor & operator/=(const OtherNumber &factor)
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:605
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2595
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< typename EnableIfScalar< Number >::type, OtherNumber >::type > operator*(const Number &factor, const Tensor< rank, dim, OtherNumber > &t)
Definition: tensor.h:1771
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2621
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:718
DEAL_II_CONSTEXPR Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:459
DEAL_II_CONSTEXPR Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2339
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
static real_type abs(const number &x)
Definition: numbers.h:627
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1416
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1452
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Definition: point.h:110
LinearAlgebra::distributed::Vector< Number > Vector
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1860
DEAL_II_CONSTEXPR void clear()
Definition: tensor.h:1536
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1884
DEAL_II_CONSTEXPR value_type & operator[](const unsigned int i)
Definition: tensor.h:1163
DEAL_II_CONSTEXPR Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2369
static ::ExceptionBase & ExcMessage(std::string arg1)
DEAL_II_CONSTEXPR Tensor & operator*=(const OtherNumber &factor)
DEAL_II_CONSTEXPR Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2313
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1718
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1842
DEAL_II_CONSTEXPR Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
DEAL_II_CONSTEXPR ProductType< T1, typename ProductType< T2, T3 >::type >::type contract3(const TensorT1< rank_1, dim, T1 > &left, const TensorT2< rank_1+rank_2, dim, T2 > &middle, const TensorT3< rank_2, dim, T3 > &right)
Definition: tensor.h:2251
#define Assert(cond, exc)
Definition: exceptions.h:1411
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1554
static DEAL_II_CONSTEXPR unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1462
Number * end_raw()
Definition: tensor.h:1226
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:452
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type operator*(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2002
static constexpr unsigned int rank
Definition: tensor.h:438
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2558
constexpr ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1664
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2, dim, typename ProductType< Number, OtherNumber >::type > outer_product(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2281
DEAL_II_CONSTEXPR Tensor operator-() const
Definition: tensor.h:1403
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1300
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:132
DEAL_II_CONSTEXPR Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1904
DEAL_II_CONSTEXPR numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1425
constexpr Tensor< 0, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< 0, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1685
Definition: tensor.h:422
Expression operator==(const Expression &lhs, const Expression &rhs)
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:1701
DEAL_II_CONSTEXPR Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:959
Number * begin_raw()
Definition: tensor.h:1206
DEAL_II_CONSTEXPR Tensor< rank_1+rank_2 - 4, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type double_contract(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2133
Expression operator-(Expression lhs, const Expression &rhs)
static ::ExceptionBase & ExcNotImplemented()
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1439
static DEAL_II_CONSTEXPR TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1515
static constexpr unsigned int n_independent_components
Definition: tensor.h:444
DEAL_II_CONSTEXPR bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1273
Expression operator!=(const Expression &lhs, const Expression &rhs)
DEAL_II_CONSTEXPR Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2412
static constexpr std::size_t memory_consumption()
Definition: tensor.h:1545
constexpr Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2397
DEAL_II_CONSTEXPR Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, OtherNumber > &src2)
Definition: tensor.h:1922
static ::ExceptionBase & ExcInternalError()
DEAL_II_CONSTEXPR Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1743