Reference documentation for deal.II version Git 497f915867 2021-09-17 22:46:48 +0200
\(\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 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
170 
173  constexpr DEAL_II_CUDA_HOST_DEV
174  Tensor(const Tensor<0, dim, Number> &other);
175 
179  constexpr DEAL_II_CUDA_HOST_DEV
180  Tensor(Tensor<0, dim, Number> &&other) noexcept;
181 #endif
182 
191  DEAL_II_DEPRECATED_EARLY
192  Number *
193  begin_raw();
194 
203  DEAL_II_DEPRECATED_EARLY
204  const Number *
205  begin_raw() const;
206 
215  DEAL_II_DEPRECATED_EARLY
216  Number *
217  end_raw();
218 
228  DEAL_II_DEPRECATED_EARLY
229  const Number *
230  end_raw() const;
231 
241  constexpr DEAL_II_CUDA_HOST_DEV
242  operator Number &();
243 
252  constexpr DEAL_II_CUDA_HOST_DEV operator const Number &() const;
253 
261  template <typename OtherNumber>
262  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
263  operator=(const Tensor<0, dim, OtherNumber> &rhs);
264 
265 #if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
266 
274  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
275  operator=(const Tensor<0, dim, Number> &rhs);
276 #endif
277 
278 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
279 
283  operator=(Tensor<0, dim, Number> &&other) noexcept;
284 #endif
285 
292  template <typename OtherNumber>
293  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
294  operator=(const OtherNumber &d);
295 
299  template <typename OtherNumber>
300  constexpr bool
301  operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
302 
306  template <typename OtherNumber>
307  constexpr bool
308  operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
309 
315  template <typename OtherNumber>
316  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
317  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
318 
324  template <typename OtherNumber>
325  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
326  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
327 
333  template <typename OtherNumber>
334  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
335  operator*=(const OtherNumber &factor);
336 
342  template <typename OtherNumber>
343  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
344  operator/=(const OtherNumber &factor);
345 
351  constexpr DEAL_II_CUDA_HOST_DEV Tensor
352  operator-() const;
353 
366  constexpr void
367  clear();
368 
374  real_type
375  norm() const;
376 
384  norm_square() const;
385 
393  template <class Iterator>
394  void
395  unroll(const Iterator begin, const Iterator end) const;
396 
402  template <class Archive>
403  void
404  serialize(Archive &ar, const unsigned int version);
405 
410  using tensor_type = Number;
411 
412 private:
416  Number value;
417 
421  template <typename Iterator>
422  Iterator
423  unroll_recursion(const Iterator current, const Iterator end) const;
424 
425  // Allow an arbitrary Tensor to access the underlying values.
426  template <int, int, typename>
427  friend class Tensor;
428 };
429 
430 
431 
505 template <int rank_, int dim, typename Number>
506 class Tensor
507 {
508 public:
509  static_assert(rank_ >= 1,
510  "Tensors must have a rank greater than or equal to one.");
511  static_assert(dim >= 0,
512  "Tensors must have a dimension greater than or equal to one.");
521  static constexpr unsigned int dimension = dim;
522 
526  static constexpr unsigned int rank = rank_;
527 
532  static constexpr unsigned int n_independent_components =
533  Tensor<rank_ - 1, dim>::n_independent_components * dim;
534 
540  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
541 
546  using array_type =
547  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
548 
555  Tensor();
556 
562  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
563  const array_type &initializer);
564 
578  template <typename ElementType, typename MemorySpace>
579  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
580  const ArrayView<ElementType, MemorySpace> &initializer);
581 
589  template <typename OtherNumber>
590  constexpr DEAL_II_CUDA_HOST_DEV
591  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
592 
596  template <typename OtherNumber>
597  constexpr Tensor(
598  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
599 
603  template <typename OtherNumber>
604  constexpr
605  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
606 
607 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
608 
611  constexpr Tensor(const Tensor<rank_, dim, Number> &);
612 
616  constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
617 #endif
618 
625  operator[](const unsigned int i);
626 
632  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
633  operator[](const unsigned int i) const;
634 
638  constexpr const Number &
639  operator[](const TableIndices<rank_> &indices) const;
640 
644  constexpr Number &
645  operator[](const TableIndices<rank_> &indices);
646 
650  Number *
651  begin_raw();
652 
656  const Number *
657  begin_raw() const;
658 
662  Number *
663  end_raw();
664 
668  const Number *
669  end_raw() const;
670 
678  template <typename OtherNumber>
679  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
680  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
681 
688  constexpr Tensor &
689  operator=(const Number &d);
690 
691 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
692 
695  constexpr Tensor<rank_, dim, Number> &
696  operator=(const Tensor<rank_, dim, Number> &);
697 
701  constexpr Tensor<rank_, dim, Number> &
702  operator=(Tensor<rank_, dim, Number> &&) noexcept;
703 #endif
704 
708  template <typename OtherNumber>
709  constexpr bool
711 
715  template <typename OtherNumber>
716  constexpr bool
718 
724  template <typename OtherNumber>
725  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
726  operator+=(const Tensor<rank_, dim, OtherNumber> &);
727 
733  template <typename OtherNumber>
734  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
735  operator-=(const Tensor<rank_, dim, OtherNumber> &);
736 
743  template <typename OtherNumber>
744  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
745  operator*=(const OtherNumber &factor);
746 
752  template <typename OtherNumber>
753  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
754  operator/=(const OtherNumber &factor);
755 
761  constexpr DEAL_II_CUDA_HOST_DEV Tensor
762  operator-() const;
763 
776  constexpr void
777  clear();
778 
788  norm() const;
789 
796  constexpr DEAL_II_CUDA_HOST_DEV
798  norm_square() const;
799 
810  template <typename OtherNumber>
811  DEAL_II_DEPRECATED_EARLY void
812  unroll(Vector<OtherNumber> &result) const;
813 
824  template <class Iterator>
825  void
826  unroll(const Iterator begin, const Iterator end) const;
827 
832  static constexpr unsigned int
833  component_to_unrolled_index(const TableIndices<rank_> &indices);
834 
840  static constexpr TableIndices<rank_>
841  unrolled_to_component_indices(const unsigned int i);
842 
847  static constexpr std::size_t
849 
855  template <class Archive>
856  void
857  serialize(Archive &ar, const unsigned int version);
858 
864 
865 private:
869  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
870  // ... avoid a compiler warning in case of dim == 0 and ensure that the
871  // array always has positive size.
872 
876  template <typename Iterator>
877  Iterator
878  unroll_recursion(const Iterator current, const Iterator end) const;
879 
886  template <typename ArrayLike, std::size_t... Indices>
887  constexpr DEAL_II_CUDA_HOST_DEV
888  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
889 
890  // Allow an arbitrary Tensor to access the underlying values.
891  template <int, int, typename>
892  friend class Tensor;
893 
894  // Point is allowed access to the coordinates. This is supposed to improve
895  // speed.
896  friend class Point<dim, Number>;
897 };
898 
899 
900 #ifndef DOXYGEN
901 namespace internal
902 {
903  // Workaround: The following 4 overloads are necessary to be able to
904  // compile the library with Apple Clang 8 and older. We should remove
905  // these overloads again when we bump the minimal required version to
906  // something later than clang-3.6 / Apple Clang 6.3.
907  template <int rank, int dim, typename T, typename U>
908  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
909  {
910  using type =
912  };
913 
914  template <int rank, int dim, typename T, typename U>
915  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
916  {
917  using type =
919  };
920 
921  template <typename T, int rank, int dim, typename U>
922  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
923  {
924  using type =
926  };
927 
928  template <int rank, int dim, typename T, typename U>
929  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
930  {
931  using type =
933  };
934  // end workaround
935 
940  template <int rank, int dim, typename T>
941  struct NumberType<Tensor<rank, dim, T>>
942  {
943  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
944  value(const Tensor<rank, dim, T> &t)
945  {
946  return t;
947  }
948 
950  value(const T &t)
951  {
953  tmp = t;
954  return tmp;
955  }
956  };
957 } // namespace internal
958 
959 
960 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
961 
962 
963 template <int dim, typename Number>
966  // Some auto-differentiable numbers need explicit
967  // zero initialization such as adtl::adouble.
968  : Tensor{0.0}
969 {}
970 
971 
972 
973 template <int dim, typename Number>
974 template <typename OtherNumber>
976 Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
977  : value(internal::NumberType<Number>::value(initializer))
978 {}
979 
980 
981 
982 template <int dim, typename Number>
983 template <typename OtherNumber>
986  : Tensor{p.value}
987 {}
988 
989 
990 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
991 template <int dim, typename Number>
994  : value{other.value}
995 {}
996 
997 
998 
999 template <int dim, typename Number>
1002  : value{std::move(other.value)}
1003 {}
1004 # endif
1005 
1006 
1007 template <int dim, typename Number>
1008 inline Number *
1010 {
1011  return std::addressof(value);
1012 }
1013 
1014 
1015 
1016 template <int dim, typename Number>
1017 inline const Number *
1019 {
1020  return std::addressof(value);
1021 }
1022 
1023 
1024 
1025 template <int dim, typename Number>
1026 inline Number *
1028 {
1030 }
1031 
1032 
1033 
1034 template <int dim, typename Number>
1035 const Number *
1037 {
1039 }
1040 
1041 
1042 
1043 template <int dim, typename Number>
1046 {
1047  // We cannot use Assert inside a CUDA kernel
1048 # ifndef __CUDA_ARCH__
1049  Assert(dim != 0,
1050  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1051 # endif
1052  return value;
1053 }
1054 
1055 
1056 template <int dim, typename Number>
1057 constexpr inline DEAL_II_ALWAYS_INLINE
1059 {
1060  // We cannot use Assert inside a CUDA kernel
1061 # ifndef __CUDA_ARCH__
1062  Assert(dim != 0,
1063  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1064 # endif
1065  return value;
1066 }
1067 
1068 
1069 
1070 template <int dim, typename Number>
1071 template <typename OtherNumber>
1072 constexpr inline DEAL_II_ALWAYS_INLINE
1075 {
1077  return *this;
1078 }
1079 
1080 
1081 # if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1082 template <int dim, typename Number>
1083 constexpr inline DEAL_II_ALWAYS_INLINE
1086 {
1087  value = p.value;
1088  return *this;
1089 }
1090 # endif
1091 
1092 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1093 template <int dim, typename Number>
1096 {
1097  value = std::move(other.value);
1098  return *this;
1099 }
1100 # endif
1101 
1102 
1103 
1104 template <int dim, typename Number>
1105 template <typename OtherNumber>
1106 constexpr inline DEAL_II_ALWAYS_INLINE
1108  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
1109 {
1111  return *this;
1112 }
1113 
1114 
1115 template <int dim, typename Number>
1116 template <typename OtherNumber>
1117 constexpr inline bool
1119 {
1120 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1121  Assert(!(std::is_same<Number, adouble>::value ||
1122  std::is_same<OtherNumber, adouble>::value),
1123  ExcMessage(
1124  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1125  "been extended to support advanced branching."));
1126 # endif
1127 
1128  return numbers::values_are_equal(value, p.value);
1129 }
1130 
1131 
1132 template <int dim, typename Number>
1133 template <typename OtherNumber>
1134 constexpr bool
1136 {
1137  return !((*this) == p);
1138 }
1139 
1140 
1141 template <int dim, typename Number>
1142 template <typename OtherNumber>
1143 constexpr inline DEAL_II_ALWAYS_INLINE
1146 {
1147  value += p.value;
1148  return *this;
1149 }
1150 
1151 
1152 template <int dim, typename Number>
1153 template <typename OtherNumber>
1154 constexpr inline DEAL_II_ALWAYS_INLINE
1157 {
1158  value -= p.value;
1159  return *this;
1160 }
1161 
1162 
1163 
1164 namespace internal
1165 {
1166  namespace ComplexWorkaround
1167  {
1168  template <typename Number, typename OtherNumber>
1169  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1170  multiply_assign_scalar(Number &val, const OtherNumber &s)
1171  {
1172  val *= s;
1173  }
1174 
1175 # ifdef __CUDA_ARCH__
1176  template <typename Number, typename OtherNumber>
1177  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1178  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1179  {
1180  printf("This function is not implemented for std::complex<Number>!\n");
1181  assert(false);
1182  }
1183 # endif
1184  } // namespace ComplexWorkaround
1185 } // namespace internal
1186 
1187 
1188 template <int dim, typename Number>
1189 template <typename OtherNumber>
1190 constexpr inline DEAL_II_ALWAYS_INLINE
1192  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1193 {
1194  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1195  return *this;
1196 }
1197 
1198 
1199 
1200 template <int dim, typename Number>
1201 template <typename OtherNumber>
1203 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1204 {
1205  value /= s;
1206  return *this;
1207 }
1208 
1209 
1210 template <int dim, typename Number>
1213 {
1214  return -value;
1215 }
1216 
1217 
1218 template <int dim, typename Number>
1221 {
1222  Assert(dim != 0,
1223  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1224  return numbers::NumberTraits<Number>::abs(value);
1225 }
1226 
1227 
1228 template <int dim, typename Number>
1232 {
1233  // We cannot use Assert inside a CUDA kernel
1234 # ifndef __CUDA_ARCH__
1235  Assert(dim != 0,
1236  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1237 # endif
1239 }
1240 
1241 
1242 
1243 template <int dim, typename Number>
1244 template <typename Iterator>
1245 Iterator
1246 Tensor<0, dim, Number>::unroll_recursion(const Iterator current,
1247  const Iterator end) const
1248 {
1249  (void)end;
1250  Assert(dim != 0,
1251  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1252  Assert(std::distance(current, end) >= 1,
1253  ExcMessage("The provided iterator range must contain at least one "
1254  "element."));
1255  *current = value;
1256  return std::next(current);
1257 }
1258 
1259 
1260 
1261 template <int dim, typename Number>
1262 constexpr inline void
1264 {
1265  // Some auto-differentiable numbers need explicit
1266  // zero initialization.
1268 }
1269 
1270 
1271 
1272 template <int dim, typename Number>
1273 template <class Iterator>
1274 inline void
1275 Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1276 {
1277  AssertDimension(std::distance(begin, end), n_independent_components);
1278  unroll_recursion(begin, end);
1279 }
1280 
1281 
1282 
1283 template <int dim, typename Number>
1284 template <class Archive>
1285 inline void
1286 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1287 {
1288  ar &value;
1289 }
1290 
1291 
1292 template <int dim, typename Number>
1294 
1295 
1296 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1297 
1298 template <int rank_, int dim, typename Number>
1299 template <typename ArrayLike, std::size_t... indices>
1301 Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1302  std::index_sequence<indices...>)
1303  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1304 {
1305  static_assert(sizeof...(indices) == dim,
1306  "dim should match the number of indices");
1307 }
1308 
1309 
1310 
1311 template <int rank_, int dim, typename Number>
1314  // We would like to use =default, but this causes compile errors with some
1315  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1316  : values{}
1317 {}
1318 
1319 
1320 
1321 template <int rank_, int dim, typename Number>
1324  : Tensor(initializer, std::make_index_sequence<dim>{})
1325 {}
1326 
1327 
1328 
1329 template <int rank_, int dim, typename Number>
1330 template <typename ElementType, typename MemorySpace>
1333  const ArrayView<ElementType, MemorySpace> &initializer)
1334 {
1336 
1337  for (unsigned int i = 0; i < n_independent_components; ++i)
1338  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1339 }
1340 
1341 
1342 
1343 template <int rank_, int dim, typename Number>
1344 template <typename OtherNumber>
1347  const Tensor<rank_, dim, OtherNumber> &initializer)
1348  : Tensor(initializer, std::make_index_sequence<dim>{})
1349 {}
1350 
1351 
1352 
1353 template <int rank_, int dim, typename Number>
1354 template <typename OtherNumber>
1355 constexpr DEAL_II_ALWAYS_INLINE
1357  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1358  : Tensor(initializer, std::make_index_sequence<dim>{})
1359 {}
1360 
1361 
1362 
1363 template <int rank_, int dim, typename Number>
1364 template <typename OtherNumber>
1366 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1367 {
1368  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
1369 }
1370 
1371 
1372 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1373 template <int rank_, int dim, typename Number>
1374 constexpr DEAL_II_ALWAYS_INLINE
1376 {
1377  for (unsigned int i = 0; i < dim; ++i)
1378  values[i] = other.values[i];
1379 }
1380 
1381 
1382 
1383 template <int rank_, int dim, typename Number>
1384 constexpr DEAL_II_ALWAYS_INLINE
1386 {
1387  for (unsigned int i = 0; i < dim; ++i)
1388  values[i] = other.values[i];
1389 }
1390 # endif
1391 
1392 
1393 namespace internal
1394 {
1395  namespace TensorSubscriptor
1396  {
1397  template <typename ArrayElementType, int dim>
1398  constexpr inline DEAL_II_ALWAYS_INLINE
1399  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1400  subscript(ArrayElementType * values,
1401  const unsigned int i,
1402  std::integral_constant<int, dim>)
1403  {
1404  // We cannot use Assert in a CUDA kernel
1405 # ifndef __CUDA_ARCH__
1406  AssertIndexRange(i, dim);
1407 # endif
1408  return values[i];
1409  }
1410 
1411  // The variables within this struct will be referenced in the next function.
1412  // It is a workaround that allows returning a reference to a static variable
1413  // while allowing constexpr evaluation of the function.
1414  // It has to be defined outside the function because constexpr functions
1415  // cannot define static variables
1416  template <typename ArrayElementType>
1417  struct Uninitialized
1418  {
1419  static ArrayElementType value;
1420  };
1421 
1422  template <typename Type>
1423  Type Uninitialized<Type>::value;
1424 
1425  template <typename ArrayElementType>
1426  constexpr inline DEAL_II_ALWAYS_INLINE
1427  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1428  subscript(ArrayElementType *,
1429  const unsigned int,
1430  std::integral_constant<int, 0>)
1431  {
1432  // We cannot use Assert in a CUDA kernel
1433 # ifndef __CUDA_ARCH__
1434  Assert(
1435  false,
1436  ExcMessage(
1437  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1438 # endif
1439  return Uninitialized<ArrayElementType>::value;
1440  }
1441  } // namespace TensorSubscriptor
1442 } // namespace internal
1443 
1444 
1445 template <int rank_, int dim, typename Number>
1448  Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1449 {
1450  return ::internal::TensorSubscriptor::subscript(
1451  values, i, std::integral_constant<int, dim>());
1452 }
1453 
1454 
1455 template <int rank_, int dim, typename Number>
1456 constexpr DEAL_II_ALWAYS_INLINE
1458  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1459 {
1460 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1461  AssertIndexRange(i, dim);
1462 # endif
1463 
1464  return values[i];
1465 }
1466 
1467 
1468 template <int rank_, int dim, typename Number>
1469 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1471 {
1472 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1473  Assert(dim != 0,
1474  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1475 # endif
1476 
1477  return TensorAccessors::extract<rank_>(*this, indices);
1478 }
1479 
1480 
1481 
1482 template <int rank_, int dim, typename Number>
1483 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1485 {
1486 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1487  Assert(dim != 0,
1488  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1489 # endif
1490 
1491  return TensorAccessors::extract<rank_>(*this, indices);
1492 }
1493 
1494 
1495 
1496 template <int rank_, int dim, typename Number>
1497 inline Number *
1499 {
1500  return std::addressof(
1501  this->operator[](this->unrolled_to_component_indices(0)));
1502 }
1503 
1504 
1505 
1506 template <int rank_, int dim, typename Number>
1507 inline const Number *
1509 {
1510  return std::addressof(
1511  this->operator[](this->unrolled_to_component_indices(0)));
1512 }
1513 
1514 
1515 
1516 template <int rank_, int dim, typename Number>
1517 inline Number *
1519 {
1521 }
1522 
1523 
1524 
1525 template <int rank_, int dim, typename Number>
1526 inline const Number *
1528 {
1530 }
1531 
1532 
1533 
1534 template <int rank_, int dim, typename Number>
1535 template <typename OtherNumber>
1538 {
1539  // The following loop could be written more concisely using std::copy, but
1540  // that function is only constexpr from C++20 on.
1541  for (unsigned int i = 0; i < dim; ++i)
1542  values[i] = t.values[i];
1543  return *this;
1544 }
1545 
1546 
1547 
1548 template <int rank_, int dim, typename Number>
1551 {
1553  (void)d;
1554 
1555  for (unsigned int i = 0; i < dim; ++i)
1557  return *this;
1558 }
1559 
1560 
1561 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1562 template <int rank_, int dim, typename Number>
1565 {
1566  for (unsigned int i = 0; i < dim; ++i)
1567  values[i] = other.values[i];
1568  return *this;
1569 }
1570 
1571 
1572 
1573 template <int rank_, int dim, typename Number>
1576  Tensor<rank_, dim, Number> &&other) noexcept
1577 {
1578  for (unsigned int i = 0; i < dim; ++i)
1579  values[i] = other.values[i];
1580  return *this;
1581 }
1582 # endif
1583 
1584 
1585 template <int rank_, int dim, typename Number>
1586 template <typename OtherNumber>
1587 constexpr inline bool
1589  const Tensor<rank_, dim, OtherNumber> &p) const
1590 {
1591  for (unsigned int i = 0; i < dim; ++i)
1592  if (values[i] != p.values[i])
1593  return false;
1594  return true;
1595 }
1596 
1597 
1598 // At some places in the library, we have Point<0> for formal reasons
1599 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1600 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1601 // in the above function that the loop end check always fails, we
1602 // implement this function here
1603 template <>
1604 template <>
1605 constexpr inline bool
1607 {
1608  return true;
1609 }
1610 
1611 
1612 template <int rank_, int dim, typename Number>
1613 template <typename OtherNumber>
1614 constexpr bool
1616  const Tensor<rank_, dim, OtherNumber> &p) const
1617 {
1618  return !((*this) == p);
1619 }
1620 
1621 
1622 template <int rank_, int dim, typename Number>
1623 template <typename OtherNumber>
1624 constexpr inline DEAL_II_ALWAYS_INLINE
1628 {
1629  for (unsigned int i = 0; i < dim; ++i)
1630  values[i] += p.values[i];
1631  return *this;
1632 }
1633 
1634 
1635 template <int rank_, int dim, typename Number>
1636 template <typename OtherNumber>
1637 constexpr inline DEAL_II_ALWAYS_INLINE
1641 {
1642  for (unsigned int i = 0; i < dim; ++i)
1643  values[i] -= p.values[i];
1644  return *this;
1645 }
1646 
1647 
1648 template <int rank_, int dim, typename Number>
1649 template <typename OtherNumber>
1650 constexpr inline DEAL_II_ALWAYS_INLINE
1652  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1653 {
1654  for (unsigned int i = 0; i < dim; ++i)
1655  values[i] *= s;
1656  return *this;
1657 }
1658 
1659 
1660 namespace internal
1661 {
1662  namespace TensorImplementation
1663  {
1664  template <int rank,
1665  int dim,
1666  typename Number,
1667  typename OtherNumber,
1668  typename std::enable_if<
1669  !std::is_integral<
1670  typename ProductType<Number, OtherNumber>::type>::value &&
1671  !std::is_same<Number, Differentiation::SD::Expression>::value,
1672  int>::type = 0>
1673  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1675  const OtherNumber &factor)
1676  {
1677  const Number inverse_factor = Number(1.) / factor;
1678  // recurse over the base objects
1679  for (unsigned int d = 0; d < dim; ++d)
1680  t[d] *= inverse_factor;
1681  }
1682 
1683 
1684  template <int rank,
1685  int dim,
1686  typename Number,
1687  typename OtherNumber,
1688  typename std::enable_if<
1689  std::is_integral<
1690  typename ProductType<Number, OtherNumber>::type>::value ||
1691  std::is_same<Number, Differentiation::SD::Expression>::value,
1692  int>::type = 0>
1693  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1695  const OtherNumber &factor)
1696  {
1697  // recurse over the base objects
1698  for (unsigned int d = 0; d < dim; ++d)
1699  t[d] /= factor;
1700  }
1701  } // namespace TensorImplementation
1702 } // namespace internal
1703 
1704 
1705 template <int rank_, int dim, typename Number>
1706 template <typename OtherNumber>
1707 constexpr inline DEAL_II_ALWAYS_INLINE
1709  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1710 {
1712  return *this;
1713 }
1714 
1715 
1716 template <int rank_, int dim, typename Number>
1717 constexpr inline DEAL_II_ALWAYS_INLINE
1720 {
1722 
1723  for (unsigned int i = 0; i < dim; ++i)
1724  tmp.values[i] = -values[i];
1725 
1726  return tmp;
1727 }
1728 
1729 
1730 template <int rank_, int dim, typename Number>
1733 {
1734  return std::sqrt(norm_square());
1735 }
1736 
1737 
1738 template <int rank_, int dim, typename Number>
1742 {
1744  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1745  for (unsigned int i = 0; i < dim; ++i)
1746  s += values[i].norm_square();
1747 
1748  return s;
1749 }
1750 
1751 
1752 
1753 template <int rank_, int dim, typename Number>
1754 template <typename OtherNumber>
1755 inline void
1756 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
1757 {
1758  unroll(result.begin(), result.end());
1759 }
1760 
1761 
1762 
1763 template <int rank_, int dim, typename Number>
1764 template <class Iterator>
1765 inline void
1766 Tensor<rank_, dim, Number>::unroll(const Iterator begin,
1767  const Iterator end) const
1768 {
1769  AssertDimension(std::distance(begin, end), n_independent_components);
1770  unroll_recursion(begin, end);
1771 }
1772 
1773 
1774 
1775 template <int rank_, int dim, typename Number>
1776 template <typename Iterator>
1777 Iterator
1778 Tensor<rank_, dim, Number>::unroll_recursion(const Iterator current,
1779  const Iterator end) const
1780 {
1781  auto next = current;
1782  for (unsigned int i = 0; i < dim; ++i)
1783  next = values[i].unroll_recursion(next, end);
1784  return next;
1785 }
1786 
1787 
1788 template <int rank_, int dim, typename Number>
1789 constexpr inline unsigned int
1791  const TableIndices<rank_> &indices)
1792 {
1793  unsigned int index = 0;
1794  for (int r = 0; r < rank_; ++r)
1795  index = index * dim + indices[r];
1796 
1797  return index;
1798 }
1799 
1800 
1801 
1802 namespace internal
1803 {
1804  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1805  // and rank=2. Make sure we don't have compiler warnings.
1806 
1807  template <int dim>
1808  inline constexpr unsigned int
1809  mod(const unsigned int x)
1810  {
1811  return x % dim;
1812  }
1813 
1814  template <>
1815  inline unsigned int
1816  mod<0>(const unsigned int x)
1817  {
1818  Assert(false, ExcInternalError());
1819  return x;
1820  }
1821 
1822  template <int dim>
1823  inline constexpr unsigned int
1824  div(const unsigned int x)
1825  {
1826  return x / dim;
1827  }
1828 
1829  template <>
1830  inline unsigned int
1831  div<0>(const unsigned int x)
1832  {
1833  Assert(false, ExcInternalError());
1834  return x;
1835  }
1836 
1837 } // namespace internal
1838 
1839 
1840 
1841 template <int rank_, int dim, typename Number>
1842 constexpr inline TableIndices<rank_>
1844 {
1845  AssertIndexRange(i, n_independent_components);
1846 
1847  TableIndices<rank_> indices;
1848 
1849  unsigned int remainder = i;
1850  for (int r = rank_ - 1; r >= 0; --r)
1851  {
1852  indices[r] = internal::mod<dim>(remainder);
1853  remainder = internal::div<dim>(remainder);
1854  }
1855  Assert(remainder == 0, ExcInternalError());
1856 
1857  return indices;
1858 }
1859 
1860 
1861 template <int rank_, int dim, typename Number>
1862 constexpr inline void
1864 {
1865  for (unsigned int i = 0; i < dim; ++i)
1867 }
1868 
1869 
1870 template <int rank_, int dim, typename Number>
1871 constexpr std::size_t
1873 {
1874  return sizeof(Tensor<rank_, dim, Number>);
1875 }
1876 
1877 
1878 template <int rank_, int dim, typename Number>
1879 template <class Archive>
1880 inline void
1881 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1882 {
1883  ar &values;
1884 }
1885 
1886 
1887 template <int rank_, int dim, typename Number>
1889 
1890 #endif // DOXYGEN
1891 
1892 /* ----------------- Non-member functions operating on tensors. ------------ */
1893 
1898 
1906 template <int rank_, int dim, typename Number>
1907 inline std::ostream &
1908 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1909 {
1910  for (unsigned int i = 0; i < dim; ++i)
1911  {
1912  out << p[i];
1913  if (i != dim - 1)
1914  out << ' ';
1915  }
1916 
1917  return out;
1918 }
1919 
1920 
1927 template <int dim, typename Number>
1928 inline std::ostream &
1929 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1930 {
1931  out << static_cast<const Number &>(p);
1932  return out;
1933 }
1934 
1935 
1937 
1941 
1942 
1953 template <int dim, typename Number, typename Other>
1956  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1957 {
1958  return object * static_cast<const Number &>(t);
1959 }
1960 
1961 
1962 
1973 template <int dim, typename Number, typename Other>
1976  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1977 {
1978  return static_cast<const Number &>(t) * object;
1979 }
1980 
1981 
1993 template <int dim, typename Number, typename OtherNumber>
1997  const Tensor<0, dim, OtherNumber> &src2)
1998 {
1999  return static_cast<const Number &>(src1) *
2000  static_cast<const OtherNumber &>(src2);
2001 }
2002 
2003 
2011 template <int dim, typename Number, typename OtherNumber>
2013  Tensor<0,
2014  dim,
2015  typename ProductType<Number,
2016  typename EnableIfScalar<OtherNumber>::type>::type>
2017  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2018 {
2019  return static_cast<const Number &>(t) / factor;
2020 }
2021 
2022 
2030 template <int dim, typename Number, typename OtherNumber>
2034  const Tensor<0, dim, OtherNumber> &q)
2035 {
2036  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2037 }
2038 
2039 
2047 template <int dim, typename Number, typename OtherNumber>
2051  const Tensor<0, dim, OtherNumber> &q)
2052 {
2053  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2054 }
2055 
2056 
2069 template <int rank, int dim, typename Number, typename OtherNumber>
2071  Tensor<rank,
2072  dim,
2073  typename ProductType<Number,
2074  typename EnableIfScalar<OtherNumber>::type>::type>
2075  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2076 {
2077  // recurse over the base objects
2079  for (unsigned int d = 0; d < dim; ++d)
2080  tt[d] = t[d] * factor;
2081  return tt;
2082 }
2083 
2084 
2097 template <int rank, int dim, typename Number, typename OtherNumber>
2099  Tensor<rank,
2100  dim,
2102  OtherNumber>::type>
2103  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2104 {
2105  // simply forward to the operator above
2106  return t * factor;
2107 }
2108 
2109 
2110 namespace internal
2111 {
2112  namespace TensorImplementation
2113  {
2114  template <int rank,
2115  int dim,
2116  typename Number,
2117  typename OtherNumber,
2118  typename std::enable_if<
2119  !std::is_integral<
2120  typename ProductType<Number, OtherNumber>::type>::value,
2121  int>::type = 0>
2125  const OtherNumber & factor)
2126  {
2128  const Number inverse_factor = Number(1.) / factor;
2129  // recurse over the base objects
2130  for (unsigned int d = 0; d < dim; ++d)
2131  tt[d] = t[d] * inverse_factor;
2132  return tt;
2133  }
2134 
2135 
2136  template <int rank,
2137  int dim,
2138  typename Number,
2139  typename OtherNumber,
2140  typename std::enable_if<
2141  std::is_integral<
2142  typename ProductType<Number, OtherNumber>::type>::value,
2143  int>::type = 0>
2147  const OtherNumber & factor)
2148  {
2150  // recurse over the base objects
2151  for (unsigned int d = 0; d < dim; ++d)
2152  tt[d] = t[d] / factor;
2153  return tt;
2154  }
2155  } // namespace TensorImplementation
2156 } // namespace internal
2157 
2158 
2168 template <int rank, int dim, typename Number, typename OtherNumber>
2170  Tensor<rank,
2171  dim,
2172  typename ProductType<Number,
2173  typename EnableIfScalar<OtherNumber>::type>::type>
2174  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2175 {
2177 }
2178 
2179 
2189 template <int rank, int dim, typename Number, typename OtherNumber>
2194 {
2196 
2197  for (unsigned int i = 0; i < dim; ++i)
2198  tmp[i] += q[i];
2199 
2200  return tmp;
2201 }
2202 
2203 
2213 template <int rank, int dim, typename Number, typename OtherNumber>
2218 {
2220 
2221  for (unsigned int i = 0; i < dim; ++i)
2222  tmp[i] -= q[i];
2223 
2224  return tmp;
2225 }
2226 
2233 template <int dim, typename Number, typename OtherNumber>
2234 inline constexpr DEAL_II_ALWAYS_INLINE
2237  const Tensor<0, dim, OtherNumber> &src2)
2238 {
2240 
2241  tmp *= src2;
2242 
2243  return tmp;
2244 }
2245 
2262 template <int rank, int dim, typename Number, typename OtherNumber>
2263 inline constexpr DEAL_II_ALWAYS_INLINE
2266  const Tensor<rank, dim, OtherNumber> &src2)
2267 {
2269 
2270  for (unsigned int i = 0; i < dim; ++i)
2271  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2273 
2274  return tmp;
2275 }
2276 
2278 
2282 
2283 
2306 template <int rank_1,
2307  int rank_2,
2308  int dim,
2309  typename Number,
2310  typename OtherNumber,
2311  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
2312 constexpr inline DEAL_II_ALWAYS_INLINE
2313  typename Tensor<rank_1 + rank_2 - 2,
2314  dim,
2315  typename ProductType<Number, OtherNumber>::type>::tensor_type
2318 {
2319  typename Tensor<rank_1 + rank_2 - 2,
2320  dim,
2321  typename ProductType<Number, OtherNumber>::type>::tensor_type
2322  result{};
2323 
2324  TensorAccessors::internal::
2325  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2326  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2327  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2328 
2329  return result;
2330 }
2331 
2332 
2361 template <int index_1,
2362  int index_2,
2363  int rank_1,
2364  int rank_2,
2365  int dim,
2366  typename Number,
2367  typename OtherNumber>
2368 constexpr inline DEAL_II_ALWAYS_INLINE
2369  typename Tensor<rank_1 + rank_2 - 2,
2370  dim,
2371  typename ProductType<Number, OtherNumber>::type>::tensor_type
2374 {
2375  Assert(0 <= index_1 && index_1 < rank_1,
2376  ExcMessage(
2377  "The specified index_1 must lie within the range [0,rank_1)"));
2378  Assert(0 <= index_2 && index_2 < rank_2,
2379  ExcMessage(
2380  "The specified index_2 must lie within the range [0,rank_2)"));
2381 
2382  using namespace TensorAccessors;
2383  using namespace TensorAccessors::internal;
2384 
2385  // Reorder index_1 to the end of src1:
2387  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2388 
2389  // Reorder index_2 to the end of src2:
2390  const ReorderedIndexView<index_2,
2391  rank_2,
2393  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2394 
2395  typename Tensor<rank_1 + rank_2 - 2,
2396  dim,
2397  typename ProductType<Number, OtherNumber>::type>::tensor_type
2398  result{};
2399  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2400  return result;
2401 }
2402 
2403 
2434 template <int index_1,
2435  int index_2,
2436  int index_3,
2437  int index_4,
2438  int rank_1,
2439  int rank_2,
2440  int dim,
2441  typename Number,
2442  typename OtherNumber>
2443 constexpr inline
2444  typename Tensor<rank_1 + rank_2 - 4,
2445  dim,
2446  typename ProductType<Number, OtherNumber>::type>::tensor_type
2449 {
2450  Assert(0 <= index_1 && index_1 < rank_1,
2451  ExcMessage(
2452  "The specified index_1 must lie within the range [0,rank_1)"));
2453  Assert(0 <= index_3 && index_3 < rank_1,
2454  ExcMessage(
2455  "The specified index_3 must lie within the range [0,rank_1)"));
2456  Assert(index_1 != index_3,
2457  ExcMessage("index_1 and index_3 must not be the same"));
2458  Assert(0 <= index_2 && index_2 < rank_2,
2459  ExcMessage(
2460  "The specified index_2 must lie within the range [0,rank_2)"));
2461  Assert(0 <= index_4 && index_4 < rank_2,
2462  ExcMessage(
2463  "The specified index_4 must lie within the range [0,rank_2)"));
2464  Assert(index_2 != index_4,
2465  ExcMessage("index_2 and index_4 must not be the same"));
2466 
2467  using namespace TensorAccessors;
2468  using namespace TensorAccessors::internal;
2469 
2470  // Reorder index_1 to the end of src1:
2472  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2473 
2474  // Reorder index_2 to the end of src2:
2476  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2477 
2478  // Now, reorder index_3 to the end of src1. We have to make sure to
2479  // preserve the original ordering: index_1 has been removed. If
2480  // index_3 > index_1, we have to use (index_3 - 1) instead:
2482  (index_3 < index_1 ? index_3 : index_3 - 1),
2483  rank_1,
2484  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2485  reord_3 =
2486  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2487  index_3 - 1,
2488  rank_1 > (reord_1);
2489 
2490  // Now, reorder index_4 to the end of src2. We have to make sure to
2491  // preserve the original ordering: index_2 has been removed. If
2492  // index_4 > index_2, we have to use (index_4 - 1) instead:
2493  ReorderedIndexView<
2494  (index_4 < index_2 ? index_4 : index_4 - 1),
2495  rank_2,
2496  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2497  reord_4 =
2498  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2499  index_4 - 1,
2500  rank_2 > (reord_2);
2501 
2502  typename Tensor<rank_1 + rank_2 - 4,
2503  dim,
2504  typename ProductType<Number, OtherNumber>::type>::tensor_type
2505  result{};
2506  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2507  return result;
2508 }
2509 
2510 
2523 template <int rank, int dim, typename Number, typename OtherNumber>
2524 constexpr inline DEAL_II_ALWAYS_INLINE
2527  const Tensor<rank, dim, OtherNumber> &right)
2528 {
2529  typename ProductType<Number, OtherNumber>::type result{};
2530  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2531  return result;
2532 }
2533 
2534 
2552 template <template <int, int, typename> class TensorT1,
2553  template <int, int, typename>
2554  class TensorT2,
2555  template <int, int, typename>
2556  class TensorT3,
2557  int rank_1,
2558  int rank_2,
2559  int dim,
2560  typename T1,
2561  typename T2,
2562  typename T3>
2563 constexpr inline DEAL_II_ALWAYS_INLINE
2565  contract3(const TensorT1<rank_1, dim, T1> & left,
2566  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2567  const TensorT3<rank_2, dim, T3> & right)
2568 {
2569  using return_type =
2571  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2572  middle,
2573  right);
2574 }
2575 
2576 
2587 template <int rank_1,
2588  int rank_2,
2589  int dim,
2590  typename Number,
2591  typename OtherNumber>
2592 constexpr inline DEAL_II_ALWAYS_INLINE
2596 {
2597  typename Tensor<rank_1 + rank_2,
2598  dim,
2599  typename ProductType<Number, OtherNumber>::type>::tensor_type
2600  result{};
2601  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2602  return result;
2603 }
2604 
2605 
2607 
2611 
2612 
2623 template <int dim, typename Number>
2626 {
2627  Assert(dim == 2, ExcInternalError());
2628 
2629  Tensor<1, dim, Number> result;
2630 
2631  result[0] = src[1];
2632  result[1] = -src[0];
2633 
2634  return result;
2635 }
2636 
2637 
2647 template <int dim, typename Number1, typename Number2>
2648 constexpr inline DEAL_II_ALWAYS_INLINE
2651  const Tensor<1, dim, Number2> &src2)
2652 {
2653  Assert(dim == 3, ExcInternalError());
2654 
2656 
2657  // avoid compiler warnings
2658  constexpr int s0 = 0 % dim;
2659  constexpr int s1 = 1 % dim;
2660  constexpr int s2 = 2 % dim;
2661 
2662  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2663  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2664  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2665 
2666  return result;
2667 }
2668 
2669 
2671 
2675 
2676 
2682 template <int dim, typename Number>
2683 constexpr inline DEAL_II_ALWAYS_INLINE Number
2685 {
2686  // Compute the determinant using the Laplace expansion of the
2687  // determinant. We expand along the last row.
2688  Number det = internal::NumberType<Number>::value(0.0);
2689 
2690  for (unsigned int k = 0; k < dim; ++k)
2691  {
2692  Tensor<2, dim - 1, Number> minor;
2693  for (unsigned int i = 0; i < dim - 1; ++i)
2694  for (unsigned int j = 0; j < dim - 1; ++j)
2695  minor[i][j] = t[i][j < k ? j : j + 1];
2696 
2697  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2698 
2699  det += t[dim - 1][k] * cofactor;
2700  }
2701 
2702  return ((dim % 2 == 0) ? 1. : -1.) * det;
2703 }
2704 
2710 template <typename Number>
2711 constexpr DEAL_II_ALWAYS_INLINE Number
2713 {
2714  return t[0][0];
2715 }
2716 
2722 template <typename Number>
2723 constexpr DEAL_II_ALWAYS_INLINE Number
2725 {
2726  // hard-coded for efficiency reasons
2727  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2728 }
2729 
2735 template <typename Number>
2736 constexpr DEAL_II_ALWAYS_INLINE Number
2738 {
2739  // hard-coded for efficiency reasons
2740  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2741  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2742  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2743  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2744  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2745  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2746  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2747 }
2748 
2749 
2756 template <int dim, typename Number>
2757 constexpr inline DEAL_II_ALWAYS_INLINE Number
2759 {
2760  Number t = d[0][0];
2761  for (unsigned int i = 1; i < dim; ++i)
2762  t += d[i][i];
2763  return t;
2764 }
2765 
2766 
2775 template <int dim, typename Number>
2776 constexpr inline Tensor<2, dim, Number>
2778 {
2779  Number return_tensor[dim][dim];
2780 
2781  // if desired, take over the
2782  // inversion of a 4x4 tensor
2783  // from the FullMatrix
2784  AssertThrow(false, ExcNotImplemented());
2785 
2786  return Tensor<2, dim, Number>(return_tensor);
2787 }
2788 
2789 
2790 #ifndef DOXYGEN
2791 
2792 template <typename Number>
2794  invert(const Tensor<2, 1, Number> &t)
2795 {
2796  Tensor<2, 1, Number> return_tensor;
2797 
2798  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2799 
2800  return return_tensor;
2801 }
2802 
2803 
2804 template <typename Number>
2806  invert(const Tensor<2, 2, Number> &t)
2807 {
2808  Tensor<2, 2, Number> return_tensor;
2809 
2810  const Number inv_det_t = internal::NumberType<Number>::value(
2811  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2812  return_tensor[0][0] = t[1][1];
2813  return_tensor[0][1] = -t[0][1];
2814  return_tensor[1][0] = -t[1][0];
2815  return_tensor[1][1] = t[0][0];
2816  return_tensor *= inv_det_t;
2817 
2818  return return_tensor;
2819 }
2820 
2821 
2822 template <typename Number>
2824  invert(const Tensor<2, 3, Number> &t)
2825 {
2826  Tensor<2, 3, Number> return_tensor;
2827 
2828  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2829  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2830  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2831  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2832  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2833  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2834  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2835  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2836  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2837  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2838  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2839  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2840  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2841  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2842  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2843  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2844  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2845  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2846  const Number inv_det_t = internal::NumberType<Number>::value(
2847  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2848  t[0][2] * return_tensor[2][0]));
2849  return_tensor *= inv_det_t;
2850 
2851  return return_tensor;
2852 }
2853 
2854 #endif /* DOXYGEN */
2855 
2856 
2862 template <int dim, typename Number>
2865 {
2867  for (unsigned int i = 0; i < dim; ++i)
2868  {
2869  tt[i][i] = t[i][i];
2870  for (unsigned int j = i + 1; j < dim; ++j)
2871  {
2872  tt[i][j] = t[j][i];
2873  tt[j][i] = t[i][j];
2874  };
2875  }
2876  return tt;
2877 }
2878 
2879 
2893 template <int dim, typename Number>
2894 constexpr Tensor<2, dim, Number>
2896 {
2897  return determinant(t) * invert(t);
2898 }
2899 
2900 
2914 template <int dim, typename Number>
2915 constexpr Tensor<2, dim, Number>
2917 {
2918  return transpose(adjugate(t));
2919 }
2920 
2921 
2985 template <int dim, typename Number>
2988 
2989 
2997 template <int dim, typename Number>
2998 inline Number
3000 {
3002  for (unsigned int j = 0; j < dim; ++j)
3003  {
3005  for (unsigned int i = 0; i < dim; ++i)
3006  sum += std::fabs(t[i][j]);
3007 
3008  if (sum > max)
3009  max = sum;
3010  }
3011 
3012  return max;
3013 }
3014 
3015 
3023 template <int dim, typename Number>
3024 inline Number
3026 {
3028  for (unsigned int i = 0; i < dim; ++i)
3029  {
3031  for (unsigned int j = 0; j < dim; ++j)
3032  sum += std::fabs(t[i][j]);
3033 
3034  if (sum > max)
3035  max = sum;
3036  }
3037 
3038  return max;
3039 }
3040 
3042 
3043 
3044 #ifndef DOXYGEN
3045 
3046 
3047 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3048 
3049 // Specialization of functions for ADOL-C number types when
3050 // the advanced branching feature is used
3051 template <int dim>
3052 inline adouble
3054 {
3056  for (unsigned int j = 0; j < dim; ++j)
3057  {
3059  for (unsigned int i = 0; i < dim; ++i)
3060  sum += std::fabs(t[i][j]);
3061 
3062  condassign(max, (sum > max), sum, max);
3063  }
3064 
3065  return max;
3066 }
3067 
3068 
3069 template <int dim>
3070 inline adouble
3072 {
3074  for (unsigned int i = 0; i < dim; ++i)
3075  {
3077  for (unsigned int j = 0; j < dim; ++j)
3078  sum += std::fabs(t[i][j]);
3079 
3080  condassign(max, (sum > max), sum, max);
3081  }
3082 
3083  return max;
3084 }
3085 
3086 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3087 
3088 
3089 #endif // DOXYGEN
3090 
3092 
3093 #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:2777
Number * begin_raw()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1653
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2526
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2916
static constexpr const T & value(const T &t)
Definition: numbers.h:703
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:925
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:2864
Iterator unroll_recursion(const Iterator current, const Iterator end) const
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:587
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1956
constexpr Tensor operator-() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2625
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2999
#define AssertIndexRange(index, range)
Definition: exceptions.h:1718
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:2650
constexpr void clear()
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3025
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:869
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:547
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
static real_type abs(const number &x)
Definition: numbers.h:609
Definition: point.h:110
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
std::size_t size() const
Definition: array_view.h:575
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:1461
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2758
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:401
VectorType::value_type * end(VectorType &V)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:98
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:540
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:2447
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
Expression fabs(const Expression &x)
static constexpr unsigned int rank
Definition: tensor.h:526
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2895
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:2565
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:2017
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:2124
Definition: tensor.h:506
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:2033
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:941
#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:532
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:2236
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:2594
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:2372
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:2684