Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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:
1841 (void)end;
1842 Assert(std::distance(begin, end) >= dim,
1843 ExcMessage(
1844 "The provided iterator range must contain at least 'dim' "
1845 "elements."));
1846 std::copy(std::begin(values), std::end(values), begin);
1847 }
1848}
1849
1850
1851
1852template <int rank_, int dim, typename Number>
1853constexpr inline unsigned int
1855 const TableIndices<rank_> &indices)
1856{
1857 unsigned int index = 0;
1858 for (int r = 0; r < rank_; ++r)
1859 index = index * dim + indices[r];
1860
1861 return index;
1862}
1863
1864
1865
1866template <int rank_, int dim, typename Number>
1867constexpr inline TableIndices<rank_>
1869{
1870 // Work-around nvcc warning
1871 unsigned int dummy = n_independent_components;
1872 AssertIndexRange(i, dummy);
1873 (void)dummy;
1874
1875 if constexpr (dim == 0)
1876 {
1877 Assert(false,
1878 ExcMessage(
1879 "A tensor with dimension 0 does not store any elements. "
1880 "There is no indexing that can address its elements."));
1881 return {};
1882 }
1883 else
1884 {
1885 TableIndices<rank_> indices;
1886
1887 unsigned int remainder = i;
1888 for (int r = rank_ - 1; r >= 0; --r)
1889 {
1890 indices[r] = remainder % dim;
1891 remainder = remainder / dim;
1892 }
1893 Assert(remainder == 0, ExcInternalError());
1894
1895 return indices;
1896 }
1897}
1898
1899
1900template <int rank_, int dim, typename Number>
1901constexpr inline void
1903{
1904 for (unsigned int i = 0; i < dim; ++i)
1905 values[i] = internal::NumberType<Number>::value(0.0);
1906}
1907
1908
1909template <int rank_, int dim, typename Number>
1910constexpr std::size_t
1912{
1913 return sizeof(Tensor<rank_, dim, Number>);
1914}
1915
1916
1917template <int rank_, int dim, typename Number>
1918template <class Archive>
1919inline void
1920Tensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
1921{
1922 if constexpr (rank_ > 1)
1923 ar &values;
1924 else
1925 ar &boost::serialization::make_array(&values[0], dim);
1926}
1927
1928
1929template <int rank_, int dim, typename Number>
1931
1932#endif // DOXYGEN
1933
1934/* ----------------- Non-member functions operating on tensors. ------------ */
1935
1948template <int rank_, int dim, typename Number>
1949inline std::ostream &
1950operator<<(std::ostream &out, const Tensor<rank_, dim, Number> &p)
1951{
1952 for (unsigned int i = 0; i < dim; ++i)
1953 {
1954 out << p[i];
1955 if (i != dim - 1)
1956 for (unsigned int j = 0; j < rank_; ++j)
1957 out << ' ';
1958 }
1959
1960 return out;
1961}
1962
1963
1970template <int dim, typename Number>
1971inline std::ostream &
1972operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
1973{
1974 out << static_cast<const Number &>(p);
1975 return out;
1976}
1977
1978
1997template <int dim, typename Number, typename Other>
2000 operator*(const Other &object, const Tensor<0, dim, Number> &t)
2001{
2002 return object * static_cast<const Number &>(t);
2003}
2004
2005
2006
2017template <int dim, typename Number, typename Other>
2020 operator*(const Tensor<0, dim, Number> &t, const Other &object)
2021{
2022 return static_cast<const Number &>(t) * object;
2023}
2024
2025
2037template <int dim, typename Number, typename OtherNumber>
2041 const Tensor<0, dim, OtherNumber> &src2)
2042{
2043 return static_cast<const Number &>(src1) *
2044 static_cast<const OtherNumber &>(src2);
2045}
2046
2047
2055template <int dim, typename Number, typename OtherNumber>
2057 Tensor<0,
2058 dim,
2059 typename ProductType<Number,
2060 typename EnableIfScalar<OtherNumber>::type>::type>
2061 operator/(const Tensor<0, dim, Number> &t, const OtherNumber &factor)
2062{
2063 return static_cast<const Number &>(t) / factor;
2064}
2065
2066
2074template <int dim, typename Number, typename OtherNumber>
2079{
2080 return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
2081}
2082
2083
2091template <int dim, typename Number, typename OtherNumber>
2096{
2097 return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
2098}
2099
2100
2113template <int rank, int dim, typename Number, typename OtherNumber>
2115 Tensor<rank,
2116 dim,
2117 typename ProductType<Number,
2118 typename EnableIfScalar<OtherNumber>::type>::type>
2119 operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2120{
2122 tt *= factor;
2123 return tt;
2124}
2125
2126
2139template <int rank, int dim, typename Number, typename OtherNumber>
2141 Tensor<rank,
2142 dim,
2144 OtherNumber>::type>
2145 operator*(const Number &factor, const Tensor<rank, dim, OtherNumber> &t)
2146{
2147 // simply forward to the operator above
2148 return t * factor;
2149}
2150
2151
2152
2162template <int rank, int dim, typename Number, typename OtherNumber>
2164 Tensor<rank,
2165 dim,
2166 typename ProductType<Number,
2167 typename EnableIfScalar<OtherNumber>::type>::type>
2168 operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
2169{
2171 tt /= factor;
2172 return tt;
2173}
2174
2175
2185template <int rank, int dim, typename Number, typename OtherNumber>
2195
2196
2206template <int rank, int dim, typename Number, typename OtherNumber>
2216
2223template <int dim, typename Number, typename OtherNumber>
2224inline constexpr DEAL_II_ALWAYS_INLINE
2227 const Tensor<0, dim, OtherNumber> &src2)
2228{
2230
2231 tmp *= src2;
2232
2233 return tmp;
2234}
2235
2252template <int rank, int dim, typename Number, typename OtherNumber>
2253inline constexpr DEAL_II_ALWAYS_INLINE
2257{
2259
2260 for (unsigned int i = 0; i < dim; ++i)
2263
2264 return tmp;
2265}
2266
2311template <int rank_1,
2312 int rank_2,
2313 int dim,
2314 typename Number,
2315 typename OtherNumber,
2316 typename = std::enable_if_t<rank_1 >= 1 && rank_2 >= 1>>
2317constexpr inline DEAL_II_ALWAYS_INLINE
2318 typename Tensor<rank_1 + rank_2 - 2,
2319 dim,
2320 typename ProductType<Number, OtherNumber>::type>::tensor_type
2323{
2324 // Treat some common cases separately. Specifically, these are the dot
2325 // product between two rank-1 tensors, and the product between a
2326 // rank-2 tensor and a rank-1 tensor. Both of these lead to a linear
2327 // loop over adjacent memory and can be dealt with efficiently; in the
2328 // latter case (rank-2 times rank-1), we implement things by deferring
2329 // to rank-1 times rank-1 dot products.
2330 if constexpr ((rank_1 == 1) && (rank_2 == 1))
2331 {
2332 // This is a dot product between two rank-1 tensors. Write it out as
2333 // a linear loop:
2334 static_assert(dim > 0, "Tensors cannot have dimension zero.");
2335 typename ProductType<Number, OtherNumber>::type sum = src1[0] * src2[0];
2336 for (unsigned int i = 1; i < dim; ++i)
2337 sum += src1[i] * src2[i];
2338
2339 return sum;
2340 }
2341 else if constexpr ((rank_1 == 2) && (rank_2 == 1))
2342 {
2343 // This is a product between a rank-2 and a rank-1 tensor. This
2344 // corresponds to taking dot products between the rows of the former
2345 // and the latter.
2346 typename Tensor<
2347 rank_1 + rank_2 - 2,
2348 dim,
2349 typename ProductType<Number, OtherNumber>::type>::tensor_type result;
2350 for (unsigned int i = 0; i < dim; ++i)
2351 result[i] += src1[i] * src2;
2352
2353 return result;
2354 }
2355 else
2356 {
2357 // Treat all of the other cases using the more general contraction
2358 // machinery.
2359 typename Tensor<
2360 rank_1 + rank_2 - 2,
2361 dim,
2362 typename ProductType<Number, OtherNumber>::type>::tensor_type result{};
2363
2364 TensorAccessors::internal::
2365 ReorderedIndexView<0, rank_2, const Tensor<rank_2, dim, OtherNumber>>
2368 src1,
2369 reordered);
2370
2371 return result;
2372 }
2373}
2374
2375
2404template <int index_1,
2405 int index_2,
2406 int rank_1,
2407 int rank_2,
2408 int dim,
2409 typename Number,
2410 typename OtherNumber>
2411constexpr inline DEAL_II_ALWAYS_INLINE
2412 typename Tensor<rank_1 + rank_2 - 2,
2413 dim,
2414 typename ProductType<Number, OtherNumber>::type>::tensor_type
2417{
2418 Assert(0 <= index_1 && index_1 < rank_1,
2419 ExcMessage(
2420 "The specified index_1 must lie within the range [0,rank_1)"));
2421 Assert(0 <= index_2 && index_2 < rank_2,
2422 ExcMessage(
2423 "The specified index_2 must lie within the range [0,rank_2)"));
2424
2425 using namespace TensorAccessors;
2426 using namespace TensorAccessors::internal;
2427
2428 // Reorder index_1 to the end of src1:
2430 reord_01 = reordered_index_view<index_1, rank_1>(src1);
2431
2432 // Reorder index_2 to the end of src2:
2433 const ReorderedIndexView<index_2,
2434 rank_2,
2436 reord_02 = reordered_index_view<index_2, rank_2>(src2);
2437
2438 typename Tensor<rank_1 + rank_2 - 2,
2439 dim,
2440 typename ProductType<Number, OtherNumber>::type>::tensor_type
2441 result{};
2442 TensorAccessors::contract<1, rank_1, rank_2, dim>(result, reord_01, reord_02);
2443 return result;
2444}
2445
2446
2477template <int index_1,
2478 int index_2,
2479 int index_3,
2480 int index_4,
2481 int rank_1,
2482 int rank_2,
2483 int dim,
2484 typename Number,
2485 typename OtherNumber>
2486constexpr inline
2487 typename Tensor<rank_1 + rank_2 - 4,
2488 dim,
2489 typename ProductType<Number, OtherNumber>::type>::tensor_type
2490 double_contract(const Tensor<rank_1, dim, Number> &src1,
2492{
2493 Assert(0 <= index_1 && index_1 < rank_1,
2494 ExcMessage(
2495 "The specified index_1 must lie within the range [0,rank_1)"));
2496 Assert(0 <= index_3 && index_3 < rank_1,
2497 ExcMessage(
2498 "The specified index_3 must lie within the range [0,rank_1)"));
2499 Assert(index_1 != index_3,
2500 ExcMessage("index_1 and index_3 must not be the same"));
2501 Assert(0 <= index_2 && index_2 < rank_2,
2502 ExcMessage(
2503 "The specified index_2 must lie within the range [0,rank_2)"));
2504 Assert(0 <= index_4 && index_4 < rank_2,
2505 ExcMessage(
2506 "The specified index_4 must lie within the range [0,rank_2)"));
2507 Assert(index_2 != index_4,
2508 ExcMessage("index_2 and index_4 must not be the same"));
2509
2510 using namespace TensorAccessors;
2511 using namespace TensorAccessors::internal;
2512
2513 // Reorder index_1 to the end of src1:
2516
2517 // Reorder index_2 to the end of src2:
2520
2521 // Now, reorder index_3 to the end of src1. We have to make sure to
2522 // preserve the original ordering: index_1 has been removed. If
2523 // index_3 > index_1, we have to use (index_3 - 1) instead:
2525 (index_3 < index_1 ? index_3 : index_3 - 1),
2526 rank_1,
2527 ReorderedIndexView<index_1, rank_1, const Tensor<rank_1, dim, Number>>>
2528 reord_3 =
2529 TensorAccessors::reordered_index_view < index_3 < index_1 ? index_3 :
2530 index_3 - 1,
2531 rank_1 > (reord_1);
2532
2533 // Now, reorder index_4 to the end of src2. We have to make sure to
2534 // preserve the original ordering: index_2 has been removed. If
2535 // index_4 > index_2, we have to use (index_4 - 1) instead:
2537 (index_4 < index_2 ? index_4 : index_4 - 1),
2538 rank_2,
2540 reord_4 =
2541 TensorAccessors::reordered_index_view < index_4 < index_2 ? index_4 :
2542 index_4 - 1,
2543 rank_2 > (reord_2);
2544
2545 typename Tensor<rank_1 + rank_2 - 4,
2546 dim,
2547 typename ProductType<Number, OtherNumber>::type>::tensor_type
2548 result{};
2549 TensorAccessors::contract<2, rank_1, rank_2, dim>(result, reord_3, reord_4);
2550 return result;
2551}
2552
2553
2566template <int rank, int dim, typename Number, typename OtherNumber>
2567constexpr inline DEAL_II_ALWAYS_INLINE
2569 scalar_product(const Tensor<rank, dim, Number> &left,
2570 const Tensor<rank, dim, OtherNumber> &right)
2571{
2574 return result;
2575}
2576
2577
2595template <template <int, int, typename> class TensorT1,
2596 template <int, int, typename>
2597 class TensorT2,
2598 template <int, int, typename>
2599 class TensorT3,
2600 int rank_1,
2601 int rank_2,
2602 int dim,
2603 typename T1,
2604 typename T2,
2605 typename T3>
2606constexpr inline DEAL_II_ALWAYS_INLINE
2608 contract3(const TensorT1<rank_1, dim, T1> &left,
2609 const TensorT2<rank_1 + rank_2, dim, T2> &middle,
2610 const TensorT3<rank_2, dim, T3> &right)
2611{
2612 using return_type =
2615 middle,
2616 right);
2617}
2618
2619
2630template <int rank_1,
2631 int rank_2,
2632 int dim,
2633 typename Number,
2634 typename OtherNumber>
2635constexpr inline DEAL_II_ALWAYS_INLINE
2639{
2640 typename Tensor<rank_1 + rank_2,
2641 dim,
2642 typename ProductType<Number, OtherNumber>::type>::tensor_type
2643 result{};
2645 return result;
2646}
2647
2648
2667template <int dim, typename Number>
2669cross_product_2d(const Tensor<1, dim, Number> &src)
2670{
2671 Assert(dim == 2, ExcInternalError());
2672
2674
2675 result[0] = src[1];
2676 result[1] = -src[0];
2677
2678 return result;
2679}
2680
2681
2691template <int dim, typename Number1, typename Number2>
2692constexpr inline DEAL_II_ALWAYS_INLINE
2694 cross_product_3d(const Tensor<1, dim, Number1> &src1,
2695 const Tensor<1, dim, Number2> &src2)
2696{
2697 Assert(dim == 3, ExcInternalError());
2698
2700
2701 // avoid compiler warnings
2702 constexpr int s0 = 0 % dim;
2703 constexpr int s1 = 1 % dim;
2704 constexpr int s2 = 2 % dim;
2705
2706 result[s0] = src1[s1] * src2[s2] - src1[s2] * src2[s1];
2707 result[s1] = src1[s2] * src2[s0] - src1[s0] * src2[s2];
2708 result[s2] = src1[s0] * src2[s1] - src1[s1] * src2[s0];
2709
2710 return result;
2711}
2712
2713
2727template <int dim, typename Number>
2728constexpr inline DEAL_II_ALWAYS_INLINE Number
2730{
2731 // Compute the determinant using the Laplace expansion of the
2732 // determinant. We expand along the last row.
2733 Number det = internal::NumberType<Number>::value(0.0);
2734
2735 for (unsigned int k = 0; k < dim; ++k)
2736 {
2737 Tensor<2, dim - 1, Number> minor;
2738 for (unsigned int i = 0; i < dim - 1; ++i)
2739 for (unsigned int j = 0; j < dim - 1; ++j)
2740 minor[i][j] = t[i][j < k ? j : j + 1];
2741
2742 const Number cofactor = ((k % 2 == 0) ? -1. : 1.) * determinant(minor);
2743
2744 det += t[dim - 1][k] * cofactor;
2745 }
2746
2747 return ((dim % 2 == 0) ? 1. : -1.) * det;
2748}
2749
2755template <typename Number>
2756constexpr DEAL_II_ALWAYS_INLINE Number
2758{
2759 return t[0][0];
2760}
2761
2767template <typename Number>
2768constexpr DEAL_II_ALWAYS_INLINE Number
2770{
2771 // hard-coded for efficiency reasons
2772 return t[0][0] * t[1][1] - t[1][0] * t[0][1];
2773}
2774
2780template <typename Number>
2781constexpr DEAL_II_ALWAYS_INLINE Number
2783{
2784 // hard-coded for efficiency reasons
2785 const Number C0 = internal::NumberType<Number>::value(t[1][1] * t[2][2]) -
2786 internal::NumberType<Number>::value(t[1][2] * t[2][1]);
2787 const Number C1 = internal::NumberType<Number>::value(t[1][2] * t[2][0]) -
2788 internal::NumberType<Number>::value(t[1][0] * t[2][2]);
2789 const Number C2 = internal::NumberType<Number>::value(t[1][0] * t[2][1]) -
2790 internal::NumberType<Number>::value(t[1][1] * t[2][0]);
2791 return t[0][0] * C0 + t[0][1] * C1 + t[0][2] * C2;
2792}
2793
2794
2801template <int dim, typename Number>
2802constexpr inline DEAL_II_ALWAYS_INLINE Number
2804{
2805 Number t = d[0][0];
2806 for (unsigned int i = 1; i < dim; ++i)
2807 t += d[i][i];
2808 return t;
2809}
2810
2811
2820template <int dim, typename Number>
2821constexpr inline Tensor<2, dim, Number>
2823{
2824 Number return_tensor[dim][dim];
2825
2826 // if desired, take over the
2827 // inversion of a 4x4 tensor
2828 // from the FullMatrix
2830
2831 return Tensor<2, dim, Number>(return_tensor);
2832}
2833
2834
2835#ifndef DOXYGEN
2836
2837template <typename Number>
2840{
2841 Tensor<2, 1, Number> return_tensor;
2842
2843 return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
2844
2845 return return_tensor;
2846}
2847
2848
2849template <typename Number>
2852{
2853 Tensor<2, 2, Number> return_tensor;
2854
2855 const Number inv_det_t = internal::NumberType<Number>::value(
2856 1.0 / (t[0][0] * t[1][1] - t[1][0] * t[0][1]));
2857 return_tensor[0][0] = t[1][1];
2858 return_tensor[0][1] = -t[0][1];
2859 return_tensor[1][0] = -t[1][0];
2860 return_tensor[1][1] = t[0][0];
2861 return_tensor *= inv_det_t;
2862
2863 return return_tensor;
2864}
2865
2866template <typename Number>
2869{
2870 Tensor<2, 3, Number> return_tensor;
2871
2872 const auto value = [](const auto &t) {
2874 };
2875
2876 return_tensor[0][0] = value(t[1][1] * t[2][2]) - value(t[1][2] * t[2][1]);
2877 return_tensor[0][1] = value(t[0][2] * t[2][1]) - value(t[0][1] * t[2][2]);
2878 return_tensor[0][2] = value(t[0][1] * t[1][2]) - value(t[0][2] * t[1][1]);
2879 return_tensor[1][0] = value(t[1][2] * t[2][0]) - value(t[1][0] * t[2][2]);
2880 return_tensor[1][1] = value(t[0][0] * t[2][2]) - value(t[0][2] * t[2][0]);
2881 return_tensor[1][2] = value(t[0][2] * t[1][0]) - value(t[0][0] * t[1][2]);
2882 return_tensor[2][0] = value(t[1][0] * t[2][1]) - value(t[1][1] * t[2][0]);
2883 return_tensor[2][1] = value(t[0][1] * t[2][0]) - value(t[0][0] * t[2][1]);
2884 return_tensor[2][2] = value(t[0][0] * t[1][1]) - value(t[0][1] * t[1][0]);
2885
2886 const Number inv_det_t =
2887 value(1.0 / (t[0][0] * return_tensor[0][0] + t[0][1] * return_tensor[1][0] +
2888 t[0][2] * return_tensor[2][0]));
2889 return_tensor *= inv_det_t;
2890
2891 return return_tensor;
2892}
2893
2894#endif /* DOXYGEN */
2895
2896
2902template <int dim, typename Number>
2905{
2907 for (unsigned int i = 0; i < dim; ++i)
2908 {
2909 tt[i][i] = t[i][i];
2910 for (unsigned int j = i + 1; j < dim; ++j)
2911 {
2912 tt[i][j] = t[j][i];
2913 tt[j][i] = t[i][j];
2914 };
2915 }
2916 return tt;
2917}
2918
2919
2933template <int dim, typename Number>
2934constexpr Tensor<2, dim, Number>
2935adjugate(const Tensor<2, dim, Number> &t)
2936{
2937 return determinant(t) * invert(t);
2938}
2939
2940
2954template <int dim, typename Number>
2955constexpr Tensor<2, dim, Number>
2956cofactor(const Tensor<2, dim, Number> &t)
2957{
2958 return transpose(adjugate(t));
2959}
2960
2961
3025template <int dim, typename Number>
3028
3029
3037template <int dim, typename Number>
3038inline Number
3040{
3041 Number max = internal::NumberType<Number>::value(0.0);
3042 for (unsigned int j = 0; j < dim; ++j)
3043 {
3044 Number sum = internal::NumberType<Number>::value(0.0);
3045 for (unsigned int i = 0; i < dim; ++i)
3046 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3047
3048 if (sum > max)
3049 max = sum;
3050 }
3051
3052 return max;
3053}
3054
3055
3063template <int dim, typename Number>
3064inline Number
3066{
3067 Number max = internal::NumberType<Number>::value(0.0);
3068 for (unsigned int i = 0; i < dim; ++i)
3069 {
3070 Number sum = internal::NumberType<Number>::value(0.0);
3071 for (unsigned int j = 0; j < dim; ++j)
3072 sum += numbers::NumberTraits<Number>::abs(t[i][j]);
3073
3074 if (sum > max)
3075 max = sum;
3076 }
3077
3078 return max;
3079}
3080
3086#ifndef DOXYGEN
3087
3088
3089# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3090
3091// Specialization of functions for ADOL-C number types when
3092// the advanced branching feature is used
3093template <int dim>
3094inline adouble
3096{
3097 adouble max = internal::NumberType<adouble>::value(0.0);
3098 for (unsigned int j = 0; j < dim; ++j)
3099 {
3100 adouble sum = internal::NumberType<adouble>::value(0.0);
3101 for (unsigned int i = 0; i < dim; ++i)
3102 sum += fabs(t[i][j]);
3103
3104 condassign(max, (sum > max), sum, max);
3105 }
3106
3107 return max;
3108}
3109
3110
3111template <int dim>
3112inline adouble
3114{
3116 for (unsigned int i = 0; i < dim; ++i)
3117 {
3119 for (unsigned int j = 0; j < dim; ++j)
3120 sum += fabs(t[i][j]);
3121
3122 condassign(max, (sum > max), sum, max);
3123 }
3124
3125 return max;
3126}
3127
3128# endif // DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
3129
3130
3131#endif // DOXYGEN
3132
3134
3135#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:2168
constexpr Tensor(const Tensor< 1, dim, Tensor< rank_ - 1, dim, OtherNumber > > &initializer)
constexpr const Number & operator[](const TableIndices< rank_ > &indices) 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:528
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:2061
Number * begin_raw()
static constexpr unsigned int rank
Definition tensor.h:490
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
Definition tensor.h:514
constexpr Tensor(const Tensor< rank_, dim, OtherNumber > &initializer)
numbers::NumberTraits< Number >::real_type norm() const
constexpr Tensor & operator-=(const Tensor< rank_, dim, OtherNumber > &)
constexpr void clear()
void unroll(const Iterator begin, const Iterator end) const
void unroll(Vector< OtherNumber > &result) const
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:2040
const Number * begin_raw() const
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:2226
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:2094
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:2188
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3065
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:2000
constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor()
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3039
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:2145
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:2077
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:2255
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:2209
constexpr ProductType< Number, Other >::type operator*(const Tensor< 0, dim, Number > &t, const Other &object)
Definition tensor.h:2020
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:2119
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:503
#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:504
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
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: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
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
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:1950
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:2226
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3065
constexpr ProductType< Other, Number >::type operator*(const Other &object, const Tensor< 0, dim, Number > &t)
Definition tensor.h:2000
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3039