Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_tensor_h
17 #define dealii_tensor_h
18 
19 #include <deal.II/base/config.h>
20 
22 #include <deal.II/base/numbers.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 <utility>
35 #include <vector>
36 
37 
39 
40 // Forward declarations:
41 #ifndef DOXYGEN
42 template <typename ElementType, typename MemorySpace>
43 class ArrayView;
44 template <int dim, typename Number>
45 class Point;
46 template <int rank_, int dim, typename Number = double>
47 class Tensor;
48 template <typename Number>
49 class Vector;
50 template <typename number>
51 class FullMatrix;
52 namespace Differentiation
53 {
54  namespace SD
55  {
56  class Expression;
57  }
58 } // namespace Differentiation
59 #endif
60 
61 
91 template <int dim, typename Number>
92 class Tensor<0, dim, Number>
93 {
94 public:
95  static_assert(dim >= 0,
96  "Tensors must have a dimension greater than or equal to one.");
97 
106  static constexpr unsigned int dimension = dim;
107 
111  static constexpr unsigned int rank = 0;
112 
116  static constexpr unsigned int n_independent_components = 1;
117 
127 
132  using value_type = Number;
133 
139  using array_type = Number;
140 
146  constexpr DEAL_II_CUDA_HOST_DEV
147  Tensor();
148 
156  template <typename OtherNumber>
157  constexpr DEAL_II_CUDA_HOST_DEV
158  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
159 
165  template <typename OtherNumber>
166  constexpr DEAL_II_CUDA_HOST_DEV
167  Tensor(const OtherNumber &initializer);
168 
169 #if __GNUC__ >= 11 || defined __INTEL_COMPILER
170 
173  constexpr DEAL_II_CUDA_HOST_DEV
174  Tensor(const Tensor<0, dim, Number> &other);
175 
180  operator=(const Tensor<0, dim, Number> &other);
181 
185  constexpr DEAL_II_CUDA_HOST_DEV
186  Tensor(Tensor<0, dim, Number> &&other) noexcept;
187 
192  operator=(Tensor<0, dim, Number> &&other) noexcept;
193 #endif
194 
198  Number *
199  begin_raw();
200 
204  const Number *
205  begin_raw() const;
206 
210  Number *
211  end_raw();
212 
217  const Number *
218  end_raw() const;
219 
229  constexpr DEAL_II_CUDA_HOST_DEV operator Number &();
230 
239  constexpr DEAL_II_CUDA_HOST_DEV operator const Number &() const;
240 
248  template <typename OtherNumber>
249  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
250  operator=(const Tensor<0, dim, OtherNumber> &rhs);
251 
258  template <typename OtherNumber>
259  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
260  operator=(const OtherNumber &d);
261 
265  template <typename OtherNumber>
266  constexpr bool
267  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
268 
272  template <typename OtherNumber>
273  constexpr bool
274  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
275 
281  template <typename OtherNumber>
282  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
283  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
284 
290  template <typename OtherNumber>
291  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
292  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
293 
299  template <typename OtherNumber>
300  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
301  operator*=(const OtherNumber &factor);
302 
308  template <typename OtherNumber>
309  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
310  operator/=(const OtherNumber &factor);
311 
317  constexpr DEAL_II_CUDA_HOST_DEV Tensor
318  operator-() const;
319 
332  constexpr void
333  clear();
334 
340  real_type
341  norm() const;
342 
350  norm_square() const;
351 
357  template <class Archive>
358  void
359  serialize(Archive &ar, const unsigned int version);
360 
365  using tensor_type = Number;
366 
367 private:
371  Number value;
372 
376  template <typename OtherNumber>
377  void
378  unroll_recursion(Vector<OtherNumber> &result,
379  unsigned int & start_index) const;
380 
381  // Allow an arbitrary Tensor to access the underlying values.
382  template <int, int, typename>
383  friend class Tensor;
384 };
385 
386 
387 
461 template <int rank_, int dim, typename Number>
462 class Tensor
463 {
464 public:
465  static_assert(rank_ >= 1,
466  "Tensors must have a rank greater than or equal to one.");
467  static_assert(dim >= 0,
468  "Tensors must have a dimension greater than or equal to one.");
477  static constexpr unsigned int dimension = dim;
478 
482  static constexpr unsigned int rank = rank_;
483 
488  static constexpr unsigned int n_independent_components =
489  Tensor<rank_ - 1, dim>::n_independent_components * dim;
490 
496  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
497 
502  using array_type =
503  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
504 
511  Tensor();
512 
518  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
519  const array_type &initializer);
520 
534  template <typename ElementType, typename MemorySpace>
535  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
536  const ArrayView<ElementType, MemorySpace> &initializer);
537 
545  template <typename OtherNumber>
546  constexpr DEAL_II_CUDA_HOST_DEV
547  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
548 
552  template <typename OtherNumber>
553  constexpr Tensor(
554  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
555 
559  template <typename OtherNumber>
560  constexpr
561  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
562 
563 #if __GNUC__ >= 11 || defined __INTEL_COMPILER
564 
567  constexpr Tensor(const Tensor<rank_, dim, Number> &);
568 
572  constexpr Tensor<rank_, dim, Number> &
573  operator=(const Tensor<rank_, dim, Number> &);
574 
578  constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
579 
583  constexpr Tensor<rank_, dim, Number> &
584  operator=(Tensor<rank_, dim, Number> &&) noexcept;
585 #endif
586 
592  constexpr DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);
593 
599  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
600  operator[](const unsigned int i) const;
601 
605  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
606 
610  constexpr Number &operator[](const TableIndices<rank_> &indices);
611 
615  Number *
616  begin_raw();
617 
621  const Number *
622  begin_raw() const;
623 
627  Number *
628  end_raw();
629 
633  const Number *
634  end_raw() const;
635 
643  template <typename OtherNumber>
644  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
645  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
646 
653  constexpr Tensor &
654  operator=(const Number &d);
655 
659  template <typename OtherNumber>
660  constexpr bool
662 
666  template <typename OtherNumber>
667  constexpr bool
669 
675  template <typename OtherNumber>
676  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
677  operator+=(const Tensor<rank_, dim, OtherNumber> &);
678 
684  template <typename OtherNumber>
685  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
686  operator-=(const Tensor<rank_, dim, OtherNumber> &);
687 
694  template <typename OtherNumber>
695  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
696  operator*=(const OtherNumber &factor);
697 
703  template <typename OtherNumber>
704  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
705  operator/=(const OtherNumber &factor);
706 
712  constexpr DEAL_II_CUDA_HOST_DEV Tensor
713  operator-() const;
714 
727  constexpr void
728  clear();
729 
739  norm() const;
740 
747  constexpr DEAL_II_CUDA_HOST_DEV
749  norm_square() const;
750 
758  template <typename OtherNumber>
759  void
760  unroll(Vector<OtherNumber> &result) const;
761 
766  static constexpr unsigned int
767  component_to_unrolled_index(const TableIndices<rank_> &indices);
768 
774  static constexpr TableIndices<rank_>
775  unrolled_to_component_indices(const unsigned int i);
776 
781  static constexpr std::size_t
783 
789  template <class Archive>
790  void
791  serialize(Archive &ar, const unsigned int version);
792 
798 
799 private:
803  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
804  // ... avoid a compiler warning in case of dim == 0 and ensure that the
805  // array always has positive size.
806 
810  template <typename OtherNumber>
811  void
812  unroll_recursion(Vector<OtherNumber> &result,
813  unsigned int & start_index) const;
814 
821  template <typename ArrayLike, std::size_t... Indices>
822  constexpr DEAL_II_CUDA_HOST_DEV
823  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
824 
825  // Allow an arbitrary Tensor to access the underlying values.
826  template <int, int, typename>
827  friend class Tensor;
828 
829  // Point is allowed access to the coordinates. This is supposed to improve
830  // speed.
831  friend class Point<dim, Number>;
832 };
833 
834 
835 #ifndef DOXYGEN
836 namespace internal
837 {
838  // Workaround: The following 4 overloads are necessary to be able to
839  // compile the library with Apple Clang 8 and older. We should remove
840  // these overloads again when we bump the minimal required version to
841  // something later than clang-3.6 / Apple Clang 6.3.
842  template <int rank, int dim, typename T, typename U>
843  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
844  {
845  using type =
847  };
848 
849  template <int rank, int dim, typename T, typename U>
850  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
851  {
852  using type =
854  };
855 
856  template <typename T, int rank, int dim, typename U>
857  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
858  {
859  using type =
861  };
862 
863  template <int rank, int dim, typename T, typename U>
864  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
865  {
866  using type =
868  };
869  // end workaround
870 
875  template <int rank, int dim, typename T>
876  struct NumberType<Tensor<rank, dim, T>>
877  {
878  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
879  value(const Tensor<rank, dim, T> &t)
880  {
881  return t;
882  }
883 
885  value(const T &t)
886  {
888  tmp = t;
889  return tmp;
890  }
891  };
892 } // namespace internal
893 
894 
895 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
896 
897 
898 template <int dim, typename Number>
901  // Some auto-differentiable numbers need explicit
902  // zero initialization such as adtl::adouble.
903  : Tensor{0.0}
904 {}
905 
906 
907 
908 template <int dim, typename Number>
909 template <typename OtherNumber>
911  Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
912  : value(internal::NumberType<Number>::value(initializer))
913 {}
914 
915 
916 
917 template <int dim, typename Number>
918 template <typename OtherNumber>
921  : Tensor{p.value}
922 {}
923 
924 
925 
926 # if __GNUC__ >= 11 || defined __INTEL_COMPILER
927 template <int dim, typename Number>
930  : value{other.value}
931 {}
932 
933 
934 
935 template <int dim, typename Number>
938 {
939  value = other.value;
940  return *this;
941 }
942 
943 
944 
945 template <int dim, typename Number>
948  : value{std::move(other.value)}
949 {}
950 
951 
952 
953 template <int dim, typename Number>
956 {
957  value = std::move(other.value);
958  return *this;
959 }
960 # endif
961 
962 
963 
964 template <int dim, typename Number>
965 inline Number *
967 {
968  return std::addressof(value);
969 }
970 
971 
972 
973 template <int dim, typename Number>
974 inline const Number *
976 {
977  return std::addressof(value);
978 }
979 
980 
981 
982 template <int dim, typename Number>
983 inline Number *
985 {
987 }
988 
989 
990 
991 template <int dim, typename Number>
992 const Number *
994 {
996 }
997 
998 
999 
1000 template <int dim, typename Number>
1001 constexpr inline DEAL_II_ALWAYS_INLINE
1003 {
1004  // We cannot use Assert inside a CUDA kernel
1005 # ifndef __CUDA_ARCH__
1006  Assert(dim != 0,
1007  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1008 # endif
1009  return value;
1010 }
1011 
1012 
1013 template <int dim, typename Number>
1014 constexpr inline DEAL_II_ALWAYS_INLINE
1016 {
1017  // We cannot use Assert inside a CUDA kernel
1018 # ifndef __CUDA_ARCH__
1019  Assert(dim != 0,
1020  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1021 # endif
1022  return value;
1023 }
1024 
1025 
1026 template <int dim, typename Number>
1027 template <typename OtherNumber>
1028 constexpr inline DEAL_II_ALWAYS_INLINE
1031 {
1033  return *this;
1034 }
1035 
1036 
1037 template <int dim, typename Number>
1038 template <typename OtherNumber>
1039 constexpr inline DEAL_II_ALWAYS_INLINE
1041  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
1042 {
1044  return *this;
1045 }
1046 
1047 
1048 template <int dim, typename Number>
1049 template <typename OtherNumber>
1050 constexpr inline bool
1052 {
1053 # if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING)
1054  Assert(!(std::is_same<Number, adouble>::value ||
1055  std::is_same<OtherNumber, adouble>::value),
1056  ExcMessage(
1057  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1058  "been extended to support advanced branching."));
1059 # endif
1060 
1061  return numbers::values_are_equal(value, p.value);
1062 }
1063 
1064 
1065 template <int dim, typename Number>
1066 template <typename OtherNumber>
1067 constexpr bool
1069 {
1070  return !((*this) == p);
1071 }
1072 
1073 
1074 template <int dim, typename Number>
1075 template <typename OtherNumber>
1076 constexpr inline DEAL_II_ALWAYS_INLINE
1079 {
1080  value += p.value;
1081  return *this;
1082 }
1083 
1084 
1085 template <int dim, typename Number>
1086 template <typename OtherNumber>
1087 constexpr inline DEAL_II_ALWAYS_INLINE
1090 {
1091  value -= p.value;
1092  return *this;
1093 }
1094 
1095 
1096 
1097 namespace internal
1098 {
1099  namespace ComplexWorkaround
1100  {
1101  template <typename Number, typename OtherNumber>
1102  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1103  multiply_assign_scalar(Number &val, const OtherNumber &s)
1104  {
1105  val *= s;
1106  }
1107 
1108 # ifdef __CUDA_ARCH__
1109  template <typename Number, typename OtherNumber>
1110  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1111  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1112  {
1113  printf("This function is not implemented for std::complex<Number>!\n");
1114  assert(false);
1115  }
1116 # endif
1117  } // namespace ComplexWorkaround
1118 } // namespace internal
1119 
1120 
1121 template <int dim, typename Number>
1122 template <typename OtherNumber>
1123 constexpr inline DEAL_II_ALWAYS_INLINE
1125  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1126 {
1127  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1128  return *this;
1129 }
1130 
1131 
1132 
1133 template <int dim, typename Number>
1134 template <typename OtherNumber>
1136 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1137 {
1138  value /= s;
1139  return *this;
1140 }
1141 
1142 
1143 template <int dim, typename Number>
1146 {
1147  return -value;
1148 }
1149 
1150 
1151 template <int dim, typename Number>
1154 {
1155  Assert(dim != 0,
1156  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1157  return numbers::NumberTraits<Number>::abs(value);
1158 }
1159 
1160 
1161 template <int dim, typename Number>
1165 {
1166  // We cannot use Assert inside a CUDA kernel
1167 # ifndef __CUDA_ARCH__
1168  Assert(dim != 0,
1169  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1170 # endif
1172 }
1173 
1174 
1175 template <int dim, typename Number>
1176 template <typename OtherNumber>
1177 inline void
1178 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1179  unsigned int & index) const
1180 {
1181  Assert(dim != 0,
1182  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1183  result[index] = value;
1184  ++index;
1185 }
1186 
1187 
1188 template <int dim, typename Number>
1189 constexpr inline void
1191 {
1192  // Some auto-differentiable numbers need explicit
1193  // zero initialization.
1195 }
1196 
1197 
1198 template <int dim, typename Number>
1199 template <class Archive>
1200 inline void
1201 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1202 {
1203  ar &value;
1204 }
1205 
1206 
1207 template <int dim, typename Number>
1209 
1210 
1211 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1212 
1213 template <int rank_, int dim, typename Number>
1214 template <typename ArrayLike, std::size_t... indices>
1216  Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1217  std::index_sequence<indices...>)
1218  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1219 {
1220  static_assert(sizeof...(indices) == dim,
1221  "dim should match the number of indices");
1222 }
1223 
1224 
1225 
1226 template <int rank_, int dim, typename Number>
1229  // We would like to use =default, but this causes compile errors with some
1230  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1231  : values{}
1232 {}
1233 
1234 
1235 
1236 template <int rank_, int dim, typename Number>
1239  : Tensor(initializer, std::make_index_sequence<dim>{})
1240 {}
1241 
1242 
1243 
1244 template <int rank_, int dim, typename Number>
1245 template <typename ElementType, typename MemorySpace>
1248  const ArrayView<ElementType, MemorySpace> &initializer)
1249 {
1251 
1252  for (unsigned int i = 0; i < n_independent_components; ++i)
1253  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1254 }
1255 
1256 
1257 
1258 template <int rank_, int dim, typename Number>
1259 template <typename OtherNumber>
1262  const Tensor<rank_, dim, OtherNumber> &initializer)
1263  : Tensor(initializer, std::make_index_sequence<dim>{})
1264 {}
1265 
1266 
1267 template <int rank_, int dim, typename Number>
1268 template <typename OtherNumber>
1269 constexpr DEAL_II_ALWAYS_INLINE
1271  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1272  : Tensor(initializer, std::make_index_sequence<dim>{})
1273 {}
1274 
1275 
1276 template <int rank_, int dim, typename Number>
1277 template <typename OtherNumber>
1279  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1280 {
1281  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1282 }
1283 
1284 
1285 # if __GNUC__ >= 11 || defined __INTEL_COMPILER
1286 template <int rank_, int dim, typename Number>
1287 constexpr DEAL_II_ALWAYS_INLINE
1289 {
1290  for (unsigned int i = 0; i < dim; ++i)
1291  values[i] = other.values[i];
1292 }
1293 
1294 
1295 template <int rank_, int dim, typename Number>
1298 {
1299  for (unsigned int i = 0; i < dim; ++i)
1300  values[i] = other.values[i];
1301  return *this;
1302 }
1303 
1304 
1305 template <int rank_, int dim, typename Number>
1306 constexpr DEAL_II_ALWAYS_INLINE
1308 {
1309  for (unsigned int i = 0; i < dim; ++i)
1310  values[i] = other.values[i];
1311 }
1312 
1313 
1314 template <int rank_, int dim, typename Number>
1317  operator=(Tensor<rank_, dim, Number> &&other) noexcept
1318 {
1319  for (unsigned int i = 0; i < dim; ++i)
1320  values[i] = other.values[i];
1321  return *this;
1322 }
1323 # endif
1324 
1325 
1326 namespace internal
1327 {
1328  namespace TensorSubscriptor
1329  {
1330  template <typename ArrayElementType, int dim>
1331  constexpr inline DEAL_II_ALWAYS_INLINE
1332  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1333  subscript(ArrayElementType * values,
1334  const unsigned int i,
1335  std::integral_constant<int, dim>)
1336  {
1337  // We cannot use Assert in a CUDA kernel
1338 # ifndef __CUDA_ARCH__
1339  AssertIndexRange(i, dim);
1340 # endif
1341  return values[i];
1342  }
1343 
1344  // The variables within this struct will be referenced in the next function.
1345  // It is a workaround that allows returning a reference to a static variable
1346  // while allowing constexpr evaluation of the function.
1347  // It has to be defined outside the function because constexpr functions
1348  // cannot define static variables
1349  template <typename ArrayElementType>
1350  struct Uninitialized
1351  {
1352  static ArrayElementType value;
1353  };
1354 
1355  template <typename Type>
1356  Type Uninitialized<Type>::value;
1357 
1358  template <typename ArrayElementType>
1359  constexpr inline DEAL_II_ALWAYS_INLINE
1360  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1361  subscript(ArrayElementType *,
1362  const unsigned int,
1363  std::integral_constant<int, 0>)
1364  {
1365  // We cannot use Assert in a CUDA kernel
1366 # ifndef __CUDA_ARCH__
1367  Assert(
1368  false,
1369  ExcMessage(
1370  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1371 # endif
1372  return Uninitialized<ArrayElementType>::value;
1373  }
1374  } // namespace TensorSubscriptor
1375 } // namespace internal
1376 
1377 
1378 template <int rank_, int dim, typename Number>
1381  operator[](const unsigned int i)
1382 {
1383  return ::internal::TensorSubscriptor::subscript(
1384  values, i, std::integral_constant<int, dim>());
1385 }
1386 
1387 
1388 template <int rank_, int dim, typename Number>
1389 constexpr DEAL_II_ALWAYS_INLINE
1391  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1392 {
1393 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1394  AssertIndexRange(i, dim);
1395 # endif
1396 
1397  return values[i];
1398 }
1399 
1400 
1401 template <int rank_, int dim, typename Number>
1402 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1404  operator[](const TableIndices<rank_> &indices) const
1405 {
1406 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1407  Assert(dim != 0,
1408  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1409 # endif
1410 
1411  return TensorAccessors::extract<rank_>(*this, indices);
1412 }
1413 
1414 
1415 
1416 template <int rank_, int dim, typename Number>
1417 constexpr inline DEAL_II_ALWAYS_INLINE Number &Tensor<rank_, dim, Number>::
1418  operator[](const TableIndices<rank_> &indices)
1419 {
1420 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1421  Assert(dim != 0,
1422  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1423 # endif
1424 
1425  return TensorAccessors::extract<rank_>(*this, indices);
1426 }
1427 
1428 
1429 
1430 template <int rank_, int dim, typename Number>
1431 inline Number *
1433 {
1434  return std::addressof(
1435  this->operator[](this->unrolled_to_component_indices(0)));
1436 }
1437 
1438 
1439 
1440 template <int rank_, int dim, typename Number>
1441 inline const Number *
1443 {
1444  return std::addressof(
1445  this->operator[](this->unrolled_to_component_indices(0)));
1446 }
1447 
1448 
1449 
1450 template <int rank_, int dim, typename Number>
1451 inline Number *
1453 {
1455 }
1456 
1457 
1458 
1459 template <int rank_, int dim, typename Number>
1460 inline const Number *
1462 {
1464 }
1465 
1466 
1467 
1468 template <int rank_, int dim, typename Number>
1469 template <typename OtherNumber>
1472 {
1473  // The following loop could be written more concisely using std::copy, but
1474  // that function is only constexpr from C++20 on.
1475  for (unsigned int i = 0; i < dim; ++i)
1476  values[i] = t.values[i];
1477  return *this;
1478 }
1479 
1480 
1481 template <int rank_, int dim, typename Number>
1484 {
1486  (void)d;
1487 
1488  for (unsigned int i = 0; i < dim; ++i)
1490  return *this;
1491 }
1492 
1493 
1494 template <int rank_, int dim, typename Number>
1495 template <typename OtherNumber>
1496 constexpr inline bool
1499 {
1500  for (unsigned int i = 0; i < dim; ++i)
1501  if (values[i] != p.values[i])
1502  return false;
1503  return true;
1504 }
1505 
1506 
1507 // At some places in the library, we have Point<0> for formal reasons
1508 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1509 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1510 // in the above function that the loop end check always fails, we
1511 // implement this function here
1512 template <>
1513 template <>
1514 constexpr inline bool
1516 {
1517  return true;
1518 }
1519 
1520 
1521 template <int rank_, int dim, typename Number>
1522 template <typename OtherNumber>
1523 constexpr bool
1526 {
1527  return !((*this) == p);
1528 }
1529 
1530 
1531 template <int rank_, int dim, typename Number>
1532 template <typename OtherNumber>
1533 constexpr inline DEAL_II_ALWAYS_INLINE
1537 {
1538  for (unsigned int i = 0; i < dim; ++i)
1539  values[i] += p.values[i];
1540  return *this;
1541 }
1542 
1543 
1544 template <int rank_, int dim, typename Number>
1545 template <typename OtherNumber>
1546 constexpr inline DEAL_II_ALWAYS_INLINE
1550 {
1551  for (unsigned int i = 0; i < dim; ++i)
1552  values[i] -= p.values[i];
1553  return *this;
1554 }
1555 
1556 
1557 template <int rank_, int dim, typename Number>
1558 template <typename OtherNumber>
1559 constexpr inline DEAL_II_ALWAYS_INLINE
1561  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1562 {
1563  for (unsigned int i = 0; i < dim; ++i)
1564  values[i] *= s;
1565  return *this;
1566 }
1567 
1568 
1569 namespace internal
1570 {
1571  namespace TensorImplementation
1572  {
1573  template <int rank,
1574  int dim,
1575  typename Number,
1576  typename OtherNumber,
1577  typename std::enable_if<
1578  !std::is_integral<
1579  typename ProductType<Number, OtherNumber>::type>::value &&
1580  !std::is_same<Number, Differentiation::SD::Expression>::value,
1581  int>::type = 0>
1582  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1584  const OtherNumber &factor)
1585  {
1586  const Number inverse_factor = Number(1.) / factor;
1587  // recurse over the base objects
1588  for (unsigned int d = 0; d < dim; ++d)
1589  t[d] *= inverse_factor;
1590  }
1591 
1592 
1593  template <int rank,
1594  int dim,
1595  typename Number,
1596  typename OtherNumber,
1597  typename std::enable_if<
1598  std::is_integral<
1599  typename ProductType<Number, OtherNumber>::type>::value ||
1600  std::is_same<Number, Differentiation::SD::Expression>::value,
1601  int>::type = 0>
1602  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1604  const OtherNumber &factor)
1605  {
1606  // recurse over the base objects
1607  for (unsigned int d = 0; d < dim; ++d)
1608  t[d] /= factor;
1609  }
1610  } // namespace TensorImplementation
1611 } // namespace internal
1612 
1613 
1614 template <int rank_, int dim, typename Number>
1615 template <typename OtherNumber>
1616 constexpr inline DEAL_II_ALWAYS_INLINE
1618  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1619 {
1621  return *this;
1622 }
1623 
1624 
1625 template <int rank_, int dim, typename Number>
1626 constexpr inline DEAL_II_ALWAYS_INLINE
1629 {
1631 
1632  for (unsigned int i = 0; i < dim; ++i)
1633  tmp.values[i] = -values[i];
1634 
1635  return tmp;
1636 }
1637 
1638 
1639 template <int rank_, int dim, typename Number>
1642 {
1643  return std::sqrt(norm_square());
1644 }
1645 
1646 
1647 template <int rank_, int dim, typename Number>
1651 {
1653  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1654  for (unsigned int i = 0; i < dim; ++i)
1655  s += values[i].norm_square();
1656 
1657  return s;
1658 }
1659 
1660 
1661 template <int rank_, int dim, typename Number>
1662 template <typename OtherNumber>
1663 inline void
1664 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1665 {
1666  AssertDimension(result.size(),
1667  (Utilities::fixed_power<rank_, unsigned int>(dim)));
1668 
1669  unsigned int index = 0;
1670  unroll_recursion(result, index);
1671 }
1672 
1673 
1674 template <int rank_, int dim, typename Number>
1675 template <typename OtherNumber>
1676 inline void
1677 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
1678  unsigned int & index) const
1679 {
1680  for (unsigned int i = 0; i < dim; ++i)
1681  values[i].unroll_recursion(result, index);
1682 }
1683 
1684 
1685 template <int rank_, int dim, typename Number>
1686 constexpr inline unsigned int
1688  const TableIndices<rank_> &indices)
1689 {
1690  unsigned int index = 0;
1691  for (int r = 0; r < rank_; ++r)
1692  index = index * dim + indices[r];
1693 
1694  return index;
1695 }
1696 
1697 
1698 
1699 namespace internal
1700 {
1701  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1702  // and rank=2. Make sure we don't have compiler warnings.
1703 
1704  template <int dim>
1705  inline constexpr unsigned int
1706  mod(const unsigned int x)
1707  {
1708  return x % dim;
1709  }
1710 
1711  template <>
1712  inline unsigned int
1713  mod<0>(const unsigned int x)
1714  {
1715  Assert(false, ExcInternalError());
1716  return x;
1717  }
1718 
1719  template <int dim>
1720  inline constexpr unsigned int
1721  div(const unsigned int x)
1722  {
1723  return x / dim;
1724  }
1725 
1726  template <>
1727  inline unsigned int
1728  div<0>(const unsigned int x)
1729  {
1730  Assert(false, ExcInternalError());
1731  return x;
1732  }
1733 
1734 } // namespace internal
1735 
1736 
1737 
1738 template <int rank_, int dim, typename Number>
1739 constexpr inline TableIndices<rank_>
1741 {
1742  AssertIndexRange(i, n_independent_components);
1743 
1744  TableIndices<rank_> indices;
1745 
1746  unsigned int remainder = i;
1747  for (int r = rank_ - 1; r >= 0; --r)
1748  {
1749  indices[r] = internal::mod<dim>(remainder);
1750  remainder = internal::div<dim>(remainder);
1751  }
1752  Assert(remainder == 0, ExcInternalError());
1753 
1754  return indices;
1755 }
1756 
1757 
1758 template <int rank_, int dim, typename Number>
1759 constexpr inline void
1761 {
1762  for (unsigned int i = 0; i < dim; ++i)
1764 }
1765 
1766 
1767 template <int rank_, int dim, typename Number>
1768 constexpr std::size_t
1770 {
1771  return sizeof(Tensor<rank_, dim, Number>);
1772 }
1773 
1774 
1775 template <int rank_, int dim, typename Number>
1776 template <class Archive>
1777 inline void
1778 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1779 {
1780  ar &values;
1781 }
1782 
1783 
1784 template <int rank_, int dim, typename Number>
1786 
1787 #endif // DOXYGEN
1788 
1789 /* ----------------- Non-member functions operating on tensors. ------------ */
1790 
1795 
1803 template <int rank_, int dim, typename Number>
1804 inline std::ostream &
1805 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1806 {
1807  for (unsigned int i = 0; i < dim; ++i)
1808  {
1809  out << p[i];
1810  if (i != dim - 1)
1811  out << ' ';
1812  }
1813 
1814  return out;
1815 }
1816 
1817 
1824 template <int dim, typename Number>
1825 inline std::ostream &
1826 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1827 {
1828  out << static_cast<const Number &>(p);
1829  return out;
1830 }
1831 
1832 
1834 
1838 
1839 
1850 template <int dim, typename Number, typename Other>
1853  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1854 {
1855  return object * static_cast<const Number &>(t);
1856 }
1857 
1858 
1859 
1870 template <int dim, typename Number, typename Other>
1873  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1874 {
1875  return static_cast<const Number &>(t) * object;
1876 }
1877 
1878 
1890 template <int dim, typename Number, typename OtherNumber>
1894  const Tensor<0, dim, OtherNumber> &src2)
1895 {
1896  return static_cast<const Number &>(src1) *
1897  static_cast<const OtherNumber &>(src2);
1898 }
1899 
1900 
1908 template <int dim, typename Number, typename OtherNumber>
1910  Tensor<0,
1911  dim,
1912  typename ProductType<Number,
1913  typename EnableIfScalar<OtherNumber>::type>::type>
1914  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
1915 {
1916  return static_cast<const Number &>(t) / factor;
1917 }
1918 
1919 
1927 template <int dim, typename Number, typename OtherNumber>
1931  const Tensor<0, dim, OtherNumber> &q)
1932 {
1933  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
1934 }
1935 
1936 
1944 template <int dim, typename Number, typename OtherNumber>
1948  const Tensor<0, dim, OtherNumber> &q)
1949 {
1950  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
1951 }
1952 
1953 
1966 template <int rank, int dim, typename Number, typename OtherNumber>
1968  Tensor<rank,
1969  dim,
1970  typename ProductType<Number,
1971  typename EnableIfScalar<OtherNumber>::type>::type>
1972  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
1973 {
1974  // recurse over the base objects
1976  for (unsigned int d = 0; d < dim; ++d)
1977  tt[d] = t[d] * factor;
1978  return tt;
1979 }
1980 
1981 
1994 template <int rank, int dim, typename Number, typename OtherNumber>
1996  Tensor<rank,
1997  dim,
1999  OtherNumber>::type>
2000  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2001 {
2002  // simply forward to the operator above
2003  return t * factor;
2004 }
2005 
2006 
2007 namespace internal
2008 {
2009  namespace TensorImplementation
2010  {
2011  template <int rank,
2012  int dim,
2013  typename Number,
2014  typename OtherNumber,
2015  typename std::enable_if<
2016  !std::is_integral<
2017  typename ProductType<Number, OtherNumber>::type>::value,
2018  int>::type = 0>
2022  const OtherNumber & factor)
2023  {
2025  const Number inverse_factor = Number(1.) / factor;
2026  // recurse over the base objects
2027  for (unsigned int d = 0; d < dim; ++d)
2028  tt[d] = t[d] * inverse_factor;
2029  return tt;
2030  }
2031 
2032 
2033  template <int rank,
2034  int dim,
2035  typename Number,
2036  typename OtherNumber,
2037  typename std::enable_if<
2038  std::is_integral<
2039  typename ProductType<Number, OtherNumber>::type>::value,
2040  int>::type = 0>
2044  const OtherNumber & factor)
2045  {
2047  // recurse over the base objects
2048  for (unsigned int d = 0; d < dim; ++d)
2049  tt[d] = t[d] / factor;
2050  return tt;
2051  }
2052  } // namespace TensorImplementation
2053 } // namespace internal
2054 
2055 
2065 template <int rank, int dim, typename Number, typename OtherNumber>
2067  Tensor<rank,
2068  dim,
2069  typename ProductType<Number,
2070  typename EnableIfScalar<OtherNumber>::type>::type>
2071  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2072 {
2074 }
2075 
2076 
2086 template <int rank, int dim, typename Number, typename OtherNumber>
2091 {
2093 
2094  for (unsigned int i = 0; i < dim; ++i)
2095  tmp[i] += q[i];
2096 
2097  return tmp;
2098 }
2099 
2100 
2110 template <int rank, int dim, typename Number, typename OtherNumber>
2115 {
2117 
2118  for (unsigned int i = 0; i < dim; ++i)
2119  tmp[i] -= q[i];
2120 
2121  return tmp;
2122 }
2123 
2130 template <int dim, typename Number, typename OtherNumber>
2131 inline constexpr DEAL_II_ALWAYS_INLINE
2134  const Tensor<0, dim, OtherNumber> &src2)
2135 {
2137 
2138  tmp *= src2;
2139 
2140  return tmp;
2141 }
2142 
2159 template <int rank, int dim, typename Number, typename OtherNumber>
2160 inline constexpr DEAL_II_ALWAYS_INLINE
2163  const Tensor<rank, dim, OtherNumber> &src2)
2164 {
2166 
2167  for (unsigned int i = 0; i < dim; ++i)
2168  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2170 
2171  return tmp;
2172 }
2173 
2175 
2179 
2180 
2203 template <int rank_1,
2204  int rank_2,
2205  int dim,
2206  typename Number,
2207  typename OtherNumber,
2208  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2209 constexpr inline DEAL_II_ALWAYS_INLINE
2210  typename Tensor<rank_1 + rank_2 - 2,
2211  dim,
2212  typename ProductType<Number, OtherNumber>::type>::tensor_type
2215 {
2216  typename Tensor<rank_1 + rank_2 - 2,
2217  dim,
2218  typename ProductType<Number, OtherNumber>::type>::tensor_type
2219  result{};
2220 
2221  TensorAccessors::internal::
2222  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2223  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2224  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2225 
2226  return result;
2227 }
2228 
2229 
2258 template <int index_1,
2259  int index_2,
2260  int rank_1,
2261  int rank_2,
2262  int dim,
2263  typename Number,
2264  typename OtherNumber>
2265 constexpr inline DEAL_II_ALWAYS_INLINE
2266  typename Tensor<rank_1 + rank_2 - 2,
2267  dim,
2268  typename ProductType<Number, OtherNumber>::type>::tensor_type
2271 {
2272  Assert(0 <= index_1 && index_1 < rank_1,
2273  ExcMessage(
2274  "The specified index_1 must lie within the range [0,rank_1)"));
2275  Assert(0 <= index_2 && index_2 < rank_2,
2276  ExcMessage(
2277  "The specified index_2 must lie within the range [0,rank_2)"));
2278 
2279  using namespace TensorAccessors;
2280  using namespace TensorAccessors::internal;
2281 
2282  // Reorder index_1 to the end of src1:
2284  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2285 
2286  // Reorder index_2 to the end of src2:
2287  const ReorderedIndexView<index_2,
2288  rank_2,
2290  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2291 
2292  typename Tensor<rank_1 + rank_2 - 2,
2293  dim,
2294  typename ProductType<Number, OtherNumber>::type>::tensor_type
2295  result{};
2296  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2297  return result;
2298 }
2299 
2300 
2331 template <int index_1,
2332  int index_2,
2333  int index_3,
2334  int index_4,
2335  int rank_1,
2336  int rank_2,
2337  int dim,
2338  typename Number,
2339  typename OtherNumber>
2340 constexpr inline
2341  typename Tensor<rank_1 + rank_2 - 4,
2342  dim,
2343  typename ProductType<Number, OtherNumber>::type>::tensor_type
2346 {
2347  Assert(0 <= index_1 && index_1 < rank_1,
2348  ExcMessage(
2349  "The specified index_1 must lie within the range [0,rank_1)"));
2350  Assert(0 <= index_3 && index_3 < rank_1,
2351  ExcMessage(
2352  "The specified index_3 must lie within the range [0,rank_1)"));
2353  Assert(index_1 != index_3,
2354  ExcMessage("index_1 and index_3 must not be the same"));
2355  Assert(0 <= index_2 && index_2 < rank_2,
2356  ExcMessage(
2357  "The specified index_2 must lie within the range [0,rank_2)"));
2358  Assert(0 <= index_4 && index_4 < rank_2,
2359  ExcMessage(
2360  "The specified index_4 must lie within the range [0,rank_2)"));
2361  Assert(index_2 != index_4,
2362  ExcMessage("index_2 and index_4 must not be the same"));
2363 
2364  using namespace TensorAccessors;
2365  using namespace TensorAccessors::internal;
2366 
2367  // Reorder index_1 to the end of src1:
2369  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2370 
2371  // Reorder index_2 to the end of src2:
2373  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2374 
2375  // Now, reorder index_3 to the end of src1. We have to make sure to
2376  // preserve the original ordering: index_1 has been removed. If
2377  // index_3 > index_1, we have to use (index_3 - 1) instead:
2379  (index_3 < index_1 ? index_3 : index_3 - 1),
2380  rank_1,
2381  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2382  reord_3 =
2383  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2384  index_3 - 1,
2385  rank_1 > (reord_1);
2386 
2387  // Now, reorder index_4 to the end of src2. We have to make sure to
2388  // preserve the original ordering: index_2 has been removed. If
2389  // index_4 > index_2, we have to use (index_4 - 1) instead:
2390  ReorderedIndexView<
2391  (index_4 < index_2 ? index_4 : index_4 - 1),
2392  rank_2,
2393  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2394  reord_4 =
2395  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2396  index_4 - 1,
2397  rank_2 > (reord_2);
2398 
2399  typename Tensor<rank_1 + rank_2 - 4,
2400  dim,
2401  typename ProductType<Number, OtherNumber>::type>::tensor_type
2402  result{};
2403  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2404  return result;
2405 }
2406 
2407 
2420 template <int rank, int dim, typename Number, typename OtherNumber>
2421 constexpr inline DEAL_II_ALWAYS_INLINE
2424  const Tensor<rank, dim, OtherNumber> &right)
2425 {
2426  typename ProductType<Number, OtherNumber>::type result{};
2427  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2428  return result;
2429 }
2430 
2431 
2449 template <template <int, int, typename> class TensorT1,
2450  template <int, int, typename> class TensorT2,
2451  template <int, int, typename> class TensorT3,
2452  int rank_1,
2453  int rank_2,
2454  int dim,
2455  typename T1,
2456  typename T2,
2457  typename T3>
2458 constexpr inline DEAL_II_ALWAYS_INLINE
2460  contract3(const TensorT1<rank_1, dim, T1> & left,
2461  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2462  const TensorT3<rank_2, dim, T3> & right)
2463 {
2464  using return_type =
2466  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2467  middle,
2468  right);
2469 }
2470 
2471 
2482 template <int rank_1,
2483  int rank_2,
2484  int dim,
2485  typename Number,
2486  typename OtherNumber>
2487 constexpr inline DEAL_II_ALWAYS_INLINE
2491 {
2492  typename Tensor<rank_1 + rank_2,
2493  dim,
2494  typename ProductType<Number, OtherNumber>::type>::tensor_type
2495  result{};
2496  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2497  return result;
2498 }
2499 
2500 
2502 
2506 
2507 
2518 template <int dim, typename Number>
2521 {
2522  Assert(dim == 2, ExcInternalError());
2523 
2524  Tensor<1, dim, Number> result;
2525 
2526  result[0] = src[1];
2527  result[1] = -src[0];
2528 
2529  return result;
2530 }
2531 
2532 
2542 template <int dim, typename Number1, typename Number2>
2543 constexpr inline DEAL_II_ALWAYS_INLINE
2546  const Tensor<1, dim, Number2> &src2)
2547 {
2548  Assert(dim == 3, ExcInternalError());
2549 
2551 
2552  // avoid compiler warnings
2553  constexpr int s0 = 0 % dim;
2554  constexpr int s1 = 1 % dim;
2555  constexpr int s2 = 2 % dim;
2556 
2557  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2558  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2559  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2560 
2561  return result;
2562 }
2563 
2564 
2566 
2570 
2571 
2577 template <int dim, typename Number>
2578 constexpr inline DEAL_II_ALWAYS_INLINE Number
2580 {
2581  // Compute the determinant using the Laplace expansion of the
2582  // determinant. We expand along the last row.
2583  Number det = internal::NumberType<Number>::value(0.0);
2584 
2585  for (unsigned int k = 0; k < dim; ++k)
2586  {
2587  Tensor<2, dim - 1, Number> minor;
2588  for (unsigned int i = 0; i < dim - 1; ++i)
2589  for (unsigned int j = 0; j < dim - 1; ++j)
2590  minor[i][j] = t[i][j < k ? j : j + 1];
2591 
2592  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2593 
2594  det += t[dim - 1][k] * cofactor;
2595  }
2596 
2597  return ((dim % 2 == 0) ? 1. : -1.) * det;
2598 }
2599 
2605 template <typename Number>
2606 constexpr DEAL_II_ALWAYS_INLINE Number
2608 {
2609  return t[0][0];
2610 }
2611 
2617 template <typename Number>
2618 constexpr DEAL_II_ALWAYS_INLINE Number
2620 {
2621  // hard-coded for efficiency reasons
2622  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2623 }
2624 
2630 template <typename Number>
2631 constexpr DEAL_II_ALWAYS_INLINE Number
2633 {
2634  // hard-coded for efficiency reasons
2635  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2636  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2637  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2638  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2639  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2640  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2641  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2642 }
2643 
2644 
2651 template <int dim, typename Number>
2652 constexpr inline DEAL_II_ALWAYS_INLINE Number
2654 {
2655  Number t = d[0][0];
2656  for (unsigned int i = 1; i < dim; ++i)
2657  t += d[i][i];
2658  return t;
2659 }
2660 
2661 
2670 template <int dim, typename Number>
2671 constexpr inline Tensor<2, dim, Number>
2673 {
2674  Number return_tensor[dim][dim];
2675 
2676  // if desired, take over the
2677  // inversion of a 4x4 tensor
2678  // from the FullMatrix
2679  AssertThrow(false, ExcNotImplemented());
2680 
2681  return Tensor<2, dim, Number>(return_tensor);
2682 }
2683 
2684 
2685 #ifndef DOXYGEN
2686 
2687 template <typename Number>
2689  invert(const Tensor<2, 1, Number> &t)
2690 {
2691  Tensor<2, 1, Number> return_tensor;
2692 
2693  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2694 
2695  return return_tensor;
2696 }
2697 
2698 
2699 template <typename Number>
2701  invert(const Tensor<2, 2, Number> &t)
2702 {
2703  Tensor<2, 2, Number> return_tensor;
2704 
2705  const Number inv_det_t = internal::NumberType<Number>::value(
2706  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2707  return_tensor[0][0] = t[1][1];
2708  return_tensor[0][1] = -t[0][1];
2709  return_tensor[1][0] = -t[1][0];
2710  return_tensor[1][1] = t[0][0];
2711  return_tensor *= inv_det_t;
2712 
2713  return return_tensor;
2714 }
2715 
2716 
2717 template <typename Number>
2719  invert(const Tensor<2, 3, Number> &t)
2720 {
2721  Tensor<2, 3, Number> return_tensor;
2722 
2723  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2724  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2725  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2726  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2727  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2728  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2729  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2730  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2731  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2732  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2733  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2734  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2735  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2736  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2737  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2738  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2739  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2740  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2741  const Number inv_det_t = internal::NumberType<Number>::value(
2742  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2743  t[0][2] * return_tensor[2][0]));
2744  return_tensor *= inv_det_t;
2745 
2746  return return_tensor;
2747 }
2748 
2749 #endif /* DOXYGEN */
2750 
2751 
2757 template <int dim, typename Number>
2760 {
2762  for (unsigned int i = 0; i < dim; ++i)
2763  {
2764  tt[i][i] = t[i][i];
2765  for (unsigned int j = i + 1; j < dim; ++j)
2766  {
2767  tt[i][j] = t[j][i];
2768  tt[j][i] = t[i][j];
2769  };
2770  }
2771  return tt;
2772 }
2773 
2774 
2788 template <int dim, typename Number>
2789 constexpr Tensor<2, dim, Number>
2791 {
2792  return determinant(t) * invert(t);
2793 }
2794 
2795 
2809 template <int dim, typename Number>
2810 constexpr Tensor<2, dim, Number>
2812 {
2813  return transpose(adjugate(t));
2814 }
2815 
2816 
2880 template <int dim, typename Number>
2883 
2884 
2892 template <int dim, typename Number>
2893 inline Number
2895 {
2897  for (unsigned int j = 0; j < dim; ++j)
2898  {
2900  for (unsigned int i = 0; i < dim; ++i)
2901  sum += std::fabs(t[i][j]);
2902 
2903  if (sum > max)
2904  max = sum;
2905  }
2906 
2907  return max;
2908 }
2909 
2910 
2918 template <int dim, typename Number>
2919 inline Number
2921 {
2923  for (unsigned int i = 0; i < dim; ++i)
2924  {
2926  for (unsigned int j = 0; j < dim; ++j)
2927  sum += std::fabs(t[i][j]);
2928 
2929  if (sum > max)
2930  max = sum;
2931  }
2932 
2933  return max;
2934 }
2935 
2937 
2938 
2939 #ifndef DOXYGEN
2940 
2941 
2942 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2943 
2944 // Specialization of functions for ADOL-C number types when
2945 // the advanced branching feature is used
2946 template <int dim>
2947 inline adouble
2949 {
2951  for (unsigned int j = 0; j < dim; ++j)
2952  {
2954  for (unsigned int i = 0; i < dim; ++i)
2955  sum += std::fabs(t[i][j]);
2956 
2957  condassign(max, (sum > max), sum, max);
2958  }
2959 
2960  return max;
2961 }
2962 
2963 
2964 template <int dim>
2965 inline adouble
2967 {
2969  for (unsigned int i = 0; i < dim; ++i)
2970  {
2972  for (unsigned int j = 0; j < dim; ++j)
2973  sum += std::fabs(t[i][j]);
2974 
2975  condassign(max, (sum > max), sum, max);
2976  }
2977 
2978  return max;
2979 }
2980 
2981 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
2982 
2983 
2984 #endif // DOXYGEN
2985 
2987 
2988 #endif
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
constexpr Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2672
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2423
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2811
static constexpr const T & value(const T &t)
Definition: numbers.h:693
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:915
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
constexpr Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2759
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
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:577
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1853
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2520
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2894
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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:2545
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2920
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:803
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
STL namespace.
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:503
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
static real_type abs(const number &x)
Definition: numbers.h:599
void unroll_recursion(Vector< OtherNumber > &result, unsigned int &start_index) const
Definition: point.h:110
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
std::size_t size() const
Definition: array_view.h:574
Number * end_raw()
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr value_type & operator[](const unsigned int i)
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
static const char T
#define Assert(cond, exc)
Definition: exceptions.h:1465
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2653
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
void serialize(Archive &ar, const unsigned int version)
constexpr bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:95
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:496
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:2344
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
Expression fabs(const Expression &x)
static constexpr unsigned int rank
Definition: tensor.h:482
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2790
constexpr Tensor & operator*=(const OtherNumber &factor)
constexpr Tensor & operator/=(const OtherNumber &factor)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static constexpr std::size_t memory_consumption()
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:2460
static const char A
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:126
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:1914
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:2021
Definition: tensor.h:462
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:1930
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:931
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
Expression operator-(Expression lhs, const Expression &rhs)
static ::ExceptionBase & ExcNotImplemented()
void unroll(Vector< OtherNumber > &result) const
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
numbers::NumberTraits< Number >::real_type norm() const
static constexpr unsigned int n_independent_components
Definition: tensor.h:488
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:2133
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:2489
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:2269
Expression operator!=(const Expression &lhs, const Expression &rhs)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
static ::ExceptionBase & ExcInternalError()
constexpr Tensor()
constexpr Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2579