Reference documentation for deal.II version GIT f6a5d312c9 2023-10-04 08:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2023 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/kokkos.h>
23 #include <deal.II/base/numbers.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 
36 
37 // Forward declarations:
38 #ifndef DOXYGEN
39 template <typename ElementType, typename MemorySpace>
40 class ArrayView;
41 
42 template <int dim, typename Number>
43 DEAL_II_CXX20_REQUIRES(dim >= 0)
44 class Point;
45 
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_HOST_DEVICE
148 
156  template <typename OtherNumber>
157  constexpr DEAL_II_HOST_DEVICE
158  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
159 
165  template <typename OtherNumber>
166  constexpr DEAL_II_HOST_DEVICE
167  Tensor(const OtherNumber &initializer);
168 
169 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
173  constexpr DEAL_II_HOST_DEVICE
174  Tensor(const Tensor<0, dim, Number> &other);
175 
179  constexpr DEAL_II_HOST_DEVICE
180  Tensor(Tensor<0, dim, Number> &&other) noexcept;
181 #endif
182 
192  Number *
194 
204  const Number *
205  begin_raw() const;
206 
216  Number *
218 
229  const Number *
230  end_raw() const;
231 
241  constexpr DEAL_II_HOST_DEVICE
242  operator Number &();
243 
252  constexpr DEAL_II_HOST_DEVICE operator const Number &() const;
253 
261  template <typename OtherNumber>
262  constexpr DEAL_II_HOST_DEVICE Tensor &
264 
265 #if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
274  constexpr DEAL_II_HOST_DEVICE Tensor &
275  operator=(const Tensor<0, dim, Number> &rhs);
276 #endif
277 
278 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
283  operator=(Tensor<0, dim, Number> &&other) noexcept;
284 #endif
285 
292  template <typename OtherNumber>
293  constexpr DEAL_II_HOST_DEVICE Tensor &
294  operator=(const OtherNumber &d) &;
295 
301  template <typename OtherNumber>
302  constexpr DEAL_II_HOST_DEVICE Tensor &
303  operator=(const OtherNumber &d) && = delete;
304 
308  template <typename OtherNumber>
309  constexpr bool
311 
315  template <typename OtherNumber>
316  constexpr bool
318 
324  template <typename OtherNumber>
325  constexpr DEAL_II_HOST_DEVICE Tensor &
327 
333  template <typename OtherNumber>
334  constexpr DEAL_II_HOST_DEVICE Tensor &
336 
342  template <typename OtherNumber>
343  constexpr DEAL_II_HOST_DEVICE Tensor &
344  operator*=(const OtherNumber &factor);
345 
351  template <typename OtherNumber>
352  constexpr DEAL_II_HOST_DEVICE Tensor &
353  operator/=(const OtherNumber &factor);
354 
360  constexpr DEAL_II_HOST_DEVICE Tensor
361  operator-() const;
362 
375  constexpr void
376  clear();
377 
383  real_type
384  norm() const;
385 
393  norm_square() const;
394 
402  template <class Iterator>
403  void
404  unroll(const Iterator begin, const Iterator end) const;
405 
411  template <class Archive>
412  void
413  serialize(Archive &ar, const unsigned int version);
414 
419  using tensor_type = Number;
420 
421 private:
425  Number value;
426 
430  template <typename Iterator>
431  Iterator
432  unroll_recursion(const Iterator current, const Iterator end) const;
433 
434  // Allow an arbitrary Tensor to access the underlying values.
435  template <int, int, typename>
436  friend class Tensor;
437 };
438 
439 
440 
514 template <int rank_, int dim, typename Number>
515 class Tensor
516 {
517 public:
518  static_assert(rank_ >= 1,
519  "Tensors must have a rank greater than or equal to one.");
520  static_assert(dim >= 0,
521  "Tensors must have a dimension greater than or equal to zero.");
530  static constexpr unsigned int dimension = dim;
531 
535  static constexpr unsigned int rank = rank_;
536 
541  static constexpr unsigned int n_independent_components =
542  Tensor<rank_ - 1, dim>::n_independent_components * dim;
543 
549  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
550 
555  using array_type =
556  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
557 
565 
571  constexpr DEAL_II_HOST_DEVICE explicit Tensor(const array_type &initializer);
572 
585  template <typename ElementType, typename MemorySpace>
586  constexpr DEAL_II_HOST_DEVICE explicit Tensor(
587  const ArrayView<ElementType, MemorySpace> &initializer);
588 
596  template <typename OtherNumber>
597  constexpr DEAL_II_HOST_DEVICE
599 
603  template <typename OtherNumber>
604  constexpr Tensor(
605  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
606 
610  template <typename OtherNumber>
611  constexpr
612  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
613 
614 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
618  constexpr Tensor(const Tensor<rank_, dim, Number> &);
619 
623  constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
624 #endif
625 
631  constexpr DEAL_II_HOST_DEVICE value_type &
632  operator[](const unsigned int i);
633 
639  constexpr DEAL_II_HOST_DEVICE const value_type &
640  operator[](const unsigned int i) const;
641 
645  constexpr const Number &
646  operator[](const TableIndices<rank_> &indices) const;
647 
651  constexpr Number &
653 
657  Number *
659 
663  const Number *
664  begin_raw() const;
665 
669  Number *
671 
675  const Number *
676  end_raw() const;
677 
685  template <typename OtherNumber>
686  constexpr DEAL_II_HOST_DEVICE Tensor &
688 
695  constexpr DEAL_II_HOST_DEVICE Tensor &
696  operator=(const Number &d) &;
697 
703  constexpr DEAL_II_HOST_DEVICE Tensor &
704  operator=(const Number &d) && = delete;
705 
706 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
710  constexpr Tensor<rank_, dim, Number> &
712 
716  constexpr Tensor<rank_, dim, Number> &
718 #endif
719 
723  template <typename OtherNumber>
724  constexpr bool
726 
730  template <typename OtherNumber>
731  constexpr bool
733 
739  template <typename OtherNumber>
740  constexpr DEAL_II_HOST_DEVICE Tensor &
742 
748  template <typename OtherNumber>
749  constexpr DEAL_II_HOST_DEVICE Tensor &
751 
758  template <typename OtherNumber>
759  constexpr DEAL_II_HOST_DEVICE Tensor &
760  operator*=(const OtherNumber &factor);
761 
767  template <typename OtherNumber>
768  constexpr DEAL_II_HOST_DEVICE Tensor &
769  operator/=(const OtherNumber &factor);
770 
776  constexpr DEAL_II_HOST_DEVICE Tensor
777  operator-() const;
778 
791  constexpr void
792  clear();
793 
803  norm() const;
804 
811  constexpr DEAL_II_HOST_DEVICE
813  norm_square() const;
814 
825  template <typename OtherNumber>
826  DEAL_II_DEPRECATED void
827  unroll(Vector<OtherNumber> &result) const;
828 
839  template <class Iterator>
840  void
841  unroll(const Iterator begin, const Iterator end) const;
842 
847  static constexpr DEAL_II_HOST_DEVICE unsigned int
849 
856  unrolled_to_component_indices(const unsigned int i);
857 
862  static constexpr std::size_t
864 
870  template <class Archive>
871  void
872  serialize(Archive &ar, const unsigned int version);
873 
879 
880 private:
884  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
885  // ... avoid a compiler warning in case of dim == 0 and ensure that the
886  // array always has positive size.
887 
891  template <typename Iterator>
892  Iterator
893  unroll_recursion(const Iterator current, const Iterator end) const;
894 
901  template <typename ArrayLike, std::size_t... Indices>
902  constexpr DEAL_II_HOST_DEVICE
903  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
904 
905  // Allow an arbitrary Tensor to access the underlying values.
906  template <int, int, typename>
907  friend class Tensor;
908 
909  // Point is allowed access to the coordinates. This is supposed to improve
910  // speed.
911  friend class Point<dim, Number>;
912 };
913 
914 
915 #ifndef DOXYGEN
916 namespace internal
917 {
918  // Workaround: The following 4 overloads are necessary to be able to
919  // compile the library with Apple Clang 8 and older. We should remove
920  // these overloads again when we bump the minimal required version to
921  // something later than clang-3.6 / Apple Clang 6.3.
922  template <int rank, int dim, typename T, typename U>
923  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
924  {
925  using type =
927  };
928 
929  template <int rank, int dim, typename T, typename U>
930  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
931  {
932  using type =
934  };
935 
936  template <typename T, int rank, int dim, typename U>
937  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
938  {
939  using type =
941  };
942 
943  template <int rank, int dim, typename T, typename U>
944  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
945  {
946  using type =
948  };
949  // end workaround
950 
955  template <int rank, int dim, typename T>
956  struct NumberType<Tensor<rank, dim, T>>
957  {
958  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const
960  value(const Tensor<rank, dim, T> &t)
961  {
962  return t;
963  }
964 
966  value(const T &t)
967  {
969  tmp = t;
970  return tmp;
971  }
972  };
973 } // namespace internal
974 
975 
976 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
977 
978 
979 template <int dim, typename Number>
982  // Some auto-differentiable numbers need explicit
983  // zero initialization such as adtl::adouble.
984  : Tensor{0.0}
985 {}
986 
987 
988 
989 template <int dim, typename Number>
990 template <typename OtherNumber>
992 Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
993  : value(internal::NumberType<Number>::value(initializer))
994 {}
995 
996 
997 
998 template <int dim, typename Number>
999 template <typename OtherNumber>
1002  : Tensor{p.value}
1003 {}
1004 
1005 
1006 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1007 template <int dim, typename Number>
1010  : value{other.value}
1011 {}
1012 
1013 
1014 
1015 template <int dim, typename Number>
1018  : value{std::move(other.value)}
1019 {}
1020 # endif
1021 
1022 
1023 template <int dim, typename Number>
1024 inline Number *
1026 {
1027  return std::addressof(value);
1028 }
1029 
1030 
1031 
1032 template <int dim, typename Number>
1033 inline const Number *
1035 {
1036  return std::addressof(value);
1037 }
1038 
1039 
1040 
1041 template <int dim, typename Number>
1042 inline Number *
1044 {
1046 }
1047 
1048 
1049 
1050 template <int dim, typename Number>
1051 const Number *
1053 {
1055 }
1056 
1057 
1058 
1059 template <int dim, typename Number>
1062 {
1063  Assert(dim != 0,
1064  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1065  return value;
1066 }
1067 
1068 
1069 template <int dim, typename Number>
1070 constexpr inline DEAL_II_ALWAYS_INLINE
1072 {
1073  Assert(dim != 0,
1074  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1075  return value;
1076 }
1077 
1078 
1079 
1080 template <int dim, typename Number>
1081 template <typename OtherNumber>
1084 {
1086  return *this;
1087 }
1088 
1089 
1090 # if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1091 template <int dim, typename Number>
1094 {
1095  value = p.value;
1096  return *this;
1097 }
1098 # endif
1099 
1100 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1101 template <int dim, typename Number>
1104 {
1105  value = std::move(other.value);
1106  return *this;
1107 }
1108 # endif
1109 
1110 
1111 
1112 template <int dim, typename Number>
1113 template <typename OtherNumber>
1115 Tensor<0, dim, Number>::operator=(const OtherNumber &d) &
1116 {
1118  return *this;
1119 }
1120 
1121 
1122 template <int dim, typename Number>
1123 template <typename OtherNumber>
1124 constexpr inline bool
1126 {
1127 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1128  Assert(!(std::is_same_v<Number, adouble> ||
1129  std::is_same_v<OtherNumber, adouble>),
1130  ExcMessage(
1131  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1132  "been extended to support advanced branching."));
1133 # endif
1134 
1135  return numbers::values_are_equal(value, p.value);
1136 }
1137 
1138 
1139 template <int dim, typename Number>
1140 template <typename OtherNumber>
1141 constexpr bool
1143 {
1144  return !((*this) == p);
1145 }
1146 
1147 
1148 template <int dim, typename Number>
1149 template <typename OtherNumber>
1152 {
1153  value += p.value;
1154  return *this;
1155 }
1156 
1157 
1158 template <int dim, typename Number>
1159 template <typename OtherNumber>
1162 {
1163  value -= p.value;
1164  return *this;
1165 }
1166 
1167 
1168 
1169 namespace internal
1170 {
1171  namespace ComplexWorkaround
1172  {
1173  template <typename Number, typename OtherNumber>
1174  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1175  multiply_assign_scalar(Number &val, const OtherNumber &s)
1176  {
1177  val *= s;
1178  }
1179 
1180  template <typename Number, typename OtherNumber>
1181  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1182  multiply_assign_scalar(std::complex<Number> &val, const OtherNumber &s)
1183  {
1184 # if KOKKOS_VERSION >= 30600
1185  KOKKOS_IF_ON_HOST((val *= s;))
1186  KOKKOS_IF_ON_DEVICE(({
1187  (void)val;
1188  (void)s;
1189  Kokkos::abort(
1190  "This function is not implemented for std::complex<Number>!\n");
1191  }))
1192 # else
1193 # ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1194  val *= s;
1195 # else
1196  (void)val;
1197  (void)s;
1198  Kokkos::abort(
1199  "This function is not implemented for std::complex<Number>!\n");
1200 # endif
1201 # endif
1202  }
1203  } // namespace ComplexWorkaround
1204 } // namespace internal
1205 
1206 
1207 template <int dim, typename Number>
1208 template <typename OtherNumber>
1210 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1211 {
1212  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1213  return *this;
1214 }
1215 
1216 
1217 
1218 template <int dim, typename Number>
1219 template <typename OtherNumber>
1221 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1222 {
1223  value /= s;
1224  return *this;
1225 }
1226 
1227 
1228 template <int dim, typename Number>
1231 {
1232  return -value;
1233 }
1234 
1235 
1236 template <int dim, typename Number>
1239 {
1240  Assert(dim != 0,
1241  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1242  return numbers::NumberTraits<Number>::abs(value);
1243 }
1244 
1245 
1246 template <int dim, typename Number>
1250 {
1251  Assert(dim != 0,
1252  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1254 }
1255 
1256 
1257 
1258 template <int dim, typename Number>
1259 template <typename Iterator>
1260 Iterator
1261 Tensor<0, dim, Number>::unroll_recursion(const Iterator current,
1262  const Iterator end) const
1263 {
1264  (void)end;
1265  Assert(dim != 0,
1266  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1267  Assert(std::distance(current, end) >= 1,
1268  ExcMessage("The provided iterator range must contain at least one "
1269  "element."));
1270  *current = value;
1271  return std::next(current);
1272 }
1273 
1274 
1275 
1276 template <int dim, typename Number>
1277 constexpr inline void
1279 {
1280  // Some auto-differentiable numbers need explicit
1281  // zero initialization.
1283 }
1284 
1285 
1286 
1287 template <int dim, typename Number>
1288 template <class Iterator>
1289 inline void
1290 Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1291 {
1294 }
1295 
1296 
1297 
1298 template <int dim, typename Number>
1299 template <class Archive>
1300 inline void
1301 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1302 {
1303  ar &value;
1304 }
1305 
1306 
1307 template <int dim, typename Number>
1309 
1310 
1311 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1312 
1313 template <int rank_, int dim, typename Number>
1314 template <typename ArrayLike, std::size_t... indices>
1316 Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1317  std::index_sequence<indices...>)
1318  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1319 {
1320  static_assert(sizeof...(indices) == dim,
1321  "dim should match the number of indices");
1322 }
1323 
1324 
1325 
1326 template <int rank_, int dim, typename Number>
1329  // We would like to use =default, but this causes compile errors with some
1330  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1331  : values{}
1332 {}
1333 
1334 
1335 
1336 template <int rank_, int dim, typename Number>
1338 Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1339  : Tensor(initializer, std::make_index_sequence<dim>{})
1340 {}
1341 
1342 
1343 
1344 template <int rank_, int dim, typename Number>
1345 template <typename ElementType, typename MemorySpace>
1348  const ArrayView<ElementType, MemorySpace> &initializer)
1349 {
1350  // make nvcc happy
1351  const int my_n_independent_components = n_independent_components;
1352  AssertDimension(initializer.size(), my_n_independent_components);
1353 
1354  for (unsigned int i = 0; i < my_n_independent_components; ++i)
1355  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1356 }
1357 
1358 
1359 
1360 template <int rank_, int dim, typename Number>
1361 template <typename OtherNumber>
1364  const Tensor<rank_, dim, OtherNumber> &initializer)
1365  : Tensor(initializer, std::make_index_sequence<dim>{})
1366 {}
1367 
1368 
1369 
1370 template <int rank_, int dim, typename Number>
1371 template <typename OtherNumber>
1372 constexpr DEAL_II_ALWAYS_INLINE
1374  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1375  : Tensor(initializer, std::make_index_sequence<dim>{})
1376 {}
1377 
1378 
1379 
1380 template <int rank_, int dim, typename Number>
1381 template <typename OtherNumber>
1383 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1384 {
1385  return Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>(values);
1386 }
1387 
1388 
1389 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1390 template <int rank_, int dim, typename Number>
1391 constexpr DEAL_II_ALWAYS_INLINE
1393 {
1394  for (unsigned int i = 0; i < dim; ++i)
1395  values[i] = other.values[i];
1396 }
1397 
1398 
1399 
1400 template <int rank_, int dim, typename Number>
1401 constexpr DEAL_II_ALWAYS_INLINE
1403 {
1404  for (unsigned int i = 0; i < dim; ++i)
1405  values[i] = other.values[i];
1406 }
1407 # endif
1408 
1409 namespace internal
1410 {
1411  namespace TensorSubscriptor
1412  {
1413  template <typename ArrayElementType, int dim>
1414  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
1415  subscript(ArrayElementType *values,
1416  const unsigned int i,
1417  std::integral_constant<int, dim>)
1418  {
1419  AssertIndexRange(i, dim);
1420  return values[i];
1421  }
1422 
1423  template <typename ArrayElementType>
1424  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
1425  subscript(ArrayElementType *dummy,
1426  const unsigned int,
1427  std::integral_constant<int, 0>)
1428  {
1429  Assert(
1430  false,
1431  ExcMessage(
1432  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1433  return *dummy;
1434  }
1435  } // namespace TensorSubscriptor
1436 } // namespace internal
1437 
1438 
1439 template <int rank_, int dim, typename Number>
1442  Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1443 {
1444  return ::internal::TensorSubscriptor::subscript(
1445  values, i, std::integral_constant<int, dim>());
1446 }
1447 
1448 
1449 template <int rank_, int dim, typename Number>
1450 constexpr DEAL_II_ALWAYS_INLINE
1452  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1453 {
1454 # if KOKKOS_VERSION < 30700
1455 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1456  AssertIndexRange(i, dim);
1457 # endif
1458 # else
1459  KOKKOS_IF_ON_HOST((AssertIndexRange(i, dim);))
1460 # endif
1461 
1462  return values[i];
1463 }
1464 
1465 
1466 template <int rank_, int dim, typename Number>
1467 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1469 {
1470 # if KOKKOS_VERSION < 30700
1471 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1472  Assert(dim != 0,
1473  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1474 # endif
1475 # else
1476  KOKKOS_IF_ON_HOST(
1477  (Assert(dim != 0,
1478  ExcMessage(
1479  "Cannot access an object of type Tensor<rank_,0,Number>"));))
1480 # endif
1481 
1482  return TensorAccessors::extract<rank_>(*this, indices);
1483 }
1484 
1485 
1486 
1487 template <int rank_, int dim, typename Number>
1488 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1490 {
1491 # if KOKKOS_VERSION < 30700
1492 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1493  Assert(dim != 0,
1494  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1495 # endif
1496 # else
1497  KOKKOS_IF_ON_HOST(
1498  (Assert(dim != 0,
1499  ExcMessage(
1500  "Cannot access an object of type Tensor<rank_,0,Number>"));))
1501 # endif
1502 
1503  return TensorAccessors::extract<rank_>(*this, indices);
1504 }
1505 
1506 
1507 
1508 template <int rank_, int dim, typename Number>
1509 inline Number *
1511 {
1512  return std::addressof(
1513  this->operator[](this->unrolled_to_component_indices(0)));
1514 }
1515 
1516 
1517 
1518 template <int rank_, int dim, typename Number>
1519 inline const Number *
1521 {
1522  return std::addressof(
1523  this->operator[](this->unrolled_to_component_indices(0)));
1524 }
1525 
1526 
1527 
1528 template <int rank_, int dim, typename Number>
1529 inline Number *
1531 {
1532  return begin_raw() + n_independent_components;
1533 }
1534 
1535 
1536 
1537 template <int rank_, int dim, typename Number>
1538 inline const Number *
1540 {
1541  return begin_raw() + n_independent_components;
1542 }
1543 
1544 
1545 
1546 template <int rank_, int dim, typename Number>
1547 template <typename OtherNumber>
1550 {
1551  // The following loop could be written more concisely using std::copy, but
1552  // that function is only constexpr from C++20 on.
1553  for (unsigned int i = 0; i < dim; ++i)
1554  values[i] = t.values[i];
1555  return *this;
1556 }
1557 
1558 
1559 
1560 template <int rank_, int dim, typename Number>
1563  Tensor<rank_, dim, Number>::operator=(const Number &d) &
1564 {
1566  (void)d;
1567 
1568  for (unsigned int i = 0; i < dim; ++i)
1570  return *this;
1571 }
1572 
1573 
1574 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1575 template <int rank_, int dim, typename Number>
1578 {
1579  for (unsigned int i = 0; i < dim; ++i)
1580  values[i] = other.values[i];
1581  return *this;
1582 }
1583 
1584 
1585 
1586 template <int rank_, int dim, typename Number>
1589  Tensor<rank_, dim, Number> &&other) noexcept
1590 {
1591  for (unsigned int i = 0; i < dim; ++i)
1592  values[i] = other.values[i];
1593  return *this;
1594 }
1595 # endif
1596 
1597 
1598 template <int rank_, int dim, typename Number>
1599 template <typename OtherNumber>
1600 constexpr inline bool
1602  const Tensor<rank_, dim, OtherNumber> &p) const
1603 {
1604  for (unsigned int i = 0; i < dim; ++i)
1605  if (values[i] != p.values[i])
1606  return false;
1607  return true;
1608 }
1609 
1610 
1611 // At some places in the library, we have Point<0> for formal reasons
1612 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1613 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1614 // in the above function that the loop end check always fails, we
1615 // implement this function here
1616 template <>
1617 template <>
1618 constexpr inline bool
1620 {
1621  return true;
1622 }
1623 
1624 
1625 template <int rank_, int dim, typename Number>
1626 template <typename OtherNumber>
1627 constexpr bool
1629  const Tensor<rank_, dim, OtherNumber> &p) const
1630 {
1631  return !((*this) == p);
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
1654 {
1655  for (unsigned int i = 0; i < dim; ++i)
1656  values[i] -= p.values[i];
1657  return *this;
1658 }
1659 
1660 
1661 template <int rank_, int dim, typename Number>
1662 template <typename OtherNumber>
1663 constexpr inline DEAL_II_ALWAYS_INLINE
1665  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1666 {
1667  for (unsigned int i = 0; i < dim; ++i)
1668  values[i] *= s;
1669  return *this;
1670 }
1671 
1672 
1673 namespace internal
1674 {
1675  namespace TensorImplementation
1676  {
1677  template <int rank,
1678  int dim,
1679  typename Number,
1680  typename OtherNumber,
1681  std::enable_if_t<
1682  !std::is_integral<
1683  typename ProductType<Number, OtherNumber>::type>::value &&
1684  !std::is_same_v<Number, Differentiation::SD::Expression>,
1685  int> = 0>
1686  constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
1688  const OtherNumber &factor)
1689  {
1690  const Number inverse_factor = Number(1.) / factor;
1691  // recurse over the base objects
1692  for (unsigned int d = 0; d < dim; ++d)
1693  t[d] *= inverse_factor;
1694  }
1695 
1696 
1697  template <int rank,
1698  int dim,
1699  typename Number,
1700  typename OtherNumber,
1701  std::enable_if_t<
1702  std::is_integral<
1703  typename ProductType<Number, OtherNumber>::type>::value ||
1704  std::is_same_v<Number, Differentiation::SD::Expression>,
1705  int> = 0>
1706  constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
1708  const OtherNumber &factor)
1709  {
1710  // recurse over the base objects
1711  for (unsigned int d = 0; d < dim; ++d)
1712  t[d] /= factor;
1713  }
1714  } // namespace TensorImplementation
1715 } // namespace internal
1716 
1717 
1718 template <int rank_, int dim, typename Number>
1719 template <typename OtherNumber>
1720 constexpr inline DEAL_II_ALWAYS_INLINE
1722  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1723 {
1725  return *this;
1726 }
1727 
1728 
1729 template <int rank_, int dim, typename Number>
1730 constexpr inline DEAL_II_ALWAYS_INLINE
1733 {
1735 
1736  for (unsigned int i = 0; i < dim; ++i)
1737  tmp.values[i] = -values[i];
1738 
1739  return tmp;
1740 }
1741 
1742 
1743 template <int rank_, int dim, typename Number>
1746 {
1747  // Make things work with AD types
1748  using std::sqrt;
1749  return sqrt(norm_square());
1750 }
1751 
1752 
1753 template <int rank_, int dim, typename Number>
1757 {
1759  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1760  for (unsigned int i = 0; i < dim; ++i)
1761  s += values[i].norm_square();
1762 
1763  return s;
1764 }
1765 
1766 
1767 
1768 template <int rank_, int dim, typename Number>
1769 template <typename OtherNumber>
1770 inline void
1772 {
1773  unroll(result.begin(), result.end());
1774 }
1775 
1776 
1777 
1778 template <int rank_, int dim, typename Number>
1779 template <class Iterator>
1780 inline void
1782  const Iterator end) const
1783 {
1784  AssertDimension(std::distance(begin, end), n_independent_components);
1785  unroll_recursion(begin, end);
1786 }
1787 
1788 
1789 
1790 template <int rank_, int dim, typename Number>
1791 template <typename Iterator>
1792 Iterator
1793 Tensor<rank_, dim, Number>::unroll_recursion(const Iterator current,
1794  const Iterator end) const
1795 {
1796  Iterator next = current;
1797  for (unsigned int i = 0; i < dim; ++i)
1798  next = values[i].unroll_recursion(next, end);
1799  return next;
1800 }
1801 
1802 
1803 template <int rank_, int dim, typename Number>
1804 constexpr inline unsigned int
1806  const TableIndices<rank_> &indices)
1807 {
1808  unsigned int index = 0;
1809  for (int r = 0; r < rank_; ++r)
1810  index = index * dim + indices[r];
1811 
1812  return index;
1813 }
1814 
1815 
1816 
1817 namespace internal
1818 {
1819  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1820  // and rank=2. Make sure we don't have compiler warnings.
1821 
1822  template <int dim>
1823  DEAL_II_HOST_DEVICE inline constexpr unsigned int
1824  mod(const unsigned int x)
1825  {
1826  return x % dim;
1827  }
1828 
1829  template <>
1830  DEAL_II_HOST_DEVICE inline unsigned int
1831  mod<0>(const unsigned int x)
1832  {
1833  Assert(false, ExcInternalError());
1834  return x;
1835  }
1836 
1837  template <int dim>
1838  DEAL_II_HOST_DEVICE inline constexpr unsigned int
1839  div(const unsigned int x)
1840  {
1841  return x / dim;
1842  }
1843 
1844  template <>
1845  DEAL_II_HOST_DEVICE inline unsigned int
1846  div<0>(const unsigned int x)
1847  {
1848  Assert(false, ExcInternalError());
1849  return x;
1850  }
1851 
1852 } // namespace internal
1853 
1854 
1855 
1856 template <int rank_, int dim, typename Number>
1857 constexpr inline TableIndices<rank_>
1859 {
1860  // Work-around nvcc warning
1861  unsigned int dummy = n_independent_components;
1862  AssertIndexRange(i, dummy);
1863  (void)dummy;
1864 
1865  TableIndices<rank_> indices;
1866 
1867  unsigned int remainder = i;
1868  for (int r = rank_ - 1; r >= 0; --r)
1869  {
1870  indices[r] = internal::mod<dim>(remainder);
1871  remainder = internal::div<dim>(remainder);
1872  }
1873  Assert(remainder == 0, ExcInternalError());
1874 
1875  return indices;
1876 }
1877 
1878 
1879 template <int rank_, int dim, typename Number>
1880 constexpr inline void
1882 {
1883  for (unsigned int i = 0; i < dim; ++i)
1885 }
1886 
1887 
1888 template <int rank_, int dim, typename Number>
1889 constexpr std::size_t
1891 {
1892  return sizeof(Tensor<rank_, dim, Number>);
1893 }
1894 
1895 
1896 template <int rank_, int dim, typename Number>
1897 template <class Archive>
1898 inline void
1899 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1900 {
1901  ar &values;
1902 }
1903 
1904 
1905 template <int rank_, int dim, typename Number>
1907 
1908 #endif // DOXYGEN
1909 
1910 /* ----------------- Non-member functions operating on tensors. ------------ */
1911 
1924 template <int rank_, int dim, typename Number>
1925 inline std::ostream &
1926 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1927 {
1928  for (unsigned int i = 0; i < dim; ++i)
1929  {
1930  out << p[i];
1931  if (i != dim - 1)
1932  out << ' ';
1933  }
1934 
1935  return out;
1936 }
1937 
1938 
1945 template <int dim, typename Number>
1946 inline std::ostream &
1947 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1948 {
1949  out << static_cast<const Number &>(p);
1950  return out;
1951 }
1952 
1953 
1972 template <int dim, typename Number, typename Other>
1975  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1976 {
1977  return object * static_cast<const Number &>(t);
1978 }
1979 
1980 
1981 
1992 template <int dim, typename Number, typename Other>
1995  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1996 {
1997  return static_cast<const Number &>(t) * object;
1998 }
1999 
2000 
2012 template <int dim, typename Number, typename OtherNumber>
2016  const Tensor<0, dim, OtherNumber> &src2)
2017 {
2018  return static_cast<const Number &>(src1) *
2019  static_cast<const OtherNumber &>(src2);
2020 }
2021 
2022 
2030 template <int dim, typename Number, typename OtherNumber>
2032  Tensor<0,
2033  dim,
2034  typename ProductType<Number,
2035  typename EnableIfScalar<OtherNumber>::type>::type>
2036  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2037 {
2038  return static_cast<const Number &>(t) / factor;
2039 }
2040 
2041 
2049 template <int dim, typename Number, typename OtherNumber>
2053  const Tensor<0, dim, OtherNumber> &q)
2054 {
2055  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2056 }
2057 
2058 
2066 template <int dim, typename Number, typename OtherNumber>
2070  const Tensor<0, dim, OtherNumber> &q)
2071 {
2072  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2073 }
2074 
2075 
2088 template <int rank, int dim, typename Number, typename OtherNumber>
2090  Tensor<rank,
2091  dim,
2092  typename ProductType<Number,
2093  typename EnableIfScalar<OtherNumber>::type>::type>
2094  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2095 {
2096  // recurse over the base objects
2098  for (unsigned int d = 0; d < dim; ++d)
2099  tt[d] = t[d] * factor;
2100  return tt;
2101 }
2102 
2103 
2116 template <int rank, int dim, typename Number, typename OtherNumber>
2118  Tensor<rank,
2119  dim,
2121  OtherNumber>::type>
2122  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2123 {
2124  // simply forward to the operator above
2125  return t * factor;
2126 }
2127 
2128 
2129 namespace internal
2130 {
2131  namespace TensorImplementation
2132  {
2133  template <int rank,
2134  int dim,
2135  typename Number,
2136  typename OtherNumber,
2137  std::enable_if_t<
2138  !std::is_integral<
2139  typename ProductType<Number, OtherNumber>::type>::value,
2140  int> = 0>
2144  const OtherNumber &factor)
2145  {
2147  const Number inverse_factor = Number(1.) / factor;
2148  // recurse over the base objects
2149  for (unsigned int d = 0; d < dim; ++d)
2150  tt[d] = t[d] * inverse_factor;
2151  return tt;
2152  }
2153 
2154 
2155  template <int rank,
2156  int dim,
2157  typename Number,
2158  typename OtherNumber,
2159  std::enable_if_t<
2160  std::is_integral<
2161  typename ProductType<Number, OtherNumber>::type>::value,
2162  int> = 0>
2166  const OtherNumber &factor)
2167  {
2169  // recurse over the base objects
2170  for (unsigned int d = 0; d < dim; ++d)
2171  tt[d] = t[d] / factor;
2172  return tt;
2173  }
2174  } // namespace TensorImplementation
2175 } // namespace internal
2176 
2177 
2187 template <int rank, int dim, typename Number, typename OtherNumber>
2189  Tensor<rank,
2190  dim,
2191  typename ProductType<Number,
2192  typename EnableIfScalar<OtherNumber>::type>::type>
2193  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2194 {
2196 }
2197 
2198 
2208 template <int rank, int dim, typename Number, typename OtherNumber>
2213 {
2215 
2216  for (unsigned int i = 0; i < dim; ++i)
2217  tmp[i] += q[i];
2218 
2219  return tmp;
2220 }
2221 
2222 
2232 template <int rank, int dim, typename Number, typename OtherNumber>
2237 {
2239 
2240  for (unsigned int i = 0; i < dim; ++i)
2241  tmp[i] -= q[i];
2242 
2243  return tmp;
2244 }
2245 
2252 template <int dim, typename Number, typename OtherNumber>
2253 inline constexpr DEAL_II_ALWAYS_INLINE
2256  const Tensor<0, dim, OtherNumber> &src2)
2257 {
2259 
2260  tmp *= src2;
2261 
2262  return tmp;
2263 }
2264 
2281 template <int rank, int dim, typename Number, typename OtherNumber>
2282 inline constexpr DEAL_II_ALWAYS_INLINE
2285  const Tensor<rank, dim, OtherNumber> &src2)
2286 {
2288 
2289  for (unsigned int i = 0; i < dim; ++i)
2290  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2292 
2293  return tmp;
2294 }
2295 
2340 template <int rank_1,
2341  int rank_2,
2342  int dim,
2343  typename Number,
2344  typename OtherNumber,
2345  typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2346 constexpr inline DEAL_II_ALWAYS_INLINE
2347  typename Tensor<rank_1 + rank_2 - 2,
2348  dim,
2349  typename ProductType<Number, OtherNumber>::type>::tensor_type
2352 {
2353  typename Tensor<rank_1 + rank_2 - 2,
2354  dim,
2355  typename ProductType<Number, OtherNumber>::type>::tensor_type
2356  result{};
2357 
2358  TensorAccessors::internal::
2359  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2360  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2361  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2362 
2363  return result;
2364 }
2365 
2366 
2395 template <int index_1,
2396  int index_2,
2397  int rank_1,
2398  int rank_2,
2399  int dim,
2400  typename Number,
2401  typename OtherNumber>
2402 constexpr inline DEAL_II_ALWAYS_INLINE
2403  typename Tensor<rank_1 + rank_2 - 2,
2404  dim,
2405  typename ProductType<Number, OtherNumber>::type>::tensor_type
2408 {
2409  Assert(0 <= index_1 && index_1 < rank_1,
2410  ExcMessage(
2411  "The specified index_1 must lie within the range [0,rank_1)"));
2412  Assert(0 <= index_2 && index_2 < rank_2,
2413  ExcMessage(
2414  "The specified index_2 must lie within the range [0,rank_2)"));
2415 
2416  using namespace TensorAccessors;
2417  using namespace TensorAccessors::internal;
2418 
2419  // Reorder index_1 to the end of src1:
2421  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2422 
2423  // Reorder index_2 to the end of src2:
2424  const ReorderedIndexView<index_2,
2425  rank_2,
2427  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2428 
2429  typename Tensor<rank_1 + rank_2 - 2,
2430  dim,
2431  typename ProductType<Number, OtherNumber>::type>::tensor_type
2432  result{};
2433  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2434  return result;
2435 }
2436 
2437 
2468 template <int index_1,
2469  int index_2,
2470  int index_3,
2471  int index_4,
2472  int rank_1,
2473  int rank_2,
2474  int dim,
2475  typename Number,
2476  typename OtherNumber>
2477 constexpr inline
2478  typename Tensor<rank_1 + rank_2 - 4,
2479  dim,
2480  typename ProductType<Number, OtherNumber>::type>::tensor_type
2481  double_contract(const Tensor<rank_1, dim, Number> &src1,
2483 {
2484  Assert(0 <= index_1 && index_1 < rank_1,
2485  ExcMessage(
2486  "The specified index_1 must lie within the range [0,rank_1)"));
2487  Assert(0 <= index_3 && index_3 < rank_1,
2488  ExcMessage(
2489  "The specified index_3 must lie within the range [0,rank_1)"));
2490  Assert(index_1 != index_3,
2491  ExcMessage("index_1 and index_3 must not be the same"));
2492  Assert(0 <= index_2 && index_2 < rank_2,
2493  ExcMessage(
2494  "The specified index_2 must lie within the range [0,rank_2)"));
2495  Assert(0 <= index_4 && index_4 < rank_2,
2496  ExcMessage(
2497  "The specified index_4 must lie within the range [0,rank_2)"));
2498  Assert(index_2 != index_4,
2499  ExcMessage("index_2 and index_4 must not be the same"));
2500 
2501  using namespace TensorAccessors;
2502  using namespace TensorAccessors::internal;
2503 
2504  // Reorder index_1 to the end of src1:
2506  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2507 
2508  // Reorder index_2 to the end of src2:
2510  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2511 
2512  // Now, reorder index_3 to the end of src1. We have to make sure to
2513  // preserve the original ordering: index_1 has been removed. If
2514  // index_3 > index_1, we have to use (index_3 - 1) instead:
2516  (index_3 < index_1 ? index_3 : index_3 - 1),
2517  rank_1,
2518  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2519  reord_3 =
2520  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2521  index_3 - 1,
2522  rank_1 > (reord_1);
2523 
2524  // Now, reorder index_4 to the end of src2. We have to make sure to
2525  // preserve the original ordering: index_2 has been removed. If
2526  // index_4 > index_2, we have to use (index_4 - 1) instead:
2528  (index_4 < index_2 ? index_4 : index_4 - 1),
2529  rank_2,
2530  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2531  reord_4 =
2532  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2533  index_4 - 1,
2534  rank_2 > (reord_2);
2535 
2536  typename Tensor<rank_1 + rank_2 - 4,
2537  dim,
2538  typename ProductType<Number, OtherNumber>::type>::tensor_type
2539  result{};
2540  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2541  return result;
2542 }
2543 
2544 
2557 template <int rank, int dim, typename Number, typename OtherNumber>
2558 constexpr inline DEAL_II_ALWAYS_INLINE
2560  scalar_product(const Tensor<rank, dim, Number> &left,
2561  const Tensor<rank, dim, OtherNumber> &right)
2562 {
2563  typename ProductType<Number, OtherNumber>::type result{};
2564  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2565  return result;
2566 }
2567 
2568 
2586 template <template <int, int, typename> class TensorT1,
2587  template <int, int, typename>
2588  class TensorT2,
2589  template <int, int, typename>
2590  class TensorT3,
2591  int rank_1,
2592  int rank_2,
2593  int dim,
2594  typename T1,
2595  typename T2,
2596  typename T3>
2597 constexpr inline DEAL_II_ALWAYS_INLINE
2599  contract3(const TensorT1<rank_1, dim, T1> &left,
2600  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2601  const TensorT3<rank_2, dim, T3> &right)
2602 {
2603  using return_type =
2605  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2606  middle,
2607  right);
2608 }
2609 
2610 
2621 template <int rank_1,
2622  int rank_2,
2623  int dim,
2624  typename Number,
2625  typename OtherNumber>
2626 constexpr inline DEAL_II_ALWAYS_INLINE
2630 {
2631  typename Tensor<rank_1 + rank_2,
2632  dim,
2633  typename ProductType<Number, OtherNumber>::type>::tensor_type
2634  result{};
2635  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2636  return result;
2637 }
2638 
2639 
2658 template <int dim, typename Number>
2660 cross_product_2d(const Tensor<1, dim, Number> &src)
2661 {
2662  Assert(dim == 2, ExcInternalError());
2663 
2664  Tensor<1, dim, Number> result;
2665 
2666  result[0] = src[1];
2667  result[1] = -src[0];
2668 
2669  return result;
2670 }
2671 
2672 
2682 template <int dim, typename Number1, typename Number2>
2683 constexpr inline DEAL_II_ALWAYS_INLINE
2685  cross_product_3d(const Tensor<1, dim, Number1> &src1,
2686  const Tensor<1, dim, Number2> &src2)
2687 {
2688  Assert(dim == 3, ExcInternalError());
2689 
2691 
2692  // avoid compiler warnings
2693  constexpr int s0 = 0 % dim;
2694  constexpr int s1 = 1 % dim;
2695  constexpr int s2 = 2 % dim;
2696 
2697  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2698  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2699  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2700 
2701  return result;
2702 }
2703 
2704 
2718 template <int dim, typename Number>
2719 constexpr inline DEAL_II_ALWAYS_INLINE Number
2721 {
2722  // Compute the determinant using the Laplace expansion of the
2723  // determinant. We expand along the last row.
2724  Number det = internal::NumberType<Number>::value(0.0);
2725 
2726  for (unsigned int k = 0; k < dim; ++k)
2727  {
2728  Tensor<2, dim - 1, Number> minor;
2729  for (unsigned int i = 0; i < dim - 1; ++i)
2730  for (unsigned int j = 0; j < dim - 1; ++j)
2731  minor[i][j] = t[i][j < k ? j : j + 1];
2732 
2733  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2734 
2735  det += t[dim - 1][k] * cofactor;
2736  }
2737 
2738  return ((dim % 2 == 0) ? 1. : -1.) * det;
2739 }
2740 
2746 template <typename Number>
2747 constexpr DEAL_II_ALWAYS_INLINE Number
2749 {
2750  return t[0][0];
2751 }
2752 
2758 template <typename Number>
2759 constexpr DEAL_II_ALWAYS_INLINE Number
2761 {
2762  // hard-coded for efficiency reasons
2763  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2764 }
2765 
2771 template <typename Number>
2772 constexpr DEAL_II_ALWAYS_INLINE Number
2774 {
2775  // hard-coded for efficiency reasons
2776  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2777  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2778  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2779  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2780  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2781  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2782  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2783 }
2784 
2785 
2792 template <int dim, typename Number>
2793 constexpr inline DEAL_II_ALWAYS_INLINE Number
2795 {
2796  Number t = d[0][0];
2797  for (unsigned int i = 1; i < dim; ++i)
2798  t += d[i][i];
2799  return t;
2800 }
2801 
2802 
2811 template <int dim, typename Number>
2812 constexpr inline Tensor<2, dim, Number>
2814 {
2815  Number return_tensor[dim][dim];
2816 
2817  // if desired, take over the
2818  // inversion of a 4x4 tensor
2819  // from the FullMatrix
2820  AssertThrow(false, ExcNotImplemented());
2821 
2822  return Tensor<2, dim, Number>(return_tensor);
2823 }
2824 
2825 
2826 #ifndef DOXYGEN
2827 
2828 template <typename Number>
2830  invert(const Tensor<2, 1, Number> &t)
2831 {
2832  Tensor<2, 1, Number> return_tensor;
2833 
2834  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2835 
2836  return return_tensor;
2837 }
2838 
2839 
2840 template <typename Number>
2842  invert(const Tensor<2, 2, Number> &t)
2843 {
2844  Tensor<2, 2, Number> return_tensor;
2845 
2846  const Number inv_det_t = internal::NumberType<Number>::value(
2847  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2848  return_tensor[0][0] = t[1][1];
2849  return_tensor[0][1] = -t[0][1];
2850  return_tensor[1][0] = -t[1][0];
2851  return_tensor[1][1] = t[0][0];
2852  return_tensor *= inv_det_t;
2853 
2854  return return_tensor;
2855 }
2856 
2857 
2858 template <typename Number>
2860  invert(const Tensor<2, 3, Number> &t)
2861 {
2862  Tensor<2, 3, Number> return_tensor;
2863 
2864  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2865  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2866  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2867  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2868  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2869  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2870  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2871  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2872  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2873  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2874  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2875  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2876  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2877  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2878  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2879  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2880  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2881  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2882  const Number inv_det_t = internal::NumberType<Number>::value(
2883  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2884  t[0][2] * return_tensor[2][0]));
2885  return_tensor *= inv_det_t;
2886 
2887  return return_tensor;
2888 }
2889 
2890 #endif /* DOXYGEN */
2891 
2892 
2898 template <int dim, typename Number>
2901 {
2903  for (unsigned int i = 0; i < dim; ++i)
2904  {
2905  tt[i][i] = t[i][i];
2906  for (unsigned int j = i + 1; j < dim; ++j)
2907  {
2908  tt[i][j] = t[j][i];
2909  tt[j][i] = t[i][j];
2910  };
2911  }
2912  return tt;
2913 }
2914 
2915 
2929 template <int dim, typename Number>
2930 constexpr Tensor<2, dim, Number>
2931 adjugate(const Tensor<2, dim, Number> &t)
2932 {
2933  return determinant(t) * invert(t);
2934 }
2935 
2936 
2950 template <int dim, typename Number>
2951 constexpr Tensor<2, dim, Number>
2952 cofactor(const Tensor<2, dim, Number> &t)
2953 {
2954  return transpose(adjugate(t));
2955 }
2956 
2957 
3021 template <int dim, typename Number>
3024 
3025 
3033 template <int dim, typename Number>
3034 inline Number
3036 {
3038  for (unsigned int j = 0; j < dim; ++j)
3039  {
3041  for (unsigned int i = 0; i < dim; ++i)
3043 
3044  if (sum > max)
3045  max = sum;
3046  }
3047 
3048  return max;
3049 }
3050 
3051 
3059 template <int dim, typename Number>
3060 inline Number
3062 {
3064  for (unsigned int i = 0; i < dim; ++i)
3065  {
3067  for (unsigned int j = 0; j < dim; ++j)
3069 
3070  if (sum > max)
3071  max = sum;
3072  }
3073 
3074  return max;
3075 }
3076 
3082 #ifndef DOXYGEN
3083 
3084 
3085 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3086 
3087 // Specialization of functions for ADOL-C number types when
3088 // the advanced branching feature is used
3089 template <int dim>
3090 inline adouble
3092 {
3094  for (unsigned int j = 0; j < dim; ++j)
3095  {
3097  for (unsigned int i = 0; i < dim; ++i)
3098  sum += fabs(t[i][j]);
3099 
3100  condassign(max, (sum > max), sum, max);
3101  }
3102 
3103  return max;
3104 }
3105 
3106 
3107 template <int dim>
3108 inline adouble
3110 {
3112  for (unsigned int i = 0; i < dim; ++i)
3113  {
3115  for (unsigned int j = 0; j < dim; ++j)
3116  sum += fabs(t[i][j]);
3117 
3118  condassign(max, (sum > max), sum, max);
3119  }
3120 
3121  return max;
3122 }
3123 
3124 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3125 
3126 
3127 #endif // DOXYGEN
3128 
3130 
3131 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
std::size_t size() const
Definition: array_view.h:573
Definition: point.h:112
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
const Number * begin_raw() const
constexpr Tensor & operator+=(const Tensor< 0, dim, OtherNumber > &rhs)
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor & operator*=(const OtherNumber &factor)
constexpr Tensor(const Tensor< 0, dim, OtherNumber > &initializer)
constexpr Tensor(const OtherNumber &initializer)
constexpr void clear()
Iterator unroll_recursion(const Iterator current, const Iterator end) const
constexpr Tensor & operator=(const OtherNumber &d) &&=delete
constexpr real_type norm_square() const
constexpr Tensor & operator-=(const Tensor< 0, dim, OtherNumber > &rhs)
const Number * end_raw() const
constexpr bool operator!=(const Tensor< 0, dim, OtherNumber > &rhs) const
real_type norm() const
constexpr Tensor & operator=(const Tensor< 0, dim, OtherNumber > &rhs)
void unroll(const Iterator begin, const Iterator end) const
constexpr bool operator==(const Tensor< 0, dim, OtherNumber > &rhs) const
typename numbers::NumberTraits< Number >::real_type real_type
Definition: tensor.h:126
constexpr Tensor & operator/=(const OtherNumber &factor)
constexpr Tensor operator-() const
constexpr Tensor & operator=(const OtherNumber &d) &
Definition: tensor.h:516
constexpr Tensor(const ArrayView< ElementType, MemorySpace > &initializer)
constexpr bool operator==(const Tensor< rank_, dim, OtherNumber > &) const
constexpr Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:2193
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:556
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
constexpr Tensor & operator*=(const OtherNumber &factor)
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
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:2036
static constexpr unsigned int rank
Definition: tensor.h:535
constexpr const Number & operator[](const TableIndices< rank_ > &indices) const
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
constexpr const value_type & operator[](const unsigned int i) const
constexpr void clear()
void unroll(const Iterator begin, const Iterator end) const
void unroll(Vector< OtherNumber > &result) const
constexpr Number & operator[](const TableIndices< rank_ > &indices)
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Iterator unroll_recursion(const Iterator current, const Iterator end) const
constexpr Tensor(const Tensor< 1, dim, Tensor< rank_ - 1, dim, OtherNumber >> &initializer)
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
constexpr ProductType< Number, OtherNumber >::type operator*(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition: tensor.h:2015
constexpr Tensor & operator=(const Number &d) &
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Number * end_raw()
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:2255
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
const Number * begin_raw() const
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:2069
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:549
friend class Tensor
Definition: tensor.h:907
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:2211
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3061
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1975
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3035
constexpr Tensor & operator/=(const OtherNumber &factor)
static constexpr unsigned int dimension
Definition: tensor.h:530
constexpr Tensor< rank, dim, typename ProductType< typename EnableIfScalar< Number >::type, OtherNumber >::type > operator*(const Number &factor, const Tensor< rank, dim, OtherNumber > &t)
Definition: tensor.h:2122
static constexpr std::size_t memory_consumption()
Number * begin_raw()
constexpr Tensor & operator=(const Number &d) &&=delete
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:2052
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< rank, dim, Number > &src1, const Tensor< rank, dim, OtherNumber > &src2)
Definition: tensor.h:2284
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &p, const Tensor< rank, dim, OtherNumber > &q)
Definition: tensor.h:2235
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:878
constexpr value_type & operator[](const unsigned int i)
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition: tensor.h:1995
constexpr Tensor(const ArrayLike &initializer, std::index_sequence< Indices... >)
const Number * end_raw() const
constexpr Tensor< rank, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:2094
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor(const array_type &initializer)
static constexpr unsigned int n_independent_components
Definition: tensor.h:541
constexpr Tensor operator-() const
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:884
numbers::NumberTraits< Number >::real_type norm() const
Definition: vector.h:110
iterator end()
iterator begin()
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:110
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Expression fabs(const Expression &x)
static const char A
static const char T
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr void contract(T1 &result, const T2 &left, const T3 &right)
constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void abort(const ExceptionBase &exc) noexcept
Definition: exceptions.cc:460
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:2143
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:937
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:921
#define DEAL_II_HOST_DEVICE
Definition: numbers.h:35
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition: numbers.h:37
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:703
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
Definition: numbers.h:594
static constexpr real_type abs_square(const number &x)
Definition: numbers.h:585
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Definition: tensor.cc:82
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3061
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3035