Reference documentation for deal.II version GIT relicensing-468-ge3ce94fd06 2024-04-23 15:40: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 constexpr DEAL_II_HOST_DEVICE
204 operator Number &();
205
214 constexpr DEAL_II_HOST_DEVICE operator const Number &() const;
215
223 template <typename OtherNumber>
224 constexpr DEAL_II_HOST_DEVICE Tensor &
226
227#if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
236 constexpr DEAL_II_HOST_DEVICE Tensor &
238#endif
239
240#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
245 operator=(Tensor<0, dim, Number> &&other) noexcept;
246#endif
247
254 template <typename OtherNumber>
255 constexpr DEAL_II_HOST_DEVICE Tensor &
256 operator=(const OtherNumber &d) &;
257
263 template <typename OtherNumber>
264 constexpr DEAL_II_HOST_DEVICE Tensor &
265 operator=(const OtherNumber &d) && = delete;
266
270 template <typename OtherNumber>
271 constexpr bool
273
277 template <typename OtherNumber>
278 constexpr bool
280
286 template <typename OtherNumber>
287 constexpr DEAL_II_HOST_DEVICE Tensor &
289
295 template <typename OtherNumber>
296 constexpr DEAL_II_HOST_DEVICE Tensor &
298
304 template <typename OtherNumber>
305 constexpr DEAL_II_HOST_DEVICE Tensor &
306 operator*=(const OtherNumber &factor);
307
313 template <typename OtherNumber>
314 constexpr DEAL_II_HOST_DEVICE Tensor &
315 operator/=(const OtherNumber &factor);
316
323 operator-() const;
324
337 constexpr void
339
346 norm() const;
347
355 norm_square() const;
356
364 template <class Iterator>
365 void
366 unroll(const Iterator begin, const Iterator end) const;
367
373 template <class Archive>
374 void
375 serialize(Archive &ar, const unsigned int version);
376
381 using tensor_type = Number;
382
383private:
387 Number value;
388
389 // Allow an arbitrary Tensor to access the underlying values.
390 template <int, int, typename>
391 friend class Tensor;
392};
393
394
395
469template <int rank_, int dim, typename Number>
471{
472public:
473 static_assert(rank_ >= 1,
474 "Tensors must have a rank greater than or equal to one.");
475 static_assert(dim >= 0,
476 "Tensors must have a dimension greater than or equal to zero.");
485 static constexpr unsigned int dimension = dim;
486
490 static constexpr unsigned int rank = rank_;
491
506 static constexpr unsigned int n_independent_components =
507 Tensor<rank_ - 1, dim>::n_independent_components * dim;
508
515 std::conditional_t<rank_ == 1, Number, Tensor<rank_ - 1, dim, Number>>;
516
528 using array_type = std::conditional_t<
529 rank_ == 1,
530 Number[(dim != 0) ? dim : 1],
531 typename Tensor<rank_ - 1, dim, Number>::array_type[(dim != 0) ? dim : 1]>;
532
540
546 constexpr DEAL_II_HOST_DEVICE explicit Tensor(const array_type &initializer);
547
560 template <typename ElementType, typename MemorySpace>
561 constexpr DEAL_II_HOST_DEVICE explicit Tensor(
562 const ArrayView<ElementType, MemorySpace> &initializer);
563
571 template <typename OtherNumber>
572 constexpr DEAL_II_HOST_DEVICE
574
578 template <typename OtherNumber>
579 constexpr Tensor(
580 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
581
585 template <typename OtherNumber>
586 constexpr
587 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const;
588
589#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
593 constexpr Tensor(const Tensor<rank_, dim, Number> &);
594
598 constexpr Tensor(Tensor<rank_, dim, Number> &&) noexcept;
599#endif
600
607 operator[](const unsigned int i);
608
614 constexpr DEAL_II_HOST_DEVICE const value_type &
615 operator[](const unsigned int i) const;
616
620 constexpr const Number &
621 operator[](const TableIndices<rank_> &indices) const;
622
626 constexpr Number &
628
632 DEAL_II_DEPRECATED_EARLY
633 Number *
635
639 DEAL_II_DEPRECATED_EARLY
640 const Number *
641 begin_raw() const;
642
646 DEAL_II_DEPRECATED_EARLY
647 Number *
649
653 DEAL_II_DEPRECATED_EARLY
654 const Number *
655 end_raw() const;
656
664 template <typename OtherNumber>
665 constexpr DEAL_II_HOST_DEVICE Tensor &
667
674 constexpr DEAL_II_HOST_DEVICE Tensor &
675 operator=(const Number &d) &;
676
682 constexpr DEAL_II_HOST_DEVICE Tensor &
683 operator=(const Number &d) && = delete;
684
685#ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
691
697#endif
698
702 template <typename OtherNumber>
703 constexpr bool
705
709 template <typename OtherNumber>
710 constexpr bool
712
718 template <typename OtherNumber>
719 constexpr DEAL_II_HOST_DEVICE Tensor &
721
727 template <typename OtherNumber>
728 constexpr DEAL_II_HOST_DEVICE Tensor &
730
737 template <typename OtherNumber>
738 constexpr DEAL_II_HOST_DEVICE Tensor &
739 operator*=(const OtherNumber &factor);
740
746 template <typename OtherNumber>
747 constexpr DEAL_II_HOST_DEVICE Tensor &
748 operator/=(const OtherNumber &factor);
749
756 operator-() const;
757
770 constexpr void
772
782 norm() const;
783
790 constexpr DEAL_II_HOST_DEVICE
792 norm_square() const;
793
804 template <typename OtherNumber>
807
818 template <class Iterator>
819 void
820 unroll(const Iterator begin, const Iterator end) const;
821
826 static constexpr DEAL_II_HOST_DEVICE unsigned int
828
835 unrolled_to_component_indices(const unsigned int i);
836
841 static constexpr std::size_t
843
849 template <class Archive>
850 void
851 serialize(Archive &ar, const unsigned int version);
852
858
859private:
865 std::conditional_t<rank_ == 1,
866 std::array<Number, dim>,
867 std::array<Tensor<rank_ - 1, dim, Number>, dim>>
869
876 template <typename ArrayLike, std::size_t... Indices>
877 constexpr DEAL_II_HOST_DEVICE
878 Tensor(const ArrayLike &initializer, std::index_sequence<Indices...>);
879
880 // Allow an arbitrary Tensor to access the underlying values.
881 template <int, int, typename>
882 friend class Tensor;
883
884 // Point is allowed access to the coordinates. This is supposed to improve
885 // speed.
886 friend class Point<dim, Number>;
887};
888
889
890#ifndef DOXYGEN
891namespace internal
892{
893 // Workaround: The following 4 overloads are necessary to be able to
894 // compile the library with Apple Clang 8 and older. We should remove
895 // these overloads again when we bump the minimal required version to
896 // something later than clang-3.6 / Apple Clang 6.3.
897 template <int rank, int dim, typename T, typename U>
898 struct ProductTypeImpl<Tensor<rank, dim, T>, std::complex<U>>
899 {
900 using type =
902 };
903
904 template <int rank, int dim, typename T, typename U>
905 struct ProductTypeImpl<Tensor<rank, dim, std::complex<T>>, std::complex<U>>
906 {
907 using type =
909 };
910
911 template <typename T, int rank, int dim, typename U>
912 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, U>>
913 {
914 using type =
916 };
917
918 template <int rank, int dim, typename T, typename U>
919 struct ProductTypeImpl<std::complex<T>, Tensor<rank, dim, std::complex<U>>>
920 {
921 using type =
923 };
924 // end workaround
925
930 template <int rank, int dim, typename T>
931 struct NumberType<Tensor<rank, dim, T>>
932 {
933 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const
936 {
937 return t;
938 }
939
941 value(const T &t)
942 {
944 tmp = t;
945 return tmp;
946 }
947 };
948} // namespace internal
949
950
951/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
952
953
954template <int dim, typename Number>
957 // Some auto-differentiable numbers need explicit
958 // zero initialization such as adtl::adouble.
959 : Tensor{0.0}
960{}
961
962
963
964template <int dim, typename Number>
965template <typename OtherNumber>
967Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
968 : value(internal::NumberType<Number>::value(initializer))
969{}
970
971
972
973template <int dim, typename Number>
974template <typename OtherNumber>
977 : Tensor{p.value}
978{}
979
980
981# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
982template <int dim, typename Number>
985 : value{other.value}
986{}
987
988
989
990template <int dim, typename Number>
993 : value{std::move(other.value)}
994{}
995# endif
996
997
998
999template <int dim, typename Number>
1002{
1003 Assert(dim != 0,
1004 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1005 return value;
1006}
1007
1008
1009template <int dim, typename Number>
1010constexpr inline DEAL_II_ALWAYS_INLINE
1012{
1013 Assert(dim != 0,
1014 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1015 return value;
1016}
1017
1018
1019
1020template <int dim, typename Number>
1021template <typename OtherNumber>
1024{
1026 return *this;
1027}
1028
1029
1030# if defined(__INTEL_COMPILER) || defined(DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG)
1031template <int dim, typename Number>
1034{
1035 value = p.value;
1036 return *this;
1037}
1038# endif
1039
1040# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1041template <int dim, typename Number>
1044{
1045 value = std::move(other.value);
1046 return *this;
1047}
1048# endif
1049
1050
1051
1052template <int dim, typename Number>
1053template <typename OtherNumber>
1055Tensor<0, dim, Number>::operator=(const OtherNumber &d) &
1056{
1058 return *this;
1059}
1060
1061
1062template <int dim, typename Number>
1063template <typename OtherNumber>
1064constexpr inline bool
1066{
1067# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1068 Assert(!(std::is_same_v<Number, adouble> ||
1069 std::is_same_v<OtherNumber, adouble>),
1070 ExcMessage(
1071 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1072 "been extended to support advanced branching."));
1073# endif
1074
1075 return numbers::values_are_equal(value, p.value);
1076}
1077
1078
1079template <int dim, typename Number>
1080template <typename OtherNumber>
1081constexpr bool
1083{
1084 return !((*this) == p);
1085}
1086
1087
1088template <int dim, typename Number>
1089template <typename OtherNumber>
1092{
1093 value += p.value;
1094 return *this;
1095}
1096
1097
1098template <int dim, typename Number>
1099template <typename OtherNumber>
1102{
1103 value -= p.value;
1104 return *this;
1105}
1106
1107
1108
1109namespace internal
1110{
1111 namespace ComplexWorkaround
1112 {
1113 template <typename Number, typename OtherNumber>
1115 multiply_assign_scalar(Number &val, const OtherNumber &s)
1116 {
1117 val *= s;
1118 }
1119
1120 template <typename Number, typename OtherNumber>
1122 multiply_assign_scalar(std::complex<Number> &val, const OtherNumber &s)
1123 {
1124# if KOKKOS_VERSION >= 30600
1125 KOKKOS_IF_ON_HOST((val *= s;))
1126 KOKKOS_IF_ON_DEVICE(({
1127 (void)val;
1128 (void)s;
1129 Kokkos::abort(
1130 "This function is not implemented for std::complex<Number>!\n");
1131 }))
1132# else
1133# ifdef KOKKOS_ACTIVE_EXECUTION_MEMORY_SPACE_HOST
1134 val *= s;
1135# else
1136 (void)val;
1137 (void)s;
1138 Kokkos::abort(
1139 "This function is not implemented for std::complex<Number>!\n");
1140# endif
1141# endif
1142 }
1143 } // namespace ComplexWorkaround
1144} // namespace internal
1145
1146
1147template <int dim, typename Number>
1148template <typename OtherNumber>
1150Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
1151{
1152 internal::ComplexWorkaround::multiply_assign_scalar(value, s);
1153 return *this;
1154}
1155
1156
1157
1158template <int dim, typename Number>
1159template <typename OtherNumber>
1161Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
1162{
1163 value /= s;
1164 return *this;
1165}
1166
1167
1168template <int dim, typename Number>
1171{
1172 return -value;
1173}
1174
1175
1176template <int dim, typename Number>
1179{
1180 Assert(dim != 0,
1181 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1183}
1184
1185
1186template <int dim, typename Number>
1190{
1191 Assert(dim != 0,
1192 ExcMessage("Cannot access an object of type Tensor<0,0,Number>"));
1194}
1195
1196
1197
1198template <int dim, typename Number>
1199constexpr inline void
1201{
1202 // Some auto-differentiable numbers need explicit
1203 // zero initialization.
1205}
1206
1207
1208
1209template <int dim, typename Number>
1210template <class Iterator>
1211inline void
1212Tensor<0, dim, Number>::unroll(const Iterator begin, const Iterator end) const
1213{
1214 (void)end;
1215 AssertDimension(std::distance(begin, end), n_independent_components);
1216 Assert(dim != 0,
1217 ExcMessage("Cannot unroll an object of type Tensor<0,0,Number>"));
1218 Assert(std::distance(begin, end) >= 1,
1219 ExcMessage("The provided iterator range must contain at least one "
1220 "element."));
1221 *begin = value;
1222}
1223
1224
1225
1226template <int dim, typename Number>
1227template <class Archive>
1228inline void
1229Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
1230{
1231 ar &value;
1232}
1233
1234
1235template <int dim, typename Number>
1237
1238
1239/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
1240
1241template <int rank_, int dim, typename Number>
1242template <typename ArrayLike, std::size_t... indices>
1244Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
1245 std::index_sequence<indices...>)
1246 // Extract from the 'initializer' a sequence of elements via template
1247 // pack evaluation. This could be as easy as
1248 // values{{ (initializer[indices])... }}
1249 // but of course in practice it is not. The challenge is that if rank>1,
1250 // we want to pass the elements initializer[indices] down to the next
1251 // lower rank tensor for evaluation unchanged. But at the rank==1 level,
1252 // we need to convert to the scalar type 'Number'. This would all be
1253 // relatively straightforward if we could rely on automatic type
1254 // conversion, but for some autodifferentiation types, the conversion
1255 // from the AD to double (i.e., the extraction of a scalar value) is
1256 // not implicit, and we need to call internal::NumberType<Number>::value() --
1257 // but as mentioned, we can only do that for rank==1.
1258 //
1259 // We can achieve all of this by dispatching to a lambda function within
1260 // which we can use a 'if constexpr'.
1261 : values{{([&initializer]() -> value_type {
1262 if constexpr (rank_ == 1)
1263 return internal::NumberType<Number>::value(initializer[indices]);
1264 else
1265 return value_type(initializer[indices]);
1266 }())...}}
1267{
1268 static_assert(sizeof...(indices) == dim,
1269 "dim should match the number of indices");
1270}
1271
1272
1273# ifdef DEAL_II_HAVE_CXX20
1274
1275template <int rank_, int dim, typename Number>
1278 : values(
1279 // In order to initialize the std::array<Number,dim>, we would need a
1280 // brace-enclosed list of length 'dim'. There is no way in C++ to create
1281 // such a list in-place, but we can come up with a lambda function that
1282 // expands such a list via template-pack expansion, and then uses this
1283 // list to initialize a std::array which it then returns.
1284 //
1285 // The trick to come up with such a lambda function is to have a function
1286 // that takes an argument that depends on a template-pack of integers.
1287 // We will call the function with an integer list of length 'dim', and
1288 // in the function itself expand that pack in a way that it serves as
1289 // a brace-enclosed list of initializers for a std::array.
1290 //
1291 // Of course, we do not want to initialize the array with the integers,
1292 // but with zeros. (Or, more correctly, a zero of the element type.)
1293 // The canonical way to do this would be using the comma operator:
1294 // (sequence_element, 0.0)
1295 // returns zero, and
1296 // (sequence, 0.0)...
1297 // returns a list of zeros of the right length. Unfortunately, some
1298 // compilers then warn that the left side of the comma expression has
1299 // no effect -- well, bummer, that was of course exactly the idea.
1300 // We could work around this by using
1301 // (sequence_element * 0.0)
1302 // instead, assuming that the compiler will optimize (known) integer
1303 // times zero to zero, and similarly for (known) integer times times
1304 // default-initialized tensor.
1305 //
1306 // But, instead of relying on compiler optimizations, a better way is
1307 // to simply have another (nested) lambda function that takes the
1308 // integer sequence element as an argument and ignores it, just
1309 // returning a zero instead.
1310 []<std::size_t... I>(
1311 const std::index_sequence<I...> &) constexpr -> decltype(values) {
1312 if constexpr (rank_ == 1)
1313 {
1314 auto get_zero_and_ignore_argument = [](int) {
1316 };
1317 return {{(get_zero_and_ignore_argument(I))...}};
1318 }
1319 else
1320 {
1321 auto get_zero_and_ignore_argument = [](int) {
1322 return Tensor<rank_ - 1, dim, Number>();
1323 };
1324 return {{(get_zero_and_ignore_argument(I))...}};
1325 }
1326 }(std::make_index_sequence<dim>()))
1327{}
1328
1329# else
1330
1331// The C++17 case works in essence the same, except that we can't use a
1332// lambda function with explicit template parameters, i.e., we can't do
1333// []<std::size_t... I>(const std::index_sequence<I...> &)
1334// as above because that's a C++20 feature. Lambda functions in C++17 can
1335// have template packs as arguments, but we need the ability to *name*
1336// that template pack (the 'I' above) and that's not possible in C++17.
1337//
1338// We work around this by moving the lambda function to a global function
1339// and using the traditional template syntax on it.
1340namespace internal
1341{
1342 namespace TensorInitialization
1343 {
1344 template <int rank, int dim, typename Number, std::size_t... I>
1345 constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
1346 make_zero_array(const std::index_sequence<I...> &)
1347 {
1348 static_assert(sizeof...(I) == dim, "This is bad.");
1349
1350 // First peel off the case dim==0. If we don't, some compilers
1351 // will warn below that we define these lambda functions but
1352 // never use them (because the expanded list has zero elements,
1353 // and the get_zero_and_ignore_argument() function is not used...)
1354 if constexpr (dim == 0)
1355 {
1356 return {};
1357 }
1358 else if constexpr (rank == 1)
1359 {
1360 auto get_zero_and_ignore_argument = [](int) {
1362 };
1363 return {{(get_zero_and_ignore_argument(I))...}};
1364 }
1365 else
1366 {
1367 auto get_zero_and_ignore_argument = [](int) {
1368 return Tensor<rank - 1, dim, Number>();
1369 };
1370 return {{(get_zero_and_ignore_argument(I))...}};
1371 }
1372 }
1373 } // namespace TensorInitialization
1374} // namespace internal
1375
1376
1377template <int rank_, int dim, typename Number>
1380 : values(internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
1381 std::make_index_sequence<dim>()))
1382{}
1383
1384
1385# endif
1386
1387
1388template <int rank_, int dim, typename Number>
1390Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
1391 : Tensor(initializer, std::make_index_sequence<dim>{})
1392{}
1393
1394
1395
1396template <int rank_, int dim, typename Number>
1397template <typename ElementType, typename MemorySpace>
1400 const ArrayView<ElementType, MemorySpace> &initializer)
1401{
1402 // make nvcc happy
1403 const int my_n_independent_components = n_independent_components;
1404 AssertDimension(initializer.size(), my_n_independent_components);
1405
1406 for (unsigned int i = 0; i < my_n_independent_components; ++i)
1407 (*this)[unrolled_to_component_indices(i)] = initializer[i];
1408}
1409
1410
1411
1412template <int rank_, int dim, typename Number>
1413template <typename OtherNumber>
1416 const Tensor<rank_, dim, OtherNumber> &initializer)
1417 : Tensor(initializer, std::make_index_sequence<dim>{})
1418{}
1419
1420
1421
1422template <int rank_, int dim, typename Number>
1423template <typename OtherNumber>
1424constexpr DEAL_II_ALWAYS_INLINE
1426 const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
1427 : Tensor(initializer, std::make_index_sequence<dim>{})
1428{}
1429
1430
1431
1432template <int rank_, int dim, typename Number>
1433template <typename OtherNumber>
1435operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
1436{
1437 Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> x;
1438 std::copy(values.begin(), values.end(), x.values.begin());
1439 return x;
1440}
1441
1442
1443# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1444template <int rank_, int dim, typename Number>
1445constexpr DEAL_II_ALWAYS_INLINE
1447 : values(other.values)
1448{}
1449
1450
1451
1452template <int rank_, int dim, typename Number>
1453constexpr DEAL_II_ALWAYS_INLINE
1455 : values(std::move(other.values))
1456{}
1457# endif
1458
1459
1460
1461template <int rank_, int dim, typename Number>
1464 Tensor<rank_, dim, Number>::operator[](const unsigned int i)
1465{
1466 Assert(dim != 0,
1467 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1468 AssertIndexRange(i, dim);
1469 DEAL_II_CXX23_ASSUME(i < dim);
1470
1471 return values[i];
1472}
1473
1474
1475template <int rank_, int dim, typename Number>
1476constexpr DEAL_II_ALWAYS_INLINE
1478 Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
1479{
1480 Assert(dim != 0,
1481 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1482 AssertIndexRange(i, dim);
1483 DEAL_II_CXX23_ASSUME(i < dim);
1484
1485 return values[i];
1486}
1487
1488
1489template <int rank_, int dim, typename Number>
1490constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1492{
1493 Assert(dim != 0,
1494 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1495
1496 return TensorAccessors::extract<rank_>(*this, indices);
1497}
1498
1499
1500
1501template <int rank_, int dim, typename Number>
1502constexpr inline DEAL_II_ALWAYS_INLINE Number &
1504{
1505 Assert(dim != 0,
1506 ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
1507
1508 return TensorAccessors::extract<rank_>(*this, indices);
1509}
1510
1511
1512
1513template <int rank_, int dim, typename Number>
1514inline Number *
1516{
1517 static_assert(rank_ == 1,
1518 "This function is only available for rank-1 tensors "
1519 "because higher-rank tensors may not store their elements "
1520 "in a contiguous array.");
1521
1522 return std::addressof(
1523 this->operator[](this->unrolled_to_component_indices(0)));
1524}
1525
1526
1527
1528template <int rank_, int dim, typename Number>
1529inline const Number *
1531{
1532 static_assert(rank_ == 1,
1533 "This function is only available for rank-1 tensors "
1534 "because higher-rank tensors may not store their elements "
1535 "in a contiguous array.");
1536
1537 return std::addressof(
1538 this->operator[](this->unrolled_to_component_indices(0)));
1539}
1540
1541
1542
1543template <int rank_, int dim, typename Number>
1544inline Number *
1546{
1547 static_assert(rank_ == 1,
1548 "This function is only available for rank-1 tensors "
1549 "because higher-rank tensors may not store their elements "
1550 "in a contiguous array.");
1551
1552 return begin_raw() + n_independent_components;
1553}
1554
1555
1556
1557template <int rank_, int dim, typename Number>
1558inline const Number *
1560{
1561 static_assert(rank_ == 1,
1562 "This function is only available for rank-1 tensors "
1563 "because higher-rank tensors may not store their elements "
1564 "in a contiguous array.");
1565
1566 return begin_raw() + n_independent_components;
1567}
1568
1569
1570
1571template <int rank_, int dim, typename Number>
1572template <typename OtherNumber>
1575{
1576 // The following loop could be written more concisely using std::copy, but
1577 // that function is only constexpr from C++20 on.
1578 for (unsigned int i = 0; i < dim; ++i)
1579 values[i] = t.values[i];
1580 return *this;
1581}
1582
1583
1584
1585template <int rank_, int dim, typename Number>
1588 Tensor<rank_, dim, Number>::operator=(const Number &d) &
1589{
1591 (void)d;
1592
1593 for (unsigned int i = 0; i < dim; ++i)
1594 values[i] = internal::NumberType<Number>::value(0.0);
1595 return *this;
1596}
1597
1598
1599# ifdef DEAL_II_DELETED_MOVE_CONSTRUCTOR_BUG
1600template <int rank_, int dim, typename Number>
1603{
1604 for (unsigned int i = 0; i < dim; ++i)
1605 values[i] = other.values[i];
1606 return *this;
1607}
1608
1609
1610
1611template <int rank_, int dim, typename Number>
1614 Tensor<rank_, dim, Number> &&other) noexcept
1615{
1616 for (unsigned int i = 0; i < dim; ++i)
1617 values[i] = other.values[i];
1618 return *this;
1619}
1620# endif
1621
1622
1623template <int rank_, int dim, typename Number>
1624template <typename OtherNumber>
1625constexpr inline bool
1627 const Tensor<rank_, dim, OtherNumber> &p) const
1628{
1629# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
1630 Assert(!(std::is_same_v<Number, adouble> ||
1631 std::is_same_v<OtherNumber, adouble>),
1632 ExcMessage(
1633 "The Tensor equality operator for ADOL-C taped numbers has not yet "
1634 "been extended to support advanced branching."));
1635# endif
1636
1637 for (unsigned int i = 0; i < dim; ++i)
1638 if (numbers::values_are_not_equal(values[i], p.values[i]))
1639 return false;
1640 return true;
1641}
1642
1643
1644// At some places in the library, we have Point<0> for formal reasons
1645// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
1646// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
1647// in the above function that the loop end check always fails, we
1648// implement this function here
1649template <>
1650template <>
1651constexpr inline bool
1653{
1654 return true;
1655}
1656
1657
1658template <int rank_, int dim, typename Number>
1659template <typename OtherNumber>
1660constexpr bool
1662 const Tensor<rank_, dim, OtherNumber> &p) const
1663{
1664 return !((*this) == p);
1665}
1666
1667
1668template <int rank_, int dim, typename Number>
1669template <typename OtherNumber>
1670constexpr inline DEAL_II_ALWAYS_INLINE
1674{
1675 for (unsigned int i = 0; i < dim; ++i)
1676 values[i] += p.values[i];
1677 return *this;
1678}
1679
1680
1681template <int rank_, int dim, typename Number>
1682template <typename OtherNumber>
1683constexpr inline DEAL_II_ALWAYS_INLINE
1687{
1688 for (unsigned int i = 0; i < dim; ++i)
1689 values[i] -= p.values[i];
1690 return *this;
1691}
1692
1693
1694template <int rank_, int dim, typename Number>
1695template <typename OtherNumber>
1696constexpr inline DEAL_II_ALWAYS_INLINE
1698 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
1699{
1700 for (unsigned int i = 0; i < dim; ++i)
1701 values[i] *= s;
1702 return *this;
1703}
1704
1705
1706
1707template <int rank_, int dim, typename Number>
1708template <typename OtherNumber>
1709constexpr inline DEAL_II_ALWAYS_INLINE
1711 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
1712{
1713 if constexpr (std::is_integral<
1714 typename ProductType<Number, OtherNumber>::type>::value ||
1715 std::is_same_v<Number, Differentiation::SD::Expression>)
1716 {
1717 // recurse over the base objects
1718 for (unsigned int d = 0; d < dim; ++d)
1719 values[d] /= s;
1720 }
1721 else
1722 {
1723 // If we can, avoid division by multiplying by the inverse of the given
1724 // factor:
1725 const Number inverse_factor = Number(1.) / s;
1726 for (unsigned int d = 0; d < dim; ++d)
1727 values[d] *= inverse_factor;
1728 }
1729
1730 return *this;
1731}
1732
1733
1734template <int rank_, int dim, typename Number>
1735constexpr inline DEAL_II_ALWAYS_INLINE
1738{
1740
1741 for (unsigned int i = 0; i < dim; ++i)
1742 tmp.values[i] = -values[i];
1743
1744 return tmp;
1745}
1746
1747
1748template <int rank_, int dim, typename Number>
1751{
1752 // Handle cases of a tensor consisting of just one number more
1753 // efficiently:
1754 if constexpr ((rank_ == 1) && (dim == 1) && std::is_arithmetic_v<Number>)
1755 {
1756 return std::abs(values[0]);
1757 }
1758 else if constexpr ((rank_ == 2) && (dim == 1) && std::is_arithmetic_v<Number>)
1759 {
1760 return std::abs(values[0][0]);
1761 }
1762 else
1763 {
1764 // Otherwise fall back to the naive algorithm of taking the square root of
1765 // the sum of squares.
1766
1767 // Make things work with AD types by letting the compiler look up
1768 // the symbol sqrt in namespace std and in the type-associated
1769 // namespaces
1770 using std::sqrt;
1771 return sqrt(norm_square());
1772 }
1773}
1774
1775
1776template <int rank_, int dim, typename Number>
1780{
1781 if constexpr (dim == 0)
1782 return internal::NumberType<
1783 typename numbers::NumberTraits<Number>::real_type>::value(0.0);
1784 else if constexpr (rank_ == 1)
1785 {
1786 // For rank-1 tensors, the square of the norm is simply the sum of
1787 // squares of the elements:
1790 for (unsigned int i = 1; i < dim; ++i)
1792
1793 return s;
1794 }
1795 else
1796 {
1797 // For higher-rank tensors, the square of the norm is the sum
1798 // of squares of sub-tensors
1800 values[0].norm_square();
1801 for (unsigned int i = 1; i < dim; ++i)
1802 s += values[i].norm_square();
1803
1804 return s;
1805 }
1806}
1807
1808
1809
1810template <int rank_, int dim, typename Number>
1811template <typename OtherNumber>
1812inline void
1814{
1815 unroll(result.begin(), result.end());
1816}
1817
1818
1819
1820template <int rank_, int dim, typename Number>
1821template <class Iterator>
1822inline void
1823Tensor<rank_, dim, Number>::unroll(const Iterator begin,
1824 const Iterator end) const
1825{
1826 if constexpr (rank_ > 1)
1827 {
1828 // For higher-rank tensors, we recurse to the sub-tensors:
1829 Iterator next = begin;
1830 for (unsigned int i = 0; i < dim; ++i)
1831 {
1832 values[i].unroll(next, end);
1833 std::advance(
1835 }
1836 }
1837 else
1838 {
1839 // For rank-1 tensors, we can simply copy the current elements from
1840 // our linear array into the output range, and then return the
1841 // iterator moved forward by 'dim' elements:
1842 (void)end;
1843 Assert(std::distance(begin, end) >= dim,
1844 ExcMessage(
1845 "The provided iterator range must contain at least 'dim' "
1846 "elements."));
1847 std::copy(std::begin(values), std::end(values), begin);
1848 }
1849}
1850
1851
1852
1853template <int rank_, int dim, typename Number>
1854constexpr inline unsigned int
1856 const TableIndices<rank_> &indices)
1857{
1858 unsigned int index = 0;
1859 for (int r = 0; r < rank_; ++r)
1860 index = index * dim + indices[r];
1861
1862 return index;
1863}
1864
1865
1866
1867template <int rank_, int dim, typename Number>
1868constexpr inline TableIndices<rank_>
1870{
1871 // Work-around nvcc warning
1872 unsigned int dummy = n_independent_components;
1873 AssertIndexRange(i, dummy);
1874 (void)dummy;
1875
1876 if constexpr (dim == 0)
1877 {
1878 Assert(false,
1879 ExcMessage(
1880 "A tensor with dimension 0 does not store any elements. "
1881 "There is no indexing that can address its elements."));
1882 return {};
1883 }
1884 else
1885 {
1886 TableIndices<rank_> indices;
1887
1888 unsigned int remainder = i;
1889 for (int r = rank_ - 1; r >= 0; --r)
1890 {
1891 indices[r] = remainder % dim;
1892 remainder = remainder / dim;
1893 }
1894 Assert(remainder == 0, ExcInternalError());
1895
1896 return indices;
1897 }
1898}
1899
1900
1901template <int rank_, int dim, typename Number>
1902constexpr inline void
1904{
1905 for (unsigned int i = 0; i < dim; ++i)
1906 values[i] = internal::NumberType<Number>::value(0.0);
1907}
1908
1909
1910template <int rank_, int dim, typename Number>
1911constexpr std::size_t
1913{
1914 return sizeof(Tensor<rank_, dim, Number>);
1915}
1916
1917
1918template <int rank_, int dim, typename Number>
1919template <class Archive>
1920inline void
1921Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1922{
1923 if constexpr (rank_ > 1)
1924 ar &values;
1925 else
1926 ar &boost::serialization::make_array(&values[0], dim);
1927}
1928
1929
1930template <int rank_, int dim, typename Number>
1932
1933#endif // DOXYGEN
1934
1935/* ----------------- Non-member functions operating on tensors. ------------ */
1936
1949template <int rank_, int dim, typename Number>
1950inline std::ostream &
1951operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1952{
1953 for (unsigned int i = 0; i < dim; ++i)
1954 {
1955 out << p[i];
1956 if (i != dim - 1)
1957 for (unsigned int j = 0; j < rank_; ++j)
1958 out << ' ';
1959 }
1960
1961 return out;
1962}
1963
1964
1971template <int dim, typename Number>
1972inline std::ostream &
1973operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1974{
1975 out << static_cast<const Number &>(p);
1976 return out;
1977}
1978
1979
1998template <int dim, typename Number, typename Other>
2001 operator*(const Other &object, const Tensor<0, dim, Number> &t)
2002{
2003 return object * static_cast<const Number &>(t);
2004}
2005
2006
2007
2018template <int dim, typename Number, typename Other>
2021 operator*(const Tensor<0, dim, Number> &t, const Other &object)
2022{
2023 return static_cast<const Number &>(t) * object;
2024}
2025
2026
2038template <int dim, typename Number, typename OtherNumber>
2042 const Tensor<0, dim, OtherNumber> &src2)
2043{
2044 return static_cast<const Number &>(src1) *
2045 static_cast<const OtherNumber &>(src2);
2046}
2047
2048
2056template <int dim, typename Number, typename OtherNumber>
2058 Tensor<0,
2059 dim,
2060 typename ProductType<Number,
2061 typename EnableIfScalar<OtherNumber>::type>::type>
2062 operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2063{
2064 return static_cast<const Number &>(t) / factor;
2065}
2066
2067
2075template <int dim, typename Number, typename OtherNumber>
2080{
2081 return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2082}
2083
2084
2092template <int dim, typename Number, typename OtherNumber>
2097{
2098 return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2099}
2100
2101
2114template <int rank, int dim, typename Number, typename OtherNumber>
2116 Tensor<rank,
2117 dim,
2118 typename ProductType<Number,
2119 typename EnableIfScalar<OtherNumber>::type>::type>
2120 operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2121{
2123 tt *= factor;
2124 return tt;
2125}
2126
2127
2140template <int rank, int dim, typename Number, typename OtherNumber>
2142 Tensor<rank,
2143 dim,
2145 OtherNumber>::type>
2146 operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2147{
2148 // simply forward to the operator above
2149 return t * factor;
2150}
2151
2152
2153
2163template <int rank, int dim, typename Number, typename OtherNumber>
2165 Tensor<rank,
2166 dim,
2167 typename ProductType<Number,
2168 typename EnableIfScalar<OtherNumber>::type>::type>
2169 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2170{
2172 tt /= factor;
2173 return tt;
2174}
2175
2176
2186template <int rank, int dim, typename Number, typename OtherNumber>
2196
2197
2207template <int rank, int dim, typename Number, typename OtherNumber>
2217
2224template <int dim, typename Number, typename OtherNumber>
2225inline constexpr DEAL_II_ALWAYS_INLINE
2228 const Tensor<0, dim, OtherNumber> &src2)
2229{
2231
2232 tmp *= src2;
2233
2234 return tmp;
2235}
2236
2253template <int rank, int dim, typename Number, typename OtherNumber>
2254inline constexpr DEAL_II_ALWAYS_INLINE
2258{
2260
2261 for (unsigned int i = 0; i < dim; ++i)
2264
2265 return tmp;
2266}
2267
2312template <int rank_1,
2313 int rank_2,
2314 int dim,
2315 typename Number,
2316 typename OtherNumber,
2317 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2318constexpr inline DEAL_II_ALWAYS_INLINE
2319 typename Tensor<rank_1 + rank_2 - 2,
2320 dim,
2321 typename ProductType<Number, OtherNumber>::type>::tensor_type
2324{
2325 // Treat some common cases separately. Specifically, these are the dot
2326 // product between two rank-1 tensors, and the product between a
2327 // rank-2 tensor and a rank-1 tensor. Both of these lead to a linear
2328 // loop over adjacent memory and can be dealt with efficiently; in the
2329 // latter case (rank-2 times rank-1), we implement things by deferring
2330 // to rank-1 times rank-1 dot products.
2331 if constexpr ((rank_1 == 1) && (rank_2 == 1))
2332 {
2333 // This is a dot product between two rank-1 tensors. Write it out as
2334 // a linear loop:
2335 static_assert(dim > 0, "Tensors cannot have dimension zero.");
2336 typename ProductType<Number, OtherNumber>::type sum = src1[0] * src2[0];
2337 for (unsigned int i = 1; i < dim; ++i)
2338 sum += src1[i] * src2[i];
2339
2340 return sum;
2341 }
2342 else if constexpr ((rank_1 == 2) && (rank_2 == 1))
2343 {
2344 // This is a product between a rank-2 and a rank-1 tensor. This
2345 // corresponds to taking dot products between the rows of the former
2346 // and the latter.
2347 typename Tensor<
2348 rank_1 + rank_2 - 2,
2349 dim,
2350 typename ProductType<Number, OtherNumber>::type>::tensor_type result;
2351 for (unsigned int i = 0; i < dim; ++i)
2352 result[i] += src1[i] * src2;
2353
2354 return result;
2355 }
2356 else
2357 {
2358 // Treat all of the other cases using the more general contraction
2359 // machinery.
2360 typename Tensor<
2361 rank_1 + rank_2 - 2,
2362 dim,
2363 typename ProductType<Number, OtherNumber>::type>::tensor_type result{};
2364
2365 TensorAccessors::internal::
2366 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2367 reordered = TensorAccessors::reordered_index_view<0, rank_2>(src2);
2368 TensorAccessors::contract<1, rank_1, rank_2, dim>(result,
2369 src1,
2370 reordered);
2371
2372 return result;
2373 }
2374}
2375
2376
2405template <int index_1,
2406 int index_2,
2407 int rank_1,
2408 int rank_2,
2409 int dim,
2410 typename Number,
2411 typename OtherNumber>
2412constexpr inline DEAL_II_ALWAYS_INLINE
2413 typename Tensor<rank_1 + rank_2 - 2,
2414 dim,
2415 typename ProductType<Number, OtherNumber>::type>::tensor_type
2418{
2419 Assert(0 <= index_1 && index_1 < rank_1,
2420 ExcMessage(
2421 "The specified index_1 must lie within the range [0,rank_1)"));
2422 Assert(0 <= index_2 && index_2 < rank_2,
2423 ExcMessage(
2424 "The specified index_2 must lie within the range [0,rank_2)"));
2425
2426 using namespace TensorAccessors;
2427 using namespace TensorAccessors::internal;
2428
2429 // Reorder index_1 to the end of src1:
2431 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2432
2433 // Reorder index_2 to the end of src2:
2434 const ReorderedIndexView<index_2,
2435 rank_2,
2437 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2438
2439 typename Tensor<rank_1 + rank_2 - 2,
2440 dim,
2441 typename ProductType<Number, OtherNumber>::type>::tensor_type
2442 result{};
2443 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2444 return result;
2445}
2446
2447
2478template <int index_1,
2479 int index_2,
2480 int index_3,
2481 int index_4,
2482 int rank_1,
2483 int rank_2,
2484 int dim,
2485 typename Number,
2486 typename OtherNumber>
2487constexpr inline
2488 typename Tensor<rank_1 + rank_2 - 4,
2489 dim,
2490 typename ProductType<Number, OtherNumber>::type>::tensor_type
2491 double_contract(const Tensor<rank_1, dim, Number> &src1,
2493{
2494 Assert(0 <= index_1 && index_1 < rank_1,
2495 ExcMessage(
2496 "The specified index_1 must lie within the range [0,rank_1)"));
2497 Assert(0 <= index_3 && index_3 < rank_1,
2498 ExcMessage(
2499 "The specified index_3 must lie within the range [0,rank_1)"));
2500 Assert(index_1 != index_3,
2501 ExcMessage("index_1 and index_3 must not be the same"));
2502 Assert(0 <= index_2 && index_2 < rank_2,
2503 ExcMessage(
2504 "The specified index_2 must lie within the range [0,rank_2)"));
2505 Assert(0 <= index_4 && index_4 < rank_2,
2506 ExcMessage(
2507 "The specified index_4 must lie within the range [0,rank_2)"));
2508 Assert(index_2 != index_4,
2509 ExcMessage("index_2 and index_4 must not be the same"));
2510
2511 using namespace TensorAccessors;
2512 using namespace TensorAccessors::internal;
2513
2514 // Reorder index_1 to the end of src1:
2516 reord_1 = TensorAccessors::reordered_index_view<index_1, rank_1>(src1);
2517
2518 // Reorder index_2 to the end of src2:
2520 reord_2 = TensorAccessors::reordered_index_view<index_2, rank_2>(src2);
2521
2522 // Now, reorder index_3 to the end of src1. We have to make sure to
2523 // preserve the original ordering: index_1 has been removed. If
2524 // index_3 > index_1, we have to use (index_3 - 1) instead:
2526 (index_3 < index_1 ? index_3 : index_3 - 1),
2527 rank_1,
2528 ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2529 reord_3 =
2530 TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2531 index_3 - 1,
2532 rank_1 > (reord_1);
2533
2534 // Now, reorder index_4 to the end of src2. We have to make sure to
2535 // preserve the original ordering: index_2 has been removed. If
2536 // index_4 > index_2, we have to use (index_4 - 1) instead:
2538 (index_4 < index_2 ? index_4 : index_4 - 1),
2539 rank_2,
2541 reord_4 =
2542 TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2543 index_4 - 1,
2544 rank_2 > (reord_2);
2545
2546 typename Tensor<rank_1 + rank_2 - 4,
2547 dim,
2548 typename ProductType<Number, OtherNumber>::type>::tensor_type
2549 result{};
2550 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2551 return result;
2552}
2553
2554
2567template <int rank, int dim, typename Number, typename OtherNumber>
2568constexpr inline DEAL_II_ALWAYS_INLINE
2570 scalar_product(const Tensor<rank, dim, Number> &left,
2571 const Tensor<rank, dim, OtherNumber> &right)
2572{
2574 TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
2575 return result;
2576}
2577
2578
2596template <template <int, int, typename> class TensorT1,
2597 template <int, int, typename>
2598 class TensorT2,
2599 template <int, int, typename>
2600 class TensorT3,
2601 int rank_1,
2602 int rank_2,
2603 int dim,
2604 typename T1,
2605 typename T2,
2606 typename T3>
2607constexpr inline DEAL_II_ALWAYS_INLINE
2609 contract3(const TensorT1<rank_1, dim, T1> &left,
2610 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2611 const TensorT3<rank_2, dim, T3> &right)
2612{
2613 using return_type =
2615 return TensorAccessors::contract3<rank_1, rank_2, dim, return_type>(left,
2616 middle,
2617 right);
2618}
2619
2620
2631template <int rank_1,
2632 int rank_2,
2633 int dim,
2634 typename Number,
2635 typename OtherNumber>
2636constexpr inline DEAL_II_ALWAYS_INLINE
2640{
2641 typename Tensor<rank_1 + rank_2,
2642 dim,
2643 typename ProductType<Number, OtherNumber>::type>::tensor_type
2644 result{};
2645 TensorAccessors::contract<0, rank_1, rank_2, dim>(result, src1, src2);
2646 return result;
2647}
2648
2649
2668template <int dim, typename Number>
2670cross_product_2d(const Tensor<1, dim, Number> &src)
2671{
2672 Assert(dim == 2, ExcInternalError());
2673
2675
2676 result[0] = src[1];
2677 result[1] = -src[0];
2678
2679 return result;
2680}
2681
2682
2692template <int dim, typename Number1, typename Number2>
2693constexpr inline DEAL_II_ALWAYS_INLINE
2695 cross_product_3d(const Tensor<1, dim, Number1> &src1,
2696 const Tensor<1, dim, Number2> &src2)
2697{
2698 Assert(dim == 3, ExcInternalError());
2699
2701
2702 // avoid compiler warnings
2703 constexpr int s0 = 0 % dim;
2704 constexpr int s1 = 1 % dim;
2705 constexpr int s2 = 2 % dim;
2706
2707 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2708 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2709 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2710
2711 return result;
2712}
2713
2714
2728template <int dim, typename Number>
2729constexpr inline DEAL_II_ALWAYS_INLINE Number
2731{
2732 // Compute the determinant using the Laplace expansion of the
2733 // determinant. We expand along the last row.
2734 Number det = internal::NumberType<Number>::value(0.0);
2735
2736 for (unsigned int k = 0; k < dim; ++k)
2737 {
2738 Tensor<2, dim - 1, Number> minor;
2739 for (unsigned int i = 0; i < dim - 1; ++i)
2740 for (unsigned int j = 0; j < dim - 1; ++j)
2741 minor[i][j] = t[i][j < k ? j : j + 1];
2742
2743 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2744
2745 det += t[dim - 1][k] * cofactor;
2746 }
2747
2748 return ((dim % 2 == 0) ? 1. : -1.) * det;
2749}
2750
2756template <typename Number>
2757constexpr DEAL_II_ALWAYS_INLINE Number
2759{
2760 return t[0][0];
2761}
2762
2768template <typename Number>
2769constexpr DEAL_II_ALWAYS_INLINE Number
2771{
2772 // hard-coded for efficiency reasons
2773 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2774}
2775
2781template <typename Number>
2782constexpr DEAL_II_ALWAYS_INLINE Number
2784{
2785 // hard-coded for efficiency reasons
2786 const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2787 internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2788 const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2789 internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2790 const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2791 internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2792 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2793}
2794
2795
2802template <int dim, typename Number>
2803constexpr inline DEAL_II_ALWAYS_INLINE Number
2805{
2806 Number t = d[0][0];
2807 for (unsigned int i = 1; i < dim; ++i)
2808 t += d[i][i];
2809 return t;
2810}
2811
2812
2821template <int dim, typename Number>
2822constexpr inline Tensor<2, dim, Number>
2824{
2825 Number return_tensor[dim][dim];
2826
2827 // if desired, take over the
2828 // inversion of a 4x4 tensor
2829 // from the FullMatrix
2831
2832 return Tensor<2, dim, Number>(return_tensor);
2833}
2834
2835
2836#ifndef DOXYGEN
2837
2838template <typename Number>
2841{
2842 Tensor<2, 1, Number> return_tensor;
2843
2844 return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2845
2846 return return_tensor;
2847}
2848
2849
2850template <typename Number>
2853{
2854 Tensor<2, 2, Number> return_tensor;
2855
2856 const Number inv_det_t = internal::NumberType<Number>::value(
2857 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2858 return_tensor[0][0] = t[1][1];
2859 return_tensor[0][1] = -t[0][1];
2860 return_tensor[1][0] = -t[1][0];
2861 return_tensor[1][1] = t[0][0];
2862 return_tensor *= inv_det_t;
2863
2864 return return_tensor;
2865}
2866
2867template <typename Number>
2870{
2871 Tensor<2, 3, Number> return_tensor;
2872
2873 const auto value = [](const auto &t) {
2875 };
2876
2877 return_tensor[0][0] = value(t[1][1] * t[2][2]) - value(t[1][2] * t[2][1]);
2878 return_tensor[0][1] = value(t[0][2] * t[2][1]) - value(t[0][1] * t[2][2]);
2879 return_tensor[0][2] = value(t[0][1] * t[1][2]) - value(t[0][2] * t[1][1]);
2880 return_tensor[1][0] = value(t[1][2] * t[2][0]) - value(t[1][0] * t[2][2]);
2881 return_tensor[1][1] = value(t[0][0] * t[2][2]) - value(t[0][2] * t[2][0]);
2882 return_tensor[1][2] = value(t[0][2] * t[1][0]) - value(t[0][0] * t[1][2]);
2883 return_tensor[2][0] = value(t[1][0] * t[2][1]) - value(t[1][1] * t[2][0]);
2884 return_tensor[2][1] = value(t[0][1] * t[2][0]) - value(t[0][0] * t[2][1]);
2885 return_tensor[2][2] = value(t[0][0] * t[1][1]) - value(t[0][1] * t[1][0]);
2886
2887 const Number inv_det_t =
2888 value(1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2889 t[0][2] * return_tensor[2][0]));
2890 return_tensor *= inv_det_t;
2891
2892 return return_tensor;
2893}
2894
2895#endif /* DOXYGEN */
2896
2897
2903template <int dim, typename Number>
2906{
2908 for (unsigned int i = 0; i < dim; ++i)
2909 {
2910 tt[i][i] = t[i][i];
2911 for (unsigned int j = i + 1; j < dim; ++j)
2912 {
2913 tt[i][j] = t[j][i];
2914 tt[j][i] = t[i][j];
2915 };
2916 }
2917 return tt;
2918}
2919
2920
2934template <int dim, typename Number>
2935constexpr Tensor<2, dim, Number>
2936adjugate(const Tensor<2, dim, Number> &t)
2937{
2938 return determinant(t) * invert(t);
2939}
2940
2941
2955template <int dim, typename Number>
2956constexpr Tensor<2, dim, Number>
2957cofactor(const Tensor<2, dim, Number> &t)
2958{
2959 return transpose(adjugate(t));
2960}
2961
2962
3026template <int dim, typename Number>
3029
3030
3038template <int dim, typename Number>
3039inline Number
3041{
3042 Number max = internal::NumberType<Number>::value(0.0);
3043 for (unsigned int j = 0; j < dim; ++j)
3044 {
3045 Number sum = internal::NumberType<Number>::value(0.0);
3046 for (unsigned int i = 0; i < dim; ++i)
3047 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3048
3049 if (sum > max)
3050 max = sum;
3051 }
3052
3053 return max;
3054}
3055
3056
3064template <int dim, typename Number>
3065inline Number
3067{
3068 Number max = internal::NumberType<Number>::value(0.0);
3069 for (unsigned int i = 0; i < dim; ++i)
3070 {
3071 Number sum = internal::NumberType<Number>::value(0.0);
3072 for (unsigned int j = 0; j < dim; ++j)
3073 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3074
3075 if (sum > max)
3076 max = sum;
3077 }
3078
3079 return max;
3080}
3081
3087#ifndef DOXYGEN
3088
3089
3090# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3091
3092// Specialization of functions for ADOL-C number types when
3093// the advanced branching feature is used
3094template <int dim>
3095inline adouble
3097{
3098 adouble max = internal::NumberType<adouble>::value(0.0);
3099 for (unsigned int j = 0; j < dim; ++j)
3100 {
3101 adouble sum = internal::NumberType<adouble>::value(0.0);
3102 for (unsigned int i = 0; i < dim; ++i)
3103 sum += fabs(t[i][j]);
3104
3105 condassign(max, (sum > max), sum, max);
3106 }
3107
3108 return max;
3109}
3110
3111
3112template <int dim>
3113inline adouble
3115{
3117 for (unsigned int i = 0; i < dim; ++i)
3118 {
3120 for (unsigned int j = 0; j < dim; ++j)
3121 sum += fabs(t[i][j]);
3122
3123 condassign(max, (sum > max), sum, max);
3124 }
3125
3126 return max;
3127}
3128
3129# endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3130
3131
3132#endif // DOXYGEN
3133
3135
3136#endif
std::size_t size() const
Definition array_view.h:684
Definition point.h:111
constexpr Tensor & operator*=(const OtherNumber &factor)
constexpr Tensor & operator=(const OtherNumber &d) &&=delete
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor & operator/=(const OtherNumber &factor)
constexpr Tensor & operator-=(const Tensor< 0, dim, OtherNumber > &rhs)
constexpr Tensor(const Tensor< 0, dim, OtherNumber > &initializer)
constexpr Tensor(const OtherNumber &initializer)
constexpr void clear()
constexpr real_type norm_square() const
constexpr bool operator!=(const Tensor< 0, dim, OtherNumber > &rhs) const
constexpr Tensor & operator=(const Tensor< 0, dim, OtherNumber > &rhs)
real_type norm() const
constexpr Tensor & operator+=(const Tensor< 0, dim, OtherNumber > &rhs)
void unroll(const Iterator begin, const Iterator end) const
constexpr bool operator==(const Tensor< 0, dim, OtherNumber > &rhs) const
typename numbers::NumberTraits< Number >::real_type real_type
Definition tensor.h: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:2169
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:2062
Number * begin_raw()
static constexpr unsigned int rank
Definition tensor.h:490
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:515
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:2041
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:531
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:2227
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:2095
friend class Tensor
Definition tensor.h:882
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:2189
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3066
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:2001
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3040
static constexpr unsigned int dimension
Definition tensor.h:485
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:2146
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:868
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:2078
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:2256
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:2210
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition tensor.h:2021
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:2120
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:506
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:990
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:1951
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:2227
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3066
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:2001
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3040