Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
symmetric_tensor.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 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_symmetric_tensor_h
16#define dealii_symmetric_tensor_h
17
18
19#include <deal.II/base/config.h>
20
26#include <deal.II/base/tensor.h>
27#include <deal.II/base/types.h>
28
29#include <array>
30
32
33// Forward declaration
34#ifndef DOXYGEN
35template <int rank, int dim, typename Number = double>
36class SymmetricTensor;
37#endif
38
49template <int dim, typename Number = double>
53
82template <int dim, typename Number = double>
86
124template <int dim, typename Number = double>
128
129template <int dim, typename Number>
132
133template <int dim, typename Number>
136
146template <int dim2, typename Number>
147DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
149
160template <int dim, typename Number>
161DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
164
178template <int dim, typename Number>
181
182
183
184namespace internal
185{
186 // Workaround: The following 4 overloads are necessary to be able to
187 // compile the library with Apple Clang 8 and older. We should remove
188 // these overloads again when we bump the minimal required version to
189 // something later than clang-3.6 / Apple Clang 6.3.
190 template <int rank, int dim, typename T, typename U>
191 struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
192 {
193 using type =
194 SymmetricTensor<rank,
195 dim,
196 std::complex<typename ProductType<T, U>::type>>;
197 };
198
199 template <int rank, int dim, typename T, typename U>
200 struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
201 std::complex<U>>
202 {
203 using type =
204 SymmetricTensor<rank,
205 dim,
206 std::complex<typename ProductType<T, U>::type>>;
207 };
208
209 template <typename T, int rank, int dim, typename U>
210 struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
211 {
212 using type =
213 SymmetricTensor<rank,
214 dim,
215 std::complex<typename ProductType<T, U>::type>>;
216 };
217
218 template <int rank, int dim, typename T, typename U>
219 struct ProductTypeImpl<std::complex<T>,
220 SymmetricTensor<rank, dim, std::complex<U>>>
221 {
222 using type =
223 SymmetricTensor<rank,
224 dim,
225 std::complex<typename ProductType<T, U>::type>>;
226 };
227 // end workaround
228
233 namespace SymmetricTensorImplementation
234 {
239 template <int rank, int dim, typename Number>
240 struct Inverse;
241 } // namespace SymmetricTensorImplementation
242
247 namespace SymmetricTensorAccessors
248 {
257 merge(const TableIndices<2> &previous_indices,
258 const unsigned int new_index,
259 const unsigned int position)
260 {
261 AssertIndexRange(position, 2);
262
263 if (position == 0)
264 return {new_index, numbers::invalid_unsigned_int};
265 else
266 return {previous_indices[0], new_index};
267 }
268
269
270
279 merge(const TableIndices<4> &previous_indices,
280 const unsigned int new_index,
281 const unsigned int position)
282 {
283 AssertIndexRange(position, 4);
284
285 switch (position)
286 {
287 case 0:
288 return {new_index,
292 case 1:
293 return {previous_indices[0],
294 new_index,
297 case 2:
298 return {previous_indices[0],
299 previous_indices[1],
300 new_index,
302 case 3:
303 return {previous_indices[0],
304 previous_indices[1],
305 previous_indices[2],
306 new_index};
307 default:
309 return {};
310 }
311 }
312
313
320 template <int rank1,
321 int rank2,
322 int dim,
323 typename Number,
324 typename OtherNumber = Number>
326 {
328 using type =
329 ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
330 };
331
332
339 template <int dim, typename Number, typename OtherNumber>
340 struct double_contraction_result<2, 2, dim, Number, OtherNumber>
341 {
343 };
344
345
346
359 template <int rank, int dim, typename Number>
361
365 template <int dim, typename Number>
366 struct StorageType<2, dim, Number>
367 {
372 static const unsigned int n_independent_components =
373 (dim * dim + dim) / 2;
374
379 };
380
381
382
386 template <int dim, typename Number>
387 struct StorageType<4, dim, Number>
388 {
394 static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
395
399 static const unsigned int n_independent_components =
400 (n_rank2_components *
402
410 };
411
412
413
418 template <int rank, int dim, bool constness, typename Number>
420
427 template <int rank, int dim, typename Number>
428 struct AccessorTypes<rank, dim, true, Number>
429 {
430 using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
431
432 using reference = const Number &;
433 };
434
441 template <int rank, int dim, typename Number>
442 struct AccessorTypes<rank, dim, false, Number>
443 {
445
446 using reference = Number &;
447 };
448
449
482 template <int rank, int dim, bool constness, int P, typename Number>
484 {
485 public:
489 using reference =
493
494 private:
516
521 constexpr DEAL_II_ALWAYS_INLINE
522 Accessor(const Accessor &) = default;
523
524 public:
529 constexpr Accessor<rank, dim, constness, P - 1, Number>
530 operator[](const unsigned int i);
531
536 constexpr Accessor<rank, dim, constness, P - 1, Number>
537 operator[](const unsigned int i) const;
538
539 private:
545
546 // Declare some other classes as friends. Make sure to work around bugs
547 // in some compilers:
548 template <int, int, typename>
549 friend class ::SymmetricTensor;
550 template <int, int, bool, int, typename>
551 friend class Accessor;
552 friend class ::SymmetricTensor<rank, dim, Number>;
553 friend class Accessor<rank, dim, constness, P + 1, Number>;
554 };
555
556
557
565 template <int rank, int dim, bool constness, typename Number>
566 class Accessor<rank, dim, constness, 1, Number>
567 {
568 public:
572 using reference =
576
577 private:
602
607 constexpr DEAL_II_ALWAYS_INLINE
608 Accessor(const Accessor &) = default;
609
610 public:
615 constexpr reference
616 operator[](const unsigned int);
617
622 constexpr reference
623 operator[](const unsigned int) const;
624
625 private:
631
632 // Declare some other classes as friends. Make sure to work around bugs
633 // in some compilers:
634 template <int, int, typename>
635 friend class ::SymmetricTensor;
636 template <int, int, bool, int, typename>
638 friend class ::SymmetricTensor<rank, dim, Number>;
639 friend class SymmetricTensorAccessors::
640 Accessor<rank, dim, constness, 2, Number>;
641 };
642 } // namespace SymmetricTensorAccessors
643} // namespace internal
644
645
646
719template <int rank_, int dim, typename Number>
721{
722public:
723 static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
724
733 static constexpr unsigned int dimension = dim;
734
738 static const unsigned int rank = rank_;
739
745 static constexpr unsigned int n_independent_components =
747 n_independent_components;
748
753 constexpr DEAL_II_ALWAYS_INLINE
754 SymmetricTensor() = default;
755
769 template <typename OtherNumber>
771
788 constexpr SymmetricTensor(const Number (&array)[n_independent_components]);
789
795 template <typename OtherNumber>
796 DEAL_II_HOST constexpr explicit SymmetricTensor(
798
805 template <typename OtherNumber>
806 DEAL_II_HOST constexpr SymmetricTensor &
808
816 constexpr SymmetricTensor &
817 operator=(const Number &d);
818
824 constexpr operator Tensor<rank_, dim, Number>() const;
825
830 constexpr bool
832
837 constexpr bool
839
843 template <typename OtherNumber>
844 DEAL_II_HOST constexpr SymmetricTensor &
846
850 template <typename OtherNumber>
851 DEAL_II_HOST constexpr SymmetricTensor &
853
858 template <typename OtherNumber>
859 DEAL_II_HOST constexpr SymmetricTensor &
860 operator*=(const OtherNumber &factor);
861
865 template <typename OtherNumber>
866 DEAL_II_HOST constexpr SymmetricTensor &
867 operator/=(const OtherNumber &factor);
868
873 constexpr SymmetricTensor
874 operator-() const;
875
928 template <typename OtherNumber>
929 DEAL_II_HOST DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
930 double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
932
937 template <typename OtherNumber>
938 DEAL_II_HOST DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
939 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
941
946 constexpr Number &
948
953 constexpr const Number &
954 operator()(const TableIndices<rank_> &indices) const;
955
961 constexpr internal::SymmetricTensorAccessors::
962 Accessor<rank_, dim, true, rank_ - 1, Number>
963 operator[](const unsigned int row) const;
964
970 constexpr internal::SymmetricTensorAccessors::
971 Accessor<rank_, dim, false, rank_ - 1, Number>
972 operator[](const unsigned int row);
973
980 constexpr const Number &
981 operator[](const TableIndices<rank_> &indices) const;
982
989 constexpr Number &
991
999 constexpr const Number &
1000 access_raw_entry(const unsigned int unrolled_index) const;
1001
1009 constexpr Number &
1010 access_raw_entry(const unsigned int unrolled_index);
1011
1023 norm() const;
1024
1032 static DEAL_II_HOST constexpr unsigned int
1034
1040 static DEAL_II_HOST constexpr TableIndices<rank_>
1041 unrolled_to_component_indices(const unsigned int i);
1042
1056 constexpr void
1058
1063 static DEAL_II_HOST constexpr std::size_t
1065
1071 template <class Archive>
1072 void
1073 serialize(Archive &ar, const unsigned int version);
1074
1075private:
1081
1085 using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
1086
1091
1092#ifndef DOXYGEN
1093
1094 // Make all other symmetric tensors friends.
1095 template <int, int, typename>
1096 friend class SymmetricTensor;
1097
1098 // Make a few more functions friends.
1099 template <int dim2, typename Number2>
1100 friend DEAL_II_HOST constexpr Number2
1102
1103 template <int dim2, typename Number2>
1104 friend DEAL_II_HOST DEAL_II_CONSTEXPR Number2
1106
1107 template <int dim2, typename Number2>
1110
1111 template <int dim2, typename Number2>
1114
1115 template <int dim2, typename Number2>
1118
1119 template <int dim2, typename Number2>
1122
1123
1124 // Make a few helper classes friends as well.
1126 Inverse<2, dim, Number>;
1127
1129 Inverse<4, dim, Number>;
1130#endif
1131};
1132
1133
1134
1135// ------------------------- inline functions ------------------------
1136
1137#ifndef DOXYGEN
1138
1139// provide declarations for static members
1140template <int rank, int dim, typename Number>
1141const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1142
1143template <int rank_, int dim, typename Number>
1144constexpr unsigned int
1145 SymmetricTensor<rank_, dim, Number>::n_independent_components;
1146
1147namespace internal
1148{
1149 namespace SymmetricTensorAccessors
1150 {
1151 template <int rank_, int dim, bool constness, int P, typename Number>
1153 Accessor<rank_, dim, constness, P, Number>::Accessor(
1154 tensor_type &tensor,
1155 const TableIndices<rank_> &previous_indices)
1156 : tensor(tensor)
1157 , previous_indices(previous_indices)
1158 {}
1159
1160
1161
1162 template <int rank_, int dim, bool constness, int P, typename Number>
1163 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1164 Accessor<rank_, dim, constness, P - 1, Number>
1165 Accessor<rank_, dim, constness, P, Number>::operator[](
1166 const unsigned int i)
1167 {
1168 return Accessor<rank_, dim, constness, P - 1, Number>(
1169 tensor, merge(previous_indices, i, rank_ - P));
1170 }
1171
1172
1173
1174 template <int rank_, int dim, bool constness, int P, typename Number>
1176 Accessor<rank_, dim, constness, P - 1, Number>
1177 Accessor<rank_, dim, constness, P, Number>::operator[](
1178 const unsigned int i) const
1179 {
1180 return Accessor<rank_, dim, constness, P - 1, Number>(
1181 tensor, merge(previous_indices, i, rank_ - P));
1182 }
1183
1184
1185
1186 template <int rank_, int dim, bool constness, typename Number>
1188 Accessor<rank_, dim, constness, 1, Number>::Accessor(
1189 tensor_type &tensor,
1190 const TableIndices<rank_> &previous_indices)
1191 : tensor(tensor)
1192 , previous_indices(previous_indices)
1193 {}
1194
1195
1196
1197 template <int rank_, int dim, bool constness, typename Number>
1198 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1199 typename Accessor<rank_, dim, constness, 1, Number>::reference
1200 Accessor<rank_, dim, constness, 1, Number>::operator[](
1201 const unsigned int i)
1202 {
1203 return tensor(merge(previous_indices, i, rank_ - 1));
1204 }
1205
1206
1207 template <int rank_, int dim, bool constness, typename Number>
1209 typename Accessor<rank_, dim, constness, 1, Number>::reference
1210 Accessor<rank_, dim, constness, 1, Number>::operator[](
1211 const unsigned int i) const
1212 {
1213 return tensor(merge(previous_indices, i, rank_ - 1));
1214 }
1215 } // namespace SymmetricTensorAccessors
1216} // namespace internal
1217
1218
1219
1220template <int rank_, int dim, typename Number>
1221template <typename OtherNumber>
1225{
1226 static_assert(rank == 2, "This function is only implemented for rank==2");
1227 for (unsigned int d = 0; d < dim; ++d)
1228 for (unsigned int e = 0; e < d; ++e)
1229 Assert(t[d][e] == t[e][d],
1230 ExcMessage("The incoming Tensor must be exactly symmetric."));
1231
1232 for (unsigned int d = 0; d < dim; ++d)
1233 data[d] = t[d][d];
1234
1235 for (unsigned int d = 0, c = 0; d < dim; ++d)
1236 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1237 data[dim + c] = t[d][e];
1238}
1239
1240
1241
1242template <int rank_, int dim, typename Number>
1243template <typename OtherNumber>
1247 : data(initializer.data)
1248{}
1249
1250
1251
1252template <int rank_, int dim, typename Number>
1253DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1255 const Number (&array)[n_independent_components])
1256 : data(
1257 *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1258{
1259 // ensure that the reinterpret_cast above actually works
1260 Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1262}
1263
1264
1265
1266template <int rank_, int dim, typename Number>
1267template <typename OtherNumber>
1268DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1272{
1273 data = t.data;
1274 return *this;
1275}
1276
1277
1278
1279template <int rank_, int dim, typename Number>
1280DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1283{
1285 ExcMessage("Only assignment with zero is allowed"));
1286 (void)d;
1287
1289
1290 return *this;
1291}
1292
1293
1294namespace internal
1295{
1296 namespace SymmetricTensorImplementation
1297 {
1298 template <int dim, typename Number>
1299 constexpr inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1300 convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1301 {
1303
1304 // diagonal entries are stored first
1305 for (unsigned int d = 0; d < dim; ++d)
1306 t[d][d] = s.access_raw_entry(d);
1307
1308 // off-diagonal entries come next, row by row
1309 for (unsigned int d = 0, c = 0; d < dim; ++d)
1310 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1311 {
1312 t[d][e] = s.access_raw_entry(dim + c);
1313 t[e][d] = s.access_raw_entry(dim + c);
1314 }
1315 return t;
1316 }
1317
1318
1319 template <int dim, typename Number>
1320 constexpr ::Tensor<4, dim, Number>
1321 convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1322 {
1323 // utilize the symmetry properties of SymmetricTensor<4,dim>
1324 // discussed in the class documentation to avoid accessing all
1325 // independent elements of the input tensor more than once
1327
1328 for (unsigned int i = 0; i < dim; ++i)
1329 for (unsigned int j = i; j < dim; ++j)
1330 for (unsigned int k = 0; k < dim; ++k)
1331 for (unsigned int l = k; l < dim; ++l)
1332 t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1333 t[TableIndices<4>(j, i, k, l)] =
1334 t[TableIndices<4>(j, i, l, k)] =
1335 st[TableIndices<4>(i, j, k, l)];
1336
1337 return t;
1338 }
1339
1340
1341 template <typename Number>
1342 struct Inverse<2, 1, Number>
1343 {
1344 constexpr static inline DEAL_II_ALWAYS_INLINE
1345 ::SymmetricTensor<2, 1, Number>
1346 value(const ::SymmetricTensor<2, 1, Number> &t)
1347 {
1349
1350 tmp[0][0] = 1.0 / t[0][0];
1351
1352 return tmp;
1353 }
1354 };
1355
1356
1357 template <typename Number>
1358 struct Inverse<2, 2, Number>
1359 {
1360 constexpr static inline DEAL_II_ALWAYS_INLINE
1361 ::SymmetricTensor<2, 2, Number>
1362 value(const ::SymmetricTensor<2, 2, Number> &t)
1363 {
1365
1366 // Sympy result: ([
1367 // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1368 // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1369 const TableIndices<2> idx_00(0, 0);
1370 const TableIndices<2> idx_01(0, 1);
1371 const TableIndices<2> idx_11(1, 1);
1372 const Number inv_det_t =
1373 1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1374 tmp[idx_00] = t[idx_11];
1375 tmp[idx_01] = -t[idx_01];
1376 tmp[idx_11] = t[idx_00];
1377 tmp *= inv_det_t;
1378
1379 return tmp;
1380 }
1381 };
1382
1383
1384 template <typename Number>
1385 struct Inverse<2, 3, Number>
1386 {
1387 constexpr static ::SymmetricTensor<2, 3, Number>
1388 value(const ::SymmetricTensor<2, 3, Number> &t)
1389 {
1391
1392 // Sympy result: ([
1393 // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1394 // 2*t01*t02*t12 - t02**2*t11),
1395 // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1396 // 2*t01*t02*t12 - t02**2*t11),
1397 // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1398 // 2*t01*t02*t12 - t02**2*t11)],
1399 // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1400 // 2*t01*t02*t12 - t02**2*t11),
1401 // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1402 // 2*t01*t02*t12 - t02**2*t11),
1403 // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1404 // 2*t01*t02*t12 + t02**2*t11)],
1405 // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1406 // 2*t01*t02*t12 - t02**2*t11),
1407 // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1408 // 2*t01*t02*t12 + t02**2*t11),
1409 // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1410 // 2*t01*t02*t12 + t02**2*t11)] ])
1411 //
1412 // =
1413 //
1414 // [ (t11*t22 - t12**2)/det_t,
1415 // (-t01*t22 + t02*t12)/det_t,
1416 // (t01*t12 - t02*t11)/det_t],
1417 // [ (-t01*t22 + t02*t12)/det_t,
1418 // (t00*t22 - t02**2)/det_t,
1419 // (-t00*t12 + t01*t02)/det_t],
1420 // [ (t01*t12 - t02*t11)/det_t,
1421 // (-t00*t12 + t01*t02)/det_t,
1422 // (t00*t11 - t01**2)/det_t] ])
1423 //
1424 // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1425 // 2*t01*t02*t12 - t02**2*t11)
1426 const TableIndices<2> idx_00(0, 0);
1427 const TableIndices<2> idx_01(0, 1);
1428 const TableIndices<2> idx_02(0, 2);
1429 const TableIndices<2> idx_11(1, 1);
1430 const TableIndices<2> idx_12(1, 2);
1431 const TableIndices<2> idx_22(2, 2);
1432 const Number inv_det_t =
1433 1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1434 t[idx_00] * t[idx_12] * t[idx_12] -
1435 t[idx_01] * t[idx_01] * t[idx_22] +
1436 2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1437 t[idx_02] * t[idx_02] * t[idx_11]);
1438 tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1439 tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1440 tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1441 tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1442 tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1443 tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1444 tmp *= inv_det_t;
1445
1446 return tmp;
1447 }
1448 };
1449
1450
1451 template <typename Number>
1452 struct Inverse<4, 1, Number>
1453 {
1454 constexpr static inline ::SymmetricTensor<4, 1, Number>
1455 value(const ::SymmetricTensor<4, 1, Number> &t)
1456 {
1458 tmp.data[0][0] = 1.0 / t.data[0][0];
1459 return tmp;
1460 }
1461 };
1462
1463
1464 template <typename Number>
1465 struct Inverse<4, 2, Number>
1466 {
1467 constexpr static inline ::SymmetricTensor<4, 2, Number>
1468 value(const ::SymmetricTensor<4, 2, Number> &t)
1469 {
1471
1472 // Inverting this tensor is a little more complicated than necessary,
1473 // since we store the data of 't' as a 3x3 matrix t.data, but the
1474 // product between a rank-4 and a rank-2 tensor is really not the
1475 // product between this matrix and the 3-vector of a rhs, but rather
1476 //
1477 // B.vec = t.data * mult * A.vec
1478 //
1479 // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1480 // capture the fact that we need to add up both the c_ij12*a_12 and the
1481 // c_ij21*a_21 terms.
1482 //
1483 // In addition, in this scheme, the identity tensor has the matrix
1484 // representation mult^-1.
1485 //
1486 // The inverse of 't' therefore has the matrix representation
1487 //
1488 // inv.data = mult^-1 * t.data^-1 * mult^-1
1489 //
1490 // in order to compute it, let's first compute the inverse of t.data and
1491 // put it into tmp.data; at the end of the function we then scale the
1492 // last row and column of the inverse by 1/2, corresponding to the left
1493 // and right multiplication with mult^-1.
1494 const Number t4 = t.data[0][0] * t.data[1][1],
1495 t6 = t.data[0][0] * t.data[1][2],
1496 t8 = t.data[0][1] * t.data[1][0],
1497 t00 = t.data[0][2] * t.data[1][0],
1498 t01 = t.data[0][1] * t.data[2][0],
1499 t04 = t.data[0][2] * t.data[2][0],
1500 t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1501 t8 * t.data[2][2] + t00 * t.data[2][1] +
1502 t01 * t.data[1][2] - t04 * t.data[1][1]);
1503 tmp.data[0][0] =
1504 (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1505 tmp.data[0][1] =
1506 -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1507 tmp.data[0][2] =
1508 -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1509 tmp.data[1][0] =
1510 -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1511 tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1512 tmp.data[1][2] = -(t6 - t00) * t07;
1513 tmp.data[2][0] =
1514 -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1515 tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1516 tmp.data[2][2] = (t4 - t8) * t07;
1517
1518 // scale last row and column as mentioned
1519 // above
1520 tmp.data[2][0] /= 2;
1521 tmp.data[2][1] /= 2;
1522 tmp.data[0][2] /= 2;
1523 tmp.data[1][2] /= 2;
1524 tmp.data[2][2] /= 4;
1525
1526 return tmp;
1527 }
1528 };
1529
1530
1531 template <typename Number>
1532 struct Inverse<4, 3, Number>
1533 {
1534 static ::SymmetricTensor<4, 3, Number>
1535 value(const ::SymmetricTensor<4, 3, Number> &t)
1536 {
1538
1539 // This function follows the exact same scheme as the 2d case, except
1540 // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1541 // Instead, we use the Gauss-Jordan algorithm implemented for
1542 // FullMatrix. For historical reasons the following code is copied from
1543 // there, with the tangential benefit that we do not need to copy the
1544 // tensor entries to and from the FullMatrix.
1545 const unsigned int N = 6;
1546
1547 // First get an estimate of the size of the elements of this matrix,
1548 // for later checks whether the pivot element is large enough, or
1549 // whether we have to fear that the matrix is not regular.
1550 Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1551 for (unsigned int i = 0; i < N; ++i)
1552 diagonal_sum += numbers::NumberTraits<Number>::abs(tmp.data[i][i]);
1553 const Number typical_diagonal_element =
1554 diagonal_sum / static_cast<double>(N);
1555 (void)typical_diagonal_element;
1556
1557 unsigned int p[N];
1558 for (unsigned int i = 0; i < N; ++i)
1559 p[i] = i;
1560
1561 for (unsigned int j = 0; j < N; ++j)
1562 {
1563 // Pivot search: search that part of the line on and right of the
1564 // diagonal for the largest element.
1565 Number max = numbers::NumberTraits<Number>::abs(tmp.data[j][j]);
1566 unsigned int r = j;
1567 for (unsigned int i = j + 1; i < N; ++i)
1568 if (numbers::NumberTraits<Number>::abs(tmp.data[i][j]) > max)
1569 {
1571 r = i;
1572 }
1573
1574 // Check whether the pivot is too small
1575 Assert(max > 1.e-16 * typical_diagonal_element,
1576 ExcMessage("This tensor seems to be noninvertible"));
1577
1578 // Row interchange
1579 if (r > j)
1580 {
1581 for (unsigned int k = 0; k < N; ++k)
1582 std::swap(tmp.data[j][k], tmp.data[r][k]);
1583
1584 std::swap(p[j], p[r]);
1585 }
1586
1587 // Transformation
1588 const Number hr = 1. / tmp.data[j][j];
1589 tmp.data[j][j] = hr;
1590 for (unsigned int k = 0; k < N; ++k)
1591 {
1592 if (k == j)
1593 continue;
1594 for (unsigned int i = 0; i < N; ++i)
1595 {
1596 if (i == j)
1597 continue;
1598 tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1599 }
1600 }
1601 for (unsigned int i = 0; i < N; ++i)
1602 {
1603 tmp.data[i][j] *= hr;
1604 tmp.data[j][i] *= -hr;
1605 }
1606 tmp.data[j][j] = hr;
1607 }
1608
1609 // Column interchange
1610 Number hv[N];
1611 for (unsigned int i = 0; i < N; ++i)
1612 {
1613 for (unsigned int k = 0; k < N; ++k)
1614 hv[p[k]] = tmp.data[i][k];
1615 for (unsigned int k = 0; k < N; ++k)
1616 tmp.data[i][k] = hv[k];
1617 }
1618
1619 // Scale rows and columns. The mult matrix
1620 // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1621 for (unsigned int i = 3; i < 6; ++i)
1622 for (unsigned int j = 0; j < 3; ++j)
1623 tmp.data[i][j] /= 2;
1624
1625 for (unsigned int i = 0; i < 3; ++i)
1626 for (unsigned int j = 3; j < 6; ++j)
1627 tmp.data[i][j] /= 2;
1628
1629 for (unsigned int i = 3; i < 6; ++i)
1630 for (unsigned int j = 3; j < 6; ++j)
1631 tmp.data[i][j] /= 4;
1632
1633 return tmp;
1634 }
1635 };
1636
1637 } // namespace SymmetricTensorImplementation
1638} // namespace internal
1639
1640
1641
1642template <int rank_, int dim, typename Number>
1645 const
1646{
1647 return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1648}
1649
1650
1651
1652template <int rank_, int dim, typename Number>
1653DEAL_II_HOST constexpr bool
1656{
1657 return data == t.data;
1658}
1659
1660
1661
1662template <int rank_, int dim, typename Number>
1663DEAL_II_HOST constexpr bool
1666{
1667 return data != t.data;
1668}
1669
1670
1671
1672template <int rank_, int dim, typename Number>
1673template <typename OtherNumber>
1674DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1678{
1679 data += t.data;
1680 return *this;
1681}
1682
1683
1684
1685template <int rank_, int dim, typename Number>
1686template <typename OtherNumber>
1687DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1691{
1692 data -= t.data;
1693 return *this;
1694}
1695
1696
1697
1698template <int rank_, int dim, typename Number>
1699template <typename OtherNumber>
1700DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1703{
1704 data *= d;
1705 return *this;
1706}
1707
1708
1709
1710template <int rank_, int dim, typename Number>
1711template <typename OtherNumber>
1712DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1715{
1716 data /= d;
1717 return *this;
1718}
1719
1720
1721
1722template <int rank_, int dim, typename Number>
1723DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1726{
1727 SymmetricTensor tmp = *this;
1728 tmp.data = -tmp.data;
1729 return tmp;
1730}
1731
1732
1733
1734template <int rank_, int dim, typename Number>
1735DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
1737{
1738 data.clear();
1739}
1740
1741
1742
1743template <int rank_, int dim, typename Number>
1744DEAL_II_HOST constexpr std::size_t
1746{
1747 // all memory consists of statically allocated memory of the current
1748 // object, no pointers
1750}
1751
1752
1753
1754namespace internal
1755{
1759 template <int dim, typename Number, typename OtherNumber = Number>
1761 typename SymmetricTensorAccessors::
1762 double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1763 perform_double_contraction(
1764 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1765 base_tensor_type &data,
1766 const typename SymmetricTensorAccessors::
1767 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1768 {
1769 using result_type = typename SymmetricTensorAccessors::
1770 double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1771
1772 switch (dim)
1773 {
1774 case 1:
1775 return data[0] * sdata[0];
1776
1777 default:
1778 // Start with the non-diagonal part. These values appear
1779 // twice in the matrix, but are only stored once. So we can
1780 // get the double-contraction sum for these elements using
1781 // only one multiplication each, and at the end multiplying
1782 // things by 2.
1783 result_type sum = data[dim] * sdata[dim];
1784 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1785 sum += data[d] * sdata[d];
1786 sum += sum; // sum *= 2
1787
1788 // Now add the contributions from the diagonal
1789 for (unsigned int d = 0; d < dim; ++d)
1790 sum += data[d] * sdata[d];
1791 return sum;
1792 }
1793 }
1794
1795
1796
1801 template <int dim, typename Number, typename OtherNumber = Number>
1803 typename SymmetricTensorAccessors::
1804 double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1805 perform_double_contraction(
1806 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1807 base_tensor_type &data,
1808 const typename SymmetricTensorAccessors::
1809 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1810 {
1811 using result_type = typename SymmetricTensorAccessors::
1812 double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1813 using value_type = typename SymmetricTensorAccessors::
1814 double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1815
1816 const unsigned int data_dim = SymmetricTensorAccessors::
1817 StorageType<2, dim, value_type>::n_independent_components;
1818 value_type tmp[data_dim]{};
1819 for (unsigned int i = 0; i < data_dim; ++i)
1820 tmp[i] =
1821 perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1822 return result_type(tmp);
1823 }
1824
1825
1826
1831 template <int dim, typename Number, typename OtherNumber = Number>
1833 typename SymmetricTensorAccessors::StorageType<
1834 2,
1835 dim,
1836 typename SymmetricTensorAccessors::
1837 double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1838 base_tensor_type
1839 perform_double_contraction(
1840 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1841 base_tensor_type &data,
1842 const typename SymmetricTensorAccessors::
1843 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1844 {
1845 using value_type = typename SymmetricTensorAccessors::
1846 double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1847 using base_tensor_type = typename SymmetricTensorAccessors::
1848 StorageType<2, dim, value_type>::base_tensor_type;
1849
1850 base_tensor_type tmp;
1851 for (unsigned int i = 0; i < tmp.dimension; ++i)
1852 {
1853 // Start with the non-diagonal part. These values appear
1854 // twice in the matrix, but are only stored once. So we can
1855 // get the double-contraction sum for these elements using
1856 // only one multiplication each, and at the end multiplying
1857 // things by 2.
1858 value_type sum = data[dim] * sdata[dim][i];
1859 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1860 sum += data[d] * sdata[d][i];
1861 sum += sum; // sum *= 2
1862
1863 // Now add the contributions from the diagonal
1864 for (unsigned int d = 0; d < dim; ++d)
1865 sum += data[d] * sdata[d][i];
1866 tmp[i] = sum;
1867 }
1868 return tmp;
1869 }
1870
1871
1872
1876 template <int dim, typename Number, typename OtherNumber = Number>
1878 typename SymmetricTensorAccessors::StorageType<
1879 4,
1880 dim,
1881 typename SymmetricTensorAccessors::
1882 double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1883 base_tensor_type
1884 perform_double_contraction(
1885 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1886 base_tensor_type &data,
1887 const typename SymmetricTensorAccessors::
1888 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1889 {
1890 using value_type = typename SymmetricTensorAccessors::
1891 double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1892 using base_tensor_type = typename SymmetricTensorAccessors::
1893 StorageType<4, dim, value_type>::base_tensor_type;
1894
1895 const unsigned int data_dim = SymmetricTensorAccessors::
1896 StorageType<2, dim, value_type>::n_independent_components;
1897 base_tensor_type tmp;
1898 for (unsigned int i = 0; i < data_dim; ++i)
1899 for (unsigned int j = 0; j < data_dim; ++j)
1900 {
1901 // Start with the non-diagonal part
1902 for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1903 tmp[i][j] += data[i][d] * sdata[d][j];
1904 tmp[i][j] += tmp[i][j]; // tmp[i][j] *= 2;
1905
1906 // Now add the contributions from the diagonal
1907 for (unsigned int d = 0; d < dim; ++d)
1908 tmp[i][j] += data[i][d] * sdata[d][j];
1909 }
1910 return tmp;
1911 }
1912
1913} // end of namespace internal
1914
1915
1916
1917template <int rank_, int dim, typename Number>
1918template <typename OtherNumber>
1920 typename internal::SymmetricTensorAccessors::
1921 double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1924{
1925 // Dispatch to functions that know the types of the involved
1926 // arguments via overloads.
1927 return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1928 s.data);
1929}
1930
1931
1932
1933template <int rank_, int dim, typename Number>
1934template <typename OtherNumber>
1936 typename internal::SymmetricTensorAccessors::
1937 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1940{
1941 typename internal::SymmetricTensorAccessors::
1942 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1943 tmp.data =
1944 internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1945 s.data);
1946 return tmp;
1947}
1948
1949
1950
1951// internal namespace to switch between the
1952// access of different tensors. There used to
1953// be explicit instantiations before for
1954// different ranks and dimensions, but since
1955// we now allow for templates on the data
1956// type, and since we cannot partially
1957// specialize the implementation, this got
1958// into a separate namespace
1959namespace internal
1960{
1961 namespace SymmetricTensorImplementation
1962 {
1963 // a function to do the unrolling from a set of indices to a
1964 // scalar index into the array in which we store the elements of
1965 // a symmetric tensor
1966 //
1967 // this function is for rank-2 tensors
1968 template <int dim>
1969 constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
1971 {
1972 AssertIndexRange(indices[0], dim);
1973 AssertIndexRange(indices[1], dim);
1974
1975 switch (dim)
1976 {
1977 case 1:
1978 {
1979 return 0;
1980 }
1981 case 2:
1982 {
1983 constexpr ::ndarray<unsigned int, 2, 2> table = {
1984 {{{0, 2}}, {{2, 1}}}};
1985 return table[indices[0]][indices[1]];
1986 }
1987 case 3:
1988 {
1989 constexpr ::ndarray<unsigned int, 3, 3> table = {
1990 {{{0, 3, 4}}, {{3, 1, 5}}, {{4, 5, 2}}}};
1991 return table[indices[0]][indices[1]];
1992 }
1993 case 4:
1994 {
1995 constexpr ::ndarray<unsigned int, 4, 4> table = {
1996 {{{0, 4, 5, 6}},
1997 {{4, 1, 7, 8}},
1998 {{5, 7, 2, 9}},
1999 {{6, 8, 9, 3}}}};
2000 return table[indices[0]][indices[1]];
2001 }
2002 default:
2003 // for the remainder, manually figure out the numbering
2004 {
2005 if (indices[0] == indices[1])
2006 return indices[0];
2007
2008 const TableIndices<2> sorted_indices(
2009 std::min(indices[0], indices[1]),
2010 std::max(indices[0], indices[1]));
2011
2012 // Here (d, e) are the row and column of the symmetric matrix and
2013 // 'dim + c' is the index into the Tensor<1, dim> actually used
2014 // for storage.
2015 unsigned int c = 0;
2016 for (unsigned int d = 0; d < dim; ++d)
2017 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2018 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2019 return dim + c;
2020
2021 // should never get here:
2023 return 0;
2024 }
2025 }
2026 }
2027
2028 // a function to do the unrolling from a set of indices to a
2029 // scalar index into the array in which we store the elements of
2030 // a symmetric tensor
2031 //
2032 // this function is for tensors of ranks not already handled
2033 // above
2034 template <int dim, int rank_>
2035 constexpr inline unsigned int
2037 {
2038 (void)indices;
2041 }
2042 } // namespace SymmetricTensorImplementation
2043
2044 template <int dim, typename Number>
2045 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2046 symmetric_tensor_access(const TableIndices<2> &indices,
2047 typename SymmetricTensorAccessors::
2048 StorageType<2, dim, Number>::base_tensor_type &data)
2049 {
2050 return data[SymmetricTensorImplementation::component_to_unrolled_index<dim>(
2051 indices)];
2052 }
2053
2054
2055
2056 template <int dim, typename Number>
2057 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2058 symmetric_tensor_access(const TableIndices<2> &indices,
2059 const typename SymmetricTensorAccessors::
2060 StorageType<2, dim, Number>::base_tensor_type &data)
2061 {
2062 return data[SymmetricTensorImplementation::component_to_unrolled_index<dim>(
2063 indices)];
2064 }
2065
2066
2067
2068 template <int dim, typename Number>
2069 constexpr inline Number &
2070 symmetric_tensor_access(const TableIndices<4> &indices,
2071 typename SymmetricTensorAccessors::
2072 StorageType<4, dim, Number>::base_tensor_type &data)
2073 {
2074 switch (dim)
2075 {
2076 case 1:
2077 return data[0][0];
2078
2079 case 2:
2080 // each entry of the tensor can be thought of as an entry in a
2081 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2082 // rank-2 tensors. this is the format in which we store rank-4
2083 // tensors. determine which position the present entry is
2084 // stored in
2085 {
2086 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2087 return data[base_index[indices[0]][indices[1]]]
2088 [base_index[indices[2]][indices[3]]];
2089 }
2090 case 3:
2091 // each entry of the tensor can be thought of as an entry in a
2092 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2093 // rank-2 tensors. this is the format in which we store rank-4
2094 // tensors. determine which position the present entry is
2095 // stored in
2096 {
2097 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2098 {3, 1, 5},
2099 {4, 5, 2}};
2100 return data[base_index[indices[0]][indices[1]]]
2101 [base_index[indices[2]][indices[3]]];
2102 }
2103
2104 default:
2106 }
2107
2108 // The code should never reach here.
2109 // We cannot return a static variable, as this class must support number
2110 // types that require no instances of the number type to be in scope during
2111 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2112 return data[0][0];
2113 }
2114
2115
2116 template <int dim, typename Number>
2117 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2118 symmetric_tensor_access(const TableIndices<4> &indices,
2119 const typename SymmetricTensorAccessors::
2120 StorageType<4, dim, Number>::base_tensor_type &data)
2121 {
2122 switch (dim)
2123 {
2124 case 1:
2125 return data[0][0];
2126
2127 case 2:
2128 // each entry of the tensor can be thought of as an entry in a
2129 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2130 // rank-2 tensors. this is the format in which we store rank-4
2131 // tensors. determine which position the present entry is
2132 // stored in
2133 {
2134 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2135 return data[base_index[indices[0]][indices[1]]]
2136 [base_index[indices[2]][indices[3]]];
2137 }
2138 case 3:
2139 // each entry of the tensor can be thought of as an entry in a
2140 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2141 // rank-2 tensors. this is the format in which we store rank-4
2142 // tensors. determine which position the present entry is
2143 // stored in
2144 {
2145 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2146 {3, 1, 5},
2147 {4, 5, 2}};
2148 return data[base_index[indices[0]][indices[1]]]
2149 [base_index[indices[2]][indices[3]]];
2150 }
2151
2152 default:
2154 }
2155
2156 // The code should never reach here.
2157 // We cannot return a static variable, as this class must support number
2158 // types that require no instances of the number type to be in scope during
2159 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2160 return data[0][0];
2161 }
2162
2163} // end of namespace internal
2164
2165
2166
2167template <int rank_, int dim, typename Number>
2168DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2170 const TableIndices<rank_> &indices)
2171{
2172 for (unsigned int r = 0; r < rank; ++r)
2173 AssertIndexRange(indices[r], dimension);
2174 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2175}
2176
2177
2178
2179template <int rank_, int dim, typename Number>
2180DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2182 const TableIndices<rank_> &indices) const
2183{
2184 for (unsigned int r = 0; r < rank; ++r)
2185 AssertIndexRange(indices[r], dimension);
2186 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2187}
2188
2189
2190
2191namespace internal
2192{
2193 namespace SymmetricTensorImplementation
2194 {
2195 template <int rank_>
2196 constexpr TableIndices<rank_>
2197 get_partially_filled_indices(const unsigned int row,
2198 const std::integral_constant<int, 2> &)
2199 {
2201 }
2202
2203
2204 template <int rank_>
2205 constexpr TableIndices<rank_>
2206 get_partially_filled_indices(const unsigned int row,
2207 const std::integral_constant<int, 4> &)
2208 {
2209 return TableIndices<rank_>(row,
2213 }
2214 } // namespace SymmetricTensorImplementation
2215} // namespace internal
2216
2217
2218template <int rank_, int dim, typename Number>
2219DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE internal::
2220 SymmetricTensorAccessors::Accessor<rank_, dim, true, rank_ - 1, Number>
2221 SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2222{
2223 return internal::SymmetricTensorAccessors::
2224 Accessor<rank_, dim, true, rank_ - 1, Number>(
2225 *this,
2226 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2227 rank_>(row, std::integral_constant<int, rank_>()));
2228}
2229
2230
2231
2232template <int rank_, int dim, typename Number>
2233DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE internal::
2234 SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2236{
2237 return internal::SymmetricTensorAccessors::
2238 Accessor<rank_, dim, false, rank_ - 1, Number>(
2239 *this,
2240 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2241 rank_>(row, std::integral_constant<int, rank_>()));
2242}
2243
2244
2245
2246template <int rank_, int dim, typename Number>
2247DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE const Number &
2249 const TableIndices<rank_> &indices) const
2250{
2251 return operator()(indices);
2252}
2253
2254
2255
2256template <int rank_, int dim, typename Number>
2257DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2259 const TableIndices<rank_> &indices)
2260{
2261 return operator()(indices);
2262}
2263
2264
2265
2266template <int rank_, int dim, typename Number>
2267DEAL_II_HOST constexpr inline const Number &
2269 const unsigned int index) const
2270{
2271 AssertIndexRange(index, n_independent_components);
2272 if constexpr (rank == 2)
2273 return data[index];
2274 else
2275 return data[decltype(data)::unrolled_to_component_indices(index)];
2276}
2277
2278
2279
2280template <int rank_, int dim, typename Number>
2281DEAL_II_HOST constexpr inline Number &
2283{
2284 AssertIndexRange(index, n_independent_components);
2285 if constexpr (rank == 2)
2286 return data[index];
2287 else
2288 return data[decltype(data)::unrolled_to_component_indices(index)];
2289}
2290
2291
2292
2293namespace internal
2294{
2295 template <int dim, typename Number>
2296 constexpr inline typename numbers::NumberTraits<Number>::real_type
2297 compute_norm(const typename SymmetricTensorAccessors::
2298 StorageType<2, dim, Number>::base_tensor_type &data)
2299 {
2300 // Make things work with AD types
2301 using std::sqrt;
2302 switch (dim)
2303 {
2304 case 1:
2306
2307 case 2:
2311
2312 case 3:
2319
2320 default:
2321 {
2322 typename numbers::NumberTraits<Number>::real_type return_value =
2324
2325 for (unsigned int d = 0; d < dim; ++d)
2326 return_value +=
2328 for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2329 return_value +=
2331
2332 return sqrt(return_value);
2333 }
2334 }
2335 }
2336
2337
2338
2339 template <int dim, typename Number>
2340 constexpr inline typename numbers::NumberTraits<Number>::real_type
2341 compute_norm(const typename SymmetricTensorAccessors::
2342 StorageType<4, dim, Number>::base_tensor_type &data)
2343 {
2344 // Make things work with AD types
2345 using std::sqrt;
2346 switch (dim)
2347 {
2348 case 1:
2350
2351 default:
2352 {
2353 typename numbers::NumberTraits<Number>::real_type return_value =
2355
2356 const unsigned int n_independent_components = data.dimension;
2357
2358 for (unsigned int i = 0; i < dim; ++i)
2359 for (unsigned int j = 0; j < dim; ++j)
2360 return_value +=
2362 for (unsigned int i = 0; i < dim; ++i)
2363 for (unsigned int j = dim; j < n_independent_components; ++j)
2364 return_value +=
2366 for (unsigned int i = dim; i < n_independent_components; ++i)
2367 for (unsigned int j = 0; j < dim; ++j)
2368 return_value +=
2370 for (unsigned int i = dim; i < n_independent_components; ++i)
2371 for (unsigned int j = dim; j < n_independent_components; ++j)
2372 return_value +=
2374
2375 return sqrt(return_value);
2376 }
2377 }
2378 }
2379
2380} // end of namespace internal
2381
2382
2383
2384template <int rank_, int dim, typename Number>
2387{
2388 return internal::compute_norm<dim, Number>(data);
2389}
2390
2391
2392
2393template <int rank_, int dim, typename Number>
2394DEAL_II_HOST constexpr unsigned int
2396 const TableIndices<rank_> &indices)
2397{
2398 return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2399 dim>(indices);
2400}
2401
2402
2403
2404namespace internal
2405{
2406 namespace SymmetricTensorImplementation
2407 {
2408 // a function to do the inverse of the unrolling from a set of
2409 // indices to a scalar index into the array in which we store
2410 // the elements of a symmetric tensor. in other words, it goes
2411 // from the scalar index into the array to a set of indices of
2412 // the tensor
2413 //
2414 // this function is for rank-2 tensors
2415 template <int dim>
2416 constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2417 unrolled_to_component_indices(const unsigned int i,
2418 const std::integral_constant<int, 2> &)
2419 {
2420 Assert(
2423 i,
2424 0,
2426 switch (dim)
2427 {
2428 case 1:
2429 {
2430 return {0, 0};
2431 }
2432
2433 case 2:
2434 {
2435 const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2436 TableIndices<2>(1, 1),
2437 TableIndices<2>(0, 1)};
2438 return table[i];
2439 }
2440
2441 case 3:
2442 {
2443 const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2444 TableIndices<2>(1, 1),
2445 TableIndices<2>(2, 2),
2446 TableIndices<2>(0, 1),
2447 TableIndices<2>(0, 2),
2448 TableIndices<2>(1, 2)};
2449 return table[i];
2450 }
2451
2452 default:
2453 if (i < dim)
2454 return {i, i};
2455
2456 for (unsigned int d = 0, c = dim; d < dim; ++d)
2457 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2458 if (c == i)
2459 return {d, e};
2460
2461 // should never get here:
2463 return {0, 0};
2464 }
2465 }
2466
2467 // a function to do the inverse of the unrolling from a set of
2468 // indices to a scalar index into the array in which we store
2469 // the elements of a symmetric tensor. in other words, it goes
2470 // from the scalar index into the array to a set of indices of
2471 // the tensor
2472 //
2473 // this function is for tensors of a rank not already handled
2474 // above
2475 template <int dim, int rank_>
2476 constexpr inline std::enable_if_t<rank_ != 2, TableIndices<rank_>>
2477 unrolled_to_component_indices(const unsigned int i,
2478 const std::integral_constant<int, rank_> &)
2479 {
2480 (void)i;
2481 Assert(
2482 (i <
2484 ExcIndexRange(i,
2485 0,
2487 n_independent_components));
2489 return TableIndices<rank_>();
2490 }
2491
2492 } // namespace SymmetricTensorImplementation
2493} // namespace internal
2494
2495template <int rank_, int dim, typename Number>
2498 const unsigned int i)
2499{
2500 return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2501 dim>(i, std::integral_constant<int, rank_>());
2502}
2503
2504
2505
2506template <int rank_, int dim, typename Number>
2507template <class Archive>
2508inline void
2509SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2510{
2511 ar &data;
2512}
2513
2514
2515#endif // DOXYGEN
2516
2517/* ----------------- Non-member functions operating on tensors. ------------ */
2518
2519
2532template <int rank_, int dim, typename Number, typename OtherNumber>
2533DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2543
2544
2557template <int rank_, int dim, typename Number, typename OtherNumber>
2558DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2568
2569
2577template <int rank_, int dim, typename Number, typename OtherNumber>
2582{
2583 return Tensor<rank_, dim, Number>(left) + right;
2584}
2585
2586
2594template <int rank_, int dim, typename Number, typename OtherNumber>
2599{
2600 return left + Tensor<rank_, dim, OtherNumber>(right);
2601}
2602
2603
2611template <int rank_, int dim, typename Number, typename OtherNumber>
2616{
2617 return Tensor<rank_, dim, Number>(left) - right;
2618}
2619
2620
2628template <int rank_, int dim, typename Number, typename OtherNumber>
2633{
2634 return left - Tensor<rank_, dim, OtherNumber>(right);
2635}
2636
2637
2638
2639template <int dim, typename Number>
2642{
2643 switch (dim)
2644 {
2645 case 1:
2646 return t.data[0];
2647 case 2:
2648 return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2649 case 3:
2650 {
2651 // in analogy to general tensors, but
2652 // there's something to be simplified for
2653 // the present case
2654 const Number tmp = t.data[3] * t.data[4] * t.data[5];
2655 return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2656 t.data[0] * t.data[5] * t.data[5] -
2657 t.data[1] * t.data[4] * t.data[4] -
2658 t.data[2] * t.data[3] * t.data[3]);
2659 }
2660 default:
2663 }
2664}
2665
2666
2667
2679template <int dim, typename Number>
2685
2686
2687
2688template <int dim, typename Number>
2689DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
2691{
2692 Number t = d.data[0];
2693 for (unsigned int i = 1; i < dim; ++i)
2694 t += d.data[i];
2695 return t;
2696}
2697
2698
2710template <int dim, typename Number>
2711DEAL_II_HOST constexpr Number
2713{
2714 return trace(t);
2715}
2716
2717
2729template <typename Number>
2730DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2735
2736
2737
2756template <typename Number>
2757DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2759{
2760 return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2761}
2762
2763
2764
2773template <typename Number>
2774DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2776{
2777 return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2778 t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2779}
2780
2781
2782
2790template <typename Number>
2791std::array<Number, 1>
2793
2794
2795
2818template <typename Number>
2819std::array<Number, 2>
2821
2822
2823
2846template <typename Number>
2847std::array<Number, 3>
2849
2850
2851
2852namespace internal
2853{
2854 namespace SymmetricTensorImplementation
2855 {
2867 template <int dim, typename Number>
2868 void
2869 tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2871 std::array<Number, dim> &d,
2872 std::array<Number, dim - 1> &e);
2873
2874
2875
2889 template <int dim, typename Number>
2890 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2891 ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2892
2893
2894
2908 template <int dim, typename Number>
2909 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2911
2912
2913
2927 template <typename Number>
2928 std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
2929 hybrid(const ::SymmetricTensor<2, 2, Number> &A);
2930
2931
2932
2947 template <typename Number>
2948 std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
2949 hybrid(const ::SymmetricTensor<2, 3, Number> &A);
2950
2955 template <int dim, typename Number>
2957 {
2958 using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
2959 bool
2960 operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
2961 {
2962 return lhs.first > rhs.first;
2963 }
2964 };
2965
2966 } // namespace SymmetricTensorImplementation
2967
2968} // namespace internal
2969
2970
2971
2972// The line below is to ensure that doxygen puts the full description
2973// of this global enumeration into the documentation
2974// See https://stackoverflow.com/a/1717984
2986{
2996 hybrid,
3014 jacobi
3015};
3016
3017
3018
3029template <int dim, typename Number>
3030std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3034
3035
3036
3045template <int rank_, int dim, typename Number>
3048{
3049 return t;
3050}
3051
3052
3053
3054template <int dim, typename Number>
3055DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3058{
3060
3061 // subtract scaled trace from the diagonal
3062 const Number tr = trace(t) * internal::NumberType<Number>::value(1.0 / dim);
3063 for (unsigned int i = 0; i < dim; ++i)
3064 tmp.data[i] -= tr;
3065
3066 return tmp;
3067}
3068
3069
3070
3071template <int dim, typename Number>
3075{
3076 // create a default constructed matrix filled with
3077 // zeros, then set the diagonal elements to one
3079 switch (dim)
3080 {
3081 case 1:
3083 break;
3084 case 2:
3085 tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3086 break;
3087 case 3:
3088 tmp.data[0] = tmp.data[1] = tmp.data[2] =
3090 break;
3091 default:
3092 for (unsigned int d = 0; d < dim; ++d)
3094 }
3095 return tmp;
3096}
3097
3098
3099
3100template <int dim, typename Number>
3103{
3105
3106 // fill the elements treating the diagonal
3107 for (unsigned int i = 0; i < dim; ++i)
3108 for (unsigned int j = 0; j < dim; ++j)
3109 tmp.data[i][j] =
3110 internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3111
3112 // then fill the ones that copy over the
3113 // non-diagonal elements. note that during
3114 // the double-contraction, we handle the
3115 // off-diagonal elements twice, so simply
3116 // copying requires a weight of 1/2
3117 for (unsigned int i = dim;
3118 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3119 n_rank2_components;
3120 ++i)
3122
3123 return tmp;
3124}
3125
3126
3127
3128template <int dim, typename Number>
3132{
3134
3135 // fill the elements treating the diagonal
3136 for (unsigned int i = 0; i < dim; ++i)
3138
3139 // then fill the ones that copy over the
3140 // non-diagonal elements. note that during
3141 // the double-contraction, we handle the
3142 // off-diagonal elements twice, so simply
3143 // copying requires a weight of 1/2
3144 for (unsigned int i = dim;
3145 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3146 n_rank2_components;
3147 ++i)
3149
3150 return tmp;
3151}
3152
3153
3154
3164template <int dim, typename Number>
3171
3172
3173
3184template <int dim, typename Number>
3191
3192
3193
3215template <int dim, typename Number>
3219{
3221
3222 // fill only the elements really needed
3223 for (unsigned int i = 0; i < dim; ++i)
3224 for (unsigned int j = i; j < dim; ++j)
3225 for (unsigned int k = 0; k < dim; ++k)
3226 for (unsigned int l = k; l < dim; ++l)
3227 tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3228
3229 return tmp;
3230}
3231
3257template <int dim, typename Number>
3258std::pair<SymmetricTensor<2, dim, Number>, SymmetricTensor<2, dim, Number>>
3260{
3261 Assert(dim <= 3, ExcNotImplemented());
3262
3263 const std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3264 eigen_system = eigenvectors(original_tensor);
3265
3266 std::pair<SymmetricTensor<2, dim, Number>, SymmetricTensor<2, dim, Number>>
3267 positive_negative_tensors;
3268
3269 auto &[positive_part_tensor, negative_part_tensor] =
3270 positive_negative_tensors;
3271
3272 positive_part_tensor = 0;
3273 for (unsigned int i = 0; i < dim; ++i)
3274 if (eigen_system[i].first > 0)
3275 positive_part_tensor += eigen_system[i].first *
3276 symmetrize(outer_product(eigen_system[i].second,
3277 eigen_system[i].second));
3278
3279 negative_part_tensor = 0;
3280 for (unsigned int i = 0; i < dim; ++i)
3281 if (eigen_system[i].first < 0)
3282 negative_part_tensor += eigen_system[i].first *
3283 symmetrize(outer_product(eigen_system[i].second,
3284 eigen_system[i].second));
3285
3286 return positive_negative_tensors;
3287}
3288
3321template <int dim, typename Number>
3322std::tuple<SymmetricTensor<2, dim, Number>,
3327 const SymmetricTensor<2, dim, Number> &original_tensor)
3328{
3329 Assert(dim <= 3, ExcNotImplemented());
3330
3331 auto heaviside_function{[](const double x) {
3332 if (std::fabs(x) < 1.0e-16)
3333 return 0.5;
3334 if (x > 0)
3335 return 1.0;
3336 else
3337 return 0.0;
3338 }};
3339
3340 std::tuple<SymmetricTensor<2, dim, Number>,
3344 positive_negative_tensors_projectors;
3345
3346 auto &[positive_part_tensor,
3347 negative_part_tensor,
3348 positive_projector,
3349 negative_projector] = positive_negative_tensors_projectors;
3350
3351 const std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3352 eigen_system = eigenvectors(original_tensor);
3353
3354 positive_part_tensor = 0;
3355 for (unsigned int i = 0; i < dim; ++i)
3356 if (eigen_system[i].first > 0)
3357 positive_part_tensor += eigen_system[i].first *
3358 symmetrize(outer_product(eigen_system[i].second,
3359 eigen_system[i].second));
3360
3361 negative_part_tensor = 0;
3362 for (unsigned int i = 0; i < dim; ++i)
3363 if (eigen_system[i].first < 0)
3364 negative_part_tensor += eigen_system[i].first *
3365 symmetrize(outer_product(eigen_system[i].second,
3366 eigen_system[i].second));
3367
3368 std::array<SymmetricTensor<2, dim, Number>, dim> M;
3369 for (unsigned int a = 0; a < dim; ++a)
3370 M[a] =
3371 symmetrize(outer_product(eigen_system[a].second, eigen_system[a].second));
3372
3373 std::array<SymmetricTensor<4, dim, Number>, dim> Q;
3374 for (unsigned int a = 0; a < dim; ++a)
3375 Q[a] = outer_product(M[a], M[a]);
3376
3377 std::array<std::array<SymmetricTensor<4, dim, Number>, dim>, dim> G;
3378 for (unsigned int a = 0; a < dim; ++a)
3379 for (unsigned int b = 0; b < dim; ++b)
3380 for (unsigned int i = 0; i < dim; ++i)
3381 for (unsigned int j = 0; j < dim; ++j)
3382 for (unsigned int k = 0; k < dim; ++k)
3383 for (unsigned int l = 0; l < dim; ++l)
3384 G[a][b][i][j][k][l] =
3385 M[a][i][k] * M[b][j][l] + M[a][i][l] * M[b][j][k];
3386
3387 // positive P
3388 positive_projector = 0;
3389 for (unsigned int a = 0; a < dim; ++a)
3390 {
3391 double lambda_a = eigen_system[a].first;
3392 positive_projector += heaviside_function(lambda_a) * Q[a];
3393 for (unsigned int b = 0; b < dim; ++b)
3394 {
3395 if (b != a)
3396 {
3397 double lambda_b = eigen_system[b].first;
3398
3399 double v_ab = 0.0;
3400 if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
3401 v_ab = (std::fmax(lambda_a, 0.0) - std::fmax(lambda_b, 0.0)) /
3402 (lambda_a - lambda_b);
3403 else
3404 v_ab = 0.5 * (heaviside_function(lambda_a) +
3405 heaviside_function(lambda_b));
3406
3407 positive_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
3408 }
3409 }
3410 }
3411
3412 // negative P
3413 negative_projector = 0;
3414 for (unsigned int a = 0; a < dim; ++a)
3415 {
3416 double lambda_a = eigen_system[a].first;
3417 negative_projector += heaviside_function(-lambda_a) * Q[a];
3418 for (unsigned int b = 0; b < dim; ++b)
3419 {
3420 if (b != a)
3421 {
3422 double lambda_b = eigen_system[b].first;
3423
3424 double v_ab = 0.0;
3425 if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
3426 v_ab = (std::fmin(lambda_a, 0.0) - std::fmin(lambda_b, 0.0)) /
3427 (lambda_a - lambda_b);
3428 else
3429 v_ab = 0.5 * (heaviside_function(-lambda_a) +
3430 heaviside_function(-lambda_b));
3431
3432 negative_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
3433 }
3434 }
3435 }
3436
3437 return positive_negative_tensors_projectors;
3438}
3439
3447template <int dim, typename Number>
3448DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3451{
3453 for (unsigned int d = 0; d < dim; ++d)
3454 result[d][d] = t[d][d];
3455
3456 const Number half = internal::NumberType<Number>::value(0.5);
3457 for (unsigned int d = 0; d < dim; ++d)
3458 for (unsigned int e = d + 1; e < dim; ++e)
3459 result[d][e] = (t[d][e] + t[e][d]) * half;
3460 return result;
3461}
3462
3463
3464
3476template <int dim, typename Number>
3477DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3479 symmetrize(const Tensor<4, dim, Number> &t, const bool major_symmetry)
3480{
3482
3483 const Number half = internal::NumberType<Number>::value(0.5);
3484
3485 // minor symmetry - A_{ijkl}=A_{jikl}=A_{ijlk}=A_{jilk}
3486 for (unsigned int i = 0; i < dim; ++i)
3487 for (unsigned int j = 0; j < dim; ++j)
3488 for (unsigned int k = 0; k < dim; ++k)
3489 for (unsigned int l = 0; l < dim; ++l)
3490 {
3491 if (i != j && k == l)
3492 {
3493 // A_{ijkk}=A_{jikk}
3494 result[i][j][k][k] = (t[i][j][k][k] + t[j][i][k][k]) * half;
3495 }
3496 else if (i == j && k != l)
3497 {
3498 // A_{iikl}=A_{iilk}
3499 result[i][i][k][l] = (t[i][i][k][l] + t[i][i][l][k]) * half;
3500 }
3501 else if (i != j && k != l)
3502 {
3503 // A_{ijkl}=A_{jilk}
3504 result[i][j][k][l] = (t[i][j][k][l] + t[j][i][k][l] +
3505 t[i][j][l][k] + t[j][i][l][k]) *
3506 half * half;
3507 }
3508 else
3509 {
3510 // A_{iijj} and A_{iiii} unchanged
3511 result[i][j][k][l] = t[i][j][k][l];
3512 }
3513 }
3514
3515 // in case major symmetry is also required
3516 if (major_symmetry)
3517 {
3518 // major symmetry - A_{ijkl}=A_{klij}
3519 for (unsigned int i = 0; i < dim; ++i)
3520 for (unsigned int j = i; j < dim; ++j)
3521 for (unsigned int k = 0; k < dim; ++k)
3522 for (unsigned int l = k; l < dim; ++l)
3523 result[i][j][k][l] = (t[i][j][k][l] + t[k][l][i][j]) * half;
3524 }
3525 return result;
3526}
3527
3528
3529
3537template <int rank_, int dim, typename Number>
3538DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3540 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3541{
3543 tt *= factor;
3544 return tt;
3545}
3546
3547
3548
3556template <int rank_, int dim, typename Number>
3558operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3559{
3560 // simply forward to the other operator
3561 return t * factor;
3562}
3563
3564
3565
3590template <int rank_, int dim, typename Number, typename OtherNumber>
3592 rank_,
3593 dim,
3594 typename ProductType<Number,
3595 typename EnableIfScalar<OtherNumber>::type>::type>
3597 const OtherNumber &factor)
3598{
3599 // form the product. we have to convert the two factors into the final
3600 // type via explicit casts because, for awkward reasons, the C++
3601 // standard committee saw it fit to not define an
3602 // operator*(float,std::complex<double>)
3603 // (as well as with switched arguments and double<->float).
3604 using product_type = typename ProductType<Number, OtherNumber>::type;
3607 return tt;
3608}
3609
3610
3611
3619template <int rank_, int dim, typename Number, typename OtherNumber>
3621 rank_,
3622 dim,
3623 typename ProductType<OtherNumber,
3624 typename EnableIfScalar<Number>::type>::type>
3625operator*(const Number &factor,
3627{
3628 // simply forward to the other operator with switched arguments
3629 return (t * factor);
3630}
3631
3632
3633
3639template <int rank_, int dim, typename Number, typename OtherNumber>
3640DEAL_II_HOST constexpr inline SymmetricTensor<
3641 rank_,
3642 dim,
3643 typename ProductType<Number,
3644 typename EnableIfScalar<OtherNumber>::type>::type>
3646 const OtherNumber &factor)
3647{
3648 using product_type = typename ProductType<Number, OtherNumber>::type;
3651 return tt;
3652}
3653
3654
3655
3662template <int rank_, int dim>
3664operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3665{
3667 tt *= factor;
3668 return tt;
3669}
3670
3671
3672
3679template <int rank_, int dim>
3681operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3682{
3684 tt *= factor;
3685 return tt;
3686}
3687
3688
3689
3695template <int rank_, int dim>
3697operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3698{
3700 tt /= factor;
3701 return tt;
3702}
3703
3713template <int dim, typename Number, typename OtherNumber>
3718{
3719 return (t1 * t2);
3720}
3721
3722
3735template <int dim, typename Number, typename OtherNumber>
3736DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3740{
3742 typename ProductType<Number, OtherNumber>::type>::value(0.0);
3743 for (unsigned int i = 0; i < dim; ++i)
3744 for (unsigned int j = 0; j < dim; ++j)
3745 s += t1[i][j] * t2[i][j];
3746 return s;
3747}
3748
3749
3762template <int dim, typename Number, typename OtherNumber>
3767{
3768 return scalar_product(t2, t1);
3769}
3770
3771
3786template <typename Number, typename OtherNumber>
3787DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
3792{
3793 tmp[0][0] = t[0][0][0][0] * s[0][0];
3794}
3795
3796
3797
3812template <typename Number, typename OtherNumber>
3813DEAL_II_HOST constexpr inline void
3818{
3819 tmp[0][0] = t[0][0][0][0] * s[0][0];
3820}
3821
3822
3823
3838template <typename Number, typename OtherNumber>
3839DEAL_II_HOST constexpr inline void
3844{
3845 const unsigned int dim = 2;
3846
3847 for (unsigned int i = 0; i < dim; ++i)
3848 for (unsigned int j = i; j < dim; ++j)
3849 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3850 2 * t[i][j][0][1] * s[0][1];
3851}
3852
3853
3854
3869template <typename Number, typename OtherNumber>
3870DEAL_II_HOST constexpr inline void
3875{
3876 const unsigned int dim = 2;
3877
3878 for (unsigned int i = 0; i < dim; ++i)
3879 for (unsigned int j = i; j < dim; ++j)
3880 tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3881 2 * s[0][1] * t[0][1][i][j];
3882}
3883
3884
3885
3900template <typename Number, typename OtherNumber>
3901DEAL_II_HOST constexpr inline void
3906{
3907 const unsigned int dim = 3;
3908
3909 for (unsigned int i = 0; i < dim; ++i)
3910 for (unsigned int j = i; j < dim; ++j)
3911 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3912 t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3913 2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3914}
3915
3916
3917
3932template <typename Number, typename OtherNumber>
3933DEAL_II_HOST constexpr inline void
3938{
3939 const unsigned int dim = 3;
3940
3941 for (unsigned int i = 0; i < dim; ++i)
3942 for (unsigned int j = i; j < dim; ++j)
3943 tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3944 s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3945 2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3946}
3947
3948
3949
3956template <int dim, typename Number, typename OtherNumber>
3957DEAL_II_HOST constexpr Tensor<1,
3958 dim,
3961 const Tensor<1, dim, OtherNumber> &src2)
3962{
3964 for (unsigned int i = 0; i < dim; ++i)
3965 {
3966 dest[i] = src1[i][0] * src2[0];
3967 for (unsigned int j = 1; j < dim; ++j)
3968 dest[i] += src1[i][j] * src2[j];
3969 }
3970 return dest;
3971}
3972
3973
3980template <int dim, typename Number, typename OtherNumber>
3981DEAL_II_HOST constexpr Tensor<1,
3982 dim,
3986{
3987 // this is easy for symmetric tensors:
3988 return src2 * src1;
3989}
3990
3991
3992
4012template <int rank_1,
4013 int rank_2,
4014 int dim,
4015 typename Number,
4016 typename OtherNumber>
4018 typename Tensor<rank_1 + rank_2 - 2,
4019 dim,
4020 typename ProductType<Number, OtherNumber>::type>::tensor_type
4026
4027
4028
4048template <int rank_1,
4049 int rank_2,
4050 int dim,
4051 typename Number,
4052 typename OtherNumber>
4054 typename Tensor<rank_1 + rank_2 - 2,
4055 dim,
4056 typename ProductType<Number, OtherNumber>::type>::tensor_type
4059{
4060 return Tensor<rank_1, dim, Number>(src1) * src2;
4061}
4062
4063
4064
4074template <int dim, typename Number>
4075inline std::ostream &
4076operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4077{
4078 // make our lives a bit simpler by outputting
4079 // the tensor through the operator for the
4080 // general Tensor class
4082
4083 for (unsigned int i = 0; i < dim; ++i)
4084 for (unsigned int j = 0; j < dim; ++j)
4085 tt[i][j] = t[i][j];
4086
4087 return out << tt;
4088}
4089
4090
4091
4101template <int dim, typename Number>
4102inline std::ostream &
4103operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4104{
4105 // make our lives a bit simpler by outputting
4106 // the tensor through the operator for the
4107 // general Tensor class
4109
4110 for (unsigned int i = 0; i < dim; ++i)
4111 for (unsigned int j = 0; j < dim; ++j)
4112 for (unsigned int k = 0; k < dim; ++k)
4113 for (unsigned int l = 0; l < dim; ++l)
4114 tt[i][j][k][l] = t[i][j][k][l];
4115
4116 return out << tt;
4117}
4118
4119
4121
4122#endif
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
DEAL_II_HOST constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, false, rank_ - 1, Number > operator[](const unsigned int row)
std::pair< SymmetricTensor< 2, dim, Number >, SymmetricTensor< 2, dim, Number > > positive_negative_split(const SymmetricTensor< 2, dim, Number > &original_tensor)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, OtherNumber > &s)
static DEAL_II_HOST constexpr std::size_t memory_consumption()
DEAL_II_HOST constexpr SymmetricTensor & operator/=(const OtherNumber &factor)
DEAL_II_HOST constexpr Tensor< rank_1+rank_2-2, dim, typenameProductType< Number, OtherNumber >::type >::tensor_type operator*(const Tensor< rank_1, dim, Number > &src1, const SymmetricTensor< rank_2, dim, OtherNumber > &src2)
DEAL_II_HOST constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const SymmetricTensor< 2, dim, Number > &src1, const Tensor< 1, dim, OtherNumber > &src2)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, typename ProductType< OtherNumber, typename EnableIfScalar< Number >::type >::type > operator*(const Number &factor, const SymmetricTensor< rank_, dim, OtherNumber > &t)
DEAL_II_HOST constexpr Number & operator[](const TableIndices< rank_ > &indices)
std::array< Number, 2 > eigenvalues(const SymmetricTensor< 2, 2, Number > &T)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
void serialize(Archive &ar, const unsigned int version)
DEAL_II_HOST constexpr Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::tuple< SymmetricTensor< 2, dim, Number >, SymmetricTensor< 2, dim, Number >, SymmetricTensor< 4, dim, Number >, SymmetricTensor< 4, dim, Number > > positive_negative_projectors(const SymmetricTensor< 2, dim, Number > &original_tensor)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
DEAL_II_HOST constexpr const Number & operator()(const TableIndices< rank_ > &indices) const
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
typename base_tensor_descriptor::base_tensor_type base_tensor_type
DEAL_II_HOST constexpr SymmetricTensor(const SymmetricTensor< rank_, dim, OtherNumber > &initializer)
DEAL_II_HOST constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
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 SymmetricTensor< rank_, dim > operator/(const SymmetricTensor< rank_, dim > &t, const double factor)
DEAL_II_HOST constexpr bool operator==(const SymmetricTensor &) const
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
DEAL_II_HOST constexpr Number & operator()(const TableIndices< rank_ > &indices)
DEAL_II_HOST constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
DEAL_II_HOST constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
DEAL_II_HOST constexpr SymmetricTensor(const Number(&array)[n_independent_components])
DEAL_II_HOST constexpr bool operator!=(const SymmetricTensor &) const
DEAL_II_HOST constexpr SymmetricTensor & operator=(const Number &d)
DEAL_II_HOST constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
DEAL_II_HOST constexpr SymmetricTensor operator-() const
DEAL_II_HOST constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, OtherNumber > &t2)
base_tensor_type data
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
DEAL_II_HOST constexpr Number & access_raw_entry(const unsigned int unrolled_index)
DEAL_II_HOST constexpr Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const Tensor< 1, dim, Number > &src1, const SymmetricTensor< 2, dim, OtherNumber > &src2)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim > operator*(const SymmetricTensor< rank_, dim > &t, const double factor)
SymmetricTensor(const Tensor< 2, dim, OtherNumber > &t)
DEAL_II_HOST constexpr void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 1, Number > &s, const SymmetricTensor< 4, 1, OtherNumber > &t)
DEAL_II_HOST constexpr Tensor< rank_1+rank_2-2, dim, typenameProductType< Number, OtherNumber >::type >::tensor_type operator*(const SymmetricTensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
DEAL_II_HOST constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
static DEAL_II_HOST constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
DEAL_II_HOST constexpr void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 3, Number > &t, const SymmetricTensor< 2, 3, OtherNumber > &s)
DEAL_II_HOST constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
DEAL_II_HOST constexpr void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
static DEAL_II_HOST constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
DEAL_II_HOST constexpr SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
DEAL_II_HOST constexpr void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, OtherNumber > &t)
DEAL_II_HOST constexpr void clear()
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
DEAL_II_HOST constexpr void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 2, Number > &s, const SymmetricTensor< 4, 2, OtherNumber > &t)
DEAL_II_HOST constexpr Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim, Number > operator*(const Number &factor, const SymmetricTensor< rank_, dim, Number > &t)
DEAL_II_HOST constexpr numbers::NumberTraits< Number >::real_type norm() const
DEAL_II_HOST constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 4, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 4, dim, OtherNumber > &s) const
DEAL_II_HOST constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
DEAL_II_HOST constexpr SymmetricTensor()=default
DEAL_II_HOST constexpr Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
DEAL_II_HOST constexpr Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > symmetrize(const Tensor< 4, dim, Number > &t, const bool major_symmetry)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
DEAL_II_HOST constexpr SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
DEAL_II_HOST constexpr const Number & operator[](const TableIndices< rank_ > &indices) const
DEAL_II_HOST constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
DEAL_II_HOST constexpr SymmetricTensor< rank_, dim > operator*(const double factor, const SymmetricTensor< rank_, dim > &t)
std::array< Number, 3 > eigenvalues(const SymmetricTensor< 2, 3, Number > &T)
DEAL_II_HOST constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
typename AccessorTypes< rank, dim, constness, Number >::reference reference
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
DEAL_II_HOST constexpr Accessor(tensor_type &tensor, const TableIndices< rank > &previous_indices)
DEAL_II_HOST constexpr reference operator[](const unsigned int)
DEAL_II_HOST constexpr reference operator[](const unsigned int) const
DEAL_II_HOST constexpr Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i) const
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
DEAL_II_HOST constexpr Accessor(tensor_type &tensor, const TableIndices< rank > &previous_indices)
DEAL_II_HOST constexpr Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i)
typename AccessorTypes< rank, dim, constness, Number >::reference reference
DEAL_II_HOST constexpr Accessor(const Accessor &)=default
#define DEAL_II_ALWAYS_INLINE
Definition config.h:161
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_CONSTEXPR
Definition config.h:268
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_HOST
Definition config.h:176
std::ostream & operator<<(std::ostream &out, const DerivativeForm< order, dim, spacedim, Number > &df)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
Definition mpi.cc:746
constexpr char N
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A, ::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim - 1 > &e)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:896
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static constexpr const T & value(const T &t)
Definition numbers.h:683
typename ProductType< Number, OtherNumber >::type value_type
std::pair< Number, Tensor< 1, dim, Number > > EigValsVecs
bool operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
static real_type abs(const number &x)
Definition numbers.h:574
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
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 > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
SymmetricTensorEigenvectorMethod
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()