Reference documentation for deal.II version GIT 6d02bd1105 2022-05-22 03:35:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symmetric_tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2021 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 
22 #include <deal.II/base/numbers.h>
25 #include <deal.II/base/tensor.h>
26 
27 #include <array>
28 
30 
31 // Forward declaration
32 #ifndef DOXYGEN
33 template <int rank, int dim, typename Number = double>
34 class SymmetricTensor;
35 #endif
36 
43 template <int dim, typename Number = double>
46 
75 template <int dim, typename Number = double>
78 
116 template <int dim, typename Number = double>
118  identity_tensor();
119 
120 template <int dim, typename Number>
123 
124 template <int dim, typename Number>
127 
128 template <int dim2, typename Number>
129 constexpr inline DEAL_II_ALWAYS_INLINE Number
131 
132 template <int dim, typename Number>
135 
136 template <int dim, typename Number>
139 
140 
141 
142 namespace 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>
316  struct StorageType;
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>
439  class Accessor
440  {
441  public:
445  using reference =
447  using tensor_type =
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 =
526  using tensor_type =
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 
667 template <int rank_, int dim, typename Number>
669 {
670 public:
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(
743  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
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
823  operator==(const SymmetricTensor &) const;
824 
828  constexpr bool
829  operator!=(const SymmetricTensor &) const;
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 
1003  constexpr typename numbers::NumberTraits<Number>::real_type
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 
1055 private:
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
1117 template <int rank, int dim, typename Number>
1118 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1119 
1120 template <int rank_, int dim, typename Number>
1121 constexpr unsigned int
1122  SymmetricTensor<rank_, dim, Number>::n_independent_components;
1123 
1124 namespace internal
1125 {
1126  namespace SymmetricTensorAccessors
1127  {
1128  template <int rank_, int dim, bool constness, int P, typename Number>
1129  constexpr DEAL_II_ALWAYS_INLINE
1131  tensor_type & tensor,
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 
1197 template <int rank_, int dim, typename Number>
1198 template <typename OtherNumber>
1199 inline DEAL_II_ALWAYS_INLINE
1201  const Tensor<2, dim, OtherNumber> &t)
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 
1219 template <int rank_, int dim, typename Number>
1220 template <typename OtherNumber>
1221 constexpr DEAL_II_ALWAYS_INLINE
1223  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1224  : data(initializer.data)
1225 {}
1226 
1227 
1228 
1229 template <int rank_, int dim, typename Number>
1230 constexpr 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),
1238  ExcInternalError());
1239 }
1240 
1241 
1242 
1243 template <int rank_, int dim, typename Number>
1244 template <typename OtherNumber>
1248 {
1249  data = t.data;
1250  return *this;
1251 }
1252 
1253 
1254 
1255 template <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 
1269 namespace 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 
1617 template <int rank_, int dim, typename Number>
1618 constexpr DEAL_II_ALWAYS_INLINE
1620  const
1621 {
1622  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1623 }
1624 
1625 
1626 
1627 template <int rank_, int dim, typename Number>
1628 constexpr bool
1630  const SymmetricTensor<rank_, dim, Number> &t) const
1631 {
1632  return data == t.data;
1633 }
1634 
1635 
1636 
1637 template <int rank_, int dim, typename Number>
1638 constexpr bool
1640  const SymmetricTensor<rank_, dim, Number> &t) const
1641 {
1642  return data != t.data;
1643 }
1644 
1645 
1646 
1647 template <int rank_, int dim, typename Number>
1648 template <typename OtherNumber>
1652 {
1653  data += t.data;
1654  return *this;
1655 }
1656 
1657 
1658 
1659 template <int rank_, int dim, typename Number>
1660 template <typename OtherNumber>
1664 {
1665  data -= t.data;
1666  return *this;
1667 }
1668 
1669 
1670 
1671 template <int rank_, int dim, typename Number>
1672 template <typename OtherNumber>
1675 {
1676  data *= d;
1677  return *this;
1678 }
1679 
1680 
1681 
1682 template <int rank_, int dim, typename Number>
1683 template <typename OtherNumber>
1686 {
1687  data /= d;
1688  return *this;
1689 }
1690 
1691 
1692 
1693 template <int rank_, int dim, typename Number>
1696 {
1697  SymmetricTensor tmp = *this;
1698  tmp.data = -tmp.data;
1699  return tmp;
1700 }
1701 
1702 
1703 
1704 template <int rank_, int dim, typename Number>
1705 constexpr inline DEAL_II_ALWAYS_INLINE void
1707 {
1708  data.clear();
1709 }
1710 
1711 
1712 
1713 template <int rank_, int dim, typename Number>
1714 constexpr std::size_t
1716 {
1717  // all memory consists of statically allocated memory of the current
1718  // object, no pointers
1719  return sizeof(SymmetricTensor<rank_, dim, Number>);
1720 }
1721 
1722 
1723 
1724 namespace 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 
1866 template <int rank_, int dim, typename Number>
1867 template <typename OtherNumber>
1872  const SymmetricTensor<2, dim, OtherNumber> &s) const
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 
1884 template <int rank_, int dim, typename Number>
1885 template <typename OtherNumber>
1889  const SymmetricTensor<4, dim, OtherNumber> &s) const
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
1909 namespace 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)),
1948  ExcInternalError());
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)),
1994  ExcInternalError());
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 
2115 template <int rank_, int dim, typename Number>
2116 constexpr 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 
2127 template <int rank_, int dim, typename Number>
2128 constexpr 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 
2139 namespace 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 
2166 template <int rank_, int dim, typename Number>
2167 constexpr 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 
2180 template <int rank_, int dim, typename Number>
2181 constexpr inline DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2182  Accessor<rank_, dim, false, rank_ - 1, Number>
2183  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
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 
2194 template <int rank_, int dim, typename Number>
2195 constexpr DEAL_II_ALWAYS_INLINE const Number &
2197  const TableIndices<rank_> &indices) const
2198 {
2199  return operator()(indices);
2200 }
2201 
2202 
2203 
2204 template <int rank_, int dim, typename Number>
2205 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2207  const TableIndices<rank_> &indices)
2208 {
2209  return operator()(indices);
2210 }
2211 
2212 
2213 
2214 template <int rank_, int dim, typename Number>
2215 inline Number *
2217 {
2218  return std::addressof(this->access_raw_entry(0));
2219 }
2220 
2221 
2222 
2223 template <int rank_, int dim, typename Number>
2224 inline const Number *
2226 {
2227  return std::addressof(this->access_raw_entry(0));
2228 }
2229 
2230 
2231 
2232 template <int rank_, int dim, typename Number>
2233 inline Number *
2235 {
2236  return begin_raw() + n_independent_components;
2237 }
2238 
2239 
2240 
2241 template <int rank_, int dim, typename Number>
2242 inline const Number *
2244 {
2245  return begin_raw() + n_independent_components;
2246 }
2247 
2248 
2249 
2250 namespace 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 
2277 template <int rank_, int dim, typename Number>
2278 constexpr 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 
2289 template <int rank_, int dim, typename Number>
2290 constexpr inline Number &
2292 {
2293  AssertIndexRange(index, n_independent_components);
2294  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2295  index)];
2296 }
2297 
2298 
2299 
2300 namespace 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 
2389 template <int rank_, int dim, typename Number>
2392 {
2393  return internal::compute_norm<dim, Number>(data);
2394 }
2395 
2396 
2397 
2398 namespace 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 
2483 template <int rank_, int dim, typename Number>
2484 constexpr unsigned int
2486  const TableIndices<rank_> &indices)
2487 {
2488  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2489  dim>(indices);
2490 }
2491 
2492 
2493 
2494 namespace 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(
2512  ExcIndexRange(
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 
2586 template <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 
2597 template <int rank_, int dim, typename Number>
2598 template <class Archive>
2599 inline void
2600 SymmetricTensor<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 
2623 template <int rank_, int dim, typename Number, typename OtherNumber>
2624 constexpr inline DEAL_II_ALWAYS_INLINE
2628 {
2630  tmp = left;
2631  tmp += right;
2632  return tmp;
2633 }
2634 
2635 
2648 template <int rank_, int dim, typename Number, typename OtherNumber>
2649 constexpr inline DEAL_II_ALWAYS_INLINE
2653 {
2655  tmp = left;
2656  tmp -= right;
2657  return tmp;
2658 }
2659 
2660 
2668 template <int rank_, int dim, typename Number, typename OtherNumber>
2669 constexpr DEAL_II_ALWAYS_INLINE
2672  const Tensor<rank_, dim, OtherNumber> & right)
2673 {
2674  return Tensor<rank_, dim, Number>(left) + right;
2675 }
2676 
2677 
2685 template <int rank_, int dim, typename Number, typename OtherNumber>
2686 constexpr DEAL_II_ALWAYS_INLINE
2690 {
2691  return left + Tensor<rank_, dim, OtherNumber>(right);
2692 }
2693 
2694 
2702 template <int rank_, int dim, typename Number, typename OtherNumber>
2703 constexpr DEAL_II_ALWAYS_INLINE
2706  const Tensor<rank_, dim, OtherNumber> & right)
2707 {
2708  return Tensor<rank_, dim, Number>(left) - right;
2709 }
2710 
2711 
2719 template <int rank_, int dim, typename Number, typename OtherNumber>
2720 constexpr DEAL_II_ALWAYS_INLINE
2724 {
2725  return left - Tensor<rank_, dim, OtherNumber>(right);
2726 }
2727 
2728 
2729 
2743 template <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 
2783 template <int dim, typename Number>
2786 {
2787  return determinant(t);
2788 }
2789 
2790 
2791 
2801 template <int dim, typename Number>
2802 constexpr 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 
2823 template <int dim, typename Number>
2824 constexpr Number
2826 {
2827  return trace(t);
2828 }
2829 
2830 
2842 template <typename Number>
2843 constexpr DEAL_II_ALWAYS_INLINE Number
2845 {
2847 }
2848 
2849 
2850 
2869 template <typename Number>
2870 constexpr DEAL_II_ALWAYS_INLINE Number
2872 {
2873  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2874 }
2875 
2876 
2877 
2886 template <typename Number>
2887 constexpr 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 
2903 template <typename Number>
2904 std::array<Number, 1>
2906 
2907 
2908 
2931 template <typename Number>
2932 std::array<Number, 2>
2934 
2935 
2936 
2959 template <typename Number>
2960 std::array<Number, 3>
2962 
2963 
2964 
2965 namespace internal
2966 {
2967  namespace SymmetricTensorImplementation
2968  {
3006  template <int dim, typename Number>
3007  void
3008  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
3009  ::Tensor<2, dim, Number> & Q,
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,
3241  jacobi
3242 };
3243 
3244 
3245 
3274 template <int dim, typename Number>
3275 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3276  std::integral_constant<int, dim>::value>
3278  const SymmetricTensorEigenvectorMethod method =
3280 
3281 
3282 
3291 template <int rank_, int dim, typename Number>
3294 {
3295  return t;
3296 }
3297 
3298 
3299 
3310 template <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 
3326 template <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 
3354 template <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)
3375  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3376 
3377  return tmp;
3378 }
3379 
3380 
3381 
3382 template <int dim, typename Number>
3385 {
3387 
3388  // fill the elements treating the diagonal
3389  for (unsigned int i = 0; i < dim; ++i)
3390  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
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)
3401  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3402 
3403  return tmp;
3404 }
3405 
3406 
3407 
3417 template <int dim, typename Number>
3420 {
3422  value(t);
3423 }
3424 
3425 
3426 
3437 template <int dim, typename Number>
3440 {
3442  value(t);
3443 }
3444 
3445 
3446 
3468 template <int dim, typename Number>
3469 constexpr 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 
3494 template <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 
3518 template <int rank_, int dim, typename Number>
3520 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3521 {
3523  tt *= factor;
3524  return tt;
3525 }
3526 
3527 
3528 
3536 template <int rank_, int dim, typename Number>
3538 operator*(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 
3571 template <int rank_, int dim, typename Number, typename OtherNumber>
3572 constexpr 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 
3601 template <int rank_, int dim, typename Number, typename OtherNumber>
3602 constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3603  rank_,
3604  dim,
3605  typename ProductType<OtherNumber,
3606  typename EnableIfScalar<Number>::type>::type>
3607 operator*(const Number & factor,
3609 {
3610  // simply forward to the other operator with switched arguments
3611  return (t * factor);
3612 }
3613 
3614 
3615 
3621 template <int rank_, int dim, typename Number, typename OtherNumber>
3622 constexpr 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 
3644 template <int rank_, int dim>
3646 operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3647 {
3649  tt *= factor;
3650  return tt;
3651 }
3652 
3653 
3654 
3661 template <int rank_, int dim>
3663 operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3664 {
3666  tt *= factor;
3667  return tt;
3668 }
3669 
3670 
3671 
3677 template <int rank_, int dim>
3678 constexpr inline SymmetricTensor<rank_, dim>
3679 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3680 {
3682  tt /= factor;
3683  return tt;
3684 }
3685 
3695 template <int dim, typename Number, typename OtherNumber>
3699 {
3700  return (t1 * t2);
3701 }
3702 
3703 
3717 template <int dim, typename Number, typename OtherNumber>
3718 constexpr 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 
3745 template <int dim, typename Number, typename OtherNumber>
3749 {
3750  return scalar_product(t2, t1);
3751 }
3752 
3753 
3768 template <typename Number, typename OtherNumber>
3769 constexpr inline DEAL_II_ALWAYS_INLINE void
3774 {
3775  tmp[0][0] = t[0][0][0][0] * s[0][0];
3776 }
3777 
3778 
3779 
3794 template <typename Number, typename OtherNumber>
3795 constexpr inline void
3800 {
3801  tmp[0][0] = t[0][0][0][0] * s[0][0];
3802 }
3803 
3804 
3805 
3820 template <typename Number, typename OtherNumber>
3821 constexpr 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 
3851 template <typename Number, typename OtherNumber>
3852 constexpr 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 
3882 template <typename Number, typename OtherNumber>
3883 constexpr 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 
3914 template <typename Number, typename OtherNumber>
3915 constexpr 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 
3938 template <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 
3957 template <int dim, typename Number, typename OtherNumber>
3961 {
3962  // this is easy for symmetric tensors:
3963  return src2 * src1;
3964 }
3965 
3966 
3967 
3987 template <int rank_1,
3988  int rank_2,
3989  int dim,
3990  typename Number,
3991  typename OtherNumber>
3992 constexpr 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 
4023 template <int rank_1,
4024  int rank_2,
4025  int dim,
4026  typename Number,
4027  typename OtherNumber>
4028 constexpr DEAL_II_ALWAYS_INLINE
4029  typename Tensor<rank_1 + rank_2 - 2,
4030  dim,
4031  typename ProductType<Number, OtherNumber>::type>::tensor_type
4033  const Tensor<rank_2, dim, OtherNumber> & src2)
4034 {
4035  return Tensor<rank_1, dim, Number>(src1) * src2;
4036 }
4037 
4038 
4039 
4049 template <int dim, typename Number>
4050 inline std::ostream &
4051 operator<<(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 
4076 template <int dim, typename Number>
4077 inline std::ostream &
4078 operator<<(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 const Number & operator()(const TableIndices< rank_ > &indices) const
constexpr bool operator==(const SymmetricTensor &) const
constexpr Number & operator()(const TableIndices< rank_ > &indices)
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 friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
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_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type operator*(const Tensor< rank_1, dim, Number > &src1, const SymmetricTensor< rank_2, dim, OtherNumber > &src2)
constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
Number * begin_raw()
constexpr const Number & operator[](const TableIndices< rank_ > &indices) const
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Number * end_raw()
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 internal::SymmetricTensorAccessors::Accessor< rank_, dim, false, rank_ - 1, Number > operator[](const unsigned int row)
constexpr SymmetricTensor & operator=(const Number &d)
constexpr Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const Tensor< 1, dim, Number > &src1, const SymmetricTensor< 2, dim, OtherNumber > &src2)
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
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)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr friend SymmetricTensor< 2, dim2, Number2 > deviator(const SymmetricTensor< 2, dim2, Number2 > &t)
constexpr numbers::NumberTraits< Number >::real_type norm() const
typename base_tensor_descriptor::base_tensor_type base_tensor_type
constexpr friend SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
constexpr bool operator!=(const SymmetricTensor &) const
const Number * end_raw() const
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)
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)
base_tensor_type data
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
const Number * begin_raw() const
SymmetricTensor(const Tensor< 2, dim, OtherNumber > &t)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr SymmetricTensor & operator/=(const OtherNumber &factor)
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
constexpr SymmetricTensor operator-() const
constexpr SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor()=default
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
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 & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor(const SymmetricTensor< rank_, dim, OtherNumber > &initializer)
constexpr SymmetricTensor< rank_, dim, Number > operator*(const Number &factor, const SymmetricTensor< rank_, dim, Number > &t)
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
constexpr friend Number2 determinant(const SymmetricTensor< 2, dim2, Number2 > &t)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr Number & access_raw_entry(const unsigned int unrolled_index)
constexpr friend 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 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 internal::SymmetricTensorAccessors::double_contraction_result< rank_, 4, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 4, dim, OtherNumber > &s) const
constexpr Tensor< rank_1+rank_2 - 2, dim, typename ProductType< Number, OtherNumber >::type >::tensor_type operator*(const SymmetricTensor< rank_1, dim, Number > &src1, const Tensor< rank_2, dim, OtherNumber > &src2)
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 & operator[](const TableIndices< rank_ > &indices)
constexpr Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
constexpr friend SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
constexpr SymmetricTensor< rank_, dim > operator/(const SymmetricTensor< rank_, dim > &t, const double factor)
std::array< Number, 3 > eigenvalues(const SymmetricTensor< 2, 3, Number > &T)
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)
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
constexpr Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i) const
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
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_CONSTEXPR
Definition: config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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 > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
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
::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()