Reference documentation for deal.II version 9.4.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\}}\)
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
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
38template <typename ElementType, typename MemorySpace>
39class ArrayView;
40template <int dim, typename Number>
41class Point;
42template <int rank_, int dim, typename Number = double>
43class Tensor;
44template <typename Number>
45class Vector;
46template <typename number>
47class FullMatrix;
48namespace Differentiation
49{
50 namespace SD
51 {
52 class Expression;
53 }
54} // namespace Differentiation
55#endif
56
57
87template <int dim, typename Number>
88class Tensor<0, dim, Number>
89{
90public:
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
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 &
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
348 operator-() const;
349
362 constexpr void
364
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
408private:
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
501template <int rank_, int dim, typename Number>
503{
504public:
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
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
693
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
758 operator-() const;
759
772 constexpr void
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
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
861private:
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
897namespace 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> &
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
959template <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
969template <int dim, typename Number>
970template <typename OtherNumber>
972Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
973 : value(internal::NumberType<Number>::value(initializer))
974{}
975
976
977
978template <int dim, typename Number>
979template <typename OtherNumber>
982 : Tensor{p.value}
983{}
984
985
986# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
987template <int dim, typename Number>
990 : value{other.value}
991{}
992
993
994
995template <int dim, typename Number>
998 : value{std::move(other.value)}
999{}
1000# endif
1001
1002
1003template <int dim, typename Number>
1004inline Number *
1006{
1007 return std::addressof(value);
1008}
1009
1010
1011
1012template <int dim, typename Number>
1013inline const Number *
1015{
1016 return std::addressof(value);
1017}
1018
1019
1020
1021template <int dim, typename Number>
1022inline Number *
1024{
1026}
1027
1028
1029
1030template <int dim, typename Number>
1031const Number *
1033{
1035}
1036
1037
1038
1039template <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
1052template <int dim, typename Number>
1053constexpr 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
1066template <int dim, typename Number>
1067template <typename OtherNumber>
1068constexpr inline DEAL_II_ALWAYS_INLINE
1071{
1073 return *this;
1074}
1075
1076
1077# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1078template <int dim, typename Number>
1079constexpr 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
1089template <int dim, typename Number>
1092{
1093 value = std::move(other.value);
1094 return *this;
1095}
1096# endif
1097
1098
1099
1100template <int dim, typename Number>
1101template <typename OtherNumber>
1102constexpr inline DEAL_II_ALWAYS_INLINE
1104 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
1105{
1107 return *this;
1108}
1109
1110
1111template <int dim, typename Number>
1112template <typename OtherNumber>
1113constexpr 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
1128template <int dim, typename Number>
1129template <typename OtherNumber>
1130constexpr bool
1132{
1133 return !((*this) == p);
1134}
1135
1136
1137template <int dim, typename Number>
1138template <typename OtherNumber>
1139constexpr inline DEAL_II_ALWAYS_INLINE
1142{
1143 value += p.value;
1144 return *this;
1145}
1146
1147
1148template <int dim, typename Number>
1149template <typename OtherNumber>
1150constexpr inline DEAL_II_ALWAYS_INLINE
1153{
1154 value -= p.value;
1155 return *this;
1156}
1157
1158
1159
1160namespace 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
1184template <int dim, typename Number>
1185template <typename OtherNumber>
1186constexpr 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
1196template <int dim, typename Number>
1197template <typename OtherNumber>
1199Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1200{
1201 value /= s;
1202 return *this;
1203}
1204
1205
1206template <int dim, typename Number>
1209{
1210 return -value;
1211}
1212
1213
1214template <int dim, typename Number>
1217{
1218 Assert(dim != 0,
1219 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1221}
1222
1223
1224template <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
1239template <int dim, typename Number>
1240template <typename Iterator>
1241Iterator
1242Tensor<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
1257template <int dim, typename Number>
1258constexpr inline void
1260{
1261 // Some auto-differentiable numbers need explicit
1262 // zero initialization.
1264}
1265
1266
1267
1268template <int dim, typename Number>
1269template <class Iterator>
1270inline void
1271Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1272{
1275}
1276
1277
1278
1279template <int dim, typename Number>
1280template <class Archive>
1281inline void
1282Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1283{
1284 ar &value;
1285}
1286
1287
1288template <int dim, typename Number>
1290
1291
1292/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1293
1294template <int rank_, int dim, typename Number>
1295template <typename ArrayLike, std::size_t... indices>
1297Tensor<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
1307template <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
1317template <int rank_, int dim, typename Number>
1319Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1320 : Tensor(initializer, std::make_index_sequence<dim>{})
1321{}
1322
1323
1324
1325template <int rank_, int dim, typename Number>
1326template <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
1339template <int rank_, int dim, typename Number>
1340template <typename OtherNumber>
1343 const Tensor<rank_, dim, OtherNumber> &initializer)
1344 : Tensor(initializer, std::make_index_sequence<dim>{})
1345{}
1346
1347
1348
1349template <int rank_, int dim, typename Number>
1350template <typename OtherNumber>
1351constexpr 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
1359template <int rank_, int dim, typename Number>
1360template <typename OtherNumber>
1362operator 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
1369template <int rank_, int dim, typename Number>
1370constexpr DEAL_II_ALWAYS_INLINE
1372{
1373 for (unsigned int i = 0; i < dim; ++i)
1374 values[i] = other.values[i];
1375}
1376
1377
1378
1379template <int rank_, int dim, typename Number>
1380constexpr 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
1389namespace 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
1441template <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
1451template <int rank_, int dim, typename Number>
1452constexpr 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
1464template <int rank_, int dim, typename Number>
1465constexpr 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
1478template <int rank_, int dim, typename Number>
1479constexpr 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
1492template <int rank_, int dim, typename Number>
1493inline Number *
1495{
1496 return std::addressof(
1497 this->operator[](this->unrolled_to_component_indices(0)));
1498}
1499
1500
1501
1502template <int rank_, int dim, typename Number>
1503inline const Number *
1505{
1506 return std::addressof(
1507 this->operator[](this->unrolled_to_component_indices(0)));
1508}
1509
1510
1511
1512template <int rank_, int dim, typename Number>
1513inline Number *
1515{
1516 return begin_raw() + n_independent_components;
1517}
1518
1519
1520
1521template <int rank_, int dim, typename Number>
1522inline const Number *
1524{
1525 return begin_raw() + n_independent_components;
1526}
1527
1528
1529
1530template <int rank_, int dim, typename Number>
1531template <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
1544template <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
1558template <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
1569template <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
1581template <int rank_, int dim, typename Number>
1582template <typename OtherNumber>
1583constexpr 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
1599template <>
1600template <>
1601constexpr inline bool
1603{
1604 return true;
1605}
1606
1607
1608template <int rank_, int dim, typename Number>
1609template <typename OtherNumber>
1610constexpr bool
1612 const Tensor<rank_, dim, OtherNumber> &p) const
1613{
1614 return !((*this) == p);
1615}
1616
1617
1618template <int rank_, int dim, typename Number>
1619template <typename OtherNumber>
1620constexpr 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
1631template <int rank_, int dim, typename Number>
1632template <typename OtherNumber>
1633constexpr 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
1644template <int rank_, int dim, typename Number>
1645template <typename OtherNumber>
1646constexpr 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
1656namespace 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
1701template <int rank_, int dim, typename Number>
1702template <typename OtherNumber>
1703constexpr inline DEAL_II_ALWAYS_INLINE
1705 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1706{
1708 return *this;
1709}
1710
1711
1712template <int rank_, int dim, typename Number>
1713constexpr 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
1726template <int rank_, int dim, typename Number>
1729{
1730 return std::sqrt(norm_square());
1731}
1732
1733
1734template <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
1749template <int rank_, int dim, typename Number>
1750template <typename OtherNumber>
1751inline void
1753{
1754 unroll(result.begin(), result.end());
1755}
1756
1757
1758
1759template <int rank_, int dim, typename Number>
1760template <class Iterator>
1761inline 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
1771template <int rank_, int dim, typename Number>
1772template <typename Iterator>
1773Iterator
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
1784template <int rank_, int dim, typename Number>
1785constexpr 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
1798namespace 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
1837template <int rank_, int dim, typename Number>
1838constexpr 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
1857template <int rank_, int dim, typename Number>
1858constexpr inline void
1860{
1861 for (unsigned int i = 0; i < dim; ++i)
1863}
1864
1865
1866template <int rank_, int dim, typename Number>
1867constexpr std::size_t
1869{
1870 return sizeof(Tensor<rank_, dim, Number>);
1871}
1872
1873
1874template <int rank_, int dim, typename Number>
1875template <class Archive>
1876inline void
1877Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1878{
1879 ar &values;
1880}
1881
1882
1883template <int rank_, int dim, typename Number>
1885
1886#endif // DOXYGEN
1887
1888/* ----------------- Non-member functions operating on tensors. ------------ */
1889
1894
1902template <int rank_, int dim, typename Number>
1903inline std::ostream &
1904operator<<(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
1923template <int dim, typename Number>
1924inline std::ostream &
1925operator<<(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
1949template <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
1969template <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
1989template <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
2007template <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
2026template <int dim, typename Number, typename OtherNumber>
2031{
2032 return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2033}
2034
2035
2043template <int dim, typename Number, typename OtherNumber>
2048{
2049 return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2050}
2051
2052
2065template <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
2093template <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
2106namespace 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<
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<
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
2164template <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
2185template <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
2209template <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
2229template <int dim, typename Number, typename OtherNumber>
2230inline constexpr DEAL_II_ALWAYS_INLINE
2233 const Tensor<0, dim, OtherNumber> &src2)
2234{
2236
2237 tmp *= src2;
2238
2239 return tmp;
2240}
2241
2258template <int rank, int dim, typename Number, typename OtherNumber>
2259inline constexpr DEAL_II_ALWAYS_INLINE
2263{
2265
2266 for (unsigned int i = 0; i < dim; ++i)
2269
2270 return tmp;
2271}
2272
2274
2278
2279
2316template <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>
2322constexpr inline DEAL_II_ALWAYS_INLINE
2323 typename Tensor<rank_1 + rank_2 - 2,
2324 dim,
2325 typename ProductType<Number, OtherNumber>::type>::tensor_type
2328{
2329 typename Tensor<rank_1 + rank_2 - 2,
2330 dim,
2331 typename ProductType<Number, OtherNumber>::type>::tensor_type
2332 result{};
2333
2334 TensorAccessors::internal::
2335 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2336 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2337 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2338
2339 return result;
2340}
2341
2342
2371template <int index_1,
2372 int index_2,
2373 int rank_1,
2374 int rank_2,
2375 int dim,
2376 typename Number,
2377 typename OtherNumber>
2378constexpr inline DEAL_II_ALWAYS_INLINE
2379 typename Tensor<rank_1 + rank_2 - 2,
2380 dim,
2381 typename ProductType<Number, OtherNumber>::type>::tensor_type
2384{
2385 Assert(0 <= index_1 && index_1 < rank_1,
2386 ExcMessage(
2387 "The specified index_1 must lie within the range [0,rank_1)"));
2388 Assert(0 <= index_2 && index_2 < rank_2,
2389 ExcMessage(
2390 "The specified index_2 must lie within the range [0,rank_2)"));
2391
2392 using namespace TensorAccessors;
2393 using namespace TensorAccessors::internal;
2394
2395 // Reorder index_1 to the end of src1:
2397 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2398
2399 // Reorder index_2 to the end of src2:
2400 const ReorderedIndexView<index_2,
2401 rank_2,
2403 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2404
2405 typename Tensor<rank_1 + rank_2 - 2,
2406 dim,
2407 typename ProductType<Number, OtherNumber>::type>::tensor_type
2408 result{};
2409 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2410 return result;
2411}
2412
2413
2444template <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>
2453constexpr 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,
2507 reord_4 =
2508 TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2509 index_4 - 1,
2510 rank_2 > (reord_2);
2511
2512 typename Tensor<rank_1 + rank_2 - 4,
2513 dim,
2514 typename ProductType<Number, OtherNumber>::type>::tensor_type
2515 result{};
2516 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2517 return result;
2518}
2519
2520
2533template <int rank, int dim, typename Number, typename OtherNumber>
2534constexpr inline DEAL_II_ALWAYS_INLINE
2537 const Tensor<rank, dim, OtherNumber> &right)
2538{
2540 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2541 return result;
2542}
2543
2544
2562template <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>
2573constexpr 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
2597template <int rank_1,
2598 int rank_2,
2599 int dim,
2600 typename Number,
2601 typename OtherNumber>
2602constexpr inline DEAL_II_ALWAYS_INLINE
2606{
2607 typename Tensor<rank_1 + rank_2,
2608 dim,
2609 typename ProductType<Number, OtherNumber>::type>::tensor_type
2610 result{};
2611 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2612 return result;
2613}
2614
2615
2617
2621
2622
2633template <int dim, typename Number>
2636{
2637 Assert(dim == 2, ExcInternalError());
2638
2640
2641 result[0] = src[1];
2642 result[1] = -src[0];
2643
2644 return result;
2645}
2646
2647
2657template <int dim, typename Number1, typename Number2>
2658constexpr 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
2692template <int dim, typename Number>
2693constexpr 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
2720template <typename Number>
2721constexpr DEAL_II_ALWAYS_INLINE Number
2723{
2724 return t[0][0];
2725}
2726
2732template <typename Number>
2733constexpr 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
2745template <typename Number>
2746constexpr 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
2766template <int dim, typename Number>
2767constexpr 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
2785template <int dim, typename Number>
2786constexpr 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
2795
2796 return Tensor<2, dim, Number>(return_tensor);
2797}
2798
2799
2800#ifndef DOXYGEN
2801
2802template <typename Number>
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
2814template <typename Number>
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
2832template <typename Number>
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
2872template <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
2903template <int dim, typename Number>
2904constexpr Tensor<2, dim, Number>
2906{
2907 return determinant(t) * invert(t);
2908}
2909
2910
2924template <int dim, typename Number>
2925constexpr Tensor<2, dim, Number>
2927{
2928 return transpose(adjugate(t));
2929}
2930
2931
2995template <int dim, typename Number>
2998
2999
3007template <int dim, typename Number>
3008inline 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
3033template <int dim, typename Number>
3034inline 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
3061template <int dim>
3062inline 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
3079template <int dim>
3080inline 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
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: point.h:111
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
constexpr Tensor & operator*=(const OtherNumber &factor)
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 Tensor & operator=(const OtherNumber &d)
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:122
constexpr Tensor operator-() const
Definition: tensor.h:503
constexpr Tensor & operator/=(const OtherNumber &factor)
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(const Tensor< 1, dim, Tensor< rank_ - 1, dim, OtherNumber > > &initializer)
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 const Number & operator[](const TableIndices< rank_ > &indices) const
constexpr ProductType< T1, typenameProductType< 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
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
Number * begin_raw()
static constexpr unsigned int rank
Definition: tensor.h:522
constexpr Number determinant(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2694
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:1992
const Number * begin_raw() const
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
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 value_type & operator[](const unsigned int i)
Number * end_raw()
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
constexpr Tensor< rank_1+rank_2-4, dim, typenameProductType< 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
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
static constexpr unsigned int dimension
Definition: tensor.h:517
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
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()
constexpr Tensor & operator=(const Number &d)
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 & operator+=(const Tensor< rank_, dim, OtherNumber > &)
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:2261
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:2212
Tensor< rank_, dim, Number > tensor_type
Definition: tensor.h:859
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
constexpr Tensor< rank_1+rank_2-2, dim, typenameProductType< 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< 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 Tensor & operator*=(const OtherNumber &factor)
constexpr const value_type & operator[](const unsigned int i) const
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
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
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
Definition: vector.h:109
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
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 * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
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
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
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