Reference documentation for deal.II version GIT b206511199 2023-01-31 13:40: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 - 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/numbers.h>
26 
27 #ifdef DEAL_II_WITH_ADOLC
28 # include <adolc/adouble.h> // Taped double
29 #endif
30 
31 #include <cmath>
32 #include <ostream>
33 
35 
36 // Forward declarations:
37 #ifndef DOXYGEN
38 template <typename ElementType, typename MemorySpace>
39 class ArrayView;
40 template <int dim, typename Number>
41 class Point;
42 template <int rank_, int dim, typename Number = double>
43 class Tensor;
44 template <typename Number>
45 class Vector;
46 template <typename number>
47 class FullMatrix;
48 namespace Differentiation
49 {
50  namespace SD
51  {
52  class Expression;
53  }
54 } // namespace Differentiation
55 #endif
56 
57 
87 template <int dim, typename Number>
88 class Tensor<0, dim, Number>
89 {
90 public:
91  static_assert(dim >= 0,
92  "Tensors must have a dimension greater than or equal to one.");
93 
102  static constexpr unsigned int dimension = dim;
103 
107  static constexpr unsigned int rank = 0;
108 
112  static constexpr unsigned int n_independent_components = 1;
113 
123 
128  using value_type = Number;
129 
135  using array_type = Number;
136 
142  constexpr DEAL_II_HOST_DEVICE
144 
152  template <typename OtherNumber>
153  constexpr DEAL_II_HOST_DEVICE
154  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
155 
161  template <typename OtherNumber>
162  constexpr DEAL_II_HOST_DEVICE
163  Tensor(const OtherNumber &initializer);
164 
165 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
169  constexpr DEAL_II_HOST_DEVICE
170  Tensor(const Tensor<0, dim, Number> &other);
171 
175  constexpr DEAL_II_HOST_DEVICE
176  Tensor(Tensor<0, dim, Number> &&other) noexcept;
177 #endif
178 
188  Number *
190 
200  const Number *
201  begin_raw() const;
202 
212  Number *
214 
225  const Number *
226  end_raw() const;
227 
237  constexpr DEAL_II_HOST_DEVICE
238  operator Number &();
239 
248  constexpr DEAL_II_HOST_DEVICE operator const Number &() const;
249 
257  template <typename OtherNumber>
258  constexpr DEAL_II_HOST_DEVICE Tensor &
260 
261 #if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
270  constexpr DEAL_II_HOST_DEVICE Tensor &
271  operator=(const Tensor<0, dim, Number> &rhs);
272 #endif
273 
274 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
279  operator=(Tensor<0, dim, Number> &&other) noexcept;
280 #endif
281 
288  template <typename OtherNumber>
289  constexpr DEAL_II_HOST_DEVICE Tensor &
290  operator=(const OtherNumber &d);
291 
295  template <typename OtherNumber>
296  constexpr bool
298 
302  template <typename OtherNumber>
303  constexpr bool
305 
311  template <typename OtherNumber>
312  constexpr DEAL_II_HOST_DEVICE Tensor &
314 
320  template <typename OtherNumber>
321  constexpr DEAL_II_HOST_DEVICE Tensor &
323 
329  template <typename OtherNumber>
330  constexpr DEAL_II_HOST_DEVICE Tensor &
331  operator*=(const OtherNumber &factor);
332 
338  template <typename OtherNumber>
339  constexpr DEAL_II_HOST_DEVICE Tensor &
340  operator/=(const OtherNumber &factor);
341 
347  constexpr DEAL_II_HOST_DEVICE Tensor
348  operator-() const;
349 
362  constexpr void
363  clear();
364 
370  real_type
371  norm() const;
372 
380  norm_square() const;
381 
389  template <class Iterator>
390  void
391  unroll(const Iterator begin, const Iterator end) const;
392 
398  template <class Archive>
399  void
400  serialize(Archive &ar, const unsigned int version);
401 
406  using tensor_type = Number;
407 
408 private:
412  Number value;
413 
417  template <typename Iterator>
418  Iterator
419  unroll_recursion(const Iterator current, const Iterator end) const;
420 
421  // Allow an arbitrary Tensor to access the underlying values.
422  template <int, int, typename>
423  friend class Tensor;
424 };
425 
426 
427 
501 template <int rank_, int dim, typename Number>
502 class Tensor
503 {
504 public:
505  static_assert(rank_ >= 1,
506  "Tensors must have a rank greater than or equal to one.");
507  static_assert(dim >= 0,
508  "Tensors must have a dimension greater than or equal to zero.");
517  static constexpr unsigned int dimension = dim;
518 
522  static constexpr unsigned int rank = rank_;
523 
528  static constexpr unsigned int n_independent_components =
529  Tensor<rank_ - 1, dim>::n_independent_components * dim;
530 
536  using value_type = typename Tensor<rank_ - 1, dim, Number>::tensor_type;
537 
542  using array_type =
543  typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1];
544 
552 
558  constexpr DEAL_II_HOST_DEVICE explicit Tensor(const array_type &initializer);
559 
573  template <typename ElementType, typename MemorySpace>
574  constexpr DEAL_II_HOST_DEVICE explicit Tensor(
575  const ArrayView<ElementType, MemorySpace> &initializer);
576 
584  template <typename OtherNumber>
585  constexpr DEAL_II_HOST_DEVICE
587 
591  template <typename OtherNumber>
592  constexpr Tensor(
593  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
594 
598  template <typename OtherNumber>
599  constexpr
600  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
601 
602 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
606  constexpr Tensor(const Tensor<rank_, dim, Number> &);
607 
611  constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
612 #endif
613 
619  constexpr DEAL_II_HOST_DEVICE value_type &
620  operator[](const unsigned int i);
621 
627  constexpr DEAL_II_HOST_DEVICE const value_type &
628  operator[](const unsigned int i) const;
629 
633  constexpr const Number &
634  operator[](const TableIndices<rank_> &indices) const;
635 
639  constexpr Number &
641 
645  Number *
647 
651  const Number *
652  begin_raw() const;
653 
657  Number *
659 
663  const Number *
664  end_raw() const;
665 
673  template <typename OtherNumber>
674  constexpr DEAL_II_HOST_DEVICE Tensor &
676 
683  constexpr Tensor &
684  operator=(const Number &d);
685 
686 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
690  constexpr Tensor<rank_, dim, Number> &
692 
696  constexpr Tensor<rank_, dim, Number> &
698 #endif
699 
703  template <typename OtherNumber>
704  constexpr bool
706 
710  template <typename OtherNumber>
711  constexpr bool
713 
719  template <typename OtherNumber>
720  constexpr DEAL_II_HOST_DEVICE Tensor &
722 
728  template <typename OtherNumber>
729  constexpr DEAL_II_HOST_DEVICE Tensor &
731 
738  template <typename OtherNumber>
739  constexpr DEAL_II_HOST_DEVICE Tensor &
740  operator*=(const OtherNumber &factor);
741 
747  template <typename OtherNumber>
748  constexpr DEAL_II_HOST_DEVICE Tensor &
749  operator/=(const OtherNumber &factor);
750 
756  constexpr DEAL_II_HOST_DEVICE Tensor
757  operator-() const;
758 
771  constexpr void
772  clear();
773 
783  norm() const;
784 
791  constexpr DEAL_II_HOST_DEVICE
793  norm_square() const;
794 
805  template <typename OtherNumber>
806  DEAL_II_DEPRECATED void
807  unroll(Vector<OtherNumber> &result) const;
808 
819  template <class Iterator>
820  void
821  unroll(const Iterator begin, const Iterator end) const;
822 
827  static constexpr unsigned int
829 
835  static constexpr TableIndices<rank_>
836  unrolled_to_component_indices(const unsigned int i);
837 
842  static constexpr std::size_t
844 
850  template <class Archive>
851  void
852  serialize(Archive &ar, const unsigned int version);
853 
859 
860 private:
864  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
865  // ... avoid a compiler warning in case of dim == 0 and ensure that the
866  // array always has positive size.
867 
871  template <typename Iterator>
872  Iterator
873  unroll_recursion(const Iterator current, const Iterator end) const;
874 
881  template <typename ArrayLike, std::size_t... Indices>
882  constexpr DEAL_II_HOST_DEVICE
883  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
884 
885  // Allow an arbitrary Tensor to access the underlying values.
886  template <int, int, typename>
887  friend class Tensor;
888 
889  // Point is allowed access to the coordinates. This is supposed to improve
890  // speed.
891  friend class Point<dim, Number>;
892 };
893 
894 
895 #ifndef DOXYGEN
896 namespace internal
897 {
898  // Workaround: The following 4 overloads are necessary to be able to
899  // compile the library with Apple Clang 8 and older. We should remove
900  // these overloads again when we bump the minimal required version to
901  // something later than clang-3.6 / Apple Clang 6.3.
902  template <int rank, int dim, typename T, typename U>
903  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
904  {
905  using type =
907  };
908 
909  template <int rank, int dim, typename T, typename U>
910  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
911  {
912  using type =
914  };
915 
916  template <typename T, int rank, int dim, typename U>
917  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
918  {
919  using type =
921  };
922 
923  template <int rank, int dim, typename T, typename U>
924  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
925  {
926  using type =
928  };
929  // end workaround
930 
935  template <int rank, int dim, typename T>
936  struct NumberType<Tensor<rank, dim, T>>
937  {
938  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
939  value(const Tensor<rank, dim, T> &t)
940  {
941  return t;
942  }
943 
945  value(const T &t)
946  {
948  tmp = t;
949  return tmp;
950  }
951  };
952 } // namespace internal
953 
954 
955 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
956 
957 
958 template <int dim, typename Number>
961  // Some auto-differentiable numbers need explicit
962  // zero initialization such as adtl::adouble.
963  : Tensor{0.0}
964 {}
965 
966 
967 
968 template <int dim, typename Number>
969 template <typename OtherNumber>
971 Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
972  : value(internal::NumberType<Number>::value(initializer))
973 {}
974 
975 
976 
977 template <int dim, typename Number>
978 template <typename OtherNumber>
981  : Tensor{p.value}
982 {}
983 
984 
985 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
986 template <int dim, typename Number>
989  : value{other.value}
990 {}
991 
992 
993 
994 template <int dim, typename Number>
997  : value{std::move(other.value)}
998 {}
999 # endif
1000 
1001 
1002 template <int dim, typename Number>
1003 inline Number *
1005 {
1006  return std::addressof(value);
1007 }
1008 
1009 
1010 
1011 template <int dim, typename Number>
1012 inline const Number *
1014 {
1015  return std::addressof(value);
1016 }
1017 
1018 
1019 
1020 template <int dim, typename Number>
1021 inline Number *
1023 {
1025 }
1026 
1027 
1028 
1029 template <int dim, typename Number>
1030 const Number *
1032 {
1034 }
1035 
1036 
1037 
1038 template <int dim, typename Number>
1041 {
1042  // We cannot use Assert inside a CUDA kernel
1043 # ifndef __CUDA_ARCH__
1044  Assert(dim != 0,
1045  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1046 # endif
1047  return value;
1048 }
1049 
1050 
1051 template <int dim, typename Number>
1052 constexpr inline DEAL_II_ALWAYS_INLINE
1054 {
1055  // We cannot use Assert inside a CUDA kernel
1056 # ifndef __CUDA_ARCH__
1057  Assert(dim != 0,
1058  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1059 # endif
1060  return value;
1061 }
1062 
1063 
1064 
1065 template <int dim, typename Number>
1066 template <typename OtherNumber>
1069 {
1071  return *this;
1072 }
1073 
1074 
1075 # if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1076 template <int dim, typename Number>
1079 {
1080  value = p.value;
1081  return *this;
1082 }
1083 # endif
1084 
1085 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1086 template <int dim, typename Number>
1089 {
1090  value = std::move(other.value);
1091  return *this;
1092 }
1093 # endif
1094 
1095 
1096 
1097 template <int dim, typename Number>
1098 template <typename OtherNumber>
1100 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
1101 {
1103  return *this;
1104 }
1105 
1106 
1107 template <int dim, typename Number>
1108 template <typename OtherNumber>
1109 constexpr inline bool
1111 {
1112 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1113  Assert(!(std::is_same<Number, adouble>::value ||
1114  std::is_same<OtherNumber, adouble>::value),
1115  ExcMessage(
1116  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1117  "been extended to support advanced branching."));
1118 # endif
1119 
1120  return numbers::values_are_equal(value, p.value);
1121 }
1122 
1123 
1124 template <int dim, typename Number>
1125 template <typename OtherNumber>
1126 constexpr bool
1128 {
1129  return !((*this) == p);
1130 }
1131 
1132 
1133 template <int dim, typename Number>
1134 template <typename OtherNumber>
1137 {
1138  value += p.value;
1139  return *this;
1140 }
1141 
1142 
1143 template <int dim, typename Number>
1144 template <typename OtherNumber>
1147 {
1148  value -= p.value;
1149  return *this;
1150 }
1151 
1152 
1153 
1154 namespace internal
1155 {
1156  namespace ComplexWorkaround
1157  {
1158  template <typename Number, typename OtherNumber>
1159  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1160  multiply_assign_scalar(Number &val, const OtherNumber &s)
1161  {
1162  val *= s;
1163  }
1164 
1165 # ifdef __CUDA_ARCH__
1166  template <typename Number, typename OtherNumber>
1167  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1168  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1169  {
1170  printf("This function is not implemented for std::complex<Number>!\n");
1171  assert(false);
1172  }
1173 # endif
1174  } // namespace ComplexWorkaround
1175 } // namespace internal
1176 
1177 
1178 template <int dim, typename Number>
1179 template <typename OtherNumber>
1181 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1182 {
1183  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1184  return *this;
1185 }
1186 
1187 
1188 
1189 template <int dim, typename Number>
1190 template <typename OtherNumber>
1192 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1193 {
1194  value /= s;
1195  return *this;
1196 }
1197 
1198 
1199 template <int dim, typename Number>
1202 {
1203  return -value;
1204 }
1205 
1206 
1207 template <int dim, typename Number>
1210 {
1211  Assert(dim != 0,
1212  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1213  return numbers::NumberTraits<Number>::abs(value);
1214 }
1215 
1216 
1217 template <int dim, typename Number>
1221 {
1222  // We cannot use Assert inside a CUDA kernel
1223 # ifndef __CUDA_ARCH__
1224  Assert(dim != 0,
1225  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1226 # endif
1228 }
1229 
1230 
1231 
1232 template <int dim, typename Number>
1233 template <typename Iterator>
1234 Iterator
1235 Tensor<0, dim, Number>::unroll_recursion(const Iterator current,
1236  const Iterator end) const
1237 {
1238  (void)end;
1239  Assert(dim != 0,
1240  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1241  Assert(std::distance(current, end) >= 1,
1242  ExcMessage("The provided iterator range must contain at least one "
1243  "element."));
1244  *current = value;
1245  return std::next(current);
1246 }
1247 
1248 
1249 
1250 template <int dim, typename Number>
1251 constexpr inline void
1253 {
1254  // Some auto-differentiable numbers need explicit
1255  // zero initialization.
1257 }
1258 
1259 
1260 
1261 template <int dim, typename Number>
1262 template <class Iterator>
1263 inline void
1264 Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1265 {
1268 }
1269 
1270 
1271 
1272 template <int dim, typename Number>
1273 template <class Archive>
1274 inline void
1275 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1276 {
1277  ar &value;
1278 }
1279 
1280 
1281 template <int dim, typename Number>
1283 
1284 
1285 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1286 
1287 template <int rank_, int dim, typename Number>
1288 template <typename ArrayLike, std::size_t... indices>
1290 Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1291  std::index_sequence<indices...>)
1292  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1293 {
1294  static_assert(sizeof...(indices) == dim,
1295  "dim should match the number of indices");
1296 }
1297 
1298 
1299 
1300 template <int rank_, int dim, typename Number>
1303  // We would like to use =default, but this causes compile errors with some
1304  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1305  : values{}
1306 {}
1307 
1308 
1309 
1310 template <int rank_, int dim, typename Number>
1312 Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1313  : Tensor(initializer, std::make_index_sequence<dim>{})
1314 {}
1315 
1316 
1317 
1318 template <int rank_, int dim, typename Number>
1319 template <typename ElementType, typename MemorySpace>
1322  const ArrayView<ElementType, MemorySpace> &initializer)
1323 {
1324  // We cannot use Assert in a CUDA kernel
1325 # ifndef __CUDA_ARCH__
1326  AssertDimension(initializer.size(), n_independent_components);
1327 # endif
1328 
1329  for (unsigned int i = 0; i < n_independent_components; ++i)
1330  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1331 }
1332 
1333 
1334 
1335 template <int rank_, int dim, typename Number>
1336 template <typename OtherNumber>
1339  const Tensor<rank_, dim, OtherNumber> &initializer)
1340  : Tensor(initializer, std::make_index_sequence<dim>{})
1341 {}
1342 
1343 
1344 
1345 template <int rank_, int dim, typename Number>
1346 template <typename OtherNumber>
1347 constexpr DEAL_II_ALWAYS_INLINE
1349  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1350  : Tensor(initializer, std::make_index_sequence<dim>{})
1351 {}
1352 
1353 
1354 
1355 template <int rank_, int dim, typename Number>
1356 template <typename OtherNumber>
1358 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1359 {
1360  return Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>(values);
1361 }
1362 
1363 
1364 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1365 template <int rank_, int dim, typename Number>
1366 constexpr DEAL_II_ALWAYS_INLINE
1368 {
1369  for (unsigned int i = 0; i < dim; ++i)
1370  values[i] = other.values[i];
1371 }
1372 
1373 
1374 
1375 template <int rank_, int dim, typename Number>
1376 constexpr DEAL_II_ALWAYS_INLINE
1378 {
1379  for (unsigned int i = 0; i < dim; ++i)
1380  values[i] = other.values[i];
1381 }
1382 # endif
1383 
1384 namespace internal
1385 {
1386  namespace TensorSubscriptor
1387  {
1388  template <typename ArrayElementType, int dim>
1389  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
1390  subscript(ArrayElementType * values,
1391  const unsigned int i,
1392  std::integral_constant<int, dim>)
1393  {
1394  // We cannot use Assert in a CUDA kernel
1395 # ifndef __CUDA_ARCH__
1396  AssertIndexRange(i, dim);
1397 # endif
1398  return values[i];
1399  }
1400 
1401  template <typename ArrayElementType>
1402  constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE ArrayElementType &
1403  subscript(ArrayElementType *dummy,
1404  const unsigned int,
1405  std::integral_constant<int, 0>)
1406  {
1407  // We cannot use Assert in a CUDA kernel
1408 # ifndef __CUDA_ARCH__
1409  Assert(
1410  false,
1411  ExcMessage(
1412  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1413 # endif
1414  return *dummy;
1415  }
1416  } // namespace TensorSubscriptor
1417 } // namespace internal
1418 
1419 
1420 template <int rank_, int dim, typename Number>
1423  Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1424 {
1425  return ::internal::TensorSubscriptor::subscript(
1426  values, i, std::integral_constant<int, dim>());
1427 }
1428 
1429 
1430 template <int rank_, int dim, typename Number>
1431 constexpr DEAL_II_ALWAYS_INLINE
1433  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1434 {
1435 # if KOKKOS_VERSION < 30700
1436 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1437  AssertIndexRange(i, dim);
1438 # endif
1439 # else
1440  KOKKOS_IF_ON_HOST((AssertIndexRange(i, dim);))
1441 # endif
1442 
1443  return values[i];
1444 }
1445 
1446 
1447 template <int rank_, int dim, typename Number>
1448 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1450 {
1451 # if KOKKOS_VERSION < 30700
1452 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1453  Assert(dim != 0,
1454  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1455 # endif
1456 # else
1457  KOKKOS_IF_ON_HOST(
1458  (Assert(dim != 0,
1459  ExcMessage(
1460  "Cannot access an object of type Tensor<rank_,0,Number>"));))
1461 # endif
1462 
1463  return TensorAccessors::extract<rank_>(*this, indices);
1464 }
1465 
1466 
1467 
1468 template <int rank_, int dim, typename Number>
1469 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1471 {
1472 # if KOKKOS_VERSION < 30700
1473 # ifdef KOKKOS_ACTIVE_MEMORY_SPACE_HOST
1474  Assert(dim != 0,
1475  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1476 # endif
1477 # else
1478  KOKKOS_IF_ON_HOST(
1479  (Assert(dim != 0,
1480  ExcMessage(
1481  "Cannot access an object of type Tensor<rank_,0,Number>"));))
1482 # endif
1483 
1484  return TensorAccessors::extract<rank_>(*this, indices);
1485 }
1486 
1487 
1488 
1489 template <int rank_, int dim, typename Number>
1490 inline Number *
1492 {
1493  return std::addressof(
1494  this->operator[](this->unrolled_to_component_indices(0)));
1495 }
1496 
1497 
1498 
1499 template <int rank_, int dim, typename Number>
1500 inline const Number *
1502 {
1503  return std::addressof(
1504  this->operator[](this->unrolled_to_component_indices(0)));
1505 }
1506 
1507 
1508 
1509 template <int rank_, int dim, typename Number>
1510 inline Number *
1512 {
1513  return begin_raw() + n_independent_components;
1514 }
1515 
1516 
1517 
1518 template <int rank_, int dim, typename Number>
1519 inline const Number *
1521 {
1522  return begin_raw() + n_independent_components;
1523 }
1524 
1525 
1526 
1527 template <int rank_, int dim, typename Number>
1528 template <typename OtherNumber>
1531 {
1532  // The following loop could be written more concisely using std::copy, but
1533  // that function is only constexpr from C++20 on.
1534  for (unsigned int i = 0; i < dim; ++i)
1535  values[i] = t.values[i];
1536  return *this;
1537 }
1538 
1539 
1540 
1541 template <int rank_, int dim, typename Number>
1544 {
1546  (void)d;
1547 
1548  for (unsigned int i = 0; i < dim; ++i)
1550  return *this;
1551 }
1552 
1553 
1554 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1555 template <int rank_, int dim, typename Number>
1558 {
1559  for (unsigned int i = 0; i < dim; ++i)
1560  values[i] = other.values[i];
1561  return *this;
1562 }
1563 
1564 
1565 
1566 template <int rank_, int dim, typename Number>
1569  Tensor<rank_, dim, Number> &&other) noexcept
1570 {
1571  for (unsigned int i = 0; i < dim; ++i)
1572  values[i] = other.values[i];
1573  return *this;
1574 }
1575 # endif
1576 
1577 
1578 template <int rank_, int dim, typename Number>
1579 template <typename OtherNumber>
1580 constexpr inline bool
1582  const Tensor<rank_, dim, OtherNumber> &p) const
1583 {
1584  for (unsigned int i = 0; i < dim; ++i)
1585  if (values[i] != p.values[i])
1586  return false;
1587  return true;
1588 }
1589 
1590 
1591 // At some places in the library, we have Point<0> for formal reasons
1592 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1593 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1594 // in the above function that the loop end check always fails, we
1595 // implement this function here
1596 template <>
1597 template <>
1598 constexpr inline bool
1600 {
1601  return true;
1602 }
1603 
1604 
1605 template <int rank_, int dim, typename Number>
1606 template <typename OtherNumber>
1607 constexpr bool
1609  const Tensor<rank_, dim, OtherNumber> &p) const
1610 {
1611  return !((*this) == p);
1612 }
1613 
1614 
1615 template <int rank_, int dim, typename Number>
1616 template <typename OtherNumber>
1617 constexpr inline DEAL_II_ALWAYS_INLINE
1621 {
1622  for (unsigned int i = 0; i < dim; ++i)
1623  values[i] += p.values[i];
1624  return *this;
1625 }
1626 
1627 
1628 template <int rank_, int dim, typename Number>
1629 template <typename OtherNumber>
1630 constexpr inline DEAL_II_ALWAYS_INLINE
1634 {
1635  for (unsigned int i = 0; i < dim; ++i)
1636  values[i] -= p.values[i];
1637  return *this;
1638 }
1639 
1640 
1641 template <int rank_, int dim, typename Number>
1642 template <typename OtherNumber>
1643 constexpr inline DEAL_II_ALWAYS_INLINE
1645  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1646 {
1647  for (unsigned int i = 0; i < dim; ++i)
1648  values[i] *= s;
1649  return *this;
1650 }
1651 
1652 
1653 namespace internal
1654 {
1655  namespace TensorImplementation
1656  {
1657  template <int rank,
1658  int dim,
1659  typename Number,
1660  typename OtherNumber,
1661  std::enable_if_t<
1662  !std::is_integral<
1663  typename ProductType<Number, OtherNumber>::type>::value &&
1664  !std::is_same<Number, Differentiation::SD::Expression>::value,
1665  int> = 0>
1666  constexpr DEAL_II_HOST_DEVICE inline DEAL_II_ALWAYS_INLINE void
1668  const OtherNumber &factor)
1669  {
1670  const Number inverse_factor = Number(1.) / factor;
1671  // recurse over the base objects
1672  for (unsigned int d = 0; d < dim; ++d)
1673  t[d] *= inverse_factor;
1674  }
1675 
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  // recurse over the base objects
1691  for (unsigned int d = 0; d < dim; ++d)
1692  t[d] /= factor;
1693  }
1694  } // namespace TensorImplementation
1695 } // namespace internal
1696 
1697 
1698 template <int rank_, int dim, typename Number>
1699 template <typename OtherNumber>
1700 constexpr inline DEAL_II_ALWAYS_INLINE
1702  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1703 {
1705  return *this;
1706 }
1707 
1708 
1709 template <int rank_, int dim, typename Number>
1710 constexpr inline DEAL_II_ALWAYS_INLINE
1713 {
1715 
1716  for (unsigned int i = 0; i < dim; ++i)
1717  tmp.values[i] = -values[i];
1718 
1719  return tmp;
1720 }
1721 
1722 
1723 template <int rank_, int dim, typename Number>
1726 {
1727  // Make things work with AD types
1728  using std::sqrt;
1729  return sqrt(norm_square());
1730 }
1731 
1732 
1733 template <int rank_, int dim, typename Number>
1737 {
1739  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1740  for (unsigned int i = 0; i < dim; ++i)
1741  s += values[i].norm_square();
1742 
1743  return s;
1744 }
1745 
1746 
1747 
1748 template <int rank_, int dim, typename Number>
1749 template <typename OtherNumber>
1750 inline void
1752 {
1753  unroll(result.begin(), result.end());
1754 }
1755 
1756 
1757 
1758 template <int rank_, int dim, typename Number>
1759 template <class Iterator>
1760 inline void
1762  const Iterator end) const
1763 {
1764  AssertDimension(std::distance(begin, end), n_independent_components);
1765  unroll_recursion(begin, end);
1766 }
1767 
1768 
1769 
1770 template <int rank_, int dim, typename Number>
1771 template <typename Iterator>
1772 Iterator
1773 Tensor<rank_, dim, Number>::unroll_recursion(const Iterator current,
1774  const Iterator end) const
1775 {
1776  auto next = current;
1777  for (unsigned int i = 0; i < dim; ++i)
1778  next = values[i].unroll_recursion(next, end);
1779  return next;
1780 }
1781 
1782 
1783 template <int rank_, int dim, typename Number>
1784 constexpr inline unsigned int
1786  const TableIndices<rank_> &indices)
1787 {
1788  unsigned int index = 0;
1789  for (int r = 0; r < rank_; ++r)
1790  index = index * dim + indices[r];
1791 
1792  return index;
1793 }
1794 
1795 
1796 
1797 namespace internal
1798 {
1799  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1800  // and rank=2. Make sure we don't have compiler warnings.
1801 
1802  template <int dim>
1803  inline constexpr unsigned int
1804  mod(const unsigned int x)
1805  {
1806  return x % dim;
1807  }
1808 
1809  template <>
1810  inline unsigned int
1811  mod<0>(const unsigned int x)
1812  {
1813  Assert(false, ExcInternalError());
1814  return x;
1815  }
1816 
1817  template <int dim>
1818  inline constexpr unsigned int
1819  div(const unsigned int x)
1820  {
1821  return x / dim;
1822  }
1823 
1824  template <>
1825  inline unsigned int
1826  div<0>(const unsigned int x)
1827  {
1828  Assert(false, ExcInternalError());
1829  return x;
1830  }
1831 
1832 } // namespace internal
1833 
1834 
1835 
1836 template <int rank_, int dim, typename Number>
1837 constexpr inline TableIndices<rank_>
1839 {
1840  AssertIndexRange(i, n_independent_components);
1841 
1842  TableIndices<rank_> indices;
1843 
1844  unsigned int remainder = i;
1845  for (int r = rank_ - 1; r >= 0; --r)
1846  {
1847  indices[r] = internal::mod<dim>(remainder);
1848  remainder = internal::div<dim>(remainder);
1849  }
1850  Assert(remainder == 0, ExcInternalError());
1851 
1852  return indices;
1853 }
1854 
1855 
1856 template <int rank_, int dim, typename Number>
1857 constexpr inline void
1859 {
1860  for (unsigned int i = 0; i < dim; ++i)
1862 }
1863 
1864 
1865 template <int rank_, int dim, typename Number>
1866 constexpr std::size_t
1868 {
1869  return sizeof(Tensor<rank_, dim, Number>);
1870 }
1871 
1872 
1873 template <int rank_, int dim, typename Number>
1874 template <class Archive>
1875 inline void
1876 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1877 {
1878  ar &values;
1879 }
1880 
1881 
1882 template <int rank_, int dim, typename Number>
1884 
1885 #endif // DOXYGEN
1886 
1887 /* ----------------- Non-member functions operating on tensors. ------------ */
1888 
1901 template <int rank_, int dim, typename Number>
1902 inline std::ostream &
1903 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1904 {
1905  for (unsigned int i = 0; i < dim; ++i)
1906  {
1907  out << p[i];
1908  if (i != dim - 1)
1909  out << ' ';
1910  }
1911 
1912  return out;
1913 }
1914 
1915 
1922 template <int dim, typename Number>
1923 inline std::ostream &
1924 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1925 {
1926  out << static_cast<const Number &>(p);
1927  return out;
1928 }
1929 
1930 
1949 template <int dim, typename Number, typename Other>
1952  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1953 {
1954  return object * static_cast<const Number &>(t);
1955 }
1956 
1957 
1958 
1969 template <int dim, typename Number, typename Other>
1972  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1973 {
1974  return static_cast<const Number &>(t) * object;
1975 }
1976 
1977 
1989 template <int dim, typename Number, typename OtherNumber>
1993  const Tensor<0, dim, OtherNumber> &src2)
1994 {
1995  return static_cast<const Number &>(src1) *
1996  static_cast<const OtherNumber &>(src2);
1997 }
1998 
1999 
2007 template <int dim, typename Number, typename OtherNumber>
2009  Tensor<0,
2010  dim,
2011  typename ProductType<Number,
2012  typename EnableIfScalar<OtherNumber>::type>::type>
2013  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2014 {
2015  return static_cast<const Number &>(t) / factor;
2016 }
2017 
2018 
2026 template <int dim, typename Number, typename OtherNumber>
2030  const Tensor<0, dim, OtherNumber> &q)
2031 {
2032  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2033 }
2034 
2035 
2043 template <int dim, typename Number, typename OtherNumber>
2047  const Tensor<0, dim, OtherNumber> &q)
2048 {
2049  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2050 }
2051 
2052 
2065 template <int rank, int dim, typename Number, typename OtherNumber>
2067  Tensor<rank,
2068  dim,
2069  typename ProductType<Number,
2070  typename EnableIfScalar<OtherNumber>::type>::type>
2071  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2072 {
2073  // recurse over the base objects
2075  for (unsigned int d = 0; d < dim; ++d)
2076  tt[d] = t[d] * factor;
2077  return tt;
2078 }
2079 
2080 
2093 template <int rank, int dim, typename Number, typename OtherNumber>
2095  Tensor<rank,
2096  dim,
2098  OtherNumber>::type>
2099  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2100 {
2101  // simply forward to the operator above
2102  return t * factor;
2103 }
2104 
2105 
2106 namespace internal
2107 {
2108  namespace TensorImplementation
2109  {
2110  template <int rank,
2111  int dim,
2112  typename Number,
2113  typename OtherNumber,
2114  std::enable_if_t<
2115  !std::is_integral<
2116  typename ProductType<Number, OtherNumber>::type>::value,
2117  int> = 0>
2121  const OtherNumber & factor)
2122  {
2124  const Number inverse_factor = Number(1.) / factor;
2125  // recurse over the base objects
2126  for (unsigned int d = 0; d < dim; ++d)
2127  tt[d] = t[d] * inverse_factor;
2128  return tt;
2129  }
2130 
2131 
2132  template <int rank,
2133  int dim,
2134  typename Number,
2135  typename OtherNumber,
2136  std::enable_if_t<
2137  std::is_integral<
2138  typename ProductType<Number, OtherNumber>::type>::value,
2139  int> = 0>
2143  const OtherNumber & factor)
2144  {
2146  // recurse over the base objects
2147  for (unsigned int d = 0; d < dim; ++d)
2148  tt[d] = t[d] / factor;
2149  return tt;
2150  }
2151  } // namespace TensorImplementation
2152 } // namespace internal
2153 
2154 
2164 template <int rank, int dim, typename Number, typename OtherNumber>
2166  Tensor<rank,
2167  dim,
2168  typename ProductType<Number,
2169  typename EnableIfScalar<OtherNumber>::type>::type>
2170  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2171 {
2173 }
2174 
2175 
2185 template <int rank, int dim, typename Number, typename OtherNumber>
2190 {
2192 
2193  for (unsigned int i = 0; i < dim; ++i)
2194  tmp[i] += q[i];
2195 
2196  return tmp;
2197 }
2198 
2199 
2209 template <int rank, int dim, typename Number, typename OtherNumber>
2214 {
2216 
2217  for (unsigned int i = 0; i < dim; ++i)
2218  tmp[i] -= q[i];
2219 
2220  return tmp;
2221 }
2222 
2229 template <int dim, typename Number, typename OtherNumber>
2230 inline constexpr DEAL_II_ALWAYS_INLINE
2233  const Tensor<0, dim, OtherNumber> &src2)
2234 {
2236 
2237  tmp *= src2;
2238 
2239  return tmp;
2240 }
2241 
2258 template <int rank, int dim, typename Number, typename OtherNumber>
2259 inline constexpr DEAL_II_ALWAYS_INLINE
2262  const Tensor<rank, dim, OtherNumber> &src2)
2263 {
2265 
2266  for (unsigned int i = 0; i < dim; ++i)
2267  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2269 
2270  return tmp;
2271 }
2272 
2317 template <int rank_1,
2318  int rank_2,
2319  int dim,
2320  typename Number,
2321  typename OtherNumber,
2322  typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2323 constexpr inline DEAL_II_ALWAYS_INLINE
2324  typename Tensor<rank_1 + rank_2 - 2,
2325  dim,
2326  typename ProductType<Number, OtherNumber>::type>::tensor_type
2329 {
2330  typename Tensor<rank_1 + rank_2 - 2,
2331  dim,
2332  typename ProductType<Number, OtherNumber>::type>::tensor_type
2333  result{};
2334 
2335  TensorAccessors::internal::
2336  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2337  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2338  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2339 
2340  return result;
2341 }
2342 
2343 
2372 template <int index_1,
2373  int index_2,
2374  int rank_1,
2375  int rank_2,
2376  int dim,
2377  typename Number,
2378  typename OtherNumber>
2379 constexpr inline DEAL_II_ALWAYS_INLINE
2380  typename Tensor<rank_1 + rank_2 - 2,
2381  dim,
2382  typename ProductType<Number, OtherNumber>::type>::tensor_type
2385 {
2386  Assert(0 <= index_1 && index_1 < rank_1,
2387  ExcMessage(
2388  "The specified index_1 must lie within the range [0,rank_1)"));
2389  Assert(0 <= index_2 && index_2 < rank_2,
2390  ExcMessage(
2391  "The specified index_2 must lie within the range [0,rank_2)"));
2392 
2393  using namespace TensorAccessors;
2394  using namespace TensorAccessors::internal;
2395 
2396  // Reorder index_1 to the end of src1:
2398  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2399 
2400  // Reorder index_2 to the end of src2:
2401  const ReorderedIndexView<index_2,
2402  rank_2,
2404  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2405 
2406  typename Tensor<rank_1 + rank_2 - 2,
2407  dim,
2408  typename ProductType<Number, OtherNumber>::type>::tensor_type
2409  result{};
2410  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2411  return result;
2412 }
2413 
2414 
2445 template <int index_1,
2446  int index_2,
2447  int index_3,
2448  int index_4,
2449  int rank_1,
2450  int rank_2,
2451  int dim,
2452  typename Number,
2453  typename OtherNumber>
2454 constexpr inline
2455  typename Tensor<rank_1 + rank_2 - 4,
2456  dim,
2457  typename ProductType<Number, OtherNumber>::type>::tensor_type
2458  double_contract(const Tensor<rank_1, dim, Number> & src1,
2460 {
2461  Assert(0 <= index_1 && index_1 < rank_1,
2462  ExcMessage(
2463  "The specified index_1 must lie within the range [0,rank_1)"));
2464  Assert(0 <= index_3 && index_3 < rank_1,
2465  ExcMessage(
2466  "The specified index_3 must lie within the range [0,rank_1)"));
2467  Assert(index_1 != index_3,
2468  ExcMessage("index_1 and index_3 must not be the same"));
2469  Assert(0 <= index_2 && index_2 < rank_2,
2470  ExcMessage(
2471  "The specified index_2 must lie within the range [0,rank_2)"));
2472  Assert(0 <= index_4 && index_4 < rank_2,
2473  ExcMessage(
2474  "The specified index_4 must lie within the range [0,rank_2)"));
2475  Assert(index_2 != index_4,
2476  ExcMessage("index_2 and index_4 must not be the same"));
2477 
2478  using namespace TensorAccessors;
2479  using namespace TensorAccessors::internal;
2480 
2481  // Reorder index_1 to the end of src1:
2483  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2484 
2485  // Reorder index_2 to the end of src2:
2487  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2488 
2489  // Now, reorder index_3 to the end of src1. We have to make sure to
2490  // preserve the original ordering: index_1 has been removed. If
2491  // index_3 > index_1, we have to use (index_3 - 1) instead:
2493  (index_3 < index_1 ? index_3 : index_3 - 1),
2494  rank_1,
2495  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2496  reord_3 =
2497  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2498  index_3 - 1,
2499  rank_1 > (reord_1);
2500 
2501  // Now, reorder index_4 to the end of src2. We have to make sure to
2502  // preserve the original ordering: index_2 has been removed. If
2503  // index_4 > index_2, we have to use (index_4 - 1) instead:
2505  (index_4 < index_2 ? index_4 : index_4 - 1),
2506  rank_2,
2507  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2508  reord_4 =
2509  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2510  index_4 - 1,
2511  rank_2 > (reord_2);
2512 
2513  typename Tensor<rank_1 + rank_2 - 4,
2514  dim,
2515  typename ProductType<Number, OtherNumber>::type>::tensor_type
2516  result{};
2517  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2518  return result;
2519 }
2520 
2521 
2534 template <int rank, int dim, typename Number, typename OtherNumber>
2535 constexpr inline DEAL_II_ALWAYS_INLINE
2537  scalar_product(const Tensor<rank, dim, Number> & left,
2538  const Tensor<rank, dim, OtherNumber> &right)
2539 {
2540  typename ProductType<Number, OtherNumber>::type result{};
2541  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2542  return result;
2543 }
2544 
2545 
2563 template <template <int, int, typename> class TensorT1,
2564  template <int, int, typename>
2565  class TensorT2,
2566  template <int, int, typename>
2567  class TensorT3,
2568  int rank_1,
2569  int rank_2,
2570  int dim,
2571  typename T1,
2572  typename T2,
2573  typename T3>
2574 constexpr inline DEAL_II_ALWAYS_INLINE
2576  contract3(const TensorT1<rank_1, dim, T1> & left,
2577  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2578  const TensorT3<rank_2, dim, T3> & right)
2579 {
2580  using return_type =
2582  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2583  middle,
2584  right);
2585 }
2586 
2587 
2598 template <int rank_1,
2599  int rank_2,
2600  int dim,
2601  typename Number,
2602  typename OtherNumber>
2603 constexpr inline DEAL_II_ALWAYS_INLINE
2607 {
2608  typename Tensor<rank_1 + rank_2,
2609  dim,
2610  typename ProductType<Number, OtherNumber>::type>::tensor_type
2611  result{};
2612  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2613  return result;
2614 }
2615 
2616 
2635 template <int dim, typename Number>
2637 cross_product_2d(const Tensor<1, dim, Number> &src)
2638 {
2639  Assert(dim == 2, ExcInternalError());
2640 
2641  Tensor<1, dim, Number> result;
2642 
2643  result[0] = src[1];
2644  result[1] = -src[0];
2645 
2646  return result;
2647 }
2648 
2649 
2659 template <int dim, typename Number1, typename Number2>
2660 constexpr inline DEAL_II_ALWAYS_INLINE
2662  cross_product_3d(const Tensor<1, dim, Number1> &src1,
2663  const Tensor<1, dim, Number2> &src2)
2664 {
2665  Assert(dim == 3, ExcInternalError());
2666 
2668 
2669  // avoid compiler warnings
2670  constexpr int s0 = 0 % dim;
2671  constexpr int s1 = 1 % dim;
2672  constexpr int s2 = 2 % dim;
2673 
2674  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2675  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2676  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2677 
2678  return result;
2679 }
2680 
2681 
2695 template <int dim, typename Number>
2696 constexpr inline DEAL_II_ALWAYS_INLINE Number
2698 {
2699  // Compute the determinant using the Laplace expansion of the
2700  // determinant. We expand along the last row.
2701  Number det = internal::NumberType<Number>::value(0.0);
2702 
2703  for (unsigned int k = 0; k < dim; ++k)
2704  {
2705  Tensor<2, dim - 1, Number> minor;
2706  for (unsigned int i = 0; i < dim - 1; ++i)
2707  for (unsigned int j = 0; j < dim - 1; ++j)
2708  minor[i][j] = t[i][j < k ? j : j + 1];
2709 
2710  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2711 
2712  det += t[dim - 1][k] * cofactor;
2713  }
2714 
2715  return ((dim % 2 == 0) ? 1. : -1.) * det;
2716 }
2717 
2723 template <typename Number>
2724 constexpr DEAL_II_ALWAYS_INLINE Number
2726 {
2727  return t[0][0];
2728 }
2729 
2735 template <typename Number>
2736 constexpr DEAL_II_ALWAYS_INLINE Number
2738 {
2739  // hard-coded for efficiency reasons
2740  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2741 }
2742 
2748 template <typename Number>
2749 constexpr DEAL_II_ALWAYS_INLINE Number
2751 {
2752  // hard-coded for efficiency reasons
2753  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2754  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2755  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2756  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2757  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2758  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2759  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2760 }
2761 
2762 
2769 template <int dim, typename Number>
2770 constexpr inline DEAL_II_ALWAYS_INLINE Number
2772 {
2773  Number t = d[0][0];
2774  for (unsigned int i = 1; i < dim; ++i)
2775  t += d[i][i];
2776  return t;
2777 }
2778 
2779 
2788 template <int dim, typename Number>
2789 constexpr inline Tensor<2, dim, Number>
2791 {
2792  Number return_tensor[dim][dim];
2793 
2794  // if desired, take over the
2795  // inversion of a 4x4 tensor
2796  // from the FullMatrix
2797  AssertThrow(false, ExcNotImplemented());
2798 
2799  return Tensor<2, dim, Number>(return_tensor);
2800 }
2801 
2802 
2803 #ifndef DOXYGEN
2804 
2805 template <typename Number>
2807  invert(const Tensor<2, 1, Number> &t)
2808 {
2809  Tensor<2, 1, Number> return_tensor;
2810 
2811  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2812 
2813  return return_tensor;
2814 }
2815 
2816 
2817 template <typename Number>
2819  invert(const Tensor<2, 2, Number> &t)
2820 {
2821  Tensor<2, 2, Number> return_tensor;
2822 
2823  const Number inv_det_t = internal::NumberType<Number>::value(
2824  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2825  return_tensor[0][0] = t[1][1];
2826  return_tensor[0][1] = -t[0][1];
2827  return_tensor[1][0] = -t[1][0];
2828  return_tensor[1][1] = t[0][0];
2829  return_tensor *= inv_det_t;
2830 
2831  return return_tensor;
2832 }
2833 
2834 
2835 template <typename Number>
2837  invert(const Tensor<2, 3, Number> &t)
2838 {
2839  Tensor<2, 3, Number> return_tensor;
2840 
2841  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2842  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2843  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2844  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2845  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2846  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2847  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2848  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2849  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2850  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2851  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2852  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2853  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2854  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2855  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2856  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2857  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2858  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2859  const Number inv_det_t = internal::NumberType<Number>::value(
2860  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2861  t[0][2] * return_tensor[2][0]));
2862  return_tensor *= inv_det_t;
2863 
2864  return return_tensor;
2865 }
2866 
2867 #endif /* DOXYGEN */
2868 
2869 
2875 template <int dim, typename Number>
2878 {
2880  for (unsigned int i = 0; i < dim; ++i)
2881  {
2882  tt[i][i] = t[i][i];
2883  for (unsigned int j = i + 1; j < dim; ++j)
2884  {
2885  tt[i][j] = t[j][i];
2886  tt[j][i] = t[i][j];
2887  };
2888  }
2889  return tt;
2890 }
2891 
2892 
2906 template <int dim, typename Number>
2907 constexpr Tensor<2, dim, Number>
2908 adjugate(const Tensor<2, dim, Number> &t)
2909 {
2910  return determinant(t) * invert(t);
2911 }
2912 
2913 
2927 template <int dim, typename Number>
2928 constexpr Tensor<2, dim, Number>
2929 cofactor(const Tensor<2, dim, Number> &t)
2930 {
2931  return transpose(adjugate(t));
2932 }
2933 
2934 
2998 template <int dim, typename Number>
3001 
3002 
3010 template <int dim, typename Number>
3011 inline Number
3013 {
3015  for (unsigned int j = 0; j < dim; ++j)
3016  {
3018  for (unsigned int i = 0; i < dim; ++i)
3020 
3021  if (sum > max)
3022  max = sum;
3023  }
3024 
3025  return max;
3026 }
3027 
3028 
3036 template <int dim, typename Number>
3037 inline Number
3039 {
3041  for (unsigned int i = 0; i < dim; ++i)
3042  {
3044  for (unsigned int j = 0; j < dim; ++j)
3046 
3047  if (sum > max)
3048  max = sum;
3049  }
3050 
3051  return max;
3052 }
3053 
3059 #ifndef DOXYGEN
3060 
3061 
3062 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3063 
3064 // Specialization of functions for ADOL-C number types when
3065 // the advanced branching feature is used
3066 template <int dim>
3067 inline adouble
3069 {
3071  for (unsigned int j = 0; j < dim; ++j)
3072  {
3074  for (unsigned int i = 0; i < dim; ++i)
3075  sum += fabs(t[i][j]);
3076 
3077  condassign(max, (sum > max), sum, max);
3078  }
3079 
3080  return max;
3081 }
3082 
3083 
3084 template <int dim>
3085 inline adouble
3087 {
3089  for (unsigned int i = 0; i < dim; ++i)
3090  {
3092  for (unsigned int j = 0; j < dim; ++j)
3093  sum += fabs(t[i][j]);
3094 
3095  condassign(max, (sum > max), sum, max);
3096  }
3097 
3098  return max;
3099 }
3100 
3101 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3102 
3103 
3104 #endif // DOXYGEN
3105 
3107 
3108 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
std::size_t size() const
Definition: array_view.h:566
Definition: point.h:111
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 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:122
constexpr Tensor & operator/=(const OtherNumber &factor)
constexpr Tensor operator-() const
constexpr Tensor & operator=(const OtherNumber &d)
Definition: tensor.h:503
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:2170
typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1] array_type
Definition: tensor.h:543
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:2013
static constexpr unsigned int rank
Definition: tensor.h:522
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:1992
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:2232
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:2046
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:536
friend class Tensor
Definition: tensor.h:887
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:2188
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3038
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1952
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3012
constexpr Tensor & operator/=(const OtherNumber &factor)
static constexpr unsigned int dimension
Definition: tensor.h:517
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:2099
static constexpr std::size_t memory_consumption()
Number * begin_raw()
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:2029
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:2261
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:2212
constexpr Tensor & operator=(const Number &d)
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:858
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:1972
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:2071
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor(const array_type &initializer)
static constexpr unsigned int n_independent_components
Definition: tensor.h:528
constexpr Tensor operator-() const
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:864
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:162
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
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:1509
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1703
#define AssertIndexRange(index, range)
Definition: exceptions.h:1768
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1619
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)
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > division_operator(const Tensor< rank, dim, Number > &t, const OtherNumber &factor)
Definition: tensor.h:2120
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:970
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:954
#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:733
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
Definition: numbers.h:624
static constexpr std::enable_if_t< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type > abs_square(const number &x)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
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:3038
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3012