Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3041-g9c4075ddf4 2025-04-07 16:30: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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 <complex>
35#include <ostream>
36#include <type_traits>
37
39
40// Forward declarations:
41#ifndef DOXYGEN
42template <typename ElementType, typename MemorySpace>
43class ArrayView;
44
45template <int dim, typename Number>
47class Point;
48
49template <int rank_, int dim, typename Number = double>
50class Tensor;
51template <typename Number>
52class Vector;
53template <typename number>
54class FullMatrix;
55namespace Differentiation
56{
57 namespace SD
58 {
59 class Expression;
60 }
61} // namespace Differentiation
62#endif
63
64
94template <int dim, typename Number>
95class Tensor<0, dim, Number>
96{
97public:
98 static_assert(dim >= 0,
99 "Tensors must have a dimension greater than or equal to one.");
100
109 static constexpr unsigned int dimension = dim;
110
114 static constexpr unsigned int rank = 0;
115
119 static constexpr unsigned int n_independent_components = 1;
120
130
135 using value_type = Number;
136
142 using array_type = Number;
143
149 constexpr DEAL_II_HOST_DEVICE
151
159 template <typename OtherNumber>
160 constexpr DEAL_II_HOST_DEVICE
162
168 template <typename OtherNumber>
169 constexpr DEAL_II_HOST_DEVICE
170 Tensor(const OtherNumber &initializer);
171
172#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
176 constexpr DEAL_II_HOST_DEVICE
177 Tensor(const Tensor<0, dim, Number> &other);
178
182 constexpr DEAL_II_HOST_DEVICE
183 Tensor(Tensor<0, dim, Number> &&other) noexcept;
184#endif
185
195 constexpr DEAL_II_HOST_DEVICE
196 operator Number &();
197
206 constexpr DEAL_II_HOST_DEVICE operator const Number &() const;
207
215 template <typename OtherNumber>
216 constexpr DEAL_II_HOST_DEVICE Tensor &
218
219#if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
228 constexpr DEAL_II_HOST_DEVICE Tensor &
230#endif
231
232#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
237 operator=(Tensor<0, dim, Number> &&other) noexcept;
238#endif
239
246 template <typename OtherNumber>
247 constexpr DEAL_II_HOST_DEVICE Tensor &
248 operator=(const OtherNumber &d) &;
249
255 template <typename OtherNumber>
256 constexpr DEAL_II_HOST_DEVICE Tensor &
257 operator=(const OtherNumber &d) && = delete;
258
262 template <typename OtherNumber>
263 constexpr bool
265
269 template <typename OtherNumber>
270 constexpr bool
272
278 template <typename OtherNumber>
279 constexpr DEAL_II_HOST_DEVICE Tensor &
281
287 template <typename OtherNumber>
288 constexpr DEAL_II_HOST_DEVICE Tensor &
290
296 template <typename OtherNumber>
297 constexpr DEAL_II_HOST_DEVICE Tensor &
298 operator*=(const OtherNumber &factor);
299
305 template <typename OtherNumber>
306 constexpr DEAL_II_HOST_DEVICE Tensor &
307 operator/=(const OtherNumber &factor);
308
315 operator-() const;
316
329 constexpr void
331
338 norm() const;
339
347 norm_square() const;
348
356 template <class Iterator>
357 void
358 unroll(const Iterator begin, const Iterator end) const;
359
365 template <class Archive>
366 void
367 serialize(Archive &ar, const unsigned int version);
368
373 using tensor_type = Number;
374
375private:
379 Number value;
380
381 // Allow an arbitrary Tensor to access the underlying values.
382 template <int, int, typename>
383 friend class Tensor;
384};
385
386
387
461template <int rank_, int dim, typename Number>
463{
464public:
465 static_assert(rank_ >= 1,
466 "Tensors must have a rank greater than or equal to one.");
467 static_assert(dim >= 0,
468 "Tensors must have a dimension greater than or equal to zero.");
477 static constexpr unsigned int dimension = dim;
478
482 static constexpr unsigned int rank = rank_;
483
498 static constexpr unsigned int n_independent_components =
499 Tensor<rank_ - 1, dim>::n_independent_components * dim;
500
507 std::conditional_t<rank_ == 1, Number, Tensor<rank_ - 1, dim, Number>>;
508
520 using array_type = std::conditional_t<
521 rank_ == 1,
522 Number[(dim != 0) ? dim : 1],
523 typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1]>;
524
532
538 constexpr DEAL_II_HOST_DEVICE explicit Tensor(const array_type &initializer);
539
552 template <typename ElementType, typename MemorySpace>
553 constexpr DEAL_II_HOST_DEVICE explicit Tensor(
554 const ArrayView<ElementType, MemorySpace> &initializer);
555
563 template <typename OtherNumber>
564 constexpr DEAL_II_HOST_DEVICE
566
570 template <typename OtherNumber>
571 constexpr Tensor(
572 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
573
577 template <typename OtherNumber>
578 constexpr
579 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
580
581#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
585 constexpr Tensor(const Tensor<rank_, dim, Number> &);
586
590 constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
591#endif
592
599 operator[](const unsigned int i);
600
606 constexpr DEAL_II_HOST_DEVICE const value_type &
607 operator[](const unsigned int i) const;
608
612 constexpr const Number &
613 operator[](const TableIndices<rank_> &indices) const;
614
618 constexpr Number &
620
625 Number *
627
632 const Number *
633 begin_raw() const;
634
639 Number *
641
646 const Number *
647 end_raw() const;
648
656 template <typename OtherNumber>
657 constexpr DEAL_II_HOST_DEVICE Tensor &
659
666 constexpr DEAL_II_HOST_DEVICE Tensor &
667 operator=(const Number &d) &;
668
674 constexpr DEAL_II_HOST_DEVICE Tensor &
675 operator=(const Number &d) && = delete;
676
677#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
683
689#endif
690
694 template <typename OtherNumber>
695 constexpr bool
697
701 template <typename OtherNumber>
702 constexpr bool
704
710 template <typename OtherNumber>
711 constexpr DEAL_II_HOST_DEVICE Tensor &
713
719 template <typename OtherNumber>
720 constexpr DEAL_II_HOST_DEVICE Tensor &
722
729 template <typename OtherNumber>
730 constexpr DEAL_II_HOST_DEVICE Tensor &
731 operator*=(const OtherNumber &factor);
732
738 template <typename OtherNumber>
739 constexpr DEAL_II_HOST_DEVICE Tensor &
740 operator/=(const OtherNumber &factor);
741
748 operator-() const;
749
762 constexpr void
764
774 norm() const;
775
782 constexpr DEAL_II_HOST_DEVICE
784 norm_square() const;
785
796 template <class Iterator>
797 void
798 unroll(const Iterator begin, const Iterator end) const;
799
804 static constexpr DEAL_II_HOST_DEVICE unsigned int
806
813 unrolled_to_component_indices(const unsigned int i);
814
819 static constexpr std::size_t
821
827 template <class Archive>
828 void
829 serialize(Archive &ar, const unsigned int version);
830
836
837private:
843#if DEAL_II_KOKKOS_VERSION_GTE(3, 7, 0)
844 std::conditional_t<rank_ == 1,
845 Kokkos::Array<Number, dim>,
846 Kokkos::Array<Tensor<rank_ - 1, dim, Number>, dim>>
847#else
848 std::conditional_t<rank_ == 1,
849 std::array<Number, dim>,
850 std::array<Tensor<rank_ - 1, dim, Number>, dim>>
851#endif
853
860 template <typename ArrayLike, std::size_t... Indices>
861 constexpr DEAL_II_HOST_DEVICE
862 Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
863
864 // Allow an arbitrary Tensor to access the underlying values.
865 template <int, int, typename>
866 friend class Tensor;
867
868 // Point is allowed access to the coordinates. This is supposed to improve
869 // speed.
870 friend class Point<dim, Number>;
871};
872
873
874#ifndef DOXYGEN
875namespace internal
876{
877 // Workaround: The following 4 overloads are necessary to be able to
878 // compile the library with Apple Clang 8 and older. We should remove
879 // these overloads again when we bump the minimal required version to
880 // something later than clang-3.6 / Apple Clang 6.3.
881 template <int rank, int dim, typename T, typename U>
882 struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
883 {
884 using type =
886 };
887
888 template <int rank, int dim, typename T, typename U>
889 struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
890 {
891 using type =
893 };
894
895 template <typename T, int rank, int dim, typename U>
896 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
897 {
898 using type =
900 };
901
902 template <int rank, int dim, typename T, typename U>
903 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
904 {
905 using type =
907 };
908 // end workaround
909
914 template <int rank, int dim, typename T>
915 struct NumberType<Tensor<rank, dim, T>>
916 {
917 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const
920 {
921 return t;
922 }
923
925 value(const T &t)
926 {
928 tmp = t;
929 return tmp;
930 }
931 };
932} // namespace internal
933
934
935/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
936
937
938template <int dim, typename Number>
941 // Some auto-differentiable numbers need explicit
942 // zero initialization such as adtl::adouble.
943 : Tensor{0.0}
944{}
945
946
947
948template <int dim, typename Number>
949template <typename OtherNumber>
951Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
952 : value(internal::NumberType<Number>::value(initializer))
953{}
954
955
956
957template <int dim, typename Number>
958template <typename OtherNumber>
961 : Tensor{p.value}
962{}
963
964
965# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
966template <int dim, typename Number>
969 : value{other.value}
970{}
971
972
973
974template <int dim, typename Number>
977 : value{std::move(other.value)}
978{}
979# endif
980
981
982
983template <int dim, typename Number>
986{
987 Assert(dim != 0,
988 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
989 return value;
990}
991
992
993template <int dim, typename Number>
994constexpr inline DEAL_II_ALWAYS_INLINE
996{
997 Assert(dim != 0,
998 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
999 return value;
1000}
1001
1002
1003
1004template <int dim, typename Number>
1005template <typename OtherNumber>
1008{
1010 return *this;
1011}
1012
1013
1014# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1015template <int dim, typename Number>
1018{
1019 value = p.value;
1020 return *this;
1021}
1022# endif
1023
1024# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1025template <int dim, typename Number>
1028{
1029 value = std::move(other.value);
1030 return *this;
1031}
1032# endif
1033
1034
1035
1036template <int dim, typename Number>
1037template <typename OtherNumber>
1039Tensor<0, dim, Number>::operator=(const OtherNumber &d) &
1040{
1042 return *this;
1043}
1044
1045
1046template <int dim, typename Number>
1047template <typename OtherNumber>
1048constexpr inline bool
1050{
1051# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1052 Assert(!(std::is_same_v<Number, adouble> ||
1053 std::is_same_v<OtherNumber, adouble>),
1054 ExcMessage(
1055 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1056 "been extended to support advanced branching."));
1057# endif
1058
1059 return numbers::values_are_equal(value, p.value);
1060}
1061
1062
1063template <int dim, typename Number>
1064template <typename OtherNumber>
1065constexpr bool
1067{
1068 return !((*this) == p);
1069}
1070
1071
1072template <int dim, typename Number>
1073template <typename OtherNumber>
1076{
1077 value += p.value;
1078 return *this;
1079}
1080
1081
1082template <int dim, typename Number>
1083template <typename OtherNumber>
1086{
1087 value -= p.value;
1088 return *this;
1089}
1090
1091
1092
1093namespace internal
1094{
1095 namespace ComplexWorkaround
1096 {
1097 template <typename Number, typename OtherNumber>
1099 multiply_assign_scalar(Number &val, const OtherNumber &s)
1100 {
1101 val *= s;
1102 }
1103
1104 template <typename Number, typename OtherNumber>
1106 multiply_assign_scalar(std::complex<Number> &val, const OtherNumber &s)
1107 {
1108# if DEAL_II_KOKKOS_VERSION_GTE(3, 6, 0)
1109 KOKKOS_IF_ON_HOST((val *= s;))
1110 KOKKOS_IF_ON_DEVICE(({
1111 (void)val;
1112 (void)s;
1113 Kokkos::abort(
1114 "This function is not implemented for std::complex<Number>!\n");
1115 }))
1116# else
1117# ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1118 val *= s;
1119# else
1120 (void)val;
1121 (void)s;
1122 Kokkos::abort(
1123 "This function is not implemented for std::complex<Number>!\n");
1124# endif
1125# endif
1126 }
1127 } // namespace ComplexWorkaround
1128} // namespace internal
1129
1130
1131template <int dim, typename Number>
1132template <typename OtherNumber>
1134Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1135{
1136 internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1137 return *this;
1138}
1139
1140
1141
1142template <int dim, typename Number>
1143template <typename OtherNumber>
1145Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1146{
1147 value /= s;
1148 return *this;
1149}
1150
1151
1152template <int dim, typename Number>
1155{
1156 return -value;
1157}
1158
1159
1160template <int dim, typename Number>
1163{
1164 Assert(dim != 0,
1165 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1167}
1168
1169
1170template <int dim, typename Number>
1174{
1175 Assert(dim != 0,
1176 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1178}
1179
1180
1181
1182template <int dim, typename Number>
1183constexpr inline void
1185{
1186 // Some auto-differentiable numbers need explicit
1187 // zero initialization.
1189}
1190
1191
1192
1193template <int dim, typename Number>
1194template <class Iterator>
1195inline void
1196Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1197{
1198 (void)end;
1199 AssertDimension(std::distance(begin, end), n_independent_components);
1200 Assert(dim != 0,
1201 ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1202 Assert(std::distance(begin, end) >= 1,
1203 ExcMessage("The provided iterator range must contain at least one "
1204 "element."));
1205 *begin = value;
1206}
1207
1208
1209
1210template <int dim, typename Number>
1211template <class Archive>
1212inline void
1213Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1214{
1215 ar &value;
1216}
1217
1218
1219template <int dim, typename Number>
1221
1222
1223/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1224
1225template <int rank_, int dim, typename Number>
1226template <typename ArrayLike, std::size_t... indices>
1228Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1229 std::index_sequence<indices...>)
1230 // Extract from the 'initializer' a sequence of elements via template
1231 // pack evaluation. This could be as easy as
1232 // values{{ (initializer[indices])... }}
1233 // but of course in practice it is not. The challenge is that if rank>1,
1234 // we want to pass the elements initializer[indices] down to the next
1235 // lower rank tensor for evaluation unchanged. But at the rank==1 level,
1236 // we need to convert to the scalar type 'Number'. This would all be
1237 // relatively straightforward if we could rely on automatic type
1238 // conversion, but for some autodifferentiation types, the conversion
1239 // from the AD to double (i.e., the extraction of a scalar value) is
1240 // not implicit, and we need to call internal::NumberType<Number>::value() --
1241 // but as mentioned, we can only do that for rank==1.
1242 //
1243 // We can achieve all of this by dispatching to a lambda function within
1244 // which we can use a 'if constexpr'.
1245 : values{{([&initializer]() -> value_type {
1246 if constexpr (rank_ == 1)
1247 return internal::NumberType<Number>::value(initializer[indices]);
1248 else
1249 return value_type(initializer[indices]);
1250 }())...}}
1251{
1252 static_assert(sizeof...(indices) == dim,
1253 "dim should match the number of indices");
1254}
1255
1256
1257# if defined(DEAL_II_HAVE_CXX20) && !defined(__NVCC__)
1258
1259template <int rank_, int dim, typename Number>
1262 : values(
1263 // In order to initialize the Kokkos::Array<Number,dim>, we would need a
1264 // brace-enclosed list of length 'dim'. There is no way in C++ to create
1265 // such a list in-place, but we can come up with a lambda function that
1266 // expands such a list via template-pack expansion, and then uses this
1267 // list to initialize a Kokkos::Array which it then returns.
1268 //
1269 // The trick to come up with such a lambda function is to have a function
1270 // that takes an argument that depends on a template-pack of integers.
1271 // We will call the function with an integer list of length 'dim', and
1272 // in the function itself expand that pack in a way that it serves as
1273 // a brace-enclosed list of initializers for a Kokkos::Array.
1274 //
1275 // Of course, we do not want to initialize the array with the integers,
1276 // but with zeros. (Or, more correctly, a zero of the element type.)
1277 // The canonical way to do this would be using the comma operator:
1278 // (sequence_element, 0.0)
1279 // returns zero, and
1280 // (sequence, 0.0)...
1281 // returns a list of zeros of the right length. Unfortunately, some
1282 // compilers then warn that the left side of the comma expression has
1283 // no effect -- well, bummer, that was of course exactly the idea.
1284 // We could work around this by using
1285 // (sequence_element * 0.0)
1286 // instead, assuming that the compiler will optimize (known) integer
1287 // times zero to zero, and similarly for (known) integer times times
1288 // default-initialized tensor.
1289 //
1290 // But, instead of relying on compiler optimizations, a better way is
1291 // to simply have another (nested) lambda function that takes the
1292 // integer sequence element as an argument and ignores it, just
1293 // returning a zero instead.
1294 []<std::size_t... I>(
1295 const std::index_sequence<I...> &) constexpr -> decltype(values) {
1296 if constexpr (dim == 0)
1297 {
1298 return {};
1299 }
1300 else if constexpr (rank_ == 1)
1301 {
1302 auto get_zero_and_ignore_argument = [](int) {
1304 };
1305 return {{(get_zero_and_ignore_argument(I))...}};
1306 }
1307 else
1308 {
1309 auto get_zero_and_ignore_argument = [](int) {
1310 return Tensor<rank_ - 1, dim, Number>();
1311 };
1312 return {{(get_zero_and_ignore_argument(I))...}};
1313 }
1314 }(std::make_index_sequence<dim>()))
1315{}
1316
1317# else
1318
1319// The C++17 case works in essence the same, except that we can't use a
1320// lambda function with explicit template parameters, i.e., we can't do
1321// []<std::size_t... I>(const std::index_sequence<I...> &)
1322// as above because that's a C++20 feature. Lambda functions in C++17 can
1323// have template packs as arguments, but we need the ability to *name*
1324// that template pack (the 'I' above) and that's not possible in C++17.
1325//
1326// We work around this by moving the lambda function to a global function
1327// and using the traditional template syntax on it.
1328namespace internal
1329{
1330 namespace TensorInitialization
1331 {
1332 template <int rank, int dim, typename Number, std::size_t... I>
1333# if DEAL_II_KOKKOS_VERSION_GTE(3, 7, 0)
1334 constexpr Kokkos::Array<typename Tensor<rank, dim, Number>::value_type, dim>
1335# else
1336 constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
1337# endif
1338 make_zero_array(const std::index_sequence<I...> &)
1339 {
1340 static_assert(sizeof...(I) == dim, "This is bad.");
1341
1342 // First peel off the case dim==0. If we don't, some compilers
1343 // will warn below that we define these lambda functions but
1344 // never use them (because the expanded list has zero elements,
1345 // and the get_zero_and_ignore_argument() function is not used...)
1346 if constexpr (dim == 0)
1347 {
1348 return {};
1349 }
1350 else if constexpr (rank == 1)
1351 {
1352 auto get_zero_and_ignore_argument = [](int) {
1354 };
1355 return {{(get_zero_and_ignore_argument(I))...}};
1356 }
1357 else
1358 {
1359 auto get_zero_and_ignore_argument = [](int) {
1360 return Tensor<rank - 1, dim, Number>();
1361 };
1362 return {{(get_zero_and_ignore_argument(I))...}};
1363 }
1364 }
1365 } // namespace TensorInitialization
1366} // namespace internal
1367
1368
1369template <int rank_, int dim, typename Number>
1372 : values(internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
1373 std::make_index_sequence<dim>()))
1374{}
1375
1376
1377# endif
1378
1379
1380template <int rank_, int dim, typename Number>
1382Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1383 : Tensor(initializer, std::make_index_sequence<dim>{})
1384{}
1385
1386
1387
1388template <int rank_, int dim, typename Number>
1389template <typename ElementType, typename MemorySpace>
1392 const ArrayView<ElementType, MemorySpace> &initializer)
1393{
1394 // make nvcc happy
1395 const int my_n_independent_components = n_independent_components;
1396 AssertDimension(initializer.size(), my_n_independent_components);
1397
1398 for (unsigned int i = 0; i < my_n_independent_components; ++i)
1399 (*this)[unrolled_to_component_indices(i)] = initializer[i];
1400}
1401
1402
1403
1404template <int rank_, int dim, typename Number>
1405template <typename OtherNumber>
1408 const Tensor<rank_, dim, OtherNumber> &initializer)
1409 : Tensor(initializer, std::make_index_sequence<dim>{})
1410{}
1411
1412
1413
1414template <int rank_, int dim, typename Number>
1415template <typename OtherNumber>
1416constexpr DEAL_II_ALWAYS_INLINE
1418 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1419 : Tensor(initializer, std::make_index_sequence<dim>{})
1420{}
1421
1422
1423
1424template <int rank_, int dim, typename Number>
1425template <typename OtherNumber>
1427operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1428{
1429 Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> x;
1430 std::copy(values.data(), values.data() + values.size(), x.values.data());
1431 return x;
1432}
1433
1434
1435# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1436template <int rank_, int dim, typename Number>
1437constexpr DEAL_II_ALWAYS_INLINE
1439 : values(other.values)
1440{}
1441
1442
1443
1444template <int rank_, int dim, typename Number>
1445constexpr DEAL_II_ALWAYS_INLINE
1447 : values(std::move(other.values))
1448{}
1449# endif
1450
1451
1452
1453template <int rank_, int dim, typename Number>
1456 Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1457{
1458 Assert(dim != 0,
1459 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1460 AssertIndexRange(i, dim);
1461 DEAL_II_CXX23_ASSUME(i < dim);
1462
1463 return values[i];
1464}
1465
1466
1467template <int rank_, int dim, typename Number>
1468constexpr DEAL_II_ALWAYS_INLINE
1470 Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1471{
1472 Assert(dim != 0,
1473 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1474 AssertIndexRange(i, dim);
1475 DEAL_II_CXX23_ASSUME(i < dim);
1476
1477 return values[i];
1478}
1479
1480
1481template <int rank_, int dim, typename Number>
1482constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1484{
1485 Assert(dim != 0,
1486 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1487
1488 return TensorAccessors::extract<rank_>(*this, indices);
1489}
1490
1491
1492
1493template <int rank_, int dim, typename Number>
1494constexpr inline DEAL_II_ALWAYS_INLINE Number &
1496{
1497 Assert(dim != 0,
1498 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1499
1500 return TensorAccessors::extract<rank_>(*this, indices);
1501}
1502
1503
1504
1505template <int rank_, int dim, typename Number>
1506inline Number *
1508{
1509 static_assert(rank_ == 1,
1510 "This function is only available for rank-1 tensors "
1511 "because higher-rank tensors may not store their elements "
1512 "in a contiguous array.");
1513
1514 return std::addressof(
1515 this->operator[](this->unrolled_to_component_indices(0)));
1516}
1517
1518
1519
1520template <int rank_, int dim, typename Number>
1521inline const Number *
1523{
1524 static_assert(rank_ == 1,
1525 "This function is only available for rank-1 tensors "
1526 "because higher-rank tensors may not store their elements "
1527 "in a contiguous array.");
1528
1529 return std::addressof(
1530 this->operator[](this->unrolled_to_component_indices(0)));
1531}
1532
1533
1534
1535template <int rank_, int dim, typename Number>
1536inline Number *
1538{
1539 static_assert(rank_ == 1,
1540 "This function is only available for rank-1 tensors "
1541 "because higher-rank tensors may not store their elements "
1542 "in a contiguous array.");
1543
1544 return begin_raw() + n_independent_components;
1545}
1546
1547
1548
1549template <int rank_, int dim, typename Number>
1550inline const Number *
1552{
1553 static_assert(rank_ == 1,
1554 "This function is only available for rank-1 tensors "
1555 "because higher-rank tensors may not store their elements "
1556 "in a contiguous array.");
1557
1558 return begin_raw() + n_independent_components;
1559}
1560
1561
1562
1563template <int rank_, int dim, typename Number>
1564template <typename OtherNumber>
1567{
1568 // The following loop could be written more concisely using std::copy, but
1569 // that function is only constexpr from C++20 on.
1570 for (unsigned int i = 0; i < dim; ++i)
1571 values[i] = t.values[i];
1572 return *this;
1573}
1574
1575
1576
1577template <int rank_, int dim, typename Number>
1580 Tensor<rank_, dim, Number>::operator=(const Number &d) &
1581{
1583 (void)d;
1584
1585 for (unsigned int i = 0; i < dim; ++i)
1586 values[i] = internal::NumberType<Number>::value(0.0);
1587 return *this;
1588}
1589
1590
1591# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1592template <int rank_, int dim, typename Number>
1595{
1596 for (unsigned int i = 0; i < dim; ++i)
1597 values[i] = other.values[i];
1598 return *this;
1599}
1600
1601
1602
1603template <int rank_, int dim, typename Number>
1606 Tensor<rank_, dim, Number> &&other) noexcept
1607{
1608 for (unsigned int i = 0; i < dim; ++i)
1609 values[i] = other.values[i];
1610 return *this;
1611}
1612# endif
1613
1614
1615template <int rank_, int dim, typename Number>
1616template <typename OtherNumber>
1617constexpr inline bool
1619 const Tensor<rank_, dim, OtherNumber> &p) const
1620{
1621# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1622 Assert(!(std::is_same_v<Number, adouble> ||
1623 std::is_same_v<OtherNumber, adouble>),
1624 ExcMessage(
1625 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1626 "been extended to support advanced branching."));
1627# endif
1628
1629 for (unsigned int i = 0; i < dim; ++i)
1630 if (numbers::values_are_not_equal(values[i], p.values[i]))
1631 return false;
1632 return true;
1633}
1634
1635
1636// At some places in the library, we have Point<0> for formal reasons
1637// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1638// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1639// in the above function that the loop end check always fails, we
1640// implement this function here
1641template <>
1642template <>
1643constexpr inline bool
1645{
1646 return true;
1647}
1648
1649
1650template <int rank_, int dim, typename Number>
1651template <typename OtherNumber>
1652constexpr bool
1654 const Tensor<rank_, dim, OtherNumber> &p) const
1655{
1656 return !((*this) == p);
1657}
1658
1659
1660template <int rank_, int dim, typename Number>
1661template <typename OtherNumber>
1662constexpr inline DEAL_II_ALWAYS_INLINE
1666{
1667 for (unsigned int i = 0; i < dim; ++i)
1668 values[i] += p.values[i];
1669 return *this;
1670}
1671
1672
1673template <int rank_, int dim, typename Number>
1674template <typename OtherNumber>
1675constexpr inline DEAL_II_ALWAYS_INLINE
1679{
1680 for (unsigned int i = 0; i < dim; ++i)
1681 values[i] -= p.values[i];
1682 return *this;
1683}
1684
1685
1686template <int rank_, int dim, typename Number>
1687template <typename OtherNumber>
1688constexpr inline DEAL_II_ALWAYS_INLINE
1690 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1691{
1692 for (unsigned int i = 0; i < dim; ++i)
1693 values[i] *= s;
1694 return *this;
1695}
1696
1697
1698
1699template <int rank_, int dim, typename Number>
1700template <typename OtherNumber>
1701constexpr inline DEAL_II_ALWAYS_INLINE
1703 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1704{
1705 if constexpr (std::is_integral_v<
1707 std::is_same_v<Number, Differentiation::SD::Expression>)
1708 {
1709 // recurse over the base objects
1710 for (unsigned int d = 0; d < dim; ++d)
1711 values[d] /= s;
1712 }
1713 else
1714 {
1715 // If we can, avoid division by multiplying by the inverse of the given
1716 // factor:
1717 const Number inverse_factor = Number(1.) / s;
1718 for (unsigned int d = 0; d < dim; ++d)
1719 values[d] *= inverse_factor;
1720 }
1721
1722 return *this;
1723}
1724
1725
1726template <int rank_, int dim, typename Number>
1727constexpr inline DEAL_II_ALWAYS_INLINE
1730{
1732
1733 for (unsigned int i = 0; i < dim; ++i)
1734 tmp.values[i] = -values[i];
1735
1736 return tmp;
1737}
1738
1739
1740template <int rank_, int dim, typename Number>
1743{
1744 // Handle cases of a tensor consisting of just one number more
1745 // efficiently:
1746 if constexpr ((rank_ == 1) && (dim == 1) && std::is_arithmetic_v<Number>)
1747 {
1748 return std::abs(values[0]);
1749 }
1750 else if constexpr ((rank_ == 2) && (dim == 1) && std::is_arithmetic_v<Number>)
1751 {
1752 return std::abs(values[0][0]);
1753 }
1754 else
1755 {
1756 // Otherwise fall back to the naive algorithm of taking the square root of
1757 // the sum of squares.
1758
1759 // Make things work with AD types by letting the compiler look up
1760 // the symbol sqrt in namespace std and in the type-associated
1761 // namespaces
1762 using std::sqrt;
1763 return sqrt(norm_square());
1764 }
1765}
1766
1767
1768template <int rank_, int dim, typename Number>
1772{
1773 if constexpr (dim == 0)
1774 return internal::NumberType<
1775 typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1776 else if constexpr (rank_ == 1)
1777 {
1778 // For rank-1 tensors, the square of the norm is simply the sum of
1779 // squares of the elements:
1782 for (unsigned int i = 1; i < dim; ++i)
1784
1785 return s;
1786 }
1787 else
1788 {
1789 // For higher-rank tensors, the square of the norm is the sum
1790 // of squares of sub-tensors
1792 values[0].norm_square();
1793 for (unsigned int i = 1; i < dim; ++i)
1794 s += values[i].norm_square();
1795
1796 return s;
1797 }
1798}
1799
1800
1801
1802template <int rank_, int dim, typename Number>
1803template <class Iterator>
1804inline void
1805Tensor<rank_, dim, Number>::unroll(const Iterator begin,
1806 const Iterator end) const
1807{
1808 if constexpr (rank_ > 1)
1809 {
1810 // For higher-rank tensors, we recurse to the sub-tensors:
1811 Iterator next = begin;
1812 for (unsigned int i = 0; i < dim; ++i)
1813 {
1814 values[i].unroll(next, end);
1815 std::advance(
1817 }
1818 }
1819 else
1820 {
1821 // For rank-1 tensors, we can simply copy the current elements from
1822 // our linear array into the output range:
1823 (void)end;
1824 Assert(std::distance(begin, end) >= dim,
1825 ExcMessage(
1826 "The provided iterator range must contain at least 'dim' "
1827 "elements."));
1828 std::copy(values.data(), values.data() + values.size(), begin);
1829 }
1830}
1831
1832
1833
1834template <int rank_, int dim, typename Number>
1835constexpr inline unsigned int
1837 const TableIndices<rank_> &indices)
1838{
1839 unsigned int index = 0;
1840 for (int r = 0; r < rank_; ++r)
1841 index = index * dim + indices[r];
1842
1843 return index;
1844}
1845
1846
1847
1848template <int rank_, int dim, typename Number>
1849constexpr inline TableIndices<rank_>
1851{
1852 // Work-around nvcc warning
1853 unsigned int dummy = n_independent_components;
1854 AssertIndexRange(i, dummy);
1855 (void)dummy;
1856
1857 if constexpr (dim == 0)
1858 {
1859 Assert(false,
1860 ExcMessage(
1861 "A tensor with dimension 0 does not store any elements. "
1862 "There is no indexing that can address its elements."));
1863 return {};
1864 }
1865 else
1866 {
1867 TableIndices<rank_> indices;
1868
1869 unsigned int remainder = i;
1870 for (int r = rank_ - 1; r >= 0; --r)
1871 {
1872 indices[r] = remainder % dim;
1873 remainder = remainder / dim;
1874 }
1875 Assert(remainder == 0, ExcInternalError());
1876
1877 return indices;
1878 }
1879}
1880
1881
1882template <int rank_, int dim, typename Number>
1883constexpr inline void
1885{
1886 for (unsigned int i = 0; i < dim; ++i)
1887 values[i] = internal::NumberType<Number>::value(0.0);
1888}
1889
1890
1891template <int rank_, int dim, typename Number>
1892constexpr std::size_t
1894{
1895 return sizeof(Tensor<rank_, dim, Number>);
1896}
1897
1898
1899template <int rank_, int dim, typename Number>
1900template <class Archive>
1901inline void
1902Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1903{
1904 for (int i = 0; i < dim; ++i)
1905 {
1906 ar &values[i];
1907 }
1908}
1909
1910
1911template <int rank_, int dim, typename Number>
1913
1914#endif // DOXYGEN
1915
1916/* ----------------- Non-member functions operating on tensors. ------------ */
1917
1930template <int rank_, int dim, typename Number>
1931inline std::ostream &
1932operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1933{
1934 for (unsigned int i = 0; i < dim; ++i)
1935 {
1936 out << p[i];
1937 if (i != dim - 1)
1938 for (unsigned int j = 0; j < rank_; ++j)
1939 out << ' ';
1940 }
1941
1942 return out;
1943}
1944
1945
1952template <int dim, typename Number>
1953inline std::ostream &
1954operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1955{
1956 out << static_cast<const Number &>(p);
1957 return out;
1958}
1959
1960
1979template <int dim, typename Number, typename Other>
1982 operator*(const Other &object, const Tensor<0, dim, Number> &t)
1983{
1984 return object * static_cast<const Number &>(t);
1985}
1986
1987
1988
1999template <int dim, typename Number, typename Other>
2002 operator*(const Tensor<0, dim, Number> &t, const Other &object)
2003{
2004 return static_cast<const Number &>(t) * object;
2005}
2006
2007
2019template <int dim, typename Number, typename OtherNumber>
2023 const Tensor<0, dim, OtherNumber> &src2)
2024{
2025 return static_cast<const Number &>(src1) *
2026 static_cast<const OtherNumber &>(src2);
2027}
2028
2029
2037template <int dim, typename Number, typename OtherNumber>
2039 Tensor<0,
2040 dim,
2041 typename ProductType<Number,
2042 typename EnableIfScalar<OtherNumber>::type>::type>
2043 operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2044{
2045 return static_cast<const Number &>(t) / factor;
2046}
2047
2048
2056template <int dim, typename Number, typename OtherNumber>
2061{
2062 return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2063}
2064
2065
2073template <int dim, typename Number, typename OtherNumber>
2078{
2079 return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2080}
2081
2082
2095template <int rank, int dim, typename Number, typename OtherNumber>
2097 Tensor<rank,
2098 dim,
2099 typename ProductType<Number,
2100 typename EnableIfScalar<OtherNumber>::type>::type>
2101 operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2102{
2104 tt *= factor;
2105 return tt;
2106}
2107
2108
2121template <int rank, int dim, typename Number, typename OtherNumber>
2123 Tensor<rank,
2124 dim,
2126 OtherNumber>::type>
2127 operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2128{
2129 // simply forward to the operator above
2130 return t * factor;
2131}
2132
2133
2134
2144template <int rank, int dim, typename Number, typename OtherNumber>
2146 Tensor<rank,
2147 dim,
2148 typename ProductType<Number,
2149 typename EnableIfScalar<OtherNumber>::type>::type>
2150 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2151{
2153 tt /= factor;
2154 return tt;
2155}
2156
2157
2167template <int rank, int dim, typename Number, typename OtherNumber>
2177
2178
2188template <int rank, int dim, typename Number, typename OtherNumber>
2198
2205template <int dim, typename Number, typename OtherNumber>
2206inline constexpr DEAL_II_ALWAYS_INLINE
2209 const Tensor<0, dim, OtherNumber> &src2)
2210{
2212
2213 tmp *= src2;
2214
2215 return tmp;
2216}
2217
2234template <int rank, int dim, typename Number, typename OtherNumber>
2235inline constexpr DEAL_II_ALWAYS_INLINE
2239{
2241
2242 for (unsigned int i = 0; i < dim; ++i)
2245
2246 return tmp;
2247}
2248
2293template <int rank_1,
2294 int rank_2,
2295 int dim,
2296 typename Number,
2297 typename OtherNumber,
2298 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2299constexpr inline DEAL_II_ALWAYS_INLINE
2300 typename Tensor<rank_1 + rank_2 - 2,
2301 dim,
2302 typename ProductType<Number, OtherNumber>::type>::tensor_type
2305{
2306 // Treat some common cases separately. Specifically, these are the dot
2307 // product between two rank-1 tensors, and the product between a
2308 // rank-2 tensor and a rank-1 tensor. Both of these lead to a linear
2309 // loop over adjacent memory and can be dealt with efficiently; in the
2310 // latter case (rank-2 times rank-1), we implement things by deferring
2311 // to rank-1 times rank-1 dot products.
2312 if constexpr ((rank_1 == 1) && (rank_2 == 1))
2313 {
2314 // This is a dot product between two rank-1 tensors. Write it out as
2315 // a linear loop:
2316 static_assert(dim > 0, "Tensors cannot have dimension zero.");
2317 typename ProductType<Number, OtherNumber>::type sum = src1[0] * src2[0];
2318 for (unsigned int i = 1; i < dim; ++i)
2319 sum += src1[i] * src2[i];
2320
2321 return sum;
2322 }
2323 else if constexpr ((rank_1 == 2) && (rank_2 == 1))
2324 {
2325 // This is a product between a rank-2 and a rank-1 tensor. This
2326 // corresponds to taking dot products between the rows of the former
2327 // and the latter.
2328 typename Tensor<
2329 rank_1 + rank_2 - 2,
2330 dim,
2331 typename ProductType<Number, OtherNumber>::type>::tensor_type result;
2332 for (unsigned int i = 0; i < dim; ++i)
2333 result[i] += src1[i] * src2;
2334
2335 return result;
2336 }
2337 else
2338 {
2339 // Treat all of the other cases using the more general contraction
2340 // machinery.
2341 typename Tensor<
2342 rank_1 + rank_2 - 2,
2343 dim,
2344 typename ProductType<Number, OtherNumber>::type>::tensor_type result{};
2345
2346 TensorAccessors::internal::
2347 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2348 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2349 TensorAccessors::contract<1, rank_1, rank_2, dim>(result,
2350 src1,
2351 reordered);
2352
2353 return result;
2354 }
2355}
2356
2357
2386template <int index_1,
2387 int index_2,
2388 int rank_1,
2389 int rank_2,
2390 int dim,
2391 typename Number,
2392 typename OtherNumber>
2393constexpr inline DEAL_II_ALWAYS_INLINE
2394 typename Tensor<rank_1 + rank_2 - 2,
2395 dim,
2396 typename ProductType<Number, OtherNumber>::type>::tensor_type
2399{
2400 Assert(0 <= index_1 && index_1 < rank_1,
2401 ExcMessage(
2402 "The specified index_1 must lie within the range [0,rank_1)"));
2403 Assert(0 <= index_2 && index_2 < rank_2,
2404 ExcMessage(
2405 "The specified index_2 must lie within the range [0,rank_2)"));
2406
2407 using namespace TensorAccessors;
2408 using namespace TensorAccessors::internal;
2409
2410 // Reorder index_1 to the end of src1:
2412 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2413
2414 // Reorder index_2 to the end of src2:
2415 const ReorderedIndexView<index_2,
2416 rank_2,
2418 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2419
2420 typename Tensor<rank_1 + rank_2 - 2,
2421 dim,
2422 typename ProductType<Number, OtherNumber>::type>::tensor_type
2423 result{};
2424 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2425 return result;
2426}
2427
2428
2459template <int index_1,
2460 int index_2,
2461 int index_3,
2462 int index_4,
2463 int rank_1,
2464 int rank_2,
2465 int dim,
2466 typename Number,
2467 typename OtherNumber>
2468constexpr inline
2469 typename Tensor<rank_1 + rank_2 - 4,
2470 dim,
2471 typename ProductType<Number, OtherNumber>::type>::tensor_type
2472 double_contract(const Tensor<rank_1, dim, Number> &src1,
2474{
2475 Assert(0 <= index_1 && index_1 < rank_1,
2476 ExcMessage(
2477 "The specified index_1 must lie within the range [0,rank_1)"));
2478 Assert(0 <= index_3 && index_3 < rank_1,
2479 ExcMessage(
2480 "The specified index_3 must lie within the range [0,rank_1)"));
2481 Assert(index_1 != index_3,
2482 ExcMessage("index_1 and index_3 must not be the same"));
2483 Assert(0 <= index_2 && index_2 < rank_2,
2484 ExcMessage(
2485 "The specified index_2 must lie within the range [0,rank_2)"));
2486 Assert(0 <= index_4 && index_4 < rank_2,
2487 ExcMessage(
2488 "The specified index_4 must lie within the range [0,rank_2)"));
2489 Assert(index_2 != index_4,
2490 ExcMessage("index_2 and index_4 must not be the same"));
2491
2492 using namespace TensorAccessors;
2493 using namespace TensorAccessors::internal;
2494
2495 // Reorder index_1 to the end of src1:
2497 reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2498
2499 // Reorder index_2 to the end of src2:
2501 reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2502
2503 // Now, reorder index_3 to the end of src1. We have to make sure to
2504 // preserve the original ordering: index_1 has been removed. If
2505 // index_3 > index_1, we have to use (index_3 - 1) instead:
2507 (index_3 < index_1 ? index_3 : index_3 - 1),
2508 rank_1,
2509 ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2510 reord_3 =
2511 TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2512 index_3 - 1,
2513 rank_1 > (reord_1);
2514
2515 // Now, reorder index_4 to the end of src2. We have to make sure to
2516 // preserve the original ordering: index_2 has been removed. If
2517 // index_4 > index_2, we have to use (index_4 - 1) instead:
2519 (index_4 < index_2 ? index_4 : index_4 - 1),
2520 rank_2,
2522 reord_4 =
2523 TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2524 index_4 - 1,
2525 rank_2 > (reord_2);
2526
2527 typename Tensor<rank_1 + rank_2 - 4,
2528 dim,
2529 typename ProductType<Number, OtherNumber>::type>::tensor_type
2530 result{};
2531 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2532 return result;
2533}
2534
2535
2548template <int rank, int dim, typename Number, typename OtherNumber>
2549constexpr inline DEAL_II_ALWAYS_INLINE
2551 scalar_product(const Tensor<rank, dim, Number> &left,
2552 const Tensor<rank, dim, OtherNumber> &right)
2553{
2555 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2556 return result;
2557}
2558
2559
2577template <template <int, int, typename> class TensorT1,
2578 template <int, int, typename>
2579 class TensorT2,
2580 template <int, int, typename>
2581 class TensorT3,
2582 int rank_1,
2583 int rank_2,
2584 int dim,
2585 typename T1,
2586 typename T2,
2587 typename T3>
2588constexpr inline DEAL_II_ALWAYS_INLINE
2590 contract3(const TensorT1<rank_1, dim, T1> &left,
2591 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2592 const TensorT3<rank_2, dim, T3> &right)
2593{
2594 using return_type =
2596 return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2597 middle,
2598 right);
2599}
2600
2601
2612template <int rank_1,
2613 int rank_2,
2614 int dim,
2615 typename Number,
2616 typename OtherNumber>
2617constexpr inline DEAL_II_ALWAYS_INLINE
2621{
2622 typename Tensor<rank_1 + rank_2,
2623 dim,
2624 typename ProductType<Number, OtherNumber>::type>::tensor_type
2625 result{};
2626 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2627 return result;
2628}
2629
2630
2649template <int dim, typename Number>
2651cross_product_2d(const Tensor<1, dim, Number> &src)
2652{
2653 Assert(dim == 2, ExcInternalError());
2654
2656
2657 result[0] = src[1];
2658 result[1] = -src[0];
2659
2660 return result;
2661}
2662
2663
2673template <int dim, typename Number1, typename Number2>
2674constexpr inline DEAL_II_ALWAYS_INLINE
2676 cross_product_3d(const Tensor<1, dim, Number1> &src1,
2677 const Tensor<1, dim, Number2> &src2)
2678{
2679 Assert(dim == 3, ExcInternalError());
2680
2682
2683 // avoid compiler warnings
2684 constexpr int s0 = 0 % dim;
2685 constexpr int s1 = 1 % dim;
2686 constexpr int s2 = 2 % dim;
2687
2688 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2689 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2690 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2691
2692 return result;
2693}
2694
2695
2709template <int dim, typename Number>
2710constexpr inline DEAL_II_ALWAYS_INLINE Number
2712{
2713 // Compute the determinant using the Laplace expansion of the
2714 // determinant. We expand along the last row.
2715 Number det = internal::NumberType<Number>::value(0.0);
2716
2717 for (unsigned int k = 0; k < dim; ++k)
2718 {
2719 Tensor<2, dim - 1, Number> minor;
2720 for (unsigned int i = 0; i < dim - 1; ++i)
2721 for (unsigned int j = 0; j < dim - 1; ++j)
2722 minor[i][j] = t[i][j < k ? j : j + 1];
2723
2724 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2725
2726 det += t[dim - 1][k] * cofactor;
2727 }
2728
2729 return ((dim % 2 == 0) ? 1. : -1.) * det;
2730}
2731
2737template <typename Number>
2738constexpr DEAL_II_ALWAYS_INLINE Number
2740{
2741 return t[0][0];
2742}
2743
2749template <typename Number>
2750constexpr DEAL_II_ALWAYS_INLINE Number
2752{
2753 // hard-coded for efficiency reasons
2754 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2755}
2756
2762template <typename Number>
2763constexpr DEAL_II_ALWAYS_INLINE Number
2765{
2766 // hard-coded for efficiency reasons
2767 const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2768 internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2769 const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2770 internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2771 const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2772 internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2773 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2774}
2775
2776
2783template <int dim, typename Number>
2784constexpr inline DEAL_II_ALWAYS_INLINE Number
2786{
2787 Number t = d[0][0];
2788 for (unsigned int i = 1; i < dim; ++i)
2789 t += d[i][i];
2790 return t;
2791}
2792
2793
2802template <int dim, typename Number>
2803constexpr inline Tensor<2, dim, Number>
2805{
2806 Number return_tensor[dim][dim];
2807
2808 // if desired, take over the
2809 // inversion of a 4x4 tensor
2810 // from the FullMatrix
2812
2813 return Tensor<2, dim, Number>(return_tensor);
2814}
2815
2816
2817#ifndef DOXYGEN
2818
2819template <typename Number>
2822{
2823 Tensor<2, 1, Number> return_tensor;
2824
2825 return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2826
2827 return return_tensor;
2828}
2829
2830
2831template <typename Number>
2834{
2835 Tensor<2, 2, Number> return_tensor;
2836
2837 const Number inv_det_t = internal::NumberType<Number>::value(
2838 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2839 return_tensor[0][0] = t[1][1];
2840 return_tensor[0][1] = -t[0][1];
2841 return_tensor[1][0] = -t[1][0];
2842 return_tensor[1][1] = t[0][0];
2843 return_tensor *= inv_det_t;
2844
2845 return return_tensor;
2846}
2847
2848template <typename Number>
2851{
2852 Tensor<2, 3, Number> return_tensor;
2853
2854 const auto value = [](const auto &t) {
2856 };
2857
2858 return_tensor[0][0] = value(t[1][1] * t[2][2]) - value(t[1][2] * t[2][1]);
2859 return_tensor[0][1] = value(t[0][2] * t[2][1]) - value(t[0][1] * t[2][2]);
2860 return_tensor[0][2] = value(t[0][1] * t[1][2]) - value(t[0][2] * t[1][1]);
2861 return_tensor[1][0] = value(t[1][2] * t[2][0]) - value(t[1][0] * t[2][2]);
2862 return_tensor[1][1] = value(t[0][0] * t[2][2]) - value(t[0][2] * t[2][0]);
2863 return_tensor[1][2] = value(t[0][2] * t[1][0]) - value(t[0][0] * t[1][2]);
2864 return_tensor[2][0] = value(t[1][0] * t[2][1]) - value(t[1][1] * t[2][0]);
2865 return_tensor[2][1] = value(t[0][1] * t[2][0]) - value(t[0][0] * t[2][1]);
2866 return_tensor[2][2] = value(t[0][0] * t[1][1]) - value(t[0][1] * t[1][0]);
2867
2868 const Number inv_det_t =
2869 value(1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2870 t[0][2] * return_tensor[2][0]));
2871 return_tensor *= inv_det_t;
2872
2873 return return_tensor;
2874}
2875
2876#endif /* DOXYGEN */
2877
2878
2884template <int dim, typename Number>
2887{
2889 for (unsigned int i = 0; i < dim; ++i)
2890 {
2891 tt[i][i] = t[i][i];
2892 for (unsigned int j = i + 1; j < dim; ++j)
2893 {
2894 tt[i][j] = t[j][i];
2895 tt[j][i] = t[i][j];
2896 };
2897 }
2898 return tt;
2899}
2900
2901
2915template <int dim, typename Number>
2916constexpr Tensor<2, dim, Number>
2917adjugate(const Tensor<2, dim, Number> &t)
2918{
2919 return determinant(t) * invert(t);
2920}
2921
2922
2936template <int dim, typename Number>
2937constexpr Tensor<2, dim, Number>
2938cofactor(const Tensor<2, dim, Number> &t)
2939{
2940 return transpose(adjugate(t));
2941}
2942
2943
3007template <int dim, typename Number>
3010
3011
3019template <int dim, typename Number>
3020inline Number
3022{
3023 Number max = internal::NumberType<Number>::value(0.0);
3024 for (unsigned int j = 0; j < dim; ++j)
3025 {
3026 Number sum = internal::NumberType<Number>::value(0.0);
3027 for (unsigned int i = 0; i < dim; ++i)
3028 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3029
3030 if (sum > max)
3031 max = sum;
3032 }
3033
3034 return max;
3035}
3036
3037
3045template <int dim, typename Number>
3046inline Number
3048{
3049 Number max = internal::NumberType<Number>::value(0.0);
3050 for (unsigned int i = 0; i < dim; ++i)
3051 {
3052 Number sum = internal::NumberType<Number>::value(0.0);
3053 for (unsigned int j = 0; j < dim; ++j)
3054 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3055
3056 if (sum > max)
3057 max = sum;
3058 }
3059
3060 return max;
3061}
3062
3068#ifndef DOXYGEN
3069
3070
3071# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3072
3073// Specialization of functions for ADOL-C number types when
3074// the advanced branching feature is used
3075template <int dim>
3076inline adouble
3078{
3079 adouble max = internal::NumberType<adouble>::value(0.0);
3080 for (unsigned int j = 0; j < dim; ++j)
3081 {
3082 adouble sum = internal::NumberType<adouble>::value(0.0);
3083 for (unsigned int i = 0; i < dim; ++i)
3084 sum += fabs(t[i][j]);
3085
3086 condassign(max, (sum > max), sum, max);
3087 }
3088
3089 return max;
3090}
3091
3092
3093template <int dim>
3094inline adouble
3096{
3098 for (unsigned int i = 0; i < dim; ++i)
3099 {
3101 for (unsigned int j = 0; j < dim; ++j)
3102 sum += fabs(t[i][j]);
3103
3104 condassign(max, (sum > max), sum, max);
3105 }
3106
3107 return max;
3108}
3109
3110# endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3111
3112
3113#endif // DOXYGEN
3114
3116
3117#endif
std::size_t size() const
Definition array_view.h:689
Definition point.h:113
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:129
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:2150
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, OtherNumber >::type > operator-(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition tensor.h:2076
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:2043
Number * begin_raw()
static constexpr unsigned int rank
Definition tensor.h:482
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:507
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:2022
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:523
constexpr bool operator!=(const Tensor< rank_, dim, OtherNumber > &) const
constexpr Tensor()
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< 0, dim, Number > &src1, const Tensor< 0, dim, OtherNumber > &src2)
Definition tensor.h:2208
constexpr value_type & operator[](const unsigned int i)
Number * end_raw()
friend class Tensor
Definition tensor.h:866
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:2170
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3047
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:1982
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3021
static constexpr unsigned int dimension
Definition tensor.h:477
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:2127
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:852
constexpr Tensor< 0, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< 0, dim, Number > &p, const Tensor< 0, dim, OtherNumber > &q)
Definition tensor.h:2059
constexpr Tensor & operator+=(const Tensor< rank_, dim, OtherNumber > &)
const Number * end_raw() const
constexpr Number & operator[](const TableIndices< rank_ > &indices)
constexpr Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > schur_product(const Tensor< rank, dim, Number > &src1, const Tensor< rank, dim, OtherNumber > &src2)
Definition tensor.h:2237
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:2191
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition tensor.h:2002
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:2101
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:498
constexpr Tensor operator-() const
constexpr Tensor & operator=(const Tensor< rank_, dim, OtherNumber > &rhs)
#define DEAL_II_ALWAYS_INLINE
Definition config.h:161
#define DEAL_II_DEPRECATED
Definition config.h:280
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:242
#define DEAL_II_CXX23_ASSUME(expr)
Definition config.h:259
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition config.h:167
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)
constexpr char T
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr internal::ReorderedIndexView< index, rank, T > reordered_index_view(T &t)
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:888
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:896
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:880
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static constexpr const T & value(const T &t)
Definition numbers.h:683
decltype(std::declval< T >() *std::declval< U >()) type
static real_type abs(const number &x)
Definition numbers.h:574
static constexpr real_type abs_square(const number &x)
Definition numbers.h:565
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:1932
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:2208
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3047
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:1982
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3021