Reference documentation for deal.II version GIT 3801df8983 2022-05-22 23: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 
187  DEAL_II_DEPRECATED_EARLY
188  Number *
190 
199  DEAL_II_DEPRECATED_EARLY
200  const Number *
201  begin_raw() const;
202 
211  DEAL_II_DEPRECATED_EARLY
212  Number *
214 
224  DEAL_II_DEPRECATED_EARLY
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_EARLY 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, Number>>(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  typename std::enable_if<
1665  !std::is_integral<
1666  typename ProductType<Number, OtherNumber>::type>::value &&
1667  !std::is_same<Number, Differentiation::SD::Expression>::value,
1668  int>::type = 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  typename std::enable_if<
1685  std::is_integral<
1686  typename ProductType<Number, OtherNumber>::type>::value ||
1687  std::is_same<Number, Differentiation::SD::Expression>::value,
1688  int>::type = 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  return std::sqrt(norm_square());
1731 }
1732 
1733 
1734 template <int rank_, int dim, typename Number>
1738 {
1740  typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1741  for (unsigned int i = 0; i < dim; ++i)
1742  s += values[i].norm_square();
1743 
1744  return s;
1745 }
1746 
1747 
1748 
1749 template <int rank_, int dim, typename Number>
1750 template <typename OtherNumber>
1751 inline void
1753 {
1754  unroll(result.begin(), result.end());
1755 }
1756 
1757 
1758 
1759 template <int rank_, int dim, typename Number>
1760 template <class Iterator>
1761 inline void
1763  const Iterator end) const
1764 {
1765  AssertDimension(std::distance(begin, end), n_independent_components);
1766  unroll_recursion(begin, end);
1767 }
1768 
1769 
1770 
1771 template <int rank_, int dim, typename Number>
1772 template <typename Iterator>
1773 Iterator
1774 Tensor<rank_, dim, Number>::unroll_recursion(const Iterator current,
1775  const Iterator end) const
1776 {
1777  auto next = current;
1778  for (unsigned int i = 0; i < dim; ++i)
1779  next = values[i].unroll_recursion(next, end);
1780  return next;
1781 }
1782 
1783 
1784 template <int rank_, int dim, typename Number>
1785 constexpr inline unsigned int
1787  const TableIndices<rank_> &indices)
1788 {
1789  unsigned int index = 0;
1790  for (int r = 0; r < rank_; ++r)
1791  index = index * dim + indices[r];
1792 
1793  return index;
1794 }
1795 
1796 
1797 
1798 namespace internal
1799 {
1800  // unrolled_to_component_indices is instantiated from DataOut for dim==0
1801  // and rank=2. Make sure we don't have compiler warnings.
1802 
1803  template <int dim>
1804  inline constexpr unsigned int
1805  mod(const unsigned int x)
1806  {
1807  return x % dim;
1808  }
1809 
1810  template <>
1811  inline unsigned int
1812  mod<0>(const unsigned int x)
1813  {
1814  Assert(false, ExcInternalError());
1815  return x;
1816  }
1817 
1818  template <int dim>
1819  inline constexpr unsigned int
1820  div(const unsigned int x)
1821  {
1822  return x / dim;
1823  }
1824 
1825  template <>
1826  inline unsigned int
1827  div<0>(const unsigned int x)
1828  {
1829  Assert(false, ExcInternalError());
1830  return x;
1831  }
1832 
1833 } // namespace internal
1834 
1835 
1836 
1837 template <int rank_, int dim, typename Number>
1838 constexpr inline TableIndices<rank_>
1840 {
1841  AssertIndexRange(i, n_independent_components);
1842 
1843  TableIndices<rank_> indices;
1844 
1845  unsigned int remainder = i;
1846  for (int r = rank_ - 1; r >= 0; --r)
1847  {
1848  indices[r] = internal::mod<dim>(remainder);
1849  remainder = internal::div<dim>(remainder);
1850  }
1851  Assert(remainder == 0, ExcInternalError());
1852 
1853  return indices;
1854 }
1855 
1856 
1857 template <int rank_, int dim, typename Number>
1858 constexpr inline void
1860 {
1861  for (unsigned int i = 0; i < dim; ++i)
1863 }
1864 
1865 
1866 template <int rank_, int dim, typename Number>
1867 constexpr std::size_t
1869 {
1870  return sizeof(Tensor<rank_, dim, Number>);
1871 }
1872 
1873 
1874 template <int rank_, int dim, typename Number>
1875 template <class Archive>
1876 inline void
1877 Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1878 {
1879  ar &values;
1880 }
1881 
1882 
1883 template <int rank_, int dim, typename Number>
1885 
1886 #endif // DOXYGEN
1887 
1888 /* ----------------- Non-member functions operating on tensors. ------------ */
1889 
1894 
1902 template <int rank_, int dim, typename Number>
1903 inline std::ostream &
1904 operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1905 {
1906  for (unsigned int i = 0; i < dim; ++i)
1907  {
1908  out << p[i];
1909  if (i != dim - 1)
1910  out << ' ';
1911  }
1912 
1913  return out;
1914 }
1915 
1916 
1923 template <int dim, typename Number>
1924 inline std::ostream &
1925 operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1926 {
1927  out << static_cast<const Number &>(p);
1928  return out;
1929 }
1930 
1931 
1933 
1937 
1938 
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  typename std::enable_if<
2115  !std::is_integral<
2116  typename ProductType<Number, OtherNumber>::type>::value,
2117  int>::type = 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  typename std::enable_if<
2137  std::is_integral<
2138  typename ProductType<Number, OtherNumber>::type>::value,
2139  int>::type = 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 
2274 
2278 
2279 
2316 template <int rank_1,
2317  int rank_2,
2318  int dim,
2319  typename Number,
2320  typename OtherNumber,
2321  typename = typename std::enable_if<rank_1 >= 1 && rank_2 >= 1>::type>
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,
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,
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,
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,
2610  result{};
2611  TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2612  return result;
2613 }
2614 
2615 
2617 
2621 
2622 
2633 template <int dim, typename Number>
2636 {
2637  Assert(dim == 2, ExcInternalError());
2638 
2639  Tensor<1, dim, Number> result;
2640 
2641  result[0] = src[1];
2642  result[1] = -src[0];
2643 
2644  return result;
2645 }
2646 
2647 
2657 template <int dim, typename Number1, typename Number2>
2658 constexpr inline DEAL_II_ALWAYS_INLINE
2661  const Tensor<1, dim, Number2> &src2)
2662 {
2663  Assert(dim == 3, ExcInternalError());
2664 
2666 
2667  // avoid compiler warnings
2668  constexpr int s0 = 0 % dim;
2669  constexpr int s1 = 1 % dim;
2670  constexpr int s2 = 2 % dim;
2671 
2672  result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2673  result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2674  result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2675 
2676  return result;
2677 }
2678 
2679 
2681 
2685 
2686 
2692 template <int dim, typename Number>
2693 constexpr inline DEAL_II_ALWAYS_INLINE Number
2695 {
2696  // Compute the determinant using the Laplace expansion of the
2697  // determinant. We expand along the last row.
2698  Number det = internal::NumberType<Number>::value(0.0);
2699 
2700  for (unsigned int k = 0; k < dim; ++k)
2701  {
2702  Tensor<2, dim - 1, Number> minor;
2703  for (unsigned int i = 0; i < dim - 1; ++i)
2704  for (unsigned int j = 0; j < dim - 1; ++j)
2705  minor[i][j] = t[i][j < k ? j : j + 1];
2706 
2707  const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2708 
2709  det += t[dim - 1][k] * cofactor;
2710  }
2711 
2712  return ((dim % 2 == 0) ? 1. : -1.) * det;
2713 }
2714 
2720 template <typename Number>
2721 constexpr DEAL_II_ALWAYS_INLINE Number
2723 {
2724  return t[0][0];
2725 }
2726 
2732 template <typename Number>
2733 constexpr DEAL_II_ALWAYS_INLINE Number
2735 {
2736  // hard-coded for efficiency reasons
2737  return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2738 }
2739 
2745 template <typename Number>
2746 constexpr DEAL_II_ALWAYS_INLINE Number
2748 {
2749  // hard-coded for efficiency reasons
2750  const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2751  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2752  const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2753  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2754  const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2755  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2756  return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2757 }
2758 
2759 
2766 template <int dim, typename Number>
2767 constexpr inline DEAL_II_ALWAYS_INLINE Number
2769 {
2770  Number t = d[0][0];
2771  for (unsigned int i = 1; i < dim; ++i)
2772  t += d[i][i];
2773  return t;
2774 }
2775 
2776 
2785 template <int dim, typename Number>
2786 constexpr inline Tensor<2, dim, Number>
2788 {
2789  Number return_tensor[dim][dim];
2790 
2791  // if desired, take over the
2792  // inversion of a 4x4 tensor
2793  // from the FullMatrix
2794  AssertThrow(false, ExcNotImplemented());
2795 
2796  return Tensor<2, dim, Number>(return_tensor);
2797 }
2798 
2799 
2800 #ifndef DOXYGEN
2801 
2802 template <typename Number>
2804  invert(const Tensor<2, 1, Number> &t)
2805 {
2806  Tensor<2, 1, Number> return_tensor;
2807 
2808  return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2809 
2810  return return_tensor;
2811 }
2812 
2813 
2814 template <typename Number>
2816  invert(const Tensor<2, 2, Number> &t)
2817 {
2818  Tensor<2, 2, Number> return_tensor;
2819 
2820  const Number inv_det_t = internal::NumberType<Number>::value(
2821  1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2822  return_tensor[0][0] = t[1][1];
2823  return_tensor[0][1] = -t[0][1];
2824  return_tensor[1][0] = -t[1][0];
2825  return_tensor[1][1] = t[0][0];
2826  return_tensor *= inv_det_t;
2827 
2828  return return_tensor;
2829 }
2830 
2831 
2832 template <typename Number>
2834  invert(const Tensor<2, 3, Number> &t)
2835 {
2836  Tensor<2, 3, Number> return_tensor;
2837 
2838  return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2839  internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2840  return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2841  internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2842  return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2843  internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2844  return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2845  internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2846  return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2847  internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2848  return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2849  internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2850  return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2851  internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2852  return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2853  internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2854  return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2855  internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2856  const Number inv_det_t = internal::NumberType<Number>::value(
2857  1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2858  t[0][2] * return_tensor[2][0]));
2859  return_tensor *= inv_det_t;
2860 
2861  return return_tensor;
2862 }
2863 
2864 #endif /* DOXYGEN */
2865 
2866 
2872 template <int dim, typename Number>
2875 {
2877  for (unsigned int i = 0; i < dim; ++i)
2878  {
2879  tt[i][i] = t[i][i];
2880  for (unsigned int j = i + 1; j < dim; ++j)
2881  {
2882  tt[i][j] = t[j][i];
2883  tt[j][i] = t[i][j];
2884  };
2885  }
2886  return tt;
2887 }
2888 
2889 
2903 template <int dim, typename Number>
2904 constexpr Tensor<2, dim, Number>
2906 {
2907  return determinant(t) * invert(t);
2908 }
2909 
2910 
2924 template <int dim, typename Number>
2925 constexpr Tensor<2, dim, Number>
2927 {
2928  return transpose(adjugate(t));
2929 }
2930 
2931 
2995 template <int dim, typename Number>
2998 
2999 
3007 template <int dim, typename Number>
3008 inline Number
3010 {
3012  for (unsigned int j = 0; j < dim; ++j)
3013  {
3015  for (unsigned int i = 0; i < dim; ++i)
3016  sum += std::fabs(t[i][j]);
3017 
3018  if (sum > max)
3019  max = sum;
3020  }
3021 
3022  return max;
3023 }
3024 
3025 
3033 template <int dim, typename Number>
3034 inline Number
3036 {
3038  for (unsigned int i = 0; i < dim; ++i)
3039  {
3041  for (unsigned int j = 0; j < dim; ++j)
3042  sum += std::fabs(t[i][j]);
3043 
3044  if (sum > max)
3045  max = sum;
3046  }
3047 
3048  return max;
3049 }
3050 
3052 
3053 
3054 #ifndef DOXYGEN
3055 
3056 
3057 # ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3058 
3059 // Specialization of functions for ADOL-C number types when
3060 // the advanced branching feature is used
3061 template <int dim>
3062 inline adouble
3064 {
3066  for (unsigned int j = 0; j < dim; ++j)
3067  {
3069  for (unsigned int i = 0; i < dim; ++i)
3070  sum += std::fabs(t[i][j]);
3071 
3072  condassign(max, (sum > max), sum, max);
3073  }
3074 
3075  return max;
3076 }
3077 
3078 
3079 template <int dim>
3080 inline adouble
3082 {
3084  for (unsigned int i = 0; i < dim; ++i)
3085  {
3087  for (unsigned int j = 0; j < dim; ++j)
3088  sum += std::fabs(t[i][j]);
3089 
3090  condassign(max, (sum > max), sum, max);
3091  }
3092 
3093  return max;
3094 }
3095 
3096 # endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3097 
3098 
3099 #endif // DOXYGEN
3100 
3102 
3103 #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
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 Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2926
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
constexpr Tensor< rank_1+rank_2, dim, typename ProductType< Number, OtherNumber >::type > outer_product(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2604
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2905
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< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2874
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:2046
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 Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2694
constexpr Tensor< rank_1+rank_2 - 4, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type double_contract(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2457
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
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
Definition: tensor.h:2536
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:2232
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:2188
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3035
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition: tensor.h:1952
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3009
constexpr Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:2768
constexpr Tensor< rank_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type contract(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2382
constexpr Tensor & operator/=(const OtherNumber &factor)
static constexpr unsigned int dimension
Definition: tensor.h:517
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2635
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()
OtherNumber::type::tensor_type operator*(const Tensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
Definition: tensor.h:2326
constexpr Number determinant(const Tensor< 2, 1, Number > &t)
Definition: tensor.h:2722
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: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: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:1972
constexpr Tensor(const ArrayLike &initializer, std::index_sequence< Indices... >)
constexpr Number determinant(const Tensor< 2, 2, Number > &t)
Definition: tensor.h:2734
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
constexpr Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2787
void serialize(Archive &ar, const unsigned int version)
constexpr Number determinant(const Tensor< 2, 3, Number > &t)
Definition: tensor.h:2747
constexpr Tensor(const array_type &initializer)
static constexpr unsigned int n_independent_components
Definition: tensor.h:528
constexpr Tensor operator-() const
constexpr ProductType< T1, typename ProductType< T2, T3 >::type >::type contract3(const TensorT1< rank_1, dim, T1 > &left, const TensorT2< rank_1+rank_2, dim, T2 > &middle, const TensorT3< rank_2, dim, T3 > &right)
Definition: tensor.h:2575
Tensor< rank_ - 1, dim, Number > values[(dim !=0) ? dim :1]
Definition: tensor.h:865
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660
numbers::NumberTraits< Number >::real_type norm() const
Definition: vector.h:109
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
iterator end()
iterator begin()
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 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:943
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:927
#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:705
decltype(std::declval< T >() *std::declval< U >()) type
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:589
static real_type abs(const number &x)
Definition: numbers.h:611
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2926
constexpr Tensor< 2, dim, Number > adjugate(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2905
constexpr Tensor< 2, dim, Number > transpose(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2874
constexpr Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2694
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3035
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3009
constexpr Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &)
Definition: tensor.h:2787