Reference documentation for deal.II version GIT a189bc2bdf 2022-12-07 02:45:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
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 
47 template <int dim, typename Number = double>
50 
79 template <int dim, typename Number = double>
82 
120 template <int dim, typename Number = double>
123 
124 template <int dim, typename Number>
127 
128 template <int dim, typename Number>
131 
141 template <int dim2, typename Number>
142 constexpr inline DEAL_II_ALWAYS_INLINE Number
144 
155 template <int dim, typename Number>
158 
172 template <int dim, typename Number>
175 
176 
177 
178 namespace internal
179 {
180  // Workaround: The following 4 overloads are necessary to be able to
181  // compile the library with Apple Clang 8 and older. We should remove
182  // these overloads again when we bump the minimal required version to
183  // something later than clang-3.6 / Apple Clang 6.3.
184  template <int rank, int dim, typename T, typename U>
185  struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
186  {
187  using type =
188  SymmetricTensor<rank,
189  dim,
190  std::complex<typename ProductType<T, U>::type>>;
191  };
192 
193  template <int rank, int dim, typename T, typename U>
194  struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
195  std::complex<U>>
196  {
197  using type =
198  SymmetricTensor<rank,
199  dim,
200  std::complex<typename ProductType<T, U>::type>>;
201  };
202 
203  template <typename T, int rank, int dim, typename U>
204  struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
205  {
206  using type =
207  SymmetricTensor<rank,
208  dim,
209  std::complex<typename ProductType<T, U>::type>>;
210  };
211 
212  template <int rank, int dim, typename T, typename U>
213  struct ProductTypeImpl<std::complex<T>,
214  SymmetricTensor<rank, dim, std::complex<U>>>
215  {
216  using type =
217  SymmetricTensor<rank,
218  dim,
219  std::complex<typename ProductType<T, U>::type>>;
220  };
221  // end workaround
222 
227  namespace SymmetricTensorImplementation
228  {
233  template <int rank, int dim, typename Number>
234  struct Inverse;
235  } // namespace SymmetricTensorImplementation
236 
241  namespace SymmetricTensorAccessors
242  {
250  merge(const TableIndices<2> &previous_indices,
251  const unsigned int new_index,
252  const unsigned int position)
253  {
254  AssertIndexRange(position, 2);
255 
256  if (position == 0)
257  return {new_index, numbers::invalid_unsigned_int};
258  else
259  return {previous_indices[0], new_index};
260  }
261 
262 
263 
271  merge(const TableIndices<4> &previous_indices,
272  const unsigned int new_index,
273  const unsigned int position)
274  {
275  AssertIndexRange(position, 4);
276 
277  switch (position)
278  {
279  case 0:
280  return {new_index,
284  case 1:
285  return {previous_indices[0],
286  new_index,
289  case 2:
290  return {previous_indices[0],
291  previous_indices[1],
292  new_index,
294  case 3:
295  return {previous_indices[0],
296  previous_indices[1],
297  previous_indices[2],
298  new_index};
299  default:
300  Assert(false, ExcInternalError());
301  return {};
302  }
303  }
304 
305 
312  template <int rank1,
313  int rank2,
314  int dim,
315  typename Number,
316  typename OtherNumber = Number>
318  {
320  using type =
321  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
322  };
323 
324 
331  template <int dim, typename Number, typename OtherNumber>
332  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
333  {
335  };
336 
337 
338 
351  template <int rank, int dim, typename Number>
352  struct StorageType;
353 
357  template <int dim, typename Number>
358  struct StorageType<2, dim, Number>
359  {
364  static const unsigned int n_independent_components =
365  (dim * dim + dim) / 2;
366 
371  };
372 
373 
374 
378  template <int dim, typename Number>
379  struct StorageType<4, dim, Number>
380  {
386  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
387 
391  static const unsigned int n_independent_components =
392  (n_rank2_components *
394 
402  };
403 
404 
405 
410  template <int rank, int dim, bool constness, typename Number>
412 
419  template <int rank, int dim, typename Number>
420  struct AccessorTypes<rank, dim, true, Number>
421  {
422  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
423 
424  using reference = Number;
425  };
426 
433  template <int rank, int dim, typename Number>
434  struct AccessorTypes<rank, dim, false, Number>
435  {
437 
438  using reference = Number &;
439  };
440 
441 
474  template <int rank, int dim, bool constness, int P, typename Number>
475  class Accessor
476  {
477  public:
481  using reference =
483  using tensor_type =
485 
486  private:
507 
511  constexpr DEAL_II_ALWAYS_INLINE
512  Accessor(const Accessor &) = default;
513 
514  public:
518  constexpr Accessor<rank, dim, constness, P - 1, Number>
519  operator[](const unsigned int i);
520 
524  constexpr Accessor<rank, dim, constness, P - 1, Number>
525  operator[](const unsigned int i) const;
526 
527  private:
533 
534  // Declare some other classes as friends. Make sure to work around bugs
535  // in some compilers:
536  template <int, int, typename>
537  friend class ::SymmetricTensor;
538  template <int, int, bool, int, typename>
539  friend class Accessor;
540  friend class ::SymmetricTensor<rank, dim, Number>;
541  friend class Accessor<rank, dim, constness, P + 1, Number>;
542  };
543 
544 
545 
553  template <int rank, int dim, bool constness, typename Number>
554  class Accessor<rank, dim, constness, 1, Number>
555  {
556  public:
560  using reference =
562  using tensor_type =
564 
565  private:
589 
593  constexpr DEAL_II_ALWAYS_INLINE
594  Accessor(const Accessor &) = default;
595 
596  public:
600  constexpr reference
601  operator[](const unsigned int);
602 
606  constexpr reference
607  operator[](const unsigned int) const;
608 
609  private:
615 
616  // Declare some other classes as friends. Make sure to work around bugs
617  // in some compilers:
618  template <int, int, typename>
619  friend class ::SymmetricTensor;
620  template <int, int, bool, int, typename>
622  friend class ::SymmetricTensor<rank, dim, Number>;
623  friend class SymmetricTensorAccessors::
624  Accessor<rank, dim, constness, 2, Number>;
625  };
626  } // namespace SymmetricTensorAccessors
627 } // namespace internal
628 
629 
630 
703 template <int rank_, int dim, typename Number>
705 {
706 public:
707  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
708 
717  static constexpr unsigned int dimension = dim;
718 
722  static const unsigned int rank = rank_;
723 
729  static constexpr unsigned int n_independent_components =
731  n_independent_components;
732 
736  constexpr DEAL_II_ALWAYS_INLINE
737  SymmetricTensor() = default;
738 
752  template <typename OtherNumber>
754 
770  constexpr SymmetricTensor(const Number (&array)[n_independent_components]);
771 
777  template <typename OtherNumber>
778  constexpr explicit SymmetricTensor(
779  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
780 
790  Number *
792 
802  const Number *
803  begin_raw() const;
804 
814  Number *
816 
827  const Number *
828  end_raw() const;
829 
836  template <typename OtherNumber>
837  constexpr SymmetricTensor &
839 
846  constexpr SymmetricTensor &
847  operator=(const Number &d);
848 
853  constexpr operator Tensor<rank_, dim, Number>() const;
854 
858  constexpr bool
859  operator==(const SymmetricTensor &) const;
860 
864  constexpr bool
865  operator!=(const SymmetricTensor &) const;
866 
870  template <typename OtherNumber>
871  constexpr SymmetricTensor &
873 
877  template <typename OtherNumber>
878  constexpr SymmetricTensor &
880 
885  template <typename OtherNumber>
886  constexpr SymmetricTensor &
887  operator*=(const OtherNumber &factor);
888 
892  template <typename OtherNumber>
893  constexpr SymmetricTensor &
894  operator/=(const OtherNumber &factor);
895 
899  constexpr SymmetricTensor
900  operator-() const;
901 
954  template <typename OtherNumber>
958 
963  template <typename OtherNumber>
967 
971  constexpr Number &
973 
977  constexpr const Number &
978  operator()(const TableIndices<rank_> &indices) const;
979 
984  constexpr internal::SymmetricTensorAccessors::
985  Accessor<rank_, dim, true, rank_ - 1, Number>
986  operator[](const unsigned int row) const;
987 
992  constexpr internal::SymmetricTensorAccessors::
993  Accessor<rank_, dim, false, rank_ - 1, Number>
994  operator[](const unsigned int row);
995 
1001  constexpr const Number &
1002  operator[](const TableIndices<rank_> &indices) const;
1003 
1009  constexpr Number &
1011 
1018  constexpr const Number &
1019  access_raw_entry(const unsigned int unrolled_index) const;
1020 
1027  constexpr Number &
1028  access_raw_entry(const unsigned int unrolled_index);
1029 
1039  constexpr typename numbers::NumberTraits<Number>::real_type
1040  norm() const;
1041 
1049  static constexpr unsigned int
1051 
1057  static constexpr TableIndices<rank_>
1058  unrolled_to_component_indices(const unsigned int i);
1059 
1072  constexpr void
1074 
1079  static constexpr std::size_t
1081 
1087  template <class Archive>
1088  void
1089  serialize(Archive &ar, const unsigned int version);
1090 
1091 private:
1097 
1101  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
1102 
1107 
1108 #ifndef DOXYGEN
1109 
1110  // Make all other symmetric tensors friends.
1111  template <int, int, typename>
1112  friend class SymmetricTensor;
1113 
1114  // Make a few more functions friends.
1115  template <int dim2, typename Number2>
1116  friend constexpr Number2
1118 
1119  template <int dim2, typename Number2>
1120  friend DEAL_II_CONSTEXPR Number2
1122 
1123  template <int dim2, typename Number2>
1124  friend constexpr SymmetricTensor<2, dim2, Number2>
1126 
1127  template <int dim2, typename Number2>
1130 
1131  template <int dim2, typename Number2>
1133  deviator_tensor();
1134 
1135  template <int dim2, typename Number2>
1137  identity_tensor();
1138 
1139 
1140  // Make a few helper classes friends as well.
1142  Inverse<2, dim, Number>;
1143 
1145  Inverse<4, dim, Number>;
1146 #endif
1147 };
1148 
1149 
1150 
1151 // ------------------------- inline functions ------------------------
1152 
1153 #ifndef DOXYGEN
1154 
1155 // provide declarations for static members
1156 template <int rank, int dim, typename Number>
1157 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1158 
1159 template <int rank_, int dim, typename Number>
1160 constexpr unsigned int
1161  SymmetricTensor<rank_, dim, Number>::n_independent_components;
1162 
1163 namespace internal
1164 {
1165  namespace SymmetricTensorAccessors
1166  {
1167  template <int rank_, int dim, bool constness, int P, typename Number>
1168  constexpr DEAL_II_ALWAYS_INLINE
1170  tensor_type & tensor,
1172  : tensor(tensor)
1174  {}
1175 
1176 
1177 
1178  template <int rank_, int dim, bool constness, int P, typename Number>
1179  constexpr inline DEAL_II_ALWAYS_INLINE
1180  Accessor<rank_, dim, constness, P - 1, Number>
1181  Accessor<rank_, dim, constness, P, Number>::operator[](
1182  const unsigned int i)
1183  {
1184  return Accessor<rank_, dim, constness, P - 1, Number>(
1185  tensor, merge(previous_indices, i, rank_ - P));
1186  }
1187 
1188 
1189 
1190  template <int rank_, int dim, bool constness, int P, typename Number>
1191  constexpr DEAL_II_ALWAYS_INLINE
1192  Accessor<rank_, dim, constness, P - 1, Number>
1193  Accessor<rank_, dim, constness, P, Number>::operator[](
1194  const unsigned int i) const
1195  {
1196  return Accessor<rank_, dim, constness, P - 1, Number>(
1197  tensor, merge(previous_indices, i, rank_ - P));
1198  }
1199 
1200 
1201 
1202  template <int rank_, int dim, bool constness, typename Number>
1203  constexpr DEAL_II_ALWAYS_INLINE
1204  Accessor<rank_, dim, constness, 1, Number>::Accessor(
1205  tensor_type & tensor,
1206  const TableIndices<rank_> &previous_indices)
1207  : tensor(tensor)
1208  , previous_indices(previous_indices)
1209  {}
1210 
1211 
1212 
1213  template <int rank_, int dim, bool constness, typename Number>
1214  constexpr inline DEAL_II_ALWAYS_INLINE
1215  typename Accessor<rank_, dim, constness, 1, Number>::reference
1216  Accessor<rank_, dim, constness, 1, Number>::operator[](
1217  const unsigned int i)
1218  {
1219  return tensor(merge(previous_indices, i, rank_ - 1));
1220  }
1221 
1222 
1223  template <int rank_, int dim, bool constness, typename Number>
1224  constexpr DEAL_II_ALWAYS_INLINE
1225  typename Accessor<rank_, dim, constness, 1, Number>::reference
1226  Accessor<rank_, dim, constness, 1, Number>::operator[](
1227  const unsigned int i) const
1228  {
1229  return tensor(merge(previous_indices, i, rank_ - 1));
1230  }
1231  } // namespace SymmetricTensorAccessors
1232 } // namespace internal
1233 
1234 
1235 
1236 template <int rank_, int dim, typename Number>
1237 template <typename OtherNumber>
1238 inline DEAL_II_ALWAYS_INLINE
1240  const Tensor<2, dim, OtherNumber> &t)
1241 {
1242  static_assert(rank == 2, "This function is only implemented for rank==2");
1243  for (unsigned int d = 0; d < dim; ++d)
1244  for (unsigned int e = 0; e < d; ++e)
1245  Assert(t[d][e] == t[e][d],
1246  ExcMessage("The incoming Tensor must be exactly symmetric."));
1247 
1248  for (unsigned int d = 0; d < dim; ++d)
1249  data[d] = t[d][d];
1250 
1251  for (unsigned int d = 0, c = 0; d < dim; ++d)
1252  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1253  data[dim + c] = t[d][e];
1254 }
1255 
1256 
1257 
1258 template <int rank_, int dim, typename Number>
1259 template <typename OtherNumber>
1260 constexpr DEAL_II_ALWAYS_INLINE
1262  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1263  : data(initializer.data)
1264 {}
1265 
1266 
1267 
1268 template <int rank_, int dim, typename Number>
1269 constexpr inline DEAL_II_ALWAYS_INLINE
1271  const Number (&array)[n_independent_components])
1272  : data(
1273  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1274 {
1275  // ensure that the reinterpret_cast above actually works
1276  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1277  ExcInternalError());
1278 }
1279 
1280 
1281 
1282 template <int rank_, int dim, typename Number>
1283 template <typename OtherNumber>
1287 {
1288  data = t.data;
1289  return *this;
1290 }
1291 
1292 
1293 
1294 template <int rank_, int dim, typename Number>
1297 {
1299  ExcMessage("Only assignment with zero is allowed"));
1300  (void)d;
1301 
1303 
1304  return *this;
1305 }
1306 
1307 
1308 namespace internal
1309 {
1310  namespace SymmetricTensorImplementation
1311  {
1312  template <int dim, typename Number>
1313  constexpr inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1314  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1315  {
1317 
1318  // diagonal entries are stored first
1319  for (unsigned int d = 0; d < dim; ++d)
1320  t[d][d] = s.access_raw_entry(d);
1321 
1322  // off-diagonal entries come next, row by row
1323  for (unsigned int d = 0, c = 0; d < dim; ++d)
1324  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1325  {
1326  t[d][e] = s.access_raw_entry(dim + c);
1327  t[e][d] = s.access_raw_entry(dim + c);
1328  }
1329  return t;
1330  }
1331 
1332 
1333  template <int dim, typename Number>
1334  constexpr ::Tensor<4, dim, Number>
1335  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1336  {
1337  // utilize the symmetry properties of SymmetricTensor<4,dim>
1338  // discussed in the class documentation to avoid accessing all
1339  // independent elements of the input tensor more than once
1341 
1342  for (unsigned int i = 0; i < dim; ++i)
1343  for (unsigned int j = i; j < dim; ++j)
1344  for (unsigned int k = 0; k < dim; ++k)
1345  for (unsigned int l = k; l < dim; ++l)
1346  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1347  t[TableIndices<4>(j, i, k, l)] =
1348  t[TableIndices<4>(j, i, l, k)] =
1349  st[TableIndices<4>(i, j, k, l)];
1350 
1351  return t;
1352  }
1353 
1354 
1355  template <typename Number>
1356  struct Inverse<2, 1, Number>
1357  {
1358  constexpr static inline DEAL_II_ALWAYS_INLINE
1359  ::SymmetricTensor<2, 1, Number>
1360  value(const ::SymmetricTensor<2, 1, Number> &t)
1361  {
1363 
1364  tmp[0][0] = 1.0 / t[0][0];
1365 
1366  return tmp;
1367  }
1368  };
1369 
1370 
1371  template <typename Number>
1372  struct Inverse<2, 2, Number>
1373  {
1374  constexpr static inline DEAL_II_ALWAYS_INLINE
1375  ::SymmetricTensor<2, 2, Number>
1376  value(const ::SymmetricTensor<2, 2, Number> &t)
1377  {
1379 
1380  // Sympy result: ([
1381  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1382  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1383  const TableIndices<2> idx_00(0, 0);
1384  const TableIndices<2> idx_01(0, 1);
1385  const TableIndices<2> idx_11(1, 1);
1386  const Number inv_det_t =
1387  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1388  tmp[idx_00] = t[idx_11];
1389  tmp[idx_01] = -t[idx_01];
1390  tmp[idx_11] = t[idx_00];
1391  tmp *= inv_det_t;
1392 
1393  return tmp;
1394  }
1395  };
1396 
1397 
1398  template <typename Number>
1399  struct Inverse<2, 3, Number>
1400  {
1401  constexpr static ::SymmetricTensor<2, 3, Number>
1402  value(const ::SymmetricTensor<2, 3, Number> &t)
1403  {
1405 
1406  // Sympy result: ([
1407  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1408  // 2*t01*t02*t12 - t02**2*t11),
1409  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1410  // 2*t01*t02*t12 - t02**2*t11),
1411  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1412  // 2*t01*t02*t12 - t02**2*t11)],
1413  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1414  // 2*t01*t02*t12 - t02**2*t11),
1415  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1416  // 2*t01*t02*t12 - t02**2*t11),
1417  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1418  // 2*t01*t02*t12 + t02**2*t11)],
1419  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1420  // 2*t01*t02*t12 - t02**2*t11),
1421  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1422  // 2*t01*t02*t12 + t02**2*t11),
1423  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1424  // 2*t01*t02*t12 + t02**2*t11)] ])
1425  //
1426  // =
1427  //
1428  // [ (t11*t22 - t12**2)/det_t,
1429  // (-t01*t22 + t02*t12)/det_t,
1430  // (t01*t12 - t02*t11)/det_t],
1431  // [ (-t01*t22 + t02*t12)/det_t,
1432  // (t00*t22 - t02**2)/det_t,
1433  // (-t00*t12 + t01*t02)/det_t],
1434  // [ (t01*t12 - t02*t11)/det_t,
1435  // (-t00*t12 + t01*t02)/det_t,
1436  // (t00*t11 - t01**2)/det_t] ])
1437  //
1438  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1439  // 2*t01*t02*t12 - t02**2*t11)
1440  const TableIndices<2> idx_00(0, 0);
1441  const TableIndices<2> idx_01(0, 1);
1442  const TableIndices<2> idx_02(0, 2);
1443  const TableIndices<2> idx_11(1, 1);
1444  const TableIndices<2> idx_12(1, 2);
1445  const TableIndices<2> idx_22(2, 2);
1446  const Number inv_det_t =
1447  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1448  t[idx_00] * t[idx_12] * t[idx_12] -
1449  t[idx_01] * t[idx_01] * t[idx_22] +
1450  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1451  t[idx_02] * t[idx_02] * t[idx_11]);
1452  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1453  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1454  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1455  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1456  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1457  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1458  tmp *= inv_det_t;
1459 
1460  return tmp;
1461  }
1462  };
1463 
1464 
1465  template <typename Number>
1466  struct Inverse<4, 1, Number>
1467  {
1468  constexpr static inline ::SymmetricTensor<4, 1, Number>
1469  value(const ::SymmetricTensor<4, 1, Number> &t)
1470  {
1472  tmp.data[0][0] = 1.0 / t.data[0][0];
1473  return tmp;
1474  }
1475  };
1476 
1477 
1478  template <typename Number>
1479  struct Inverse<4, 2, Number>
1480  {
1481  constexpr static inline ::SymmetricTensor<4, 2, Number>
1482  value(const ::SymmetricTensor<4, 2, Number> &t)
1483  {
1485 
1486  // Inverting this tensor is a little more complicated than necessary,
1487  // since we store the data of 't' as a 3x3 matrix t.data, but the
1488  // product between a rank-4 and a rank-2 tensor is really not the
1489  // product between this matrix and the 3-vector of a rhs, but rather
1490  //
1491  // B.vec = t.data * mult * A.vec
1492  //
1493  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1494  // capture the fact that we need to add up both the c_ij12*a_12 and the
1495  // c_ij21*a_21 terms.
1496  //
1497  // In addition, in this scheme, the identity tensor has the matrix
1498  // representation mult^-1.
1499  //
1500  // The inverse of 't' therefore has the matrix representation
1501  //
1502  // inv.data = mult^-1 * t.data^-1 * mult^-1
1503  //
1504  // in order to compute it, let's first compute the inverse of t.data and
1505  // put it into tmp.data; at the end of the function we then scale the
1506  // last row and column of the inverse by 1/2, corresponding to the left
1507  // and right multiplication with mult^-1.
1508  const Number t4 = t.data[0][0] * t.data[1][1],
1509  t6 = t.data[0][0] * t.data[1][2],
1510  t8 = t.data[0][1] * t.data[1][0],
1511  t00 = t.data[0][2] * t.data[1][0],
1512  t01 = t.data[0][1] * t.data[2][0],
1513  t04 = t.data[0][2] * t.data[2][0],
1514  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1515  t8 * t.data[2][2] + t00 * t.data[2][1] +
1516  t01 * t.data[1][2] - t04 * t.data[1][1]);
1517  tmp.data[0][0] =
1518  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1519  tmp.data[0][1] =
1520  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1521  tmp.data[0][2] =
1522  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1523  tmp.data[1][0] =
1524  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1525  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1526  tmp.data[1][2] = -(t6 - t00) * t07;
1527  tmp.data[2][0] =
1528  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1529  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1530  tmp.data[2][2] = (t4 - t8) * t07;
1531 
1532  // scale last row and column as mentioned
1533  // above
1534  tmp.data[2][0] /= 2;
1535  tmp.data[2][1] /= 2;
1536  tmp.data[0][2] /= 2;
1537  tmp.data[1][2] /= 2;
1538  tmp.data[2][2] /= 4;
1539 
1540  return tmp;
1541  }
1542  };
1543 
1544 
1545  template <typename Number>
1546  struct Inverse<4, 3, Number>
1547  {
1548  static ::SymmetricTensor<4, 3, Number>
1549  value(const ::SymmetricTensor<4, 3, Number> &t)
1550  {
1552 
1553  // This function follows the exact same scheme as the 2d case, except
1554  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1555  // Instead, we use the Gauss-Jordan algorithm implemented for
1556  // FullMatrix. For historical reasons the following code is copied from
1557  // there, with the tangential benefit that we do not need to copy the
1558  // tensor entries to and from the FullMatrix.
1559  const unsigned int N = 6;
1560 
1561  // First get an estimate of the size of the elements of this matrix,
1562  // for later checks whether the pivot element is large enough, or
1563  // whether we have to fear that the matrix is not regular.
1564  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1565  for (unsigned int i = 0; i < N; ++i)
1566  diagonal_sum += numbers::NumberTraits<Number>::abs(tmp.data[i][i]);
1567  const Number typical_diagonal_element =
1568  diagonal_sum / static_cast<double>(N);
1569  (void)typical_diagonal_element;
1570 
1571  unsigned int p[N];
1572  for (unsigned int i = 0; i < N; ++i)
1573  p[i] = i;
1574 
1575  for (unsigned int j = 0; j < N; ++j)
1576  {
1577  // Pivot search: search that part of the line on and right of the
1578  // diagonal for the largest element.
1579  Number max = numbers::NumberTraits<Number>::abs(tmp.data[j][j]);
1580  unsigned int r = j;
1581  for (unsigned int i = j + 1; i < N; ++i)
1582  if (numbers::NumberTraits<Number>::abs(tmp.data[i][j]) > max)
1583  {
1585  r = i;
1586  }
1587 
1588  // Check whether the pivot is too small
1589  Assert(max > 1.e-16 * typical_diagonal_element,
1590  ExcMessage("This tensor seems to be noninvertible"));
1591 
1592  // Row interchange
1593  if (r > j)
1594  {
1595  for (unsigned int k = 0; k < N; ++k)
1596  std::swap(tmp.data[j][k], tmp.data[r][k]);
1597 
1598  std::swap(p[j], p[r]);
1599  }
1600 
1601  // Transformation
1602  const Number hr = 1. / tmp.data[j][j];
1603  tmp.data[j][j] = hr;
1604  for (unsigned int k = 0; k < N; ++k)
1605  {
1606  if (k == j)
1607  continue;
1608  for (unsigned int i = 0; i < N; ++i)
1609  {
1610  if (i == j)
1611  continue;
1612  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1613  }
1614  }
1615  for (unsigned int i = 0; i < N; ++i)
1616  {
1617  tmp.data[i][j] *= hr;
1618  tmp.data[j][i] *= -hr;
1619  }
1620  tmp.data[j][j] = hr;
1621  }
1622 
1623  // Column interchange
1624  Number hv[N];
1625  for (unsigned int i = 0; i < N; ++i)
1626  {
1627  for (unsigned int k = 0; k < N; ++k)
1628  hv[p[k]] = tmp.data[i][k];
1629  for (unsigned int k = 0; k < N; ++k)
1630  tmp.data[i][k] = hv[k];
1631  }
1632 
1633  // Scale rows and columns. The mult matrix
1634  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1635  for (unsigned int i = 3; i < 6; ++i)
1636  for (unsigned int j = 0; j < 3; ++j)
1637  tmp.data[i][j] /= 2;
1638 
1639  for (unsigned int i = 0; i < 3; ++i)
1640  for (unsigned int j = 3; j < 6; ++j)
1641  tmp.data[i][j] /= 2;
1642 
1643  for (unsigned int i = 3; i < 6; ++i)
1644  for (unsigned int j = 3; j < 6; ++j)
1645  tmp.data[i][j] /= 4;
1646 
1647  return tmp;
1648  }
1649  };
1650 
1651  } // namespace SymmetricTensorImplementation
1652 } // namespace internal
1653 
1654 
1655 
1656 template <int rank_, int dim, typename Number>
1657 constexpr DEAL_II_ALWAYS_INLINE
1659  const
1660 {
1661  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1662 }
1663 
1664 
1665 
1666 template <int rank_, int dim, typename Number>
1667 constexpr bool
1669  const SymmetricTensor<rank_, dim, Number> &t) const
1670 {
1671  return data == t.data;
1672 }
1673 
1674 
1675 
1676 template <int rank_, int dim, typename Number>
1677 constexpr bool
1679  const SymmetricTensor<rank_, dim, Number> &t) const
1680 {
1681  return data != t.data;
1682 }
1683 
1684 
1685 
1686 template <int rank_, int dim, typename Number>
1687 template <typename OtherNumber>
1691 {
1692  data += t.data;
1693  return *this;
1694 }
1695 
1696 
1697 
1698 template <int rank_, int dim, typename Number>
1699 template <typename OtherNumber>
1703 {
1704  data -= t.data;
1705  return *this;
1706 }
1707 
1708 
1709 
1710 template <int rank_, int dim, typename Number>
1711 template <typename OtherNumber>
1714 {
1715  data *= d;
1716  return *this;
1717 }
1718 
1719 
1720 
1721 template <int rank_, int dim, typename Number>
1722 template <typename OtherNumber>
1725 {
1726  data /= d;
1727  return *this;
1728 }
1729 
1730 
1731 
1732 template <int rank_, int dim, typename Number>
1735 {
1736  SymmetricTensor tmp = *this;
1737  tmp.data = -tmp.data;
1738  return tmp;
1739 }
1740 
1741 
1742 
1743 template <int rank_, int dim, typename Number>
1744 constexpr inline DEAL_II_ALWAYS_INLINE void
1746 {
1747  data.clear();
1748 }
1749 
1750 
1751 
1752 template <int rank_, int dim, typename Number>
1753 constexpr std::size_t
1755 {
1756  // all memory consists of statically allocated memory of the current
1757  // object, no pointers
1758  return sizeof(SymmetricTensor<rank_, dim, Number>);
1759 }
1760 
1761 
1762 
1763 namespace internal
1764 {
1765  template <int dim, typename Number, typename OtherNumber = Number>
1769  perform_double_contraction(
1770  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1771  base_tensor_type &data,
1772  const typename SymmetricTensorAccessors::
1773  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1774  {
1775  using result_type = typename SymmetricTensorAccessors::
1777 
1778  switch (dim)
1779  {
1780  case 1:
1781  return data[0] * sdata[0];
1782  default:
1783  // Start with the non-diagonal part to avoid some multiplications by
1784  // 2.
1785 
1786  result_type sum = data[dim] * sdata[dim];
1787  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1788  sum += data[d] * sdata[d];
1789  sum += sum; // sum = sum * 2.;
1790 
1791  // Now add the contributions from the diagonal
1792  for (unsigned int d = 0; d < dim; ++d)
1793  sum += data[d] * sdata[d];
1794  return sum;
1795  }
1796  }
1797 
1798 
1799 
1800  template <int dim, typename Number, typename OtherNumber = Number>
1804  perform_double_contraction(
1805  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1806  base_tensor_type &data,
1807  const typename SymmetricTensorAccessors::
1808  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1809  {
1810  using result_type = typename SymmetricTensorAccessors::
1812  using value_type = typename SymmetricTensorAccessors::
1814 
1815  const unsigned int data_dim = SymmetricTensorAccessors::
1816  StorageType<2, dim, value_type>::n_independent_components;
1817  value_type tmp[data_dim]{};
1818  for (unsigned int i = 0; i < data_dim; ++i)
1819  tmp[i] =
1820  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1821  return result_type(tmp);
1822  }
1823 
1824 
1825 
1826  template <int dim, typename Number, typename OtherNumber = Number>
1828  typename SymmetricTensorAccessors::StorageType<
1829  2,
1830  dim,
1833  base_tensor_type
1834  perform_double_contraction(
1835  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1836  base_tensor_type &data,
1837  const typename SymmetricTensorAccessors::
1838  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1839  {
1840  using value_type = typename SymmetricTensorAccessors::
1842  using base_tensor_type = typename SymmetricTensorAccessors::
1843  StorageType<2, dim, value_type>::base_tensor_type;
1844 
1845  base_tensor_type tmp;
1846  for (unsigned int i = 0; i < tmp.dimension; ++i)
1847  {
1848  // Start with the non-diagonal part
1849  value_type sum = data[dim] * sdata[dim][i];
1850  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1851  sum += data[d] * sdata[d][i];
1852  sum += sum; // sum = sum * 2.;
1853 
1854  // Now add the contributions from the diagonal
1855  for (unsigned int d = 0; d < dim; ++d)
1856  sum += data[d] * sdata[d][i];
1857  tmp[i] = sum;
1858  }
1859  return tmp;
1860  }
1861 
1862 
1863 
1864  template <int dim, typename Number, typename OtherNumber = Number>
1866  typename SymmetricTensorAccessors::StorageType<
1867  4,
1868  dim,
1871  base_tensor_type
1872  perform_double_contraction(
1873  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1874  base_tensor_type &data,
1875  const typename SymmetricTensorAccessors::
1876  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1877  {
1878  using value_type = typename SymmetricTensorAccessors::
1880  using base_tensor_type = typename SymmetricTensorAccessors::
1881  StorageType<4, dim, value_type>::base_tensor_type;
1882 
1883  const unsigned int data_dim = SymmetricTensorAccessors::
1884  StorageType<2, dim, value_type>::n_independent_components;
1885  base_tensor_type tmp;
1886  for (unsigned int i = 0; i < data_dim; ++i)
1887  for (unsigned int j = 0; j < data_dim; ++j)
1888  {
1889  // Start with the non-diagonal part
1890  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1891  tmp[i][j] += data[i][d] * sdata[d][j];
1892  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1893 
1894  // Now add the contributions from the diagonal
1895  for (unsigned int d = 0; d < dim; ++d)
1896  tmp[i][j] += data[i][d] * sdata[d][j];
1897  }
1898  return tmp;
1899  }
1900 
1901 } // end of namespace internal
1902 
1903 
1904 
1905 template <int rank_, int dim, typename Number>
1906 template <typename OtherNumber>
1911  const SymmetricTensor<2, dim, OtherNumber> &s) const
1912 {
1913  // need to have two different function calls
1914  // because a scalar and rank-2 tensor are not
1915  // the same data type (see internal function
1916  // above)
1917  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1918  s.data);
1919 }
1920 
1921 
1922 
1923 template <int rank_, int dim, typename Number>
1924 template <typename OtherNumber>
1928  const SymmetricTensor<4, dim, OtherNumber> &s) const
1929 {
1932  tmp.data =
1933  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1934  s.data);
1935  return tmp;
1936 }
1937 
1938 
1939 
1940 // internal namespace to switch between the
1941 // access of different tensors. There used to
1942 // be explicit instantiations before for
1943 // different ranks and dimensions, but since
1944 // we now allow for templates on the data
1945 // type, and since we cannot partially
1946 // specialize the implementation, this got
1947 // into a separate namespace
1948 namespace internal
1949 {
1950  // The variables within this struct will be referenced in the next functions.
1951  // It is a workaround that allows returning a reference to a static variable
1952  // while allowing constexpr evaluation of the function.
1953  // It has to be defined outside the function because constexpr functions
1954  // cannot define static variables.
1955  // A similar struct has also been defined in tensor.h
1956  template <typename Type>
1957  struct Uninitialized
1958  {
1959  static Type value;
1960  };
1961 
1962  template <typename Type>
1963  Type Uninitialized<Type>::value;
1964 
1965  template <int dim, typename Number>
1966  constexpr inline DEAL_II_ALWAYS_INLINE Number &
1967  symmetric_tensor_access(const TableIndices<2> &indices,
1968  typename SymmetricTensorAccessors::
1969  StorageType<2, dim, Number>::base_tensor_type &data)
1970  {
1971  // 1d is very simple and done first
1972  if (dim == 1)
1973  return data[0];
1974 
1975  // first treat the main diagonal elements, which are stored consecutively
1976  // at the beginning
1977  if (indices[0] == indices[1])
1978  return data[indices[0]];
1979 
1980  // the rest is messier and requires a few switches.
1981  switch (dim)
1982  {
1983  case 2:
1984  // at least for the 2x2 case it is reasonably simple
1985  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1986  ((indices[0] == 0) && (indices[1] == 1)),
1987  ExcInternalError());
1988  return data[2];
1989 
1990  default:
1991  // to do the rest, sort our indices before comparing
1992  {
1993  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1994  std::max(indices[0], indices[1]));
1995  for (unsigned int d = 0, c = 0; d < dim; ++d)
1996  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1997  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1998  return data[dim + c];
1999  Assert(false, ExcInternalError());
2000  }
2001  }
2002 
2003  // The code should never reach there.
2004  // Returns a dummy reference to a dummy variable just to make the
2005  // compiler happy.
2006  return Uninitialized<Number>::value;
2007  }
2008 
2009 
2010 
2011  template <int dim, typename Number>
2012  constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2013  symmetric_tensor_access(const TableIndices<2> &indices,
2014  const typename SymmetricTensorAccessors::
2015  StorageType<2, dim, Number>::base_tensor_type &data)
2016  {
2017  // 1d is very simple and done first
2018  if (dim == 1)
2019  return data[0];
2020 
2021  // first treat the main diagonal elements, which are stored consecutively
2022  // at the beginning
2023  if (indices[0] == indices[1])
2024  return data[indices[0]];
2025 
2026  // the rest is messier and requires a few switches.
2027  switch (dim)
2028  {
2029  case 2:
2030  // at least for the 2x2 case it is reasonably simple
2031  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
2032  ((indices[0] == 0) && (indices[1] == 1)),
2033  ExcInternalError());
2034  return data[2];
2035 
2036  default:
2037  // to do the rest, sort our indices before comparing
2038  {
2039  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2040  std::max(indices[0], indices[1]));
2041  for (unsigned int d = 0, c = 0; d < dim; ++d)
2042  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2043  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2044  return data[dim + c];
2045  Assert(false, ExcInternalError());
2046  }
2047  }
2048 
2049  // The code should never reach there.
2050  // Returns a dummy reference to a dummy variable just to make the
2051  // compiler happy.
2052  return Uninitialized<Number>::value;
2053  }
2054 
2055 
2056 
2057  template <int dim, typename Number>
2058  constexpr inline Number &
2059  symmetric_tensor_access(const TableIndices<4> &indices,
2060  typename SymmetricTensorAccessors::
2061  StorageType<4, dim, Number>::base_tensor_type &data)
2062  {
2063  switch (dim)
2064  {
2065  case 1:
2066  return data[0][0];
2067 
2068  case 2:
2069  // each entry of the tensor can be thought of as an entry in a
2070  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2071  // rank-2 tensors. this is the format in which we store rank-4
2072  // tensors. determine which position the present entry is
2073  // stored in
2074  {
2075  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2076  return data[base_index[indices[0]][indices[1]]]
2077  [base_index[indices[2]][indices[3]]];
2078  }
2079  case 3:
2080  // each entry of the tensor can be thought of as an entry in a
2081  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2082  // rank-2 tensors. this is the format in which we store rank-4
2083  // tensors. determine which position the present entry is
2084  // stored in
2085  {
2086  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2087  {3, 1, 5},
2088  {4, 5, 2}};
2089  return data[base_index[indices[0]][indices[1]]]
2090  [base_index[indices[2]][indices[3]]];
2091  }
2092 
2093  default:
2094  Assert(false, ExcNotImplemented());
2095  }
2096 
2097  // The code should never reach there.
2098  // Returns a dummy reference to a dummy variable just to make the
2099  // compiler happy.
2100  return Uninitialized<Number>::value;
2101  }
2102 
2103 
2104  template <int dim, typename Number>
2105  constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2106  symmetric_tensor_access(const TableIndices<4> &indices,
2107  const typename SymmetricTensorAccessors::
2108  StorageType<4, dim, Number>::base_tensor_type &data)
2109  {
2110  switch (dim)
2111  {
2112  case 1:
2113  return data[0][0];
2114 
2115  case 2:
2116  // each entry of the tensor can be thought of as an entry in a
2117  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2118  // rank-2 tensors. this is the format in which we store rank-4
2119  // tensors. determine which position the present entry is
2120  // stored in
2121  {
2122  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2123  return data[base_index[indices[0]][indices[1]]]
2124  [base_index[indices[2]][indices[3]]];
2125  }
2126  case 3:
2127  // each entry of the tensor can be thought of as an entry in a
2128  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2129  // rank-2 tensors. this is the format in which we store rank-4
2130  // tensors. determine which position the present entry is
2131  // stored in
2132  {
2133  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2134  {3, 1, 5},
2135  {4, 5, 2}};
2136  return data[base_index[indices[0]][indices[1]]]
2137  [base_index[indices[2]][indices[3]]];
2138  }
2139 
2140  default:
2141  Assert(false, ExcNotImplemented());
2142  }
2143 
2144  // The code should never reach there.
2145  // Returns a dummy reference to a dummy variable just to make the
2146  // compiler happy.
2147  return Uninitialized<Number>::value;
2148  }
2149 
2150 } // end of namespace internal
2151 
2152 
2153 
2154 template <int rank_, int dim, typename Number>
2155 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2157  const TableIndices<rank_> &indices)
2158 {
2159  for (unsigned int r = 0; r < rank; ++r)
2160  AssertIndexRange(indices[r], dimension);
2161  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2162 }
2163 
2164 
2165 
2166 template <int rank_, int dim, typename Number>
2167 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2169  const TableIndices<rank_> &indices) const
2170 {
2171  for (unsigned int r = 0; r < rank; ++r)
2172  AssertIndexRange(indices[r], dimension);
2173  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2174 }
2175 
2176 
2177 
2178 namespace internal
2179 {
2180  namespace SymmetricTensorImplementation
2181  {
2182  template <int rank_>
2183  constexpr TableIndices<rank_>
2184  get_partially_filled_indices(const unsigned int row,
2185  const std::integral_constant<int, 2> &)
2186  {
2188  }
2189 
2190 
2191  template <int rank_>
2192  constexpr TableIndices<rank_>
2193  get_partially_filled_indices(const unsigned int row,
2194  const std::integral_constant<int, 4> &)
2195  {
2196  return TableIndices<rank_>(row,
2200  }
2201  } // namespace SymmetricTensorImplementation
2202 } // namespace internal
2203 
2204 
2205 template <int rank_, int dim, typename Number>
2206 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2207  Accessor<rank_, dim, true, rank_ - 1, Number>
2208  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2209 {
2210  return internal::SymmetricTensorAccessors::
2211  Accessor<rank_, dim, true, rank_ - 1, Number>(
2212  *this,
2213  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2214  rank_>(row, std::integral_constant<int, rank_>()));
2215 }
2216 
2217 
2218 
2219 template <int rank_, int dim, typename Number>
2220 constexpr inline DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2221  Accessor<rank_, dim, false, rank_ - 1, Number>
2222  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2223 {
2224  return internal::SymmetricTensorAccessors::
2225  Accessor<rank_, dim, false, rank_ - 1, Number>(
2226  *this,
2227  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2228  rank_>(row, std::integral_constant<int, rank_>()));
2229 }
2230 
2231 
2232 
2233 template <int rank_, int dim, typename Number>
2234 constexpr DEAL_II_ALWAYS_INLINE const Number &
2236  const TableIndices<rank_> &indices) const
2237 {
2238  return operator()(indices);
2239 }
2240 
2241 
2242 
2243 template <int rank_, int dim, typename Number>
2244 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2246  const TableIndices<rank_> &indices)
2247 {
2248  return operator()(indices);
2249 }
2250 
2251 
2252 
2253 template <int rank_, int dim, typename Number>
2254 inline Number *
2256 {
2257  return std::addressof(this->access_raw_entry(0));
2258 }
2259 
2260 
2261 
2262 template <int rank_, int dim, typename Number>
2263 inline const Number *
2265 {
2266  return std::addressof(this->access_raw_entry(0));
2267 }
2268 
2269 
2270 
2271 template <int rank_, int dim, typename Number>
2272 inline Number *
2274 {
2275  return begin_raw() + n_independent_components;
2276 }
2277 
2278 
2279 
2280 template <int rank_, int dim, typename Number>
2281 inline const Number *
2283 {
2284  return begin_raw() + n_independent_components;
2285 }
2286 
2287 
2288 
2289 namespace internal
2290 {
2291  namespace SymmetricTensorImplementation
2292  {
2293  template <int dim, typename Number>
2294  constexpr unsigned int
2295  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2296  const unsigned int index)
2297  {
2298  return index;
2299  }
2300 
2301 
2302  template <int dim, typename Number>
2303  constexpr ::TableIndices<2>
2304  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2305  const unsigned int index)
2306  {
2309  }
2310 
2311  } // namespace SymmetricTensorImplementation
2312 } // namespace internal
2313 
2314 
2315 
2316 template <int rank_, int dim, typename Number>
2317 constexpr inline const Number &
2319  const unsigned int index) const
2320 {
2321  AssertIndexRange(index, n_independent_components);
2322  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2323  index)];
2324 }
2325 
2326 
2327 
2328 template <int rank_, int dim, typename Number>
2329 constexpr inline Number &
2331 {
2332  AssertIndexRange(index, n_independent_components);
2333  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2334  index)];
2335 }
2336 
2337 
2338 
2339 namespace internal
2340 {
2341  template <int dim, typename Number>
2342  constexpr inline typename numbers::NumberTraits<Number>::real_type
2343  compute_norm(const typename SymmetricTensorAccessors::
2344  StorageType<2, dim, Number>::base_tensor_type &data)
2345  {
2346  // Make things work with AD types
2347  using std::sqrt;
2348  switch (dim)
2349  {
2350  case 1:
2351  return numbers::NumberTraits<Number>::abs(data[0]);
2352 
2353  case 2:
2357 
2358  case 3:
2365 
2366  default:
2367  {
2368  typename numbers::NumberTraits<Number>::real_type return_value =
2370 
2371  for (unsigned int d = 0; d < dim; ++d)
2372  return_value +=
2374  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2375  return_value +=
2377 
2378  return sqrt(return_value);
2379  }
2380  }
2381  }
2382 
2383 
2384 
2385  template <int dim, typename Number>
2386  constexpr inline typename numbers::NumberTraits<Number>::real_type
2387  compute_norm(const typename SymmetricTensorAccessors::
2388  StorageType<4, dim, Number>::base_tensor_type &data)
2389  {
2390  // Make things work with AD types
2391  using std::sqrt;
2392  switch (dim)
2393  {
2394  case 1:
2395  return numbers::NumberTraits<Number>::abs(data[0][0]);
2396 
2397  default:
2398  {
2399  typename numbers::NumberTraits<Number>::real_type return_value =
2401 
2402  const unsigned int n_independent_components = data.dimension;
2403 
2404  for (unsigned int i = 0; i < dim; ++i)
2405  for (unsigned int j = 0; j < dim; ++j)
2406  return_value +=
2408  for (unsigned int i = 0; i < dim; ++i)
2409  for (unsigned int j = dim; j < n_independent_components; ++j)
2410  return_value +=
2412  for (unsigned int i = dim; i < n_independent_components; ++i)
2413  for (unsigned int j = 0; j < dim; ++j)
2414  return_value +=
2416  for (unsigned int i = dim; i < n_independent_components; ++i)
2417  for (unsigned int j = dim; j < n_independent_components; ++j)
2418  return_value +=
2420 
2421  return sqrt(return_value);
2422  }
2423  }
2424  }
2425 
2426 } // end of namespace internal
2427 
2428 
2429 
2430 template <int rank_, int dim, typename Number>
2433 {
2434  return internal::compute_norm<dim, Number>(data);
2435 }
2436 
2437 
2438 
2439 namespace internal
2440 {
2441  namespace SymmetricTensorImplementation
2442  {
2443  // a function to do the unrolling from a set of indices to a
2444  // scalar index into the array in which we store the elements of
2445  // a symmetric tensor
2446  //
2447  // this function is for rank-2 tensors
2448  template <int dim>
2449  constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
2451  {
2452  AssertIndexRange(indices[0], dim);
2453  AssertIndexRange(indices[1], dim);
2454 
2455  switch (dim)
2456  {
2457  case 1:
2458  {
2459  return 0;
2460  }
2461 
2462  case 2:
2463  {
2464  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2465  return table[indices[0]][indices[1]];
2466  }
2467 
2468  case 3:
2469  {
2470  constexpr unsigned int table[3][3] = {{0, 3, 4},
2471  {3, 1, 5},
2472  {4, 5, 2}};
2473  return table[indices[0]][indices[1]];
2474  }
2475 
2476  case 4:
2477  {
2478  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2479  {4, 1, 7, 8},
2480  {5, 7, 2, 9},
2481  {6, 8, 9, 3}};
2482  return table[indices[0]][indices[1]];
2483  }
2484 
2485  default:
2486  // for the remainder, manually figure out the numbering
2487  {
2488  if (indices[0] == indices[1])
2489  return indices[0];
2490 
2491  TableIndices<2> sorted_indices(indices);
2492  sorted_indices.sort();
2493 
2494  for (unsigned int d = 0, c = 0; d < dim; ++d)
2495  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2496  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2497  return dim + c;
2498 
2499  // should never get here:
2500  Assert(false, ExcInternalError());
2501  return 0;
2502  }
2503  }
2504  }
2505 
2506  // a function to do the unrolling from a set of indices to a
2507  // scalar index into the array in which we store the elements of
2508  // a symmetric tensor
2509  //
2510  // this function is for tensors of ranks not already handled
2511  // above
2512  template <int dim, int rank_>
2513  constexpr inline unsigned int
2515  {
2516  (void)indices;
2517  Assert(false, ExcNotImplemented());
2519  }
2520  } // namespace SymmetricTensorImplementation
2521 } // namespace internal
2522 
2523 
2524 template <int rank_, int dim, typename Number>
2525 constexpr unsigned int
2527  const TableIndices<rank_> &indices)
2528 {
2529  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2530  dim>(indices);
2531 }
2532 
2533 
2534 
2535 namespace internal
2536 {
2537  namespace SymmetricTensorImplementation
2538  {
2539  // a function to do the inverse of the unrolling from a set of
2540  // indices to a scalar index into the array in which we store
2541  // the elements of a symmetric tensor. in other words, it goes
2542  // from the scalar index into the array to a set of indices of
2543  // the tensor
2544  //
2545  // this function is for rank-2 tensors
2546  template <int dim>
2547  constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2548  unrolled_to_component_indices(const unsigned int i,
2549  const std::integral_constant<int, 2> &)
2550  {
2551  Assert(
2553  ExcIndexRange(
2554  i,
2555  0,
2557  switch (dim)
2558  {
2559  case 1:
2560  {
2561  return {0, 0};
2562  }
2563 
2564  case 2:
2565  {
2566  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2567  TableIndices<2>(1, 1),
2568  TableIndices<2>(0, 1)};
2569  return table[i];
2570  }
2571 
2572  case 3:
2573  {
2574  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2575  TableIndices<2>(1, 1),
2576  TableIndices<2>(2, 2),
2577  TableIndices<2>(0, 1),
2578  TableIndices<2>(0, 2),
2579  TableIndices<2>(1, 2)};
2580  return table[i];
2581  }
2582 
2583  default:
2584  if (i < dim)
2585  return {i, i};
2586 
2587  for (unsigned int d = 0, c = dim; d < dim; ++d)
2588  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2589  if (c == i)
2590  return {d, e};
2591 
2592  // should never get here:
2593  Assert(false, ExcInternalError());
2594  return {0, 0};
2595  }
2596  }
2597 
2598  // a function to do the inverse of the unrolling from a set of
2599  // indices to a scalar index into the array in which we store
2600  // the elements of a symmetric tensor. in other words, it goes
2601  // from the scalar index into the array to a set of indices of
2602  // the tensor
2603  //
2604  // this function is for tensors of a rank not already handled
2605  // above
2606  template <int dim, int rank_>
2607  constexpr inline std::enable_if_t<rank_ != 2, TableIndices<rank_>>
2608  unrolled_to_component_indices(const unsigned int i,
2609  const std::integral_constant<int, rank_> &)
2610  {
2611  (void)i;
2612  Assert(
2613  (i <
2615  ExcIndexRange(i,
2616  0,
2618  n_independent_components));
2619  Assert(false, ExcNotImplemented());
2620  return TableIndices<rank_>();
2621  }
2622 
2623  } // namespace SymmetricTensorImplementation
2624 } // namespace internal
2625 
2626 template <int rank_, int dim, typename Number>
2629  const unsigned int i)
2630 {
2631  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2632  dim>(i, std::integral_constant<int, rank_>());
2633 }
2634 
2635 
2636 
2637 template <int rank_, int dim, typename Number>
2638 template <class Archive>
2639 inline void
2640 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2641 {
2642  ar &data;
2643 }
2644 
2645 
2646 #endif // DOXYGEN
2647 
2648 /* ----------------- Non-member functions operating on tensors. ------------ */
2649 
2650 
2663 template <int rank_, int dim, typename Number, typename OtherNumber>
2664 constexpr inline DEAL_II_ALWAYS_INLINE
2668 {
2670  tmp = left;
2671  tmp += right;
2672  return tmp;
2673 }
2674 
2675 
2688 template <int rank_, int dim, typename Number, typename OtherNumber>
2689 constexpr inline DEAL_II_ALWAYS_INLINE
2693 {
2695  tmp = left;
2696  tmp -= right;
2697  return tmp;
2698 }
2699 
2700 
2708 template <int rank_, int dim, typename Number, typename OtherNumber>
2709 constexpr DEAL_II_ALWAYS_INLINE
2712  const Tensor<rank_, dim, OtherNumber> & right)
2713 {
2714  return Tensor<rank_, dim, Number>(left) + right;
2715 }
2716 
2717 
2725 template <int rank_, int dim, typename Number, typename OtherNumber>
2726 constexpr DEAL_II_ALWAYS_INLINE
2730 {
2731  return left + Tensor<rank_, dim, OtherNumber>(right);
2732 }
2733 
2734 
2742 template <int rank_, int dim, typename Number, typename OtherNumber>
2743 constexpr DEAL_II_ALWAYS_INLINE
2746  const Tensor<rank_, dim, OtherNumber> & right)
2747 {
2748  return Tensor<rank_, dim, Number>(left) - right;
2749 }
2750 
2751 
2759 template <int rank_, int dim, typename Number, typename OtherNumber>
2760 constexpr DEAL_II_ALWAYS_INLINE
2764 {
2765  return left - Tensor<rank_, dim, OtherNumber>(right);
2766 }
2767 
2768 
2769 
2770 template <int dim, typename Number>
2773 {
2774  switch (dim)
2775  {
2776  case 1:
2777  return t.data[0];
2778  case 2:
2779  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2780  case 3:
2781  {
2782  // in analogy to general tensors, but
2783  // there's something to be simplified for
2784  // the present case
2785  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2786  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2787  t.data[0] * t.data[5] * t.data[5] -
2788  t.data[1] * t.data[4] * t.data[4] -
2789  t.data[2] * t.data[3] * t.data[3]);
2790  }
2791  default:
2792  Assert(false, ExcNotImplemented());
2794  }
2795 }
2796 
2797 
2798 
2810 template <int dim, typename Number>
2813 {
2814  return determinant(t);
2815 }
2816 
2817 
2818 
2819 template <int dim, typename Number>
2820 constexpr inline DEAL_II_ALWAYS_INLINE Number
2822 {
2823  Number t = d.data[0];
2824  for (unsigned int i = 1; i < dim; ++i)
2825  t += d.data[i];
2826  return t;
2827 }
2828 
2829 
2841 template <int dim, typename Number>
2842 constexpr Number
2844 {
2845  return trace(t);
2846 }
2847 
2848 
2860 template <typename Number>
2861 constexpr DEAL_II_ALWAYS_INLINE Number
2863 {
2865 }
2866 
2867 
2868 
2887 template <typename Number>
2888 constexpr DEAL_II_ALWAYS_INLINE Number
2890 {
2891  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2892 }
2893 
2894 
2895 
2904 template <typename Number>
2905 constexpr DEAL_II_ALWAYS_INLINE Number
2907 {
2908  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2909  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2910 }
2911 
2912 
2913 
2921 template <typename Number>
2922 std::array<Number, 1>
2924 
2925 
2926 
2949 template <typename Number>
2950 std::array<Number, 2>
2952 
2953 
2954 
2977 template <typename Number>
2978 std::array<Number, 3>
2980 
2981 
2982 
2983 namespace internal
2984 {
2985  namespace SymmetricTensorImplementation
2986  {
3024  template <int dim, typename Number>
3025  void
3026  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
3027  ::Tensor<2, dim, Number> & Q,
3028  std::array<Number, dim> & d,
3029  std::array<Number, dim - 1> & e);
3030 
3031 
3032 
3072  template <int dim, typename Number>
3073  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3074  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
3075 
3076 
3077 
3117  template <int dim, typename Number>
3118  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3120 
3121 
3122 
3136  template <typename Number>
3137  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3138  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3139 
3140 
3141 
3174  template <typename Number>
3175  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3176  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3177 
3182  template <int dim, typename Number>
3184  {
3185  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3186  bool
3187  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3188  {
3189  return lhs.first > rhs.first;
3190  }
3191  };
3192 
3193  } // namespace SymmetricTensorImplementation
3194 
3195 } // namespace internal
3196 
3197 
3198 
3199 // The line below is to ensure that doxygen puts the full description
3200 // of this global enumeration into the documentation
3201 // See https://stackoverflow.com/a/1717984
3231 {
3241  hybrid,
3259  jacobi
3260 };
3261 
3262 
3263 
3292 template <int dim, typename Number>
3293 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3294  std::integral_constant<int, dim>::value>
3296  const SymmetricTensorEigenvectorMethod method =
3298 
3299 
3300 
3309 template <int rank_, int dim, typename Number>
3312 {
3313  return t;
3314 }
3315 
3316 
3317 
3318 template <int dim, typename Number>
3321 {
3323 
3324  // subtract scaled trace from the diagonal
3325  const Number tr = trace(t) / dim;
3326  for (unsigned int i = 0; i < dim; ++i)
3327  tmp.data[i] -= tr;
3328 
3329  return tmp;
3330 }
3331 
3332 
3333 
3334 template <int dim, typename Number>
3337 {
3338  // create a default constructed matrix filled with
3339  // zeros, then set the diagonal elements to one
3341  switch (dim)
3342  {
3343  case 1:
3345  break;
3346  case 2:
3347  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3348  break;
3349  case 3:
3350  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3352  break;
3353  default:
3354  for (unsigned int d = 0; d < dim; ++d)
3356  }
3357  return tmp;
3358 }
3359 
3360 
3361 
3362 template <int dim, typename Number>
3365 {
3367 
3368  // fill the elements treating the diagonal
3369  for (unsigned int i = 0; i < dim; ++i)
3370  for (unsigned int j = 0; j < dim; ++j)
3371  tmp.data[i][j] =
3372  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3373 
3374  // then fill the ones that copy over the
3375  // non-diagonal elements. note that during
3376  // the double-contraction, we handle the
3377  // off-diagonal elements twice, so simply
3378  // copying requires a weight of 1/2
3379  for (unsigned int i = dim;
3380  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3381  n_rank2_components;
3382  ++i)
3383  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3384 
3385  return tmp;
3386 }
3387 
3388 
3389 
3390 template <int dim, typename Number>
3393 {
3395 
3396  // fill the elements treating the diagonal
3397  for (unsigned int i = 0; i < dim; ++i)
3398  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3399 
3400  // then fill the ones that copy over the
3401  // non-diagonal elements. note that during
3402  // the double-contraction, we handle the
3403  // off-diagonal elements twice, so simply
3404  // copying requires a weight of 1/2
3405  for (unsigned int i = dim;
3406  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3407  n_rank2_components;
3408  ++i)
3409  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3410 
3411  return tmp;
3412 }
3413 
3414 
3415 
3425 template <int dim, typename Number>
3428 {
3430  value(t);
3431 }
3432 
3433 
3434 
3445 template <int dim, typename Number>
3448 {
3450  value(t);
3451 }
3452 
3453 
3454 
3476 template <int dim, typename Number>
3477 constexpr inline SymmetricTensor<4, dim, Number>
3480 {
3482 
3483  // fill only the elements really needed
3484  for (unsigned int i = 0; i < dim; ++i)
3485  for (unsigned int j = i; j < dim; ++j)
3486  for (unsigned int k = 0; k < dim; ++k)
3487  for (unsigned int l = k; l < dim; ++l)
3488  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3489 
3490  return tmp;
3491 }
3492 
3493 
3494 
3502 template <int dim, typename Number>
3505 {
3507  for (unsigned int d = 0; d < dim; ++d)
3508  result[d][d] = t[d][d];
3509 
3510  const Number half = internal::NumberType<Number>::value(0.5);
3511  for (unsigned int d = 0; d < dim; ++d)
3512  for (unsigned int e = d + 1; e < dim; ++e)
3513  result[d][e] = (t[d][e] + t[e][d]) * half;
3514  return result;
3515 }
3516 
3517 
3518 
3526 template <int rank_, int dim, typename Number>
3528 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3529 {
3531  tt *= factor;
3532  return tt;
3533 }
3534 
3535 
3536 
3544 template <int rank_, int dim, typename Number>
3546 operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3547 {
3548  // simply forward to the other operator
3549  return t * factor;
3550 }
3551 
3552 
3553 
3578 template <int rank_, int dim, typename Number, typename OtherNumber>
3579 constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3580  rank_,
3581  dim,
3582  typename ProductType<Number,
3583  typename EnableIfScalar<OtherNumber>::type>::type>
3585  const OtherNumber & factor)
3586 {
3587  // form the product. we have to convert the two factors into the final
3588  // type via explicit casts because, for awkward reasons, the C++
3589  // standard committee saw it fit to not define an
3590  // operator*(float,std::complex<double>)
3591  // (as well as with switched arguments and double<->float).
3592  using product_type = typename ProductType<Number, OtherNumber>::type;
3595  return tt;
3596 }
3597 
3598 
3599 
3607 template <int rank_, int dim, typename Number, typename OtherNumber>
3608 constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3609  rank_,
3610  dim,
3611  typename ProductType<OtherNumber,
3612  typename EnableIfScalar<Number>::type>::type>
3613 operator*(const Number & factor,
3615 {
3616  // simply forward to the other operator with switched arguments
3617  return (t * factor);
3618 }
3619 
3620 
3621 
3627 template <int rank_, int dim, typename Number, typename OtherNumber>
3628 constexpr inline SymmetricTensor<
3629  rank_,
3630  dim,
3631  typename ProductType<Number,
3632  typename EnableIfScalar<OtherNumber>::type>::type>
3634  const OtherNumber & factor)
3635 {
3636  using product_type = typename ProductType<Number, OtherNumber>::type;
3639  return tt;
3640 }
3641 
3642 
3643 
3650 template <int rank_, int dim>
3652 operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3653 {
3655  tt *= factor;
3656  return tt;
3657 }
3658 
3659 
3660 
3667 template <int rank_, int dim>
3669 operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3670 {
3672  tt *= factor;
3673  return tt;
3674 }
3675 
3676 
3677 
3683 template <int rank_, int dim>
3684 constexpr inline SymmetricTensor<rank_, dim>
3685 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3686 {
3688  tt /= factor;
3689  return tt;
3690 }
3691 
3701 template <int dim, typename Number, typename OtherNumber>
3705 {
3706  return (t1 * t2);
3707 }
3708 
3709 
3722 template <int dim, typename Number, typename OtherNumber>
3723 constexpr inline DEAL_II_ALWAYS_INLINE
3726  const Tensor<2, dim, OtherNumber> & t2)
3727 {
3729  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3730  for (unsigned int i = 0; i < dim; ++i)
3731  for (unsigned int j = 0; j < dim; ++j)
3732  s += t1[i][j] * t2[i][j];
3733  return s;
3734 }
3735 
3736 
3749 template <int dim, typename Number, typename OtherNumber>
3753 {
3754  return scalar_product(t2, t1);
3755 }
3756 
3757 
3772 template <typename Number, typename OtherNumber>
3773 constexpr inline DEAL_II_ALWAYS_INLINE void
3778 {
3779  tmp[0][0] = t[0][0][0][0] * s[0][0];
3780 }
3781 
3782 
3783 
3798 template <typename Number, typename OtherNumber>
3799 constexpr inline void
3804 {
3805  tmp[0][0] = t[0][0][0][0] * s[0][0];
3806 }
3807 
3808 
3809 
3824 template <typename Number, typename OtherNumber>
3825 constexpr inline void
3830 {
3831  const unsigned int dim = 2;
3832 
3833  for (unsigned int i = 0; i < dim; ++i)
3834  for (unsigned int j = i; j < dim; ++j)
3835  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3836  2 * t[i][j][0][1] * s[0][1];
3837 }
3838 
3839 
3840 
3855 template <typename Number, typename OtherNumber>
3856 constexpr inline void
3861 {
3862  const unsigned int dim = 2;
3863 
3864  for (unsigned int i = 0; i < dim; ++i)
3865  for (unsigned int j = i; j < dim; ++j)
3866  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3867  2 * s[0][1] * t[0][1][i][j];
3868 }
3869 
3870 
3871 
3886 template <typename Number, typename OtherNumber>
3887 constexpr inline void
3892 {
3893  const unsigned int dim = 3;
3894 
3895  for (unsigned int i = 0; i < dim; ++i)
3896  for (unsigned int j = i; j < dim; ++j)
3897  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3898  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3899  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3900 }
3901 
3902 
3903 
3918 template <typename Number, typename OtherNumber>
3919 constexpr inline void
3924 {
3925  const unsigned int dim = 3;
3926 
3927  for (unsigned int i = 0; i < dim; ++i)
3928  for (unsigned int j = i; j < dim; ++j)
3929  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3930  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3931  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3932 }
3933 
3934 
3935 
3942 template <int dim, typename Number, typename OtherNumber>
3945  const Tensor<1, dim, OtherNumber> & src2)
3946 {
3948  for (unsigned int i = 0; i < dim; ++i)
3949  for (unsigned int j = 0; j < dim; ++j)
3950  dest[i] += src1[i][j] * src2[j];
3951  return dest;
3952 }
3953 
3954 
3961 template <int dim, typename Number, typename OtherNumber>
3965 {
3966  // this is easy for symmetric tensors:
3967  return src2 * src1;
3968 }
3969 
3970 
3971 
3991 template <int rank_1,
3992  int rank_2,
3993  int dim,
3994  typename Number,
3995  typename OtherNumber>
3996 constexpr DEAL_II_ALWAYS_INLINE
3997  typename Tensor<rank_1 + rank_2 - 2,
3998  dim,
3999  typename ProductType<Number, OtherNumber>::type>::tensor_type
4002 {
4003  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4004 }
4005 
4006 
4007 
4027 template <int rank_1,
4028  int rank_2,
4029  int dim,
4030  typename Number,
4031  typename OtherNumber>
4032 constexpr DEAL_II_ALWAYS_INLINE
4033  typename Tensor<rank_1 + rank_2 - 2,
4034  dim,
4035  typename ProductType<Number, OtherNumber>::type>::tensor_type
4037  const Tensor<rank_2, dim, OtherNumber> & src2)
4038 {
4039  return Tensor<rank_1, dim, Number>(src1) * src2;
4040 }
4041 
4042 
4043 
4053 template <int dim, typename Number>
4054 inline std::ostream &
4055 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4056 {
4057  // make our lives a bit simpler by outputting
4058  // the tensor through the operator for the
4059  // general Tensor class
4061 
4062  for (unsigned int i = 0; i < dim; ++i)
4063  for (unsigned int j = 0; j < dim; ++j)
4064  tt[i][j] = t[i][j];
4065 
4066  return out << tt;
4067 }
4068 
4069 
4070 
4080 template <int dim, typename Number>
4081 inline std::ostream &
4082 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4083 {
4084  // make our lives a bit simpler by outputting
4085  // the tensor through the operator for the
4086  // general Tensor class
4088 
4089  for (unsigned int i = 0; i < dim; ++i)
4090  for (unsigned int j = 0; j < dim; ++j)
4091  for (unsigned int k = 0; k < dim; ++k)
4092  for (unsigned int l = 0; l < dim; ++l)
4093  tt[i][j][k][l] = t[i][j][k][l];
4094 
4095  return out << tt;
4096 }
4097 
4098 
4100 
4101 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
constexpr const Number & operator()(const TableIndices< rank_ > &indices) const
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
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 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 Number determinant(const SymmetricTensor< 2, dim, Number > &)
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)
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)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr numbers::NumberTraits< Number >::real_type norm() const
typename base_tensor_descriptor::base_tensor_type base_tensor_type
constexpr SymmetricTensor< 4, dim, Number > 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< 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)
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 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 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 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 SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
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 SymmetricTensor< 4, dim, Number > 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 SymmetricTensor< 2, dim, Number > 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_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_CONSTEXPR
Definition: config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
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:970
static const unsigned int invalid_unsigned_int
Definition: types.h:206
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:733
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:624
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
SymmetricTensorEigenvectorMethod
constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()