Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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#ifdef DEAL_II_WITH_ADOLC
28# include <adolc/adouble.h> // Taped double
29#endif
30
31// boost::serialization::make_array used to be in array.hpp, but was
32// moved to a different file in BOOST 1.64
33#include <boost/version.hpp>
34#if BOOST_VERSION >= 106400
35# include <boost/serialization/array_wrapper.hpp>
36#else
37# include <boost/serialization/array.hpp>
38#endif
39
40
41#include <array>
42#include <cmath>
43#include <ostream>
44#include <type_traits>
45
47
48// Forward declarations:
49#ifndef DOXYGEN
50template <typename ElementType, typename MemorySpace>
51class ArrayView;
52
53template <int dim, typename Number>
55class Point;
56
57template <int rank_, int dim, typename Number = double>
58class Tensor;
59template <typename Number>
60class Vector;
61template <typename number>
62class FullMatrix;
63namespace Differentiation
64{
65 namespace SD
66 {
67 class Expression;
68 }
69} // namespace Differentiation
70#endif
71
72
102template <int dim, typename Number>
103class Tensor<0, dim, Number>
104{
105public:
106 static_assert(dim >= 0,
107 "Tensors must have a dimension greater than or equal to one.");
108
117 static constexpr unsigned int dimension = dim;
118
122 static constexpr unsigned int rank = 0;
123
127 static constexpr unsigned int n_independent_components = 1;
128
138
143 using value_type = Number;
144
150 using array_type = Number;
151
157 constexpr DEAL_II_HOST_DEVICE
159
167 template <typename OtherNumber>
168 constexpr DEAL_II_HOST_DEVICE
170
176 template <typename OtherNumber>
177 constexpr DEAL_II_HOST_DEVICE
178 Tensor(const OtherNumber &initializer);
179
180#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
184 constexpr DEAL_II_HOST_DEVICE
185 Tensor(const Tensor<0, dim, Number> &other);
186
190 constexpr DEAL_II_HOST_DEVICE
191 Tensor(Tensor<0, dim, Number> &&other) noexcept;
192#endif
193
203 Number *
205
215 const Number *
216 begin_raw() const;
217
227 Number *
229
240 const Number *
241 end_raw() const;
242
252 constexpr DEAL_II_HOST_DEVICE
253 operator Number &();
254
263 constexpr DEAL_II_HOST_DEVICE operator const Number &() const;
264
272 template <typename OtherNumber>
273 constexpr DEAL_II_HOST_DEVICE Tensor &
275
276#if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
285 constexpr DEAL_II_HOST_DEVICE Tensor &
287#endif
288
289#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
294 operator=(Tensor<0, dim, Number> &&other) noexcept;
295#endif
296
303 template <typename OtherNumber>
304 constexpr DEAL_II_HOST_DEVICE Tensor &
305 operator=(const OtherNumber &d) &;
306
312 template <typename OtherNumber>
313 constexpr DEAL_II_HOST_DEVICE Tensor &
314 operator=(const OtherNumber &d) && = delete;
315
319 template <typename OtherNumber>
320 constexpr bool
322
326 template <typename OtherNumber>
327 constexpr bool
329
335 template <typename OtherNumber>
336 constexpr DEAL_II_HOST_DEVICE Tensor &
338
344 template <typename OtherNumber>
345 constexpr DEAL_II_HOST_DEVICE Tensor &
347
353 template <typename OtherNumber>
354 constexpr DEAL_II_HOST_DEVICE Tensor &
355 operator*=(const OtherNumber &factor);
356
362 template <typename OtherNumber>
363 constexpr DEAL_II_HOST_DEVICE Tensor &
364 operator/=(const OtherNumber &factor);
365
372 operator-() const;
373
386 constexpr void
388
395 norm() const;
396
404 norm_square() const;
405
413 template <class Iterator>
414 void
415 unroll(const Iterator begin, const Iterator end) const;
416
422 template <class Archive>
423 void
424 serialize(Archive &ar, const unsigned int version);
425
430 using tensor_type = Number;
431
432private:
436 Number value;
437
438 // Allow an arbitrary Tensor to access the underlying values.
439 template <int, int, typename>
440 friend class Tensor;
441};
442
443
444
518template <int rank_, int dim, typename Number>
520{
521public:
522 static_assert(rank_ >= 1,
523 "Tensors must have a rank greater than or equal to one.");
524 static_assert(dim >= 0,
525 "Tensors must have a dimension greater than or equal to zero.");
534 static constexpr unsigned int dimension = dim;
535
539 static constexpr unsigned int rank = rank_;
540
545 static constexpr unsigned int n_independent_components =
546 Tensor<rank_ - 1, dim>::n_independent_components * dim;
547
554 std::conditional_t<rank_ == 1, Number, Tensor<rank_ - 1, dim, Number>>;
555
567 using array_type = std::conditional_t<
568 rank_ == 1,
569 Number[(dim != 0) ? dim : 1],
570 typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1]>;
571
579
585 constexpr DEAL_II_HOST_DEVICE explicit Tensor(const array_type &initializer);
586
599 template <typename ElementType, typename MemorySpace>
600 constexpr DEAL_II_HOST_DEVICE explicit Tensor(
601 const ArrayView<ElementType, MemorySpace> &initializer);
602
610 template <typename OtherNumber>
611 constexpr DEAL_II_HOST_DEVICE
613
617 template <typename OtherNumber>
618 constexpr Tensor(
619 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
620
624 template <typename OtherNumber>
625 constexpr
626 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
627
628#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
632 constexpr Tensor(const Tensor<rank_, dim, Number> &);
633
637 constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
638#endif
639
646 operator[](const unsigned int i);
647
653 constexpr DEAL_II_HOST_DEVICE const value_type &
654 operator[](const unsigned int i) const;
655
659 constexpr const Number &
660 operator[](const TableIndices<rank_> &indices) const;
661
665 constexpr Number &
667
671 Number *
673
677 const Number *
678 begin_raw() const;
679
683 Number *
685
689 const Number *
690 end_raw() const;
691
699 template <typename OtherNumber>
700 constexpr DEAL_II_HOST_DEVICE Tensor &
702
709 constexpr DEAL_II_HOST_DEVICE Tensor &
710 operator=(const Number &d) &;
711
717 constexpr DEAL_II_HOST_DEVICE Tensor &
718 operator=(const Number &d) && = delete;
719
720#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
726
732#endif
733
737 template <typename OtherNumber>
738 constexpr bool
740
744 template <typename OtherNumber>
745 constexpr bool
747
753 template <typename OtherNumber>
754 constexpr DEAL_II_HOST_DEVICE Tensor &
756
762 template <typename OtherNumber>
763 constexpr DEAL_II_HOST_DEVICE Tensor &
765
772 template <typename OtherNumber>
773 constexpr DEAL_II_HOST_DEVICE Tensor &
774 operator*=(const OtherNumber &factor);
775
781 template <typename OtherNumber>
782 constexpr DEAL_II_HOST_DEVICE Tensor &
783 operator/=(const OtherNumber &factor);
784
791 operator-() const;
792
805 constexpr void
807
817 norm() const;
818
825 constexpr DEAL_II_HOST_DEVICE
827 norm_square() const;
828
839 template <typename OtherNumber>
842
853 template <class Iterator>
854 void
855 unroll(const Iterator begin, const Iterator end) const;
856
861 static constexpr DEAL_II_HOST_DEVICE unsigned int
863
870 unrolled_to_component_indices(const unsigned int i);
871
876 static constexpr std::size_t
878
884 template <class Archive>
885 void
886 serialize(Archive &ar, const unsigned int version);
887
893
894private:
900 std::conditional_t<rank_ == 1,
901 std::array<Number, dim>,
902 std::array<Tensor<rank_ - 1, dim, Number>, dim>>
904
911 template <typename ArrayLike, std::size_t... Indices>
912 constexpr DEAL_II_HOST_DEVICE
913 Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
914
915 // Allow an arbitrary Tensor to access the underlying values.
916 template <int, int, typename>
917 friend class Tensor;
918
919 // Point is allowed access to the coordinates. This is supposed to improve
920 // speed.
921 friend class Point<dim, Number>;
922};
923
924
925#ifndef DOXYGEN
926namespace internal
927{
928 // Workaround: The following 4 overloads are necessary to be able to
929 // compile the library with Apple Clang 8 and older. We should remove
930 // these overloads again when we bump the minimal required version to
931 // something later than clang-3.6 / Apple Clang 6.3.
932 template <int rank, int dim, typename T, typename U>
933 struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
934 {
935 using type =
937 };
938
939 template <int rank, int dim, typename T, typename U>
940 struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
941 {
942 using type =
944 };
945
946 template <typename T, int rank, int dim, typename U>
947 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
948 {
949 using type =
951 };
952
953 template <int rank, int dim, typename T, typename U>
954 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
955 {
956 using type =
958 };
959 // end workaround
960
965 template <int rank, int dim, typename T>
966 struct NumberType<Tensor<rank, dim, T>>
967 {
968 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const
971 {
972 return t;
973 }
974
976 value(const T &t)
977 {
979 tmp = t;
980 return tmp;
981 }
982 };
983} // namespace internal
984
985
986/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
987
988
989template <int dim, typename Number>
992 // Some auto-differentiable numbers need explicit
993 // zero initialization such as adtl::adouble.
994 : Tensor{0.0}
995{}
996
997
998
999template <int dim, typename Number>
1000template <typename OtherNumber>
1002Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
1003 : value(internal::NumberType<Number>::value(initializer))
1004{}
1005
1006
1007
1008template <int dim, typename Number>
1009template <typename OtherNumber>
1012 : Tensor{p.value}
1013{}
1014
1015
1016# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1017template <int dim, typename Number>
1020 : value{other.value}
1021{}
1022
1023
1024
1025template <int dim, typename Number>
1028 : value{std::move(other.value)}
1029{}
1030# endif
1031
1032
1033template <int dim, typename Number>
1034inline Number *
1036{
1037 return std::addressof(value);
1038}
1039
1040
1041
1042template <int dim, typename Number>
1043inline const Number *
1045{
1046 return std::addressof(value);
1047}
1048
1049
1050
1051template <int dim, typename Number>
1052inline Number *
1054{
1055 return begin_raw() + n_independent_components;
1056}
1057
1058
1059
1060template <int dim, typename Number>
1061const Number *
1063{
1064 return begin_raw() + n_independent_components;
1065}
1066
1067
1068
1069template <int dim, typename Number>
1072{
1073 Assert(dim != 0,
1074 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1075 return value;
1076}
1077
1078
1079template <int dim, typename Number>
1080constexpr inline DEAL_II_ALWAYS_INLINE
1082{
1083 Assert(dim != 0,
1084 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1085 return value;
1086}
1087
1088
1089
1090template <int dim, typename Number>
1091template <typename OtherNumber>
1094{
1096 return *this;
1097}
1098
1099
1100# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1101template <int dim, typename Number>
1104{
1105 value = p.value;
1106 return *this;
1107}
1108# endif
1109
1110# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1111template <int dim, typename Number>
1114{
1115 value = std::move(other.value);
1116 return *this;
1117}
1118# endif
1119
1120
1121
1122template <int dim, typename Number>
1123template <typename OtherNumber>
1125Tensor<0, dim, Number>::operator=(const OtherNumber &d) &
1126{
1128 return *this;
1129}
1130
1131
1132template <int dim, typename Number>
1133template <typename OtherNumber>
1134constexpr inline bool
1136{
1137# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1138 Assert(!(std::is_same_v<Number, adouble> ||
1139 std::is_same_v<OtherNumber, adouble>),
1140 ExcMessage(
1141 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1142 "been extended to support advanced branching."));
1143# endif
1144
1145 return numbers::values_are_equal(value, p.value);
1146}
1147
1148
1149template <int dim, typename Number>
1150template <typename OtherNumber>
1151constexpr bool
1153{
1154 return !((*this) == p);
1155}
1156
1157
1158template <int dim, typename Number>
1159template <typename OtherNumber>
1162{
1163 value += p.value;
1164 return *this;
1165}
1166
1167
1168template <int dim, typename Number>
1169template <typename OtherNumber>
1172{
1173 value -= p.value;
1174 return *this;
1175}
1176
1177
1178
1179namespace internal
1180{
1181 namespace ComplexWorkaround
1182 {
1183 template <typename Number, typename OtherNumber>
1185 multiply_assign_scalar(Number &val, const OtherNumber &s)
1186 {
1187 val *= s;
1188 }
1189
1190 template <typename Number, typename OtherNumber>
1192 multiply_assign_scalar(std::complex<Number> &val, const OtherNumber &s)
1193 {
1194# if KOKKOS_VERSION >= 30600
1195 KOKKOS_IF_ON_HOST((val *= s;))
1196 KOKKOS_IF_ON_DEVICE(({
1197 (void)val;
1198 (void)s;
1199 Kokkos::abort(
1200 "This function is not implemented for std::complex<Number>!\n");
1201 }))
1202# else
1203# ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1204 val *= s;
1205# else
1206 (void)val;
1207 (void)s;
1208 Kokkos::abort(
1209 "This function is not implemented for std::complex<Number>!\n");
1210# endif
1211# endif
1212 }
1213 } // namespace ComplexWorkaround
1214} // namespace internal
1215
1216
1217template <int dim, typename Number>
1218template <typename OtherNumber>
1220Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1221{
1222 internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1223 return *this;
1224}
1225
1226
1227
1228template <int dim, typename Number>
1229template <typename OtherNumber>
1231Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1232{
1233 value /= s;
1234 return *this;
1235}
1236
1237
1238template <int dim, typename Number>
1241{
1242 return -value;
1243}
1244
1245
1246template <int dim, typename Number>
1249{
1250 Assert(dim != 0,
1251 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1253}
1254
1255
1256template <int dim, typename Number>
1260{
1261 Assert(dim != 0,
1262 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1264}
1265
1266
1267
1268template <int dim, typename Number>
1269constexpr inline void
1271{
1272 // Some auto-differentiable numbers need explicit
1273 // zero initialization.
1275}
1276
1277
1278
1279template <int dim, typename Number>
1280template <class Iterator>
1281inline void
1282Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1283{
1284 (void)end;
1285 AssertDimension(std::distance(begin, end), n_independent_components);
1286 Assert(dim != 0,
1287 ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1288 Assert(std::distance(begin, end) >= 1,
1289 ExcMessage("The provided iterator range must contain at least one "
1290 "element."));
1291 *begin = value;
1292}
1293
1294
1295
1296template <int dim, typename Number>
1297template <class Archive>
1298inline void
1299Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1300{
1301 ar &value;
1302}
1303
1304
1305template <int dim, typename Number>
1307
1308
1309/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1310
1311template <int rank_, int dim, typename Number>
1312template <typename ArrayLike, std::size_t... indices>
1314Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1315 std::index_sequence<indices...>)
1316 // Extract from the 'initializer' a sequence of elements via template
1317 // pack evaluation. This could be as easy as
1318 // values{{ (initializer[indices])... }}
1319 // but of course in practice it is not. The challenge is that if rank>1,
1320 // we want to pass the elements initializer[indices] down to the next
1321 // lower rank tensor for evaluation unchanged. But at the rank==1 level,
1322 // we need to convert to the scalar type 'Number'. This would all be
1323 // relatively straightforward if we could rely on automatic type
1324 // conversion, but for some autodifferentiation types, the conversion
1325 // from the AD to double (i.e., the extraction of a scalar value) is
1326 // not implicit, and we need to call internal::NumberType<Number>::value() --
1327 // but as mentioned, we can only do that for rank==1.
1328 //
1329 // We can achieve all of this by dispatching to a lambda function within
1330 // which we can use a 'if constexpr'.
1331 : values{{([&initializer]() -> value_type {
1332 if constexpr (rank_ == 1)
1333 return internal::NumberType<Number>::value(initializer[indices]);
1334 else
1335 return value_type(initializer[indices]);
1336 }())...}}
1337{
1338 static_assert(sizeof...(indices) == dim,
1339 "dim should match the number of indices");
1340}
1341
1342
1343# ifdef DEAL_II_HAVE_CXX20
1344
1345template <int rank_, int dim, typename Number>
1348 : values(
1349 // In order to initialize the std::array<Number,dim>, we would need a
1350 // brace-enclosed list of length 'dim'. There is no way in C++ to create
1351 // such a list in-place, but we can come up with a lambda function that
1352 // expands such a list via template-pack expansion, and then uses this
1353 // list to initialize a std::array which it then returns.
1354 //
1355 // The trick to come up with such a lambda function is to have a function
1356 // that takes an argument that depends on a template-pack of integers.
1357 // We will call the function with an integer list of length 'dim', and
1358 // in the function itself expand that pack in a way that it serves as
1359 // a brace-enclosed list of initializers for a std::array.
1360 //
1361 // Of course, we do not want to initialize the array with the integers,
1362 // but with zeros. (Or, more correctly, a zero of the element type.)
1363 // The canonical way to do this would be using the comma operator:
1364 // (sequence_element, 0.0)
1365 // returns zero, and
1366 // (sequence, 0.0)...
1367 // returns a list of zeros of the right length. Unfortunately, some
1368 // compilers then warn that the left side of the comma expression has
1369 // no effect -- well, bummer, that was of course exactly the idea.
1370 // We could work around this by using
1371 // (sequence_element * 0.0)
1372 // instead, assuming that the compiler will optimize (known) integer
1373 // times zero to zero, and similarly for (known) integer times times
1374 // default-initialized tensor.
1375 //
1376 // But, instead of relying on compiler optimizations, a better way is
1377 // to simply have another (nested) lambda function that takes the
1378 // integer sequence element as an argument and ignores it, just
1379 // returning a zero instead.
1380 []<std::size_t... I>(
1381 const std::index_sequence<I...> &) constexpr -> decltype(values) {
1382 if constexpr (rank_ == 1)
1383 {
1384 auto get_zero_and_ignore_argument = [](int) {
1386 };
1387 return {{(get_zero_and_ignore_argument(I))...}};
1388 }
1389 else
1390 {
1391 auto get_zero_and_ignore_argument = [](int) {
1392 return Tensor<rank_ - 1, dim, Number>();
1393 };
1394 return {{(get_zero_and_ignore_argument(I))...}};
1395 }
1396 }(std::make_index_sequence<dim>()))
1397{}
1398
1399# else
1400
1401// The C++17 case works in essence the same, except that we can't use a
1402// lambda function with explicit template parameters, i.e., we can't do
1403// []<std::size_t... I>(const std::index_sequence<I...> &)
1404// as above because that's a C++20 feature. Lambda functions in C++17 can
1405// have template packs as arguments, but we need the ability to *name*
1406// that template pack (the 'I' above) and that's not possible in C++17.
1407//
1408// We work around this by moving the lambda function to a global function
1409// and using the traditional template syntax on it.
1410namespace internal
1411{
1412 namespace TensorInitialization
1413 {
1414 template <int rank, int dim, typename Number, std::size_t... I>
1415 constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
1416 make_zero_array(const std::index_sequence<I...> &)
1417 {
1418 static_assert(sizeof...(I) == dim, "This is bad.");
1419
1420 // First peel off the case dim==0. If we don't, some compilers
1421 // will warn below that we define these lambda functions but
1422 // never use them (because the expanded list has zero elements,
1423 // and the get_zero_and_ignore_argument() function is not used...)
1424 if constexpr (dim == 0)
1425 {
1426 return {};
1427 }
1428 else if constexpr (rank == 1)
1429 {
1430 auto get_zero_and_ignore_argument = [](int) {
1432 };
1433 return {{(get_zero_and_ignore_argument(I))...}};
1434 }
1435 else
1436 {
1437 auto get_zero_and_ignore_argument = [](int) {
1438 return Tensor<rank - 1, dim, Number>();
1439 };
1440 return {{(get_zero_and_ignore_argument(I))...}};
1441 }
1442 }
1443 } // namespace TensorInitialization
1444} // namespace internal
1445
1446
1447template <int rank_, int dim, typename Number>
1450 : values(internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
1451 std::make_index_sequence<dim>()))
1452{}
1453
1454
1455# endif
1456
1457
1458template <int rank_, int dim, typename Number>
1460Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1461 : Tensor(initializer, std::make_index_sequence<dim>{})
1462{}
1463
1464
1465
1466template <int rank_, int dim, typename Number>
1467template <typename ElementType, typename MemorySpace>
1470 const ArrayView<ElementType, MemorySpace> &initializer)
1471{
1472 // make nvcc happy
1473 const int my_n_independent_components = n_independent_components;
1474 AssertDimension(initializer.size(), my_n_independent_components);
1475
1476 for (unsigned int i = 0; i < my_n_independent_components; ++i)
1477 (*this)[unrolled_to_component_indices(i)] = initializer[i];
1478}
1479
1480
1481
1482template <int rank_, int dim, typename Number>
1483template <typename OtherNumber>
1486 const Tensor<rank_, dim, OtherNumber> &initializer)
1487 : Tensor(initializer, std::make_index_sequence<dim>{})
1488{}
1489
1490
1491
1492template <int rank_, int dim, typename Number>
1493template <typename OtherNumber>
1494constexpr DEAL_II_ALWAYS_INLINE
1496 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1497 : Tensor(initializer, std::make_index_sequence<dim>{})
1498{}
1499
1500
1501
1502template <int rank_, int dim, typename Number>
1503template <typename OtherNumber>
1505operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1506{
1507 Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> x;
1508 std::copy(values.begin(), values.end(), x.values.begin());
1509 return x;
1510}
1511
1512
1513# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1514template <int rank_, int dim, typename Number>
1515constexpr DEAL_II_ALWAYS_INLINE
1517 : values(other.values)
1518{}
1519
1520
1521
1522template <int rank_, int dim, typename Number>
1523constexpr DEAL_II_ALWAYS_INLINE
1525 : values(std::move(other.values))
1526{}
1527# endif
1528
1529
1530
1531template <int rank_, int dim, typename Number>
1534 Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1535{
1536 Assert(dim != 0,
1537 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1538 AssertIndexRange(i, dim);
1539 DEAL_II_CXX23_ASSUME(i < dim);
1540
1541 return values[i];
1542}
1543
1544
1545template <int rank_, int dim, typename Number>
1546constexpr DEAL_II_ALWAYS_INLINE
1548 Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1549{
1550 Assert(dim != 0,
1551 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1552 AssertIndexRange(i, dim);
1553 DEAL_II_CXX23_ASSUME(i < dim);
1554
1555 return values[i];
1556}
1557
1558
1559template <int rank_, int dim, typename Number>
1560constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1562{
1563 Assert(dim != 0,
1564 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1565
1566 return TensorAccessors::extract<rank_>(*this, indices);
1567}
1568
1569
1570
1571template <int rank_, int dim, typename Number>
1572constexpr inline DEAL_II_ALWAYS_INLINE Number &
1574{
1575 Assert(dim != 0,
1576 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1577
1578 return TensorAccessors::extract<rank_>(*this, indices);
1579}
1580
1581
1582
1583template <int rank_, int dim, typename Number>
1584inline Number *
1586{
1587 return std::addressof(
1588 this->operator[](this->unrolled_to_component_indices(0)));
1589}
1590
1591
1592
1593template <int rank_, int dim, typename Number>
1594inline const Number *
1596{
1597 return std::addressof(
1598 this->operator[](this->unrolled_to_component_indices(0)));
1599}
1600
1601
1602
1603template <int rank_, int dim, typename Number>
1604inline Number *
1606{
1607 return begin_raw() + n_independent_components;
1608}
1609
1610
1611
1612template <int rank_, int dim, typename Number>
1613inline const Number *
1615{
1616 return begin_raw() + n_independent_components;
1617}
1618
1619
1620
1621template <int rank_, int dim, typename Number>
1622template <typename OtherNumber>
1625{
1626 // The following loop could be written more concisely using std::copy, but
1627 // that function is only constexpr from C++20 on.
1628 for (unsigned int i = 0; i < dim; ++i)
1629 values[i] = t.values[i];
1630 return *this;
1631}
1632
1633
1634
1635template <int rank_, int dim, typename Number>
1638 Tensor<rank_, dim, Number>::operator=(const Number &d) &
1639{
1641 (void)d;
1642
1643 for (unsigned int i = 0; i < dim; ++i)
1644 values[i] = internal::NumberType<Number>::value(0.0);
1645 return *this;
1646}
1647
1648
1649# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1650template <int rank_, int dim, typename Number>
1653{
1654 for (unsigned int i = 0; i < dim; ++i)
1655 values[i] = other.values[i];
1656 return *this;
1657}
1658
1659
1660
1661template <int rank_, int dim, typename Number>
1664 Tensor<rank_, dim, Number> &&other) noexcept
1665{
1666 for (unsigned int i = 0; i < dim; ++i)
1667 values[i] = other.values[i];
1668 return *this;
1669}
1670# endif
1671
1672
1673template <int rank_, int dim, typename Number>
1674template <typename OtherNumber>
1675constexpr inline bool
1677 const Tensor<rank_, dim, OtherNumber> &p) const
1678{
1679# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1680 Assert(!(std::is_same_v<Number, adouble> ||
1681 std::is_same_v<OtherNumber, adouble>),
1682 ExcMessage(
1683 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1684 "been extended to support advanced branching."));
1685# endif
1686
1687 for (unsigned int i = 0; i < dim; ++i)
1688 if (numbers::values_are_not_equal(values[i], p.values[i]))
1689 return false;
1690 return true;
1691}
1692
1693
1694// At some places in the library, we have Point<0> for formal reasons
1695// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1696// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1697// in the above function that the loop end check always fails, we
1698// implement this function here
1699template <>
1700template <>
1701constexpr inline bool
1703{
1704 return true;
1705}
1706
1707
1708template <int rank_, int dim, typename Number>
1709template <typename OtherNumber>
1710constexpr bool
1712 const Tensor<rank_, dim, OtherNumber> &p) const
1713{
1714 return !((*this) == p);
1715}
1716
1717
1718template <int rank_, int dim, typename Number>
1719template <typename OtherNumber>
1720constexpr inline DEAL_II_ALWAYS_INLINE
1724{
1725 for (unsigned int i = 0; i < dim; ++i)
1726 values[i] += p.values[i];
1727 return *this;
1728}
1729
1730
1731template <int rank_, int dim, typename Number>
1732template <typename OtherNumber>
1733constexpr inline DEAL_II_ALWAYS_INLINE
1737{
1738 for (unsigned int i = 0; i < dim; ++i)
1739 values[i] -= p.values[i];
1740 return *this;
1741}
1742
1743
1744template <int rank_, int dim, typename Number>
1745template <typename OtherNumber>
1746constexpr inline DEAL_II_ALWAYS_INLINE
1748 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1749{
1750 for (unsigned int i = 0; i < dim; ++i)
1751 values[i] *= s;
1752 return *this;
1753}
1754
1755
1756
1757template <int rank_, int dim, typename Number>
1758template <typename OtherNumber>
1759constexpr inline DEAL_II_ALWAYS_INLINE
1761 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1762{
1763 if constexpr (std::is_integral<
1764 typename ProductType<Number, OtherNumber>::type>::value ||
1765 std::is_same_v<Number, Differentiation::SD::Expression>)
1766 {
1767 // recurse over the base objects
1768 for (unsigned int d = 0; d < dim; ++d)
1769 values[d] /= s;
1770 }
1771 else
1772 {
1773 // If we can, avoid division by multiplying by the inverse of the given
1774 // factor:
1775 const Number inverse_factor = Number(1.) / s;
1776 for (unsigned int d = 0; d < dim; ++d)
1777 values[d] *= inverse_factor;
1778 }
1779
1780 return *this;
1781}
1782
1783
1784template <int rank_, int dim, typename Number>
1785constexpr inline DEAL_II_ALWAYS_INLINE
1788{
1790
1791 for (unsigned int i = 0; i < dim; ++i)
1792 tmp.values[i] = -values[i];
1793
1794 return tmp;
1795}
1796
1797
1798template <int rank_, int dim, typename Number>
1801{
1802 // Handle cases of a tensor consisting of just one number more
1803 // efficiently:
1804 if constexpr ((rank_ == 1) && (dim == 1) && std::is_arithmetic_v<Number>)
1805 {
1806 return std::abs(values[0]);
1807 }
1808 else if constexpr ((rank_ == 2) && (dim == 1) && std::is_arithmetic_v<Number>)
1809 {
1810 return std::abs(values[0][0]);
1811 }
1812 else
1813 {
1814 // Otherwise fall back to the naive algorithm of taking the square root of
1815 // the sum of squares.
1816
1817 // Make things work with AD types by letting the compiler look up
1818 // the symbol sqrt in namespace std and in the type-associated
1819 // namespaces
1820 using std::sqrt;
1821 return sqrt(norm_square());
1822 }
1823}
1824
1825
1826template <int rank_, int dim, typename Number>
1830{
1831 if constexpr (dim == 0)
1832 return internal::NumberType<
1833 typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1834 else if constexpr (rank_ == 1)
1835 {
1836 // For rank-1 tensors, the square of the norm is simply the sum of
1837 // squares of the elements:
1840 for (unsigned int i = 1; i < dim; ++i)
1842
1843 return s;
1844 }
1845 else
1846 {
1847 // For higher-rank tensors, the square of the norm is the sum
1848 // of squares of sub-tensors
1850 values[0].norm_square();
1851 for (unsigned int i = 1; i < dim; ++i)
1852 s += values[i].norm_square();
1853
1854 return s;
1855 }
1856}
1857
1858
1859
1860template <int rank_, int dim, typename Number>
1861template <typename OtherNumber>
1862inline void
1864{
1865 unroll(result.begin(), result.end());
1866}
1867
1868
1869
1870template <int rank_, int dim, typename Number>
1871template <class Iterator>
1872inline void
1873Tensor<rank_, dim, Number>::unroll(const Iterator begin,
1874 const Iterator end) const
1875{
1876 if constexpr (rank_ > 1)
1877 {
1878 // For higher-rank tensors, we recurse to the sub-tensors:
1879 Iterator next = begin;
1880 for (unsigned int i = 0; i < dim; ++i)
1881 {
1882 values[i].unroll(next, end);
1883 std::advance(
1885 }
1886 }
1887 else
1888 {
1889 // For rank-1 tensors, we can simply copy the current elements from
1890 // our linear array into the output range, and then return the
1891 // iterator moved forward by 'dim' elements:
1892 (void)end;
1893 Assert(std::distance(begin, end) >= dim,
1894 ExcMessage(
1895 "The provided iterator range must contain at least 'dim' "
1896 "elements."));
1897 std::copy(std::begin(values), std::end(values), begin);
1898 }
1899}
1900
1901
1902
1903template <int rank_, int dim, typename Number>
1904constexpr inline unsigned int
1906 const TableIndices<rank_> &indices)
1907{
1908 unsigned int index = 0;
1909 for (int r = 0; r < rank_; ++r)
1910 index = index * dim + indices[r];
1911
1912 return index;
1913}
1914
1915
1916
1917template <int rank_, int dim, typename Number>
1918constexpr inline TableIndices<rank_>
1920{
1921 // Work-around nvcc warning
1922 unsigned int dummy = n_independent_components;
1923 AssertIndexRange(i, dummy);
1924 (void)dummy;
1925
1926 if constexpr (dim == 0)
1927 {
1928 Assert(false,
1929 ExcMessage(
1930 "A tensor with dimension 0 does not store any elements. "
1931 "There is no indexing that can address its elements."));
1932 return {};
1933 }
1934 else
1935 {
1936 TableIndices<rank_> indices;
1937
1938 unsigned int remainder = i;
1939 for (int r = rank_ - 1; r >= 0; --r)
1940 {
1941 indices[r] = remainder % dim;
1942 remainder = remainder / dim;
1943 }
1944 Assert(remainder == 0, ExcInternalError());
1945
1946 return indices;
1947 }
1948}
1949
1950
1951template <int rank_, int dim, typename Number>
1952constexpr inline void
1954{
1955 for (unsigned int i = 0; i < dim; ++i)
1956 values[i] = internal::NumberType<Number>::value(0.0);
1957}
1958
1959
1960template <int rank_, int dim, typename Number>
1961constexpr std::size_t
1963{
1964 return sizeof(Tensor<rank_, dim, Number>);
1965}
1966
1967
1968template <int rank_, int dim, typename Number>
1969template <class Archive>
1970inline void
1971Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1972{
1973 if constexpr (rank_ > 1)
1974 ar &values;
1975 else
1976 ar &boost::serialization::make_array(&values[0], dim);
1977}
1978
1979
1980template <int rank_, int dim, typename Number>
1982
1983#endif // DOXYGEN
1984
1985/* ----------------- Non-member functions operating on tensors. ------------ */
1986
1999template <int rank_, int dim, typename Number>
2000inline std::ostream &
2001operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
2002{
2003 for (unsigned int i = 0; i < dim; ++i)
2004 {
2005 out << p[i];
2006 if (i != dim - 1)
2007 for (unsigned int j = 0; j < rank_; ++j)
2008 out << ' ';
2009 }
2010
2011 return out;
2012}
2013
2014
2021template <int dim, typename Number>
2022inline std::ostream &
2023operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
2024{
2025 out << static_cast<const Number &>(p);
2026 return out;
2027}
2028
2029
2048template <int dim, typename Number, typename Other>
2051 operator*(const Other &object, const Tensor<0, dim, Number> &t)
2052{
2053 return object * static_cast<const Number &>(t);
2054}
2055
2056
2057
2068template <int dim, typename Number, typename Other>
2071 operator*(const Tensor<0, dim, Number> &t, const Other &object)
2072{
2073 return static_cast<const Number &>(t) * object;
2074}
2075
2076
2088template <int dim, typename Number, typename OtherNumber>
2092 const Tensor<0, dim, OtherNumber> &src2)
2093{
2094 return static_cast<const Number &>(src1) *
2095 static_cast<const OtherNumber &>(src2);
2096}
2097
2098
2106template <int dim, typename Number, typename OtherNumber>
2108 Tensor<0,
2109 dim,
2110 typename ProductType<Number,
2111 typename EnableIfScalar<OtherNumber>::type>::type>
2112 operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2113{
2114 return static_cast<const Number &>(t) / factor;
2115}
2116
2117
2125template <int dim, typename Number, typename OtherNumber>
2130{
2131 return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2132}
2133
2134
2142template <int dim, typename Number, typename OtherNumber>
2147{
2148 return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2149}
2150
2151
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{
2172 // recurse over the base objects
2174 for (unsigned int d = 0; d < dim; ++d)
2175 tt[d] = t[d] * factor;
2176 return tt;
2177}
2178
2179
2192template <int rank, int dim, typename Number, typename OtherNumber>
2194 Tensor<rank,
2195 dim,
2197 OtherNumber>::type>
2198 operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2199{
2200 // simply forward to the operator above
2201 return t * factor;
2202}
2203
2204
2205
2215template <int rank, int dim, typename Number, typename OtherNumber>
2217 Tensor<rank,
2218 dim,
2219 typename ProductType<Number,
2220 typename EnableIfScalar<OtherNumber>::type>::type>
2221 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2222{
2224 if constexpr (std::is_integral<
2226 {
2227 // recurse over the base objects
2228 for (unsigned int d = 0; d < dim; ++d)
2229 tt[d] = t[d] / factor;
2230 }
2231 else
2232 {
2233 const Number inverse_factor = Number(1.) / factor;
2234 // recurse over the base objects
2235 for (unsigned int d = 0; d < dim; ++d)
2236 tt[d] = t[d] * inverse_factor;
2237 }
2238 return tt;
2239}
2240
2241
2251template <int rank, int dim, typename Number, typename OtherNumber>
2256{
2258
2259 for (unsigned int i = 0; i < dim; ++i)
2260 tmp[i] += q[i];
2261
2262 return tmp;
2263}
2264
2265
2275template <int rank, int dim, typename Number, typename OtherNumber>
2280{
2282
2283 for (unsigned int i = 0; i < dim; ++i)
2284 tmp[i] -= q[i];
2285
2286 return tmp;
2287}
2288
2295template <int dim, typename Number, typename OtherNumber>
2296inline constexpr DEAL_II_ALWAYS_INLINE
2299 const Tensor<0, dim, OtherNumber> &src2)
2300{
2302
2303 tmp *= src2;
2304
2305 return tmp;
2306}
2307
2324template <int rank, int dim, typename Number, typename OtherNumber>
2325inline constexpr DEAL_II_ALWAYS_INLINE
2329{
2331
2332 for (unsigned int i = 0; i < dim; ++i)
2335
2336 return tmp;
2337}
2338
2383template <int rank_1,
2384 int rank_2,
2385 int dim,
2386 typename Number,
2387 typename OtherNumber,
2388 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2389constexpr inline DEAL_II_ALWAYS_INLINE
2390 typename Tensor<rank_1 + rank_2 - 2,
2391 dim,
2392 typename ProductType<Number, OtherNumber>::type>::tensor_type
2395{
2396 typename Tensor<rank_1 + rank_2 - 2,
2397 dim,
2398 typename ProductType<Number, OtherNumber>::type>::tensor_type
2399 result{};
2400
2401 TensorAccessors::internal::
2402 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2403 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2404 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, src1, reordered);
2405
2406 return result;
2407}
2408
2409
2438template <int index_1,
2439 int index_2,
2440 int rank_1,
2441 int rank_2,
2442 int dim,
2443 typename Number,
2444 typename OtherNumber>
2445constexpr inline DEAL_II_ALWAYS_INLINE
2446 typename Tensor<rank_1 + rank_2 - 2,
2447 dim,
2448 typename ProductType<Number, OtherNumber>::type>::tensor_type
2451{
2452 Assert(0 <= index_1 && index_1 < rank_1,
2453 ExcMessage(
2454 "The specified index_1 must lie within the range [0,rank_1)"));
2455 Assert(0 <= index_2 && index_2 < rank_2,
2456 ExcMessage(
2457 "The specified index_2 must lie within the range [0,rank_2)"));
2458
2459 using namespace TensorAccessors;
2460 using namespace TensorAccessors::internal;
2461
2462 // Reorder index_1 to the end of src1:
2464 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2465
2466 // Reorder index_2 to the end of src2:
2467 const ReorderedIndexView<index_2,
2468 rank_2,
2470 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2471
2472 typename Tensor<rank_1 + rank_2 - 2,
2473 dim,
2474 typename ProductType<Number, OtherNumber>::type>::tensor_type
2475 result{};
2476 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2477 return result;
2478}
2479
2480
2511template <int index_1,
2512 int index_2,
2513 int index_3,
2514 int index_4,
2515 int rank_1,
2516 int rank_2,
2517 int dim,
2518 typename Number,
2519 typename OtherNumber>
2520constexpr inline
2521 typename Tensor<rank_1 + rank_2 - 4,
2522 dim,
2523 typename ProductType<Number, OtherNumber>::type>::tensor_type
2524 double_contract(const Tensor<rank_1, dim, Number> &src1,
2526{
2527 Assert(0 <= index_1 && index_1 < rank_1,
2528 ExcMessage(
2529 "The specified index_1 must lie within the range [0,rank_1)"));
2530 Assert(0 <= index_3 && index_3 < rank_1,
2531 ExcMessage(
2532 "The specified index_3 must lie within the range [0,rank_1)"));
2533 Assert(index_1 != index_3,
2534 ExcMessage("index_1 and index_3 must not be the same"));
2535 Assert(0 <= index_2 && index_2 < rank_2,
2536 ExcMessage(
2537 "The specified index_2 must lie within the range [0,rank_2)"));
2538 Assert(0 <= index_4 && index_4 < rank_2,
2539 ExcMessage(
2540 "The specified index_4 must lie within the range [0,rank_2)"));
2541 Assert(index_2 != index_4,
2542 ExcMessage("index_2 and index_4 must not be the same"));
2543
2544 using namespace TensorAccessors;
2545 using namespace TensorAccessors::internal;
2546
2547 // Reorder index_1 to the end of src1:
2549 reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2550
2551 // Reorder index_2 to the end of src2:
2553 reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2554
2555 // Now, reorder index_3 to the end of src1. We have to make sure to
2556 // preserve the original ordering: index_1 has been removed. If
2557 // index_3 > index_1, we have to use (index_3 - 1) instead:
2559 (index_3 < index_1 ? index_3 : index_3 - 1),
2560 rank_1,
2561 ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2562 reord_3 =
2563 TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2564 index_3 - 1,
2565 rank_1 > (reord_1);
2566
2567 // Now, reorder index_4 to the end of src2. We have to make sure to
2568 // preserve the original ordering: index_2 has been removed. If
2569 // index_4 > index_2, we have to use (index_4 - 1) instead:
2571 (index_4 < index_2 ? index_4 : index_4 - 1),
2572 rank_2,
2574 reord_4 =
2575 TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2576 index_4 - 1,
2577 rank_2 > (reord_2);
2578
2579 typename Tensor<rank_1 + rank_2 - 4,
2580 dim,
2581 typename ProductType<Number, OtherNumber>::type>::tensor_type
2582 result{};
2583 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2584 return result;
2585}
2586
2587
2600template <int rank, int dim, typename Number, typename OtherNumber>
2601constexpr inline DEAL_II_ALWAYS_INLINE
2603 scalar_product(const Tensor<rank, dim, Number> &left,
2604 const Tensor<rank, dim, OtherNumber> &right)
2605{
2607 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2608 return result;
2609}
2610
2611
2629template <template <int, int, typename> class TensorT1,
2630 template <int, int, typename>
2631 class TensorT2,
2632 template <int, int, typename>
2633 class TensorT3,
2634 int rank_1,
2635 int rank_2,
2636 int dim,
2637 typename T1,
2638 typename T2,
2639 typename T3>
2640constexpr inline DEAL_II_ALWAYS_INLINE
2642 contract3(const TensorT1<rank_1, dim, T1> &left,
2643 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2644 const TensorT3<rank_2, dim, T3> &right)
2645{
2646 using return_type =
2648 return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2649 middle,
2650 right);
2651}
2652
2653
2664template <int rank_1,
2665 int rank_2,
2666 int dim,
2667 typename Number,
2668 typename OtherNumber>
2669constexpr inline DEAL_II_ALWAYS_INLINE
2673{
2674 typename Tensor<rank_1 + rank_2,
2675 dim,
2676 typename ProductType<Number, OtherNumber>::type>::tensor_type
2677 result{};
2678 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2679 return result;
2680}
2681
2682
2701template <int dim, typename Number>
2703cross_product_2d(const Tensor<1, dim, Number> &src)
2704{
2705 Assert(dim == 2, ExcInternalError());
2706
2708
2709 result[0] = src[1];
2710 result[1] = -src[0];
2711
2712 return result;
2713}
2714
2715
2725template <int dim, typename Number1, typename Number2>
2726constexpr inline DEAL_II_ALWAYS_INLINE
2728 cross_product_3d(const Tensor<1, dim, Number1> &src1,
2729 const Tensor<1, dim, Number2> &src2)
2730{
2731 Assert(dim == 3, ExcInternalError());
2732
2734
2735 // avoid compiler warnings
2736 constexpr int s0 = 0 % dim;
2737 constexpr int s1 = 1 % dim;
2738 constexpr int s2 = 2 % dim;
2739
2740 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2741 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2742 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2743
2744 return result;
2745}
2746
2747
2761template <int dim, typename Number>
2762constexpr inline DEAL_II_ALWAYS_INLINE Number
2764{
2765 // Compute the determinant using the Laplace expansion of the
2766 // determinant. We expand along the last row.
2767 Number det = internal::NumberType<Number>::value(0.0);
2768
2769 for (unsigned int k = 0; k < dim; ++k)
2770 {
2771 Tensor<2, dim - 1, Number> minor;
2772 for (unsigned int i = 0; i < dim - 1; ++i)
2773 for (unsigned int j = 0; j < dim - 1; ++j)
2774 minor[i][j] = t[i][j < k ? j : j + 1];
2775
2776 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2777
2778 det += t[dim - 1][k] * cofactor;
2779 }
2780
2781 return ((dim % 2 == 0) ? 1. : -1.) * det;
2782}
2783
2789template <typename Number>
2790constexpr DEAL_II_ALWAYS_INLINE Number
2792{
2793 return t[0][0];
2794}
2795
2801template <typename Number>
2802constexpr DEAL_II_ALWAYS_INLINE Number
2804{
2805 // hard-coded for efficiency reasons
2806 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2807}
2808
2814template <typename Number>
2815constexpr DEAL_II_ALWAYS_INLINE Number
2817{
2818 // hard-coded for efficiency reasons
2819 const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2820 internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2821 const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2822 internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2823 const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2824 internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2825 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2826}
2827
2828
2835template <int dim, typename Number>
2836constexpr inline DEAL_II_ALWAYS_INLINE Number
2838{
2839 Number t = d[0][0];
2840 for (unsigned int i = 1; i < dim; ++i)
2841 t += d[i][i];
2842 return t;
2843}
2844
2845
2854template <int dim, typename Number>
2855constexpr inline Tensor<2, dim, Number>
2857{
2858 Number return_tensor[dim][dim];
2859
2860 // if desired, take over the
2861 // inversion of a 4x4 tensor
2862 // from the FullMatrix
2864
2865 return Tensor<2, dim, Number>(return_tensor);
2866}
2867
2868
2869#ifndef DOXYGEN
2870
2871template <typename Number>
2874{
2875 Tensor<2, 1, Number> return_tensor;
2876
2877 return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2878
2879 return return_tensor;
2880}
2881
2882
2883template <typename Number>
2886{
2887 Tensor<2, 2, Number> return_tensor;
2888
2889 const Number inv_det_t = internal::NumberType<Number>::value(
2890 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2891 return_tensor[0][0] = t[1][1];
2892 return_tensor[0][1] = -t[0][1];
2893 return_tensor[1][0] = -t[1][0];
2894 return_tensor[1][1] = t[0][0];
2895 return_tensor *= inv_det_t;
2896
2897 return return_tensor;
2898}
2899
2900
2901template <typename Number>
2904{
2905 Tensor<2, 3, Number> return_tensor;
2906
2907 return_tensor[0][0] = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2908 internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2909 return_tensor[0][1] = internal::NumberType<Number>::value(t[0][2] * t[2][1]) -
2910 internal::NumberType<Number>::value(t[0][1] * t[2][2]);
2911 return_tensor[0][2] = internal::NumberType<Number>::value(t[0][1] * t[1][2]) -
2912 internal::NumberType<Number>::value(t[0][2] * t[1][1]);
2913 return_tensor[1][0] = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2914 internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2915 return_tensor[1][1] = internal::NumberType<Number>::value(t[0][0] * t[2][2]) -
2916 internal::NumberType<Number>::value(t[0][2] * t[2][0]);
2917 return_tensor[1][2] = internal::NumberType<Number>::value(t[0][2] * t[1][0]) -
2918 internal::NumberType<Number>::value(t[0][0] * t[1][2]);
2919 return_tensor[2][0] = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2920 internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2921 return_tensor[2][1] = internal::NumberType<Number>::value(t[0][1] * t[2][0]) -
2922 internal::NumberType<Number>::value(t[0][0] * t[2][1]);
2923 return_tensor[2][2] = internal::NumberType<Number>::value(t[0][0] * t[1][1]) -
2924 internal::NumberType<Number>::value(t[0][1] * t[1][0]);
2925 const Number inv_det_t = internal::NumberType<Number>::value(
2926 1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2927 t[0][2] * return_tensor[2][0]));
2928 return_tensor *= inv_det_t;
2929
2930 return return_tensor;
2931}
2932
2933#endif /* DOXYGEN */
2934
2935
2941template <int dim, typename Number>
2944{
2946 for (unsigned int i = 0; i < dim; ++i)
2947 {
2948 tt[i][i] = t[i][i];
2949 for (unsigned int j = i + 1; j < dim; ++j)
2950 {
2951 tt[i][j] = t[j][i];
2952 tt[j][i] = t[i][j];
2953 };
2954 }
2955 return tt;
2956}
2957
2958
2972template <int dim, typename Number>
2973constexpr Tensor<2, dim, Number>
2974adjugate(const Tensor<2, dim, Number> &t)
2975{
2976 return determinant(t) * invert(t);
2977}
2978
2979
2993template <int dim, typename Number>
2994constexpr Tensor<2, dim, Number>
2995cofactor(const Tensor<2, dim, Number> &t)
2996{
2997 return transpose(adjugate(t));
2998}
2999
3000
3064template <int dim, typename Number>
3067
3068
3076template <int dim, typename Number>
3077inline Number
3079{
3080 Number max = internal::NumberType<Number>::value(0.0);
3081 for (unsigned int j = 0; j < dim; ++j)
3082 {
3083 Number sum = internal::NumberType<Number>::value(0.0);
3084 for (unsigned int i = 0; i < dim; ++i)
3085 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3086
3087 if (sum > max)
3088 max = sum;
3089 }
3090
3091 return max;
3092}
3093
3094
3102template <int dim, typename Number>
3103inline Number
3105{
3106 Number max = internal::NumberType<Number>::value(0.0);
3107 for (unsigned int i = 0; i < dim; ++i)
3108 {
3109 Number sum = internal::NumberType<Number>::value(0.0);
3110 for (unsigned int j = 0; j < dim; ++j)
3111 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3112
3113 if (sum > max)
3114 max = sum;
3115 }
3116
3117 return max;
3118}
3119
3125#ifndef DOXYGEN
3126
3127
3128# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3129
3130// Specialization of functions for ADOL-C number types when
3131// the advanced branching feature is used
3132template <int dim>
3133inline adouble
3135{
3136 adouble max = internal::NumberType<adouble>::value(0.0);
3137 for (unsigned int j = 0; j < dim; ++j)
3138 {
3139 adouble sum = internal::NumberType<adouble>::value(0.0);
3140 for (unsigned int i = 0; i < dim; ++i)
3141 sum += fabs(t[i][j]);
3142
3143 condassign(max, (sum > max), sum, max);
3144 }
3145
3146 return max;
3147}
3148
3149
3150template <int dim>
3151inline adouble
3153{
3155 for (unsigned int i = 0; i < dim; ++i)
3156 {
3158 for (unsigned int j = 0; j < dim; ++j)
3159 sum += fabs(t[i][j]);
3160
3161 condassign(max, (sum > max), sum, max);
3162 }
3163
3164 return max;
3165}
3166
3167# endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3168
3169
3170#endif // DOXYGEN
3171
3173
3174#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)
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()
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:137
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:2221
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:2112
Number * begin_raw()
static constexpr unsigned int rank
Definition tensor.h:539
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:554
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
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:2091
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:570
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:2298
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:2145
friend class Tensor
Definition tensor.h:917
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:2254
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3104
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:2051
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3078
static constexpr unsigned int dimension
Definition tensor.h:534
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:2198
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:903
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:2128
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:2327
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:2278
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition tensor.h:2071
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:2170
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:545
constexpr Tensor operator-() const
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
iterator end()
iterator begin()
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_DEPRECATED
Definition config.h:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_CXX23_ASSUME(expr)
Definition config.h:194
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
constexpr void contract(T1 &result, const T2 &left, const T3 &right)
constexpr T1 contract3(const T2 &left, const T3 &middle, const T4 &right)
T sum(const T &t, const MPI_Comm mpi_communicator)
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:928
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:936
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:920
const Iterator & begin
Definition parallel.h:609
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:34
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition numbers.h:36
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:702
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
Definition numbers.h:593
static constexpr real_type abs_square(const number &x)
Definition numbers.h:584
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Definition tensor.cc:81
std::ostream & operator<<(std::ostream &out, const Tensor< rank_, dim, Number > &p)
Definition tensor.h:2001
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:2298
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3104
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:2051
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3078