Reference documentation for deal.II version GIT relicensing-822-gf1642d4ea5 2024-06-12 14:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
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
25#include <deal.II/base/tensor.h>
26
27#include <array>
28
30
31// Forward declaration
32#ifndef DOXYGEN
33template <int rank, int dim, typename Number = double>
34class SymmetricTensor;
35#endif
36
47template <int dim, typename Number = double>
51
80template <int dim, typename Number = double>
84
122template <int dim, typename Number = double>
126
127template <int dim, typename Number>
130
131template <int dim, typename Number>
134
144template <int dim2, typename Number>
145DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
147
158template <int dim, typename Number>
159DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
162
176template <int dim, typename Number>
179
180
181
182namespace internal
183{
184 // Workaround: The following 4 overloads are necessary to be able to
185 // compile the library with Apple Clang 8 and older. We should remove
186 // these overloads again when we bump the minimal required version to
187 // something later than clang-3.6 / Apple Clang 6.3.
188 template <int rank, int dim, typename T, typename U>
189 struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
190 {
191 using type =
192 SymmetricTensor<rank,
193 dim,
194 std::complex<typename ProductType<T, U>::type>>;
195 };
196
197 template <int rank, int dim, typename T, typename U>
198 struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
199 std::complex<U>>
200 {
201 using type =
202 SymmetricTensor<rank,
203 dim,
204 std::complex<typename ProductType<T, U>::type>>;
205 };
206
207 template <typename T, int rank, int dim, typename U>
208 struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
209 {
210 using type =
211 SymmetricTensor<rank,
212 dim,
213 std::complex<typename ProductType<T, U>::type>>;
214 };
215
216 template <int rank, int dim, typename T, typename U>
217 struct ProductTypeImpl<std::complex<T>,
218 SymmetricTensor<rank, dim, std::complex<U>>>
219 {
220 using type =
221 SymmetricTensor<rank,
222 dim,
223 std::complex<typename ProductType<T, U>::type>>;
224 };
225 // end workaround
226
231 namespace SymmetricTensorImplementation
232 {
237 template <int rank, int dim, typename Number>
238 struct Inverse;
239 } // namespace SymmetricTensorImplementation
240
245 namespace SymmetricTensorAccessors
246 {
255 merge(const TableIndices<2> &previous_indices,
256 const unsigned int new_index,
257 const unsigned int position)
258 {
259 AssertIndexRange(position, 2);
260
261 if (position == 0)
263 else
264 return {previous_indices[0], new_index};
265 }
266
267
268
277 merge(const TableIndices<4> &previous_indices,
278 const unsigned int new_index,
279 const unsigned int position)
280 {
281 AssertIndexRange(position, 4);
282
283 switch (position)
284 {
285 case 0:
286 return {new_index,
290 case 1:
291 return {previous_indices[0],
292 new_index,
295 case 2:
296 return {previous_indices[0],
297 previous_indices[1],
298 new_index,
300 case 3:
301 return {previous_indices[0],
302 previous_indices[1],
303 previous_indices[2],
304 new_index};
305 default:
307 return {};
308 }
309 }
310
311
318 template <int rank1,
319 int rank2,
320 int dim,
321 typename Number,
322 typename OtherNumber = Number>
329
330
337 template <int dim, typename Number, typename OtherNumber>
338 struct double_contraction_result<2, 2, dim, Number, OtherNumber>
339 {
341 };
342
343
344
357 template <int rank, int dim, typename Number>
359
363 template <int dim, typename Number>
364 struct StorageType<2, dim, Number>
365 {
370 static const unsigned int n_independent_components =
371 (dim * dim + dim) / 2;
372
377 };
378
379
380
384 template <int dim, typename Number>
385 struct StorageType<4, dim, Number>
386 {
392 static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
393
397 static const unsigned int n_independent_components =
398 (n_rank2_components *
400
408 };
409
410
411
416 template <int rank, int dim, bool constness, typename Number>
418
425 template <int rank, int dim, typename Number>
426 struct AccessorTypes<rank, dim, true, Number>
427 {
428 using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
429
430 using reference = const Number &;
431 };
432
439 template <int rank, int dim, typename Number>
440 struct AccessorTypes<rank, dim, false, Number>
441 {
443
444 using reference = Number &;
445 };
446
447
480 template <int rank, int dim, bool constness, int P, typename Number>
482 {
483 public:
487 using reference =
491
492 private:
514
519 constexpr DEAL_II_ALWAYS_INLINE
520 Accessor(const Accessor &) = default;
521
522 public:
527 constexpr Accessor<rank, dim, constness, P - 1, Number>
528 operator[](const unsigned int i);
529
534 constexpr Accessor<rank, dim, constness, P - 1, Number>
535 operator[](const unsigned int i) const;
536
537 private:
543
544 // Declare some other classes as friends. Make sure to work around bugs
545 // in some compilers:
546 template <int, int, typename>
547 friend class ::SymmetricTensor;
548 template <int, int, bool, int, typename>
549 friend class Accessor;
550 friend class ::SymmetricTensor<rank, dim, Number>;
551 friend class Accessor<rank, dim, constness, P + 1, Number>;
552 };
553
554
555
563 template <int rank, int dim, bool constness, typename Number>
564 class Accessor<rank, dim, constness, 1, Number>
565 {
566 public:
570 using reference =
574
575 private:
600
605 constexpr DEAL_II_ALWAYS_INLINE
606 Accessor(const Accessor &) = default;
607
608 public:
613 constexpr reference
614 operator[](const unsigned int);
615
620 constexpr reference
621 operator[](const unsigned int) const;
622
623 private:
629
630 // Declare some other classes as friends. Make sure to work around bugs
631 // in some compilers:
632 template <int, int, typename>
633 friend class ::SymmetricTensor;
634 template <int, int, bool, int, typename>
636 friend class ::SymmetricTensor<rank, dim, Number>;
637 friend class SymmetricTensorAccessors::
638 Accessor<rank, dim, constness, 2, Number>;
639 };
640 } // namespace SymmetricTensorAccessors
641} // namespace internal
642
643
644
717template <int rank_, int dim, typename Number>
719{
720public:
721 static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
722
731 static constexpr unsigned int dimension = dim;
732
736 static const unsigned int rank = rank_;
737
743 static constexpr unsigned int n_independent_components =
746
751 constexpr DEAL_II_ALWAYS_INLINE
752 SymmetricTensor() = default;
753
767 template <typename OtherNumber>
769
786 constexpr SymmetricTensor(const Number (&array)[n_independent_components]);
787
793 template <typename OtherNumber>
794 DEAL_II_HOST constexpr explicit SymmetricTensor(
796
806 Number *
808
818 const Number *
819 begin_raw() const;
820
830 Number *
832
843 const Number *
844 end_raw() const;
845
852 template <typename OtherNumber>
853 DEAL_II_HOST constexpr SymmetricTensor &
855
863 constexpr SymmetricTensor &
864 operator=(const Number &d);
865
871 constexpr operator Tensor<rank_, dim, Number>() const;
872
877 constexpr bool
879
884 constexpr bool
886
890 template <typename OtherNumber>
891 DEAL_II_HOST constexpr SymmetricTensor &
893
897 template <typename OtherNumber>
898 DEAL_II_HOST constexpr SymmetricTensor &
900
905 template <typename OtherNumber>
906 DEAL_II_HOST constexpr SymmetricTensor &
907 operator*=(const OtherNumber &factor);
908
912 template <typename OtherNumber>
913 DEAL_II_HOST constexpr SymmetricTensor &
914 operator/=(const OtherNumber &factor);
915
920 constexpr SymmetricTensor
921 operator-() const;
922
975 template <typename OtherNumber>
976 DEAL_II_HOST DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
977 double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
979
984 template <typename OtherNumber>
985 DEAL_II_HOST DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
986 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
988
993 constexpr Number &
995
1000 constexpr const Number &
1001 operator()(const TableIndices<rank_> &indices) const;
1002
1008 constexpr internal::SymmetricTensorAccessors::
1009 Accessor<rank_, dim, true, rank_ - 1, Number>
1010 operator[](const unsigned int row) const;
1011
1017 constexpr internal::SymmetricTensorAccessors::
1018 Accessor<rank_, dim, false, rank_ - 1, Number>
1019 operator[](const unsigned int row);
1020
1027 constexpr const Number &
1028 operator[](const TableIndices<rank_> &indices) const;
1029
1036 constexpr Number &
1038
1046 constexpr const Number &
1047 access_raw_entry(const unsigned int unrolled_index) const;
1048
1056 constexpr Number &
1058
1070 norm() const;
1071
1079 static DEAL_II_HOST constexpr unsigned int
1081
1087 static DEAL_II_HOST constexpr TableIndices<rank_>
1088 unrolled_to_component_indices(const unsigned int i);
1089
1103 constexpr void
1105
1110 static DEAL_II_HOST constexpr std::size_t
1112
1118 template <class Archive>
1119 void
1120 serialize(Archive &ar, const unsigned int version);
1121
1122private:
1128
1132 using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
1133
1138
1139#ifndef DOXYGEN
1140
1141 // Make all other symmetric tensors friends.
1142 template <int, int, typename>
1143 friend class SymmetricTensor;
1144
1145 // Make a few more functions friends.
1146 template <int dim2, typename Number2>
1147 friend DEAL_II_HOST constexpr Number2
1149
1150 template <int dim2, typename Number2>
1151 friend DEAL_II_HOST DEAL_II_CONSTEXPR Number2
1153
1154 template <int dim2, typename Number2>
1157
1158 template <int dim2, typename Number2>
1161
1162 template <int dim2, typename Number2>
1165
1166 template <int dim2, typename Number2>
1169
1170
1171 // Make a few helper classes friends as well.
1173 Inverse<2, dim, Number>;
1174
1176 Inverse<4, dim, Number>;
1177#endif
1178};
1179
1180
1181
1182// ------------------------- inline functions ------------------------
1183
1184#ifndef DOXYGEN
1185
1186// provide declarations for static members
1187template <int rank, int dim, typename Number>
1188const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1189
1190template <int rank_, int dim, typename Number>
1191constexpr unsigned int
1192 SymmetricTensor<rank_, dim, Number>::n_independent_components;
1193
1195{
1196 namespace SymmetricTensorAccessors
1197 {
1198 template <int rank_, int dim, bool constness, int P, typename Number>
1200 Accessor<rank_, dim, constness, P, Number>::Accessor(
1201 tensor_type &tensor,
1202 const TableIndices<rank_> &previous_indices)
1203 : tensor(tensor)
1204 , previous_indices(previous_indices)
1205 {}
1206
1207
1208
1209 template <int rank_, int dim, bool constness, int P, typename Number>
1210 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1211 Accessor<rank_, dim, constness, P - 1, Number>
1212 Accessor<rank_, dim, constness, P, Number>::operator[](
1213 const unsigned int i)
1214 {
1215 return Accessor<rank_, dim, constness, P - 1, Number>(
1216 tensor, merge(previous_indices, i, rank_ - P));
1217 }
1218
1219
1220
1221 template <int rank_, int dim, bool constness, int P, typename Number>
1223 Accessor<rank_, dim, constness, P - 1, Number>
1224 Accessor<rank_, dim, constness, P, Number>::operator[](
1225 const unsigned int i) const
1226 {
1227 return Accessor<rank_, dim, constness, P - 1, Number>(
1228 tensor, merge(previous_indices, i, rank_ - P));
1229 }
1230
1231
1232
1233 template <int rank_, int dim, bool constness, typename Number>
1235 Accessor<rank_, dim, constness, 1, Number>::Accessor(
1236 tensor_type &tensor,
1237 const TableIndices<rank_> &previous_indices)
1238 : tensor(tensor)
1239 , previous_indices(previous_indices)
1240 {}
1241
1242
1243
1244 template <int rank_, int dim, bool constness, typename Number>
1245 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1246 typename Accessor<rank_, dim, constness, 1, Number>::reference
1247 Accessor<rank_, dim, constness, 1, Number>::operator[](
1248 const unsigned int i)
1249 {
1250 return tensor(merge(previous_indices, i, rank_ - 1));
1251 }
1252
1253
1254 template <int rank_, int dim, bool constness, typename Number>
1256 typename Accessor<rank_, dim, constness, 1, Number>::reference
1257 Accessor<rank_, dim, constness, 1, Number>::operator[](
1258 const unsigned int i) const
1259 {
1260 return tensor(merge(previous_indices, i, rank_ - 1));
1261 }
1262 } // namespace SymmetricTensorAccessors
1263} // namespace internal
1264
1265
1266
1267template <int rank_, int dim, typename Number>
1268template <typename OtherNumber>
1272{
1273 static_assert(rank == 2, "This function is only implemented for rank==2");
1274 for (unsigned int d = 0; d < dim; ++d)
1275 for (unsigned int e = 0; e < d; ++e)
1276 Assert(t[d][e] == t[e][d],
1277 ExcMessage("The incoming Tensor must be exactly symmetric."));
1278
1279 for (unsigned int d = 0; d < dim; ++d)
1280 data[d] = t[d][d];
1281
1282 for (unsigned int d = 0, c = 0; d < dim; ++d)
1283 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1284 data[dim + c] = t[d][e];
1285}
1286
1287
1288
1289template <int rank_, int dim, typename Number>
1290template <typename OtherNumber>
1294 : data(initializer.data)
1295{}
1296
1297
1298
1299template <int rank_, int dim, typename Number>
1300DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1302 const Number (&array)[n_independent_components])
1303 : data(
1304 *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1305{
1306 // ensure that the reinterpret_cast above actually works
1307 Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1309}
1310
1311
1312
1313template <int rank_, int dim, typename Number>
1314template <typename OtherNumber>
1315DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1319{
1320 data = t.data;
1321 return *this;
1322}
1323
1324
1325
1326template <int rank_, int dim, typename Number>
1327DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1330{
1332 ExcMessage("Only assignment with zero is allowed"));
1333 (void)d;
1334
1336
1337 return *this;
1338}
1339
1340
1341namespace internal
1342{
1343 namespace SymmetricTensorImplementation
1344 {
1345 template <int dim, typename Number>
1346 constexpr inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1348 {
1350
1351 // diagonal entries are stored first
1352 for (unsigned int d = 0; d < dim; ++d)
1353 t[d][d] = s.access_raw_entry(d);
1354
1355 // off-diagonal entries come next, row by row
1356 for (unsigned int d = 0, c = 0; d < dim; ++d)
1357 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1358 {
1359 t[d][e] = s.access_raw_entry(dim + c);
1360 t[e][d] = s.access_raw_entry(dim + c);
1361 }
1362 return t;
1363 }
1364
1365
1366 template <int dim, typename Number>
1369 {
1370 // utilize the symmetry properties of SymmetricTensor<4,dim>
1371 // discussed in the class documentation to avoid accessing all
1372 // independent elements of the input tensor more than once
1374
1375 for (unsigned int i = 0; i < dim; ++i)
1376 for (unsigned int j = i; j < dim; ++j)
1377 for (unsigned int k = 0; k < dim; ++k)
1378 for (unsigned int l = k; l < dim; ++l)
1379 t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1380 t[TableIndices<4>(j, i, k, l)] =
1381 t[TableIndices<4>(j, i, l, k)] =
1382 st[TableIndices<4>(i, j, k, l)];
1383
1384 return t;
1385 }
1386
1387
1388 template <typename Number>
1389 struct Inverse<2, 1, Number>
1390 {
1391 constexpr static inline DEAL_II_ALWAYS_INLINE
1394 {
1396
1397 tmp[0][0] = 1.0 / t[0][0];
1398
1399 return tmp;
1400 }
1401 };
1402
1403
1404 template <typename Number>
1405 struct Inverse<2, 2, Number>
1406 {
1407 constexpr static inline DEAL_II_ALWAYS_INLINE
1410 {
1412
1413 // Sympy result: ([
1414 // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1415 // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1416 const TableIndices<2> idx_00(0, 0);
1417 const TableIndices<2> idx_01(0, 1);
1418 const TableIndices<2> idx_11(1, 1);
1419 const Number inv_det_t =
1420 1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1421 tmp[idx_00] = t[idx_11];
1422 tmp[idx_01] = -t[idx_01];
1423 tmp[idx_11] = t[idx_00];
1424 tmp *= inv_det_t;
1425
1426 return tmp;
1427 }
1428 };
1429
1430
1431 template <typename Number>
1432 struct Inverse<2, 3, Number>
1433 {
1436 {
1438
1439 // Sympy result: ([
1440 // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1441 // 2*t01*t02*t12 - t02**2*t11),
1442 // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1443 // 2*t01*t02*t12 - t02**2*t11),
1444 // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1445 // 2*t01*t02*t12 - t02**2*t11)],
1446 // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1447 // 2*t01*t02*t12 - t02**2*t11),
1448 // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1449 // 2*t01*t02*t12 - t02**2*t11),
1450 // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1451 // 2*t01*t02*t12 + t02**2*t11)],
1452 // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1453 // 2*t01*t02*t12 - t02**2*t11),
1454 // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1455 // 2*t01*t02*t12 + t02**2*t11),
1456 // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1457 // 2*t01*t02*t12 + t02**2*t11)] ])
1458 //
1459 // =
1460 //
1461 // [ (t11*t22 - t12**2)/det_t,
1462 // (-t01*t22 + t02*t12)/det_t,
1463 // (t01*t12 - t02*t11)/det_t],
1464 // [ (-t01*t22 + t02*t12)/det_t,
1465 // (t00*t22 - t02**2)/det_t,
1466 // (-t00*t12 + t01*t02)/det_t],
1467 // [ (t01*t12 - t02*t11)/det_t,
1468 // (-t00*t12 + t01*t02)/det_t,
1469 // (t00*t11 - t01**2)/det_t] ])
1470 //
1471 // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1472 // 2*t01*t02*t12 - t02**2*t11)
1473 const TableIndices<2> idx_00(0, 0);
1474 const TableIndices<2> idx_01(0, 1);
1475 const TableIndices<2> idx_02(0, 2);
1476 const TableIndices<2> idx_11(1, 1);
1477 const TableIndices<2> idx_12(1, 2);
1478 const TableIndices<2> idx_22(2, 2);
1479 const Number inv_det_t =
1480 1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1481 t[idx_00] * t[idx_12] * t[idx_12] -
1482 t[idx_01] * t[idx_01] * t[idx_22] +
1483 2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1484 t[idx_02] * t[idx_02] * t[idx_11]);
1485 tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1486 tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1487 tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1488 tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1489 tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1490 tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1491 tmp *= inv_det_t;
1492
1493 return tmp;
1494 }
1495 };
1496
1497
1498 template <typename Number>
1499 struct Inverse<4, 1, Number>
1500 {
1503 {
1505 tmp.data[0][0] = 1.0 / t.data[0][0];
1506 return tmp;
1507 }
1508 };
1509
1510
1511 template <typename Number>
1512 struct Inverse<4, 2, Number>
1513 {
1516 {
1518
1519 // Inverting this tensor is a little more complicated than necessary,
1520 // since we store the data of 't' as a 3x3 matrix t.data, but the
1521 // product between a rank-4 and a rank-2 tensor is really not the
1522 // product between this matrix and the 3-vector of a rhs, but rather
1523 //
1524 // B.vec = t.data * mult * A.vec
1525 //
1526 // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1527 // capture the fact that we need to add up both the c_ij12*a_12 and the
1528 // c_ij21*a_21 terms.
1529 //
1530 // In addition, in this scheme, the identity tensor has the matrix
1531 // representation mult^-1.
1532 //
1533 // The inverse of 't' therefore has the matrix representation
1534 //
1535 // inv.data = mult^-1 * t.data^-1 * mult^-1
1536 //
1537 // in order to compute it, let's first compute the inverse of t.data and
1538 // put it into tmp.data; at the end of the function we then scale the
1539 // last row and column of the inverse by 1/2, corresponding to the left
1540 // and right multiplication with mult^-1.
1541 const Number t4 = t.data[0][0] * t.data[1][1],
1542 t6 = t.data[0][0] * t.data[1][2],
1543 t8 = t.data[0][1] * t.data[1][0],
1544 t00 = t.data[0][2] * t.data[1][0],
1545 t01 = t.data[0][1] * t.data[2][0],
1546 t04 = t.data[0][2] * t.data[2][0],
1547 t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1548 t8 * t.data[2][2] + t00 * t.data[2][1] +
1549 t01 * t.data[1][2] - t04 * t.data[1][1]);
1550 tmp.data[0][0] =
1551 (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1552 tmp.data[0][1] =
1553 -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1554 tmp.data[0][2] =
1555 -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1556 tmp.data[1][0] =
1557 -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1558 tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1559 tmp.data[1][2] = -(t6 - t00) * t07;
1560 tmp.data[2][0] =
1561 -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1562 tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1563 tmp.data[2][2] = (t4 - t8) * t07;
1564
1565 // scale last row and column as mentioned
1566 // above
1567 tmp.data[2][0] /= 2;
1568 tmp.data[2][1] /= 2;
1569 tmp.data[0][2] /= 2;
1570 tmp.data[1][2] /= 2;
1571 tmp.data[2][2] /= 4;
1572
1573 return tmp;
1574 }
1575 };
1576
1577
1578 template <typename Number>
1579 struct Inverse<4, 3, Number>
1580 {
1583 {
1585
1586 // This function follows the exact same scheme as the 2d case, except
1587 // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1588 // Instead, we use the Gauss-Jordan algorithm implemented for
1589 // FullMatrix. For historical reasons the following code is copied from
1590 // there, with the tangential benefit that we do not need to copy the
1591 // tensor entries to and from the FullMatrix.
1592 const unsigned int N = 6;
1593
1594 // First get an estimate of the size of the elements of this matrix,
1595 // for later checks whether the pivot element is large enough, or
1596 // whether we have to fear that the matrix is not regular.
1598 for (unsigned int i = 0; i < N; ++i)
1600 const Number typical_diagonal_element =
1601 diagonal_sum / static_cast<double>(N);
1603
1604 unsigned int p[N];
1605 for (unsigned int i = 0; i < N; ++i)
1606 p[i] = i;
1607
1608 for (unsigned int j = 0; j < N; ++j)
1609 {
1610 // Pivot search: search that part of the line on and right of the
1611 // diagonal for the largest element.
1612 Number max = numbers::NumberTraits<Number>::abs(tmp.data[j][j]);
1613 unsigned int r = j;
1614 for (unsigned int i = j + 1; i < N; ++i)
1615 if (numbers::NumberTraits<Number>::abs(tmp.data[i][j]) > max)
1616 {
1618 r = i;
1619 }
1620
1621 // Check whether the pivot is too small
1622 Assert(max > 1.e-16 * typical_diagonal_element,
1623 ExcMessage("This tensor seems to be noninvertible"));
1624
1625 // Row interchange
1626 if (r > j)
1627 {
1628 for (unsigned int k = 0; k < N; ++k)
1629 std::swap(tmp.data[j][k], tmp.data[r][k]);
1630
1631 std::swap(p[j], p[r]);
1632 }
1633
1634 // Transformation
1635 const Number hr = 1. / tmp.data[j][j];
1636 tmp.data[j][j] = hr;
1637 for (unsigned int k = 0; k < N; ++k)
1638 {
1639 if (k == j)
1640 continue;
1641 for (unsigned int i = 0; i < N; ++i)
1642 {
1643 if (i == j)
1644 continue;
1645 tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1646 }
1647 }
1648 for (unsigned int i = 0; i < N; ++i)
1649 {
1650 tmp.data[i][j] *= hr;
1651 tmp.data[j][i] *= -hr;
1652 }
1653 tmp.data[j][j] = hr;
1654 }
1655
1656 // Column interchange
1657 Number hv[N];
1658 for (unsigned int i = 0; i < N; ++i)
1659 {
1660 for (unsigned int k = 0; k < N; ++k)
1661 hv[p[k]] = tmp.data[i][k];
1662 for (unsigned int k = 0; k < N; ++k)
1663 tmp.data[i][k] = hv[k];
1664 }
1665
1666 // Scale rows and columns. The mult matrix
1667 // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1668 for (unsigned int i = 3; i < 6; ++i)
1669 for (unsigned int j = 0; j < 3; ++j)
1670 tmp.data[i][j] /= 2;
1671
1672 for (unsigned int i = 0; i < 3; ++i)
1673 for (unsigned int j = 3; j < 6; ++j)
1674 tmp.data[i][j] /= 2;
1675
1676 for (unsigned int i = 3; i < 6; ++i)
1677 for (unsigned int j = 3; j < 6; ++j)
1678 tmp.data[i][j] /= 4;
1679
1680 return tmp;
1681 }
1682 };
1683
1684 } // namespace SymmetricTensorImplementation
1685} // namespace internal
1686
1687
1688
1689template <int rank_, int dim, typename Number>
1692 const
1693{
1694 return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1695}
1696
1697
1698
1699template <int rank_, int dim, typename Number>
1700DEAL_II_HOST constexpr bool
1703{
1704 return data == t.data;
1705}
1706
1707
1708
1709template <int rank_, int dim, typename Number>
1710DEAL_II_HOST constexpr bool
1713{
1714 return data != t.data;
1715}
1716
1717
1718
1719template <int rank_, int dim, typename Number>
1720template <typename OtherNumber>
1721DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1725{
1726 data += t.data;
1727 return *this;
1728}
1729
1730
1731
1732template <int rank_, int dim, typename Number>
1733template <typename OtherNumber>
1734DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1738{
1739 data -= t.data;
1740 return *this;
1741}
1742
1743
1744
1745template <int rank_, int dim, typename Number>
1746template <typename OtherNumber>
1747DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1750{
1751 data *= d;
1752 return *this;
1753}
1754
1755
1756
1757template <int rank_, int dim, typename Number>
1758template <typename OtherNumber>
1759DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1762{
1763 data /= d;
1764 return *this;
1765}
1766
1767
1768
1769template <int rank_, int dim, typename Number>
1770DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1773{
1774 SymmetricTensor tmp = *this;
1775 tmp.data = -tmp.data;
1776 return tmp;
1777}
1778
1779
1780
1781template <int rank_, int dim, typename Number>
1782DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
1784{
1785 data.clear();
1786}
1787
1788
1789
1790template <int rank_, int dim, typename Number>
1791DEAL_II_HOST constexpr std::size_t
1793{
1794 // all memory consists of statically allocated memory of the current
1795 // object, no pointers
1797}
1798
1799
1800
1801namespace internal
1802{
1806 template <int dim, typename Number, typename OtherNumber = Number>
1808 typename SymmetricTensorAccessors::
1809 double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1811 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1812 base_tensor_type &data,
1813 const typename SymmetricTensorAccessors::
1814 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1815 {
1816 using result_type = typename SymmetricTensorAccessors::
1817 double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1818
1819 switch (dim)
1820 {
1821 case 1:
1822 return data[0] * sdata[0];
1823
1824 default:
1825 // Start with the non-diagonal part. These values appear
1826 // twice in the matrix, but are only stored once. So we can
1827 // get the double-contraction sum for these elements using
1828 // only one multiplication each, and at the end multiplying
1829 // things by 2.
1830 result_type sum = data[dim] * sdata[dim];
1831 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1832 sum += data[d] * sdata[d];
1833 sum += sum; // sum *= 2
1834
1835 // Now add the contributions from the diagonal
1836 for (unsigned int d = 0; d < dim; ++d)
1837 sum += data[d] * sdata[d];
1838 return sum;
1839 }
1840 }
1841
1842
1843
1848 template <int dim, typename Number, typename OtherNumber = Number>
1850 typename SymmetricTensorAccessors::
1851 double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1853 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1854 base_tensor_type &data,
1855 const typename SymmetricTensorAccessors::
1856 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1857 {
1858 using result_type = typename SymmetricTensorAccessors::
1859 double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1860 using value_type = typename SymmetricTensorAccessors::
1861 double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1862
1863 const unsigned int data_dim = SymmetricTensorAccessors::
1864 StorageType<2, dim, value_type>::n_independent_components;
1865 value_type tmp[data_dim]{};
1866 for (unsigned int i = 0; i < data_dim; ++i)
1867 tmp[i] =
1869 return result_type(tmp);
1870 }
1871
1872
1873
1878 template <int dim, typename Number, typename OtherNumber = Number>
1880 typename SymmetricTensorAccessors::StorageType<
1881 2,
1882 dim,
1883 typename SymmetricTensorAccessors::
1884 double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1887 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1888 base_tensor_type &data,
1889 const typename SymmetricTensorAccessors::
1890 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1891 {
1892 using value_type = typename SymmetricTensorAccessors::
1893 double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1894 using base_tensor_type = typename SymmetricTensorAccessors::
1895 StorageType<2, dim, value_type>::base_tensor_type;
1896
1897 base_tensor_type tmp;
1898 for (unsigned int i = 0; i < tmp.dimension; ++i)
1899 {
1900 // Start with the non-diagonal part. These values appear
1901 // twice in the matrix, but are only stored once. So we can
1902 // get the double-contraction sum for these elements using
1903 // only one multiplication each, and at the end multiplying
1904 // things by 2.
1905 value_type sum = data[dim] * sdata[dim][i];
1906 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1907 sum += data[d] * sdata[d][i];
1908 sum += sum; // sum *= 2
1909
1910 // Now add the contributions from the diagonal
1911 for (unsigned int d = 0; d < dim; ++d)
1912 sum += data[d] * sdata[d][i];
1913 tmp[i] = sum;
1914 }
1915 return tmp;
1916 }
1917
1918
1919
1923 template <int dim, typename Number, typename OtherNumber = Number>
1925 typename SymmetricTensorAccessors::StorageType<
1926 4,
1927 dim,
1928 typename SymmetricTensorAccessors::
1929 double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1932 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1933 base_tensor_type &data,
1934 const typename SymmetricTensorAccessors::
1935 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1936 {
1937 using value_type = typename SymmetricTensorAccessors::
1938 double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1939 using base_tensor_type = typename SymmetricTensorAccessors::
1940 StorageType<4, dim, value_type>::base_tensor_type;
1941
1942 const unsigned int data_dim = SymmetricTensorAccessors::
1943 StorageType<2, dim, value_type>::n_independent_components;
1944 base_tensor_type tmp;
1945 for (unsigned int i = 0; i < data_dim; ++i)
1946 for (unsigned int j = 0; j < data_dim; ++j)
1947 {
1948 // Start with the non-diagonal part
1949 for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1950 tmp[i][j] += data[i][d] * sdata[d][j];
1951 tmp[i][j] += tmp[i][j]; // tmp[i][j] *= 2;
1952
1953 // Now add the contributions from the diagonal
1954 for (unsigned int d = 0; d < dim; ++d)
1955 tmp[i][j] += data[i][d] * sdata[d][j];
1956 }
1957 return tmp;
1958 }
1959
1960} // end of namespace internal
1961
1962
1963
1964template <int rank_, int dim, typename Number>
1965template <typename OtherNumber>
1967 typename internal::SymmetricTensorAccessors::
1968 double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1971{
1972 // Dispatch to functions that know the types of the involved
1973 // arguments via overloads.
1974 return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1975 s.data);
1976}
1977
1978
1979
1980template <int rank_, int dim, typename Number>
1981template <typename OtherNumber>
1983 typename internal::SymmetricTensorAccessors::
1984 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1987{
1988 typename internal::SymmetricTensorAccessors::
1989 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1990 tmp.data =
1991 internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1992 s.data);
1993 return tmp;
1994}
1995
1996
1997
1998// internal namespace to switch between the
1999// access of different tensors. There used to
2000// be explicit instantiations before for
2001// different ranks and dimensions, but since
2002// we now allow for templates on the data
2003// type, and since we cannot partially
2004// specialize the implementation, this got
2005// into a separate namespace
2006namespace internal
2007{
2008 namespace SymmetricTensorImplementation
2009 {
2010 // a function to do the unrolling from a set of indices to a
2011 // scalar index into the array in which we store the elements of
2012 // a symmetric tensor
2013 //
2014 // this function is for rank-2 tensors
2015 template <int dim>
2016 constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
2018 {
2019 AssertIndexRange(indices[0], dim);
2020 AssertIndexRange(indices[1], dim);
2021
2022 switch (dim)
2023 {
2024 case 1:
2025 {
2026 return 0;
2027 }
2028 case 2:
2029 {
2031 {{{0, 2}}, {{2, 1}}}};
2032 return table[indices[0]][indices[1]];
2033 }
2034 case 3:
2035 {
2037 {{{0, 3, 4}}, {{3, 1, 5}}, {{4, 5, 2}}}};
2038 return table[indices[0]][indices[1]];
2039 }
2040 case 4:
2041 {
2043 {{{0, 4, 5, 6}},
2044 {{4, 1, 7, 8}},
2045 {{5, 7, 2, 9}},
2046 {{6, 8, 9, 3}}}};
2047 return table[indices[0]][indices[1]];
2048 }
2049 default:
2050 // for the remainder, manually figure out the numbering
2051 {
2052 if (indices[0] == indices[1])
2053 return indices[0];
2054
2056 std::min(indices[0], indices[1]),
2057 std::max(indices[0], indices[1]));
2058
2059 // Here (d, e) are the row and column of the symmetric matrix and
2060 // 'dim + c' is the index into the Tensor<1, dim> actually used
2061 // for storage.
2062 unsigned int c = 0;
2063 for (unsigned int d = 0; d < dim; ++d)
2064 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2065 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2066 return dim + c;
2067
2068 // should never get here:
2070 return 0;
2071 }
2072 }
2073 }
2074
2075 // a function to do the unrolling from a set of indices to a
2076 // scalar index into the array in which we store the elements of
2077 // a symmetric tensor
2078 //
2079 // this function is for tensors of ranks not already handled
2080 // above
2081 template <int dim, int rank_>
2082 constexpr inline unsigned int
2084 {
2085 (void)indices;
2088 }
2089 } // namespace SymmetricTensorImplementation
2090
2091 template <int dim, typename Number>
2092 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2094 typename SymmetricTensorAccessors::
2095 StorageType<2, dim, Number>::base_tensor_type &data)
2096 {
2098 indices)];
2099 }
2100
2101
2102
2103 template <int dim, typename Number>
2104 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2106 const typename SymmetricTensorAccessors::
2107 StorageType<2, dim, Number>::base_tensor_type &data)
2108 {
2110 indices)];
2111 }
2112
2113
2114
2115 template <int dim, typename Number>
2116 constexpr inline Number &
2118 typename SymmetricTensorAccessors::
2119 StorageType<4, dim, Number>::base_tensor_type &data)
2120 {
2121 switch (dim)
2122 {
2123 case 1:
2124 return data[0][0];
2125
2126 case 2:
2127 // each entry of the tensor can be thought of as an entry in a
2128 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2129 // rank-2 tensors. this is the format in which we store rank-4
2130 // tensors. determine which position the present entry is
2131 // stored in
2132 {
2133 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2134 return data[base_index[indices[0]][indices[1]]]
2135 [base_index[indices[2]][indices[3]]];
2136 }
2137 case 3:
2138 // each entry of the tensor can be thought of as an entry in a
2139 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2140 // rank-2 tensors. this is the format in which we store rank-4
2141 // tensors. determine which position the present entry is
2142 // stored in
2143 {
2144 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2145 {3, 1, 5},
2146 {4, 5, 2}};
2147 return data[base_index[indices[0]][indices[1]]]
2148 [base_index[indices[2]][indices[3]]];
2149 }
2150
2151 default:
2153 }
2154
2155 // The code should never reach here.
2156 // We cannot return a static variable, as this class must support number
2157 // types that require no instances of the number type to be in scope during
2158 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2159 return data[0][0];
2160 }
2161
2162
2163 template <int dim, typename Number>
2164 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2166 const typename SymmetricTensorAccessors::
2167 StorageType<4, dim, Number>::base_tensor_type &data)
2168 {
2169 switch (dim)
2170 {
2171 case 1:
2172 return data[0][0];
2173
2174 case 2:
2175 // each entry of the tensor can be thought of as an entry in a
2176 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2177 // rank-2 tensors. this is the format in which we store rank-4
2178 // tensors. determine which position the present entry is
2179 // stored in
2180 {
2181 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2182 return data[base_index[indices[0]][indices[1]]]
2183 [base_index[indices[2]][indices[3]]];
2184 }
2185 case 3:
2186 // each entry of the tensor can be thought of as an entry in a
2187 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2188 // rank-2 tensors. this is the format in which we store rank-4
2189 // tensors. determine which position the present entry is
2190 // stored in
2191 {
2192 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2193 {3, 1, 5},
2194 {4, 5, 2}};
2195 return data[base_index[indices[0]][indices[1]]]
2196 [base_index[indices[2]][indices[3]]];
2197 }
2198
2199 default:
2201 }
2202
2203 // The code should never reach here.
2204 // We cannot return a static variable, as this class must support number
2205 // types that require no instances of the number type to be in scope during
2206 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2207 return data[0][0];
2208 }
2209
2210} // end of namespace internal
2211
2212
2213
2214template <int rank_, int dim, typename Number>
2215DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2217 const TableIndices<rank_> &indices)
2218{
2219 for (unsigned int r = 0; r < rank; ++r)
2220 AssertIndexRange(indices[r], dimension);
2221 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2222}
2223
2224
2225
2226template <int rank_, int dim, typename Number>
2227DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2229 const TableIndices<rank_> &indices) const
2230{
2231 for (unsigned int r = 0; r < rank; ++r)
2232 AssertIndexRange(indices[r], dimension);
2233 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2234}
2235
2236
2237
2238namespace internal
2239{
2240 namespace SymmetricTensorImplementation
2241 {
2242 template <int rank_>
2243 constexpr TableIndices<rank_>
2244 get_partially_filled_indices(const unsigned int row,
2245 const std::integral_constant<int, 2> &)
2246 {
2248 }
2249
2250
2251 template <int rank_>
2252 constexpr TableIndices<rank_>
2253 get_partially_filled_indices(const unsigned int row,
2254 const std::integral_constant<int, 4> &)
2255 {
2256 return TableIndices<rank_>(row,
2260 }
2261 } // namespace SymmetricTensorImplementation
2262} // namespace internal
2263
2264
2265template <int rank_, int dim, typename Number>
2266DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE internal::
2267 SymmetricTensorAccessors::Accessor<rank_, dim, true, rank_ - 1, Number>
2268 SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2269{
2270 return internal::SymmetricTensorAccessors::
2271 Accessor<rank_, dim, true, rank_ - 1, Number>(
2272 *this,
2273 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2274 rank_>(row, std::integral_constant<int, rank_>()));
2275}
2276
2277
2278
2279template <int rank_, int dim, typename Number>
2280DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE internal::
2281 SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2283{
2284 return internal::SymmetricTensorAccessors::
2285 Accessor<rank_, dim, false, rank_ - 1, Number>(
2286 *this,
2287 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2288 rank_>(row, std::integral_constant<int, rank_>()));
2289}
2290
2291
2292
2293template <int rank_, int dim, typename Number>
2294DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE const Number &
2296 const TableIndices<rank_> &indices) const
2297{
2298 return operator()(indices);
2299}
2300
2301
2302
2303template <int rank_, int dim, typename Number>
2304DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2306 const TableIndices<rank_> &indices)
2307{
2308 return operator()(indices);
2309}
2310
2311
2312
2313template <int rank_, int dim, typename Number>
2314inline Number *
2316{
2317 return std::addressof(this->access_raw_entry(0));
2318}
2319
2320
2321
2322template <int rank_, int dim, typename Number>
2323inline const Number *
2325{
2326 return std::addressof(this->access_raw_entry(0));
2327}
2328
2329
2330
2331template <int rank_, int dim, typename Number>
2332inline Number *
2334{
2336}
2337
2338
2339
2340template <int rank_, int dim, typename Number>
2341inline const Number *
2343{
2345}
2346
2347
2348
2349template <int rank_, int dim, typename Number>
2350DEAL_II_HOST constexpr inline const Number &
2352 const unsigned int index) const
2353{
2355 if constexpr (rank == 2)
2356 return data[index];
2357 else
2358 return data[decltype(data)::unrolled_to_component_indices(index)];
2359}
2360
2361
2362
2363template <int rank_, int dim, typename Number>
2364DEAL_II_HOST constexpr inline Number &
2366{
2368 if constexpr (rank == 2)
2369 return data[index];
2370 else
2371 return data[decltype(data)::unrolled_to_component_indices(index)];
2372}
2373
2374
2375
2376namespace internal
2377{
2378 template <int dim, typename Number>
2379 constexpr inline typename numbers::NumberTraits<Number>::real_type
2380 compute_norm(const typename SymmetricTensorAccessors::
2381 StorageType<2, dim, Number>::base_tensor_type &data)
2382 {
2383 // Make things work with AD types
2384 using std::sqrt;
2385 switch (dim)
2386 {
2387 case 1:
2388 return numbers::NumberTraits<Number>::abs(data[0]);
2389
2390 case 2:
2394
2395 case 3:
2402
2403 default:
2404 {
2405 typename numbers::NumberTraits<Number>::real_type return_value =
2407
2408 for (unsigned int d = 0; d < dim; ++d)
2409 return_value +=
2411 for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2412 return_value +=
2414
2415 return sqrt(return_value);
2416 }
2417 }
2418 }
2419
2420
2421
2422 template <int dim, typename Number>
2423 constexpr inline typename numbers::NumberTraits<Number>::real_type
2424 compute_norm(const typename SymmetricTensorAccessors::
2425 StorageType<4, dim, Number>::base_tensor_type &data)
2426 {
2427 // Make things work with AD types
2428 using std::sqrt;
2429 switch (dim)
2430 {
2431 case 1:
2432 return numbers::NumberTraits<Number>::abs(data[0][0]);
2433
2434 default:
2435 {
2436 typename numbers::NumberTraits<Number>::real_type return_value =
2438
2439 const unsigned int n_independent_components = data.dimension;
2440
2441 for (unsigned int i = 0; i < dim; ++i)
2442 for (unsigned int j = 0; j < dim; ++j)
2443 return_value +=
2445 for (unsigned int i = 0; i < dim; ++i)
2446 for (unsigned int j = dim; j < n_independent_components; ++j)
2447 return_value +=
2449 for (unsigned int i = dim; i < n_independent_components; ++i)
2450 for (unsigned int j = 0; j < dim; ++j)
2451 return_value +=
2453 for (unsigned int i = dim; i < n_independent_components; ++i)
2454 for (unsigned int j = dim; j < n_independent_components; ++j)
2455 return_value +=
2457
2458 return sqrt(return_value);
2459 }
2460 }
2461 }
2462
2463} // end of namespace internal
2464
2465
2466
2467template <int rank_, int dim, typename Number>
2470{
2471 return internal::compute_norm<dim, Number>(data);
2472}
2473
2474
2475
2476template <int rank_, int dim, typename Number>
2477DEAL_II_HOST constexpr unsigned int
2479 const TableIndices<rank_> &indices)
2480{
2481 return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2482 dim>(indices);
2483}
2484
2485
2486
2487namespace internal
2488{
2489 namespace SymmetricTensorImplementation
2490 {
2491 // a function to do the inverse of the unrolling from a set of
2492 // indices to a scalar index into the array in which we store
2493 // the elements of a symmetric tensor. in other words, it goes
2494 // from the scalar index into the array to a set of indices of
2495 // the tensor
2496 //
2497 // this function is for rank-2 tensors
2498 template <int dim>
2499 constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2500 unrolled_to_component_indices(const unsigned int i,
2501 const std::integral_constant<int, 2> &)
2502 {
2503 Assert(
2506 i,
2507 0,
2509 switch (dim)
2510 {
2511 case 1:
2512 {
2513 return {0, 0};
2514 }
2515
2516 case 2:
2517 {
2518 const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2519 TableIndices<2>(1, 1),
2520 TableIndices<2>(0, 1)};
2521 return table[i];
2522 }
2523
2524 case 3:
2525 {
2526 const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2527 TableIndices<2>(1, 1),
2528 TableIndices<2>(2, 2),
2529 TableIndices<2>(0, 1),
2530 TableIndices<2>(0, 2),
2531 TableIndices<2>(1, 2)};
2532 return table[i];
2533 }
2534
2535 default:
2536 if (i < dim)
2537 return {i, i};
2538
2539 for (unsigned int d = 0, c = dim; d < dim; ++d)
2540 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2541 if (c == i)
2542 return {d, e};
2543
2544 // should never get here:
2546 return {0, 0};
2547 }
2548 }
2549
2550 // a function to do the inverse of the unrolling from a set of
2551 // indices to a scalar index into the array in which we store
2552 // the elements of a symmetric tensor. in other words, it goes
2553 // from the scalar index into the array to a set of indices of
2554 // the tensor
2555 //
2556 // this function is for tensors of a rank not already handled
2557 // above
2558 template <int dim, int rank_>
2559 constexpr inline std::enable_if_t<rank_ != 2, TableIndices<rank_>>
2560 unrolled_to_component_indices(const unsigned int i,
2561 const std::integral_constant<int, rank_> &)
2562 {
2563 (void)i;
2564 Assert(
2565 (i <
2567 ExcIndexRange(i,
2568 0,
2570 n_independent_components));
2572 return TableIndices<rank_>();
2573 }
2574
2575 } // namespace SymmetricTensorImplementation
2576} // namespace internal
2577
2578template <int rank_, int dim, typename Number>
2581 const unsigned int i)
2582{
2583 return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2584 dim>(i, std::integral_constant<int, rank_>());
2585}
2586
2587
2588
2589template <int rank_, int dim, typename Number>
2590template <class Archive>
2591inline void
2593{
2594 ar &data;
2595}
2596
2597
2598#endif // DOXYGEN
2599
2600/* ----------------- Non-member functions operating on tensors. ------------ */
2601
2602
2615template <int rank_, int dim, typename Number, typename OtherNumber>
2616DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2626
2627
2640template <int rank_, int dim, typename Number, typename OtherNumber>
2641DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2651
2652
2660template <int rank_, int dim, typename Number, typename OtherNumber>
2665{
2666 return Tensor<rank_, dim, Number>(left) + right;
2667}
2668
2669
2677template <int rank_, int dim, typename Number, typename OtherNumber>
2682{
2683 return left + Tensor<rank_, dim, OtherNumber>(right);
2684}
2685
2686
2694template <int rank_, int dim, typename Number, typename OtherNumber>
2699{
2700 return Tensor<rank_, dim, Number>(left) - right;
2701}
2702
2703
2711template <int rank_, int dim, typename Number, typename OtherNumber>
2716{
2717 return left - Tensor<rank_, dim, OtherNumber>(right);
2718}
2719
2720
2721
2722template <int dim, typename Number>
2725{
2726 switch (dim)
2727 {
2728 case 1:
2729 return t.data[0];
2730 case 2:
2731 return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2732 case 3:
2733 {
2734 // in analogy to general tensors, but
2735 // there's something to be simplified for
2736 // the present case
2737 const Number tmp = t.data[3] * t.data[4] * t.data[5];
2738 return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2739 t.data[0] * t.data[5] * t.data[5] -
2740 t.data[1] * t.data[4] * t.data[4] -
2741 t.data[2] * t.data[3] * t.data[3]);
2742 }
2743 default:
2746 }
2747}
2748
2749
2750
2762template <int dim, typename Number>
2768
2769
2770
2771template <int dim, typename Number>
2772DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
2774{
2775 Number t = d.data[0];
2776 for (unsigned int i = 1; i < dim; ++i)
2777 t += d.data[i];
2778 return t;
2779}
2780
2781
2793template <int dim, typename Number>
2794DEAL_II_HOST constexpr Number
2796{
2797 return trace(t);
2798}
2799
2800
2812template <typename Number>
2813DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2818
2819
2820
2839template <typename Number>
2840DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2842{
2843 return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2844}
2845
2846
2847
2856template <typename Number>
2857DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2859{
2860 return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2861 t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2862}
2863
2864
2865
2873template <typename Number>
2874std::array<Number, 1>
2876
2877
2878
2901template <typename Number>
2902std::array<Number, 2>
2904
2905
2906
2929template <typename Number>
2930std::array<Number, 3>
2932
2933
2934
2935namespace internal
2936{
2937 namespace SymmetricTensorImplementation
2938 {
2950 template <int dim, typename Number>
2951 void
2954 std::array<Number, dim> &d,
2955 std::array<Number, dim - 1> &e);
2956
2957
2958
2972 template <int dim, typename Number>
2973 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2975
2976
2977
2991 template <int dim, typename Number>
2992 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2994
2995
2996
3010 template <typename Number>
3011 std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3013
3014
3015
3030 template <typename Number>
3031 std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3033
3038 template <int dim, typename Number>
3040 {
3041 using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3042 bool
3044 {
3045 return lhs.first > rhs.first;
3046 }
3047 };
3048
3049 } // namespace SymmetricTensorImplementation
3050
3051} // namespace internal
3052
3053
3054
3055// The line below is to ensure that doxygen puts the full description
3056// of this global enumeration into the documentation
3057// See https://stackoverflow.com/a/1717984
3069{
3079 hybrid,
3097 jacobi
3098};
3099
3100
3101
3112template <int dim, typename Number>
3113std::array<std::pair<Number, Tensor<1, dim, Number>>,
3114 std::integral_constant<int, dim>::value>
3118
3119
3120
3129template <int rank_, int dim, typename Number>
3132{
3133 return t;
3134}
3135
3136
3137
3138template <int dim, typename Number>
3139DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3142{
3144
3145 // subtract scaled trace from the diagonal
3146 const Number tr = trace(t) * internal::NumberType<Number>::value(1.0 / dim);
3147 for (unsigned int i = 0; i < dim; ++i)
3148 tmp.data[i] -= tr;
3149
3150 return tmp;
3151}
3152
3153
3154
3155template <int dim, typename Number>
3159{
3160 // create a default constructed matrix filled with
3161 // zeros, then set the diagonal elements to one
3163 switch (dim)
3164 {
3165 case 1:
3166 tmp.data[0] = internal::NumberType<Number>::value(1.);
3167 break;
3168 case 2:
3169 tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3170 break;
3171 case 3:
3172 tmp.data[0] = tmp.data[1] = tmp.data[2] =
3174 break;
3175 default:
3176 for (unsigned int d = 0; d < dim; ++d)
3177 tmp.data[d] = internal::NumberType<Number>::value(1.);
3178 }
3179 return tmp;
3180}
3181
3182
3183
3184template <int dim, typename Number>
3187{
3189
3190 // fill the elements treating the diagonal
3191 for (unsigned int i = 0; i < dim; ++i)
3192 for (unsigned int j = 0; j < dim; ++j)
3193 tmp.data[i][j] =
3194 internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3195
3196 // then fill the ones that copy over the
3197 // non-diagonal elements. note that during
3198 // the double-contraction, we handle the
3199 // off-diagonal elements twice, so simply
3200 // copying requires a weight of 1/2
3201 for (unsigned int i = dim;
3204 ++i)
3205 tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3206
3207 return tmp;
3208}
3209
3210
3211
3212template <int dim, typename Number>
3216{
3218
3219 // fill the elements treating the diagonal
3220 for (unsigned int i = 0; i < dim; ++i)
3221 tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3222
3223 // then fill the ones that copy over the
3224 // non-diagonal elements. note that during
3225 // the double-contraction, we handle the
3226 // off-diagonal elements twice, so simply
3227 // copying requires a weight of 1/2
3228 for (unsigned int i = dim;
3231 ++i)
3232 tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3233
3234 return tmp;
3235}
3236
3237
3238
3248template <int dim, typename Number>
3255
3256
3257
3268template <int dim, typename Number>
3275
3276
3277
3299template <int dim, typename Number>
3303{
3305
3306 // fill only the elements really needed
3307 for (unsigned int i = 0; i < dim; ++i)
3308 for (unsigned int j = i; j < dim; ++j)
3309 for (unsigned int k = 0; k < dim; ++k)
3310 for (unsigned int l = k; l < dim; ++l)
3311 tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3312
3313 return tmp;
3314}
3315
3341template <int dim, typename Number>
3342std::pair<SymmetricTensor<2, dim, Number>, SymmetricTensor<2, dim, Number>>
3344{
3345 Assert(dim <= 3, ExcNotImplemented());
3346
3347 const std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3349
3350 std::pair<SymmetricTensor<2, dim, Number>, SymmetricTensor<2, dim, Number>>
3352
3355
3357 for (unsigned int i = 0; i < dim; ++i)
3358 if (eigen_system[i].first > 0)
3361 eigen_system[i].second));
3362
3364 for (unsigned int i = 0; i < dim; ++i)
3365 if (eigen_system[i].first < 0)
3368 eigen_system[i].second));
3369
3371}
3372
3405template <int dim, typename Number>
3406std::tuple<SymmetricTensor<2, dim, Number>,
3412{
3413 Assert(dim <= 3, ExcNotImplemented());
3414
3415 auto heaviside_function{[](const double x) {
3416 if (std::fabs(x) < 1.0e-16)
3417 return 0.5;
3418 if (x > 0)
3419 return 1.0;
3420 else
3421 return 0.0;
3422 }};
3423
3424 std::tuple<SymmetricTensor<2, dim, Number>,
3429
3430 auto &[positive_part_tensor,
3434
3435 const std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3437
3439 for (unsigned int i = 0; i < dim; ++i)
3440 if (eigen_system[i].first > 0)
3443 eigen_system[i].second));
3444
3446 for (unsigned int i = 0; i < dim; ++i)
3447 if (eigen_system[i].first < 0)
3450 eigen_system[i].second));
3451
3452 std::array<SymmetricTensor<2, dim, Number>, dim> M;
3453 for (unsigned int a = 0; a < dim; ++a)
3454 M[a] =
3456
3457 std::array<SymmetricTensor<4, dim, Number>, dim> Q;
3458 for (unsigned int a = 0; a < dim; ++a)
3459 Q[a] = outer_product(M[a], M[a]);
3460
3461 std::array<std::array<SymmetricTensor<4, dim, Number>, dim>, dim> G;
3462 for (unsigned int a = 0; a < dim; ++a)
3463 for (unsigned int b = 0; b < dim; ++b)
3464 for (unsigned int i = 0; i < dim; ++i)
3465 for (unsigned int j = 0; j < dim; ++j)
3466 for (unsigned int k = 0; k < dim; ++k)
3467 for (unsigned int l = 0; l < dim; ++l)
3468 G[a][b][i][j][k][l] =
3469 M[a][i][k] * M[b][j][l] + M[a][i][l] * M[b][j][k];
3470
3471 // positive P
3473 for (unsigned int a = 0; a < dim; ++a)
3474 {
3475 double lambda_a = eigen_system[a].first;
3477 for (unsigned int b = 0; b < dim; ++b)
3478 {
3479 if (b != a)
3480 {
3481 double lambda_b = eigen_system[b].first;
3482
3483 double v_ab = 0.0;
3484 if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
3485 v_ab = (std::fmax(lambda_a, 0.0) - std::fmax(lambda_b, 0.0)) /
3486 (lambda_a - lambda_b);
3487 else
3488 v_ab = 0.5 * (heaviside_function(lambda_a) +
3490
3491 positive_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
3492 }
3493 }
3494 }
3495
3496 // negative P
3498 for (unsigned int a = 0; a < dim; ++a)
3499 {
3500 double lambda_a = eigen_system[a].first;
3502 for (unsigned int b = 0; b < dim; ++b)
3503 {
3504 if (b != a)
3505 {
3506 double lambda_b = eigen_system[b].first;
3507
3508 double v_ab = 0.0;
3509 if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
3510 v_ab = (std::fmin(lambda_a, 0.0) - std::fmin(lambda_b, 0.0)) /
3511 (lambda_a - lambda_b);
3512 else
3513 v_ab = 0.5 * (heaviside_function(-lambda_a) +
3515
3516 negative_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
3517 }
3518 }
3519 }
3520
3522}
3523
3531template <int dim, typename Number>
3532DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3535{
3537 for (unsigned int d = 0; d < dim; ++d)
3538 result[d][d] = t[d][d];
3539
3540 const Number half = internal::NumberType<Number>::value(0.5);
3541 for (unsigned int d = 0; d < dim; ++d)
3542 for (unsigned int e = d + 1; e < dim; ++e)
3543 result[d][e] = (t[d][e] + t[e][d]) * half;
3544 return result;
3545}
3546
3547
3548
3556template <int rank_, int dim, typename Number>
3557DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3559 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3560{
3562 tt *= factor;
3563 return tt;
3564}
3565
3566
3567
3575template <int rank_, int dim, typename Number>
3577operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3578{
3579 // simply forward to the other operator
3580 return t * factor;
3581}
3582
3583
3584
3609template <int rank_, int dim, typename Number, typename OtherNumber>
3611 rank_,
3612 dim,
3613 typename ProductType<Number,
3614 typename EnableIfScalar<OtherNumber>::type>::type>
3616 const OtherNumber &factor)
3617{
3618 // form the product. we have to convert the two factors into the final
3619 // type via explicit casts because, for awkward reasons, the C++
3620 // standard committee saw it fit to not define an
3621 // operator*(float,std::complex<double>)
3622 // (as well as with switched arguments and double<->float).
3626 return tt;
3627}
3628
3629
3630
3638template <int rank_, int dim, typename Number, typename OtherNumber>
3640 rank_,
3641 dim,
3642 typename ProductType<OtherNumber,
3643 typename EnableIfScalar<Number>::type>::type>
3644operator*(const Number &factor,
3646{
3647 // simply forward to the other operator with switched arguments
3648 return (t * factor);
3649}
3650
3651
3652
3658template <int rank_, int dim, typename Number, typename OtherNumber>
3659DEAL_II_HOST constexpr inline SymmetricTensor<
3660 rank_,
3661 dim,
3662 typename ProductType<Number,
3663 typename EnableIfScalar<OtherNumber>::type>::type>
3672
3673
3674
3681template <int rank_, int dim>
3683operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3684{
3686 tt *= factor;
3687 return tt;
3688}
3689
3690
3691
3698template <int rank_, int dim>
3700operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3701{
3703 tt *= factor;
3704 return tt;
3705}
3706
3707
3708
3714template <int rank_, int dim>
3716operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3717{
3719 tt /= factor;
3720 return tt;
3721}
3722
3732template <int dim, typename Number, typename OtherNumber>
3740
3741
3754template <int dim, typename Number, typename OtherNumber>
3755DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3759{
3761 typename ProductType<Number, OtherNumber>::type>::value(0.0);
3762 for (unsigned int i = 0; i < dim; ++i)
3763 for (unsigned int j = 0; j < dim; ++j)
3764 s += t1[i][j] * t2[i][j];
3765 return s;
3766}
3767
3768
3781template <int dim, typename Number, typename OtherNumber>
3786{
3787 return scalar_product(t2, t1);
3788}
3789
3790
3805template <typename Number, typename OtherNumber>
3806DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
3811{
3812 tmp[0][0] = t[0][0][0][0] * s[0][0];
3813}
3814
3815
3816
3831template <typename Number, typename OtherNumber>
3832DEAL_II_HOST constexpr inline void
3837{
3838 tmp[0][0] = t[0][0][0][0] * s[0][0];
3839}
3840
3841
3842
3857template <typename Number, typename OtherNumber>
3858DEAL_II_HOST constexpr inline void
3863{
3864 const unsigned int dim = 2;
3865
3866 for (unsigned int i = 0; i < dim; ++i)
3867 for (unsigned int j = i; j < dim; ++j)
3868 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3869 2 * t[i][j][0][1] * s[0][1];
3870}
3871
3872
3873
3888template <typename Number, typename OtherNumber>
3889DEAL_II_HOST constexpr inline void
3894{
3895 const unsigned int dim = 2;
3896
3897 for (unsigned int i = 0; i < dim; ++i)
3898 for (unsigned int j = i; j < dim; ++j)
3899 tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3900 2 * s[0][1] * t[0][1][i][j];
3901}
3902
3903
3904
3919template <typename Number, typename OtherNumber>
3920DEAL_II_HOST constexpr inline void
3925{
3926 const unsigned int dim = 3;
3927
3928 for (unsigned int i = 0; i < dim; ++i)
3929 for (unsigned int j = i; j < dim; ++j)
3930 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3931 t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3932 2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3933}
3934
3935
3936
3951template <typename Number, typename OtherNumber>
3952DEAL_II_HOST constexpr inline void
3957{
3958 const unsigned int dim = 3;
3959
3960 for (unsigned int i = 0; i < dim; ++i)
3961 for (unsigned int j = i; j < dim; ++j)
3962 tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3963 s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3964 2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3965}
3966
3967
3968
3975template <int dim, typename Number, typename OtherNumber>
3976DEAL_II_HOST constexpr Tensor<1,
3977 dim,
3981{
3983 for (unsigned int i = 0; i < dim; ++i)
3984 {
3985 dest[i] = src1[i][0] * src2[0];
3986 for (unsigned int j = 1; j < dim; ++j)
3987 dest[i] += src1[i][j] * src2[j];
3988 }
3989 return dest;
3990}
3991
3992
3999template <int dim, typename Number, typename OtherNumber>
4000DEAL_II_HOST constexpr Tensor<1,
4001 dim,
4005{
4006 // this is easy for symmetric tensors:
4007 return src2 * src1;
4008}
4009
4010
4011
4031template <int rank_1,
4032 int rank_2,
4033 int dim,
4034 typename Number,
4035 typename OtherNumber>
4037 typename Tensor<rank_1 + rank_2 - 2,
4038 dim,
4039 typename ProductType<Number, OtherNumber>::type>::tensor_type
4045
4046
4047
4067template <int rank_1,
4068 int rank_2,
4069 int dim,
4070 typename Number,
4071 typename OtherNumber>
4073 typename Tensor<rank_1 + rank_2 - 2,
4074 dim,
4075 typename ProductType<Number, OtherNumber>::type>::tensor_type
4081
4082
4083
4093template <int dim, typename Number>
4094inline std::ostream &
4095operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4096{
4097 // make our lives a bit simpler by outputting
4098 // the tensor through the operator for the
4099 // general Tensor class
4101
4102 for (unsigned int i = 0; i < dim; ++i)
4103 for (unsigned int j = 0; j < dim; ++j)
4104 tt[i][j] = t[i][j];
4105
4106 return out << tt;
4107}
4108
4109
4110
4120template <int dim, typename Number>
4121inline std::ostream &
4122operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4123{
4124 // make our lives a bit simpler by outputting
4125 // the tensor through the operator for the
4126 // general Tensor class
4128
4129 for (unsigned int i = 0; i < dim; ++i)
4130 for (unsigned int j = 0; j < dim; ++j)
4131 for (unsigned int k = 0; k < dim; ++k)
4132 for (unsigned int l = 0; l < dim; ++l)
4133 tt[i][j][k][l] = t[i][j][k][l];
4134
4135 return out << tt;
4136}
4137
4138
4140
4141#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)
const Number * begin_raw() const
const Number * end_raw() const
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
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)
Number * begin_raw()
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)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
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< 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)
Number * end_raw()
DEAL_II_HOST constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
Number * begin_raw()
static constexpr unsigned int rank
Definition tensor.h:490
constexpr void clear()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
friend class Tensor
Definition tensor.h:882
static constexpr unsigned int dimension
Definition tensor.h:485
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
static constexpr unsigned int n_independent_components
Definition tensor.h:506
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:109
#define DEAL_II_DEPRECATED
Definition config.h:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_CONSTEXPR
Definition config.h:236
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
std::ostream & operator<<(std::ostream &out, const DerivativeForm< order, dim, spacedim, Number > &df)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
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)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
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 bool value_is_zero(const Number &value)
Definition numbers.h:936
static const unsigned int invalid_unsigned_int
Definition types.h:220
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 > &)
#define DEAL_II_HOST
Definition numbers.h:45
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
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:593
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 > &)
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()
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
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()