Reference documentation for deal.II version GIT 85919f3e70 2023-05-28 07:10:01+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 - 2022 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 
586  template <typename ElementType, typename MemorySpace>
587  constexpr DEAL_II_HOST_DEVICE explicit Tensor(
588  const ArrayView<ElementType, MemorySpace> &initializer);
589 
597  template <typename OtherNumber>
598  constexpr DEAL_II_HOST_DEVICE
600 
604  template <typename OtherNumber>
605  constexpr Tensor(
606  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
607 
611  template <typename OtherNumber>
612  constexpr
613  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
614 
615 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
619  constexpr Tensor(const Tensor<rank_, dim, Number> &);
620 
624  constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
625 #endif
626 
632  constexpr DEAL_II_HOST_DEVICE value_type &
633  operator[](const unsigned int i);
634 
640  constexpr DEAL_II_HOST_DEVICE const value_type &
641  operator[](const unsigned int i) const;
642 
646  constexpr const Number &
647  operator[](const TableIndices<rank_> &indices) const;
648 
652  constexpr Number &
654 
658  Number *
660 
664  const Number *
665  begin_raw() const;
666 
670  Number *
672 
676  const Number *
677  end_raw() const;
678 
686  template <typename OtherNumber>
687  constexpr DEAL_II_HOST_DEVICE Tensor &
689 
696  constexpr Tensor &
697  operator=(const Number &d) &;
698 
704  constexpr DEAL_II_HOST_DEVICE Tensor &
705  operator=(const Number &d) && = delete;
706 
707 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
711  constexpr Tensor<rank_, dim, Number> &
713 
717  constexpr Tensor<rank_, dim, Number> &
719 #endif
720 
724  template <typename OtherNumber>
725  constexpr bool
727 
731  template <typename OtherNumber>
732  constexpr bool
734 
740  template <typename OtherNumber>
741  constexpr DEAL_II_HOST_DEVICE Tensor &
743 
749  template <typename OtherNumber>
750  constexpr DEAL_II_HOST_DEVICE Tensor &
752 
759  template <typename OtherNumber>
760  constexpr DEAL_II_HOST_DEVICE Tensor &
761  operator*=(const OtherNumber &factor);
762 
768  template <typename OtherNumber>
769  constexpr DEAL_II_HOST_DEVICE Tensor &
770  operator/=(const OtherNumber &factor);
771 
777  constexpr DEAL_II_HOST_DEVICE Tensor
778  operator-() const;
779 
792  constexpr void
793  clear();
794 
804  norm() const;
805 
812  constexpr DEAL_II_HOST_DEVICE
814  norm_square() const;
815 
826  template <typename OtherNumber>
827  DEAL_II_DEPRECATED void
828  unroll(Vector<OtherNumber> &result) const;
829 
840  template <class Iterator>
841  void
842  unroll(const Iterator begin, const Iterator end) const;
843 
848  static constexpr unsigned int
850 
856  static constexpr TableIndices<rank_>
857  unrolled_to_component_indices(const unsigned int i);
858 
863  static constexpr std::size_t
865 
871  template <class Archive>
872  void
873  serialize(Archive &ar, const unsigned int version);
874 
880 
881 private:
885  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
886  // ... avoid a compiler warning in case of dim == 0 and ensure that the
887  // array always has positive size.
888 
892  template <typename Iterator>
893  Iterator
894  unroll_recursion(const Iterator current, const Iterator end) const;
895 
902  template <typename ArrayLike, std::size_t... Indices>
903  constexpr DEAL_II_HOST_DEVICE
904  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
905 
906  // Allow an arbitrary Tensor to access the underlying values.
907  template <int, int, typename>
908  friend class Tensor;
909 
910  // Point is allowed access to the coordinates. This is supposed to improve
911  // speed.
912  friend class Point<dim, Number>;
913 };
914 
915 
916 #ifndef DOXYGEN
917 namespace internal
918 {
919  // Workaround: The following 4 overloads are necessary to be able to
920  // compile the library with Apple Clang 8 and older. We should remove
921  // these overloads again when we bump the minimal required version to
922  // something later than clang-3.6 / Apple Clang 6.3.
923  template <int rank, int dim, typename T, typename U>
924  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
925  {
926  using type =
928  };
929 
930  template <int rank, int dim, typename T, typename U>
931  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
932  {
933  using type =
935  };
936 
937  template <typename T, int rank, int dim, typename U>
938  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
939  {
940  using type =
942  };
943 
944  template <int rank, int dim, typename T, typename U>
945  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
946  {
947  using type =
949  };
950  // end workaround
951 
956  template <int rank, int dim, typename T>
957  struct NumberType<Tensor<rank, dim, T>>
958  {
959  static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const
961  value(const Tensor<rank, dim, T> &t)
962  {
963  return t;
964  }
965 
967  value(const T &t)
968  {
970  tmp = t;
971  return tmp;
972  }
973  };
974 } // namespace internal
975 
976 
977 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
978 
979 
980 template <int dim, typename Number>
983  // Some auto-differentiable numbers need explicit
984  // zero initialization such as adtl::adouble.
985  : Tensor{0.0}
986 {}
987 
988 
989 
990 template <int dim, typename Number>
991 template <typename OtherNumber>
993 Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
994  : value(internal::NumberType<Number>::value(initializer))
995 {}
996 
997 
998 
999 template <int dim, typename Number>
1000 template <typename OtherNumber>
1003  : Tensor{p.value}
1004 {}
1005 
1006 
1007 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1008 template <int dim, typename Number>
1011  : value{other.value}
1012 {}
1013 
1014 
1015 
1016 template <int dim, typename Number>
1019  : value{std::move(other.value)}
1020 {}
1021 # endif
1022 
1023 
1024 template <int dim, typename Number>
1025 inline Number *
1027 {
1028  return std::addressof(value);
1029 }
1030 
1031 
1032 
1033 template <int dim, typename Number>
1034 inline const Number *
1036 {
1037  return std::addressof(value);
1038 }
1039 
1040 
1041 
1042 template <int dim, typename Number>
1043 inline Number *
1045 {
1047 }
1048 
1049 
1050 
1051 template <int dim, typename Number>
1052 const Number *
1054 {
1056 }
1057 
1058 
1059 
1060 template <int dim, typename Number>
1063 {
1064  Assert(dim != 0,
1065  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1066  return value;
1067 }
1068 
1069 
1070 template <int dim, typename Number>
1071 constexpr inline DEAL_II_ALWAYS_INLINE
1073 {
1074  Assert(dim != 0,
1075  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1076  return value;
1077 }
1078 
1079 
1080 
1081 template <int dim, typename Number>
1082 template <typename OtherNumber>
1085 {
1087  return *this;
1088 }
1089 
1090 
1091 # if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1092 template <int dim, typename Number>
1095 {
1096  value = p.value;
1097  return *this;
1098 }
1099 # endif
1100 
1101 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1102 template <int dim, typename Number>
1105 {
1106  value = std::move(other.value);
1107  return *this;
1108 }
1109 # endif
1110 
1111 
1112 
1113 template <int dim, typename Number>
1114 template <typename OtherNumber>
1116 Tensor<0, dim, Number>::operator=(const OtherNumber &d) &
1117 {
1119  return *this;
1120 }
1121 
1122 
1123 template <int dim, typename Number>
1124 template <typename OtherNumber>
1125 constexpr inline bool
1127 {
1128 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1129  Assert(!(std::is_same<Number, adouble>::value ||
1130  std::is_same<OtherNumber, adouble>::value),
1131  ExcMessage(
1132  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1133  "been extended to support advanced branching."));
1134 # endif
1135 
1136  return numbers::values_are_equal(value, p.value);
1137 }
1138 
1139 
1140 template <int dim, typename Number>
1141 template <typename OtherNumber>
1142 constexpr bool
1144 {
1145  return !((*this) == p);
1146 }
1147 
1148 
1149 template <int dim, typename Number>
1150 template <typename OtherNumber>
1153 {
1154  value += p.value;
1155  return *this;
1156 }
1157 
1158 
1159 template <int dim, typename Number>
1160 template <typename OtherNumber>
1163 {
1164  value -= p.value;
1165  return *this;
1166 }
1167 
1168 
1169 
1170 namespace internal
1171 {
1172  namespace ComplexWorkaround
1173  {
1174  template <typename Number, typename OtherNumber>
1175  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1176  multiply_assign_scalar(Number &val, const OtherNumber &s)
1177  {
1178  val *= s;
1179  }
1180 
1181  template <typename Number, typename OtherNumber>
1182  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1183  multiply_assign_scalar(std::complex<Number> &val, const OtherNumber &s)
1184  {
1185 # if KOKKOS_VERSION >= 30600
1186  KOKKOS_IF_ON_HOST((val *= s;))
1187  KOKKOS_IF_ON_DEVICE(({
1188  (void)val;
1189  (void)s;
1190  Kokkos::abort(
1191  "This function is not implemented for std::complex<Number>!\n");
1192  }))
1193 # else
1194 # ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1195  val *= s;
1196 # else
1197  (void)val;
1198  (void)s;
1199  Kokkos::abort(
1200  "This function is not implemented for std::complex<Number>!\n");
1201 # endif
1202 # endif
1203  }
1204  } // namespace ComplexWorkaround
1205 } // namespace internal
1206 
1207 
1208 template <int dim, typename Number>
1209 template <typename OtherNumber>
1211 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1212 {
1213  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1214  return *this;
1215 }
1216 
1217 
1218 
1219 template <int dim, typename Number>
1220 template <typename OtherNumber>
1222 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1223 {
1224  value /= s;
1225  return *this;
1226 }
1227 
1228 
1229 template <int dim, typename Number>
1232 {
1233  return -value;
1234 }
1235 
1236 
1237 template <int dim, typename Number>
1240 {
1241  Assert(dim != 0,
1242  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1243  return numbers::NumberTraits<Number>::abs(value);
1244 }
1245 
1246 
1247 template <int dim, typename Number>
1251 {
1252  Assert(dim != 0,
1253  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1255 }
1256 
1257 
1258 
1259 template <int dim, typename Number>
1260 template <typename Iterator>
1261 Iterator
1262 Tensor<0, dim, Number>::unroll_recursion(const Iterator current,
1263  const Iterator end) const
1264 {
1265  (void)end;
1266  Assert(dim != 0,
1267  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1268  Assert(std::distance(current, end) >= 1,
1269  ExcMessage("The provided iterator range must contain at least one "
1270  "element."));
1271  *current = value;
1272  return std::next(current);
1273 }
1274 
1275 
1276 
1277 template <int dim, typename Number>
1278 constexpr inline void
1280 {
1281  // Some auto-differentiable numbers need explicit
1282  // zero initialization.
1284 }
1285 
1286 
1287 
1288 template <int dim, typename Number>
1289 template <class Iterator>
1290 inline void
1291 Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1292 {
1295 }
1296 
1297 
1298 
1299 template <int dim, typename Number>
1300 template <class Archive>
1301 inline void
1302 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1303 {
1304  ar &value;
1305 }
1306 
1307 
1308 template <int dim, typename Number>
1310 
1311 
1312 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1313 
1314 template <int rank_, int dim, typename Number>
1315 template <typename ArrayLike, std::size_t... indices>
1317 Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1318  std::index_sequence<indices...>)
1319  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1320 {
1321  static_assert(sizeof...(indices) == dim,
1322  "dim should match the number of indices");
1323 }
1324 
1325 
1326 
1327 template <int rank_, int dim, typename Number>
1330  // We would like to use =default, but this causes compile errors with some
1331  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1332  : values{}
1333 {}
1334 
1335 
1336 
1337 template <int rank_, int dim, typename Number>
1339 Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1340  : Tensor(initializer, std::make_index_sequence<dim>{})
1341 {}
1342 
1343 
1344 
1345 template <int rank_, int dim, typename Number>
1346 template <typename ElementType, typename MemorySpace>
1349  const ArrayView<ElementType, MemorySpace> &initializer)
1350 {
1351  // make nvcc happy
1352  const int my_n_independent_components = n_independent_components;
1353  AssertDimension(initializer.size(), my_n_independent_components);
1354 
1355  for (unsigned int i = 0; i < my_n_independent_components; ++i)
1356  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1357 }
1358 
1359 
1360 
1361 template <int rank_, int dim, typename Number>
1362 template <typename OtherNumber>
1365  const Tensor<rank_, dim, OtherNumber> &initializer)
1366  : Tensor(initializer, std::make_index_sequence<dim>{})
1367 {}
1368 
1369 
1370 
1371 template <int rank_, int dim, typename Number>
1372 template <typename OtherNumber>
1373 constexpr DEAL_II_ALWAYS_INLINE
1375  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1376  : Tensor(initializer, std::make_index_sequence<dim>{})
1377 {}
1378 
1379 
1380 
1381 template <int rank_, int dim, typename Number>
1382 template <typename OtherNumber>
1384 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1385 {
1386  return Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>(values);
1387 }
1388 
1389 
1390 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1391 template <int rank_, int dim, typename Number>
1392 constexpr DEAL_II_ALWAYS_INLINE
1394 {
1395  for (unsigned int i = 0; i < dim; ++i)
1396  values[i] = other.values[i];
1397 }
1398 
1399 
1400 
1401 template <int rank_, int dim, typename Number>
1402 constexpr DEAL_II_ALWAYS_INLINE
1404 {
1405  for (unsigned int i = 0; i < dim; ++i)
1406  values[i] = other.values[i];
1407 }
1408 # endif
1409 
1410 namespace internal
1411 {
1412  namespace TensorSubscriptor
1413  {
1414  template <typename ArrayElementType, int dim>
1415  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
1416  subscript(ArrayElementType * values,
1417  const unsigned int i,
1418  std::integral_constant<int, dim>)
1419  {
1420  AssertIndexRange(i, dim);
1421  return values[i];
1422  }
1423 
1424  template <typename ArrayElementType>
1425  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
1426  subscript(ArrayElementType *dummy,
1427  const unsigned int,
1428  std::integral_constant<int, 0>)
1429  {
1430  Assert(
1431  false,
1432  ExcMessage(
1433  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1434  return *dummy;
1435  }
1436  } // namespace TensorSubscriptor
1437 } // namespace internal
1438 
1439 
1440 template <int rank_, int dim, typename Number>
1443  Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1444 {
1445  return ::internal::TensorSubscriptor::subscript(
1446  values, i, std::integral_constant<int, dim>());
1447 }
1448 
1449 
1450 template <int rank_, int dim, typename Number>
1451 constexpr DEAL_II_ALWAYS_INLINE
1453  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1454 {
1455 # if KOKKOS_VERSION < 30700
1456 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1457  AssertIndexRange(i, dim);
1458 # endif
1459 # else
1460  KOKKOS_IF_ON_HOST((AssertIndexRange(i, dim);))
1461 # endif
1462 
1463  return values[i];
1464 }
1465 
1466 
1467 template <int rank_, int dim, typename Number>
1468 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1470 {
1471 # if KOKKOS_VERSION < 30700
1472 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1473  Assert(dim != 0,
1474  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1475 # endif
1476 # else
1477  KOKKOS_IF_ON_HOST(
1478  (Assert(dim != 0,
1479  ExcMessage(
1480  "Cannot access an object of type Tensor<rank_,0,Number>"));))
1481 # endif
1482 
1483  return TensorAccessors::extract<rank_>(*this, indices);
1484 }
1485 
1486 
1487 
1488 template <int rank_, int dim, typename Number>
1489 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1491 {
1492 # if KOKKOS_VERSION < 30700
1493 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1494  Assert(dim != 0,
1495  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1496 # endif
1497 # else
1498  KOKKOS_IF_ON_HOST(
1499  (Assert(dim != 0,
1500  ExcMessage(
1501  "Cannot access an object of type Tensor<rank_,0,Number>"));))
1502 # endif
1503 
1504  return TensorAccessors::extract<rank_>(*this, indices);
1505 }
1506 
1507 
1508 
1509 template <int rank_, int dim, typename Number>
1510 inline Number *
1512 {
1513  return std::addressof(
1514  this->operator[](this->unrolled_to_component_indices(0)));
1515 }
1516 
1517 
1518 
1519 template <int rank_, int dim, typename Number>
1520 inline const Number *
1522 {
1523  return std::addressof(
1524  this->operator[](this->unrolled_to_component_indices(0)));
1525 }
1526 
1527 
1528 
1529 template <int rank_, int dim, typename Number>
1530 inline Number *
1532 {
1533  return begin_raw() + n_independent_components;
1534 }
1535 
1536 
1537 
1538 template <int rank_, int dim, typename Number>
1539 inline const Number *
1541 {
1542  return begin_raw() + n_independent_components;
1543 }
1544 
1545 
1546 
1547 template <int rank_, int dim, typename Number>
1548 template <typename OtherNumber>
1551 {
1552  // The following loop could be written more concisely using std::copy, but
1553  // that function is only constexpr from C++20 on.
1554  for (unsigned int i = 0; i < dim; ++i)
1555  values[i] = t.values[i];
1556  return *this;
1557 }
1558 
1559 
1560 
1561 template <int rank_, int dim, typename Number>
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<Number, Differentiation::SD::Expression>::value,
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<Number, Differentiation::SD::Expression>::value,
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  auto 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  inline constexpr unsigned int
1824  mod(const unsigned int x)
1825  {
1826  return x % dim;
1827  }
1828 
1829  template <>
1830  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  inline constexpr unsigned int
1839  div(const unsigned int x)
1840  {
1841  return x / dim;
1842  }
1843 
1844  template <>
1845  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  AssertIndexRange(i, n_independent_components);
1861 
1862  TableIndices<rank_> indices;
1863 
1864  unsigned int remainder = i;
1865  for (int r = rank_ - 1; r >= 0; --r)
1866  {
1867  indices[r] = internal::mod<dim>(remainder);
1868  remainder = internal::div<dim>(remainder);
1869  }
1870  Assert(remainder == 0, ExcInternalError());
1871 
1872  return indices;
1873 }
1874 
1875 
1876 template <int rank_, int dim, typename Number>
1877 constexpr inline void
1879 {
1880  for (unsigned int i = 0; i < dim; ++i)
1882 }
1883 
1884 
1885 template <int rank_, int dim, typename Number>
1886 constexpr std::size_t
1888 {
1889  return sizeof(Tensor<rank_, dim, Number>);
1890 }
1891 
1892 
1893 template <int rank_, int dim, typename Number>
1894 template <class Archive>
1895 inline void
1896 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1897 {
1898  ar &values;
1899 }
1900 
1901 
1902 template <int rank_, int dim, typename Number>
1904 
1905 #endif // DOXYGEN
1906 
1907 /* ----------------- Non-member functions operating on tensors. ------------ */
1908 
1921 template <int rank_, int dim, typename Number>
1922 inline std::ostream &
1923 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1924 {
1925  for (unsigned int i = 0; i < dim; ++i)
1926  {
1927  out << p[i];
1928  if (i != dim - 1)
1929  out << ' ';
1930  }
1931 
1932  return out;
1933 }
1934 
1935 
1942 template <int dim, typename Number>
1943 inline std::ostream &
1944 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1945 {
1946  out << static_cast<const Number &>(p);
1947  return out;
1948 }
1949 
1950 
1969 template <int dim, typename Number, typename Other>
1972  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1973 {
1974  return object * static_cast<const Number &>(t);
1975 }
1976 
1977 
1978 
1989 template <int dim, typename Number, typename Other>
1992  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1993 {
1994  return static_cast<const Number &>(t) * object;
1995 }
1996 
1997 
2009 template <int dim, typename Number, typename OtherNumber>
2013  const Tensor<0, dim, OtherNumber> &src2)
2014 {
2015  return static_cast<const Number &>(src1) *
2016  static_cast<const OtherNumber &>(src2);
2017 }
2018 
2019 
2027 template <int dim, typename Number, typename OtherNumber>
2029  Tensor<0,
2030  dim,
2031  typename ProductType<Number,
2032  typename EnableIfScalar<OtherNumber>::type>::type>
2033  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2034 {
2035  return static_cast<const Number &>(t) / factor;
2036 }
2037 
2038 
2046 template <int dim, typename Number, typename OtherNumber>
2050  const Tensor<0, dim, OtherNumber> &q)
2051 {
2052  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2053 }
2054 
2055 
2063 template <int dim, typename Number, typename OtherNumber>
2067  const Tensor<0, dim, OtherNumber> &q)
2068 {
2069  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2070 }
2071 
2072 
2085 template <int rank, int dim, typename Number, typename OtherNumber>
2087  Tensor<rank,
2088  dim,
2089  typename ProductType<Number,
2090  typename EnableIfScalar<OtherNumber>::type>::type>
2091  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2092 {
2093  // recurse over the base objects
2095  for (unsigned int d = 0; d < dim; ++d)
2096  tt[d] = t[d] * factor;
2097  return tt;
2098 }
2099 
2100 
2113 template <int rank, int dim, typename Number, typename OtherNumber>
2115  Tensor<rank,
2116  dim,
2118  OtherNumber>::type>
2119  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2120 {
2121  // simply forward to the operator above
2122  return t * factor;
2123 }
2124 
2125 
2126 namespace internal
2127 {
2128  namespace TensorImplementation
2129  {
2130  template <int rank,
2131  int dim,
2132  typename Number,
2133  typename OtherNumber,
2134  std::enable_if_t<
2135  !std::is_integral<
2136  typename ProductType<Number, OtherNumber>::type>::value,
2137  int> = 0>
2141  const OtherNumber & factor)
2142  {
2144  const Number inverse_factor = Number(1.) / factor;
2145  // recurse over the base objects
2146  for (unsigned int d = 0; d < dim; ++d)
2147  tt[d] = t[d] * inverse_factor;
2148  return tt;
2149  }
2150 
2151 
2152  template <int rank,
2153  int dim,
2154  typename Number,
2155  typename OtherNumber,
2156  std::enable_if_t<
2157  std::is_integral<
2158  typename ProductType<Number, OtherNumber>::type>::value,
2159  int> = 0>
2163  const OtherNumber & factor)
2164  {
2166  // recurse over the base objects
2167  for (unsigned int d = 0; d < dim; ++d)
2168  tt[d] = t[d] / factor;
2169  return tt;
2170  }
2171  } // namespace TensorImplementation
2172 } // namespace internal
2173 
2174 
2184 template <int rank, int dim, typename Number, typename OtherNumber>
2186  Tensor<rank,
2187  dim,
2188  typename ProductType<Number,
2189  typename EnableIfScalar<OtherNumber>::type>::type>
2190  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2191 {
2193 }
2194 
2195 
2205 template <int rank, int dim, typename Number, typename OtherNumber>
2210 {
2212 
2213  for (unsigned int i = 0; i < dim; ++i)
2214  tmp[i] += q[i];
2215 
2216  return tmp;
2217 }
2218 
2219 
2229 template <int rank, int dim, typename Number, typename OtherNumber>
2234 {
2236 
2237  for (unsigned int i = 0; i < dim; ++i)
2238  tmp[i] -= q[i];
2239 
2240  return tmp;
2241 }
2242 
2249 template <int dim, typename Number, typename OtherNumber>
2250 inline constexpr DEAL_II_ALWAYS_INLINE
2253  const Tensor<0, dim, OtherNumber> &src2)
2254 {
2256 
2257  tmp *= src2;
2258 
2259  return tmp;
2260 }
2261 
2278 template <int rank, int dim, typename Number, typename OtherNumber>
2279 inline constexpr DEAL_II_ALWAYS_INLINE
2282  const Tensor<rank, dim, OtherNumber> &src2)
2283 {
2285 
2286  for (unsigned int i = 0; i < dim; ++i)
2287  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2289 
2290  return tmp;
2291 }
2292 
2337 template <int rank_1,
2338  int rank_2,
2339  int dim,
2340  typename Number,
2341  typename OtherNumber,
2342  typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2343 constexpr inline DEAL_II_ALWAYS_INLINE
2344  typename Tensor<rank_1 + rank_2 - 2,
2345  dim,
2346  typename ProductType<Number, OtherNumber>::type>::tensor_type
2349 {
2350  typename Tensor<rank_1 + rank_2 - 2,
2351  dim,
2352  typename ProductType<Number, OtherNumber>::type>::tensor_type
2353  result{};
2354 
2355  TensorAccessors::internal::
2356  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2357  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2358  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2359 
2360  return result;
2361 }
2362 
2363 
2392 template <int index_1,
2393  int index_2,
2394  int rank_1,
2395  int rank_2,
2396  int dim,
2397  typename Number,
2398  typename OtherNumber>
2399 constexpr inline DEAL_II_ALWAYS_INLINE
2400  typename Tensor<rank_1 + rank_2 - 2,
2401  dim,
2402  typename ProductType<Number, OtherNumber>::type>::tensor_type
2405 {
2406  Assert(0 <= index_1 && index_1 < rank_1,
2407  ExcMessage(
2408  "The specified index_1 must lie within the range [0,rank_1)"));
2409  Assert(0 <= index_2 && index_2 < rank_2,
2410  ExcMessage(
2411  "The specified index_2 must lie within the range [0,rank_2)"));
2412 
2413  using namespace TensorAccessors;
2414  using namespace TensorAccessors::internal;
2415 
2416  // Reorder index_1 to the end of src1:
2418  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2419 
2420  // Reorder index_2 to the end of src2:
2421  const ReorderedIndexView<index_2,
2422  rank_2,
2424  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2425 
2426  typename Tensor<rank_1 + rank_2 - 2,
2427  dim,
2428  typename ProductType<Number, OtherNumber>::type>::tensor_type
2429  result{};
2430  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2431  return result;
2432 }
2433 
2434 
2465 template <int index_1,
2466  int index_2,
2467  int index_3,
2468  int index_4,
2469  int rank_1,
2470  int rank_2,
2471  int dim,
2472  typename Number,
2473  typename OtherNumber>
2474 constexpr inline
2475  typename Tensor<rank_1 + rank_2 - 4,
2476  dim,
2477  typename ProductType<Number, OtherNumber>::type>::tensor_type
2478  double_contract(const Tensor<rank_1, dim, Number> & src1,
2480 {
2481  Assert(0 <= index_1 && index_1 < rank_1,
2482  ExcMessage(
2483  "The specified index_1 must lie within the range [0,rank_1)"));
2484  Assert(0 <= index_3 && index_3 < rank_1,
2485  ExcMessage(
2486  "The specified index_3 must lie within the range [0,rank_1)"));
2487  Assert(index_1 != index_3,
2488  ExcMessage("index_1 and index_3 must not be the same"));
2489  Assert(0 <= index_2 && index_2 < rank_2,
2490  ExcMessage(
2491  "The specified index_2 must lie within the range [0,rank_2)"));
2492  Assert(0 <= index_4 && index_4 < rank_2,
2493  ExcMessage(
2494  "The specified index_4 must lie within the range [0,rank_2)"));
2495  Assert(index_2 != index_4,
2496  ExcMessage("index_2 and index_4 must not be the same"));
2497 
2498  using namespace TensorAccessors;
2499  using namespace TensorAccessors::internal;
2500 
2501  // Reorder index_1 to the end of src1:
2503  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2504 
2505  // Reorder index_2 to the end of src2:
2507  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2508 
2509  // Now, reorder index_3 to the end of src1. We have to make sure to
2510  // preserve the original ordering: index_1 has been removed. If
2511  // index_3 > index_1, we have to use (index_3 - 1) instead:
2513  (index_3 < index_1 ? index_3 : index_3 - 1),
2514  rank_1,
2515  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2516  reord_3 =
2517  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2518  index_3 - 1,
2519  rank_1 > (reord_1);
2520 
2521  // Now, reorder index_4 to the end of src2. We have to make sure to
2522  // preserve the original ordering: index_2 has been removed. If
2523  // index_4 > index_2, we have to use (index_4 - 1) instead:
2525  (index_4 < index_2 ? index_4 : index_4 - 1),
2526  rank_2,
2527  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2528  reord_4 =
2529  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2530  index_4 - 1,
2531  rank_2 > (reord_2);
2532 
2533  typename Tensor<rank_1 + rank_2 - 4,
2534  dim,
2535  typename ProductType<Number, OtherNumber>::type>::tensor_type
2536  result{};
2537  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2538  return result;
2539 }
2540 
2541 
2554 template <int rank, int dim, typename Number, typename OtherNumber>
2555 constexpr inline DEAL_II_ALWAYS_INLINE
2557  scalar_product(const Tensor<rank, dim, Number> & left,
2558  const Tensor<rank, dim, OtherNumber> &right)
2559 {
2560  typename ProductType<Number, OtherNumber>::type result{};
2561  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2562  return result;
2563 }
2564 
2565 
2583 template <template <int, int, typename> class TensorT1,
2584  template <int, int, typename>
2585  class TensorT2,
2586  template <int, int, typename>
2587  class TensorT3,
2588  int rank_1,
2589  int rank_2,
2590  int dim,
2591  typename T1,
2592  typename T2,
2593  typename T3>
2594 constexpr inline DEAL_II_ALWAYS_INLINE
2596  contract3(const TensorT1<rank_1, dim, T1> & left,
2597  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2598  const TensorT3<rank_2, dim, T3> & right)
2599 {
2600  using return_type =
2602  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2603  middle,
2604  right);
2605 }
2606 
2607 
2618 template <int rank_1,
2619  int rank_2,
2620  int dim,
2621  typename Number,
2622  typename OtherNumber>
2623 constexpr inline DEAL_II_ALWAYS_INLINE
2627 {
2628  typename Tensor<rank_1 + rank_2,
2629  dim,
2630  typename ProductType<Number, OtherNumber>::type>::tensor_type
2631  result{};
2632  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2633  return result;
2634 }
2635 
2636 
2655 template <int dim, typename Number>
2657 cross_product_2d(const Tensor<1, dim, Number> &src)
2658 {
2659  Assert(dim == 2, ExcInternalError());
2660 
2661  Tensor<1, dim, Number> result;
2662 
2663  result[0] = src[1];
2664  result[1] = -src[0];
2665 
2666  return result;
2667 }
2668 
2669 
2679 template <int dim, typename Number1, typename Number2>
2680 constexpr inline DEAL_II_ALWAYS_INLINE
2682  cross_product_3d(const Tensor<1, dim, Number1> &src1,
2683  const Tensor<1, dim, Number2> &src2)
2684 {
2685  Assert(dim == 3, ExcInternalError());
2686 
2688 
2689  // avoid compiler warnings
2690  constexpr int s0 = 0 % dim;
2691  constexpr int s1 = 1 % dim;
2692  constexpr int s2 = 2 % dim;
2693 
2694  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2695  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2696  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2697 
2698  return result;
2699 }
2700 
2701 
2715 template <int dim, typename Number>
2716 constexpr inline DEAL_II_ALWAYS_INLINE Number
2718 {
2719  // Compute the determinant using the Laplace expansion of the
2720  // determinant. We expand along the last row.
2721  Number det = internal::NumberType<Number>::value(0.0);
2722 
2723  for (unsigned int k = 0; k < dim; ++k)
2724  {
2725  Tensor<2, dim - 1, Number> minor;
2726  for (unsigned int i = 0; i < dim - 1; ++i)
2727  for (unsigned int j = 0; j < dim - 1; ++j)
2728  minor[i][j] = t[i][j < k ? j : j + 1];
2729 
2730  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2731 
2732  det += t[dim - 1][k] * cofactor;
2733  }
2734 
2735  return ((dim % 2 == 0) ? 1. : -1.) * det;
2736 }
2737 
2743 template <typename Number>
2744 constexpr DEAL_II_ALWAYS_INLINE Number
2746 {
2747  return t[0][0];
2748 }
2749 
2755 template <typename Number>
2756 constexpr DEAL_II_ALWAYS_INLINE Number
2758 {
2759  // hard-coded for efficiency reasons
2760  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2761 }
2762 
2768 template <typename Number>
2769 constexpr DEAL_II_ALWAYS_INLINE Number
2771 {
2772  // hard-coded for efficiency reasons
2773  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2774  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2775  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2776  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2777  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2778  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2779  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2780 }
2781 
2782 
2789 template <int dim, typename Number>
2790 constexpr inline DEAL_II_ALWAYS_INLINE Number
2792 {
2793  Number t = d[0][0];
2794  for (unsigned int i = 1; i < dim; ++i)
2795  t += d[i][i];
2796  return t;
2797 }
2798 
2799 
2808 template <int dim, typename Number>
2809 constexpr inline Tensor<2, dim, Number>
2811 {
2812  Number return_tensor[dim][dim];
2813 
2814  // if desired, take over the
2815  // inversion of a 4x4 tensor
2816  // from the FullMatrix
2817  AssertThrow(false, ExcNotImplemented());
2818 
2819  return Tensor<2, dim, Number>(return_tensor);
2820 }
2821 
2822 
2823 #ifndef DOXYGEN
2824 
2825 template <typename Number>
2827  invert(const Tensor<2, 1, Number> &t)
2828 {
2829  Tensor<2, 1, Number> return_tensor;
2830 
2831  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2832 
2833  return return_tensor;
2834 }
2835 
2836 
2837 template <typename Number>
2839  invert(const Tensor<2, 2, Number> &t)
2840 {
2841  Tensor<2, 2, Number> return_tensor;
2842 
2843  const Number inv_det_t = internal::NumberType<Number>::value(
2844  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2845  return_tensor[0][0] = t[1][1];
2846  return_tensor[0][1] = -t[0][1];
2847  return_tensor[1][0] = -t[1][0];
2848  return_tensor[1][1] = t[0][0];
2849  return_tensor *= inv_det_t;
2850 
2851  return return_tensor;
2852 }
2853 
2854 
2855 template <typename Number>
2857  invert(const Tensor<2, 3, Number> &t)
2858 {
2859  Tensor<2, 3, Number> return_tensor;
2860 
2861  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2862  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2863  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2864  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2865  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2866  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2867  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2868  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2869  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2870  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2871  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2872  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2873  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2874  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2875  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2876  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2877  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2878  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2879  const Number inv_det_t = internal::NumberType<Number>::value(
2880  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2881  t[0][2] * return_tensor[2][0]));
2882  return_tensor *= inv_det_t;
2883 
2884  return return_tensor;
2885 }
2886 
2887 #endif /* DOXYGEN */
2888 
2889 
2895 template <int dim, typename Number>
2898 {
2900  for (unsigned int i = 0; i < dim; ++i)
2901  {
2902  tt[i][i] = t[i][i];
2903  for (unsigned int j = i + 1; j < dim; ++j)
2904  {
2905  tt[i][j] = t[j][i];
2906  tt[j][i] = t[i][j];
2907  };
2908  }
2909  return tt;
2910 }
2911 
2912 
2926 template <int dim, typename Number>
2927 constexpr Tensor<2, dim, Number>
2928 adjugate(const Tensor<2, dim, Number> &t)
2929 {
2930  return determinant(t) * invert(t);
2931 }
2932 
2933 
2947 template <int dim, typename Number>
2948 constexpr Tensor<2, dim, Number>
2949 cofactor(const Tensor<2, dim, Number> &t)
2950 {
2951  return transpose(adjugate(t));
2952 }
2953 
2954 
3018 template <int dim, typename Number>
3021 
3022 
3030 template <int dim, typename Number>
3031 inline Number
3033 {
3035  for (unsigned int j = 0; j < dim; ++j)
3036  {
3038  for (unsigned int i = 0; i < dim; ++i)
3040 
3041  if (sum > max)
3042  max = sum;
3043  }
3044 
3045  return max;
3046 }
3047 
3048 
3056 template <int dim, typename Number>
3057 inline Number
3059 {
3061  for (unsigned int i = 0; i < dim; ++i)
3062  {
3064  for (unsigned int j = 0; j < dim; ++j)
3066 
3067  if (sum > max)
3068  max = sum;
3069  }
3070 
3071  return max;
3072 }
3073 
3079 #ifndef DOXYGEN
3080 
3081 
3082 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3083 
3084 // Specialization of functions for ADOL-C number types when
3085 // the advanced branching feature is used
3086 template <int dim>
3087 inline adouble
3089 {
3091  for (unsigned int j = 0; j < dim; ++j)
3092  {
3094  for (unsigned int i = 0; i < dim; ++i)
3095  sum += fabs(t[i][j]);
3096 
3097  condassign(max, (sum > max), sum, max);
3098  }
3099 
3100  return max;
3101 }
3102 
3103 
3104 template <int dim>
3105 inline adouble
3107 {
3109  for (unsigned int i = 0; i < dim; ++i)
3110  {
3112  for (unsigned int j = 0; j < dim; ++j)
3113  sum += fabs(t[i][j]);
3114 
3115  condassign(max, (sum > max), sum, max);
3116  }
3117 
3118  return max;
3119 }
3120 
3121 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3122 
3123 
3124 #endif // DOXYGEN
3125 
3127 
3128 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
std::size_t size() const
Definition: array_view.h:576
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:2190
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:2033
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:2012
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:2252
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:2066
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:549
friend class Tensor
Definition: tensor.h:908
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:2208
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3058
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1972
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3032
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:2119
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:2049
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:2281
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:2232
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:879
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:1992
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:2091
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:885
numbers::NumberTraits< Number >::real_type norm() const
Definition: vector.h:109
iterator end()
iterator begin()
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:108
#define DEAL_II_DEPRECATED
Definition: config.h:174
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:162
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
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:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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:2140
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:939
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:923
#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< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:702
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
Definition: numbers.h:593
static constexpr real_type abs_square(const number &x)
Definition: numbers.h:584
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:3058
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3032