Reference documentation for deal.II version GIT 6113e8d4d8 2022-09-29 15:30: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_CUDA_HOST_DEV
144 
152  template <typename OtherNumber>
153  constexpr DEAL_II_CUDA_HOST_DEV
154  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
155 
161  template <typename OtherNumber>
162  constexpr DEAL_II_CUDA_HOST_DEV
163  Tensor(const OtherNumber &initializer);
164 
165 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
169  constexpr DEAL_II_CUDA_HOST_DEV
170  Tensor(const Tensor<0, dim, Number> &other);
171 
175  constexpr DEAL_II_CUDA_HOST_DEV
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_CUDA_HOST_DEV
238  operator Number &();
239 
248  constexpr DEAL_II_CUDA_HOST_DEV operator const Number &() const;
249 
257  template <typename OtherNumber>
258  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
260 
261 #if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
270  constexpr DEAL_II_CUDA_HOST_DEV 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_CUDA_HOST_DEV 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_CUDA_HOST_DEV Tensor &
314 
320  template <typename OtherNumber>
321  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
323 
329  template <typename OtherNumber>
330  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
331  operator*=(const OtherNumber &factor);
332 
338  template <typename OtherNumber>
339  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
340  operator/=(const OtherNumber &factor);
341 
347  constexpr DEAL_II_CUDA_HOST_DEV 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_CUDA_HOST_DEV explicit Tensor(
559  const array_type &initializer);
560 
574  template <typename ElementType, typename MemorySpace>
575  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
576  const ArrayView<ElementType, MemorySpace> &initializer);
577 
585  template <typename OtherNumber>
586  constexpr DEAL_II_CUDA_HOST_DEV
588 
592  template <typename OtherNumber>
593  constexpr Tensor(
594  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
595 
599  template <typename OtherNumber>
600  constexpr
601  operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
602 
603 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
607  constexpr Tensor(const Tensor<rank_, dim, Number> &);
608 
612  constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
613 #endif
614 
621  operator[](const unsigned int i);
622 
628  constexpr DEAL_II_CUDA_HOST_DEV const value_type &
629  operator[](const unsigned int i) const;
630 
634  constexpr const Number &
635  operator[](const TableIndices<rank_> &indices) const;
636 
640  constexpr Number &
642 
646  Number *
648 
652  const Number *
653  begin_raw() const;
654 
658  Number *
660 
664  const Number *
665  end_raw() const;
666 
674  template <typename OtherNumber>
675  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
677 
684  constexpr Tensor &
685  operator=(const Number &d);
686 
687 #ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
691  constexpr Tensor<rank_, dim, Number> &
693 
697  constexpr Tensor<rank_, dim, Number> &
699 #endif
700 
704  template <typename OtherNumber>
705  constexpr bool
707 
711  template <typename OtherNumber>
712  constexpr bool
714 
720  template <typename OtherNumber>
721  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
723 
729  template <typename OtherNumber>
730  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
732 
739  template <typename OtherNumber>
740  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
741  operator*=(const OtherNumber &factor);
742 
748  template <typename OtherNumber>
749  constexpr DEAL_II_CUDA_HOST_DEV Tensor &
750  operator/=(const OtherNumber &factor);
751 
757  constexpr DEAL_II_CUDA_HOST_DEV Tensor
758  operator-() const;
759 
772  constexpr void
773  clear();
774 
784  norm() const;
785 
792  constexpr DEAL_II_CUDA_HOST_DEV
794  norm_square() const;
795 
806  template <typename OtherNumber>
807  DEAL_II_DEPRECATED void
808  unroll(Vector<OtherNumber> &result) const;
809 
820  template <class Iterator>
821  void
822  unroll(const Iterator begin, const Iterator end) const;
823 
828  static constexpr unsigned int
830 
836  static constexpr TableIndices<rank_>
837  unrolled_to_component_indices(const unsigned int i);
838 
843  static constexpr std::size_t
845 
851  template <class Archive>
852  void
853  serialize(Archive &ar, const unsigned int version);
854 
860 
861 private:
865  Tensor<rank_ - 1, dim, Number> values[(dim != 0) ? dim : 1];
866  // ... avoid a compiler warning in case of dim == 0 and ensure that the
867  // array always has positive size.
868 
872  template <typename Iterator>
873  Iterator
874  unroll_recursion(const Iterator current, const Iterator end) const;
875 
882  template <typename ArrayLike, std::size_t... Indices>
883  constexpr DEAL_II_CUDA_HOST_DEV
884  Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
885 
886  // Allow an arbitrary Tensor to access the underlying values.
887  template <int, int, typename>
888  friend class Tensor;
889 
890  // Point is allowed access to the coordinates. This is supposed to improve
891  // speed.
892  friend class Point<dim, Number>;
893 };
894 
895 
896 #ifndef DOXYGEN
897 namespace internal
898 {
899  // Workaround: The following 4 overloads are necessary to be able to
900  // compile the library with Apple Clang 8 and older. We should remove
901  // these overloads again when we bump the minimal required version to
902  // something later than clang-3.6 / Apple Clang 6.3.
903  template <int rank, int dim, typename T, typename U>
904  struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
905  {
906  using type =
908  };
909 
910  template <int rank, int dim, typename T, typename U>
911  struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
912  {
913  using type =
915  };
916 
917  template <typename T, int rank, int dim, typename U>
918  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
919  {
920  using type =
922  };
923 
924  template <int rank, int dim, typename T, typename U>
925  struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
926  {
927  using type =
929  };
930  // end workaround
931 
936  template <int rank, int dim, typename T>
937  struct NumberType<Tensor<rank, dim, T>>
938  {
939  static constexpr DEAL_II_ALWAYS_INLINE const Tensor<rank, dim, T> &
940  value(const Tensor<rank, dim, T> &t)
941  {
942  return t;
943  }
944 
946  value(const T &t)
947  {
949  tmp = t;
950  return tmp;
951  }
952  };
953 } // namespace internal
954 
955 
956 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
957 
958 
959 template <int dim, typename Number>
962  // Some auto-differentiable numbers need explicit
963  // zero initialization such as adtl::adouble.
964  : Tensor{0.0}
965 {}
966 
967 
968 
969 template <int dim, typename Number>
970 template <typename OtherNumber>
972 Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
973  : value(internal::NumberType<Number>::value(initializer))
974 {}
975 
976 
977 
978 template <int dim, typename Number>
979 template <typename OtherNumber>
982  : Tensor{p.value}
983 {}
984 
985 
986 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
987 template <int dim, typename Number>
990  : value{other.value}
991 {}
992 
993 
994 
995 template <int dim, typename Number>
998  : value{std::move(other.value)}
999 {}
1000 # endif
1001 
1002 
1003 template <int dim, typename Number>
1004 inline Number *
1006 {
1007  return std::addressof(value);
1008 }
1009 
1010 
1011 
1012 template <int dim, typename Number>
1013 inline const Number *
1015 {
1016  return std::addressof(value);
1017 }
1018 
1019 
1020 
1021 template <int dim, typename Number>
1022 inline Number *
1024 {
1026 }
1027 
1028 
1029 
1030 template <int dim, typename Number>
1031 const Number *
1033 {
1035 }
1036 
1037 
1038 
1039 template <int dim, typename Number>
1042 {
1043  // We cannot use Assert inside a CUDA kernel
1044 # ifndef __CUDA_ARCH__
1045  Assert(dim != 0,
1046  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1047 # endif
1048  return value;
1049 }
1050 
1051 
1052 template <int dim, typename Number>
1053 constexpr inline DEAL_II_ALWAYS_INLINE
1055 {
1056  // We cannot use Assert inside a CUDA kernel
1057 # ifndef __CUDA_ARCH__
1058  Assert(dim != 0,
1059  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1060 # endif
1061  return value;
1062 }
1063 
1064 
1065 
1066 template <int dim, typename Number>
1067 template <typename OtherNumber>
1068 constexpr inline DEAL_II_ALWAYS_INLINE
1071 {
1073  return *this;
1074 }
1075 
1076 
1077 # if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1078 template <int dim, typename Number>
1079 constexpr inline DEAL_II_ALWAYS_INLINE
1082 {
1083  value = p.value;
1084  return *this;
1085 }
1086 # endif
1087 
1088 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1089 template <int dim, typename Number>
1092 {
1093  value = std::move(other.value);
1094  return *this;
1095 }
1096 # endif
1097 
1098 
1099 
1100 template <int dim, typename Number>
1101 template <typename OtherNumber>
1102 constexpr inline DEAL_II_ALWAYS_INLINE
1104  Tensor<0, dim, Number>::operator=(const OtherNumber &d)
1105 {
1107  return *this;
1108 }
1109 
1110 
1111 template <int dim, typename Number>
1112 template <typename OtherNumber>
1113 constexpr inline bool
1115 {
1116 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1117  Assert(!(std::is_same<Number, adouble>::value ||
1118  std::is_same<OtherNumber, adouble>::value),
1119  ExcMessage(
1120  "The Tensor equality operator for ADOL-C taped numbers has not yet "
1121  "been extended to support advanced branching."));
1122 # endif
1123 
1124  return numbers::values_are_equal(value, p.value);
1125 }
1126 
1127 
1128 template <int dim, typename Number>
1129 template <typename OtherNumber>
1130 constexpr bool
1132 {
1133  return !((*this) == p);
1134 }
1135 
1136 
1137 template <int dim, typename Number>
1138 template <typename OtherNumber>
1139 constexpr inline DEAL_II_ALWAYS_INLINE
1142 {
1143  value += p.value;
1144  return *this;
1145 }
1146 
1147 
1148 template <int dim, typename Number>
1149 template <typename OtherNumber>
1150 constexpr inline DEAL_II_ALWAYS_INLINE
1153 {
1154  value -= p.value;
1155  return *this;
1156 }
1157 
1158 
1159 
1160 namespace internal
1161 {
1162  namespace ComplexWorkaround
1163  {
1164  template <typename Number, typename OtherNumber>
1165  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1166  multiply_assign_scalar(Number &val, const OtherNumber &s)
1167  {
1168  val *= s;
1169  }
1170 
1171 # ifdef __CUDA_ARCH__
1172  template <typename Number, typename OtherNumber>
1173  constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV void
1174  multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
1175  {
1176  printf("This function is not implemented for std::complex<Number>!\n");
1177  assert(false);
1178  }
1179 # endif
1180  } // namespace ComplexWorkaround
1181 } // namespace internal
1182 
1183 
1184 template <int dim, typename Number>
1185 template <typename OtherNumber>
1186 constexpr inline DEAL_II_ALWAYS_INLINE
1188  Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1189 {
1190  internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1191  return *this;
1192 }
1193 
1194 
1195 
1196 template <int dim, typename Number>
1197 template <typename OtherNumber>
1199 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1200 {
1201  value /= s;
1202  return *this;
1203 }
1204 
1205 
1206 template <int dim, typename Number>
1209 {
1210  return -value;
1211 }
1212 
1213 
1214 template <int dim, typename Number>
1217 {
1218  Assert(dim != 0,
1219  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1220  return numbers::NumberTraits<Number>::abs(value);
1221 }
1222 
1223 
1224 template <int dim, typename Number>
1228 {
1229  // We cannot use Assert inside a CUDA kernel
1230 # ifndef __CUDA_ARCH__
1231  Assert(dim != 0,
1232  ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1233 # endif
1235 }
1236 
1237 
1238 
1239 template <int dim, typename Number>
1240 template <typename Iterator>
1241 Iterator
1242 Tensor<0, dim, Number>::unroll_recursion(const Iterator current,
1243  const Iterator end) const
1244 {
1245  (void)end;
1246  Assert(dim != 0,
1247  ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1248  Assert(std::distance(current, end) >= 1,
1249  ExcMessage("The provided iterator range must contain at least one "
1250  "element."));
1251  *current = value;
1252  return std::next(current);
1253 }
1254 
1255 
1256 
1257 template <int dim, typename Number>
1258 constexpr inline void
1260 {
1261  // Some auto-differentiable numbers need explicit
1262  // zero initialization.
1264 }
1265 
1266 
1267 
1268 template <int dim, typename Number>
1269 template <class Iterator>
1270 inline void
1271 Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1272 {
1275 }
1276 
1277 
1278 
1279 template <int dim, typename Number>
1280 template <class Archive>
1281 inline void
1282 Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1283 {
1284  ar &value;
1285 }
1286 
1287 
1288 template <int dim, typename Number>
1290 
1291 
1292 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1293 
1294 template <int rank_, int dim, typename Number>
1295 template <typename ArrayLike, std::size_t... indices>
1297 Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1298  std::index_sequence<indices...>)
1299  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
1300 {
1301  static_assert(sizeof...(indices) == dim,
1302  "dim should match the number of indices");
1303 }
1304 
1305 
1306 
1307 template <int rank_, int dim, typename Number>
1310  // We would like to use =default, but this causes compile errors with some
1311  // MSVC versions and internal compiler errors with -O1 in gcc 5.4.
1312  : values{}
1313 {}
1314 
1315 
1316 
1317 template <int rank_, int dim, typename Number>
1319 Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1320  : Tensor(initializer, std::make_index_sequence<dim>{})
1321 {}
1322 
1323 
1324 
1325 template <int rank_, int dim, typename Number>
1326 template <typename ElementType, typename MemorySpace>
1329  const ArrayView<ElementType, MemorySpace> &initializer)
1330 {
1331  AssertDimension(initializer.size(), n_independent_components);
1332 
1333  for (unsigned int i = 0; i < n_independent_components; ++i)
1334  (*this)[unrolled_to_component_indices(i)] = initializer[i];
1335 }
1336 
1337 
1338 
1339 template <int rank_, int dim, typename Number>
1340 template <typename OtherNumber>
1343  const Tensor<rank_, dim, OtherNumber> &initializer)
1344  : Tensor(initializer, std::make_index_sequence<dim>{})
1345 {}
1346 
1347 
1348 
1349 template <int rank_, int dim, typename Number>
1350 template <typename OtherNumber>
1351 constexpr DEAL_II_ALWAYS_INLINE
1353  const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1354  : Tensor(initializer, std::make_index_sequence<dim>{})
1355 {}
1356 
1357 
1358 
1359 template <int rank_, int dim, typename Number>
1360 template <typename OtherNumber>
1362 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1363 {
1364  return Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>(values);
1365 }
1366 
1367 
1368 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1369 template <int rank_, int dim, typename Number>
1370 constexpr DEAL_II_ALWAYS_INLINE
1372 {
1373  for (unsigned int i = 0; i < dim; ++i)
1374  values[i] = other.values[i];
1375 }
1376 
1377 
1378 
1379 template <int rank_, int dim, typename Number>
1380 constexpr DEAL_II_ALWAYS_INLINE
1382 {
1383  for (unsigned int i = 0; i < dim; ++i)
1384  values[i] = other.values[i];
1385 }
1386 # endif
1387 
1388 
1389 namespace internal
1390 {
1391  namespace TensorSubscriptor
1392  {
1393  template <typename ArrayElementType, int dim>
1394  constexpr inline DEAL_II_ALWAYS_INLINE
1395  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1396  subscript(ArrayElementType * values,
1397  const unsigned int i,
1398  std::integral_constant<int, dim>)
1399  {
1400  // We cannot use Assert in a CUDA kernel
1401 # ifndef __CUDA_ARCH__
1402  AssertIndexRange(i, dim);
1403 # endif
1404  return values[i];
1405  }
1406 
1407  // The variables within this struct will be referenced in the next function.
1408  // It is a workaround that allows returning a reference to a static variable
1409  // while allowing constexpr evaluation of the function.
1410  // It has to be defined outside the function because constexpr functions
1411  // cannot define static variables
1412  template <typename ArrayElementType>
1413  struct Uninitialized
1414  {
1415  static ArrayElementType value;
1416  };
1417 
1418  template <typename Type>
1419  Type Uninitialized<Type>::value;
1420 
1421  template <typename ArrayElementType>
1422  constexpr inline DEAL_II_ALWAYS_INLINE
1423  DEAL_II_CUDA_HOST_DEV ArrayElementType &
1424  subscript(ArrayElementType *,
1425  const unsigned int,
1426  std::integral_constant<int, 0>)
1427  {
1428  // We cannot use Assert in a CUDA kernel
1429 # ifndef __CUDA_ARCH__
1430  Assert(
1431  false,
1432  ExcMessage(
1433  "Cannot access elements of an object of type Tensor<rank,0,Number>."));
1434 # endif
1435  return Uninitialized<ArrayElementType>::value;
1436  }
1437  } // namespace TensorSubscriptor
1438 } // namespace internal
1439 
1440 
1441 template <int rank_, int dim, typename Number>
1444  Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1445 {
1446  return ::internal::TensorSubscriptor::subscript(
1447  values, i, std::integral_constant<int, dim>());
1448 }
1449 
1450 
1451 template <int rank_, int dim, typename Number>
1452 constexpr DEAL_II_ALWAYS_INLINE
1454  Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1455 {
1456 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1457  AssertIndexRange(i, dim);
1458 # endif
1459 
1460  return values[i];
1461 }
1462 
1463 
1464 template <int rank_, int dim, typename Number>
1465 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1467 {
1468 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1469  Assert(dim != 0,
1470  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1471 # endif
1472 
1473  return TensorAccessors::extract<rank_>(*this, indices);
1474 }
1475 
1476 
1477 
1478 template <int rank_, int dim, typename Number>
1479 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1481 {
1482 # ifndef DEAL_II_COMPILER_CUDA_AWARE
1483  Assert(dim != 0,
1484  ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1485 # endif
1486 
1487  return TensorAccessors::extract<rank_>(*this, indices);
1488 }
1489 
1490 
1491 
1492 template <int rank_, int dim, typename Number>
1493 inline Number *
1495 {
1496  return std::addressof(
1497  this->operator[](this->unrolled_to_component_indices(0)));
1498 }
1499 
1500 
1501 
1502 template <int rank_, int dim, typename Number>
1503 inline const Number *
1505 {
1506  return std::addressof(
1507  this->operator[](this->unrolled_to_component_indices(0)));
1508 }
1509 
1510 
1511 
1512 template <int rank_, int dim, typename Number>
1513 inline Number *
1515 {
1516  return begin_raw() + n_independent_components;
1517 }
1518 
1519 
1520 
1521 template <int rank_, int dim, typename Number>
1522 inline const Number *
1524 {
1525  return begin_raw() + n_independent_components;
1526 }
1527 
1528 
1529 
1530 template <int rank_, int dim, typename Number>
1531 template <typename OtherNumber>
1534 {
1535  // The following loop could be written more concisely using std::copy, but
1536  // that function is only constexpr from C++20 on.
1537  for (unsigned int i = 0; i < dim; ++i)
1538  values[i] = t.values[i];
1539  return *this;
1540 }
1541 
1542 
1543 
1544 template <int rank_, int dim, typename Number>
1547 {
1549  (void)d;
1550 
1551  for (unsigned int i = 0; i < dim; ++i)
1553  return *this;
1554 }
1555 
1556 
1557 # ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1558 template <int rank_, int dim, typename Number>
1561 {
1562  for (unsigned int i = 0; i < dim; ++i)
1563  values[i] = other.values[i];
1564  return *this;
1565 }
1566 
1567 
1568 
1569 template <int rank_, int dim, typename Number>
1572  Tensor<rank_, dim, Number> &&other) noexcept
1573 {
1574  for (unsigned int i = 0; i < dim; ++i)
1575  values[i] = other.values[i];
1576  return *this;
1577 }
1578 # endif
1579 
1580 
1581 template <int rank_, int dim, typename Number>
1582 template <typename OtherNumber>
1583 constexpr inline bool
1585  const Tensor<rank_, dim, OtherNumber> &p) const
1586 {
1587  for (unsigned int i = 0; i < dim; ++i)
1588  if (values[i] != p.values[i])
1589  return false;
1590  return true;
1591 }
1592 
1593 
1594 // At some places in the library, we have Point<0> for formal reasons
1595 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1596 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1597 // in the above function that the loop end check always fails, we
1598 // implement this function here
1599 template <>
1600 template <>
1601 constexpr inline bool
1603 {
1604  return true;
1605 }
1606 
1607 
1608 template <int rank_, int dim, typename Number>
1609 template <typename OtherNumber>
1610 constexpr bool
1612  const Tensor<rank_, dim, OtherNumber> &p) const
1613 {
1614  return !((*this) == p);
1615 }
1616 
1617 
1618 template <int rank_, int dim, typename Number>
1619 template <typename OtherNumber>
1620 constexpr inline DEAL_II_ALWAYS_INLINE
1624 {
1625  for (unsigned int i = 0; i < dim; ++i)
1626  values[i] += p.values[i];
1627  return *this;
1628 }
1629 
1630 
1631 template <int rank_, int dim, typename Number>
1632 template <typename OtherNumber>
1633 constexpr inline DEAL_II_ALWAYS_INLINE
1637 {
1638  for (unsigned int i = 0; i < dim; ++i)
1639  values[i] -= p.values[i];
1640  return *this;
1641 }
1642 
1643 
1644 template <int rank_, int dim, typename Number>
1645 template <typename OtherNumber>
1646 constexpr inline DEAL_II_ALWAYS_INLINE
1648  Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1649 {
1650  for (unsigned int i = 0; i < dim; ++i)
1651  values[i] *= s;
1652  return *this;
1653 }
1654 
1655 
1656 namespace internal
1657 {
1658  namespace TensorImplementation
1659  {
1660  template <int rank,
1661  int dim,
1662  typename Number,
1663  typename OtherNumber,
1664  std::enable_if_t<
1665  !std::is_integral<
1666  typename ProductType<Number, OtherNumber>::type>::value &&
1667  !std::is_same<Number, Differentiation::SD::Expression>::value,
1668  int> = 0>
1669  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1671  const OtherNumber &factor)
1672  {
1673  const Number inverse_factor = Number(1.) / factor;
1674  // recurse over the base objects
1675  for (unsigned int d = 0; d < dim; ++d)
1676  t[d] *= inverse_factor;
1677  }
1678 
1679 
1680  template <int rank,
1681  int dim,
1682  typename Number,
1683  typename OtherNumber,
1684  std::enable_if_t<
1685  std::is_integral<
1686  typename ProductType<Number, OtherNumber>::type>::value ||
1687  std::is_same<Number, Differentiation::SD::Expression>::value,
1688  int> = 0>
1689  constexpr DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
1691  const OtherNumber &factor)
1692  {
1693  // recurse over the base objects
1694  for (unsigned int d = 0; d < dim; ++d)
1695  t[d] /= factor;
1696  }
1697  } // namespace TensorImplementation
1698 } // namespace internal
1699 
1700 
1701 template <int rank_, int dim, typename Number>
1702 template <typename OtherNumber>
1703 constexpr inline DEAL_II_ALWAYS_INLINE
1705  Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1706 {
1708  return *this;
1709 }
1710 
1711 
1712 template <int rank_, int dim, typename Number>
1713 constexpr inline DEAL_II_ALWAYS_INLINE
1716 {
1718 
1719  for (unsigned int i = 0; i < dim; ++i)
1720  tmp.values[i] = -values[i];
1721 
1722  return tmp;
1723 }
1724 
1725 
1726 template <int rank_, int dim, typename Number>
1729 {
1730  // Make things work with AD types
1731  using std::sqrt;
1732  return sqrt(norm_square());
1733 }
1734 
1735 
1736 template <int rank_, int dim, typename Number>
1740 {
1742  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1743  for (unsigned int i = 0; i < dim; ++i)
1744  s += values[i].norm_square();
1745 
1746  return s;
1747 }
1748 
1749 
1750 
1751 template <int rank_, int dim, typename Number>
1752 template <typename OtherNumber>
1753 inline void
1755 {
1756  unroll(result.begin(), result.end());
1757 }
1758 
1759 
1760 
1761 template <int rank_, int dim, typename Number>
1762 template <class Iterator>
1763 inline void
1765  const Iterator end) const
1766 {
1767  AssertDimension(std::distance(begin, end), n_independent_components);
1768  unroll_recursion(begin, end);
1769 }
1770 
1771 
1772 
1773 template <int rank_, int dim, typename Number>
1774 template <typename Iterator>
1775 Iterator
1776 Tensor<rank_, dim, Number>::unroll_recursion(const Iterator current,
1777  const Iterator end) const
1778 {
1779  auto next = current;
1780  for (unsigned int i = 0; i < dim; ++i)
1781  next = values[i].unroll_recursion(next, end);
1782  return next;
1783 }
1784 
1785 
1786 template <int rank_, int dim, typename Number>
1787 constexpr inline unsigned int
1789  const TableIndices<rank_> &indices)
1790 {
1791  unsigned int index = 0;
1792  for (int r = 0; r < rank_; ++r)
1793  index = index * dim + indices[r];
1794 
1795  return index;
1796 }
1797 
1798 
1799 
1800 namespace internal
1801 {
1802  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1803  // and rank=2. Make sure we don't have compiler warnings.
1804 
1805  template <int dim>
1806  inline constexpr unsigned int
1807  mod(const unsigned int x)
1808  {
1809  return x % dim;
1810  }
1811 
1812  template <>
1813  inline unsigned int
1814  mod<0>(const unsigned int x)
1815  {
1816  Assert(false, ExcInternalError());
1817  return x;
1818  }
1819 
1820  template <int dim>
1821  inline constexpr unsigned int
1822  div(const unsigned int x)
1823  {
1824  return x / dim;
1825  }
1826 
1827  template <>
1828  inline unsigned int
1829  div<0>(const unsigned int x)
1830  {
1831  Assert(false, ExcInternalError());
1832  return x;
1833  }
1834 
1835 } // namespace internal
1836 
1837 
1838 
1839 template <int rank_, int dim, typename Number>
1840 constexpr inline TableIndices<rank_>
1842 {
1843  AssertIndexRange(i, n_independent_components);
1844 
1845  TableIndices<rank_> indices;
1846 
1847  unsigned int remainder = i;
1848  for (int r = rank_ - 1; r >= 0; --r)
1849  {
1850  indices[r] = internal::mod<dim>(remainder);
1851  remainder = internal::div<dim>(remainder);
1852  }
1853  Assert(remainder == 0, ExcInternalError());
1854 
1855  return indices;
1856 }
1857 
1858 
1859 template <int rank_, int dim, typename Number>
1860 constexpr inline void
1862 {
1863  for (unsigned int i = 0; i < dim; ++i)
1865 }
1866 
1867 
1868 template <int rank_, int dim, typename Number>
1869 constexpr std::size_t
1871 {
1872  return sizeof(Tensor<rank_, dim, Number>);
1873 }
1874 
1875 
1876 template <int rank_, int dim, typename Number>
1877 template <class Archive>
1878 inline void
1879 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1880 {
1881  ar &values;
1882 }
1883 
1884 
1885 template <int rank_, int dim, typename Number>
1887 
1888 #endif // DOXYGEN
1889 
1890 /* ----------------- Non-member functions operating on tensors. ------------ */
1891 
1904 template <int rank_, int dim, typename Number>
1905 inline std::ostream &
1906 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1907 {
1908  for (unsigned int i = 0; i < dim; ++i)
1909  {
1910  out << p[i];
1911  if (i != dim - 1)
1912  out << ' ';
1913  }
1914 
1915  return out;
1916 }
1917 
1918 
1925 template <int dim, typename Number>
1926 inline std::ostream &
1927 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1928 {
1929  out << static_cast<const Number &>(p);
1930  return out;
1931 }
1932 
1933 
1950 template <int dim, typename Number, typename Other>
1953  operator*(const Other &object, const Tensor<0, dim, Number> &t)
1954 {
1955  return object * static_cast<const Number &>(t);
1956 }
1957 
1958 
1959 
1970 template <int dim, typename Number, typename Other>
1973  operator*(const Tensor<0, dim, Number> &t, const Other &object)
1974 {
1975  return static_cast<const Number &>(t) * object;
1976 }
1977 
1978 
1990 template <int dim, typename Number, typename OtherNumber>
1994  const Tensor<0, dim, OtherNumber> &src2)
1995 {
1996  return static_cast<const Number &>(src1) *
1997  static_cast<const OtherNumber &>(src2);
1998 }
1999 
2000 
2008 template <int dim, typename Number, typename OtherNumber>
2010  Tensor<0,
2011  dim,
2012  typename ProductType<Number,
2013  typename EnableIfScalar<OtherNumber>::type>::type>
2014  operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2015 {
2016  return static_cast<const Number &>(t) / factor;
2017 }
2018 
2019 
2027 template <int dim, typename Number, typename OtherNumber>
2031  const Tensor<0, dim, OtherNumber> &q)
2032 {
2033  return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2034 }
2035 
2036 
2044 template <int dim, typename Number, typename OtherNumber>
2048  const Tensor<0, dim, OtherNumber> &q)
2049 {
2050  return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2051 }
2052 
2053 
2066 template <int rank, int dim, typename Number, typename OtherNumber>
2068  Tensor<rank,
2069  dim,
2070  typename ProductType<Number,
2071  typename EnableIfScalar<OtherNumber>::type>::type>
2072  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2073 {
2074  // recurse over the base objects
2076  for (unsigned int d = 0; d < dim; ++d)
2077  tt[d] = t[d] * factor;
2078  return tt;
2079 }
2080 
2081 
2094 template <int rank, int dim, typename Number, typename OtherNumber>
2096  Tensor<rank,
2097  dim,
2099  OtherNumber>::type>
2100  operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2101 {
2102  // simply forward to the operator above
2103  return t * factor;
2104 }
2105 
2106 
2107 namespace internal
2108 {
2109  namespace TensorImplementation
2110  {
2111  template <int rank,
2112  int dim,
2113  typename Number,
2114  typename OtherNumber,
2115  std::enable_if_t<
2116  !std::is_integral<
2117  typename ProductType<Number, OtherNumber>::type>::value,
2118  int> = 0>
2122  const OtherNumber & factor)
2123  {
2125  const Number inverse_factor = Number(1.) / factor;
2126  // recurse over the base objects
2127  for (unsigned int d = 0; d < dim; ++d)
2128  tt[d] = t[d] * inverse_factor;
2129  return tt;
2130  }
2131 
2132 
2133  template <int rank,
2134  int dim,
2135  typename Number,
2136  typename OtherNumber,
2137  std::enable_if_t<
2138  std::is_integral<
2139  typename ProductType<Number, OtherNumber>::type>::value,
2140  int> = 0>
2144  const OtherNumber & factor)
2145  {
2147  // recurse over the base objects
2148  for (unsigned int d = 0; d < dim; ++d)
2149  tt[d] = t[d] / factor;
2150  return tt;
2151  }
2152  } // namespace TensorImplementation
2153 } // namespace internal
2154 
2155 
2165 template <int rank, int dim, typename Number, typename OtherNumber>
2167  Tensor<rank,
2168  dim,
2169  typename ProductType<Number,
2170  typename EnableIfScalar<OtherNumber>::type>::type>
2171  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2172 {
2174 }
2175 
2176 
2186 template <int rank, int dim, typename Number, typename OtherNumber>
2191 {
2193 
2194  for (unsigned int i = 0; i < dim; ++i)
2195  tmp[i] += q[i];
2196 
2197  return tmp;
2198 }
2199 
2200 
2210 template <int rank, int dim, typename Number, typename OtherNumber>
2215 {
2217 
2218  for (unsigned int i = 0; i < dim; ++i)
2219  tmp[i] -= q[i];
2220 
2221  return tmp;
2222 }
2223 
2230 template <int dim, typename Number, typename OtherNumber>
2231 inline constexpr DEAL_II_ALWAYS_INLINE
2234  const Tensor<0, dim, OtherNumber> &src2)
2235 {
2237 
2238  tmp *= src2;
2239 
2240  return tmp;
2241 }
2242 
2259 template <int rank, int dim, typename Number, typename OtherNumber>
2260 inline constexpr DEAL_II_ALWAYS_INLINE
2263  const Tensor<rank, dim, OtherNumber> &src2)
2264 {
2266 
2267  for (unsigned int i = 0; i < dim; ++i)
2268  tmp[i] = schur_product(Tensor<rank - 1, dim, Number>(src1[i]),
2270 
2271  return tmp;
2272 }
2273 
2316 template <int rank_1,
2317  int rank_2,
2318  int dim,
2319  typename Number,
2320  typename OtherNumber,
2321  typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2322 constexpr inline DEAL_II_ALWAYS_INLINE
2323  typename Tensor<rank_1 + rank_2 - 2,
2324  dim,
2325  typename ProductType<Number, OtherNumber>::type>::tensor_type
2328 {
2329  typename Tensor<rank_1 + rank_2 - 2,
2330  dim,
2331  typename ProductType<Number, OtherNumber>::type>::tensor_type
2332  result{};
2333 
2334  TensorAccessors::internal::
2335  ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2336  reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2337  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2338 
2339  return result;
2340 }
2341 
2342 
2371 template <int index_1,
2372  int index_2,
2373  int rank_1,
2374  int rank_2,
2375  int dim,
2376  typename Number,
2377  typename OtherNumber>
2378 constexpr inline DEAL_II_ALWAYS_INLINE
2379  typename Tensor<rank_1 + rank_2 - 2,
2380  dim,
2381  typename ProductType<Number, OtherNumber>::type>::tensor_type
2384 {
2385  Assert(0 <= index_1 && index_1 < rank_1,
2386  ExcMessage(
2387  "The specified index_1 must lie within the range [0,rank_1)"));
2388  Assert(0 <= index_2 && index_2 < rank_2,
2389  ExcMessage(
2390  "The specified index_2 must lie within the range [0,rank_2)"));
2391 
2392  using namespace TensorAccessors;
2393  using namespace TensorAccessors::internal;
2394 
2395  // Reorder index_1 to the end of src1:
2397  reord_01 = reordered_index_view<index_1, rank_1>(src1);
2398 
2399  // Reorder index_2 to the end of src2:
2400  const ReorderedIndexView<index_2,
2401  rank_2,
2403  reord_02 = reordered_index_view<index_2, rank_2>(src2);
2404 
2405  typename Tensor<rank_1 + rank_2 - 2,
2406  dim,
2407  typename ProductType<Number, OtherNumber>::type>::tensor_type
2408  result{};
2409  TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2410  return result;
2411 }
2412 
2413 
2444 template <int index_1,
2445  int index_2,
2446  int index_3,
2447  int index_4,
2448  int rank_1,
2449  int rank_2,
2450  int dim,
2451  typename Number,
2452  typename OtherNumber>
2453 constexpr inline
2454  typename Tensor<rank_1 + rank_2 - 4,
2455  dim,
2456  typename ProductType<Number, OtherNumber>::type>::tensor_type
2459 {
2460  Assert(0 <= index_1 && index_1 < rank_1,
2461  ExcMessage(
2462  "The specified index_1 must lie within the range [0,rank_1)"));
2463  Assert(0 <= index_3 && index_3 < rank_1,
2464  ExcMessage(
2465  "The specified index_3 must lie within the range [0,rank_1)"));
2466  Assert(index_1 != index_3,
2467  ExcMessage("index_1 and index_3 must not be the same"));
2468  Assert(0 <= index_2 && index_2 < rank_2,
2469  ExcMessage(
2470  "The specified index_2 must lie within the range [0,rank_2)"));
2471  Assert(0 <= index_4 && index_4 < rank_2,
2472  ExcMessage(
2473  "The specified index_4 must lie within the range [0,rank_2)"));
2474  Assert(index_2 != index_4,
2475  ExcMessage("index_2 and index_4 must not be the same"));
2476 
2477  using namespace TensorAccessors;
2478  using namespace TensorAccessors::internal;
2479 
2480  // Reorder index_1 to the end of src1:
2482  reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2483 
2484  // Reorder index_2 to the end of src2:
2486  reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2487 
2488  // Now, reorder index_3 to the end of src1. We have to make sure to
2489  // preserve the original ordering: index_1 has been removed. If
2490  // index_3 > index_1, we have to use (index_3 - 1) instead:
2492  (index_3 < index_1 ? index_3 : index_3 - 1),
2493  rank_1,
2494  ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2495  reord_3 =
2496  TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2497  index_3 - 1,
2498  rank_1 > (reord_1);
2499 
2500  // Now, reorder index_4 to the end of src2. We have to make sure to
2501  // preserve the original ordering: index_2 has been removed. If
2502  // index_4 > index_2, we have to use (index_4 - 1) instead:
2504  (index_4 < index_2 ? index_4 : index_4 - 1),
2505  rank_2,
2506  ReorderedIndexView<index_2, rank_2, const Tensor<rank_2, dim, OtherNumber>>>
2507  reord_4 =
2508  TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2509  index_4 - 1,
2510  rank_2 > (reord_2);
2511 
2512  typename Tensor<rank_1 + rank_2 - 4,
2513  dim,
2514  typename ProductType<Number, OtherNumber>::type>::tensor_type
2515  result{};
2516  TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2517  return result;
2518 }
2519 
2520 
2533 template <int rank, int dim, typename Number, typename OtherNumber>
2534 constexpr inline DEAL_II_ALWAYS_INLINE
2537  const Tensor<rank, dim, OtherNumber> &right)
2538 {
2539  typename ProductType<Number, OtherNumber>::type result{};
2540  TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2541  return result;
2542 }
2543 
2544 
2562 template <template <int, int, typename> class TensorT1,
2563  template <int, int, typename>
2564  class TensorT2,
2565  template <int, int, typename>
2566  class TensorT3,
2567  int rank_1,
2568  int rank_2,
2569  int dim,
2570  typename T1,
2571  typename T2,
2572  typename T3>
2573 constexpr inline DEAL_II_ALWAYS_INLINE
2575  contract3(const TensorT1<rank_1, dim, T1> & left,
2576  const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2577  const TensorT3<rank_2, dim, T3> & right)
2578 {
2579  using return_type =
2581  return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2582  middle,
2583  right);
2584 }
2585 
2586 
2597 template <int rank_1,
2598  int rank_2,
2599  int dim,
2600  typename Number,
2601  typename OtherNumber>
2602 constexpr inline DEAL_II_ALWAYS_INLINE
2606 {
2607  typename Tensor<rank_1 + rank_2,
2608  dim,
2609  typename ProductType<Number, OtherNumber>::type>::tensor_type
2610  result{};
2611  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2612  return result;
2613 }
2614 
2615 
2632 template <int dim, typename Number>
2634 cross_product_2d(const Tensor<1, dim, Number> &src)
2635 {
2636  Assert(dim == 2, ExcInternalError());
2637 
2638  Tensor<1, dim, Number> result;
2639 
2640  result[0] = src[1];
2641  result[1] = -src[0];
2642 
2643  return result;
2644 }
2645 
2646 
2656 template <int dim, typename Number1, typename Number2>
2657 constexpr inline DEAL_II_ALWAYS_INLINE
2659  cross_product_3d(const Tensor<1, dim, Number1> &src1,
2660  const Tensor<1, dim, Number2> &src2)
2661 {
2662  Assert(dim == 3, ExcInternalError());
2663 
2665 
2666  // avoid compiler warnings
2667  constexpr int s0 = 0 % dim;
2668  constexpr int s1 = 1 % dim;
2669  constexpr int s2 = 2 % dim;
2670 
2671  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2672  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2673  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2674 
2675  return result;
2676 }
2677 
2678 
2690 template <int dim, typename Number>
2691 constexpr inline DEAL_II_ALWAYS_INLINE Number
2693 {
2694  // Compute the determinant using the Laplace expansion of the
2695  // determinant. We expand along the last row.
2696  Number det = internal::NumberType<Number>::value(0.0);
2697 
2698  for (unsigned int k = 0; k < dim; ++k)
2699  {
2700  Tensor<2, dim - 1, Number> minor;
2701  for (unsigned int i = 0; i < dim - 1; ++i)
2702  for (unsigned int j = 0; j < dim - 1; ++j)
2703  minor[i][j] = t[i][j < k ? j : j + 1];
2704 
2705  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2706 
2707  det += t[dim - 1][k] * cofactor;
2708  }
2709 
2710  return ((dim % 2 == 0) ? 1. : -1.) * det;
2711 }
2712 
2718 template <typename Number>
2719 constexpr DEAL_II_ALWAYS_INLINE Number
2721 {
2722  return t[0][0];
2723 }
2724 
2730 template <typename Number>
2731 constexpr DEAL_II_ALWAYS_INLINE Number
2733 {
2734  // hard-coded for efficiency reasons
2735  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2736 }
2737 
2743 template <typename Number>
2744 constexpr DEAL_II_ALWAYS_INLINE Number
2746 {
2747  // hard-coded for efficiency reasons
2748  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2749  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2750  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2751  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2752  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2753  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2754  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2755 }
2756 
2757 
2764 template <int dim, typename Number>
2765 constexpr inline DEAL_II_ALWAYS_INLINE Number
2767 {
2768  Number t = d[0][0];
2769  for (unsigned int i = 1; i < dim; ++i)
2770  t += d[i][i];
2771  return t;
2772 }
2773 
2774 
2783 template <int dim, typename Number>
2784 constexpr inline Tensor<2, dim, Number>
2786 {
2787  Number return_tensor[dim][dim];
2788 
2789  // if desired, take over the
2790  // inversion of a 4x4 tensor
2791  // from the FullMatrix
2792  AssertThrow(false, ExcNotImplemented());
2793 
2794  return Tensor<2, dim, Number>(return_tensor);
2795 }
2796 
2797 
2798 #ifndef DOXYGEN
2799 
2800 template <typename Number>
2802  invert(const Tensor<2, 1, Number> &t)
2803 {
2804  Tensor<2, 1, Number> return_tensor;
2805 
2806  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2807 
2808  return return_tensor;
2809 }
2810 
2811 
2812 template <typename Number>
2814  invert(const Tensor<2, 2, Number> &t)
2815 {
2816  Tensor<2, 2, Number> return_tensor;
2817 
2818  const Number inv_det_t = internal::NumberType<Number>::value(
2819  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2820  return_tensor[0][0] = t[1][1];
2821  return_tensor[0][1] = -t[0][1];
2822  return_tensor[1][0] = -t[1][0];
2823  return_tensor[1][1] = t[0][0];
2824  return_tensor *= inv_det_t;
2825 
2826  return return_tensor;
2827 }
2828 
2829 
2830 template <typename Number>
2832  invert(const Tensor<2, 3, Number> &t)
2833 {
2834  Tensor<2, 3, Number> return_tensor;
2835 
2836  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2837  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2838  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2839  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2840  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2841  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2842  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2843  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2844  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2845  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2846  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2847  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2848  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2849  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2850  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2851  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2852  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2853  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2854  const Number inv_det_t = internal::NumberType<Number>::value(
2855  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2856  t[0][2] * return_tensor[2][0]));
2857  return_tensor *= inv_det_t;
2858 
2859  return return_tensor;
2860 }
2861 
2862 #endif /* DOXYGEN */
2863 
2864 
2870 template <int dim, typename Number>
2873 {
2875  for (unsigned int i = 0; i < dim; ++i)
2876  {
2877  tt[i][i] = t[i][i];
2878  for (unsigned int j = i + 1; j < dim; ++j)
2879  {
2880  tt[i][j] = t[j][i];
2881  tt[j][i] = t[i][j];
2882  };
2883  }
2884  return tt;
2885 }
2886 
2887 
2901 template <int dim, typename Number>
2902 constexpr Tensor<2, dim, Number>
2903 adjugate(const Tensor<2, dim, Number> &t)
2904 {
2905  return determinant(t) * invert(t);
2906 }
2907 
2908 
2922 template <int dim, typename Number>
2923 constexpr Tensor<2, dim, Number>
2924 cofactor(const Tensor<2, dim, Number> &t)
2925 {
2926  return transpose(adjugate(t));
2927 }
2928 
2929 
2993 template <int dim, typename Number>
2996 
2997 
3005 template <int dim, typename Number>
3006 inline Number
3008 {
3010  for (unsigned int j = 0; j < dim; ++j)
3011  {
3013  for (unsigned int i = 0; i < dim; ++i)
3015 
3016  if (sum > max)
3017  max = sum;
3018  }
3019 
3020  return max;
3021 }
3022 
3023 
3031 template <int dim, typename Number>
3032 inline Number
3034 {
3036  for (unsigned int i = 0; i < dim; ++i)
3037  {
3039  for (unsigned int j = 0; j < dim; ++j)
3041 
3042  if (sum > max)
3043  max = sum;
3044  }
3045 
3046  return max;
3047 }
3048 
3052 #ifndef DOXYGEN
3053 
3054 
3055 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3056 
3057 // Specialization of functions for ADOL-C number types when
3058 // the advanced branching feature is used
3059 template <int dim>
3060 inline adouble
3062 {
3064  for (unsigned int j = 0; j < dim; ++j)
3065  {
3067  for (unsigned int i = 0; i < dim; ++i)
3068  sum += std::fabs(t[i][j]);
3069 
3070  condassign(max, (sum > max), sum, max);
3071  }
3072 
3073  return max;
3074 }
3075 
3076 
3077 template <int dim>
3078 inline adouble
3080 {
3082  for (unsigned int i = 0; i < dim; ++i)
3083  {
3085  for (unsigned int j = 0; j < dim; ++j)
3086  sum += std::fabs(t[i][j]);
3087 
3088  condassign(max, (sum > max), sum, max);
3089  }
3090 
3091  return max;
3092 }
3093 
3094 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3095 
3096 
3097 #endif // DOXYGEN
3098 
3100 
3101 #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: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:2171
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, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:2047
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:2014
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:1993
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
Number * end_raw()
constexpr Tensor()
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:2233
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
const Number * begin_raw() const
typename Tensor< rank_ - 1, dim, Number >::tensor_type value_type
Definition: tensor.h:536
friend class Tensor
Definition: tensor.h:888
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:2189
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3033
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1953
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3007
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:2100
static constexpr std::size_t memory_consumption()
Number * begin_raw()
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition: tensor.h:2030
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:2262
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:2213
constexpr Tensor & operator=(const Number &d)
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:859
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:1973
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:2072
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:865
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:102
#define DEAL_II_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
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:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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:2121
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_CUDA_HOST_DEV
Definition: numbers.h:34
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static constexpr 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 void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
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)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Definition: tensor.cc:71
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3033
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3007