Reference documentation for deal.II version 9.1.1
\(\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/table_indices.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/tensor_accessors.h>
26 #include <deal.II/base/utilities.h>
27 
28 #ifdef DEAL_II_WITH_ADOLC
29 # include <adolc/adouble.h> // Taped double
30 #endif
31 
32 #include <cmath>
33 #include <ostream>
34 #include <vector>
35 
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 // Forward declarations:
40 
41 template <int dim, typename Number>
42 class Point;
43 template <int rank_, int dim, typename Number = double>
44 class Tensor;
45 template <typename Number>
46 class Vector;
47 template <typename Number>
48 class VectorizedArray;
49 
50 #ifndef DOXYGEN
51 // Overload invalid tensor types of negative rank that come up during
52 // overload resolution of operator* and related contraction variants.
53 template <int dim, typename Number>
54 class Tensor<-2, dim, Number>
55 {};
56 
57 template <int dim, typename Number>
58 class Tensor<-1, dim, Number>
59 {};
60 #endif /* DOXYGEN */
61 
62 
93 template <int dim, typename Number>
94 class Tensor<0, dim, Number>
95 {
96 public:
105  static constexpr unsigned int dimension = dim;
106 
110  static constexpr unsigned int rank = 0;
111 
115  static constexpr unsigned int n_independent_components = 1;
116 
126 
131  using value_type = Number;
132 
138  using array_type = Number;
139 
145  DEAL_II_CUDA_HOST_DEV
146  Tensor();
147 
153  template <typename OtherNumber>
154  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
155 
159  template <typename OtherNumber>
160  Tensor(const OtherNumber &initializer);
161 
165  Number *
166  begin_raw();
167 
171  constexpr const Number *
172  begin_raw() const;
173 
177  Number *
178  end_raw();
179 
184  constexpr const Number *
185  end_raw() const;
186 
196  DEAL_II_CUDA_HOST_DEV operator Number &();
197 
206  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
207 
213  template <typename OtherNumber>
214  Tensor &
216 
217 #ifdef __INTEL_COMPILER
218 
224  Tensor &
225  operator=(const Tensor<0, dim, Number> &rhs);
226 #endif
227 
232  template <typename OtherNumber>
233  Tensor &
234  operator=(const OtherNumber &d);
235 
239  template <typename OtherNumber>
240  bool
241  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
242 
246  template <typename OtherNumber>
247  constexpr bool
248  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
249 
253  template <typename OtherNumber>
254  Tensor &
256 
260  template <typename OtherNumber>
261  Tensor &
263 
269  template <typename OtherNumber>
270  DEAL_II_CUDA_HOST_DEV Tensor &
271  operator*=(const OtherNumber &factor);
272 
276  template <typename OtherNumber>
277  Tensor &
278  operator/=(const OtherNumber &factor);
279 
283  constexpr Tensor
284  operator-() const;
285 
298  void
299  clear();
300 
306  real_type
307  norm() const;
308 
315  DEAL_II_CUDA_HOST_DEV real_type
316  norm_square() const;
317 
322  template <class Archive>
323  void
324  serialize(Archive &ar, const unsigned int version);
325 
330  using tensor_type = Number;
331 
332 private:
336  Number value;
337 
341  template <typename OtherNumber>
342  void
343  unroll_recursion(Vector<OtherNumber> &result,
344  unsigned int & start_index) const;
345 
349  template <int, int, typename>
350  friend class Tensor;
351 };
352 
353 
354 
396 template <int rank_, int dim, typename Number>
397 class Tensor
398 {
399 public:
408  static constexpr unsigned int dimension = dim;
409 
413  static constexpr unsigned int rank = rank_;
414 
419  static constexpr unsigned int n_independent_components =
420  Tensor<rank_ - 1, dim>::n_independent_components * dim;
421 
427  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
428 
433  using array_type =
434  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
435 
441  constexpr DEAL_II_CUDA_HOST_DEV
442  Tensor() = default;
443 
447  explicit Tensor(const array_type &initializer);
448 
454  template <typename OtherNumber>
455  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
456 
460  template <typename OtherNumber>
461  Tensor(
462  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
463 
467  template <typename OtherNumber>
468  constexpr
469  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
470 
476  DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);
477 
483  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
484  operator[](const unsigned int i) const;
485 
489  const Number &operator[](const TableIndices<rank_> &indices) const;
490 
494  Number &operator[](const TableIndices<rank_> &indices);
495 
499  Number *
500  begin_raw();
501 
505  constexpr const Number *
506  begin_raw() const;
507 
511  Number *
512  end_raw();
513 
517  constexpr const Number *
518  end_raw() const;
519 
525  template <typename OtherNumber>
526  Tensor &
528 
535  Tensor &
536  operator=(const Number &d);
537 
541  template <typename OtherNumber>
542  bool
544 
548  template <typename OtherNumber>
549  constexpr bool
551 
555  template <typename OtherNumber>
556  Tensor &
558 
562  template <typename OtherNumber>
563  Tensor &
565 
572  template <typename OtherNumber>
573  DEAL_II_CUDA_HOST_DEV Tensor &
574  operator*=(const OtherNumber &factor);
575 
579  template <typename OtherNumber>
580  Tensor &
581  operator/=(const OtherNumber &factor);
582 
586  Tensor
587  operator-() const;
588 
601  void
602  clear();
603 
611  norm() const;
612 
619  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
620  norm_square() const;
621 
629  template <typename OtherNumber>
630  void
631  unroll(Vector<OtherNumber> &result) const;
632 
637  static unsigned int
639 
644  static TableIndices<rank_>
645  unrolled_to_component_indices(const unsigned int i);
646 
651  static constexpr std::size_t
653 
658  template <class Archive>
659  void
660  serialize(Archive &ar, const unsigned int version);
661 
667 
668 private:
672  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
673  // ... avoid a compiler warning in case of dim == 0 and ensure that the
674  // array always has positive size.
675 
679  template <typename OtherNumber>
680  void
681  unroll_recursion(Vector<OtherNumber> &result,
682  unsigned int & start_index) const;
683 
687  template <int, int, typename>
688  friend class Tensor;
689 
694  friend class Point<dim, Number>;
695 };
696 
697 
698 namespace internal
699 {
713  template <int rank, int dim, typename T>
714  struct NumberType<Tensor<rank, dim, T>>
715  {
716  static constexpr const Tensor<rank, dim, T> &
717  value(const Tensor<rank, dim, T> &t)
718  {
719  return t;
720  }
721 
722  static Tensor<rank, dim, T>
723  value(const T &t)
724  {
726  tmp = t;
727  return tmp;
728  }
729  };
730 
731  template <int rank, int dim, typename T>
732  struct NumberType<Tensor<rank, dim, VectorizedArray<T>>>
733  {
734  static constexpr const Tensor<rank, dim, VectorizedArray<T>> &
735  value(const Tensor<rank, dim, VectorizedArray<T>> &t)
736  {
737  return t;
738  }
739 
741  value(const T &t)
742  {
745  return tmp;
746  }
747 
749  value(const VectorizedArray<T> &t)
750  {
752  tmp = t;
753  return tmp;
754  }
755  };
756 } // namespace internal
757 
758 
759 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
760 
761 
762 template <int dim, typename Number>
763 DEAL_II_CUDA_HOST_DEV inline Tensor<0, dim, Number>::Tensor()
764  // Some auto-differentiable numbers need explicit
765  // zero initialization.
766  : value(internal::NumberType<Number>::value(0.0))
767 {}
768 
769 
770 
771 template <int dim, typename Number>
772 template <typename OtherNumber>
773 inline Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
774 {
775  value = internal::NumberType<Number>::value(initializer);
776 }
777 
778 
779 
780 template <int dim, typename Number>
781 template <typename OtherNumber>
783 {
784  value = p.value;
785 }
786 
787 
788 
789 template <int dim, typename Number>
790 inline Number *
792 {
793  return std::addressof(value);
794 }
795 
796 
797 
798 template <int dim, typename Number>
799 constexpr const Number *
801 {
802  return std::addressof(value);
803 }
804 
805 
806 
807 template <int dim, typename Number>
808 inline Number *
810 {
812 }
813 
814 
815 
816 template <int dim, typename Number>
817 constexpr const Number *
819 {
821 }
822 
823 
824 
825 template <int dim, typename Number>
826 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
827 {
828  // We cannot use Assert inside a CUDA kernel
829 #ifndef __CUDA_ARCH__
830  Assert(dim != 0,
831  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
832 #endif
833  return value;
834 }
835 
836 
837 template <int dim, typename Number>
838 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
839  operator const Number &() const
840 {
841  // We cannot use Assert inside a CUDA kernel
842 #ifndef __CUDA_ARCH__
843  Assert(dim != 0,
844  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
845 #endif
846  return value;
847 }
848 
849 
850 template <int dim, typename Number>
851 template <typename OtherNumber>
852 inline Tensor<0, dim, Number> &
854 {
856  return *this;
857 }
858 
859 
860 #ifdef __INTEL_COMPILER
861 template <int dim, typename Number>
862 inline Tensor<0, dim, Number> &
864 {
865  value = p.value;
866  return *this;
867 }
868 #endif
869 
870 
871 template <int dim, typename Number>
872 template <typename OtherNumber>
873 inline Tensor<0, dim, Number> &
874 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
875 {
877  return *this;
878 }
879 
880 
881 template <int dim, typename Number>
882 template <typename OtherNumber>
883 inline bool
885 {
886 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
887  Assert(!(std::is_same<Number, adouble>::value ||
888  std::is_same<OtherNumber, adouble>::value),
889  ExcMessage(
890  "The Tensor equality operator for ADOL-C taped numbers has not yet "
891  "been extended to support advanced branching."));
892 #endif
893 
894  return numbers::values_are_equal(value, p.value);
895 }
896 
897 
898 template <int dim, typename Number>
899 template <typename OtherNumber>
900 constexpr bool
902 {
903  return !((*this) == p);
904 }
905 
906 
907 template <int dim, typename Number>
908 template <typename OtherNumber>
909 inline Tensor<0, dim, Number> &
911 {
912  value += p.value;
913  return *this;
914 }
915 
916 
917 template <int dim, typename Number>
918 template <typename OtherNumber>
919 inline Tensor<0, dim, Number> &
921 {
922  value -= p.value;
923  return *this;
924 }
925 
926 
927 
928 namespace internal
929 {
930  namespace ComplexWorkaround
931  {
932  template <typename Number, typename OtherNumber>
933  inline DEAL_II_CUDA_HOST_DEV void
934  multiply_assign_scalar(Number &val, const OtherNumber &s)
935  {
936  val *= s;
937  }
938 
939 #ifdef __CUDA_ARCH__
940  template <typename Number, typename OtherNumber>
941  inline DEAL_II_CUDA_HOST_DEV void
942  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
943  {
944  printf("This function is not implemented for std::complex<Number>!\n");
945  assert(false);
946  }
947 #endif
948  } // namespace ComplexWorkaround
949 } // namespace internal
950 
951 
952 template <int dim, typename Number>
953 template <typename OtherNumber>
954 inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
955 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
956 {
957  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
958  return *this;
959 }
960 
961 
962 
963 template <int dim, typename Number>
964 template <typename OtherNumber>
965 inline Tensor<0, dim, Number> &
966 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
967 {
968  value /= s;
969  return *this;
970 }
971 
972 
973 template <int dim, typename Number>
974 constexpr Tensor<0, dim, Number>
976 {
977  return -value;
978 }
979 
980 
981 template <int dim, typename Number>
982 inline typename Tensor<0, dim, Number>::real_type
984 {
985  Assert(dim != 0,
986  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
988 }
989 
990 
991 template <int dim, typename Number>
992 inline typename Tensor<0, dim, Number>::real_type DEAL_II_CUDA_HOST_DEV
994 {
995  // We cannot use Assert inside a CUDA kernel
996 #ifndef __CUDA_ARCH__
997  Assert(dim != 0,
998  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
999 #endif
1001 }
1002 
1003 
1004 template <int dim, typename Number>
1005 template <typename OtherNumber>
1006 inline void
1007 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1008  unsigned int & index) const
1009 {
1010  Assert(dim != 0,
1011  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1012  result[index] = value;
1013  ++index;
1014 }
1015 
1016 
1017 template <int dim, typename Number>
1018 inline void
1020 {
1021  // Some auto-differentiable numbers need explicit
1022  // zero initialization.
1024 }
1025 
1026 
1027 template <int dim, typename Number>
1028 template <class Archive>
1029 inline void
1030 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1031 {
1032  ar &value;
1033 }
1034 
1035 
1036 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1037 
1038 
1039 template <int rank_, int dim, typename Number>
1040 inline DEAL_II_ALWAYS_INLINE
1042 {
1043  for (unsigned int i = 0; i < dim; ++i)
1044  values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
1045 }
1046 
1047 
1048 template <int rank_, int dim, typename Number>
1049 template <typename OtherNumber>
1050 inline DEAL_II_ALWAYS_INLINE
1052  const Tensor<rank_, dim, OtherNumber> &initializer)
1053 {
1054  for (unsigned int i = 0; i != dim; ++i)
1055  values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
1056 }
1057 
1058 
1059 template <int rank_, int dim, typename Number>
1060 template <typename OtherNumber>
1061 inline DEAL_II_ALWAYS_INLINE
1063  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1064 {
1065  for (unsigned int i = 0; i < dim; ++i)
1066  values[i] = initializer[i];
1067 }
1068 
1069 
1070 template <int rank_, int dim, typename Number>
1071 template <typename OtherNumber>
1072 constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
1073  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1074 {
1075  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1076 }
1077 
1078 
1079 
1080 namespace internal
1081 {
1082  namespace TensorSubscriptor
1083  {
1084  template <typename ArrayElementType, int dim>
1085  inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV ArrayElementType &
1086  subscript(ArrayElementType * values,
1087  const unsigned int i,
1088  std::integral_constant<int, dim>)
1089  {
1090  // We cannot use Assert in a CUDA kernel
1091 #ifndef __CUDA_ARCH__
1092  Assert(i < dim, ExcIndexRange(i, 0, dim));
1093 #endif
1094  return values[i];
1095  }
1096 
1097 
1098  template <typename ArrayElementType>
1099  ArrayElementType &
1100  subscript(ArrayElementType *,
1101  const unsigned int,
1102  std::integral_constant<int, 0>)
1103  {
1104  Assert(
1105  false,
1106  ExcMessage(
1107  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1108  static ArrayElementType t;
1109  return t;
1110  }
1111  } // namespace TensorSubscriptor
1112 } // namespace internal
1113 
1114 
1115 template <int rank_, int dim, typename Number>
1116 inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
1118  operator[](const unsigned int i)
1119 {
1120  return ::internal::TensorSubscriptor::subscript(
1121  values, i, std::integral_constant<int, dim>());
1122 }
1123 
1124 
1125 template <int rank_, int dim, typename Number>
1126 constexpr DEAL_II_ALWAYS_INLINE
1127  DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
1128  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1129 {
1130  return ::internal::TensorSubscriptor::subscript(
1131  values, i, std::integral_constant<int, dim>());
1132 }
1133 
1134 
1135 template <int rank_, int dim, typename Number>
1136 inline const Number &Tensor<rank_, dim, Number>::
1137  operator[](const TableIndices<rank_> &indices) const
1138 {
1139  Assert(dim != 0,
1140  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1141 
1142  return TensorAccessors::extract<rank_>(*this, indices);
1143 }
1144 
1145 
1146 
1147 template <int rank_, int dim, typename Number>
1148 inline Number &Tensor<rank_, dim, Number>::
1150 {
1151  Assert(dim != 0,
1152  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1153 
1154  return TensorAccessors::extract<rank_>(*this, indices);
1155 }
1156 
1157 
1158 
1159 template <int rank_, int dim, typename Number>
1160 inline Number *
1162 {
1163  return std::addressof(
1164  this->operator[](this->unrolled_to_component_indices(0)));
1165 }
1166 
1167 
1168 
1169 template <int rank_, int dim, typename Number>
1170 constexpr const Number *
1172 {
1173  return std::addressof(
1174  this->operator[](this->unrolled_to_component_indices(0)));
1175 }
1176 
1177 
1178 
1179 template <int rank_, int dim, typename Number>
1180 inline Number *
1182 {
1183  return begin_raw() + n_independent_components;
1184 }
1185 
1186 
1187 
1188 template <int rank_, int dim, typename Number>
1189 constexpr const Number *
1191 {
1192  return begin_raw() + n_independent_components;
1193 }
1194 
1195 
1196 
1197 template <int rank_, int dim, typename Number>
1198 template <typename OtherNumber>
1199 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1201 {
1202  if (dim > 0)
1203  std::copy(&t.values[0], &t.values[0] + dim, &values[0]);
1204  return *this;
1205 }
1206 
1207 
1208 template <int rank_, int dim, typename Number>
1209 inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
1211 {
1213  ExcMessage("Only assignment with zero is allowed"));
1214  (void)d;
1215 
1216  for (unsigned int i = 0; i < dim; ++i)
1217  values[i] = internal::NumberType<Number>::value(0.0);
1218  return *this;
1219 }
1220 
1221 
1222 template <int rank_, int dim, typename Number>
1223 template <typename OtherNumber>
1224 inline bool
1227 {
1228  for (unsigned int i = 0; i < dim; ++i)
1229  if (values[i] != p.values[i])
1230  return false;
1231  return true;
1232 }
1233 
1234 
1235 // At some places in the library, we have Point<0> for formal reasons
1236 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1237 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1238 // in the above function that the loop end check always fails, we
1239 // implement this function here
1240 template <>
1241 template <>
1242 inline bool
1244 {
1245  return true;
1246 }
1247 
1248 
1249 template <int rank_, int dim, typename Number>
1250 template <typename OtherNumber>
1251 constexpr bool
1254 {
1255  return !((*this) == p);
1256 }
1257 
1258 
1259 template <int rank_, int dim, typename Number>
1260 template <typename OtherNumber>
1263 {
1264  for (unsigned int i = 0; i < dim; ++i)
1265  values[i] += p.values[i];
1266  return *this;
1267 }
1268 
1269 
1270 template <int rank_, int dim, typename Number>
1271 template <typename OtherNumber>
1274 {
1275  for (unsigned int i = 0; i < dim; ++i)
1276  values[i] -= p.values[i];
1277  return *this;
1278 }
1279 
1280 
1281 template <int rank_, int dim, typename Number>
1282 template <typename OtherNumber>
1283 inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
1284 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1285 {
1286  for (unsigned int i = 0; i < dim; ++i)
1287  values[i] *= s;
1288  return *this;
1289 }
1290 
1291 
1292 template <int rank_, int dim, typename Number>
1293 template <typename OtherNumber>
1295 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1296 {
1297  for (unsigned int i = 0; i < dim; ++i)
1298  values[i] /= s;
1299  return *this;
1300 }
1301 
1302 
1303 template <int rank_, int dim, typename Number>
1306 {
1308 
1309  for (unsigned int i = 0; i < dim; ++i)
1310  tmp.values[i] = -values[i];
1311 
1312  return tmp;
1313 }
1314 
1315 
1316 template <int rank_, int dim, typename Number>
1319 {
1320  return std::sqrt(norm_square());
1321 }
1322 
1323 
1324 template <int rank_, int dim, typename Number>
1325 inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
1327 {
1329  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1330  for (unsigned int i = 0; i < dim; ++i)
1331  s += values[i].norm_square();
1332 
1333  return s;
1334 }
1335 
1336 
1337 template <int rank_, int dim, typename Number>
1338 template <typename OtherNumber>
1339 inline void
1340 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1341 {
1342  AssertDimension(result.size(),
1343  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1344 
1345  unsigned int index = 0;
1346  unroll_recursion(result, index);
1347 }
1348 
1349 
1350 template <int rank_, int dim, typename Number>
1351 template <typename OtherNumber>
1352 inline void
1354  unsigned int & index) const
1355 {
1356  for (unsigned int i = 0; i < dim; ++i)
1357  values[i].unroll_recursion(result, index);
1358 }
1359 
1360 
1361 template <int rank_, int dim, typename Number>
1362 inline unsigned int
1364  const TableIndices<rank_> &indices)
1365 {
1366  unsigned int index = 0;
1367  for (int r = 0; r < rank_; ++r)
1368  index = index * dim + indices[r];
1369 
1370  return index;
1371 }
1372 
1373 
1374 
1375 namespace internal
1376 {
1377  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1378  // and rank=2. Make sure we don't have compiler warnings.
1379 
1380  template <int dim>
1381  inline unsigned int
1382  mod(const unsigned int x)
1383  {
1384  return x % dim;
1385  }
1386 
1387  template <>
1388  inline unsigned int
1389  mod<0>(const unsigned int x)
1390  {
1391  Assert(false, ExcInternalError());
1392  return x;
1393  }
1394 
1395  template <int dim>
1396  inline unsigned int
1397  div(const unsigned int x)
1398  {
1399  return x / dim;
1400  }
1401 
1402  template <>
1403  inline unsigned int
1404  div<0>(const unsigned int x)
1405  {
1406  Assert(false, ExcInternalError());
1407  return x;
1408  }
1409 
1410 } // namespace internal
1411 
1412 
1413 
1414 template <int rank_, int dim, typename Number>
1415 inline TableIndices<rank_>
1417 {
1418  Assert(i < n_independent_components,
1419  ExcIndexRange(i, 0, n_independent_components));
1420 
1421  TableIndices<rank_> indices;
1422 
1423  unsigned int remainder = i;
1424  for (int r = rank_ - 1; r >= 0; --r)
1425  {
1426  indices[r] = internal::mod<dim>(remainder);
1427  remainder = internal::div<dim>(remainder);
1428  }
1429  Assert(remainder == 0, ExcInternalError());
1430 
1431  return indices;
1432 }
1433 
1434 
1435 template <int rank_, int dim, typename Number>
1436 inline void
1438 {
1439  for (unsigned int i = 0; i < dim; ++i)
1440  values[i] = internal::NumberType<Number>::value(0.0);
1441 }
1442 
1443 
1444 template <int rank_, int dim, typename Number>
1445 constexpr std::size_t
1447 {
1448  return sizeof(Tensor<rank_, dim, Number>);
1449 }
1450 
1451 
1452 template <int rank_, int dim, typename Number>
1453 template <class Archive>
1454 inline void
1455 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1456 {
1457  ar &values;
1458 }
1459 
1460 
1461 /* ----------------- Non-member functions operating on tensors. ------------ */
1462 
1467 
1475 template <int rank_, int dim, typename Number>
1476 inline std::ostream &
1477 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1478 {
1479  for (unsigned int i = 0; i < dim; ++i)
1480  {
1481  out << p[i];
1482  if (i != dim - 1)
1483  out << ' ';
1484  }
1485 
1486  return out;
1487 }
1488 
1489 
1496 template <int dim, typename Number>
1497 inline std::ostream &
1498 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1499 {
1500  out << static_cast<const Number &>(p);
1501  return out;
1502 }
1503 
1504 
1506 
1510 
1511 
1520 template <int dim, typename Number, typename Other>
1521 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Other, Number>::type
1522 operator*(const Other &object, const Tensor<0, dim, Number> &t)
1523 {
1524  return object * static_cast<const Number &>(t);
1525 }
1526 
1527 
1528 
1537 template <int dim, typename Number, typename Other>
1538 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, Other>::type
1539 operator*(const Tensor<0, dim, Number> &t, const Other &object)
1540 {
1541  return static_cast<const Number &>(t) * object;
1542 }
1543 
1544 
1554 template <int dim, typename Number, typename OtherNumber>
1555 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
1557  const Tensor<0, dim, OtherNumber> &src2)
1558 {
1559  return static_cast<const Number &>(src1) *
1560  static_cast<const OtherNumber &>(src2);
1561 }
1562 
1563 
1569 template <int dim, typename Number, typename OtherNumber>
1570 constexpr DEAL_II_ALWAYS_INLINE
1571  Tensor<0,
1572  dim,
1573  typename ProductType<Number,
1574  typename EnableIfScalar<OtherNumber>::type>::type>
1575  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1576 {
1577  return static_cast<const Number &>(t) / factor;
1578 }
1579 
1580 
1586 template <int dim, typename Number, typename OtherNumber>
1587 constexpr DEAL_II_ALWAYS_INLINE
1590  const Tensor<0, dim, OtherNumber> &q)
1591 {
1592  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1593 }
1594 
1595 
1601 template <int dim, typename Number, typename OtherNumber>
1602 constexpr DEAL_II_ALWAYS_INLINE
1605  const Tensor<0, dim, OtherNumber> &q)
1606 {
1607  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1608 }
1609 
1610 
1621 template <int rank, int dim, typename Number, typename OtherNumber>
1622 inline DEAL_II_ALWAYS_INLINE
1623  Tensor<rank,
1624  dim,
1625  typename ProductType<Number,
1626  typename EnableIfScalar<OtherNumber>::type>::type>
1627  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1628 {
1629  // recurse over the base objects
1631  for (unsigned int d = 0; d < dim; ++d)
1632  tt[d] = t[d] * factor;
1633  return tt;
1634 }
1635 
1636 
1647 template <int rank, int dim, typename Number, typename OtherNumber>
1648 constexpr DEAL_II_ALWAYS_INLINE
1649  Tensor<rank,
1650  dim,
1652  OtherNumber>::type>
1653  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
1654 {
1655  // simply forward to the operator above
1656  return t * factor;
1657 }
1658 
1659 
1667 template <int rank, int dim, typename Number, typename OtherNumber>
1668 inline Tensor<
1669  rank,
1670  dim,
1671  typename ProductType<Number,
1672  typename EnableIfScalar<OtherNumber>::type>::type>
1673 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1674 {
1675  // recurse over the base objects
1677  for (unsigned int d = 0; d < dim; ++d)
1678  tt[d] = t[d] / factor;
1679  return tt;
1680 }
1681 
1682 
1690 template <int rank, int dim, typename Number, typename OtherNumber>
1691 inline DEAL_II_ALWAYS_INLINE
1695 {
1697 
1698  for (unsigned int i = 0; i < dim; ++i)
1699  tmp[i] += q[i];
1700 
1701  return tmp;
1702 }
1703 
1704 
1712 template <int rank, int dim, typename Number, typename OtherNumber>
1713 inline DEAL_II_ALWAYS_INLINE
1717 {
1719 
1720  for (unsigned int i = 0; i < dim; ++i)
1721  tmp[i] -= q[i];
1722 
1723  return tmp;
1724 }
1725 
1726 
1728 
1732 
1733 
1757 template <int rank_1,
1758  int rank_2,
1759  int dim,
1760  typename Number,
1761  typename OtherNumber>
1762 inline DEAL_II_ALWAYS_INLINE
1763  typename Tensor<rank_1 + rank_2 - 2,
1764  dim,
1765  typename ProductType<Number, OtherNumber>::type>::tensor_type
1768 {
1769  typename Tensor<rank_1 + rank_2 - 2,
1770  dim,
1771  typename ProductType<Number, OtherNumber>::type>::tensor_type
1772  result;
1773 
1774  TensorAccessors::internal::
1775  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1776  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
1777  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
1778 
1779  return result;
1780 }
1781 
1782 
1812 template <int index_1,
1813  int index_2,
1814  int rank_1,
1815  int rank_2,
1816  int dim,
1817  typename Number,
1818  typename OtherNumber>
1819 inline DEAL_II_ALWAYS_INLINE
1820  typename Tensor<rank_1 + rank_2 - 2,
1821  dim,
1822  typename ProductType<Number, OtherNumber>::type>::tensor_type
1825 {
1826  Assert(0 <= index_1 && index_1 < rank_1,
1827  ExcMessage(
1828  "The specified index_1 must lie within the range [0,rank_1)"));
1829  Assert(0 <= index_2 && index_2 < rank_2,
1830  ExcMessage(
1831  "The specified index_2 must lie within the range [0,rank_2)"));
1832 
1833  using namespace TensorAccessors;
1834  using namespace TensorAccessors::internal;
1835 
1836  // Reorder index_1 to the end of src1:
1837  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
1838  reord_01 = reordered_index_view<index_1, rank_1>(src1);
1839 
1840  // Reorder index_2 to the end of src2:
1841  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1842  reord_02 = reordered_index_view<index_2, rank_2>(src2);
1843 
1844  typename Tensor<rank_1 + rank_2 - 2,
1845  dim,
1846  typename ProductType<Number, OtherNumber>::type>::tensor_type
1847  result;
1848  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
1849  return result;
1850 }
1851 
1852 
1884 template <int index_1,
1885  int index_2,
1886  int index_3,
1887  int index_4,
1888  int rank_1,
1889  int rank_2,
1890  int dim,
1891  typename Number,
1892  typename OtherNumber>
1893 inline
1894  typename Tensor<rank_1 + rank_2 - 4,
1895  dim,
1896  typename ProductType<Number, OtherNumber>::type>::tensor_type
1899 {
1900  Assert(0 <= index_1 && index_1 < rank_1,
1901  ExcMessage(
1902  "The specified index_1 must lie within the range [0,rank_1)"));
1903  Assert(0 <= index_3 && index_3 < rank_1,
1904  ExcMessage(
1905  "The specified index_3 must lie within the range [0,rank_1)"));
1906  Assert(index_1 != index_3,
1907  ExcMessage("index_1 and index_3 must not be the same"));
1908  Assert(0 <= index_2 && index_2 < rank_2,
1909  ExcMessage(
1910  "The specified index_2 must lie within the range [0,rank_2)"));
1911  Assert(0 <= index_4 && index_4 < rank_2,
1912  ExcMessage(
1913  "The specified index_4 must lie within the range [0,rank_2)"));
1914  Assert(index_2 != index_4,
1915  ExcMessage("index_2 and index_4 must not be the same"));
1916 
1917  using namespace TensorAccessors;
1918  using namespace TensorAccessors::internal;
1919 
1920  // Reorder index_1 to the end of src1:
1921  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>
1922  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
1923 
1924  // Reorder index_2 to the end of src2:
1925  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>
1926  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
1927 
1928  // Now, reorder index_3 to the end of src1. We have to make sure to
1929  // preserve the original ordering: index_1 has been removed. If
1930  // index_3 > index_1, we have to use (index_3 - 1) instead:
1931  ReorderedIndexView<
1932  (index_3 < index_1 ? index_3 : index_3 - 1),
1933  rank_1,
1934  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
1935  reord_3 =
1936  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
1937  index_3 - 1,
1938  rank_1 > (reord_1);
1939 
1940  // Now, reorder index_4 to the end of src2. We have to make sure to
1941  // preserve the original ordering: index_2 has been removed. If
1942  // index_4 > index_2, we have to use (index_4 - 1) instead:
1943  ReorderedIndexView<
1944  (index_4 < index_2 ? index_4 : index_4 - 1),
1945  rank_2,
1946  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
1947  reord_4 =
1948  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
1949  index_4 - 1,
1950  rank_2 > (reord_2);
1951 
1952  typename Tensor<rank_1 + rank_2 - 4,
1953  dim,
1954  typename ProductType<Number, OtherNumber>::type>::tensor_type
1955  result;
1956  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
1957  return result;
1958 }
1959 
1960 
1974 template <int rank, int dim, typename Number, typename OtherNumber>
1975 inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
1977  const Tensor<rank, dim, OtherNumber> &right)
1978 {
1979  typename ProductType<Number, OtherNumber>::type result;
1980  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
1981  return result;
1982 }
1983 
1984 
2003 template <template <int, int, typename> class TensorT1,
2004  template <int, int, typename> class TensorT2,
2005  template <int, int, typename> class TensorT3,
2006  int rank_1,
2007  int rank_2,
2008  int dim,
2009  typename T1,
2010  typename T2,
2011  typename T3>
2013 contract3(const TensorT1<rank_1, dim, T1> & left,
2014  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2015  const TensorT3<rank_2, dim, T3> & right)
2016 {
2017  using return_type =
2019  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2020  middle,
2021  right);
2022 }
2023 
2024 
2036 template <int rank_1,
2037  int rank_2,
2038  int dim,
2039  typename Number,
2040  typename OtherNumber>
2041 inline DEAL_II_ALWAYS_INLINE
2045 {
2046  typename Tensor<rank_1 + rank_2,
2047  dim,
2048  typename ProductType<Number, OtherNumber>::type>::tensor_type
2049  result;
2050  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2051  return result;
2052 }
2053 
2054 
2056 
2060 
2061 
2073 template <int dim, typename Number>
2074 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2076 {
2077  Assert(dim == 2, ExcInternalError());
2078 
2079  Tensor<1, dim, Number> result;
2080 
2081  result[0] = src[1];
2082  result[1] = -src[0];
2083 
2084  return result;
2085 }
2086 
2087 
2098 template <int dim, typename Number>
2099 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
2101  const Tensor<1, dim, Number> &src2)
2102 {
2103  Assert(dim == 3, ExcInternalError());
2104 
2105  Tensor<1, dim, Number> result;
2106 
2107  result[0] = src1[1] * src2[2] - src1[2] * src2[1];
2108  result[1] = src1[2] * src2[0] - src1[0] * src2[2];
2109  result[2] = src1[0] * src2[1] - src1[1] * src2[0];
2110 
2111  return result;
2112 }
2113 
2114 
2116 
2120 
2121 
2128 template <int dim, typename Number>
2129 inline Number
2131 {
2132  // Compute the determinant using the Laplace expansion of the
2133  // determinant. We expand along the last row.
2134  Number det = internal::NumberType<Number>::value(0.0);
2135 
2136  for (unsigned int k = 0; k < dim; ++k)
2137  {
2138  Tensor<2, dim - 1, Number> minor;
2139  for (unsigned int i = 0; i < dim - 1; ++i)
2140  for (unsigned int j = 0; j < dim - 1; ++j)
2141  minor[i][j] = t[i][j < k ? j : j + 1];
2142 
2143  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2144 
2145  det += t[dim - 1][k] * cofactor;
2146  }
2147 
2148  return ((dim % 2 == 0) ? 1. : -1.) * det;
2149 }
2150 
2156 template <typename Number>
2157 constexpr Number
2159 {
2160  return t[0][0];
2161 }
2162 
2163 
2171 template <int dim, typename Number>
2172 inline DEAL_II_ALWAYS_INLINE Number
2174 {
2175  Number t = d[0][0];
2176  for (unsigned int i = 1; i < dim; ++i)
2177  t += d[i][i];
2178  return t;
2179 }
2180 
2181 
2191 template <int dim, typename Number>
2194 {
2195  Number return_tensor[dim][dim];
2196 
2197  // if desired, take over the
2198  // inversion of a 4x4 tensor
2199  // from the FullMatrix
2200  AssertThrow(false, ExcNotImplemented());
2201 
2202  return Tensor<2, dim, Number>(return_tensor);
2203 }
2204 
2205 
2206 #ifndef DOXYGEN
2207 
2208 template <typename Number>
2209 inline Tensor<2, 1, Number>
2210 invert(const Tensor<2, 1, Number> &t)
2211 {
2212  Number return_tensor[1][1];
2213 
2214  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2215 
2216  return Tensor<2, 1, Number>(return_tensor);
2217 }
2218 
2219 
2220 template <typename Number>
2221 inline Tensor<2, 2, Number>
2222 invert(const Tensor<2, 2, Number> &t)
2223 {
2224  Tensor<2, 2, Number> return_tensor;
2225 
2226  // this is Maple output,
2227  // thus a bit unstructured
2228  const Number inv_det_t = internal::NumberType<Number>::value(
2229  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2230  return_tensor[0][0] = t[1][1];
2231  return_tensor[0][1] = -t[0][1];
2232  return_tensor[1][0] = -t[1][0];
2233  return_tensor[1][1] = t[0][0];
2234  return_tensor *= inv_det_t;
2235 
2236  return return_tensor;
2237 }
2238 
2239 
2240 template <typename Number>
2241 inline Tensor<2, 3, Number>
2242 invert(const Tensor<2, 3, Number> &t)
2243 {
2244  Tensor<2, 3, Number> return_tensor;
2245 
2246  const Number t4 = internal::NumberType<Number>::value(t[0][0] * t[1][1]),
2247  t6 = internal::NumberType<Number>::value(t[0][0] * t[1][2]),
2248  t8 = internal::NumberType<Number>::value(t[0][1] * t[1][0]),
2249  t00 = internal::NumberType<Number>::value(t[0][2] * t[1][0]),
2250  t01 = internal::NumberType<Number>::value(t[0][1] * t[2][0]),
2251  t04 = internal::NumberType<Number>::value(t[0][2] * t[2][0]),
2253  1.0 / (t4 * t[2][2] - t6 * t[2][1] - t8 * t[2][2] +
2254  t00 * t[2][1] + t01 * t[1][2] - t04 * t[1][1]));
2255  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2256  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2257  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2258  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2259  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2260  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2261  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2262  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2263  return_tensor[1][1] =
2264  internal::NumberType<Number>::value(t[0][0] * t[2][2]) - t04;
2265  return_tensor[1][2] = t00 - t6;
2266  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2267  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2268  return_tensor[2][1] =
2269  t01 - internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2270  return_tensor[2][2] = internal::NumberType<Number>::value(t4 - t8);
2271  return_tensor *= inv_det_t;
2272 
2273  return return_tensor;
2274 }
2275 
2276 #endif /* DOXYGEN */
2277 
2278 
2285 template <int dim, typename Number>
2286 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
2288 {
2290  for (unsigned int i = 0; i < dim; ++i)
2291  {
2292  tt[i][i] = t[i][i];
2293  for (unsigned int j = i + 1; j < dim; ++j)
2294  {
2295  tt[i][j] = t[j][i];
2296  tt[j][i] = t[i][j];
2297  };
2298  }
2299  return tt;
2300 }
2301 
2302 
2317 template <int dim, typename Number>
2318 constexpr Tensor<2, dim, Number>
2320 {
2321  return determinant(t) * invert(t);
2322 }
2323 
2324 
2339 template <int dim, typename Number>
2340 constexpr Tensor<2, dim, Number>
2342 {
2343  return transpose(adjugate(t));
2344 }
2345 
2346 
2354 template <int dim, typename Number>
2355 inline Number
2357 {
2358  Number max = internal::NumberType<Number>::value(0.0);
2359  for (unsigned int j = 0; j < dim; ++j)
2360  {
2361  Number sum = internal::NumberType<Number>::value(0.0);
2362  for (unsigned int i = 0; i < dim; ++i)
2363  sum += std::fabs(t[i][j]);
2364 
2365  if (sum > max)
2366  max = sum;
2367  }
2368 
2369  return max;
2370 }
2371 
2372 
2380 template <int dim, typename Number>
2381 inline Number
2383 {
2384  Number max = internal::NumberType<Number>::value(0.0);
2385  for (unsigned int i = 0; i < dim; ++i)
2386  {
2387  Number sum = internal::NumberType<Number>::value(0.0);
2388  for (unsigned int j = 0; j < dim; ++j)
2389  sum += std::fabs(t[i][j]);
2390 
2391  if (sum > max)
2392  max = sum;
2393  }
2394 
2395  return max;
2396 }
2397 
2399 
2400 
2401 #ifndef DOXYGEN
2402 
2403 
2404 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2405 
2406 // Specialization of functions for ADOL-C number types when
2407 // the advanced branching feature is used
2408 template <int dim>
2409 inline adouble
2410 l1_norm(const Tensor<2, dim, adouble> &t)
2411 {
2412  adouble max = internal::NumberType<adouble>::value(0.0);
2413  for (unsigned int j = 0; j < dim; ++j)
2414  {
2415  adouble sum = internal::NumberType<adouble>::value(0.0);
2416  for (unsigned int i = 0; i < dim; ++i)
2417  sum += std::fabs(t[i][j]);
2418 
2419  condassign(max, (sum > max), sum, max);
2420  }
2421 
2422  return max;
2423 }
2424 
2425 
2426 template <int dim>
2427 inline adouble
2428 linfty_norm(const Tensor<2, dim, adouble> &t)
2429 {
2431  for (unsigned int i = 0; i < dim; ++i)
2432  {
2434  for (unsigned int j = 0; j < dim; ++j)
2435  sum += std::fabs(t[i][j]);
2436 
2437  condassign(max, (sum > max), sum, max);
2438  }
2439 
2440  return max;
2441 }
2442 
2443 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2444 
2445 
2446 #endif // DOXYGEN
2447 
2448 DEAL_II_NAMESPACE_CLOSE
2449 
2450 // include deprecated non-member functions operating on Tensor
2451 #include <deal.II/base/tensor_deprecated.h>
2452 
2453 #endif
constexpr Tensor()=default
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:1766
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Number determinant(const SymmetricTensor< 2, dim, Number > &)
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:1363
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2341
Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1673
bool value_is_zero(const Number &value)
Definition: numbers.h:876
bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1226
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:513
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1522
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2356
Tensor & operator/=(const OtherNumber &factor)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2382
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:1897
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:672
typename Tensor< rank_ - 1, dim, VectorizedArray< Number > >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:434
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static real_type abs(const number &x)
Definition: numbers.h:535
internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1318
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: tensor.h:1353
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Definition: point.h:110
bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:860
LinearAlgebra::distributed::Vector< Number > Vector
Tensor & operator*=(const OtherNumber &factor)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:666
Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2287
static ::ExceptionBase & ExcMessage(std::string arg1)
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:2013
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:1604
T sum(const T &t, const MPI_Comm &mpi_communicator)
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2193
Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2075
#define Assert(cond, exc)
Definition: exceptions.h:1407
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:1976
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:1653
void serialize(Archive &ar, const unsigned int version)
Definition: tensor.h:1455
Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:1627
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1326
Number * end_raw()
Definition: tensor.h:1181
typename Tensor< rank_ - 1, dim, VectorizedArray< Number > >::tensor_type value_type
Definition: tensor.h:427
value_type & operator[](const unsigned int i)
Definition: tensor.h:1118
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1693
Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
static constexpr unsigned int rank
Definition: tensor.h:413
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2319
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1416
constexpr ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:1556
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:2100
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Definition: tensor.h:1253
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:125
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:1575
Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2130
Definition: mpi.h:90
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition: tensor.h:1539
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:1715
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:1589
Number * begin_raw()
Definition: tensor.h:1161
static ::ExceptionBase & ExcNotImplemented()
static constexpr unsigned int dimension
Definition: tensor.h:408
void unroll(Vector< OtherNumber > &result) const
Definition: tensor.h:1340
Tensor operator-() const
Definition: tensor.h:1305
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:1823
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:2043
Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
static constexpr unsigned int n_independent_components
Definition: tensor.h:419
void clear()
Definition: tensor.h:1437
T max(const T &t, const MPI_Comm &mpi_communicator)
static constexpr std::size_t memory_consumption()
Definition: tensor.h:1446
constexpr Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2158
static ::ExceptionBase & ExcInternalError()