Reference documentation for deal.II version GIT d77e5ebb0a 2023-01-27 22:35: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  template <int dim, typename Number>
1951  constexpr inline DEAL_II_ALWAYS_INLINE Number &
1952  symmetric_tensor_access(const TableIndices<2> &indices,
1953  typename SymmetricTensorAccessors::
1954  StorageType<2, dim, Number>::base_tensor_type &data)
1955  {
1956  // 1d is very simple and done first
1957  if (dim == 1)
1958  return data[0];
1959 
1960  // first treat the main diagonal elements, which are stored consecutively
1961  // at the beginning
1962  if (indices[0] == indices[1])
1963  return data[indices[0]];
1964 
1965  // the rest is messier and requires a few switches.
1966  switch (dim)
1967  {
1968  case 2:
1969  // at least for the 2x2 case it is reasonably simple
1970  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1971  ((indices[0] == 0) && (indices[1] == 1)),
1972  ExcInternalError());
1973  return data[2];
1974 
1975  default:
1976  // to do the rest, sort our indices before comparing
1977  {
1978  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1979  std::max(indices[0], indices[1]));
1980  for (unsigned int d = 0, c = 0; d < dim; ++d)
1981  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1982  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1983  return data[dim + c];
1984  Assert(false, ExcInternalError());
1985  }
1986  }
1987 
1988  // The code should never reach here.
1989  // We cannot return a static variable, as this class must support number
1990  // types that require no instances of the number type to be in scope during
1991  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
1992  return data[0];
1993  }
1994 
1995 
1996 
1997  template <int dim, typename Number>
1998  constexpr inline DEAL_II_ALWAYS_INLINE const Number &
1999  symmetric_tensor_access(const TableIndices<2> &indices,
2000  const typename SymmetricTensorAccessors::
2001  StorageType<2, dim, Number>::base_tensor_type &data)
2002  {
2003  // 1d is very simple and done first
2004  if (dim == 1)
2005  return data[0];
2006 
2007  // first treat the main diagonal elements, which are stored consecutively
2008  // at the beginning
2009  if (indices[0] == indices[1])
2010  return data[indices[0]];
2011 
2012  // the rest is messier and requires a few switches.
2013  switch (dim)
2014  {
2015  case 2:
2016  // at least for the 2x2 case it is reasonably simple
2017  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
2018  ((indices[0] == 0) && (indices[1] == 1)),
2019  ExcInternalError());
2020  return data[2];
2021 
2022  default:
2023  // to do the rest, sort our indices before comparing
2024  {
2025  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2026  std::max(indices[0], indices[1]));
2027  for (unsigned int d = 0, c = 0; d < dim; ++d)
2028  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2029  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2030  return data[dim + c];
2031  Assert(false, ExcInternalError());
2032  }
2033  }
2034 
2035  // The code should never reach here.
2036  // We cannot return a static variable, as this class must support number
2037  // types that require no instances of the number type to be in scope during
2038  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2039  return data[0];
2040  }
2041 
2042 
2043 
2044  template <int dim, typename Number>
2045  constexpr inline Number &
2046  symmetric_tensor_access(const TableIndices<4> &indices,
2047  typename SymmetricTensorAccessors::
2048  StorageType<4, dim, Number>::base_tensor_type &data)
2049  {
2050  switch (dim)
2051  {
2052  case 1:
2053  return data[0][0];
2054 
2055  case 2:
2056  // each entry of the tensor can be thought of as an entry in a
2057  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2058  // rank-2 tensors. this is the format in which we store rank-4
2059  // tensors. determine which position the present entry is
2060  // stored in
2061  {
2062  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2063  return data[base_index[indices[0]][indices[1]]]
2064  [base_index[indices[2]][indices[3]]];
2065  }
2066  case 3:
2067  // each entry of the tensor can be thought of as an entry in a
2068  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2069  // rank-2 tensors. this is the format in which we store rank-4
2070  // tensors. determine which position the present entry is
2071  // stored in
2072  {
2073  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2074  {3, 1, 5},
2075  {4, 5, 2}};
2076  return data[base_index[indices[0]][indices[1]]]
2077  [base_index[indices[2]][indices[3]]];
2078  }
2079 
2080  default:
2081  Assert(false, ExcNotImplemented());
2082  }
2083 
2084  // The code should never reach here.
2085  // We cannot return a static variable, as this class must support number
2086  // types that require no instances of the number type to be in scope during
2087  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2088  return data[0][0];
2089  }
2090 
2091 
2092  template <int dim, typename Number>
2093  constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2094  symmetric_tensor_access(const TableIndices<4> &indices,
2095  const typename SymmetricTensorAccessors::
2096  StorageType<4, dim, Number>::base_tensor_type &data)
2097  {
2098  switch (dim)
2099  {
2100  case 1:
2101  return data[0][0];
2102 
2103  case 2:
2104  // each entry of the tensor can be thought of as an entry in a
2105  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2106  // rank-2 tensors. this is the format in which we store rank-4
2107  // tensors. determine which position the present entry is
2108  // stored in
2109  {
2110  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2111  return data[base_index[indices[0]][indices[1]]]
2112  [base_index[indices[2]][indices[3]]];
2113  }
2114  case 3:
2115  // each entry of the tensor can be thought of as an entry in a
2116  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2117  // rank-2 tensors. this is the format in which we store rank-4
2118  // tensors. determine which position the present entry is
2119  // stored in
2120  {
2121  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2122  {3, 1, 5},
2123  {4, 5, 2}};
2124  return data[base_index[indices[0]][indices[1]]]
2125  [base_index[indices[2]][indices[3]]];
2126  }
2127 
2128  default:
2129  Assert(false, ExcNotImplemented());
2130  }
2131 
2132  // The code should never reach here.
2133  // We cannot return a static variable, as this class must support number
2134  // types that require no instances of the number type to be in scope during
2135  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2136  return data[0][0];
2137  }
2138 
2139 } // end of namespace internal
2140 
2141 
2142 
2143 template <int rank_, int dim, typename Number>
2144 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2146  const TableIndices<rank_> &indices)
2147 {
2148  for (unsigned int r = 0; r < rank; ++r)
2149  AssertIndexRange(indices[r], dimension);
2150  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2151 }
2152 
2153 
2154 
2155 template <int rank_, int dim, typename Number>
2156 constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2158  const TableIndices<rank_> &indices) const
2159 {
2160  for (unsigned int r = 0; r < rank; ++r)
2161  AssertIndexRange(indices[r], dimension);
2162  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2163 }
2164 
2165 
2166 
2167 namespace internal
2168 {
2169  namespace SymmetricTensorImplementation
2170  {
2171  template <int rank_>
2172  constexpr TableIndices<rank_>
2173  get_partially_filled_indices(const unsigned int row,
2174  const std::integral_constant<int, 2> &)
2175  {
2177  }
2178 
2179 
2180  template <int rank_>
2181  constexpr TableIndices<rank_>
2182  get_partially_filled_indices(const unsigned int row,
2183  const std::integral_constant<int, 4> &)
2184  {
2185  return TableIndices<rank_>(row,
2189  }
2190  } // namespace SymmetricTensorImplementation
2191 } // namespace internal
2192 
2193 
2194 template <int rank_, int dim, typename Number>
2195 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2196  Accessor<rank_, dim, true, rank_ - 1, Number>
2197  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2198 {
2199  return internal::SymmetricTensorAccessors::
2200  Accessor<rank_, dim, true, rank_ - 1, Number>(
2201  *this,
2202  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2203  rank_>(row, std::integral_constant<int, rank_>()));
2204 }
2205 
2206 
2207 
2208 template <int rank_, int dim, typename Number>
2209 constexpr inline DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2210  Accessor<rank_, dim, false, rank_ - 1, Number>
2211  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2212 {
2213  return internal::SymmetricTensorAccessors::
2214  Accessor<rank_, dim, false, rank_ - 1, Number>(
2215  *this,
2216  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2217  rank_>(row, std::integral_constant<int, rank_>()));
2218 }
2219 
2220 
2221 
2222 template <int rank_, int dim, typename Number>
2223 constexpr DEAL_II_ALWAYS_INLINE const Number &
2225  const TableIndices<rank_> &indices) const
2226 {
2227  return operator()(indices);
2228 }
2229 
2230 
2231 
2232 template <int rank_, int dim, typename Number>
2233 constexpr inline DEAL_II_ALWAYS_INLINE Number &
2235  const TableIndices<rank_> &indices)
2236 {
2237  return operator()(indices);
2238 }
2239 
2240 
2241 
2242 template <int rank_, int dim, typename Number>
2243 inline Number *
2245 {
2246  return std::addressof(this->access_raw_entry(0));
2247 }
2248 
2249 
2250 
2251 template <int rank_, int dim, typename Number>
2252 inline const Number *
2254 {
2255  return std::addressof(this->access_raw_entry(0));
2256 }
2257 
2258 
2259 
2260 template <int rank_, int dim, typename Number>
2261 inline Number *
2263 {
2264  return begin_raw() + n_independent_components;
2265 }
2266 
2267 
2268 
2269 template <int rank_, int dim, typename Number>
2270 inline const Number *
2272 {
2273  return begin_raw() + n_independent_components;
2274 }
2275 
2276 
2277 
2278 namespace internal
2279 {
2280  namespace SymmetricTensorImplementation
2281  {
2282  template <int dim, typename Number>
2283  constexpr unsigned int
2284  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2285  const unsigned int index)
2286  {
2287  return index;
2288  }
2289 
2290 
2291  template <int dim, typename Number>
2292  constexpr ::TableIndices<2>
2293  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2294  const unsigned int index)
2295  {
2298  }
2299 
2300  } // namespace SymmetricTensorImplementation
2301 } // namespace internal
2302 
2303 
2304 
2305 template <int rank_, int dim, typename Number>
2306 constexpr inline const Number &
2308  const unsigned int index) const
2309 {
2310  AssertIndexRange(index, n_independent_components);
2311  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2312  index)];
2313 }
2314 
2315 
2316 
2317 template <int rank_, int dim, typename Number>
2318 constexpr inline Number &
2320 {
2321  AssertIndexRange(index, n_independent_components);
2322  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2323  index)];
2324 }
2325 
2326 
2327 
2328 namespace internal
2329 {
2330  template <int dim, typename Number>
2331  constexpr inline typename numbers::NumberTraits<Number>::real_type
2332  compute_norm(const typename SymmetricTensorAccessors::
2333  StorageType<2, dim, Number>::base_tensor_type &data)
2334  {
2335  // Make things work with AD types
2336  using std::sqrt;
2337  switch (dim)
2338  {
2339  case 1:
2340  return numbers::NumberTraits<Number>::abs(data[0]);
2341 
2342  case 2:
2346 
2347  case 3:
2354 
2355  default:
2356  {
2357  typename numbers::NumberTraits<Number>::real_type return_value =
2359 
2360  for (unsigned int d = 0; d < dim; ++d)
2361  return_value +=
2363  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2364  return_value +=
2366 
2367  return sqrt(return_value);
2368  }
2369  }
2370  }
2371 
2372 
2373 
2374  template <int dim, typename Number>
2375  constexpr inline typename numbers::NumberTraits<Number>::real_type
2376  compute_norm(const typename SymmetricTensorAccessors::
2377  StorageType<4, dim, Number>::base_tensor_type &data)
2378  {
2379  // Make things work with AD types
2380  using std::sqrt;
2381  switch (dim)
2382  {
2383  case 1:
2384  return numbers::NumberTraits<Number>::abs(data[0][0]);
2385 
2386  default:
2387  {
2388  typename numbers::NumberTraits<Number>::real_type return_value =
2390 
2391  const unsigned int n_independent_components = data.dimension;
2392 
2393  for (unsigned int i = 0; i < dim; ++i)
2394  for (unsigned int j = 0; j < dim; ++j)
2395  return_value +=
2397  for (unsigned int i = 0; i < dim; ++i)
2398  for (unsigned int j = dim; j < n_independent_components; ++j)
2399  return_value +=
2401  for (unsigned int i = dim; i < n_independent_components; ++i)
2402  for (unsigned int j = 0; j < dim; ++j)
2403  return_value +=
2405  for (unsigned int i = dim; i < n_independent_components; ++i)
2406  for (unsigned int j = dim; j < n_independent_components; ++j)
2407  return_value +=
2409 
2410  return sqrt(return_value);
2411  }
2412  }
2413  }
2414 
2415 } // end of namespace internal
2416 
2417 
2418 
2419 template <int rank_, int dim, typename Number>
2422 {
2423  return internal::compute_norm<dim, Number>(data);
2424 }
2425 
2426 
2427 
2428 namespace internal
2429 {
2430  namespace SymmetricTensorImplementation
2431  {
2432  // a function to do the unrolling from a set of indices to a
2433  // scalar index into the array in which we store the elements of
2434  // a symmetric tensor
2435  //
2436  // this function is for rank-2 tensors
2437  template <int dim>
2438  constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
2440  {
2441  AssertIndexRange(indices[0], dim);
2442  AssertIndexRange(indices[1], dim);
2443 
2444  switch (dim)
2445  {
2446  case 1:
2447  {
2448  return 0;
2449  }
2450 
2451  case 2:
2452  {
2453  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2454  return table[indices[0]][indices[1]];
2455  }
2456 
2457  case 3:
2458  {
2459  constexpr unsigned int table[3][3] = {{0, 3, 4},
2460  {3, 1, 5},
2461  {4, 5, 2}};
2462  return table[indices[0]][indices[1]];
2463  }
2464 
2465  case 4:
2466  {
2467  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2468  {4, 1, 7, 8},
2469  {5, 7, 2, 9},
2470  {6, 8, 9, 3}};
2471  return table[indices[0]][indices[1]];
2472  }
2473 
2474  default:
2475  // for the remainder, manually figure out the numbering
2476  {
2477  if (indices[0] == indices[1])
2478  return indices[0];
2479 
2480  TableIndices<2> sorted_indices(indices);
2481  sorted_indices.sort();
2482 
2483  for (unsigned int d = 0, c = 0; d < dim; ++d)
2484  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2485  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2486  return dim + c;
2487 
2488  // should never get here:
2489  Assert(false, ExcInternalError());
2490  return 0;
2491  }
2492  }
2493  }
2494 
2495  // a function to do the unrolling from a set of indices to a
2496  // scalar index into the array in which we store the elements of
2497  // a symmetric tensor
2498  //
2499  // this function is for tensors of ranks not already handled
2500  // above
2501  template <int dim, int rank_>
2502  constexpr inline unsigned int
2504  {
2505  (void)indices;
2506  Assert(false, ExcNotImplemented());
2508  }
2509  } // namespace SymmetricTensorImplementation
2510 } // namespace internal
2511 
2512 
2513 template <int rank_, int dim, typename Number>
2514 constexpr unsigned int
2516  const TableIndices<rank_> &indices)
2517 {
2518  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2519  dim>(indices);
2520 }
2521 
2522 
2523 
2524 namespace internal
2525 {
2526  namespace SymmetricTensorImplementation
2527  {
2528  // a function to do the inverse of the unrolling from a set of
2529  // indices to a scalar index into the array in which we store
2530  // the elements of a symmetric tensor. in other words, it goes
2531  // from the scalar index into the array to a set of indices of
2532  // the tensor
2533  //
2534  // this function is for rank-2 tensors
2535  template <int dim>
2536  constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2537  unrolled_to_component_indices(const unsigned int i,
2538  const std::integral_constant<int, 2> &)
2539  {
2540  Assert(
2542  ExcIndexRange(
2543  i,
2544  0,
2546  switch (dim)
2547  {
2548  case 1:
2549  {
2550  return {0, 0};
2551  }
2552 
2553  case 2:
2554  {
2555  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2556  TableIndices<2>(1, 1),
2557  TableIndices<2>(0, 1)};
2558  return table[i];
2559  }
2560 
2561  case 3:
2562  {
2563  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2564  TableIndices<2>(1, 1),
2565  TableIndices<2>(2, 2),
2566  TableIndices<2>(0, 1),
2567  TableIndices<2>(0, 2),
2568  TableIndices<2>(1, 2)};
2569  return table[i];
2570  }
2571 
2572  default:
2573  if (i < dim)
2574  return {i, i};
2575 
2576  for (unsigned int d = 0, c = dim; d < dim; ++d)
2577  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2578  if (c == i)
2579  return {d, e};
2580 
2581  // should never get here:
2582  Assert(false, ExcInternalError());
2583  return {0, 0};
2584  }
2585  }
2586 
2587  // a function to do the inverse of the unrolling from a set of
2588  // indices to a scalar index into the array in which we store
2589  // the elements of a symmetric tensor. in other words, it goes
2590  // from the scalar index into the array to a set of indices of
2591  // the tensor
2592  //
2593  // this function is for tensors of a rank not already handled
2594  // above
2595  template <int dim, int rank_>
2596  constexpr inline std::enable_if_t<rank_ != 2, TableIndices<rank_>>
2597  unrolled_to_component_indices(const unsigned int i,
2598  const std::integral_constant<int, rank_> &)
2599  {
2600  (void)i;
2601  Assert(
2602  (i <
2604  ExcIndexRange(i,
2605  0,
2607  n_independent_components));
2608  Assert(false, ExcNotImplemented());
2609  return TableIndices<rank_>();
2610  }
2611 
2612  } // namespace SymmetricTensorImplementation
2613 } // namespace internal
2614 
2615 template <int rank_, int dim, typename Number>
2618  const unsigned int i)
2619 {
2620  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2621  dim>(i, std::integral_constant<int, rank_>());
2622 }
2623 
2624 
2625 
2626 template <int rank_, int dim, typename Number>
2627 template <class Archive>
2628 inline void
2629 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2630 {
2631  ar &data;
2632 }
2633 
2634 
2635 #endif // DOXYGEN
2636 
2637 /* ----------------- Non-member functions operating on tensors. ------------ */
2638 
2639 
2652 template <int rank_, int dim, typename Number, typename OtherNumber>
2653 constexpr inline DEAL_II_ALWAYS_INLINE
2657 {
2659  tmp = left;
2660  tmp += right;
2661  return tmp;
2662 }
2663 
2664 
2677 template <int rank_, int dim, typename Number, typename OtherNumber>
2678 constexpr inline DEAL_II_ALWAYS_INLINE
2682 {
2684  tmp = left;
2685  tmp -= right;
2686  return tmp;
2687 }
2688 
2689 
2697 template <int rank_, int dim, typename Number, typename OtherNumber>
2698 constexpr DEAL_II_ALWAYS_INLINE
2701  const Tensor<rank_, dim, OtherNumber> & right)
2702 {
2703  return Tensor<rank_, dim, Number>(left) + right;
2704 }
2705 
2706 
2714 template <int rank_, int dim, typename Number, typename OtherNumber>
2715 constexpr DEAL_II_ALWAYS_INLINE
2719 {
2720  return left + Tensor<rank_, dim, OtherNumber>(right);
2721 }
2722 
2723 
2731 template <int rank_, int dim, typename Number, typename OtherNumber>
2732 constexpr DEAL_II_ALWAYS_INLINE
2735  const Tensor<rank_, dim, OtherNumber> & right)
2736 {
2737  return Tensor<rank_, dim, Number>(left) - right;
2738 }
2739 
2740 
2748 template <int rank_, int dim, typename Number, typename OtherNumber>
2749 constexpr DEAL_II_ALWAYS_INLINE
2753 {
2754  return left - Tensor<rank_, dim, OtherNumber>(right);
2755 }
2756 
2757 
2758 
2759 template <int dim, typename Number>
2762 {
2763  switch (dim)
2764  {
2765  case 1:
2766  return t.data[0];
2767  case 2:
2768  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2769  case 3:
2770  {
2771  // in analogy to general tensors, but
2772  // there's something to be simplified for
2773  // the present case
2774  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2775  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2776  t.data[0] * t.data[5] * t.data[5] -
2777  t.data[1] * t.data[4] * t.data[4] -
2778  t.data[2] * t.data[3] * t.data[3]);
2779  }
2780  default:
2781  Assert(false, ExcNotImplemented());
2783  }
2784 }
2785 
2786 
2787 
2799 template <int dim, typename Number>
2802 {
2803  return determinant(t);
2804 }
2805 
2806 
2807 
2808 template <int dim, typename Number>
2809 constexpr inline DEAL_II_ALWAYS_INLINE Number
2811 {
2812  Number t = d.data[0];
2813  for (unsigned int i = 1; i < dim; ++i)
2814  t += d.data[i];
2815  return t;
2816 }
2817 
2818 
2830 template <int dim, typename Number>
2831 constexpr Number
2833 {
2834  return trace(t);
2835 }
2836 
2837 
2849 template <typename Number>
2850 constexpr DEAL_II_ALWAYS_INLINE Number
2852 {
2854 }
2855 
2856 
2857 
2876 template <typename Number>
2877 constexpr DEAL_II_ALWAYS_INLINE Number
2879 {
2880  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2881 }
2882 
2883 
2884 
2893 template <typename Number>
2894 constexpr DEAL_II_ALWAYS_INLINE Number
2896 {
2897  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2898  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2899 }
2900 
2901 
2902 
2910 template <typename Number>
2911 std::array<Number, 1>
2913 
2914 
2915 
2938 template <typename Number>
2939 std::array<Number, 2>
2941 
2942 
2943 
2966 template <typename Number>
2967 std::array<Number, 3>
2969 
2970 
2971 
2972 namespace internal
2973 {
2974  namespace SymmetricTensorImplementation
2975  {
3013  template <int dim, typename Number>
3014  void
3015  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
3016  ::Tensor<2, dim, Number> & Q,
3017  std::array<Number, dim> & d,
3018  std::array<Number, dim - 1> & e);
3019 
3020 
3021 
3061  template <int dim, typename Number>
3062  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3063  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
3064 
3065 
3066 
3106  template <int dim, typename Number>
3107  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3109 
3110 
3111 
3125  template <typename Number>
3126  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3127  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3128 
3129 
3130 
3163  template <typename Number>
3164  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3165  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3166 
3171  template <int dim, typename Number>
3173  {
3174  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3175  bool
3176  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3177  {
3178  return lhs.first > rhs.first;
3179  }
3180  };
3181 
3182  } // namespace SymmetricTensorImplementation
3183 
3184 } // namespace internal
3185 
3186 
3187 
3188 // The line below is to ensure that doxygen puts the full description
3189 // of this global enumeration into the documentation
3190 // See https://stackoverflow.com/a/1717984
3220 {
3230  hybrid,
3248  jacobi
3249 };
3250 
3251 
3252 
3281 template <int dim, typename Number>
3282 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3283  std::integral_constant<int, dim>::value>
3285  const SymmetricTensorEigenvectorMethod method =
3287 
3288 
3289 
3298 template <int rank_, int dim, typename Number>
3301 {
3302  return t;
3303 }
3304 
3305 
3306 
3307 template <int dim, typename Number>
3310 {
3312 
3313  // subtract scaled trace from the diagonal
3314  const Number tr = trace(t) / dim;
3315  for (unsigned int i = 0; i < dim; ++i)
3316  tmp.data[i] -= tr;
3317 
3318  return tmp;
3319 }
3320 
3321 
3322 
3323 template <int dim, typename Number>
3326 {
3327  // create a default constructed matrix filled with
3328  // zeros, then set the diagonal elements to one
3330  switch (dim)
3331  {
3332  case 1:
3334  break;
3335  case 2:
3336  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3337  break;
3338  case 3:
3339  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3341  break;
3342  default:
3343  for (unsigned int d = 0; d < dim; ++d)
3345  }
3346  return tmp;
3347 }
3348 
3349 
3350 
3351 template <int dim, typename Number>
3354 {
3356 
3357  // fill the elements treating the diagonal
3358  for (unsigned int i = 0; i < dim; ++i)
3359  for (unsigned int j = 0; j < dim; ++j)
3360  tmp.data[i][j] =
3361  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3362 
3363  // then fill the ones that copy over the
3364  // non-diagonal elements. note that during
3365  // the double-contraction, we handle the
3366  // off-diagonal elements twice, so simply
3367  // copying requires a weight of 1/2
3368  for (unsigned int i = dim;
3369  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3370  n_rank2_components;
3371  ++i)
3372  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3373 
3374  return tmp;
3375 }
3376 
3377 
3378 
3379 template <int dim, typename Number>
3382 {
3384 
3385  // fill the elements treating the diagonal
3386  for (unsigned int i = 0; i < dim; ++i)
3387  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3388 
3389  // then fill the ones that copy over the
3390  // non-diagonal elements. note that during
3391  // the double-contraction, we handle the
3392  // off-diagonal elements twice, so simply
3393  // copying requires a weight of 1/2
3394  for (unsigned int i = dim;
3395  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3396  n_rank2_components;
3397  ++i)
3398  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3399 
3400  return tmp;
3401 }
3402 
3403 
3404 
3414 template <int dim, typename Number>
3417 {
3419  value(t);
3420 }
3421 
3422 
3423 
3434 template <int dim, typename Number>
3437 {
3439  value(t);
3440 }
3441 
3442 
3443 
3465 template <int dim, typename Number>
3466 constexpr inline SymmetricTensor<4, dim, Number>
3469 {
3471 
3472  // fill only the elements really needed
3473  for (unsigned int i = 0; i < dim; ++i)
3474  for (unsigned int j = i; j < dim; ++j)
3475  for (unsigned int k = 0; k < dim; ++k)
3476  for (unsigned int l = k; l < dim; ++l)
3477  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3478 
3479  return tmp;
3480 }
3481 
3482 
3483 
3491 template <int dim, typename Number>
3494 {
3496  for (unsigned int d = 0; d < dim; ++d)
3497  result[d][d] = t[d][d];
3498 
3499  const Number half = internal::NumberType<Number>::value(0.5);
3500  for (unsigned int d = 0; d < dim; ++d)
3501  for (unsigned int e = d + 1; e < dim; ++e)
3502  result[d][e] = (t[d][e] + t[e][d]) * half;
3503  return result;
3504 }
3505 
3506 
3507 
3515 template <int rank_, int dim, typename Number>
3517 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3518 {
3520  tt *= factor;
3521  return tt;
3522 }
3523 
3524 
3525 
3533 template <int rank_, int dim, typename Number>
3535 operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3536 {
3537  // simply forward to the other operator
3538  return t * factor;
3539 }
3540 
3541 
3542 
3567 template <int rank_, int dim, typename Number, typename OtherNumber>
3568 constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3569  rank_,
3570  dim,
3571  typename ProductType<Number,
3572  typename EnableIfScalar<OtherNumber>::type>::type>
3574  const OtherNumber & factor)
3575 {
3576  // form the product. we have to convert the two factors into the final
3577  // type via explicit casts because, for awkward reasons, the C++
3578  // standard committee saw it fit to not define an
3579  // operator*(float,std::complex<double>)
3580  // (as well as with switched arguments and double<->float).
3581  using product_type = typename ProductType<Number, OtherNumber>::type;
3584  return tt;
3585 }
3586 
3587 
3588 
3596 template <int rank_, int dim, typename Number, typename OtherNumber>
3597 constexpr inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3598  rank_,
3599  dim,
3600  typename ProductType<OtherNumber,
3601  typename EnableIfScalar<Number>::type>::type>
3602 operator*(const Number & factor,
3604 {
3605  // simply forward to the other operator with switched arguments
3606  return (t * factor);
3607 }
3608 
3609 
3610 
3616 template <int rank_, int dim, typename Number, typename OtherNumber>
3617 constexpr inline SymmetricTensor<
3618  rank_,
3619  dim,
3620  typename ProductType<Number,
3621  typename EnableIfScalar<OtherNumber>::type>::type>
3623  const OtherNumber & factor)
3624 {
3625  using product_type = typename ProductType<Number, OtherNumber>::type;
3628  return tt;
3629 }
3630 
3631 
3632 
3639 template <int rank_, int dim>
3641 operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3642 {
3644  tt *= factor;
3645  return tt;
3646 }
3647 
3648 
3649 
3656 template <int rank_, int dim>
3658 operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3659 {
3661  tt *= factor;
3662  return tt;
3663 }
3664 
3665 
3666 
3672 template <int rank_, int dim>
3673 constexpr inline SymmetricTensor<rank_, dim>
3674 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3675 {
3677  tt /= factor;
3678  return tt;
3679 }
3680 
3690 template <int dim, typename Number, typename OtherNumber>
3694 {
3695  return (t1 * t2);
3696 }
3697 
3698 
3711 template <int dim, typename Number, typename OtherNumber>
3712 constexpr inline DEAL_II_ALWAYS_INLINE
3715  const Tensor<2, dim, OtherNumber> & t2)
3716 {
3718  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3719  for (unsigned int i = 0; i < dim; ++i)
3720  for (unsigned int j = 0; j < dim; ++j)
3721  s += t1[i][j] * t2[i][j];
3722  return s;
3723 }
3724 
3725 
3738 template <int dim, typename Number, typename OtherNumber>
3742 {
3743  return scalar_product(t2, t1);
3744 }
3745 
3746 
3761 template <typename Number, typename OtherNumber>
3762 constexpr inline DEAL_II_ALWAYS_INLINE void
3767 {
3768  tmp[0][0] = t[0][0][0][0] * s[0][0];
3769 }
3770 
3771 
3772 
3787 template <typename Number, typename OtherNumber>
3788 constexpr inline void
3793 {
3794  tmp[0][0] = t[0][0][0][0] * s[0][0];
3795 }
3796 
3797 
3798 
3813 template <typename Number, typename OtherNumber>
3814 constexpr inline void
3819 {
3820  const unsigned int dim = 2;
3821 
3822  for (unsigned int i = 0; i < dim; ++i)
3823  for (unsigned int j = i; j < dim; ++j)
3824  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3825  2 * t[i][j][0][1] * s[0][1];
3826 }
3827 
3828 
3829 
3844 template <typename Number, typename OtherNumber>
3845 constexpr inline void
3850 {
3851  const unsigned int dim = 2;
3852 
3853  for (unsigned int i = 0; i < dim; ++i)
3854  for (unsigned int j = i; j < dim; ++j)
3855  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3856  2 * s[0][1] * t[0][1][i][j];
3857 }
3858 
3859 
3860 
3875 template <typename Number, typename OtherNumber>
3876 constexpr inline void
3881 {
3882  const unsigned int dim = 3;
3883 
3884  for (unsigned int i = 0; i < dim; ++i)
3885  for (unsigned int j = i; j < dim; ++j)
3886  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3887  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3888  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3889 }
3890 
3891 
3892 
3907 template <typename Number, typename OtherNumber>
3908 constexpr inline void
3913 {
3914  const unsigned int dim = 3;
3915 
3916  for (unsigned int i = 0; i < dim; ++i)
3917  for (unsigned int j = i; j < dim; ++j)
3918  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3919  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3920  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3921 }
3922 
3923 
3924 
3931 template <int dim, typename Number, typename OtherNumber>
3934  const Tensor<1, dim, OtherNumber> & src2)
3935 {
3937  for (unsigned int i = 0; i < dim; ++i)
3938  for (unsigned int j = 0; j < dim; ++j)
3939  dest[i] += src1[i][j] * src2[j];
3940  return dest;
3941 }
3942 
3943 
3950 template <int dim, typename Number, typename OtherNumber>
3954 {
3955  // this is easy for symmetric tensors:
3956  return src2 * src1;
3957 }
3958 
3959 
3960 
3980 template <int rank_1,
3981  int rank_2,
3982  int dim,
3983  typename Number,
3984  typename OtherNumber>
3985 constexpr DEAL_II_ALWAYS_INLINE
3986  typename Tensor<rank_1 + rank_2 - 2,
3987  dim,
3988  typename ProductType<Number, OtherNumber>::type>::tensor_type
3991 {
3992  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
3993 }
3994 
3995 
3996 
4016 template <int rank_1,
4017  int rank_2,
4018  int dim,
4019  typename Number,
4020  typename OtherNumber>
4021 constexpr DEAL_II_ALWAYS_INLINE
4022  typename Tensor<rank_1 + rank_2 - 2,
4023  dim,
4024  typename ProductType<Number, OtherNumber>::type>::tensor_type
4026  const Tensor<rank_2, dim, OtherNumber> & src2)
4027 {
4028  return Tensor<rank_1, dim, Number>(src1) * src2;
4029 }
4030 
4031 
4032 
4042 template <int dim, typename Number>
4043 inline std::ostream &
4044 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4045 {
4046  // make our lives a bit simpler by outputting
4047  // the tensor through the operator for the
4048  // general Tensor class
4050 
4051  for (unsigned int i = 0; i < dim; ++i)
4052  for (unsigned int j = 0; j < dim; ++j)
4053  tt[i][j] = t[i][j];
4054 
4055  return out << tt;
4056 }
4057 
4058 
4059 
4069 template <int dim, typename Number>
4070 inline std::ostream &
4071 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4072 {
4073  // make our lives a bit simpler by outputting
4074  // the tensor through the operator for the
4075  // general Tensor class
4077 
4078  for (unsigned int i = 0; i < dim; ++i)
4079  for (unsigned int j = 0; j < dim; ++j)
4080  for (unsigned int k = 0; k < dim; ++k)
4081  for (unsigned int l = 0; l < dim; ++l)
4082  tt[i][j][k][l] = t[i][j][k][l];
4083 
4084  return out << tt;
4085 }
4086 
4087 
4089 
4090 #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:108
#define DEAL_II_DEPRECATED
Definition: config.h:162
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_CONSTEXPR
Definition: config.h:175
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1509
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1768
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< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
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()