Reference documentation for deal.II version 9.4.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\}}\)
symmetric_tensor.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2005 - 2022 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
43template <int dim, typename Number = double>
46
75template <int dim, typename Number = double>
78
116template <int dim, typename Number = double>
119
120template <int dim, typename Number>
123
124template <int dim, typename Number>
127
128template <int dim2, typename Number>
129constexpr inline DEAL_II_ALWAYS_INLINE Number
131
132template <int dim, typename Number>
135
136template <int dim, typename Number>
139
140
141
142namespace internal
143{
144 // Workaround: The following 4 overloads are necessary to be able to
145 // compile the library with Apple Clang 8 and older. We should remove
146 // these overloads again when we bump the minimal required version to
147 // something later than clang-3.6 / Apple Clang 6.3.
148 template <int rank, int dim, typename T, typename U>
149 struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
150 {
151 using type =
152 SymmetricTensor<rank,
153 dim,
154 std::complex<typename ProductType<T, U>::type>>;
155 };
156
157 template <int rank, int dim, typename T, typename U>
158 struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
159 std::complex<U>>
160 {
161 using type =
162 SymmetricTensor<rank,
163 dim,
164 std::complex<typename ProductType<T, U>::type>>;
165 };
166
167 template <typename T, int rank, int dim, typename U>
168 struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
169 {
170 using type =
171 SymmetricTensor<rank,
172 dim,
173 std::complex<typename ProductType<T, U>::type>>;
174 };
175
176 template <int rank, int dim, typename T, typename U>
177 struct ProductTypeImpl<std::complex<T>,
178 SymmetricTensor<rank, dim, std::complex<U>>>
179 {
180 using type =
181 SymmetricTensor<rank,
182 dim,
183 std::complex<typename ProductType<T, U>::type>>;
184 };
185 // end workaround
186
191 namespace SymmetricTensorImplementation
192 {
197 template <int rank, int dim, typename Number>
198 struct Inverse;
199 } // namespace SymmetricTensorImplementation
200
205 namespace SymmetricTensorAccessors
206 {
214 merge(const TableIndices<2> &previous_indices,
215 const unsigned int new_index,
216 const unsigned int position)
217 {
218 AssertIndexRange(position, 2);
219
220 if (position == 0)
221 return {new_index, numbers::invalid_unsigned_int};
222 else
223 return {previous_indices[0], new_index};
224 }
225
226
227
235 merge(const TableIndices<4> &previous_indices,
236 const unsigned int new_index,
237 const unsigned int position)
238 {
239 AssertIndexRange(position, 4);
240
241 switch (position)
242 {
243 case 0:
244 return {new_index,
248 case 1:
249 return {previous_indices[0],
250 new_index,
253 case 2:
254 return {previous_indices[0],
255 previous_indices[1],
256 new_index,
258 case 3:
259 return {previous_indices[0],
260 previous_indices[1],
261 previous_indices[2],
262 new_index};
263 default:
264 Assert(false, ExcInternalError());
265 return {};
266 }
267 }
268
269
276 template <int rank1,
277 int rank2,
278 int dim,
279 typename Number,
280 typename OtherNumber = Number>
282 {
284 using type =
285 ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
286 };
287
288
295 template <int dim, typename Number, typename OtherNumber>
296 struct double_contraction_result<2, 2, dim, Number, OtherNumber>
297 {
299 };
300
301
302
315 template <int rank, int dim, typename Number>
317
321 template <int dim, typename Number>
322 struct StorageType<2, dim, Number>
323 {
328 static const unsigned int n_independent_components =
329 (dim * dim + dim) / 2;
330
335 };
336
337
338
342 template <int dim, typename Number>
343 struct StorageType<4, dim, Number>
344 {
350 static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
351
355 static const unsigned int n_independent_components =
356 (n_rank2_components *
358
366 };
367
368
369
374 template <int rank, int dim, bool constness, typename Number>
376
383 template <int rank, int dim, typename Number>
384 struct AccessorTypes<rank, dim, true, Number>
385 {
386 using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
387
388 using reference = Number;
389 };
390
397 template <int rank, int dim, typename Number>
398 struct AccessorTypes<rank, dim, false, Number>
399 {
401
402 using reference = Number &;
403 };
404
405
438 template <int rank, int dim, bool constness, int P, typename Number>
440 {
441 public:
445 using reference =
449
450 private:
471
475 constexpr DEAL_II_ALWAYS_INLINE
476 Accessor(const Accessor &) = default;
477
478 public:
482 constexpr Accessor<rank, dim, constness, P - 1, Number>
483 operator[](const unsigned int i);
484
488 constexpr Accessor<rank, dim, constness, P - 1, Number>
489 operator[](const unsigned int i) const;
490
491 private:
497
498 // Declare some other classes as friends. Make sure to work around bugs
499 // in some compilers:
500 template <int, int, typename>
501 friend class ::SymmetricTensor;
502 template <int, int, bool, int, typename>
503 friend class Accessor;
504 friend class ::SymmetricTensor<rank, dim, Number>;
505 friend class Accessor<rank, dim, constness, P + 1, Number>;
506 };
507
508
509
517 template <int rank, int dim, bool constness, typename Number>
518 class Accessor<rank, dim, constness, 1, Number>
519 {
520 public:
524 using reference =
528
529 private:
553
557 constexpr DEAL_II_ALWAYS_INLINE
558 Accessor(const Accessor &) = default;
559
560 public:
564 constexpr reference
565 operator[](const unsigned int);
566
570 constexpr reference
571 operator[](const unsigned int) const;
572
573 private:
579
580 // Declare some other classes as friends. Make sure to work around bugs
581 // in some compilers:
582 template <int, int, typename>
583 friend class ::SymmetricTensor;
584 template <int, int, bool, int, typename>
586 friend class ::SymmetricTensor<rank, dim, Number>;
587 friend class SymmetricTensorAccessors::
588 Accessor<rank, dim, constness, 2, Number>;
589 };
590 } // namespace SymmetricTensorAccessors
591} // namespace internal
592
593
594
667template <int rank_, int dim, typename Number>
669{
670public:
671 static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
672
681 static constexpr unsigned int dimension = dim;
682
686 static const unsigned int rank = rank_;
687
693 static constexpr unsigned int n_independent_components =
695 n_independent_components;
696
700 constexpr DEAL_II_ALWAYS_INLINE
701 SymmetricTensor() = default;
702
716 template <typename OtherNumber>
718
734 constexpr SymmetricTensor(const Number (&array)[n_independent_components]);
735
741 template <typename OtherNumber>
742 constexpr explicit SymmetricTensor(
744
753 DEAL_II_DEPRECATED_EARLY
754 Number *
756
765 DEAL_II_DEPRECATED_EARLY
766 const Number *
767 begin_raw() const;
768
777 DEAL_II_DEPRECATED_EARLY
778 Number *
780
790 DEAL_II_DEPRECATED_EARLY
791 const Number *
792 end_raw() const;
793
800 template <typename OtherNumber>
801 constexpr SymmetricTensor &
803
810 constexpr SymmetricTensor &
811 operator=(const Number &d);
812
817 constexpr operator Tensor<rank_, dim, Number>() const;
818
822 constexpr bool
824
828 constexpr bool
830
834 template <typename OtherNumber>
835 constexpr SymmetricTensor &
837
841 template <typename OtherNumber>
842 constexpr SymmetricTensor &
844
849 template <typename OtherNumber>
850 constexpr SymmetricTensor &
851 operator*=(const OtherNumber &factor);
852
856 template <typename OtherNumber>
857 constexpr SymmetricTensor &
858 operator/=(const OtherNumber &factor);
859
863 constexpr SymmetricTensor
864 operator-() const;
865
918 template <typename OtherNumber>
922
927 template <typename OtherNumber>
931
935 constexpr Number &
937
941 constexpr const Number &
942 operator()(const TableIndices<rank_> &indices) const;
943
948 constexpr internal::SymmetricTensorAccessors::
949 Accessor<rank_, dim, true, rank_ - 1, Number>
950 operator[](const unsigned int row) const;
951
956 constexpr internal::SymmetricTensorAccessors::
957 Accessor<rank_, dim, false, rank_ - 1, Number>
958 operator[](const unsigned int row);
959
965 constexpr const Number &
966 operator[](const TableIndices<rank_> &indices) const;
967
973 constexpr Number &
975
982 constexpr const Number &
983 access_raw_entry(const unsigned int unrolled_index) const;
984
991 constexpr Number &
992 access_raw_entry(const unsigned int unrolled_index);
993
1004 norm() const;
1005
1013 static constexpr unsigned int
1015
1021 static constexpr TableIndices<rank_>
1022 unrolled_to_component_indices(const unsigned int i);
1023
1036 constexpr void
1038
1043 static constexpr std::size_t
1045
1051 template <class Archive>
1052 void
1053 serialize(Archive &ar, const unsigned int version);
1054
1055private:
1061
1065 using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
1066
1071
1072 // Make all other symmetric tensors friends.
1073 template <int, int, typename>
1074 friend class SymmetricTensor;
1075
1076 // Make a few more functions friends.
1077 template <int dim2, typename Number2>
1078 friend constexpr Number2
1080
1081 template <int dim2, typename Number2>
1082 friend DEAL_II_CONSTEXPR Number2
1084
1085 template <int dim2, typename Number2>
1086 friend constexpr SymmetricTensor<2, dim2, Number2>
1088
1089 template <int dim2, typename Number2>
1092
1093 template <int dim2, typename Number2>
1096
1097 template <int dim2, typename Number2>
1100
1101
1102 // Make a few helper classes friends as well.
1104 Inverse<2, dim, Number>;
1105
1107 Inverse<4, dim, Number>;
1108};
1109
1110
1111
1112// ------------------------- inline functions ------------------------
1113
1114#ifndef DOXYGEN
1115
1116// provide declarations for static members
1117template <int rank, int dim, typename Number>
1118const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1119
1120template <int rank_, int dim, typename Number>
1121constexpr unsigned int
1122 SymmetricTensor<rank_, dim, Number>::n_independent_components;
1123
1124namespace internal
1125{
1126 namespace SymmetricTensorAccessors
1127 {
1128 template <int rank_, int dim, bool constness, int P, typename Number>
1129 constexpr DEAL_II_ALWAYS_INLINE
1133 : tensor(tensor)
1135 {}
1136
1137
1138
1139 template <int rank_, int dim, bool constness, int P, typename Number>
1140 constexpr inline DEAL_II_ALWAYS_INLINE
1141 Accessor<rank_, dim, constness, P - 1, Number>
1142 Accessor<rank_, dim, constness, P, Number>::operator[](
1143 const unsigned int i)
1144 {
1145 return Accessor<rank_, dim, constness, P - 1, Number>(
1146 tensor, merge(previous_indices, i, rank_ - P));
1147 }
1148
1149
1150
1151 template <int rank_, int dim, bool constness, int P, typename Number>
1152 constexpr DEAL_II_ALWAYS_INLINE
1153 Accessor<rank_, dim, constness, P - 1, Number>
1154 Accessor<rank_, dim, constness, P, Number>::operator[](
1155 const unsigned int i) const
1156 {
1157 return Accessor<rank_, dim, constness, P - 1, Number>(
1158 tensor, merge(previous_indices, i, rank_ - P));
1159 }
1160
1161
1162
1163 template <int rank_, int dim, bool constness, typename Number>
1164 constexpr DEAL_II_ALWAYS_INLINE
1165 Accessor<rank_, dim, constness, 1, Number>::Accessor(
1166 tensor_type & tensor,
1167 const TableIndices<rank_> &previous_indices)
1168 : tensor(tensor)
1169 , previous_indices(previous_indices)
1170 {}
1171
1172
1173
1174 template <int rank_, int dim, bool constness, typename Number>
1175 constexpr inline DEAL_II_ALWAYS_INLINE
1176 typename Accessor<rank_, dim, constness, 1, Number>::reference
1177 Accessor<rank_, dim, constness, 1, Number>::operator[](
1178 const unsigned int i)
1179 {
1180 return tensor(merge(previous_indices, i, rank_ - 1));
1181 }
1182
1183
1184 template <int rank_, int dim, bool constness, typename Number>
1185 constexpr DEAL_II_ALWAYS_INLINE
1186 typename Accessor<rank_, dim, constness, 1, Number>::reference
1187 Accessor<rank_, dim, constness, 1, Number>::operator[](
1188 const unsigned int i) const
1189 {
1190 return tensor(merge(previous_indices, i, rank_ - 1));
1191 }
1192 } // namespace SymmetricTensorAccessors
1193} // namespace internal
1194
1195
1196
1197template <int rank_, int dim, typename Number>
1198template <typename OtherNumber>
1202{
1203 static_assert(rank == 2, "This function is only implemented for rank==2");
1204 for (unsigned int d = 0; d < dim; ++d)
1205 for (unsigned int e = 0; e < d; ++e)
1206 Assert(t[d][e] == t[e][d],
1207 ExcMessage("The incoming Tensor must be exactly symmetric."));
1208
1209 for (unsigned int d = 0; d < dim; ++d)
1210 data[d] = t[d][d];
1211
1212 for (unsigned int d = 0, c = 0; d < dim; ++d)
1213 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1214 data[dim + c] = t[d][e];
1215}
1216
1217
1218
1219template <int rank_, int dim, typename Number>
1220template <typename OtherNumber>
1221constexpr DEAL_II_ALWAYS_INLINE
1224 : data(initializer.data)
1225{}
1226
1227
1228
1229template <int rank_, int dim, typename Number>
1230constexpr inline DEAL_II_ALWAYS_INLINE
1232 const Number (&array)[n_independent_components])
1233 : data(
1234 *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1235{
1236 // ensure that the reinterpret_cast above actually works
1237 Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1239}
1240
1241
1242
1243template <int rank_, int dim, typename Number>
1244template <typename OtherNumber>
1248{
1249 data = t.data;
1250 return *this;
1251}
1252
1253
1254
1255template <int rank_, int dim, typename Number>
1258{
1260 ExcMessage("Only assignment with zero is allowed"));
1261 (void)d;
1262
1264
1265 return *this;
1266}
1267
1268
1269namespace internal
1270{
1271 namespace SymmetricTensorImplementation
1272 {
1273 template <int dim, typename Number>
1274 constexpr inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1275 convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1276 {
1278
1279 // diagonal entries are stored first
1280 for (unsigned int d = 0; d < dim; ++d)
1281 t[d][d] = s.access_raw_entry(d);
1282
1283 // off-diagonal entries come next, row by row
1284 for (unsigned int d = 0, c = 0; d < dim; ++d)
1285 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1286 {
1287 t[d][e] = s.access_raw_entry(dim + c);
1288 t[e][d] = s.access_raw_entry(dim + c);
1289 }
1290 return t;
1291 }
1292
1293
1294 template <int dim, typename Number>
1295 constexpr ::Tensor<4, dim, Number>
1296 convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1297 {
1298 // utilize the symmetry properties of SymmetricTensor<4,dim>
1299 // discussed in the class documentation to avoid accessing all
1300 // independent elements of the input tensor more than once
1302
1303 for (unsigned int i = 0; i < dim; ++i)
1304 for (unsigned int j = i; j < dim; ++j)
1305 for (unsigned int k = 0; k < dim; ++k)
1306 for (unsigned int l = k; l < dim; ++l)
1307 t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1308 t[TableIndices<4>(j, i, k, l)] =
1309 t[TableIndices<4>(j, i, l, k)] =
1310 st[TableIndices<4>(i, j, k, l)];
1311
1312 return t;
1313 }
1314
1315
1316 template <typename Number>
1317 struct Inverse<2, 1, Number>
1318 {
1319 constexpr static inline DEAL_II_ALWAYS_INLINE
1320 ::SymmetricTensor<2, 1, Number>
1321 value(const ::SymmetricTensor<2, 1, Number> &t)
1322 {
1324
1325 tmp[0][0] = 1.0 / t[0][0];
1326
1327 return tmp;
1328 }
1329 };
1330
1331
1332 template <typename Number>
1333 struct Inverse<2, 2, Number>
1334 {
1335 constexpr static inline DEAL_II_ALWAYS_INLINE
1336 ::SymmetricTensor<2, 2, Number>
1337 value(const ::SymmetricTensor<2, 2, Number> &t)
1338 {
1340
1341 // Sympy result: ([
1342 // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1343 // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1344 const TableIndices<2> idx_00(0, 0);
1345 const TableIndices<2> idx_01(0, 1);
1346 const TableIndices<2> idx_11(1, 1);
1347 const Number inv_det_t =
1348 1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1349 tmp[idx_00] = t[idx_11];
1350 tmp[idx_01] = -t[idx_01];
1351 tmp[idx_11] = t[idx_00];
1352 tmp *= inv_det_t;
1353
1354 return tmp;
1355 }
1356 };
1357
1358
1359 template <typename Number>
1360 struct Inverse<2, 3, Number>
1361 {
1362 constexpr static ::SymmetricTensor<2, 3, Number>
1363 value(const ::SymmetricTensor<2, 3, Number> &t)
1364 {
1366
1367 // Sympy result: ([
1368 // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1369 // 2*t01*t02*t12 - t02**2*t11),
1370 // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1371 // 2*t01*t02*t12 - t02**2*t11),
1372 // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1373 // 2*t01*t02*t12 - t02**2*t11)],
1374 // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1375 // 2*t01*t02*t12 - t02**2*t11),
1376 // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1377 // 2*t01*t02*t12 - t02**2*t11),
1378 // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1379 // 2*t01*t02*t12 + t02**2*t11)],
1380 // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1381 // 2*t01*t02*t12 - t02**2*t11),
1382 // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1383 // 2*t01*t02*t12 + t02**2*t11),
1384 // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1385 // 2*t01*t02*t12 + t02**2*t11)] ])
1386 //
1387 // =
1388 //
1389 // [ (t11*t22 - t12**2)/det_t,
1390 // (-t01*t22 + t02*t12)/det_t,
1391 // (t01*t12 - t02*t11)/det_t],
1392 // [ (-t01*t22 + t02*t12)/det_t,
1393 // (t00*t22 - t02**2)/det_t,
1394 // (-t00*t12 + t01*t02)/det_t],
1395 // [ (t01*t12 - t02*t11)/det_t,
1396 // (-t00*t12 + t01*t02)/det_t,
1397 // (t00*t11 - t01**2)/det_t] ])
1398 //
1399 // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1400 // 2*t01*t02*t12 - t02**2*t11)
1401 const TableIndices<2> idx_00(0, 0);
1402 const TableIndices<2> idx_01(0, 1);
1403 const TableIndices<2> idx_02(0, 2);
1404 const TableIndices<2> idx_11(1, 1);
1405 const TableIndices<2> idx_12(1, 2);
1406 const TableIndices<2> idx_22(2, 2);
1407 const Number inv_det_t =
1408 1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1409 t[idx_00] * t[idx_12] * t[idx_12] -
1410 t[idx_01] * t[idx_01] * t[idx_22] +
1411 2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1412 t[idx_02] * t[idx_02] * t[idx_11]);
1413 tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1414 tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1415 tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1416 tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1417 tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1418 tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1419 tmp *= inv_det_t;
1420
1421 return tmp;
1422 }
1423 };
1424
1425
1426 template <typename Number>
1427 struct Inverse<4, 1, Number>
1428 {
1429 constexpr static inline ::SymmetricTensor<4, 1, Number>
1430 value(const ::SymmetricTensor<4, 1, Number> &t)
1431 {
1433 tmp.data[0][0] = 1.0 / t.data[0][0];
1434 return tmp;
1435 }
1436 };
1437
1438
1439 template <typename Number>
1440 struct Inverse<4, 2, Number>
1441 {
1442 constexpr static inline ::SymmetricTensor<4, 2, Number>
1443 value(const ::SymmetricTensor<4, 2, Number> &t)
1444 {
1446
1447 // Inverting this tensor is a little more complicated than necessary,
1448 // since we store the data of 't' as a 3x3 matrix t.data, but the
1449 // product between a rank-4 and a rank-2 tensor is really not the
1450 // product between this matrix and the 3-vector of a rhs, but rather
1451 //
1452 // B.vec = t.data * mult * A.vec
1453 //
1454 // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1455 // capture the fact that we need to add up both the c_ij12*a_12 and the
1456 // c_ij21*a_21 terms.
1457 //
1458 // In addition, in this scheme, the identity tensor has the matrix
1459 // representation mult^-1.
1460 //
1461 // The inverse of 't' therefore has the matrix representation
1462 //
1463 // inv.data = mult^-1 * t.data^-1 * mult^-1
1464 //
1465 // in order to compute it, let's first compute the inverse of t.data and
1466 // put it into tmp.data; at the end of the function we then scale the
1467 // last row and column of the inverse by 1/2, corresponding to the left
1468 // and right multiplication with mult^-1.
1469 const Number t4 = t.data[0][0] * t.data[1][1],
1470 t6 = t.data[0][0] * t.data[1][2],
1471 t8 = t.data[0][1] * t.data[1][0],
1472 t00 = t.data[0][2] * t.data[1][0],
1473 t01 = t.data[0][1] * t.data[2][0],
1474 t04 = t.data[0][2] * t.data[2][0],
1475 t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1476 t8 * t.data[2][2] + t00 * t.data[2][1] +
1477 t01 * t.data[1][2] - t04 * t.data[1][1]);
1478 tmp.data[0][0] =
1479 (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1480 tmp.data[0][1] =
1481 -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1482 tmp.data[0][2] =
1483 -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1484 tmp.data[1][0] =
1485 -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1486 tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1487 tmp.data[1][2] = -(t6 - t00) * t07;
1488 tmp.data[2][0] =
1489 -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1490 tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1491 tmp.data[2][2] = (t4 - t8) * t07;
1492
1493 // scale last row and column as mentioned
1494 // above
1495 tmp.data[2][0] /= 2;
1496 tmp.data[2][1] /= 2;
1497 tmp.data[0][2] /= 2;
1498 tmp.data[1][2] /= 2;
1499 tmp.data[2][2] /= 4;
1500
1501 return tmp;
1502 }
1503 };
1504
1505
1506 template <typename Number>
1507 struct Inverse<4, 3, Number>
1508 {
1509 static ::SymmetricTensor<4, 3, Number>
1510 value(const ::SymmetricTensor<4, 3, Number> &t)
1511 {
1513
1514 // This function follows the exact same scheme as the 2d case, except
1515 // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1516 // Instead, we use the Gauss-Jordan algorithm implemented for
1517 // FullMatrix. For historical reasons the following code is copied from
1518 // there, with the tangential benefit that we do not need to copy the
1519 // tensor entries to and from the FullMatrix.
1520 const unsigned int N = 6;
1521
1522 // First get an estimate of the size of the elements of this matrix,
1523 // for later checks whether the pivot element is large enough, or
1524 // whether we have to fear that the matrix is not regular.
1525 Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1526 for (unsigned int i = 0; i < N; ++i)
1527 diagonal_sum += std::fabs(tmp.data[i][i]);
1528 const Number typical_diagonal_element =
1529 diagonal_sum / static_cast<double>(N);
1530 (void)typical_diagonal_element;
1531
1532 unsigned int p[N];
1533 for (unsigned int i = 0; i < N; ++i)
1534 p[i] = i;
1535
1536 for (unsigned int j = 0; j < N; ++j)
1537 {
1538 // Pivot search: search that part of the line on and right of the
1539 // diagonal for the largest element.
1540 Number max = std::fabs(tmp.data[j][j]);
1541 unsigned int r = j;
1542 for (unsigned int i = j + 1; i < N; ++i)
1543 if (std::fabs(tmp.data[i][j]) > max)
1544 {
1545 max = std::fabs(tmp.data[i][j]);
1546 r = i;
1547 }
1548
1549 // Check whether the pivot is too small
1550 Assert(max > 1.e-16 * typical_diagonal_element,
1551 ExcMessage("This tensor seems to be noninvertible"));
1552
1553 // Row interchange
1554 if (r > j)
1555 {
1556 for (unsigned int k = 0; k < N; ++k)
1557 std::swap(tmp.data[j][k], tmp.data[r][k]);
1558
1559 std::swap(p[j], p[r]);
1560 }
1561
1562 // Transformation
1563 const Number hr = 1. / tmp.data[j][j];
1564 tmp.data[j][j] = hr;
1565 for (unsigned int k = 0; k < N; ++k)
1566 {
1567 if (k == j)
1568 continue;
1569 for (unsigned int i = 0; i < N; ++i)
1570 {
1571 if (i == j)
1572 continue;
1573 tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1574 }
1575 }
1576 for (unsigned int i = 0; i < N; ++i)
1577 {
1578 tmp.data[i][j] *= hr;
1579 tmp.data[j][i] *= -hr;
1580 }
1581 tmp.data[j][j] = hr;
1582 }
1583
1584 // Column interchange
1585 Number hv[N];
1586 for (unsigned int i = 0; i < N; ++i)
1587 {
1588 for (unsigned int k = 0; k < N; ++k)
1589 hv[p[k]] = tmp.data[i][k];
1590 for (unsigned int k = 0; k < N; ++k)
1591 tmp.data[i][k] = hv[k];
1592 }
1593
1594 // Scale rows and columns. The mult matrix
1595 // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1596 for (unsigned int i = 3; i < 6; ++i)
1597 for (unsigned int j = 0; j < 3; ++j)
1598 tmp.data[i][j] /= 2;
1599
1600 for (unsigned int i = 0; i < 3; ++i)
1601 for (unsigned int j = 3; j < 6; ++j)
1602 tmp.data[i][j] /= 2;
1603
1604 for (unsigned int i = 3; i < 6; ++i)
1605 for (unsigned int j = 3; j < 6; ++j)
1606 tmp.data[i][j] /= 4;
1607
1608 return tmp;
1609 }
1610 };
1611
1612 } // namespace SymmetricTensorImplementation
1613} // namespace internal
1614
1615
1616
1617template <int rank_, int dim, typename Number>
1618constexpr DEAL_II_ALWAYS_INLINE
1620 const
1621{
1622 return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1623}
1624
1625
1626
1627template <int rank_, int dim, typename Number>
1628constexpr bool
1631{
1632 return data == t.data;
1633}
1634
1635
1636
1637template <int rank_, int dim, typename Number>
1638constexpr bool
1641{
1642 return data != t.data;
1643}
1644
1645
1646
1647template <int rank_, int dim, typename Number>
1648template <typename OtherNumber>
1652{
1653 data += t.data;
1654 return *this;
1655}
1656
1657
1658
1659template <int rank_, int dim, typename Number>
1660template <typename OtherNumber>
1664{
1665 data -= t.data;
1666 return *this;
1667}
1668
1669
1670
1671template <int rank_, int dim, typename Number>
1672template <typename OtherNumber>
1675{
1676 data *= d;
1677 return *this;
1678}
1679
1680
1681
1682template <int rank_, int dim, typename Number>
1683template <typename OtherNumber>
1686{
1687 data /= d;
1688 return *this;
1689}
1690
1691
1692
1693template <int rank_, int dim, typename Number>
1696{
1697 SymmetricTensor tmp = *this;
1698 tmp.data = -tmp.data;
1699 return tmp;
1700}
1701
1702
1703
1704template <int rank_, int dim, typename Number>
1705constexpr inline DEAL_II_ALWAYS_INLINE void
1707{
1708 data.clear();
1709}
1710
1711
1712
1713template <int rank_, int dim, typename Number>
1714constexpr std::size_t
1716{
1717 // all memory consists of statically allocated memory of the current
1718 // object, no pointers
1720}
1721
1722
1723
1724namespace internal
1725{
1726 template <int dim, typename Number, typename OtherNumber = Number>
1730 perform_double_contraction(
1731 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1732 base_tensor_type &data,
1733 const typename SymmetricTensorAccessors::
1734 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1735 {
1736 using result_type = typename SymmetricTensorAccessors::
1738
1739 switch (dim)
1740 {
1741 case 1:
1742 return data[0] * sdata[0];
1743 default:
1744 // Start with the non-diagonal part to avoid some multiplications by
1745 // 2.
1746
1747 result_type sum = data[dim] * sdata[dim];
1748 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1749 sum += data[d] * sdata[d];
1750 sum += sum; // sum = sum * 2.;
1751
1752 // Now add the contributions from the diagonal
1753 for (unsigned int d = 0; d < dim; ++d)
1754 sum += data[d] * sdata[d];
1755 return sum;
1756 }
1757 }
1758
1759
1760
1761 template <int dim, typename Number, typename OtherNumber = Number>
1765 perform_double_contraction(
1766 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1767 base_tensor_type &data,
1768 const typename SymmetricTensorAccessors::
1769 StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1770 {
1771 using result_type = typename SymmetricTensorAccessors::
1773 using value_type = typename SymmetricTensorAccessors::
1775
1776 const unsigned int data_dim = SymmetricTensorAccessors::
1777 StorageType<2, dim, value_type>::n_independent_components;
1778 value_type tmp[data_dim]{};
1779 for (unsigned int i = 0; i < data_dim; ++i)
1780 tmp[i] =
1781 perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1782 return result_type(tmp);
1783 }
1784
1785
1786
1787 template <int dim, typename Number, typename OtherNumber = Number>
1789 typename SymmetricTensorAccessors::StorageType<
1790 2,
1791 dim,
1794 base_tensor_type
1795 perform_double_contraction(
1796 const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1797 base_tensor_type &data,
1798 const typename SymmetricTensorAccessors::
1799 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1800 {
1801 using value_type = typename SymmetricTensorAccessors::
1803 using base_tensor_type = typename SymmetricTensorAccessors::
1804 StorageType<2, dim, value_type>::base_tensor_type;
1805
1806 base_tensor_type tmp;
1807 for (unsigned int i = 0; i < tmp.dimension; ++i)
1808 {
1809 // Start with the non-diagonal part
1810 value_type sum = data[dim] * sdata[dim][i];
1811 for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1812 sum += data[d] * sdata[d][i];
1813 sum += sum; // sum = sum * 2.;
1814
1815 // Now add the contributions from the diagonal
1816 for (unsigned int d = 0; d < dim; ++d)
1817 sum += data[d] * sdata[d][i];
1818 tmp[i] = sum;
1819 }
1820 return tmp;
1821 }
1822
1823
1824
1825 template <int dim, typename Number, typename OtherNumber = Number>
1827 typename SymmetricTensorAccessors::StorageType<
1828 4,
1829 dim,
1832 base_tensor_type
1833 perform_double_contraction(
1834 const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1835 base_tensor_type &data,
1836 const typename SymmetricTensorAccessors::
1837 StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1838 {
1839 using value_type = typename SymmetricTensorAccessors::
1841 using base_tensor_type = typename SymmetricTensorAccessors::
1842 StorageType<4, dim, value_type>::base_tensor_type;
1843
1844 const unsigned int data_dim = SymmetricTensorAccessors::
1845 StorageType<2, dim, value_type>::n_independent_components;
1846 base_tensor_type tmp;
1847 for (unsigned int i = 0; i < data_dim; ++i)
1848 for (unsigned int j = 0; j < data_dim; ++j)
1849 {
1850 // Start with the non-diagonal part
1851 for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1852 tmp[i][j] += data[i][d] * sdata[d][j];
1853 tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1854
1855 // Now add the contributions from the diagonal
1856 for (unsigned int d = 0; d < dim; ++d)
1857 tmp[i][j] += data[i][d] * sdata[d][j];
1858 }
1859 return tmp;
1860 }
1861
1862} // end of namespace internal
1863
1864
1865
1866template <int rank_, int dim, typename Number>
1867template <typename OtherNumber>
1873{
1874 // need to have two different function calls
1875 // because a scalar and rank-2 tensor are not
1876 // the same data type (see internal function
1877 // above)
1878 return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1879 s.data);
1880}
1881
1882
1883
1884template <int rank_, int dim, typename Number>
1885template <typename OtherNumber>
1890{
1893 tmp.data =
1894 internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1895 s.data);
1896 return tmp;
1897}
1898
1899
1900
1901// internal namespace to switch between the
1902// access of different tensors. There used to
1903// be explicit instantiations before for
1904// different ranks and dimensions, but since
1905// we now allow for templates on the data
1906// type, and since we cannot partially
1907// specialize the implementation, this got
1908// into a separate namespace
1909namespace internal
1910{
1911 // The variables within this struct will be referenced in the next functions.
1912 // It is a workaround that allows returning a reference to a static variable
1913 // while allowing constexpr evaluation of the function.
1914 // It has to be defined outside the function because constexpr functions
1915 // cannot define static variables.
1916 // A similar struct has also been defined in tensor.h
1917 template <typename Type>
1918 struct Uninitialized
1919 {
1920 static Type value;
1921 };
1922
1923 template <typename Type>
1924 Type Uninitialized<Type>::value;
1925
1926 template <int dim, typename Number>
1927 constexpr inline DEAL_II_ALWAYS_INLINE Number &
1928 symmetric_tensor_access(const TableIndices<2> &indices,
1929 typename SymmetricTensorAccessors::
1930 StorageType<2, dim, Number>::base_tensor_type &data)
1931 {
1932 // 1d is very simple and done first
1933 if (dim == 1)
1934 return data[0];
1935
1936 // first treat the main diagonal elements, which are stored consecutively
1937 // at the beginning
1938 if (indices[0] == indices[1])
1939 return data[indices[0]];
1940
1941 // the rest is messier and requires a few switches.
1942 switch (dim)
1943 {
1944 case 2:
1945 // at least for the 2x2 case it is reasonably simple
1946 Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1947 ((indices[0] == 0) && (indices[1] == 1)),
1949 return data[2];
1950
1951 default:
1952 // to do the rest, sort our indices before comparing
1953 {
1954 TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1955 std::max(indices[0], indices[1]));
1956 for (unsigned int d = 0, c = 0; d < dim; ++d)
1957 for (unsigned int e = d + 1; e < dim; ++e, ++c)
1958 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1959 return data[dim + c];
1960 Assert(false, ExcInternalError());
1961 }
1962 }
1963
1964 // The code should never reach there.
1965 // Returns a dummy reference to a dummy variable just to make the
1966 // compiler happy.
1967 return Uninitialized<Number>::value;
1968 }
1969
1970
1971
1972 template <int dim, typename Number>
1973 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1974 symmetric_tensor_access(const TableIndices<2> &indices,
1975 const typename SymmetricTensorAccessors::
1976 StorageType<2, dim, Number>::base_tensor_type &data)
1977 {
1978 // 1d is very simple and done first
1979 if (dim == 1)
1980 return data[0];
1981
1982 // first treat the main diagonal elements, which are stored consecutively
1983 // at the beginning
1984 if (indices[0] == indices[1])
1985 return data[indices[0]];
1986
1987 // the rest is messier and requires a few switches.
1988 switch (dim)
1989 {
1990 case 2:
1991 // at least for the 2x2 case it is reasonably simple
1992 Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1993 ((indices[0] == 0) && (indices[1] == 1)),
1995 return data[2];
1996
1997 default:
1998 // to do the rest, sort our indices before comparing
1999 {
2000 TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2001 std::max(indices[0], indices[1]));
2002 for (unsigned int d = 0, c = 0; d < dim; ++d)
2003 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2004 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2005 return data[dim + c];
2006 Assert(false, ExcInternalError());
2007 }
2008 }
2009
2010 // The code should never reach there.
2011 // Returns a dummy reference to a dummy variable just to make the
2012 // compiler happy.
2013 return Uninitialized<Number>::value;
2014 }
2015
2016
2017
2018 template <int dim, typename Number>
2019 constexpr inline Number &
2020 symmetric_tensor_access(const TableIndices<4> &indices,
2021 typename SymmetricTensorAccessors::
2022 StorageType<4, dim, Number>::base_tensor_type &data)
2023 {
2024 switch (dim)
2025 {
2026 case 1:
2027 return data[0][0];
2028
2029 case 2:
2030 // each entry of the tensor can be thought of as an entry in a
2031 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2032 // rank-2 tensors. this is the format in which we store rank-4
2033 // tensors. determine which position the present entry is
2034 // stored in
2035 {
2036 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2037 return data[base_index[indices[0]][indices[1]]]
2038 [base_index[indices[2]][indices[3]]];
2039 }
2040 case 3:
2041 // each entry of the tensor can be thought of as an entry in a
2042 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2043 // rank-2 tensors. this is the format in which we store rank-4
2044 // tensors. determine which position the present entry is
2045 // stored in
2046 {
2047 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2048 {3, 1, 5},
2049 {4, 5, 2}};
2050 return data[base_index[indices[0]][indices[1]]]
2051 [base_index[indices[2]][indices[3]]];
2052 }
2053
2054 default:
2055 Assert(false, ExcNotImplemented());
2056 }
2057
2058 // The code should never reach there.
2059 // Returns a dummy reference to a dummy variable just to make the
2060 // compiler happy.
2061 return Uninitialized<Number>::value;
2062 }
2063
2064
2065 template <int dim, typename Number>
2066 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2067 symmetric_tensor_access(const TableIndices<4> &indices,
2068 const typename SymmetricTensorAccessors::
2069 StorageType<4, dim, Number>::base_tensor_type &data)
2070 {
2071 switch (dim)
2072 {
2073 case 1:
2074 return data[0][0];
2075
2076 case 2:
2077 // each entry of the tensor can be thought of as an entry in a
2078 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2079 // rank-2 tensors. this is the format in which we store rank-4
2080 // tensors. determine which position the present entry is
2081 // stored in
2082 {
2083 constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2084 return data[base_index[indices[0]][indices[1]]]
2085 [base_index[indices[2]][indices[3]]];
2086 }
2087 case 3:
2088 // each entry of the tensor can be thought of as an entry in a
2089 // matrix that maps the rolled-out rank-2 tensors into rolled-out
2090 // rank-2 tensors. this is the format in which we store rank-4
2091 // tensors. determine which position the present entry is
2092 // stored in
2093 {
2094 constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2095 {3, 1, 5},
2096 {4, 5, 2}};
2097 return data[base_index[indices[0]][indices[1]]]
2098 [base_index[indices[2]][indices[3]]];
2099 }
2100
2101 default:
2102 Assert(false, ExcNotImplemented());
2103 }
2104
2105 // The code should never reach there.
2106 // Returns a dummy reference to a dummy variable just to make the
2107 // compiler happy.
2108 return Uninitialized<Number>::value;
2109 }
2110
2111} // end of namespace internal
2112
2113
2114
2115template <int rank_, int dim, typename Number>
2116constexpr inline DEAL_II_ALWAYS_INLINE Number &
2118 const TableIndices<rank_> &indices)
2119{
2120 for (unsigned int r = 0; r < rank; ++r)
2121 AssertIndexRange(indices[r], dimension);
2122 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2123}
2124
2125
2126
2127template <int rank_, int dim, typename Number>
2128constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2130 const TableIndices<rank_> &indices) const
2131{
2132 for (unsigned int r = 0; r < rank; ++r)
2133 AssertIndexRange(indices[r], dimension);
2134 return internal::symmetric_tensor_access<dim, Number>(indices, data);
2135}
2136
2137
2138
2139namespace internal
2140{
2141 namespace SymmetricTensorImplementation
2142 {
2143 template <int rank_>
2144 constexpr TableIndices<rank_>
2145 get_partially_filled_indices(const unsigned int row,
2146 const std::integral_constant<int, 2> &)
2147 {
2149 }
2150
2151
2152 template <int rank_>
2153 constexpr TableIndices<rank_>
2154 get_partially_filled_indices(const unsigned int row,
2155 const std::integral_constant<int, 4> &)
2156 {
2157 return TableIndices<rank_>(row,
2161 }
2162 } // namespace SymmetricTensorImplementation
2163} // namespace internal
2164
2165
2166template <int rank_, int dim, typename Number>
2167constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2168 Accessor<rank_, dim, true, rank_ - 1, Number>
2169 SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2170{
2171 return internal::SymmetricTensorAccessors::
2172 Accessor<rank_, dim, true, rank_ - 1, Number>(
2173 *this,
2174 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2175 rank_>(row, std::integral_constant<int, rank_>()));
2176}
2177
2178
2179
2180template <int rank_, int dim, typename Number>
2181constexpr inline DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2182 Accessor<rank_, dim, false, rank_ - 1, Number>
2184{
2185 return internal::SymmetricTensorAccessors::
2186 Accessor<rank_, dim, false, rank_ - 1, Number>(
2187 *this,
2188 internal::SymmetricTensorImplementation::get_partially_filled_indices<
2189 rank_>(row, std::integral_constant<int, rank_>()));
2190}
2191
2192
2193
2194template <int rank_, int dim, typename Number>
2195constexpr DEAL_II_ALWAYS_INLINE const Number &
2197 const TableIndices<rank_> &indices) const
2198{
2199 return operator()(indices);
2200}
2201
2202
2203
2204template <int rank_, int dim, typename Number>
2205constexpr inline DEAL_II_ALWAYS_INLINE Number &
2207 const TableIndices<rank_> &indices)
2208{
2209 return operator()(indices);
2210}
2211
2212
2213
2214template <int rank_, int dim, typename Number>
2215inline Number *
2217{
2218 return std::addressof(this->access_raw_entry(0));
2219}
2220
2221
2222
2223template <int rank_, int dim, typename Number>
2224inline const Number *
2226{
2227 return std::addressof(this->access_raw_entry(0));
2228}
2229
2230
2231
2232template <int rank_, int dim, typename Number>
2233inline Number *
2235{
2236 return begin_raw() + n_independent_components;
2237}
2238
2239
2240
2241template <int rank_, int dim, typename Number>
2242inline const Number *
2244{
2245 return begin_raw() + n_independent_components;
2246}
2247
2248
2249
2250namespace internal
2251{
2252 namespace SymmetricTensorImplementation
2253 {
2254 template <int dim, typename Number>
2255 constexpr unsigned int
2256 entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2257 const unsigned int index)
2258 {
2259 return index;
2260 }
2261
2262
2263 template <int dim, typename Number>
2264 constexpr ::TableIndices<2>
2265 entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2266 const unsigned int index)
2267 {
2270 }
2271
2272 } // namespace SymmetricTensorImplementation
2273} // namespace internal
2274
2275
2276
2277template <int rank_, int dim, typename Number>
2278constexpr inline const Number &
2280 const unsigned int index) const
2281{
2282 AssertIndexRange(index, n_independent_components);
2283 return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2284 index)];
2285}
2286
2287
2288
2289template <int rank_, int dim, typename Number>
2290constexpr inline Number &
2292{
2293 AssertIndexRange(index, n_independent_components);
2294 return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2295 index)];
2296}
2297
2298
2299
2300namespace internal
2301{
2302 template <int dim, typename Number>
2303 constexpr inline typename numbers::NumberTraits<Number>::real_type
2304 compute_norm(const typename SymmetricTensorAccessors::
2305 StorageType<2, dim, Number>::base_tensor_type &data)
2306 {
2307 switch (dim)
2308 {
2309 case 1:
2310 return numbers::NumberTraits<Number>::abs(data[0]);
2311
2312 case 2:
2313 return std::sqrt(
2317
2318 case 3:
2319 return std::sqrt(
2326
2327 default:
2328 {
2329 typename numbers::NumberTraits<Number>::real_type return_value =
2331
2332 for (unsigned int d = 0; d < dim; ++d)
2333 return_value +=
2335 for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2336 return_value +=
2338
2339 return std::sqrt(return_value);
2340 }
2341 }
2342 }
2343
2344
2345
2346 template <int dim, typename Number>
2347 constexpr inline typename numbers::NumberTraits<Number>::real_type
2348 compute_norm(const typename SymmetricTensorAccessors::
2349 StorageType<4, dim, Number>::base_tensor_type &data)
2350 {
2351 switch (dim)
2352 {
2353 case 1:
2354 return numbers::NumberTraits<Number>::abs(data[0][0]);
2355
2356 default:
2357 {
2358 typename numbers::NumberTraits<Number>::real_type return_value =
2360
2361 const unsigned int n_independent_components = data.dimension;
2362
2363 for (unsigned int i = 0; i < dim; ++i)
2364 for (unsigned int j = 0; j < dim; ++j)
2365 return_value +=
2367 for (unsigned int i = 0; i < dim; ++i)
2368 for (unsigned int j = dim; j < n_independent_components; ++j)
2369 return_value +=
2371 for (unsigned int i = dim; i < n_independent_components; ++i)
2372 for (unsigned int j = 0; j < dim; ++j)
2373 return_value +=
2375 for (unsigned int i = dim; i < n_independent_components; ++i)
2376 for (unsigned int j = dim; j < n_independent_components; ++j)
2377 return_value +=
2379
2380 return std::sqrt(return_value);
2381 }
2382 }
2383 }
2384
2385} // end of namespace internal
2386
2387
2388
2389template <int rank_, int dim, typename Number>
2392{
2393 return internal::compute_norm<dim, Number>(data);
2394}
2395
2396
2397
2398namespace internal
2399{
2400 namespace SymmetricTensorImplementation
2401 {
2402 // a function to do the unrolling from a set of indices to a
2403 // scalar index into the array in which we store the elements of
2404 // a symmetric tensor
2405 //
2406 // this function is for rank-2 tensors
2407 template <int dim>
2408 constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
2410 {
2411 AssertIndexRange(indices[0], dim);
2412 AssertIndexRange(indices[1], dim);
2413
2414 switch (dim)
2415 {
2416 case 1:
2417 {
2418 return 0;
2419 }
2420
2421 case 2:
2422 {
2423 constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2424 return table[indices[0]][indices[1]];
2425 }
2426
2427 case 3:
2428 {
2429 constexpr unsigned int table[3][3] = {{0, 3, 4},
2430 {3, 1, 5},
2431 {4, 5, 2}};
2432 return table[indices[0]][indices[1]];
2433 }
2434
2435 case 4:
2436 {
2437 constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2438 {4, 1, 7, 8},
2439 {5, 7, 2, 9},
2440 {6, 8, 9, 3}};
2441 return table[indices[0]][indices[1]];
2442 }
2443
2444 default:
2445 // for the remainder, manually figure out the numbering
2446 {
2447 if (indices[0] == indices[1])
2448 return indices[0];
2449
2450 TableIndices<2> sorted_indices(indices);
2451 sorted_indices.sort();
2452
2453 for (unsigned int d = 0, c = 0; d < dim; ++d)
2454 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2455 if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2456 return dim + c;
2457
2458 // should never get here:
2459 Assert(false, ExcInternalError());
2460 return 0;
2461 }
2462 }
2463 }
2464
2465 // a function to do the unrolling from a set of indices to a
2466 // scalar index into the array in which we store the elements of
2467 // a symmetric tensor
2468 //
2469 // this function is for tensors of ranks not already handled
2470 // above
2471 template <int dim, int rank_>
2472 constexpr inline unsigned int
2474 {
2475 (void)indices;
2476 Assert(false, ExcNotImplemented());
2478 }
2479 } // namespace SymmetricTensorImplementation
2480} // namespace internal
2481
2482
2483template <int rank_, int dim, typename Number>
2484constexpr unsigned int
2486 const TableIndices<rank_> &indices)
2487{
2488 return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2489 dim>(indices);
2490}
2491
2492
2493
2494namespace internal
2495{
2496 namespace SymmetricTensorImplementation
2497 {
2498 // a function to do the inverse of the unrolling from a set of
2499 // indices to a scalar index into the array in which we store
2500 // the elements of a symmetric tensor. in other words, it goes
2501 // from the scalar index into the array to a set of indices of
2502 // the tensor
2503 //
2504 // this function is for rank-2 tensors
2505 template <int dim>
2506 constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2507 unrolled_to_component_indices(const unsigned int i,
2508 const std::integral_constant<int, 2> &)
2509 {
2510 Assert(
2513 i,
2514 0,
2516 switch (dim)
2517 {
2518 case 1:
2519 {
2520 return {0, 0};
2521 }
2522
2523 case 2:
2524 {
2525 const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2526 TableIndices<2>(1, 1),
2527 TableIndices<2>(0, 1)};
2528 return table[i];
2529 }
2530
2531 case 3:
2532 {
2533 const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2534 TableIndices<2>(1, 1),
2535 TableIndices<2>(2, 2),
2536 TableIndices<2>(0, 1),
2537 TableIndices<2>(0, 2),
2538 TableIndices<2>(1, 2)};
2539 return table[i];
2540 }
2541
2542 default:
2543 if (i < dim)
2544 return {i, i};
2545
2546 for (unsigned int d = 0, c = dim; d < dim; ++d)
2547 for (unsigned int e = d + 1; e < dim; ++e, ++c)
2548 if (c == i)
2549 return {d, e};
2550
2551 // should never get here:
2552 Assert(false, ExcInternalError());
2553 return {0, 0};
2554 }
2555 }
2556
2557 // a function to do the inverse of the unrolling from a set of
2558 // indices to a scalar index into the array in which we store
2559 // the elements of a symmetric tensor. in other words, it goes
2560 // from the scalar index into the array to a set of indices of
2561 // the tensor
2562 //
2563 // this function is for tensors of a rank not already handled
2564 // above
2565 template <int dim, int rank_>
2566 constexpr inline
2567 typename std::enable_if<rank_ != 2, TableIndices<rank_>>::type
2568 unrolled_to_component_indices(const unsigned int i,
2569 const std::integral_constant<int, rank_> &)
2570 {
2571 (void)i;
2572 Assert(
2573 (i <
2575 ExcIndexRange(i,
2576 0,
2578 n_independent_components));
2579 Assert(false, ExcNotImplemented());
2580 return TableIndices<rank_>();
2581 }
2582
2583 } // namespace SymmetricTensorImplementation
2584} // namespace internal
2585
2586template <int rank_, int dim, typename Number>
2589 const unsigned int i)
2590{
2591 return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2592 dim>(i, std::integral_constant<int, rank_>());
2593}
2594
2595
2596
2597template <int rank_, int dim, typename Number>
2598template <class Archive>
2599inline void
2600SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2601{
2602 ar &data;
2603}
2604
2605
2606#endif // DOXYGEN
2607
2608/* ----------------- Non-member functions operating on tensors. ------------ */
2609
2610
2623template <int rank_, int dim, typename Number, typename OtherNumber>
2624constexpr inline DEAL_II_ALWAYS_INLINE
2628{
2630 tmp = left;
2631 tmp += right;
2632 return tmp;
2633}
2634
2635
2648template <int rank_, int dim, typename Number, typename OtherNumber>
2649constexpr inline DEAL_II_ALWAYS_INLINE
2653{
2655 tmp = left;
2656 tmp -= right;
2657 return tmp;
2658}
2659
2660
2668template <int rank_, int dim, typename Number, typename OtherNumber>
2669constexpr DEAL_II_ALWAYS_INLINE
2672 const Tensor<rank_, dim, OtherNumber> & right)
2673{
2674 return Tensor<rank_, dim, Number>(left) + right;
2675}
2676
2677
2685template <int rank_, int dim, typename Number, typename OtherNumber>
2686constexpr DEAL_II_ALWAYS_INLINE
2690{
2691 return left + Tensor<rank_, dim, OtherNumber>(right);
2692}
2693
2694
2702template <int rank_, int dim, typename Number, typename OtherNumber>
2703constexpr DEAL_II_ALWAYS_INLINE
2706 const Tensor<rank_, dim, OtherNumber> & right)
2707{
2708 return Tensor<rank_, dim, Number>(left) - right;
2709}
2710
2711
2719template <int rank_, int dim, typename Number, typename OtherNumber>
2720constexpr DEAL_II_ALWAYS_INLINE
2724{
2725 return left - Tensor<rank_, dim, OtherNumber>(right);
2726}
2727
2728
2729
2743template <int dim, typename Number>
2746{
2747 switch (dim)
2748 {
2749 case 1:
2750 return t.data[0];
2751 case 2:
2752 return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2753 case 3:
2754 {
2755 // in analogy to general tensors, but
2756 // there's something to be simplified for
2757 // the present case
2758 const Number tmp = t.data[3] * t.data[4] * t.data[5];
2759 return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2760 t.data[0] * t.data[5] * t.data[5] -
2761 t.data[1] * t.data[4] * t.data[4] -
2762 t.data[2] * t.data[3] * t.data[3]);
2763 }
2764 default:
2765 Assert(false, ExcNotImplemented());
2767 }
2768}
2769
2770
2771
2783template <int dim, typename Number>
2786{
2787 return determinant(t);
2788}
2789
2790
2791
2801template <int dim, typename Number>
2802constexpr inline DEAL_II_ALWAYS_INLINE Number
2804{
2805 Number t = d.data[0];
2806 for (unsigned int i = 1; i < dim; ++i)
2807 t += d.data[i];
2808 return t;
2809}
2810
2811
2823template <int dim, typename Number>
2824constexpr Number
2826{
2827 return trace(t);
2828}
2829
2830
2842template <typename Number>
2843constexpr DEAL_II_ALWAYS_INLINE Number
2845{
2847}
2848
2849
2850
2869template <typename Number>
2870constexpr DEAL_II_ALWAYS_INLINE Number
2872{
2873 return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2874}
2875
2876
2877
2886template <typename Number>
2887constexpr DEAL_II_ALWAYS_INLINE Number
2889{
2890 return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2891 t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2892}
2893
2894
2895
2903template <typename Number>
2904std::array<Number, 1>
2906
2907
2908
2931template <typename Number>
2932std::array<Number, 2>
2934
2935
2936
2959template <typename Number>
2960std::array<Number, 3>
2962
2963
2964
2965namespace internal
2966{
2967 namespace SymmetricTensorImplementation
2968 {
3006 template <int dim, typename Number>
3007 void
3008 tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
3010 std::array<Number, dim> & d,
3011 std::array<Number, dim - 1> & e);
3012
3013
3014
3054 template <int dim, typename Number>
3055 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3056 ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
3057
3058
3059
3099 template <int dim, typename Number>
3100 std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3102
3103
3104
3118 template <typename Number>
3119 std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3120 hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3121
3122
3123
3156 template <typename Number>
3157 std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3158 hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3159
3164 template <int dim, typename Number>
3166 {
3167 using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3168 bool
3169 operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3170 {
3171 return lhs.first > rhs.first;
3172 }
3173 };
3174
3175 } // namespace SymmetricTensorImplementation
3176
3177} // namespace internal
3178
3179
3180
3181// The line below is to ensure that doxygen puts the full description
3182// of this global enumeration into the documentation
3183// See https://stackoverflow.com/a/1717984
3213{
3223 hybrid,
3233 ql_implicit_shifts,
3241 jacobi
3242};
3243
3244
3245
3274template <int dim, typename Number>
3275std::array<std::pair<Number, Tensor<1, dim, Number>>,
3276 std::integral_constant<int, dim>::value>
3280
3281
3282
3291template <int rank_, int dim, typename Number>
3294{
3295 return t;
3296}
3297
3298
3299
3310template <int dim, typename Number>
3313{
3315
3316 // subtract scaled trace from the diagonal
3317 const Number tr = trace(t) / dim;
3318 for (unsigned int i = 0; i < dim; ++i)
3319 tmp.data[i] -= tr;
3320
3321 return tmp;
3322}
3323
3324
3325
3326template <int dim, typename Number>
3329{
3330 // create a default constructed matrix filled with
3331 // zeros, then set the diagonal elements to one
3333 switch (dim)
3334 {
3335 case 1:
3337 break;
3338 case 2:
3339 tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3340 break;
3341 case 3:
3342 tmp.data[0] = tmp.data[1] = tmp.data[2] =
3344 break;
3345 default:
3346 for (unsigned int d = 0; d < dim; ++d)
3348 }
3349 return tmp;
3350}
3351
3352
3353
3354template <int dim, typename Number>
3357{
3359
3360 // fill the elements treating the diagonal
3361 for (unsigned int i = 0; i < dim; ++i)
3362 for (unsigned int j = 0; j < dim; ++j)
3363 tmp.data[i][j] =
3364 internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3365
3366 // then fill the ones that copy over the
3367 // non-diagonal elements. note that during
3368 // the double-contraction, we handle the
3369 // off-diagonal elements twice, so simply
3370 // copying requires a weight of 1/2
3371 for (unsigned int i = dim;
3372 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3373 n_rank2_components;
3374 ++i)
3376
3377 return tmp;
3378}
3379
3380
3381
3382template <int dim, typename Number>
3385{
3387
3388 // fill the elements treating the diagonal
3389 for (unsigned int i = 0; i < dim; ++i)
3391
3392 // then fill the ones that copy over the
3393 // non-diagonal elements. note that during
3394 // the double-contraction, we handle the
3395 // off-diagonal elements twice, so simply
3396 // copying requires a weight of 1/2
3397 for (unsigned int i = dim;
3398 i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3399 n_rank2_components;
3400 ++i)
3402
3403 return tmp;
3404}
3405
3406
3407
3417template <int dim, typename Number>
3420{
3422 value(t);
3423}
3424
3425
3426
3437template <int dim, typename Number>
3440{
3442 value(t);
3443}
3444
3445
3446
3468template <int dim, typename Number>
3469constexpr inline SymmetricTensor<4, dim, Number>
3472{
3474
3475 // fill only the elements really needed
3476 for (unsigned int i = 0; i < dim; ++i)
3477 for (unsigned int j = i; j < dim; ++j)
3478 for (unsigned int k = 0; k < dim; ++k)
3479 for (unsigned int l = k; l < dim; ++l)
3480 tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3481
3482 return tmp;
3483}
3484
3485
3486
3494template <int dim, typename Number>
3497{
3499 for (unsigned int d = 0; d < dim; ++d)
3500 result[d][d] = t[d][d];
3501
3502 const Number half = internal::NumberType<Number>::value(0.5);
3503 for (unsigned int d = 0; d < dim; ++d)
3504 for (unsigned int e = d + 1; e < dim; ++e)
3505 result[d][e] = (t[d][e] + t[e][d]) * half;
3506 return result;
3507}
3508
3509
3510
3518template <int rank_, int dim, typename Number>
3520operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3521{
3523 tt *= factor;
3524 return tt;
3525}
3526
3527
3528
3536template <int rank_, int dim, typename Number>
3538operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3539{
3540 // simply forward to the other operator
3541 return t * factor;
3542}
3543
3544
3545
3571template <int rank_, int dim, typename Number, typename OtherNumber>
3572constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3573 rank_,
3574 dim,
3575 typename ProductType<Number,
3576 typename EnableIfScalar<OtherNumber>::type>::type>
3578 const OtherNumber & factor)
3579{
3580 // form the product. we have to convert the two factors into the final
3581 // type via explicit casts because, for awkward reasons, the C++
3582 // standard committee saw it fit to not define an
3583 // operator*(float,std::complex<double>)
3584 // (as well as with switched arguments and double<->float).
3585 using product_type = typename ProductType<Number, OtherNumber>::type;
3588 return tt;
3589}
3590
3591
3592
3601template <int rank_, int dim, typename Number, typename OtherNumber>
3602constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3603 rank_,
3604 dim,
3605 typename ProductType<OtherNumber,
3606 typename EnableIfScalar<Number>::type>::type>
3607operator*(const Number & factor,
3609{
3610 // simply forward to the other operator with switched arguments
3611 return (t * factor);
3612}
3613
3614
3615
3621template <int rank_, int dim, typename Number, typename OtherNumber>
3622constexpr inline SymmetricTensor<
3623 rank_,
3624 dim,
3625 typename ProductType<Number,
3626 typename EnableIfScalar<OtherNumber>::type>::type>
3628 const OtherNumber & factor)
3629{
3630 using product_type = typename ProductType<Number, OtherNumber>::type;
3633 return tt;
3634}
3635
3636
3637
3644template <int rank_, int dim>
3646operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3647{
3649 tt *= factor;
3650 return tt;
3651}
3652
3653
3654
3661template <int rank_, int dim>
3663operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3664{
3666 tt *= factor;
3667 return tt;
3668}
3669
3670
3671
3677template <int rank_, int dim>
3678constexpr inline SymmetricTensor<rank_, dim>
3679operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3680{
3682 tt /= factor;
3683 return tt;
3684}
3685
3695template <int dim, typename Number, typename OtherNumber>
3699{
3700 return (t1 * t2);
3701}
3702
3703
3717template <int dim, typename Number, typename OtherNumber>
3718constexpr inline DEAL_II_ALWAYS_INLINE
3721 const Tensor<2, dim, OtherNumber> & t2)
3722{
3724 typename ProductType<Number, OtherNumber>::type>::value(0.0);
3725 for (unsigned int i = 0; i < dim; ++i)
3726 for (unsigned int j = 0; j < dim; ++j)
3727 s += t1[i][j] * t2[i][j];
3728 return s;
3729}
3730
3731
3745template <int dim, typename Number, typename OtherNumber>
3749{
3750 return scalar_product(t2, t1);
3751}
3752
3753
3768template <typename Number, typename OtherNumber>
3769constexpr inline DEAL_II_ALWAYS_INLINE void
3774{
3775 tmp[0][0] = t[0][0][0][0] * s[0][0];
3776}
3777
3778
3779
3794template <typename Number, typename OtherNumber>
3795constexpr inline void
3800{
3801 tmp[0][0] = t[0][0][0][0] * s[0][0];
3802}
3803
3804
3805
3820template <typename Number, typename OtherNumber>
3821constexpr inline void
3826{
3827 const unsigned int dim = 2;
3828
3829 for (unsigned int i = 0; i < dim; ++i)
3830 for (unsigned int j = i; j < dim; ++j)
3831 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3832 2 * t[i][j][0][1] * s[0][1];
3833}
3834
3835
3836
3851template <typename Number, typename OtherNumber>
3852constexpr inline void
3857{
3858 const unsigned int dim = 2;
3859
3860 for (unsigned int i = 0; i < dim; ++i)
3861 for (unsigned int j = i; j < dim; ++j)
3862 tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3863 2 * s[0][1] * t[0][1][i][j];
3864}
3865
3866
3867
3882template <typename Number, typename OtherNumber>
3883constexpr inline void
3888{
3889 const unsigned int dim = 3;
3890
3891 for (unsigned int i = 0; i < dim; ++i)
3892 for (unsigned int j = i; j < dim; ++j)
3893 tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3894 t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3895 2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3896}
3897
3898
3899
3914template <typename Number, typename OtherNumber>
3915constexpr inline void
3920{
3921 const unsigned int dim = 3;
3922
3923 for (unsigned int i = 0; i < dim; ++i)
3924 for (unsigned int j = i; j < dim; ++j)
3925 tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3926 s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3927 2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3928}
3929
3930
3931
3938template <int dim, typename Number, typename OtherNumber>
3941 const Tensor<1, dim, OtherNumber> & src2)
3942{
3944 for (unsigned int i = 0; i < dim; ++i)
3945 for (unsigned int j = 0; j < dim; ++j)
3946 dest[i] += src1[i][j] * src2[j];
3947 return dest;
3948}
3949
3950
3957template <int dim, typename Number, typename OtherNumber>
3961{
3962 // this is easy for symmetric tensors:
3963 return src2 * src1;
3964}
3965
3966
3967
3987template <int rank_1,
3988 int rank_2,
3989 int dim,
3990 typename Number,
3991 typename OtherNumber>
3992constexpr DEAL_II_ALWAYS_INLINE
3993 typename Tensor<rank_1 + rank_2 - 2,
3994 dim,
3995 typename ProductType<Number, OtherNumber>::type>::tensor_type
3998{
3999 return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4000}
4001
4002
4003
4023template <int rank_1,
4024 int rank_2,
4025 int dim,
4026 typename Number,
4027 typename OtherNumber>
4028constexpr DEAL_II_ALWAYS_INLINE
4029 typename Tensor<rank_1 + rank_2 - 2,
4030 dim,
4031 typename ProductType<Number, OtherNumber>::type>::tensor_type
4034{
4035 return Tensor<rank_1, dim, Number>(src1) * src2;
4036}
4037
4038
4039
4049template <int dim, typename Number>
4050inline std::ostream &
4051operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4052{
4053 // make our lives a bit simpler by outputting
4054 // the tensor through the operator for the
4055 // general Tensor class
4057
4058 for (unsigned int i = 0; i < dim; ++i)
4059 for (unsigned int j = 0; j < dim; ++j)
4060 tt[i][j] = t[i][j];
4061
4062 return out << tt;
4063}
4064
4065
4066
4076template <int dim, typename Number>
4077inline std::ostream &
4078operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4079{
4080 // make our lives a bit simpler by outputting
4081 // the tensor through the operator for the
4082 // general Tensor class
4084
4085 for (unsigned int i = 0; i < dim; ++i)
4086 for (unsigned int j = 0; j < dim; ++j)
4087 for (unsigned int k = 0; k < dim; ++k)
4088 for (unsigned int l = 0; l < dim; ++l)
4089 tt[i][j][k][l] = t[i][j][k][l];
4090
4091 return out << tt;
4092}
4093
4094
4096
4097#endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
constexpr bool operator==(const SymmetricTensor &) const
constexpr Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
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)
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)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
constexpr SymmetricTensor< rank_, dim > operator*(const SymmetricTensor< rank_, dim > &t, const double factor)
constexpr Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const SymmetricTensor< 2, dim, Number > &src1, const Tensor< 1, dim, OtherNumber > &src2)
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr SymmetricTensor(const Number(&array)[n_independent_components])
constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
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 constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
constexpr Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const Tensor< 1, dim, Number > &src1, const SymmetricTensor< 2, dim, OtherNumber > &src2)
std::array< Number, 2 > eigenvalues(const SymmetricTensor< 2, 2, Number > &T)
friend class SymmetricTensor
void serialize(Archive &ar, const unsigned int version)
constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
friend constexpr Number2 determinant(const SymmetricTensor< 2, dim2, Number2 > &t)
const Number * begin_raw() const
const Number * end_raw() const
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, false, rank_ - 1, Number > operator[](const unsigned int row)
typename base_tensor_descriptor::base_tensor_type base_tensor_type
constexpr const Number & operator()(const TableIndices< rank_ > &indices) const
friend constexpr SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
constexpr bool operator!=(const SymmetricTensor &) const
constexpr Number & operator[](const TableIndices< rank_ > &indices)
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)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
static constexpr std::size_t memory_consumption()
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, OtherNumber > &t2)
friend constexpr SymmetricTensor< 2, dim2, Number2 > deviator(const SymmetricTensor< 2, dim2, Number2 > &t)
friend constexpr Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
Number * begin_raw()
base_tensor_type data
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
constexpr const Number & operator[](const TableIndices< rank_ > &indices) const
constexpr SymmetricTensor & operator=(const Number &d)
SymmetricTensor(const Tensor< 2, dim, OtherNumber > &t)
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr Number & access_raw_entry(const unsigned int unrolled_index)
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
constexpr numbers::NumberTraits< Number >::real_type norm() const
constexpr SymmetricTensor operator-() const
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
constexpr SymmetricTensor()=default
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)
constexpr Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor(const SymmetricTensor< rank_, dim, OtherNumber > &initializer)
constexpr SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor< rank_, dim, Number > operator*(const Number &factor, const SymmetricTensor< rank_, dim, Number > &t)
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)
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)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
friend constexpr SymmetricTensor< 4, dim2, Number2 > identity_tensor()
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)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
constexpr SymmetricTensor< rank_, dim > operator*(const double factor, const SymmetricTensor< rank_, dim > &t)
constexpr SymmetricTensor & operator/=(const OtherNumber &factor)
constexpr void clear()
constexpr Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
constexpr SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
friend constexpr SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 4, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 4, dim, OtherNumber > &s) const
constexpr SymmetricTensor< rank_, dim > operator/(const SymmetricTensor< rank_, dim > &t, const double factor)
std::array< Number, 3 > eigenvalues(const SymmetricTensor< 2, 3, Number > &T)
Number * end_raw()
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
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)
constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
constexpr Number & operator()(const TableIndices< rank_ > &indices)
Definition: tensor.h:503
typename AccessorTypes< rank, dim, constness, Number >::reference reference
constexpr Accessor(tensor_type &tensor, const TableIndices< rank > &previous_indices)
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
constexpr Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i)
constexpr Accessor(tensor_type &tensor, const TableIndices< rank > &previous_indices)
constexpr Accessor(const Accessor &)=default
typename AccessorTypes< rank, dim, constness, Number >::reference reference
constexpr Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i) const
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_CONSTEXPR
Definition: config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression fabs(const Expression &x)
static const char A
static const char T
static const char N
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
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)
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
constexpr TableIndices< 4 > merge(const TableIndices< 4 > &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:943
static const unsigned int invalid_unsigned_int
Definition: types.h:201
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 > &)
constexpr SymmetricTensor< rank_, dim, typename ProductType< OtherNumber, typename EnableIfScalar< Number >::type >::type > operator*(const Number &factor, const SymmetricTensor< rank_, dim, OtherNumber > &t)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static constexpr const T & value(const T &t)
Definition: numbers.h:705
typename ProductType< Number, OtherNumber >::type value_type
::SymmetricTensor< rank1+rank2 - 4, dim, value_type > 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:611
constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
SymmetricTensorEigenvectorMethod
constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()