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