Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
symmetric_tensor.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2005 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_symmetric_tensor_h
17#define dealii_symmetric_tensor_h
18
19
20#include <deal.II/base/config.h>
21
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)
262 return {new_index, numbers::invalid_unsigned_int};
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:
306 Assert(false, ExcInternalError());
307 return {};
308 }
309 }
310
311
318 template <int rank1,
319 int rank2,
320 int dim,
321 typename Number,
322 typename OtherNumber = Number>
324 {
326 using type =
327 ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
328 };
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 =
745 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 &
1057 access_raw_entry(const unsigned int unrolled_index);
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
1194namespace internal
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>
1347 convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
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>
1367 constexpr ::Tensor<4, dim, Number>
1368 convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
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
1392 ::SymmetricTensor<2, 1, Number>
1393 value(const ::SymmetricTensor<2, 1, Number> &t)
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
1408 ::SymmetricTensor<2, 2, Number>
1409 value(const ::SymmetricTensor<2, 2, Number> &t)
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 {
1434 constexpr static ::SymmetricTensor<2, 3, Number>
1435 value(const ::SymmetricTensor<2, 3, Number> &t)
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 {
1501 constexpr static inline ::SymmetricTensor<4, 1, Number>
1502 value(const ::SymmetricTensor<4, 1, Number> &t)
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 {
1514 constexpr static inline ::SymmetricTensor<4, 2, Number>
1515 value(const ::SymmetricTensor<4, 2, Number> &t)
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 {
1581 static ::SymmetricTensor<4, 3, Number>
1582 value(const ::SymmetricTensor<4, 3, Number> &t)
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.
1597 Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1598 for (unsigned int i = 0; i < N; ++i)
1599 diagonal_sum += numbers::NumberTraits<Number>::abs(tmp.data[i][i]);
1600 const Number typical_diagonal_element =
1601 diagonal_sum / static_cast<double>(N);
1602 (void)typical_diagonal_element;
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{
1803 template <int dim, typename Number, typename OtherNumber = Number>
1805 typename SymmetricTensorAccessors::
1806 double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1807 perform_double_contraction(
1808 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1809 base_tensor_type &data,
1810 const typename SymmetricTensorAccessors::
1811 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1812 {
1813 using result_type = typename SymmetricTensorAccessors::
1814 double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1815
1816 switch (dim)
1817 {
1818 case 1:
1819 return data[0] * sdata[0];
1820 default:
1821 // Start with the non-diagonal part to avoid some multiplications by
1822 // 2.
1823
1824 result_type sum = data[dim] * sdata[dim];
1825 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1826 sum += data[d] * sdata[d];
1827 sum += sum; // sum = sum * 2.;
1828
1829 // Now add the contributions from the diagonal
1830 for (unsigned int d = 0; d < dim; ++d)
1831 sum += data[d] * sdata[d];
1832 return sum;
1833 }
1834 }
1835
1836
1837
1838 template <int dim, typename Number, typename OtherNumber = Number>
1840 typename SymmetricTensorAccessors::
1841 double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1842 perform_double_contraction(
1843 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1844 base_tensor_type &data,
1845 const typename SymmetricTensorAccessors::
1846 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1847 {
1848 using result_type = typename SymmetricTensorAccessors::
1849 double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1850 using value_type = typename SymmetricTensorAccessors::
1851 double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1852
1853 const unsigned int data_dim = SymmetricTensorAccessors::
1854 StorageType<2, dim, value_type>::n_independent_components;
1855 value_type tmp[data_dim]{};
1856 for (unsigned int i = 0; i < data_dim; ++i)
1857 tmp[i] =
1858 perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1859 return result_type(tmp);
1860 }
1861
1862
1863
1864 template <int dim, typename Number, typename OtherNumber = Number>
1866 typename SymmetricTensorAccessors::StorageType<
1867 2,
1868 dim,
1869 typename SymmetricTensorAccessors::
1870 double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1871 base_tensor_type
1872 perform_double_contraction(
1873 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1874 base_tensor_type &data,
1875 const typename SymmetricTensorAccessors::
1876 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1877 {
1878 using value_type = typename SymmetricTensorAccessors::
1879 double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1880 using base_tensor_type = typename SymmetricTensorAccessors::
1881 StorageType<2, dim, value_type>::base_tensor_type;
1882
1883 base_tensor_type tmp;
1884 for (unsigned int i = 0; i < tmp.dimension; ++i)
1885 {
1886 // Start with the non-diagonal part
1887 value_type sum = data[dim] * sdata[dim][i];
1888 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1889 sum += data[d] * sdata[d][i];
1890 sum += sum; // sum = sum * 2.;
1891
1892 // Now add the contributions from the diagonal
1893 for (unsigned int d = 0; d < dim; ++d)
1894 sum += data[d] * sdata[d][i];
1895 tmp[i] = sum;
1896 }
1897 return tmp;
1898 }
1899
1900
1901
1902 template <int dim, typename Number, typename OtherNumber = Number>
1904 typename SymmetricTensorAccessors::StorageType<
1905 4,
1906 dim,
1907 typename SymmetricTensorAccessors::
1908 double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1909 base_tensor_type
1910 perform_double_contraction(
1911 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1912 base_tensor_type &data,
1913 const typename SymmetricTensorAccessors::
1914 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1915 {
1916 using value_type = typename SymmetricTensorAccessors::
1917 double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1918 using base_tensor_type = typename SymmetricTensorAccessors::
1919 StorageType<4, dim, value_type>::base_tensor_type;
1920
1921 const unsigned int data_dim = SymmetricTensorAccessors::
1922 StorageType<2, dim, value_type>::n_independent_components;
1923 base_tensor_type tmp;
1924 for (unsigned int i = 0; i < data_dim; ++i)
1925 for (unsigned int j = 0; j < data_dim; ++j)
1926 {
1927 // Start with the non-diagonal part
1928 for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1929 tmp[i][j] += data[i][d] * sdata[d][j];
1930 tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1931
1932 // Now add the contributions from the diagonal
1933 for (unsigned int d = 0; d < dim; ++d)
1934 tmp[i][j] += data[i][d] * sdata[d][j];
1935 }
1936 return tmp;
1937 }
1938
1939} // end of namespace internal
1940
1941
1942
1943template <int rank_, int dim, typename Number>
1944template <typename OtherNumber>
1946 typename internal::SymmetricTensorAccessors::
1947 double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1950{
1951 // need to have two different function calls
1952 // because a scalar and rank-2 tensor are not
1953 // the same data type (see internal function
1954 // above)
1955 return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1956 s.data);
1957}
1958
1959
1960
1961template <int rank_, int dim, typename Number>
1962template <typename OtherNumber>
1964 typename internal::SymmetricTensorAccessors::
1965 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1968{
1969 typename internal::SymmetricTensorAccessors::
1970 double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1971 tmp.data =
1972 internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1973 s.data);
1974 return tmp;
1975}
1976
1977
1978
1979// internal namespace to switch between the
1980// access of different tensors. There used to
1981// be explicit instantiations before for
1982// different ranks and dimensions, but since
1983// we now allow for templates on the data
1984// type, and since we cannot partially
1985// specialize the implementation, this got
1986// into a separate namespace
1987namespace internal
1988{
1989 template <int dim, typename Number>
1990 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1991 symmetric_tensor_access(const TableIndices<2> &indices,
1992 typename SymmetricTensorAccessors::
1993 StorageType<2, dim, Number>::base_tensor_type &data)
1994 {
1995 // 1d is very simple and done first
1996 if (dim == 1)
1997 return data[0];
1998
1999 // first treat the main diagonal elements, which are stored consecutively
2000 // at the beginning
2001 if (indices[0] == indices[1])
2002 return data[indices[0]];
2003
2004 // the rest is messier and requires a few switches.
2005 switch (dim)
2006 {
2007 case 2:
2008 // at least for the 2x2 case it is reasonably simple
2009 Assert(((indices[0] == 1) && (indices[1] == 0)) ||
2010 ((indices[0] == 0) && (indices[1] == 1)),
2012 return data[2];
2013
2014 default:
2015 // to do the rest, sort our indices before comparing
2016 {
2017 TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2018 std::max(indices[0], indices[1]));
2019 for (unsigned int d = 0, c = 0; d < dim; ++d)
2020 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2021 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2022 return data[dim + c];
2023 Assert(false, ExcInternalError());
2024 }
2025 }
2026
2027 // The code should never reach here.
2028 // We cannot return a static variable, as this class must support number
2029 // types that require no instances of the number type to be in scope during
2030 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2031 return data[0];
2032 }
2033
2034
2035
2036 template <int dim, typename Number>
2037 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2038 symmetric_tensor_access(const TableIndices<2> &indices,
2039 const typename SymmetricTensorAccessors::
2040 StorageType<2, dim, Number>::base_tensor_type &data)
2041 {
2042 // 1d is very simple and done first
2043 if (dim == 1)
2044 return data[0];
2045
2046 // first treat the main diagonal elements, which are stored consecutively
2047 // at the beginning
2048 if (indices[0] == indices[1])
2049 return data[indices[0]];
2050
2051 // the rest is messier and requires a few switches.
2052 switch (dim)
2053 {
2054 case 2:
2055 // at least for the 2x2 case it is reasonably simple
2056 Assert(((indices[0] == 1) && (indices[1] == 0)) ||
2057 ((indices[0] == 0) && (indices[1] == 1)),
2059 return data[2];
2060
2061 default:
2062 // to do the rest, sort our indices before comparing
2063 {
2064 TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2065 std::max(indices[0], indices[1]));
2066 for (unsigned int d = 0, c = 0; d < dim; ++d)
2067 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2068 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2069 return data[dim + c];
2070 Assert(false, ExcInternalError());
2071 }
2072 }
2073
2074 // The code should never reach here.
2075 // We cannot return a static variable, as this class must support number
2076 // types that require no instances of the number type to be in scope during
2077 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2078 return data[0];
2079 }
2080
2081
2082
2083 template <int dim, typename Number>
2084 constexpr inline Number &
2085 symmetric_tensor_access(const TableIndices<4> &indices,
2086 typename SymmetricTensorAccessors::
2087 StorageType<4, dim, Number>::base_tensor_type &data)
2088 {
2089 switch (dim)
2090 {
2091 case 1:
2092 return data[0][0];
2093
2094 case 2:
2095 // each entry of the tensor can be thought of as an entry in a
2096 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2097 // rank-2 tensors. this is the format in which we store rank-4
2098 // tensors. determine which position the present entry is
2099 // stored in
2100 {
2101 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2102 return data[base_index[indices[0]][indices[1]]]
2103 [base_index[indices[2]][indices[3]]];
2104 }
2105 case 3:
2106 // each entry of the tensor can be thought of as an entry in a
2107 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2108 // rank-2 tensors. this is the format in which we store rank-4
2109 // tensors. determine which position the present entry is
2110 // stored in
2111 {
2112 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2113 {3, 1, 5},
2114 {4, 5, 2}};
2115 return data[base_index[indices[0]][indices[1]]]
2116 [base_index[indices[2]][indices[3]]];
2117 }
2118
2119 default:
2120 Assert(false, ExcNotImplemented());
2121 }
2122
2123 // The code should never reach here.
2124 // We cannot return a static variable, as this class must support number
2125 // types that require no instances of the number type to be in scope during
2126 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2127 return data[0][0];
2128 }
2129
2130
2131 template <int dim, typename Number>
2132 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2133 symmetric_tensor_access(const TableIndices<4> &indices,
2134 const typename SymmetricTensorAccessors::
2135 StorageType<4, dim, Number>::base_tensor_type &data)
2136 {
2137 switch (dim)
2138 {
2139 case 1:
2140 return data[0][0];
2141
2142 case 2:
2143 // each entry of the tensor can be thought of as an entry in a
2144 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2145 // rank-2 tensors. this is the format in which we store rank-4
2146 // tensors. determine which position the present entry is
2147 // stored in
2148 {
2149 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2150 return data[base_index[indices[0]][indices[1]]]
2151 [base_index[indices[2]][indices[3]]];
2152 }
2153 case 3:
2154 // each entry of the tensor can be thought of as an entry in a
2155 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2156 // rank-2 tensors. this is the format in which we store rank-4
2157 // tensors. determine which position the present entry is
2158 // stored in
2159 {
2160 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2161 {3, 1, 5},
2162 {4, 5, 2}};
2163 return data[base_index[indices[0]][indices[1]]]
2164 [base_index[indices[2]][indices[3]]];
2165 }
2166
2167 default:
2168 Assert(false, ExcNotImplemented());
2169 }
2170
2171 // The code should never reach here.
2172 // We cannot return a static variable, as this class must support number
2173 // types that require no instances of the number type to be in scope during
2174 // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2175 return data[0][0];
2176 }
2177
2178} // end of namespace internal
2179
2180
2181
2182template <int rank_, int dim, typename Number>
2183DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2185 const TableIndices<rank_> &indices)
2186{
2187 for (unsigned int r = 0; r < rank; ++r)
2188 AssertIndexRange(indices[r], dimension);
2189 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2190}
2191
2192
2193
2194template <int rank_, int dim, typename Number>
2195DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2197 const TableIndices<rank_> &indices) const
2198{
2199 for (unsigned int r = 0; r < rank; ++r)
2200 AssertIndexRange(indices[r], dimension);
2201 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2202}
2203
2204
2205
2206namespace internal
2207{
2208 namespace SymmetricTensorImplementation
2209 {
2210 template <int rank_>
2211 constexpr TableIndices<rank_>
2212 get_partially_filled_indices(const unsigned int row,
2213 const std::integral_constant<int, 2> &)
2214 {
2216 }
2217
2218
2219 template <int rank_>
2220 constexpr TableIndices<rank_>
2221 get_partially_filled_indices(const unsigned int row,
2222 const std::integral_constant<int, 4> &)
2223 {
2224 return TableIndices<rank_>(row,
2228 }
2229 } // namespace SymmetricTensorImplementation
2230} // namespace internal
2231
2232
2233template <int rank_, int dim, typename Number>
2234DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE internal::
2235 SymmetricTensorAccessors::Accessor<rank_, dim, true, rank_ - 1, Number>
2236 SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2237{
2238 return internal::SymmetricTensorAccessors::
2239 Accessor<rank_, dim, true, rank_ - 1, Number>(
2240 *this,
2241 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2242 rank_>(row, std::integral_constant<int, rank_>()));
2243}
2244
2245
2246
2247template <int rank_, int dim, typename Number>
2248DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE internal::
2249 SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2251{
2252 return internal::SymmetricTensorAccessors::
2253 Accessor<rank_, dim, false, rank_ - 1, Number>(
2254 *this,
2255 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2256 rank_>(row, std::integral_constant<int, rank_>()));
2257}
2258
2259
2260
2261template <int rank_, int dim, typename Number>
2262DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE const Number &
2264 const TableIndices<rank_> &indices) const
2265{
2266 return operator()(indices);
2267}
2268
2269
2270
2271template <int rank_, int dim, typename Number>
2272DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2274 const TableIndices<rank_> &indices)
2275{
2276 return operator()(indices);
2277}
2278
2279
2280
2281template <int rank_, int dim, typename Number>
2282inline Number *
2284{
2285 return std::addressof(this->access_raw_entry(0));
2286}
2287
2288
2289
2290template <int rank_, int dim, typename Number>
2291inline const Number *
2293{
2294 return std::addressof(this->access_raw_entry(0));
2295}
2296
2297
2298
2299template <int rank_, int dim, typename Number>
2300inline Number *
2302{
2303 return begin_raw() + n_independent_components;
2304}
2305
2306
2307
2308template <int rank_, int dim, typename Number>
2309inline const Number *
2311{
2312 return begin_raw() + n_independent_components;
2313}
2314
2315
2316
2317namespace internal
2318{
2319 namespace SymmetricTensorImplementation
2320 {
2321 template <int dim, typename Number>
2322 constexpr unsigned int
2323 entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2324 const unsigned int index)
2325 {
2326 return index;
2327 }
2328
2329
2330 template <int dim, typename Number>
2331 constexpr ::TableIndices<2>
2332 entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2333 const unsigned int index)
2334 {
2337 }
2338
2339 } // namespace SymmetricTensorImplementation
2340} // namespace internal
2341
2342
2343
2344template <int rank_, int dim, typename Number>
2345DEAL_II_HOST constexpr inline const Number &
2347 const unsigned int index) const
2348{
2349 AssertIndexRange(index, n_independent_components);
2350 return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2351 index)];
2352}
2353
2354
2355
2356template <int rank_, int dim, typename Number>
2357DEAL_II_HOST constexpr inline Number &
2359{
2360 AssertIndexRange(index, n_independent_components);
2361 return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2362 index)];
2363}
2364
2365
2366
2367namespace internal
2368{
2369 template <int dim, typename Number>
2370 constexpr inline typename numbers::NumberTraits<Number>::real_type
2371 compute_norm(const typename SymmetricTensorAccessors::
2372 StorageType<2, dim, Number>::base_tensor_type &data)
2373 {
2374 // Make things work with AD types
2375 using std::sqrt;
2376 switch (dim)
2377 {
2378 case 1:
2379 return numbers::NumberTraits<Number>::abs(data[0]);
2380
2381 case 2:
2385
2386 case 3:
2393
2394 default:
2395 {
2396 typename numbers::NumberTraits<Number>::real_type return_value =
2398
2399 for (unsigned int d = 0; d < dim; ++d)
2400 return_value +=
2402 for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2403 return_value +=
2405
2406 return sqrt(return_value);
2407 }
2408 }
2409 }
2410
2411
2412
2413 template <int dim, typename Number>
2414 constexpr inline typename numbers::NumberTraits<Number>::real_type
2415 compute_norm(const typename SymmetricTensorAccessors::
2416 StorageType<4, dim, Number>::base_tensor_type &data)
2417 {
2418 // Make things work with AD types
2419 using std::sqrt;
2420 switch (dim)
2421 {
2422 case 1:
2423 return numbers::NumberTraits<Number>::abs(data[0][0]);
2424
2425 default:
2426 {
2427 typename numbers::NumberTraits<Number>::real_type return_value =
2429
2430 const unsigned int n_independent_components = data.dimension;
2431
2432 for (unsigned int i = 0; i < dim; ++i)
2433 for (unsigned int j = 0; j < dim; ++j)
2434 return_value +=
2436 for (unsigned int i = 0; i < dim; ++i)
2437 for (unsigned int j = dim; j < n_independent_components; ++j)
2438 return_value +=
2440 for (unsigned int i = dim; i < n_independent_components; ++i)
2441 for (unsigned int j = 0; j < dim; ++j)
2442 return_value +=
2444 for (unsigned int i = dim; i < n_independent_components; ++i)
2445 for (unsigned int j = dim; j < n_independent_components; ++j)
2446 return_value +=
2448
2449 return sqrt(return_value);
2450 }
2451 }
2452 }
2453
2454} // end of namespace internal
2455
2456
2457
2458template <int rank_, int dim, typename Number>
2461{
2462 return internal::compute_norm<dim, Number>(data);
2463}
2464
2465
2466
2467namespace internal
2468{
2469 namespace SymmetricTensorImplementation
2470 {
2471 // a function to do the unrolling from a set of indices to a
2472 // scalar index into the array in which we store the elements of
2473 // a symmetric tensor
2474 //
2475 // this function is for rank-2 tensors
2476 template <int dim>
2477 constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
2479 {
2480 AssertIndexRange(indices[0], dim);
2481 AssertIndexRange(indices[1], dim);
2482
2483 switch (dim)
2484 {
2485 case 1:
2486 {
2487 return 0;
2488 }
2489
2490 case 2:
2491 {
2492 constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2493 return table[indices[0]][indices[1]];
2494 }
2495
2496 case 3:
2497 {
2498 constexpr unsigned int table[3][3] = {{0, 3, 4},
2499 {3, 1, 5},
2500 {4, 5, 2}};
2501 return table[indices[0]][indices[1]];
2502 }
2503
2504 case 4:
2505 {
2506 constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2507 {4, 1, 7, 8},
2508 {5, 7, 2, 9},
2509 {6, 8, 9, 3}};
2510 return table[indices[0]][indices[1]];
2511 }
2512
2513 default:
2514 // for the remainder, manually figure out the numbering
2515 {
2516 if (indices[0] == indices[1])
2517 return indices[0];
2518
2519 TableIndices<2> sorted_indices(indices);
2520 sorted_indices.sort();
2521
2522 for (unsigned int d = 0, c = 0; d < dim; ++d)
2523 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2524 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2525 return dim + c;
2526
2527 // should never get here:
2528 Assert(false, ExcInternalError());
2529 return 0;
2530 }
2531 }
2532 }
2533
2534 // a function to do the unrolling from a set of indices to a
2535 // scalar index into the array in which we store the elements of
2536 // a symmetric tensor
2537 //
2538 // this function is for tensors of ranks not already handled
2539 // above
2540 template <int dim, int rank_>
2541 constexpr inline unsigned int
2543 {
2544 (void)indices;
2545 Assert(false, ExcNotImplemented());
2547 }
2548 } // namespace SymmetricTensorImplementation
2549} // namespace internal
2550
2551
2552template <int rank_, int dim, typename Number>
2553DEAL_II_HOST constexpr unsigned int
2555 const TableIndices<rank_> &indices)
2556{
2557 return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2558 dim>(indices);
2559}
2560
2561
2562
2563namespace internal
2564{
2565 namespace SymmetricTensorImplementation
2566 {
2567 // a function to do the inverse of the unrolling from a set of
2568 // indices to a scalar index into the array in which we store
2569 // the elements of a symmetric tensor. in other words, it goes
2570 // from the scalar index into the array to a set of indices of
2571 // the tensor
2572 //
2573 // this function is for rank-2 tensors
2574 template <int dim>
2575 constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2576 unrolled_to_component_indices(const unsigned int i,
2577 const std::integral_constant<int, 2> &)
2578 {
2579 Assert(
2582 i,
2583 0,
2585 switch (dim)
2586 {
2587 case 1:
2588 {
2589 return {0, 0};
2590 }
2591
2592 case 2:
2593 {
2594 const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2595 TableIndices<2>(1, 1),
2596 TableIndices<2>(0, 1)};
2597 return table[i];
2598 }
2599
2600 case 3:
2601 {
2602 const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2603 TableIndices<2>(1, 1),
2604 TableIndices<2>(2, 2),
2605 TableIndices<2>(0, 1),
2606 TableIndices<2>(0, 2),
2607 TableIndices<2>(1, 2)};
2608 return table[i];
2609 }
2610
2611 default:
2612 if (i < dim)
2613 return {i, i};
2614
2615 for (unsigned int d = 0, c = dim; d < dim; ++d)
2616 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2617 if (c == i)
2618 return {d, e};
2619
2620 // should never get here:
2621 Assert(false, ExcInternalError());
2622 return {0, 0};
2623 }
2624 }
2625
2626 // a function to do the inverse of the unrolling from a set of
2627 // indices to a scalar index into the array in which we store
2628 // the elements of a symmetric tensor. in other words, it goes
2629 // from the scalar index into the array to a set of indices of
2630 // the tensor
2631 //
2632 // this function is for tensors of a rank not already handled
2633 // above
2634 template <int dim, int rank_>
2635 constexpr inline std::enable_if_t<rank_ != 2, TableIndices<rank_>>
2636 unrolled_to_component_indices(const unsigned int i,
2637 const std::integral_constant<int, rank_> &)
2638 {
2639 (void)i;
2640 Assert(
2641 (i <
2643 ExcIndexRange(i,
2644 0,
2646 n_independent_components));
2647 Assert(false, ExcNotImplemented());
2648 return TableIndices<rank_>();
2649 }
2650
2651 } // namespace SymmetricTensorImplementation
2652} // namespace internal
2653
2654template <int rank_, int dim, typename Number>
2657 const unsigned int i)
2658{
2659 return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2660 dim>(i, std::integral_constant<int, rank_>());
2661}
2662
2663
2664
2665template <int rank_, int dim, typename Number>
2666template <class Archive>
2667inline void
2668SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2669{
2670 ar &data;
2671}
2672
2673
2674#endif // DOXYGEN
2675
2676/* ----------------- Non-member functions operating on tensors. ------------ */
2677
2678
2691template <int rank_, int dim, typename Number, typename OtherNumber>
2692DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2696{
2698 tmp = left;
2699 tmp += right;
2700 return tmp;
2701}
2702
2703
2716template <int rank_, int dim, typename Number, typename OtherNumber>
2717DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2721{
2723 tmp = left;
2724 tmp -= right;
2725 return tmp;
2726}
2727
2728
2736template <int rank_, int dim, typename Number, typename OtherNumber>
2740 const Tensor<rank_, dim, OtherNumber> & right)
2741{
2742 return Tensor<rank_, dim, Number>(left) + right;
2743}
2744
2745
2753template <int rank_, int dim, typename Number, typename OtherNumber>
2758{
2759 return left + Tensor<rank_, dim, OtherNumber>(right);
2760}
2761
2762
2770template <int rank_, int dim, typename Number, typename OtherNumber>
2774 const Tensor<rank_, dim, OtherNumber> & right)
2775{
2776 return Tensor<rank_, dim, Number>(left) - right;
2777}
2778
2779
2787template <int rank_, int dim, typename Number, typename OtherNumber>
2792{
2793 return left - Tensor<rank_, dim, OtherNumber>(right);
2794}
2795
2796
2797
2798template <int dim, typename Number>
2801{
2802 switch (dim)
2803 {
2804 case 1:
2805 return t.data[0];
2806 case 2:
2807 return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2808 case 3:
2809 {
2810 // in analogy to general tensors, but
2811 // there's something to be simplified for
2812 // the present case
2813 const Number tmp = t.data[3] * t.data[4] * t.data[5];
2814 return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2815 t.data[0] * t.data[5] * t.data[5] -
2816 t.data[1] * t.data[4] * t.data[4] -
2817 t.data[2] * t.data[3] * t.data[3]);
2818 }
2819 default:
2820 Assert(false, ExcNotImplemented());
2822 }
2823}
2824
2825
2826
2838template <int dim, typename Number>
2841{
2842 return determinant(t);
2843}
2844
2845
2846
2847template <int dim, typename Number>
2848DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
2850{
2851 Number t = d.data[0];
2852 for (unsigned int i = 1; i < dim; ++i)
2853 t += d.data[i];
2854 return t;
2855}
2856
2857
2869template <int dim, typename Number>
2870DEAL_II_HOST constexpr Number
2872{
2873 return trace(t);
2874}
2875
2876
2888template <typename Number>
2889DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2891{
2893}
2894
2895
2896
2915template <typename Number>
2916DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2918{
2919 return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2920}
2921
2922
2923
2932template <typename Number>
2933DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2935{
2936 return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2937 t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2938}
2939
2940
2941
2949template <typename Number>
2950std::array<Number, 1>
2952
2953
2954
2977template <typename Number>
2978std::array<Number, 2>
2980
2981
2982
3005template <typename Number>
3006std::array<Number, 3>
3008
3009
3010
3011namespace internal
3012{
3013 namespace SymmetricTensorImplementation
3014 {
3052 template <int dim, typename Number>
3053 void
3054 tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
3056 std::array<Number, dim> & d,
3057 std::array<Number, dim - 1> & e);
3058
3059
3060
3100 template <int dim, typename Number>
3101 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3102 ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
3103
3104
3105
3145 template <int dim, typename Number>
3146 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3148
3149
3150
3164 template <typename Number>
3165 std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3166 hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3167
3168
3169
3202 template <typename Number>
3203 std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3204 hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3205
3210 template <int dim, typename Number>
3212 {
3213 using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3214 bool
3215 operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3216 {
3217 return lhs.first > rhs.first;
3218 }
3219 };
3220
3221 } // namespace SymmetricTensorImplementation
3222
3223} // namespace internal
3224
3225
3226
3227// The line below is to ensure that doxygen puts the full description
3228// of this global enumeration into the documentation
3229// See https://stackoverflow.com/a/1717984
3259{
3269 hybrid,
3287 jacobi
3288};
3289
3290
3291
3320template <int dim, typename Number>
3321std::array<std::pair<Number, Tensor<1, dim, Number>>,
3322 std::integral_constant<int, dim>::value>
3326
3327
3328
3337template <int rank_, int dim, typename Number>
3340{
3341 return t;
3342}
3343
3344
3345
3346template <int dim, typename Number>
3347DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3350{
3352
3353 // subtract scaled trace from the diagonal
3354 const Number tr = trace(t) * internal::NumberType<Number>::value(1.0 / dim);
3355 for (unsigned int i = 0; i < dim; ++i)
3356 tmp.data[i] -= tr;
3357
3358 return tmp;
3359}
3360
3361
3362
3363template <int dim, typename Number>
3367{
3368 // create a default constructed matrix filled with
3369 // zeros, then set the diagonal elements to one
3371 switch (dim)
3372 {
3373 case 1:
3375 break;
3376 case 2:
3377 tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3378 break;
3379 case 3:
3380 tmp.data[0] = tmp.data[1] = tmp.data[2] =
3382 break;
3383 default:
3384 for (unsigned int d = 0; d < dim; ++d)
3386 }
3387 return tmp;
3388}
3389
3390
3391
3392template <int dim, typename Number>
3395{
3397
3398 // fill the elements treating the diagonal
3399 for (unsigned int i = 0; i < dim; ++i)
3400 for (unsigned int j = 0; j < dim; ++j)
3401 tmp.data[i][j] =
3402 internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3403
3404 // then fill the ones that copy over the
3405 // non-diagonal elements. note that during
3406 // the double-contraction, we handle the
3407 // off-diagonal elements twice, so simply
3408 // copying requires a weight of 1/2
3409 for (unsigned int i = dim;
3410 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3411 n_rank2_components;
3412 ++i)
3414
3415 return tmp;
3416}
3417
3418
3419
3420template <int dim, typename Number>
3424{
3426
3427 // fill the elements treating the diagonal
3428 for (unsigned int i = 0; i < dim; ++i)
3430
3431 // then fill the ones that copy over the
3432 // non-diagonal elements. note that during
3433 // the double-contraction, we handle the
3434 // off-diagonal elements twice, so simply
3435 // copying requires a weight of 1/2
3436 for (unsigned int i = dim;
3437 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3438 n_rank2_components;
3439 ++i)
3441
3442 return tmp;
3443}
3444
3445
3446
3456template <int dim, typename Number>
3459{
3461 value(t);
3462}
3463
3464
3465
3476template <int dim, typename Number>
3479{
3481 value(t);
3482}
3483
3484
3485
3507template <int dim, typename Number>
3511{
3513
3514 // fill only the elements really needed
3515 for (unsigned int i = 0; i < dim; ++i)
3516 for (unsigned int j = i; j < dim; ++j)
3517 for (unsigned int k = 0; k < dim; ++k)
3518 for (unsigned int l = k; l < dim; ++l)
3519 tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3520
3521 return tmp;
3522}
3523
3524
3525
3533template <int dim, typename Number>
3534DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3537{
3539 for (unsigned int d = 0; d < dim; ++d)
3540 result[d][d] = t[d][d];
3541
3542 const Number half = internal::NumberType<Number>::value(0.5);
3543 for (unsigned int d = 0; d < dim; ++d)
3544 for (unsigned int e = d + 1; e < dim; ++e)
3545 result[d][e] = (t[d][e] + t[e][d]) * half;
3546 return result;
3547}
3548
3549
3550
3558template <int rank_, int dim, typename Number>
3559DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3561 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3562{
3564 tt *= factor;
3565 return tt;
3566}
3567
3568
3569
3577template <int rank_, int dim, typename Number>
3579operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3580{
3581 // simply forward to the other operator
3582 return t * factor;
3583}
3584
3585
3586
3611template <int rank_, int dim, typename Number, typename OtherNumber>
3613 rank_,
3614 dim,
3615 typename ProductType<Number,
3616 typename EnableIfScalar<OtherNumber>::type>::type>
3618 const OtherNumber & factor)
3619{
3620 // form the product. we have to convert the two factors into the final
3621 // type via explicit casts because, for awkward reasons, the C++
3622 // standard committee saw it fit to not define an
3623 // operator*(float,std::complex<double>)
3624 // (as well as with switched arguments and double<->float).
3625 using product_type = typename ProductType<Number, OtherNumber>::type;
3628 return tt;
3629}
3630
3631
3632
3640template <int rank_, int dim, typename Number, typename OtherNumber>
3642 rank_,
3643 dim,
3644 typename ProductType<OtherNumber,
3645 typename EnableIfScalar<Number>::type>::type>
3646operator*(const Number & factor,
3648{
3649 // simply forward to the other operator with switched arguments
3650 return (t * factor);
3651}
3652
3653
3654
3660template <int rank_, int dim, typename Number, typename OtherNumber>
3661DEAL_II_HOST constexpr inline SymmetricTensor<
3662 rank_,
3663 dim,
3664 typename ProductType<Number,
3665 typename EnableIfScalar<OtherNumber>::type>::type>
3667 const OtherNumber & factor)
3668{
3669 using product_type = typename ProductType<Number, OtherNumber>::type;
3672 return tt;
3673}
3674
3675
3676
3683template <int rank_, int dim>
3685operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3686{
3688 tt *= factor;
3689 return tt;
3690}
3691
3692
3693
3700template <int rank_, int dim>
3702operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3703{
3705 tt *= factor;
3706 return tt;
3707}
3708
3709
3710
3716template <int rank_, int dim>
3718operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3719{
3721 tt /= factor;
3722 return tt;
3723}
3724
3734template <int dim, typename Number, typename OtherNumber>
3739{
3740 return (t1 * t2);
3741}
3742
3743
3756template <int dim, typename Number, typename OtherNumber>
3757DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3760 const Tensor<2, dim, OtherNumber> & t2)
3761{
3763 typename ProductType<Number, OtherNumber>::type>::value(0.0);
3764 for (unsigned int i = 0; i < dim; ++i)
3765 for (unsigned int j = 0; j < dim; ++j)
3766 s += t1[i][j] * t2[i][j];
3767 return s;
3768}
3769
3770
3783template <int dim, typename Number, typename OtherNumber>
3788{
3789 return scalar_product(t2, t1);
3790}
3791
3792
3807template <typename Number, typename OtherNumber>
3808DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
3813{
3814 tmp[0][0] = t[0][0][0][0] * s[0][0];
3815}
3816
3817
3818
3833template <typename Number, typename OtherNumber>
3834DEAL_II_HOST constexpr inline void
3839{
3840 tmp[0][0] = t[0][0][0][0] * s[0][0];
3841}
3842
3843
3844
3859template <typename Number, typename OtherNumber>
3860DEAL_II_HOST constexpr inline void
3865{
3866 const unsigned int dim = 2;
3867
3868 for (unsigned int i = 0; i < dim; ++i)
3869 for (unsigned int j = i; j < dim; ++j)
3870 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3871 2 * t[i][j][0][1] * s[0][1];
3872}
3873
3874
3875
3890template <typename Number, typename OtherNumber>
3891DEAL_II_HOST constexpr inline void
3896{
3897 const unsigned int dim = 2;
3898
3899 for (unsigned int i = 0; i < dim; ++i)
3900 for (unsigned int j = i; j < dim; ++j)
3901 tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3902 2 * s[0][1] * t[0][1][i][j];
3903}
3904
3905
3906
3921template <typename Number, typename OtherNumber>
3922DEAL_II_HOST constexpr inline void
3927{
3928 const unsigned int dim = 3;
3929
3930 for (unsigned int i = 0; i < dim; ++i)
3931 for (unsigned int j = i; j < dim; ++j)
3932 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3933 t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3934 2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3935}
3936
3937
3938
3953template <typename Number, typename OtherNumber>
3954DEAL_II_HOST constexpr inline void
3959{
3960 const unsigned int dim = 3;
3961
3962 for (unsigned int i = 0; i < dim; ++i)
3963 for (unsigned int j = i; j < dim; ++j)
3964 tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3965 s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3966 2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3967}
3968
3969
3970
3977template <int dim, typename Number, typename OtherNumber>
3978DEAL_II_HOST constexpr Tensor<1,
3979 dim,
3982 const Tensor<1, dim, OtherNumber> & src2)
3983{
3985 for (unsigned int i = 0; i < dim; ++i)
3986 for (unsigned int j = 0; j < dim; ++j)
3987 dest[i] += src1[i][j] * src2[j];
3988 return dest;
3989}
3990
3991
3998template <int dim, typename Number, typename OtherNumber>
3999DEAL_II_HOST constexpr Tensor<1,
4000 dim,
4004{
4005 // this is easy for symmetric tensors:
4006 return src2 * src1;
4007}
4008
4009
4010
4030template <int rank_1,
4031 int rank_2,
4032 int dim,
4033 typename Number,
4034 typename OtherNumber>
4036 typename Tensor<rank_1 + rank_2 - 2,
4037 dim,
4038 typename ProductType<Number, OtherNumber>::type>::tensor_type
4041{
4042 return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4043}
4044
4045
4046
4066template <int rank_1,
4067 int rank_2,
4068 int dim,
4069 typename Number,
4070 typename OtherNumber>
4072 typename Tensor<rank_1 + rank_2 - 2,
4073 dim,
4074 typename ProductType<Number, OtherNumber>::type>::tensor_type
4077{
4078 return Tensor<rank_1, dim, Number>(src1) * src2;
4079}
4080
4081
4082
4092template <int dim, typename Number>
4093inline std::ostream &
4094operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4095{
4096 // make our lives a bit simpler by outputting
4097 // the tensor through the operator for the
4098 // general Tensor class
4100
4101 for (unsigned int i = 0; i < dim; ++i)
4102 for (unsigned int j = 0; j < dim; ++j)
4103 tt[i][j] = t[i][j];
4104
4105 return out << tt;
4106}
4107
4108
4109
4119template <int dim, typename Number>
4120inline std::ostream &
4121operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4122{
4123 // make our lives a bit simpler by outputting
4124 // the tensor through the operator for the
4125 // general Tensor class
4127
4128 for (unsigned int i = 0; i < dim; ++i)
4129 for (unsigned int j = 0; j < dim; ++j)
4130 for (unsigned int k = 0; k < dim; ++k)
4131 for (unsigned int l = 0; l < dim; ++l)
4132 tt[i][j][k][l] = t[i][j][k][l];
4133
4134 return out << tt;
4135}
4136
4137
4139
4140#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)
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)
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
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:106
#define DEAL_II_DEPRECATED
Definition config.h:172
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_CONSTEXPR
Definition config.h:185
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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)
LogStream & operator<<(LogStream &log, const T &t)
Definition logstream.h:406
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:939
static const unsigned int invalid_unsigned_int
Definition types.h:213
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:47
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::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 > invert(const SymmetricTensor< 2, dim, Number > &)
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()