deal.II version GIT relicensing-2069-gd7d5dbb18d 2024-11-01 22:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tensor_h
16#define dealii_tensor_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/kokkos.h>
26
27#include <Kokkos_Array.hpp>
28
29#ifdef DEAL_II_WITH_ADOLC
30# include <adolc/adouble.h> // Taped double
31#endif
32
33#include <cmath>
34#include <ostream>
35#include <type_traits>
36
38
39// Forward declarations:
40#ifndef DOXYGEN
41template <typename ElementType, typename MemorySpace>
42class ArrayView;
43
44template <int dim, typename Number>
46class Point;
47
48template <int rank_, int dim, typename Number = double>
49class Tensor;
50template <typename Number>
51class Vector;
52template <typename number>
53class FullMatrix;
54namespace Differentiation
55{
56 namespace SD
57 {
58 class Expression;
59 }
60} // namespace Differentiation
61#endif
62
63
93template <int dim, typename Number>
94class Tensor<0, dim, Number>
95{
96public:
97 static_assert(dim >= 0,
98 "Tensors must have a dimension greater than or equal to one.");
99
108 static constexpr unsigned int dimension = dim;
109
113 static constexpr unsigned int rank = 0;
114
118 static constexpr unsigned int n_independent_components = 1;
119
129
134 using value_type = Number;
135
141 using array_type = Number;
142
148 constexpr DEAL_II_HOST_DEVICE
150
158 template <typename OtherNumber>
159 constexpr DEAL_II_HOST_DEVICE
161
167 template <typename OtherNumber>
168 constexpr DEAL_II_HOST_DEVICE
169 Tensor(const OtherNumber &initializer);
170
171#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
175 constexpr DEAL_II_HOST_DEVICE
176 Tensor(const Tensor<0, dim, Number> &other);
177
181 constexpr DEAL_II_HOST_DEVICE
182 Tensor(Tensor<0, dim, Number> &&other) noexcept;
183#endif
184
194 constexpr DEAL_II_HOST_DEVICE
195 operator Number &();
196
205 constexpr DEAL_II_HOST_DEVICE operator const Number &() const;
206
214 template <typename OtherNumber>
215 constexpr DEAL_II_HOST_DEVICE Tensor &
217
218#if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
227 constexpr DEAL_II_HOST_DEVICE Tensor &
229#endif
230
231#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
236 operator=(Tensor<0, dim, Number> &&other) noexcept;
237#endif
238
245 template <typename OtherNumber>
246 constexpr DEAL_II_HOST_DEVICE Tensor &
247 operator=(const OtherNumber &d) &;
248
254 template <typename OtherNumber>
255 constexpr DEAL_II_HOST_DEVICE Tensor &
256 operator=(const OtherNumber &d) && = delete;
257
261 template <typename OtherNumber>
262 constexpr bool
264
268 template <typename OtherNumber>
269 constexpr bool
271
277 template <typename OtherNumber>
278 constexpr DEAL_II_HOST_DEVICE Tensor &
280
286 template <typename OtherNumber>
287 constexpr DEAL_II_HOST_DEVICE Tensor &
289
295 template <typename OtherNumber>
296 constexpr DEAL_II_HOST_DEVICE Tensor &
297 operator*=(const OtherNumber &factor);
298
304 template <typename OtherNumber>
305 constexpr DEAL_II_HOST_DEVICE Tensor &
306 operator/=(const OtherNumber &factor);
307
314 operator-() const;
315
328 constexpr void
330
337 norm() const;
338
346 norm_square() const;
347
355 template <class Iterator>
356 void
357 unroll(const Iterator begin, const Iterator end) const;
358
364 template <class Archive>
365 void
366 serialize(Archive &ar, const unsigned int version);
367
372 using tensor_type = Number;
373
374private:
378 Number value;
379
380 // Allow an arbitrary Tensor to access the underlying values.
381 template <int, int, typename>
382 friend class Tensor;
383};
384
385
386
460template <int rank_, int dim, typename Number>
462{
463public:
464 static_assert(rank_ >= 1,
465 "Tensors must have a rank greater than or equal to one.");
466 static_assert(dim >= 0,
467 "Tensors must have a dimension greater than or equal to zero.");
476 static constexpr unsigned int dimension = dim;
477
481 static constexpr unsigned int rank = rank_;
482
497 static constexpr unsigned int n_independent_components =
498 Tensor<rank_ - 1, dim>::n_independent_components * dim;
499
506 std::conditional_t<rank_ == 1, Number, Tensor<rank_ - 1, dim, Number>>;
507
519 using array_type = std::conditional_t<
520 rank_ == 1,
521 Number[(dim != 0) ? dim : 1],
522 typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1]>;
523
531
537 constexpr DEAL_II_HOST_DEVICE explicit Tensor(const array_type &initializer);
538
551 template <typename ElementType, typename MemorySpace>
552 constexpr DEAL_II_HOST_DEVICE explicit Tensor(
553 const ArrayView<ElementType, MemorySpace> &initializer);
554
562 template <typename OtherNumber>
563 constexpr DEAL_II_HOST_DEVICE
565
569 template <typename OtherNumber>
570 constexpr Tensor(
571 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
572
576 template <typename OtherNumber>
577 constexpr
578 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
579
580#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
584 constexpr Tensor(const Tensor<rank_, dim, Number> &);
585
589 constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
590#endif
591
598 operator[](const unsigned int i);
599
605 constexpr DEAL_II_HOST_DEVICE const value_type &
606 operator[](const unsigned int i) const;
607
611 constexpr const Number &
612 operator[](const TableIndices<rank_> &indices) const;
613
617 constexpr Number &
619
624 Number *
626
631 const Number *
632 begin_raw() const;
633
638 Number *
640
645 const Number *
646 end_raw() const;
647
655 template <typename OtherNumber>
656 constexpr DEAL_II_HOST_DEVICE Tensor &
658
665 constexpr DEAL_II_HOST_DEVICE Tensor &
666 operator=(const Number &d) &;
667
673 constexpr DEAL_II_HOST_DEVICE Tensor &
674 operator=(const Number &d) && = delete;
675
676#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
682
688#endif
689
693 template <typename OtherNumber>
694 constexpr bool
696
700 template <typename OtherNumber>
701 constexpr bool
703
709 template <typename OtherNumber>
710 constexpr DEAL_II_HOST_DEVICE Tensor &
712
718 template <typename OtherNumber>
719 constexpr DEAL_II_HOST_DEVICE Tensor &
721
728 template <typename OtherNumber>
729 constexpr DEAL_II_HOST_DEVICE Tensor &
730 operator*=(const OtherNumber &factor);
731
737 template <typename OtherNumber>
738 constexpr DEAL_II_HOST_DEVICE Tensor &
739 operator/=(const OtherNumber &factor);
740
747 operator-() const;
748
761 constexpr void
763
773 norm() const;
774
781 constexpr DEAL_II_HOST_DEVICE
783 norm_square() const;
784
795 template <class Iterator>
796 void
797 unroll(const Iterator begin, const Iterator end) const;
798
803 static constexpr DEAL_II_HOST_DEVICE unsigned int
805
812 unrolled_to_component_indices(const unsigned int i);
813
818 static constexpr std::size_t
820
826 template <class Archive>
827 void
828 serialize(Archive &ar, const unsigned int version);
829
835
836private:
842#if KOKKOS_VERSION >= 30700
843 std::conditional_t<rank_ == 1,
844 Kokkos::Array<Number, dim>,
845 Kokkos::Array<Tensor<rank_ - 1, dim, Number>, dim>>
846#else
847 std::conditional_t<rank_ == 1,
848 std::array<Number, dim>,
849 std::array<Tensor<rank_ - 1, dim, Number>, dim>>
850#endif
852
859 template <typename ArrayLike, std::size_t... Indices>
860 constexpr DEAL_II_HOST_DEVICE
861 Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
862
863 // Allow an arbitrary Tensor to access the underlying values.
864 template <int, int, typename>
865 friend class Tensor;
866
867 // Point is allowed access to the coordinates. This is supposed to improve
868 // speed.
869 friend class Point<dim, Number>;
870};
871
872
873#ifndef DOXYGEN
874namespace internal
875{
876 // Workaround: The following 4 overloads are necessary to be able to
877 // compile the library with Apple Clang 8 and older. We should remove
878 // these overloads again when we bump the minimal required version to
879 // something later than clang-3.6 / Apple Clang 6.3.
880 template <int rank, int dim, typename T, typename U>
881 struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
882 {
883 using type =
885 };
886
887 template <int rank, int dim, typename T, typename U>
888 struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
889 {
890 using type =
892 };
893
894 template <typename T, int rank, int dim, typename U>
895 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
896 {
897 using type =
899 };
900
901 template <int rank, int dim, typename T, typename U>
902 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
903 {
904 using type =
906 };
907 // end workaround
908
913 template <int rank, int dim, typename T>
914 struct NumberType<Tensor<rank, dim, T>>
915 {
916 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const
919 {
920 return t;
921 }
922
924 value(const T &t)
925 {
927 tmp = t;
928 return tmp;
929 }
930 };
931} // namespace internal
932
933
934/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
935
936
937template <int dim, typename Number>
940 // Some auto-differentiable numbers need explicit
941 // zero initialization such as adtl::adouble.
942 : Tensor{0.0}
943{}
944
945
946
947template <int dim, typename Number>
948template <typename OtherNumber>
950Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
951 : value(internal::NumberType<Number>::value(initializer))
952{}
953
954
955
956template <int dim, typename Number>
957template <typename OtherNumber>
960 : Tensor{p.value}
961{}
962
963
964# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
965template <int dim, typename Number>
968 : value{other.value}
969{}
970
971
972
973template <int dim, typename Number>
976 : value{std::move(other.value)}
977{}
978# endif
979
980
981
982template <int dim, typename Number>
985{
986 Assert(dim != 0,
987 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
988 return value;
989}
990
991
992template <int dim, typename Number>
993constexpr inline DEAL_II_ALWAYS_INLINE
995{
996 Assert(dim != 0,
997 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
998 return value;
999}
1000
1001
1002
1003template <int dim, typename Number>
1004template <typename OtherNumber>
1007{
1009 return *this;
1010}
1011
1012
1013# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1014template <int dim, typename Number>
1017{
1018 value = p.value;
1019 return *this;
1020}
1021# endif
1022
1023# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1024template <int dim, typename Number>
1027{
1028 value = std::move(other.value);
1029 return *this;
1030}
1031# endif
1032
1033
1034
1035template <int dim, typename Number>
1036template <typename OtherNumber>
1038Tensor<0, dim, Number>::operator=(const OtherNumber &d) &
1039{
1041 return *this;
1042}
1043
1044
1045template <int dim, typename Number>
1046template <typename OtherNumber>
1047constexpr inline bool
1049{
1050# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1051 Assert(!(std::is_same_v<Number, adouble> ||
1052 std::is_same_v<OtherNumber, adouble>),
1053 ExcMessage(
1054 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1055 "been extended to support advanced branching."));
1056# endif
1057
1058 return numbers::values_are_equal(value, p.value);
1059}
1060
1061
1062template <int dim, typename Number>
1063template <typename OtherNumber>
1064constexpr bool
1066{
1067 return !((*this) == p);
1068}
1069
1070
1071template <int dim, typename Number>
1072template <typename OtherNumber>
1075{
1076 value += p.value;
1077 return *this;
1078}
1079
1080
1081template <int dim, typename Number>
1082template <typename OtherNumber>
1085{
1086 value -= p.value;
1087 return *this;
1088}
1089
1090
1091
1092namespace internal
1093{
1094 namespace ComplexWorkaround
1095 {
1096 template <typename Number, typename OtherNumber>
1098 multiply_assign_scalar(Number &val, const OtherNumber &s)
1099 {
1100 val *= s;
1101 }
1102
1103 template <typename Number, typename OtherNumber>
1105 multiply_assign_scalar(std::complex<Number> &val, const OtherNumber &s)
1106 {
1107# if KOKKOS_VERSION >= 30600
1108 KOKKOS_IF_ON_HOST((val *= s;))
1109 KOKKOS_IF_ON_DEVICE(({
1110 (void)val;
1111 (void)s;
1112 Kokkos::abort(
1113 "This function is not implemented for std::complex<Number>!\n");
1114 }))
1115# else
1116# ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1117 val *= s;
1118# else
1119 (void)val;
1120 (void)s;
1121 Kokkos::abort(
1122 "This function is not implemented for std::complex<Number>!\n");
1123# endif
1124# endif
1125 }
1126 } // namespace ComplexWorkaround
1127} // namespace internal
1128
1129
1130template <int dim, typename Number>
1131template <typename OtherNumber>
1133Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1134{
1135 internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1136 return *this;
1137}
1138
1139
1140
1141template <int dim, typename Number>
1142template <typename OtherNumber>
1144Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1145{
1146 value /= s;
1147 return *this;
1148}
1149
1150
1151template <int dim, typename Number>
1154{
1155 return -value;
1156}
1157
1158
1159template <int dim, typename Number>
1162{
1163 Assert(dim != 0,
1164 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1166}
1167
1168
1169template <int dim, typename Number>
1173{
1174 Assert(dim != 0,
1175 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1177}
1178
1179
1180
1181template <int dim, typename Number>
1182constexpr inline void
1184{
1185 // Some auto-differentiable numbers need explicit
1186 // zero initialization.
1188}
1189
1190
1191
1192template <int dim, typename Number>
1193template <class Iterator>
1194inline void
1195Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1196{
1197 (void)end;
1198 AssertDimension(std::distance(begin, end), n_independent_components);
1199 Assert(dim != 0,
1200 ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1201 Assert(std::distance(begin, end) >= 1,
1202 ExcMessage("The provided iterator range must contain at least one "
1203 "element."));
1204 *begin = value;
1205}
1206
1207
1208
1209template <int dim, typename Number>
1210template <class Archive>
1211inline void
1212Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1213{
1214 ar &value;
1215}
1216
1217
1218template <int dim, typename Number>
1220
1221
1222/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1223
1224template <int rank_, int dim, typename Number>
1225template <typename ArrayLike, std::size_t... indices>
1227Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1228 std::index_sequence<indices...>)
1229 // Extract from the 'initializer' a sequence of elements via template
1230 // pack evaluation. This could be as easy as
1231 // values{{ (initializer[indices])... }}
1232 // but of course in practice it is not. The challenge is that if rank>1,
1233 // we want to pass the elements initializer[indices] down to the next
1234 // lower rank tensor for evaluation unchanged. But at the rank==1 level,
1235 // we need to convert to the scalar type 'Number'. This would all be
1236 // relatively straightforward if we could rely on automatic type
1237 // conversion, but for some autodifferentiation types, the conversion
1238 // from the AD to double (i.e., the extraction of a scalar value) is
1239 // not implicit, and we need to call internal::NumberType<Number>::value() --
1240 // but as mentioned, we can only do that for rank==1.
1241 //
1242 // We can achieve all of this by dispatching to a lambda function within
1243 // which we can use a 'if constexpr'.
1244 : values{{([&initializer]() -> value_type {
1245 if constexpr (rank_ == 1)
1246 return internal::NumberType<Number>::value(initializer[indices]);
1247 else
1248 return value_type(initializer[indices]);
1249 }())...}}
1250{
1251 static_assert(sizeof...(indices) == dim,
1252 "dim should match the number of indices");
1253}
1254
1255
1256# ifdef DEAL_II_HAVE_CXX20
1257
1258template <int rank_, int dim, typename Number>
1261 : values(
1262 // In order to initialize the Kokkos::Array<Number,dim>, we would need a
1263 // brace-enclosed list of length 'dim'. There is no way in C++ to create
1264 // such a list in-place, but we can come up with a lambda function that
1265 // expands such a list via template-pack expansion, and then uses this
1266 // list to initialize a Kokkos::Array which it then returns.
1267 //
1268 // The trick to come up with such a lambda function is to have a function
1269 // that takes an argument that depends on a template-pack of integers.
1270 // We will call the function with an integer list of length 'dim', and
1271 // in the function itself expand that pack in a way that it serves as
1272 // a brace-enclosed list of initializers for a Kokkos::Array.
1273 //
1274 // Of course, we do not want to initialize the array with the integers,
1275 // but with zeros. (Or, more correctly, a zero of the element type.)
1276 // The canonical way to do this would be using the comma operator:
1277 // (sequence_element, 0.0)
1278 // returns zero, and
1279 // (sequence, 0.0)...
1280 // returns a list of zeros of the right length. Unfortunately, some
1281 // compilers then warn that the left side of the comma expression has
1282 // no effect -- well, bummer, that was of course exactly the idea.
1283 // We could work around this by using
1284 // (sequence_element * 0.0)
1285 // instead, assuming that the compiler will optimize (known) integer
1286 // times zero to zero, and similarly for (known) integer times times
1287 // default-initialized tensor.
1288 //
1289 // But, instead of relying on compiler optimizations, a better way is
1290 // to simply have another (nested) lambda function that takes the
1291 // integer sequence element as an argument and ignores it, just
1292 // returning a zero instead.
1293 []<std::size_t... I>(
1294 const std::index_sequence<I...> &) constexpr -> decltype(values) {
1295 if constexpr (dim == 0)
1296 {
1297 return {};
1298 }
1299 else if constexpr (rank_ == 1)
1300 {
1301 auto get_zero_and_ignore_argument = [](int) {
1303 };
1304 return {{(get_zero_and_ignore_argument(I))...}};
1305 }
1306 else
1307 {
1308 auto get_zero_and_ignore_argument = [](int) {
1309 return Tensor<rank_ - 1, dim, Number>();
1310 };
1311 return {{(get_zero_and_ignore_argument(I))...}};
1312 }
1313 }(std::make_index_sequence<dim>()))
1314{}
1315
1316# else
1317
1318// The C++17 case works in essence the same, except that we can't use a
1319// lambda function with explicit template parameters, i.e., we can't do
1320// []<std::size_t... I>(const std::index_sequence<I...> &)
1321// as above because that's a C++20 feature. Lambda functions in C++17 can
1322// have template packs as arguments, but we need the ability to *name*
1323// that template pack (the 'I' above) and that's not possible in C++17.
1324//
1325// We work around this by moving the lambda function to a global function
1326// and using the traditional template syntax on it.
1327namespace internal
1328{
1329 namespace TensorInitialization
1330 {
1331 template <int rank, int dim, typename Number, std::size_t... I>
1332# if KOKKOS_VERSION >= 30700
1333 constexpr Kokkos::Array<typename Tensor<rank, dim, Number>::value_type, dim>
1334# else
1335 constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
1336# endif
1337 make_zero_array(const std::index_sequence<I...> &)
1338 {
1339 static_assert(sizeof...(I) == dim, "This is bad.");
1340
1341 // First peel off the case dim==0. If we don't, some compilers
1342 // will warn below that we define these lambda functions but
1343 // never use them (because the expanded list has zero elements,
1344 // and the get_zero_and_ignore_argument() function is not used...)
1345 if constexpr (dim == 0)
1346 {
1347 return {};
1348 }
1349 else if constexpr (rank == 1)
1350 {
1351 auto get_zero_and_ignore_argument = [](int) {
1353 };
1354 return {{(get_zero_and_ignore_argument(I))...}};
1355 }
1356 else
1357 {
1358 auto get_zero_and_ignore_argument = [](int) {
1359 return Tensor<rank - 1, dim, Number>();
1360 };
1361 return {{(get_zero_and_ignore_argument(I))...}};
1362 }
1363 }
1364 } // namespace TensorInitialization
1365} // namespace internal
1366
1367
1368template <int rank_, int dim, typename Number>
1371 : values(internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
1372 std::make_index_sequence<dim>()))
1373{}
1374
1375
1376# endif
1377
1378
1379template <int rank_, int dim, typename Number>
1381Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1382 : Tensor(initializer, std::make_index_sequence<dim>{})
1383{}
1384
1385
1386
1387template <int rank_, int dim, typename Number>
1388template <typename ElementType, typename MemorySpace>
1391 const ArrayView<ElementType, MemorySpace> &initializer)
1392{
1393 // make nvcc happy
1394 const int my_n_independent_components = n_independent_components;
1395 AssertDimension(initializer.size(), my_n_independent_components);
1396
1397 for (unsigned int i = 0; i < my_n_independent_components; ++i)
1398 (*this)[unrolled_to_component_indices(i)] = initializer[i];
1399}
1400
1401
1402
1403template <int rank_, int dim, typename Number>
1404template <typename OtherNumber>
1407 const Tensor<rank_, dim, OtherNumber> &initializer)
1408 : Tensor(initializer, std::make_index_sequence<dim>{})
1409{}
1410
1411
1412
1413template <int rank_, int dim, typename Number>
1414template <typename OtherNumber>
1415constexpr DEAL_II_ALWAYS_INLINE
1417 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1418 : Tensor(initializer, std::make_index_sequence<dim>{})
1419{}
1420
1421
1422
1423template <int rank_, int dim, typename Number>
1424template <typename OtherNumber>
1426operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1427{
1428 Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> x;
1429 std::copy(values.begin(), values.end(), x.values.begin());
1430 return x;
1431}
1432
1433
1434# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1435template <int rank_, int dim, typename Number>
1436constexpr DEAL_II_ALWAYS_INLINE
1438 : values(other.values)
1439{}
1440
1441
1442
1443template <int rank_, int dim, typename Number>
1444constexpr DEAL_II_ALWAYS_INLINE
1446 : values(std::move(other.values))
1447{}
1448# endif
1449
1450
1451
1452template <int rank_, int dim, typename Number>
1455 Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1456{
1457 Assert(dim != 0,
1458 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1459 AssertIndexRange(i, dim);
1460 DEAL_II_CXX23_ASSUME(i < dim);
1461
1462 return values[i];
1463}
1464
1465
1466template <int rank_, int dim, typename Number>
1467constexpr DEAL_II_ALWAYS_INLINE
1469 Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1470{
1471 Assert(dim != 0,
1472 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1473 AssertIndexRange(i, dim);
1474 DEAL_II_CXX23_ASSUME(i < dim);
1475
1476 return values[i];
1477}
1478
1479
1480template <int rank_, int dim, typename Number>
1481constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1483{
1484 Assert(dim != 0,
1485 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1486
1487 return TensorAccessors::extract<rank_>(*this, indices);
1488}
1489
1490
1491
1492template <int rank_, int dim, typename Number>
1493constexpr inline DEAL_II_ALWAYS_INLINE Number &
1495{
1496 Assert(dim != 0,
1497 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1498
1499 return TensorAccessors::extract<rank_>(*this, indices);
1500}
1501
1502
1503
1504template <int rank_, int dim, typename Number>
1505inline Number *
1507{
1508 static_assert(rank_ == 1,
1509 "This function is only available for rank-1 tensors "
1510 "because higher-rank tensors may not store their elements "
1511 "in a contiguous array.");
1512
1513 return std::addressof(
1514 this->operator[](this->unrolled_to_component_indices(0)));
1515}
1516
1517
1518
1519template <int rank_, int dim, typename Number>
1520inline const Number *
1522{
1523 static_assert(rank_ == 1,
1524 "This function is only available for rank-1 tensors "
1525 "because higher-rank tensors may not store their elements "
1526 "in a contiguous array.");
1527
1528 return std::addressof(
1529 this->operator[](this->unrolled_to_component_indices(0)));
1530}
1531
1532
1533
1534template <int rank_, int dim, typename Number>
1535inline Number *
1537{
1538 static_assert(rank_ == 1,
1539 "This function is only available for rank-1 tensors "
1540 "because higher-rank tensors may not store their elements "
1541 "in a contiguous array.");
1542
1543 return begin_raw() + n_independent_components;
1544}
1545
1546
1547
1548template <int rank_, int dim, typename Number>
1549inline const Number *
1551{
1552 static_assert(rank_ == 1,
1553 "This function is only available for rank-1 tensors "
1554 "because higher-rank tensors may not store their elements "
1555 "in a contiguous array.");
1556
1557 return begin_raw() + n_independent_components;
1558}
1559
1560
1561
1562template <int rank_, int dim, typename Number>
1563template <typename OtherNumber>
1566{
1567 // The following loop could be written more concisely using std::copy, but
1568 // that function is only constexpr from C++20 on.
1569 for (unsigned int i = 0; i < dim; ++i)
1570 values[i] = t.values[i];
1571 return *this;
1572}
1573
1574
1575
1576template <int rank_, int dim, typename Number>
1579 Tensor<rank_, dim, Number>::operator=(const Number &d) &
1580{
1582 (void)d;
1583
1584 for (unsigned int i = 0; i < dim; ++i)
1585 values[i] = internal::NumberType<Number>::value(0.0);
1586 return *this;
1587}
1588
1589
1590# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1591template <int rank_, int dim, typename Number>
1594{
1595 for (unsigned int i = 0; i < dim; ++i)
1596 values[i] = other.values[i];
1597 return *this;
1598}
1599
1600
1601
1602template <int rank_, int dim, typename Number>
1605 Tensor<rank_, dim, Number> &&other) noexcept
1606{
1607 for (unsigned int i = 0; i < dim; ++i)
1608 values[i] = other.values[i];
1609 return *this;
1610}
1611# endif
1612
1613
1614template <int rank_, int dim, typename Number>
1615template <typename OtherNumber>
1616constexpr inline bool
1618 const Tensor<rank_, dim, OtherNumber> &p) const
1619{
1620# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1621 Assert(!(std::is_same_v<Number, adouble> ||
1622 std::is_same_v<OtherNumber, adouble>),
1623 ExcMessage(
1624 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1625 "been extended to support advanced branching."));
1626# endif
1627
1628 for (unsigned int i = 0; i < dim; ++i)
1629 if (numbers::values_are_not_equal(values[i], p.values[i]))
1630 return false;
1631 return true;
1632}
1633
1634
1635// At some places in the library, we have Point<0> for formal reasons
1636// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1637// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1638// in the above function that the loop end check always fails, we
1639// implement this function here
1640template <>
1641template <>
1642constexpr inline bool
1644{
1645 return true;
1646}
1647
1648
1649template <int rank_, int dim, typename Number>
1650template <typename OtherNumber>
1651constexpr bool
1653 const Tensor<rank_, dim, OtherNumber> &p) const
1654{
1655 return !((*this) == p);
1656}
1657
1658
1659template <int rank_, int dim, typename Number>
1660template <typename OtherNumber>
1661constexpr inline DEAL_II_ALWAYS_INLINE
1665{
1666 for (unsigned int i = 0; i < dim; ++i)
1667 values[i] += p.values[i];
1668 return *this;
1669}
1670
1671
1672template <int rank_, int dim, typename Number>
1673template <typename OtherNumber>
1674constexpr inline DEAL_II_ALWAYS_INLINE
1678{
1679 for (unsigned int i = 0; i < dim; ++i)
1680 values[i] -= p.values[i];
1681 return *this;
1682}
1683
1684
1685template <int rank_, int dim, typename Number>
1686template <typename OtherNumber>
1687constexpr inline DEAL_II_ALWAYS_INLINE
1689 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1690{
1691 for (unsigned int i = 0; i < dim; ++i)
1692 values[i] *= s;
1693 return *this;
1694}
1695
1696
1697
1698template <int rank_, int dim, typename Number>
1699template <typename OtherNumber>
1700constexpr inline DEAL_II_ALWAYS_INLINE
1702 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1703{
1704 if constexpr (std::is_integral<
1705 typename ProductType<Number, OtherNumber>::type>::value ||
1706 std::is_same_v<Number, Differentiation::SD::Expression>)
1707 {
1708 // recurse over the base objects
1709 for (unsigned int d = 0; d < dim; ++d)
1710 values[d] /= s;
1711 }
1712 else
1713 {
1714 // If we can, avoid division by multiplying by the inverse of the given
1715 // factor:
1716 const Number inverse_factor = Number(1.) / s;
1717 for (unsigned int d = 0; d < dim; ++d)
1718 values[d] *= inverse_factor;
1719 }
1720
1721 return *this;
1722}
1723
1724
1725template <int rank_, int dim, typename Number>
1726constexpr inline DEAL_II_ALWAYS_INLINE
1729{
1731
1732 for (unsigned int i = 0; i < dim; ++i)
1733 tmp.values[i] = -values[i];
1734
1735 return tmp;
1736}
1737
1738
1739template <int rank_, int dim, typename Number>
1742{
1743 // Handle cases of a tensor consisting of just one number more
1744 // efficiently:
1745 if constexpr ((rank_ == 1) && (dim == 1) && std::is_arithmetic_v<Number>)
1746 {
1747 return std::abs(values[0]);
1748 }
1749 else if constexpr ((rank_ == 2) && (dim == 1) && std::is_arithmetic_v<Number>)
1750 {
1751 return std::abs(values[0][0]);
1752 }
1753 else
1754 {
1755 // Otherwise fall back to the naive algorithm of taking the square root of
1756 // the sum of squares.
1757
1758 // Make things work with AD types by letting the compiler look up
1759 // the symbol sqrt in namespace std and in the type-associated
1760 // namespaces
1761 using std::sqrt;
1762 return sqrt(norm_square());
1763 }
1764}
1765
1766
1767template <int rank_, int dim, typename Number>
1771{
1772 if constexpr (dim == 0)
1773 return internal::NumberType<
1774 typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1775 else if constexpr (rank_ == 1)
1776 {
1777 // For rank-1 tensors, the square of the norm is simply the sum of
1778 // squares of the elements:
1781 for (unsigned int i = 1; i < dim; ++i)
1783
1784 return s;
1785 }
1786 else
1787 {
1788 // For higher-rank tensors, the square of the norm is the sum
1789 // of squares of sub-tensors
1791 values[0].norm_square();
1792 for (unsigned int i = 1; i < dim; ++i)
1793 s += values[i].norm_square();
1794
1795 return s;
1796 }
1797}
1798
1799
1800
1801template <int rank_, int dim, typename Number>
1802template <class Iterator>
1803inline void
1804Tensor<rank_, dim, Number>::unroll(const Iterator begin,
1805 const Iterator end) const
1806{
1807 if constexpr (rank_ > 1)
1808 {
1809 // For higher-rank tensors, we recurse to the sub-tensors:
1810 Iterator next = begin;
1811 for (unsigned int i = 0; i < dim; ++i)
1812 {
1813 values[i].unroll(next, end);
1814 std::advance(
1816 }
1817 }
1818 else
1819 {
1820 // For rank-1 tensors, we can simply copy the current elements from
1821 // our linear array into the output range:
1822 (void)end;
1823 Assert(std::distance(begin, end) >= dim,
1824 ExcMessage(
1825 "The provided iterator range must contain at least 'dim' "
1826 "elements."));
1827 std::copy(values.data(), values.data() + values.size(), begin);
1828 }
1829}
1830
1831
1832
1833template <int rank_, int dim, typename Number>
1834constexpr inline unsigned int
1836 const TableIndices<rank_> &indices)
1837{
1838 unsigned int index = 0;
1839 for (int r = 0; r < rank_; ++r)
1840 index = index * dim + indices[r];
1841
1842 return index;
1843}
1844
1845
1846
1847template <int rank_, int dim, typename Number>
1848constexpr inline TableIndices<rank_>
1850{
1851 // Work-around nvcc warning
1852 unsigned int dummy = n_independent_components;
1853 AssertIndexRange(i, dummy);
1854 (void)dummy;
1855
1856 if constexpr (dim == 0)
1857 {
1858 Assert(false,
1859 ExcMessage(
1860 "A tensor with dimension 0 does not store any elements. "
1861 "There is no indexing that can address its elements."));
1862 return {};
1863 }
1864 else
1865 {
1866 TableIndices<rank_> indices;
1867
1868 unsigned int remainder = i;
1869 for (int r = rank_ - 1; r >= 0; --r)
1870 {
1871 indices[r] = remainder % dim;
1872 remainder = remainder / dim;
1873 }
1874 Assert(remainder == 0, ExcInternalError());
1875
1876 return indices;
1877 }
1878}
1879
1880
1881template <int rank_, int dim, typename Number>
1882constexpr inline void
1884{
1885 for (unsigned int i = 0; i < dim; ++i)
1886 values[i] = internal::NumberType<Number>::value(0.0);
1887}
1888
1889
1890template <int rank_, int dim, typename Number>
1891constexpr std::size_t
1893{
1894 return sizeof(Tensor<rank_, dim, Number>);
1895}
1896
1897
1898template <int rank_, int dim, typename Number>
1899template <class Archive>
1900inline void
1901Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1902{
1903 for (int i = 0; i < dim; ++i)
1904 {
1905 ar &values[i];
1906 }
1907}
1908
1909
1910template <int rank_, int dim, typename Number>
1912
1913#endif // DOXYGEN
1914
1915/* ----------------- Non-member functions operating on tensors. ------------ */
1916
1929template <int rank_, int dim, typename Number>
1930inline std::ostream &
1931operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1932{
1933 for (unsigned int i = 0; i < dim; ++i)
1934 {
1935 out << p[i];
1936 if (i != dim - 1)
1937 for (unsigned int j = 0; j < rank_; ++j)
1938 out << ' ';
1939 }
1940
1941 return out;
1942}
1943
1944
1951template <int dim, typename Number>
1952inline std::ostream &
1953operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1954{
1955 out << static_cast<const Number &>(p);
1956 return out;
1957}
1958
1959
1978template <int dim, typename Number, typename Other>
1981 operator*(const Other &object, const Tensor<0, dim, Number> &t)
1982{
1983 return object * static_cast<const Number &>(t);
1984}
1985
1986
1987
1998template <int dim, typename Number, typename Other>
2001 operator*(const Tensor<0, dim, Number> &t, const Other &object)
2002{
2003 return static_cast<const Number &>(t) * object;
2004}
2005
2006
2018template <int dim, typename Number, typename OtherNumber>
2022 const Tensor<0, dim, OtherNumber> &src2)
2023{
2024 return static_cast<const Number &>(src1) *
2025 static_cast<const OtherNumber &>(src2);
2026}
2027
2028
2036template <int dim, typename Number, typename OtherNumber>
2038 Tensor<0,
2039 dim,
2040 typename ProductType<Number,
2041 typename EnableIfScalar<OtherNumber>::type>::type>
2042 operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2043{
2044 return static_cast<const Number &>(t) / factor;
2045}
2046
2047
2055template <int dim, typename Number, typename OtherNumber>
2060{
2061 return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2062}
2063
2064
2072template <int dim, typename Number, typename OtherNumber>
2077{
2078 return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2079}
2080
2081
2094template <int rank, int dim, typename Number, typename OtherNumber>
2096 Tensor<rank,
2097 dim,
2098 typename ProductType<Number,
2099 typename EnableIfScalar<OtherNumber>::type>::type>
2100 operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2101{
2103 tt *= factor;
2104 return tt;
2105}
2106
2107
2120template <int rank, int dim, typename Number, typename OtherNumber>
2122 Tensor<rank,
2123 dim,
2125 OtherNumber>::type>
2126 operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2127{
2128 // simply forward to the operator above
2129 return t * factor;
2130}
2131
2132
2133
2143template <int rank, int dim, typename Number, typename OtherNumber>
2145 Tensor<rank,
2146 dim,
2147 typename ProductType<Number,
2148 typename EnableIfScalar<OtherNumber>::type>::type>
2149 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2150{
2152 tt /= factor;
2153 return tt;
2154}
2155
2156
2166template <int rank, int dim, typename Number, typename OtherNumber>
2176
2177
2187template <int rank, int dim, typename Number, typename OtherNumber>
2197
2204template <int dim, typename Number, typename OtherNumber>
2205inline constexpr DEAL_II_ALWAYS_INLINE
2208 const Tensor<0, dim, OtherNumber> &src2)
2209{
2211
2212 tmp *= src2;
2213
2214 return tmp;
2215}
2216
2233template <int rank, int dim, typename Number, typename OtherNumber>
2234inline constexpr DEAL_II_ALWAYS_INLINE
2238{
2240
2241 for (unsigned int i = 0; i < dim; ++i)
2244
2245 return tmp;
2246}
2247
2292template <int rank_1,
2293 int rank_2,
2294 int dim,
2295 typename Number,
2296 typename OtherNumber,
2297 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2298constexpr inline DEAL_II_ALWAYS_INLINE
2299 typename Tensor<rank_1 + rank_2 - 2,
2300 dim,
2301 typename ProductType<Number, OtherNumber>::type>::tensor_type
2304{
2305 // Treat some common cases separately. Specifically, these are the dot
2306 // product between two rank-1 tensors, and the product between a
2307 // rank-2 tensor and a rank-1 tensor. Both of these lead to a linear
2308 // loop over adjacent memory and can be dealt with efficiently; in the
2309 // latter case (rank-2 times rank-1), we implement things by deferring
2310 // to rank-1 times rank-1 dot products.
2311 if constexpr ((rank_1 == 1) && (rank_2 == 1))
2312 {
2313 // This is a dot product between two rank-1 tensors. Write it out as
2314 // a linear loop:
2315 static_assert(dim > 0, "Tensors cannot have dimension zero.");
2316 typename ProductType<Number, OtherNumber>::type sum = src1[0] * src2[0];
2317 for (unsigned int i = 1; i < dim; ++i)
2318 sum += src1[i] * src2[i];
2319
2320 return sum;
2321 }
2322 else if constexpr ((rank_1 == 2) && (rank_2 == 1))
2323 {
2324 // This is a product between a rank-2 and a rank-1 tensor. This
2325 // corresponds to taking dot products between the rows of the former
2326 // and the latter.
2327 typename Tensor<
2328 rank_1 + rank_2 - 2,
2329 dim,
2330 typename ProductType<Number, OtherNumber>::type>::tensor_type result;
2331 for (unsigned int i = 0; i < dim; ++i)
2332 result[i] += src1[i] * src2;
2333
2334 return result;
2335 }
2336 else
2337 {
2338 // Treat all of the other cases using the more general contraction
2339 // machinery.
2340 typename Tensor<
2341 rank_1 + rank_2 - 2,
2342 dim,
2343 typename ProductType<Number, OtherNumber>::type>::tensor_type result{};
2344
2345 TensorAccessors::internal::
2346 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2347 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2348 TensorAccessors::contract<1, rank_1, rank_2, dim>(result,
2349 src1,
2350 reordered);
2351
2352 return result;
2353 }
2354}
2355
2356
2385template <int index_1,
2386 int index_2,
2387 int rank_1,
2388 int rank_2,
2389 int dim,
2390 typename Number,
2391 typename OtherNumber>
2392constexpr inline DEAL_II_ALWAYS_INLINE
2393 typename Tensor<rank_1 + rank_2 - 2,
2394 dim,
2395 typename ProductType<Number, OtherNumber>::type>::tensor_type
2398{
2399 Assert(0 <= index_1 && index_1 < rank_1,
2400 ExcMessage(
2401 "The specified index_1 must lie within the range [0,rank_1)"));
2402 Assert(0 <= index_2 && index_2 < rank_2,
2403 ExcMessage(
2404 "The specified index_2 must lie within the range [0,rank_2)"));
2405
2406 using namespace TensorAccessors;
2407 using namespace TensorAccessors::internal;
2408
2409 // Reorder index_1 to the end of src1:
2411 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2412
2413 // Reorder index_2 to the end of src2:
2414 const ReorderedIndexView<index_2,
2415 rank_2,
2417 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2418
2419 typename Tensor<rank_1 + rank_2 - 2,
2420 dim,
2421 typename ProductType<Number, OtherNumber>::type>::tensor_type
2422 result{};
2423 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2424 return result;
2425}
2426
2427
2458template <int index_1,
2459 int index_2,
2460 int index_3,
2461 int index_4,
2462 int rank_1,
2463 int rank_2,
2464 int dim,
2465 typename Number,
2466 typename OtherNumber>
2467constexpr inline
2468 typename Tensor<rank_1 + rank_2 - 4,
2469 dim,
2470 typename ProductType<Number, OtherNumber>::type>::tensor_type
2471 double_contract(const Tensor<rank_1, dim, Number> &src1,
2473{
2474 Assert(0 <= index_1 && index_1 < rank_1,
2475 ExcMessage(
2476 "The specified index_1 must lie within the range [0,rank_1)"));
2477 Assert(0 <= index_3 && index_3 < rank_1,
2478 ExcMessage(
2479 "The specified index_3 must lie within the range [0,rank_1)"));
2480 Assert(index_1 != index_3,
2481 ExcMessage("index_1 and index_3 must not be the same"));
2482 Assert(0 <= index_2 && index_2 < rank_2,
2483 ExcMessage(
2484 "The specified index_2 must lie within the range [0,rank_2)"));
2485 Assert(0 <= index_4 && index_4 < rank_2,
2486 ExcMessage(
2487 "The specified index_4 must lie within the range [0,rank_2)"));
2488 Assert(index_2 != index_4,
2489 ExcMessage("index_2 and index_4 must not be the same"));
2490
2491 using namespace TensorAccessors;
2492 using namespace TensorAccessors::internal;
2493
2494 // Reorder index_1 to the end of src1:
2496 reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2497
2498 // Reorder index_2 to the end of src2:
2500 reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2501
2502 // Now, reorder index_3 to the end of src1. We have to make sure to
2503 // preserve the original ordering: index_1 has been removed. If
2504 // index_3 > index_1, we have to use (index_3 - 1) instead:
2506 (index_3 < index_1 ? index_3 : index_3 - 1),
2507 rank_1,
2508 ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2509 reord_3 =
2510 TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2511 index_3 - 1,
2512 rank_1 > (reord_1);
2513
2514 // Now, reorder index_4 to the end of src2. We have to make sure to
2515 // preserve the original ordering: index_2 has been removed. If
2516 // index_4 > index_2, we have to use (index_4 - 1) instead:
2518 (index_4 < index_2 ? index_4 : index_4 - 1),
2519 rank_2,
2521 reord_4 =
2522 TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2523 index_4 - 1,
2524 rank_2 > (reord_2);
2525
2526 typename Tensor<rank_1 + rank_2 - 4,
2527 dim,
2528 typename ProductType<Number, OtherNumber>::type>::tensor_type
2529 result{};
2530 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2531 return result;
2532}
2533
2534
2547template <int rank, int dim, typename Number, typename OtherNumber>
2548constexpr inline DEAL_II_ALWAYS_INLINE
2550 scalar_product(const Tensor<rank, dim, Number> &left,
2551 const Tensor<rank, dim, OtherNumber> &right)
2552{
2554 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2555 return result;
2556}
2557
2558
2576template <template <int, int, typename> class TensorT1,
2577 template <int, int, typename>
2578 class TensorT2,
2579 template <int, int, typename>
2580 class TensorT3,
2581 int rank_1,
2582 int rank_2,
2583 int dim,
2584 typename T1,
2585 typename T2,
2586 typename T3>
2587constexpr inline DEAL_II_ALWAYS_INLINE
2589 contract3(const TensorT1<rank_1, dim, T1> &left,
2590 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2591 const TensorT3<rank_2, dim, T3> &right)
2592{
2593 using return_type =
2595 return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2596 middle,
2597 right);
2598}
2599
2600
2611template <int rank_1,
2612 int rank_2,
2613 int dim,
2614 typename Number,
2615 typename OtherNumber>
2616constexpr inline DEAL_II_ALWAYS_INLINE
2620{
2621 typename Tensor<rank_1 + rank_2,
2622 dim,
2623 typename ProductType<Number, OtherNumber>::type>::tensor_type
2624 result{};
2625 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2626 return result;
2627}
2628
2629
2648template <int dim, typename Number>
2650cross_product_2d(const Tensor<1, dim, Number> &src)
2651{
2652 Assert(dim == 2, ExcInternalError());
2653
2655
2656 result[0] = src[1];
2657 result[1] = -src[0];
2658
2659 return result;
2660}
2661
2662
2672template <int dim, typename Number1, typename Number2>
2673constexpr inline DEAL_II_ALWAYS_INLINE
2675 cross_product_3d(const Tensor<1, dim, Number1> &src1,
2676 const Tensor<1, dim, Number2> &src2)
2677{
2678 Assert(dim == 3, ExcInternalError());
2679
2681
2682 // avoid compiler warnings
2683 constexpr int s0 = 0 % dim;
2684 constexpr int s1 = 1 % dim;
2685 constexpr int s2 = 2 % dim;
2686
2687 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2688 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2689 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2690
2691 return result;
2692}
2693
2694
2708template <int dim, typename Number>
2709constexpr inline DEAL_II_ALWAYS_INLINE Number
2711{
2712 // Compute the determinant using the Laplace expansion of the
2713 // determinant. We expand along the last row.
2714 Number det = internal::NumberType<Number>::value(0.0);
2715
2716 for (unsigned int k = 0; k < dim; ++k)
2717 {
2718 Tensor<2, dim - 1, Number> minor;
2719 for (unsigned int i = 0; i < dim - 1; ++i)
2720 for (unsigned int j = 0; j < dim - 1; ++j)
2721 minor[i][j] = t[i][j < k ? j : j + 1];
2722
2723 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2724
2725 det += t[dim - 1][k] * cofactor;
2726 }
2727
2728 return ((dim % 2 == 0) ? 1. : -1.) * det;
2729}
2730
2736template <typename Number>
2737constexpr DEAL_II_ALWAYS_INLINE Number
2739{
2740 return t[0][0];
2741}
2742
2748template <typename Number>
2749constexpr DEAL_II_ALWAYS_INLINE Number
2751{
2752 // hard-coded for efficiency reasons
2753 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2754}
2755
2761template <typename Number>
2762constexpr DEAL_II_ALWAYS_INLINE Number
2764{
2765 // hard-coded for efficiency reasons
2766 const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2767 internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2768 const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2769 internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2770 const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2771 internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2772 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2773}
2774
2775
2782template <int dim, typename Number>
2783constexpr inline DEAL_II_ALWAYS_INLINE Number
2785{
2786 Number t = d[0][0];
2787 for (unsigned int i = 1; i < dim; ++i)
2788 t += d[i][i];
2789 return t;
2790}
2791
2792
2801template <int dim, typename Number>
2802constexpr inline Tensor<2, dim, Number>
2804{
2805 Number return_tensor[dim][dim];
2806
2807 // if desired, take over the
2808 // inversion of a 4x4 tensor
2809 // from the FullMatrix
2811
2812 return Tensor<2, dim, Number>(return_tensor);
2813}
2814
2815
2816#ifndef DOXYGEN
2817
2818template <typename Number>
2821{
2822 Tensor<2, 1, Number> return_tensor;
2823
2824 return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2825
2826 return return_tensor;
2827}
2828
2829
2830template <typename Number>
2833{
2834 Tensor<2, 2, Number> return_tensor;
2835
2836 const Number inv_det_t = internal::NumberType<Number>::value(
2837 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2838 return_tensor[0][0] = t[1][1];
2839 return_tensor[0][1] = -t[0][1];
2840 return_tensor[1][0] = -t[1][0];
2841 return_tensor[1][1] = t[0][0];
2842 return_tensor *= inv_det_t;
2843
2844 return return_tensor;
2845}
2846
2847template <typename Number>
2850{
2851 Tensor<2, 3, Number> return_tensor;
2852
2853 const auto value = [](const auto &t) {
2855 };
2856
2857 return_tensor[0][0] = value(t[1][1] * t[2][2]) - value(t[1][2] * t[2][1]);
2858 return_tensor[0][1] = value(t[0][2] * t[2][1]) - value(t[0][1] * t[2][2]);
2859 return_tensor[0][2] = value(t[0][1] * t[1][2]) - value(t[0][2] * t[1][1]);
2860 return_tensor[1][0] = value(t[1][2] * t[2][0]) - value(t[1][0] * t[2][2]);
2861 return_tensor[1][1] = value(t[0][0] * t[2][2]) - value(t[0][2] * t[2][0]);
2862 return_tensor[1][2] = value(t[0][2] * t[1][0]) - value(t[0][0] * t[1][2]);
2863 return_tensor[2][0] = value(t[1][0] * t[2][1]) - value(t[1][1] * t[2][0]);
2864 return_tensor[2][1] = value(t[0][1] * t[2][0]) - value(t[0][0] * t[2][1]);
2865 return_tensor[2][2] = value(t[0][0] * t[1][1]) - value(t[0][1] * t[1][0]);
2866
2867 const Number inv_det_t =
2868 value(1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2869 t[0][2] * return_tensor[2][0]));
2870 return_tensor *= inv_det_t;
2871
2872 return return_tensor;
2873}
2874
2875#endif /* DOXYGEN */
2876
2877
2883template <int dim, typename Number>
2886{
2888 for (unsigned int i = 0; i < dim; ++i)
2889 {
2890 tt[i][i] = t[i][i];
2891 for (unsigned int j = i + 1; j < dim; ++j)
2892 {
2893 tt[i][j] = t[j][i];
2894 tt[j][i] = t[i][j];
2895 };
2896 }
2897 return tt;
2898}
2899
2900
2914template <int dim, typename Number>
2915constexpr Tensor<2, dim, Number>
2916adjugate(const Tensor<2, dim, Number> &t)
2917{
2918 return determinant(t) * invert(t);
2919}
2920
2921
2935template <int dim, typename Number>
2936constexpr Tensor<2, dim, Number>
2937cofactor(const Tensor<2, dim, Number> &t)
2938{
2939 return transpose(adjugate(t));
2940}
2941
2942
3006template <int dim, typename Number>
3009
3010
3018template <int dim, typename Number>
3019inline Number
3021{
3022 Number max = internal::NumberType<Number>::value(0.0);
3023 for (unsigned int j = 0; j < dim; ++j)
3024 {
3025 Number sum = internal::NumberType<Number>::value(0.0);
3026 for (unsigned int i = 0; i < dim; ++i)
3027 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3028
3029 if (sum > max)
3030 max = sum;
3031 }
3032
3033 return max;
3034}
3035
3036
3044template <int dim, typename Number>
3045inline Number
3047{
3048 Number max = internal::NumberType<Number>::value(0.0);
3049 for (unsigned int i = 0; i < dim; ++i)
3050 {
3051 Number sum = internal::NumberType<Number>::value(0.0);
3052 for (unsigned int j = 0; j < dim; ++j)
3053 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3054
3055 if (sum > max)
3056 max = sum;
3057 }
3058
3059 return max;
3060}
3061
3067#ifndef DOXYGEN
3068
3069
3070# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3071
3072// Specialization of functions for ADOL-C number types when
3073// the advanced branching feature is used
3074template <int dim>
3075inline adouble
3077{
3078 adouble max = internal::NumberType<adouble>::value(0.0);
3079 for (unsigned int j = 0; j < dim; ++j)
3080 {
3081 adouble sum = internal::NumberType<adouble>::value(0.0);
3082 for (unsigned int i = 0; i < dim; ++i)
3083 sum += fabs(t[i][j]);
3084
3085 condassign(max, (sum > max), sum, max);
3086 }
3087
3088 return max;
3089}
3090
3091
3092template <int dim>
3093inline adouble
3095{
3097 for (unsigned int i = 0; i < dim; ++i)
3098 {
3100 for (unsigned int j = 0; j < dim; ++j)
3101 sum += fabs(t[i][j]);
3102
3103 condassign(max, (sum > max), sum, max);
3104 }
3105
3106 return max;
3107}
3108
3109# endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3110
3111
3112#endif // DOXYGEN
3113
3115
3116#endif
std::size_t size() const
Definition array_view.h:684
Definition point.h:111
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)
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()
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)
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:128
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:2149
constexpr Tensor(const Tensor< 1, dim, Tensor< rank_ - 1, dim, OtherNumber > > &initializer)
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:2042
Number * begin_raw()
static constexpr unsigned int rank
Definition tensor.h:481
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:506
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
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:2021
const Number * begin_raw() const
std::conditional_t< rank_==1, Number[(dim !=0) ? dim :1], typename Tensor< rank_ - 1, dim, Number >::array_type[(dim !=0) ? dim :1]> array_type
Definition tensor.h:522
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:2207
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:2075
friend class Tensor
Definition tensor.h:865
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:2169
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3046
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:1981
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3020
static constexpr unsigned int dimension
Definition tensor.h:476
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:2126
static constexpr std::size_t memory_consumption()
constexpr Tensor & operator=(const Number &d) &
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
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:2058
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:2236
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:2190
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition tensor.h:2001
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:2100
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:497
constexpr Tensor operator-() const
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_DEPRECATED
Definition config.h:205
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:175
#define DEAL_II_CXX23_ASSUME(expr)
Definition config.h:192
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:902
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:910
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:894
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
#define DEAL_II_HOST_DEVICE
Definition numbers.h:30
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition numbers.h:31
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:697
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
Definition numbers.h:588
static constexpr real_type abs_square(const number &x)
Definition numbers.h:579
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:81
std::ostream & operator<<(std::ostream &out, const Tensor< rank_, dim, Number > &p)
Definition tensor.h:1931
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:2207
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3046
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:1981
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3020