Reference documentation for deal.II version GIT 6a72d26406 2023-06-07 13:05:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
symmetric_tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 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>
51 
80 template <int dim, typename Number = double>
84 
122 template <int dim, typename Number = double>
126 
127 template <int dim, typename Number>
130 
131 template <int dim, typename Number>
134 
144 template <int dim2, typename Number>
145 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
147 
158 template <int dim, typename Number>
159 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
162 
176 template <int dim, typename Number>
179 
180 
181 
182 namespace internal
183 {
184  // Workaround: The following 4 overloads are necessary to be able to
185  // compile the library with Apple Clang 8 and older. We should remove
186  // these overloads again when we bump the minimal required version to
187  // something later than clang-3.6 / Apple Clang 6.3.
188  template <int rank, int dim, typename T, typename U>
189  struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
190  {
191  using type =
192  SymmetricTensor<rank,
193  dim,
194  std::complex<typename ProductType<T, U>::type>>;
195  };
196 
197  template <int rank, int dim, typename T, typename U>
198  struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
199  std::complex<U>>
200  {
201  using type =
202  SymmetricTensor<rank,
203  dim,
204  std::complex<typename ProductType<T, U>::type>>;
205  };
206 
207  template <typename T, int rank, int dim, typename U>
208  struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
209  {
210  using type =
211  SymmetricTensor<rank,
212  dim,
213  std::complex<typename ProductType<T, U>::type>>;
214  };
215 
216  template <int rank, int dim, typename T, typename U>
217  struct ProductTypeImpl<std::complex<T>,
218  SymmetricTensor<rank, dim, std::complex<U>>>
219  {
220  using type =
221  SymmetricTensor<rank,
222  dim,
223  std::complex<typename ProductType<T, U>::type>>;
224  };
225  // end workaround
226 
231  namespace SymmetricTensorImplementation
232  {
237  template <int rank, int dim, typename Number>
238  struct Inverse;
239  } // namespace SymmetricTensorImplementation
240 
245  namespace SymmetricTensorAccessors
246  {
255  merge(const TableIndices<2> &previous_indices,
256  const unsigned int new_index,
257  const unsigned int position)
258  {
259  AssertIndexRange(position, 2);
260 
261  if (position == 0)
262  return {new_index, numbers::invalid_unsigned_int};
263  else
264  return {previous_indices[0], new_index};
265  }
266 
267 
268 
277  merge(const TableIndices<4> &previous_indices,
278  const unsigned int new_index,
279  const unsigned int position)
280  {
281  AssertIndexRange(position, 4);
282 
283  switch (position)
284  {
285  case 0:
286  return {new_index,
290  case 1:
291  return {previous_indices[0],
292  new_index,
295  case 2:
296  return {previous_indices[0],
297  previous_indices[1],
298  new_index,
300  case 3:
301  return {previous_indices[0],
302  previous_indices[1],
303  previous_indices[2],
304  new_index};
305  default:
306  Assert(false, ExcInternalError());
307  return {};
308  }
309  }
310 
311 
318  template <int rank1,
319  int rank2,
320  int dim,
321  typename Number,
322  typename OtherNumber = Number>
324  {
326  using type =
327  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
328  };
329 
330 
337  template <int dim, typename Number, typename OtherNumber>
338  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
339  {
341  };
342 
343 
344 
357  template <int rank, int dim, typename Number>
358  struct StorageType;
359 
363  template <int dim, typename Number>
364  struct StorageType<2, dim, Number>
365  {
370  static const unsigned int n_independent_components =
371  (dim * dim + dim) / 2;
372 
377  };
378 
379 
380 
384  template <int dim, typename Number>
385  struct StorageType<4, dim, Number>
386  {
392  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
393 
397  static const unsigned int n_independent_components =
398  (n_rank2_components *
400 
408  };
409 
410 
411 
416  template <int rank, int dim, bool constness, typename Number>
418 
425  template <int rank, int dim, typename Number>
426  struct AccessorTypes<rank, dim, true, Number>
427  {
428  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
429 
430  using reference = const Number &;
431  };
432 
439  template <int rank, int dim, typename Number>
440  struct AccessorTypes<rank, dim, false, Number>
441  {
443 
444  using reference = Number &;
445  };
446 
447 
480  template <int rank, int dim, bool constness, int P, typename Number>
481  class Accessor
482  {
483  public:
487  using reference =
489  using tensor_type =
491 
492  private:
514 
519  constexpr DEAL_II_ALWAYS_INLINE
520  Accessor(const Accessor &) = default;
521 
522  public:
527  constexpr Accessor<rank, dim, constness, P - 1, Number>
528  operator[](const unsigned int i);
529 
534  constexpr Accessor<rank, dim, constness, P - 1, Number>
535  operator[](const unsigned int i) const;
536 
537  private:
543 
544  // Declare some other classes as friends. Make sure to work around bugs
545  // in some compilers:
546  template <int, int, typename>
547  friend class ::SymmetricTensor;
548  template <int, int, bool, int, typename>
549  friend class Accessor;
550  friend class ::SymmetricTensor<rank, dim, Number>;
551  friend class Accessor<rank, dim, constness, P + 1, Number>;
552  };
553 
554 
555 
563  template <int rank, int dim, bool constness, typename Number>
564  class Accessor<rank, dim, constness, 1, Number>
565  {
566  public:
570  using reference =
572  using tensor_type =
574 
575  private:
600 
605  constexpr DEAL_II_ALWAYS_INLINE
606  Accessor(const Accessor &) = default;
607 
608  public:
613  constexpr reference
614  operator[](const unsigned int);
615 
620  constexpr reference
621  operator[](const unsigned int) const;
622 
623  private:
629 
630  // Declare some other classes as friends. Make sure to work around bugs
631  // in some compilers:
632  template <int, int, typename>
633  friend class ::SymmetricTensor;
634  template <int, int, bool, int, typename>
636  friend class ::SymmetricTensor<rank, dim, Number>;
637  friend class SymmetricTensorAccessors::
638  Accessor<rank, dim, constness, 2, Number>;
639  };
640  } // namespace SymmetricTensorAccessors
641 } // namespace internal
642 
643 
644 
717 template <int rank_, int dim, typename Number>
719 {
720 public:
721  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
722 
731  static constexpr unsigned int dimension = dim;
732 
736  static const unsigned int rank = rank_;
737 
743  static constexpr unsigned int n_independent_components =
745  n_independent_components;
746 
751  constexpr DEAL_II_ALWAYS_INLINE
752  SymmetricTensor() = default;
753 
767  template <typename OtherNumber>
769 
786  constexpr SymmetricTensor(const Number (&array)[n_independent_components]);
787 
793  template <typename OtherNumber>
794  DEAL_II_HOST constexpr explicit SymmetricTensor(
795  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
796 
806  Number *
808 
818  const Number *
819  begin_raw() const;
820 
830  Number *
832 
843  const Number *
844  end_raw() const;
845 
852  template <typename OtherNumber>
853  DEAL_II_HOST constexpr SymmetricTensor &
855 
863  constexpr SymmetricTensor &
864  operator=(const Number &d);
865 
871  constexpr operator Tensor<rank_, dim, Number>() const;
872 
877  constexpr bool
878  operator==(const SymmetricTensor &) const;
879 
884  constexpr bool
885  operator!=(const SymmetricTensor &) const;
886 
890  template <typename OtherNumber>
891  DEAL_II_HOST constexpr SymmetricTensor &
893 
897  template <typename OtherNumber>
898  DEAL_II_HOST constexpr SymmetricTensor &
900 
905  template <typename OtherNumber>
906  DEAL_II_HOST constexpr SymmetricTensor &
907  operator*=(const OtherNumber &factor);
908 
912  template <typename OtherNumber>
913  DEAL_II_HOST constexpr SymmetricTensor &
914  operator/=(const OtherNumber &factor);
915 
920  constexpr SymmetricTensor
921  operator-() const;
922 
975  template <typename OtherNumber>
979 
984  template <typename OtherNumber>
988 
993  constexpr Number &
995 
1000  constexpr const Number &
1001  operator()(const TableIndices<rank_> &indices) const;
1002 
1007  DEAL_II_HOST
1008  constexpr internal::SymmetricTensorAccessors::
1009  Accessor<rank_, dim, true, rank_ - 1, Number>
1010  operator[](const unsigned int row) const;
1011 
1016  DEAL_II_HOST
1017  constexpr internal::SymmetricTensorAccessors::
1018  Accessor<rank_, dim, false, rank_ - 1, Number>
1019  operator[](const unsigned int row);
1020 
1026  DEAL_II_HOST
1027  constexpr const Number &
1028  operator[](const TableIndices<rank_> &indices) const;
1029 
1035  DEAL_II_HOST
1036  constexpr Number &
1038 
1045  DEAL_II_HOST
1046  constexpr const Number &
1047  access_raw_entry(const unsigned int unrolled_index) const;
1048 
1055  DEAL_II_HOST
1056  constexpr Number &
1057  access_raw_entry(const unsigned int unrolled_index);
1058 
1068  DEAL_II_HOST
1069  constexpr typename numbers::NumberTraits<Number>::real_type
1070  norm() const;
1071 
1079  static DEAL_II_HOST constexpr unsigned int
1081 
1087  static DEAL_II_HOST constexpr TableIndices<rank_>
1088  unrolled_to_component_indices(const unsigned int i);
1089 
1102  DEAL_II_HOST
1103  constexpr void
1105 
1110  static DEAL_II_HOST constexpr std::size_t
1112 
1118  template <class Archive>
1119  void
1120  serialize(Archive &ar, const unsigned int version);
1121 
1122 private:
1128 
1132  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
1133 
1138 
1139 #ifndef DOXYGEN
1140 
1141  // Make all other symmetric tensors friends.
1142  template <int, int, typename>
1143  friend class SymmetricTensor;
1144 
1145  // Make a few more functions friends.
1146  template <int dim2, typename Number2>
1147  friend DEAL_II_HOST constexpr Number2
1149 
1150  template <int dim2, typename Number2>
1151  friend DEAL_II_HOST DEAL_II_CONSTEXPR Number2
1153 
1154  template <int dim2, typename Number2>
1157 
1158  template <int dim2, typename Number2>
1161 
1162  template <int dim2, typename Number2>
1164  deviator_tensor();
1165 
1166  template <int dim2, typename Number2>
1168  identity_tensor();
1169 
1170 
1171  // Make a few helper classes friends as well.
1173  Inverse<2, dim, Number>;
1174 
1176  Inverse<4, dim, Number>;
1177 #endif
1178 };
1179 
1180 
1181 
1182 // ------------------------- inline functions ------------------------
1183 
1184 #ifndef DOXYGEN
1185 
1186 // provide declarations for static members
1187 template <int rank, int dim, typename Number>
1188 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1189 
1190 template <int rank_, int dim, typename Number>
1191 constexpr unsigned int
1192  SymmetricTensor<rank_, dim, Number>::n_independent_components;
1193 
1194 namespace internal
1195 {
1196  namespace SymmetricTensorAccessors
1197  {
1198  template <int rank_, int dim, bool constness, int P, typename Number>
1201  tensor_type & tensor,
1203  : tensor(tensor)
1205  {}
1206 
1207 
1208 
1209  template <int rank_, int dim, bool constness, int P, typename Number>
1210  DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1211  Accessor<rank_, dim, constness, P - 1, Number>
1212  Accessor<rank_, dim, constness, P, Number>::operator[](
1213  const unsigned int i)
1214  {
1215  return Accessor<rank_, dim, constness, P - 1, Number>(
1216  tensor, merge(previous_indices, i, rank_ - P));
1217  }
1218 
1219 
1220 
1221  template <int rank_, int dim, bool constness, int P, typename Number>
1223  Accessor<rank_, dim, constness, P - 1, Number>
1224  Accessor<rank_, dim, constness, P, Number>::operator[](
1225  const unsigned int i) const
1226  {
1227  return Accessor<rank_, dim, constness, P - 1, Number>(
1228  tensor, merge(previous_indices, i, rank_ - P));
1229  }
1230 
1231 
1232 
1233  template <int rank_, int dim, bool constness, typename Number>
1235  Accessor<rank_, dim, constness, 1, Number>::Accessor(
1236  tensor_type & tensor,
1237  const TableIndices<rank_> &previous_indices)
1238  : tensor(tensor)
1239  , previous_indices(previous_indices)
1240  {}
1241 
1242 
1243 
1244  template <int rank_, int dim, bool constness, typename Number>
1245  DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1246  typename Accessor<rank_, dim, constness, 1, Number>::reference
1247  Accessor<rank_, dim, constness, 1, Number>::operator[](
1248  const unsigned int i)
1249  {
1250  return tensor(merge(previous_indices, i, rank_ - 1));
1251  }
1252 
1253 
1254  template <int rank_, int dim, bool constness, typename Number>
1256  typename Accessor<rank_, dim, constness, 1, Number>::reference
1257  Accessor<rank_, dim, constness, 1, Number>::operator[](
1258  const unsigned int i) const
1259  {
1260  return tensor(merge(previous_indices, i, rank_ - 1));
1261  }
1262  } // namespace SymmetricTensorAccessors
1263 } // namespace internal
1264 
1265 
1266 
1267 template <int rank_, int dim, typename Number>
1268 template <typename OtherNumber>
1269 inline DEAL_II_ALWAYS_INLINE
1271  const Tensor<2, dim, OtherNumber> &t)
1272 {
1273  static_assert(rank == 2, "This function is only implemented for rank==2");
1274  for (unsigned int d = 0; d < dim; ++d)
1275  for (unsigned int e = 0; e < d; ++e)
1276  Assert(t[d][e] == t[e][d],
1277  ExcMessage("The incoming Tensor must be exactly symmetric."));
1278 
1279  for (unsigned int d = 0; d < dim; ++d)
1280  data[d] = t[d][d];
1281 
1282  for (unsigned int d = 0, c = 0; d < dim; ++d)
1283  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1284  data[dim + c] = t[d][e];
1285 }
1286 
1287 
1288 
1289 template <int rank_, int dim, typename Number>
1290 template <typename OtherNumber>
1293  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1294  : data(initializer.data)
1295 {}
1296 
1297 
1298 
1299 template <int rank_, int dim, typename Number>
1300 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1302  const Number (&array)[n_independent_components])
1303  : data(
1304  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1305 {
1306  // ensure that the reinterpret_cast above actually works
1307  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1308  ExcInternalError());
1309 }
1310 
1311 
1312 
1313 template <int rank_, int dim, typename Number>
1314 template <typename OtherNumber>
1315 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1319 {
1320  data = t.data;
1321  return *this;
1322 }
1323 
1324 
1325 
1326 template <int rank_, int dim, typename Number>
1327 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1330 {
1332  ExcMessage("Only assignment with zero is allowed"));
1333  (void)d;
1334 
1336 
1337  return *this;
1338 }
1339 
1340 
1341 namespace internal
1342 {
1343  namespace SymmetricTensorImplementation
1344  {
1345  template <int dim, typename Number>
1346  constexpr inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1347  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1348  {
1350 
1351  // diagonal entries are stored first
1352  for (unsigned int d = 0; d < dim; ++d)
1353  t[d][d] = s.access_raw_entry(d);
1354 
1355  // off-diagonal entries come next, row by row
1356  for (unsigned int d = 0, c = 0; d < dim; ++d)
1357  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1358  {
1359  t[d][e] = s.access_raw_entry(dim + c);
1360  t[e][d] = s.access_raw_entry(dim + c);
1361  }
1362  return t;
1363  }
1364 
1365 
1366  template <int dim, typename Number>
1367  constexpr ::Tensor<4, dim, Number>
1368  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1369  {
1370  // utilize the symmetry properties of SymmetricTensor<4,dim>
1371  // discussed in the class documentation to avoid accessing all
1372  // independent elements of the input tensor more than once
1374 
1375  for (unsigned int i = 0; i < dim; ++i)
1376  for (unsigned int j = i; j < dim; ++j)
1377  for (unsigned int k = 0; k < dim; ++k)
1378  for (unsigned int l = k; l < dim; ++l)
1379  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1380  t[TableIndices<4>(j, i, k, l)] =
1381  t[TableIndices<4>(j, i, l, k)] =
1382  st[TableIndices<4>(i, j, k, l)];
1383 
1384  return t;
1385  }
1386 
1387 
1388  template <typename Number>
1389  struct Inverse<2, 1, Number>
1390  {
1391  constexpr static inline DEAL_II_ALWAYS_INLINE
1392  ::SymmetricTensor<2, 1, Number>
1393  value(const ::SymmetricTensor<2, 1, Number> &t)
1394  {
1396 
1397  tmp[0][0] = 1.0 / t[0][0];
1398 
1399  return tmp;
1400  }
1401  };
1402 
1403 
1404  template <typename Number>
1405  struct Inverse<2, 2, Number>
1406  {
1407  constexpr static inline DEAL_II_ALWAYS_INLINE
1408  ::SymmetricTensor<2, 2, Number>
1409  value(const ::SymmetricTensor<2, 2, Number> &t)
1410  {
1412 
1413  // Sympy result: ([
1414  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1415  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1416  const TableIndices<2> idx_00(0, 0);
1417  const TableIndices<2> idx_01(0, 1);
1418  const TableIndices<2> idx_11(1, 1);
1419  const Number inv_det_t =
1420  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1421  tmp[idx_00] = t[idx_11];
1422  tmp[idx_01] = -t[idx_01];
1423  tmp[idx_11] = t[idx_00];
1424  tmp *= inv_det_t;
1425 
1426  return tmp;
1427  }
1428  };
1429 
1430 
1431  template <typename Number>
1432  struct Inverse<2, 3, Number>
1433  {
1434  constexpr static ::SymmetricTensor<2, 3, Number>
1435  value(const ::SymmetricTensor<2, 3, Number> &t)
1436  {
1438 
1439  // Sympy result: ([
1440  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1441  // 2*t01*t02*t12 - t02**2*t11),
1442  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1443  // 2*t01*t02*t12 - t02**2*t11),
1444  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1445  // 2*t01*t02*t12 - t02**2*t11)],
1446  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1447  // 2*t01*t02*t12 - t02**2*t11),
1448  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1449  // 2*t01*t02*t12 - t02**2*t11),
1450  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1451  // 2*t01*t02*t12 + t02**2*t11)],
1452  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1453  // 2*t01*t02*t12 - t02**2*t11),
1454  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1455  // 2*t01*t02*t12 + t02**2*t11),
1456  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1457  // 2*t01*t02*t12 + t02**2*t11)] ])
1458  //
1459  // =
1460  //
1461  // [ (t11*t22 - t12**2)/det_t,
1462  // (-t01*t22 + t02*t12)/det_t,
1463  // (t01*t12 - t02*t11)/det_t],
1464  // [ (-t01*t22 + t02*t12)/det_t,
1465  // (t00*t22 - t02**2)/det_t,
1466  // (-t00*t12 + t01*t02)/det_t],
1467  // [ (t01*t12 - t02*t11)/det_t,
1468  // (-t00*t12 + t01*t02)/det_t,
1469  // (t00*t11 - t01**2)/det_t] ])
1470  //
1471  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1472  // 2*t01*t02*t12 - t02**2*t11)
1473  const TableIndices<2> idx_00(0, 0);
1474  const TableIndices<2> idx_01(0, 1);
1475  const TableIndices<2> idx_02(0, 2);
1476  const TableIndices<2> idx_11(1, 1);
1477  const TableIndices<2> idx_12(1, 2);
1478  const TableIndices<2> idx_22(2, 2);
1479  const Number inv_det_t =
1480  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1481  t[idx_00] * t[idx_12] * t[idx_12] -
1482  t[idx_01] * t[idx_01] * t[idx_22] +
1483  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1484  t[idx_02] * t[idx_02] * t[idx_11]);
1485  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1486  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1487  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1488  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1489  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1490  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1491  tmp *= inv_det_t;
1492 
1493  return tmp;
1494  }
1495  };
1496 
1497 
1498  template <typename Number>
1499  struct Inverse<4, 1, Number>
1500  {
1501  constexpr static inline ::SymmetricTensor<4, 1, Number>
1502  value(const ::SymmetricTensor<4, 1, Number> &t)
1503  {
1505  tmp.data[0][0] = 1.0 / t.data[0][0];
1506  return tmp;
1507  }
1508  };
1509 
1510 
1511  template <typename Number>
1512  struct Inverse<4, 2, Number>
1513  {
1514  constexpr static inline ::SymmetricTensor<4, 2, Number>
1515  value(const ::SymmetricTensor<4, 2, Number> &t)
1516  {
1518 
1519  // Inverting this tensor is a little more complicated than necessary,
1520  // since we store the data of 't' as a 3x3 matrix t.data, but the
1521  // product between a rank-4 and a rank-2 tensor is really not the
1522  // product between this matrix and the 3-vector of a rhs, but rather
1523  //
1524  // B.vec = t.data * mult * A.vec
1525  //
1526  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1527  // capture the fact that we need to add up both the c_ij12*a_12 and the
1528  // c_ij21*a_21 terms.
1529  //
1530  // In addition, in this scheme, the identity tensor has the matrix
1531  // representation mult^-1.
1532  //
1533  // The inverse of 't' therefore has the matrix representation
1534  //
1535  // inv.data = mult^-1 * t.data^-1 * mult^-1
1536  //
1537  // in order to compute it, let's first compute the inverse of t.data and
1538  // put it into tmp.data; at the end of the function we then scale the
1539  // last row and column of the inverse by 1/2, corresponding to the left
1540  // and right multiplication with mult^-1.
1541  const Number t4 = t.data[0][0] * t.data[1][1],
1542  t6 = t.data[0][0] * t.data[1][2],
1543  t8 = t.data[0][1] * t.data[1][0],
1544  t00 = t.data[0][2] * t.data[1][0],
1545  t01 = t.data[0][1] * t.data[2][0],
1546  t04 = t.data[0][2] * t.data[2][0],
1547  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1548  t8 * t.data[2][2] + t00 * t.data[2][1] +
1549  t01 * t.data[1][2] - t04 * t.data[1][1]);
1550  tmp.data[0][0] =
1551  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1552  tmp.data[0][1] =
1553  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1554  tmp.data[0][2] =
1555  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1556  tmp.data[1][0] =
1557  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1558  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1559  tmp.data[1][2] = -(t6 - t00) * t07;
1560  tmp.data[2][0] =
1561  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1562  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1563  tmp.data[2][2] = (t4 - t8) * t07;
1564 
1565  // scale last row and column as mentioned
1566  // above
1567  tmp.data[2][0] /= 2;
1568  tmp.data[2][1] /= 2;
1569  tmp.data[0][2] /= 2;
1570  tmp.data[1][2] /= 2;
1571  tmp.data[2][2] /= 4;
1572 
1573  return tmp;
1574  }
1575  };
1576 
1577 
1578  template <typename Number>
1579  struct Inverse<4, 3, Number>
1580  {
1581  static ::SymmetricTensor<4, 3, Number>
1582  value(const ::SymmetricTensor<4, 3, Number> &t)
1583  {
1585 
1586  // This function follows the exact same scheme as the 2d case, except
1587  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1588  // Instead, we use the Gauss-Jordan algorithm implemented for
1589  // FullMatrix. For historical reasons the following code is copied from
1590  // there, with the tangential benefit that we do not need to copy the
1591  // tensor entries to and from the FullMatrix.
1592  const unsigned int N = 6;
1593 
1594  // First get an estimate of the size of the elements of this matrix,
1595  // for later checks whether the pivot element is large enough, or
1596  // whether we have to fear that the matrix is not regular.
1597  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1598  for (unsigned int i = 0; i < N; ++i)
1599  diagonal_sum += numbers::NumberTraits<Number>::abs(tmp.data[i][i]);
1600  const Number typical_diagonal_element =
1601  diagonal_sum / static_cast<double>(N);
1602  (void)typical_diagonal_element;
1603 
1604  unsigned int p[N];
1605  for (unsigned int i = 0; i < N; ++i)
1606  p[i] = i;
1607 
1608  for (unsigned int j = 0; j < N; ++j)
1609  {
1610  // Pivot search: search that part of the line on and right of the
1611  // diagonal for the largest element.
1612  Number max = numbers::NumberTraits<Number>::abs(tmp.data[j][j]);
1613  unsigned int r = j;
1614  for (unsigned int i = j + 1; i < N; ++i)
1615  if (numbers::NumberTraits<Number>::abs(tmp.data[i][j]) > max)
1616  {
1618  r = i;
1619  }
1620 
1621  // Check whether the pivot is too small
1622  Assert(max > 1.e-16 * typical_diagonal_element,
1623  ExcMessage("This tensor seems to be noninvertible"));
1624 
1625  // Row interchange
1626  if (r > j)
1627  {
1628  for (unsigned int k = 0; k < N; ++k)
1629  std::swap(tmp.data[j][k], tmp.data[r][k]);
1630 
1631  std::swap(p[j], p[r]);
1632  }
1633 
1634  // Transformation
1635  const Number hr = 1. / tmp.data[j][j];
1636  tmp.data[j][j] = hr;
1637  for (unsigned int k = 0; k < N; ++k)
1638  {
1639  if (k == j)
1640  continue;
1641  for (unsigned int i = 0; i < N; ++i)
1642  {
1643  if (i == j)
1644  continue;
1645  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1646  }
1647  }
1648  for (unsigned int i = 0; i < N; ++i)
1649  {
1650  tmp.data[i][j] *= hr;
1651  tmp.data[j][i] *= -hr;
1652  }
1653  tmp.data[j][j] = hr;
1654  }
1655 
1656  // Column interchange
1657  Number hv[N];
1658  for (unsigned int i = 0; i < N; ++i)
1659  {
1660  for (unsigned int k = 0; k < N; ++k)
1661  hv[p[k]] = tmp.data[i][k];
1662  for (unsigned int k = 0; k < N; ++k)
1663  tmp.data[i][k] = hv[k];
1664  }
1665 
1666  // Scale rows and columns. The mult matrix
1667  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1668  for (unsigned int i = 3; i < 6; ++i)
1669  for (unsigned int j = 0; j < 3; ++j)
1670  tmp.data[i][j] /= 2;
1671 
1672  for (unsigned int i = 0; i < 3; ++i)
1673  for (unsigned int j = 3; j < 6; ++j)
1674  tmp.data[i][j] /= 2;
1675 
1676  for (unsigned int i = 3; i < 6; ++i)
1677  for (unsigned int j = 3; j < 6; ++j)
1678  tmp.data[i][j] /= 4;
1679 
1680  return tmp;
1681  }
1682  };
1683 
1684  } // namespace SymmetricTensorImplementation
1685 } // namespace internal
1686 
1687 
1688 
1689 template <int rank_, int dim, typename Number>
1692  const
1693 {
1694  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1695 }
1696 
1697 
1698 
1699 template <int rank_, int dim, typename Number>
1700 DEAL_II_HOST constexpr bool
1702  const SymmetricTensor<rank_, dim, Number> &t) const
1703 {
1704  return data == t.data;
1705 }
1706 
1707 
1708 
1709 template <int rank_, int dim, typename Number>
1710 DEAL_II_HOST constexpr bool
1712  const SymmetricTensor<rank_, dim, Number> &t) const
1713 {
1714  return data != t.data;
1715 }
1716 
1717 
1718 
1719 template <int rank_, int dim, typename Number>
1720 template <typename OtherNumber>
1721 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1725 {
1726  data += t.data;
1727  return *this;
1728 }
1729 
1730 
1731 
1732 template <int rank_, int dim, typename Number>
1733 template <typename OtherNumber>
1734 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1738 {
1739  data -= t.data;
1740  return *this;
1741 }
1742 
1743 
1744 
1745 template <int rank_, int dim, typename Number>
1746 template <typename OtherNumber>
1747 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1750 {
1751  data *= d;
1752  return *this;
1753 }
1754 
1755 
1756 
1757 template <int rank_, int dim, typename Number>
1758 template <typename OtherNumber>
1759 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1762 {
1763  data /= d;
1764  return *this;
1765 }
1766 
1767 
1768 
1769 template <int rank_, int dim, typename Number>
1770 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
1773 {
1774  SymmetricTensor tmp = *this;
1775  tmp.data = -tmp.data;
1776  return tmp;
1777 }
1778 
1779 
1780 
1781 template <int rank_, int dim, typename Number>
1782 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
1784 {
1785  data.clear();
1786 }
1787 
1788 
1789 
1790 template <int rank_, int dim, typename Number>
1791 DEAL_II_HOST constexpr std::size_t
1793 {
1794  // all memory consists of statically allocated memory of the current
1795  // object, no pointers
1796  return sizeof(SymmetricTensor<rank_, dim, Number>);
1797 }
1798 
1799 
1800 
1801 namespace internal
1802 {
1803  template <int dim, typename Number, typename OtherNumber = Number>
1807  perform_double_contraction(
1808  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1809  base_tensor_type &data,
1810  const typename SymmetricTensorAccessors::
1811  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1812  {
1813  using result_type = typename SymmetricTensorAccessors::
1815 
1816  switch (dim)
1817  {
1818  case 1:
1819  return data[0] * sdata[0];
1820  default:
1821  // Start with the non-diagonal part to avoid some multiplications by
1822  // 2.
1823 
1824  result_type sum = data[dim] * sdata[dim];
1825  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1826  sum += data[d] * sdata[d];
1827  sum += sum; // sum = sum * 2.;
1828 
1829  // Now add the contributions from the diagonal
1830  for (unsigned int d = 0; d < dim; ++d)
1831  sum += data[d] * sdata[d];
1832  return sum;
1833  }
1834  }
1835 
1836 
1837 
1838  template <int dim, typename Number, typename OtherNumber = Number>
1842  perform_double_contraction(
1843  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1844  base_tensor_type &data,
1845  const typename SymmetricTensorAccessors::
1846  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1847  {
1848  using result_type = typename SymmetricTensorAccessors::
1850  using value_type = typename SymmetricTensorAccessors::
1852 
1853  const unsigned int data_dim = SymmetricTensorAccessors::
1854  StorageType<2, dim, value_type>::n_independent_components;
1855  value_type tmp[data_dim]{};
1856  for (unsigned int i = 0; i < data_dim; ++i)
1857  tmp[i] =
1858  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1859  return result_type(tmp);
1860  }
1861 
1862 
1863 
1864  template <int dim, typename Number, typename OtherNumber = Number>
1866  typename SymmetricTensorAccessors::StorageType<
1867  2,
1868  dim,
1871  base_tensor_type
1872  perform_double_contraction(
1873  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1874  base_tensor_type &data,
1875  const typename SymmetricTensorAccessors::
1876  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1877  {
1878  using value_type = typename SymmetricTensorAccessors::
1880  using base_tensor_type = typename SymmetricTensorAccessors::
1881  StorageType<2, dim, value_type>::base_tensor_type;
1882 
1883  base_tensor_type tmp;
1884  for (unsigned int i = 0; i < tmp.dimension; ++i)
1885  {
1886  // Start with the non-diagonal part
1887  value_type sum = data[dim] * sdata[dim][i];
1888  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1889  sum += data[d] * sdata[d][i];
1890  sum += sum; // sum = sum * 2.;
1891 
1892  // Now add the contributions from the diagonal
1893  for (unsigned int d = 0; d < dim; ++d)
1894  sum += data[d] * sdata[d][i];
1895  tmp[i] = sum;
1896  }
1897  return tmp;
1898  }
1899 
1900 
1901 
1902  template <int dim, typename Number, typename OtherNumber = Number>
1904  typename SymmetricTensorAccessors::StorageType<
1905  4,
1906  dim,
1909  base_tensor_type
1910  perform_double_contraction(
1911  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1912  base_tensor_type &data,
1913  const typename SymmetricTensorAccessors::
1914  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1915  {
1916  using value_type = typename SymmetricTensorAccessors::
1918  using base_tensor_type = typename SymmetricTensorAccessors::
1919  StorageType<4, dim, value_type>::base_tensor_type;
1920 
1921  const unsigned int data_dim = SymmetricTensorAccessors::
1922  StorageType<2, dim, value_type>::n_independent_components;
1923  base_tensor_type tmp;
1924  for (unsigned int i = 0; i < data_dim; ++i)
1925  for (unsigned int j = 0; j < data_dim; ++j)
1926  {
1927  // Start with the non-diagonal part
1928  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1929  tmp[i][j] += data[i][d] * sdata[d][j];
1930  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1931 
1932  // Now add the contributions from the diagonal
1933  for (unsigned int d = 0; d < dim; ++d)
1934  tmp[i][j] += data[i][d] * sdata[d][j];
1935  }
1936  return tmp;
1937  }
1938 
1939 } // end of namespace internal
1940 
1941 
1942 
1943 template <int rank_, int dim, typename Number>
1944 template <typename OtherNumber>
1949  const SymmetricTensor<2, dim, OtherNumber> &s) const
1950 {
1951  // need to have two different function calls
1952  // because a scalar and rank-2 tensor are not
1953  // the same data type (see internal function
1954  // above)
1955  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1956  s.data);
1957 }
1958 
1959 
1960 
1961 template <int rank_, int dim, typename Number>
1962 template <typename OtherNumber>
1967  const SymmetricTensor<4, dim, OtherNumber> &s) const
1968 {
1971  tmp.data =
1972  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1973  s.data);
1974  return tmp;
1975 }
1976 
1977 
1978 
1979 // internal namespace to switch between the
1980 // access of different tensors. There used to
1981 // be explicit instantiations before for
1982 // different ranks and dimensions, but since
1983 // we now allow for templates on the data
1984 // type, and since we cannot partially
1985 // specialize the implementation, this got
1986 // into a separate namespace
1987 namespace internal
1988 {
1989  template <int dim, typename Number>
1990  constexpr inline DEAL_II_ALWAYS_INLINE Number &
1991  symmetric_tensor_access(const TableIndices<2> &indices,
1992  typename SymmetricTensorAccessors::
1993  StorageType<2, dim, Number>::base_tensor_type &data)
1994  {
1995  // 1d is very simple and done first
1996  if (dim == 1)
1997  return data[0];
1998 
1999  // first treat the main diagonal elements, which are stored consecutively
2000  // at the beginning
2001  if (indices[0] == indices[1])
2002  return data[indices[0]];
2003 
2004  // the rest is messier and requires a few switches.
2005  switch (dim)
2006  {
2007  case 2:
2008  // at least for the 2x2 case it is reasonably simple
2009  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
2010  ((indices[0] == 0) && (indices[1] == 1)),
2011  ExcInternalError());
2012  return data[2];
2013 
2014  default:
2015  // to do the rest, sort our indices before comparing
2016  {
2017  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2018  std::max(indices[0], indices[1]));
2019  for (unsigned int d = 0, c = 0; d < dim; ++d)
2020  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2021  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2022  return data[dim + c];
2023  Assert(false, ExcInternalError());
2024  }
2025  }
2026 
2027  // The code should never reach here.
2028  // We cannot return a static variable, as this class must support number
2029  // types that require no instances of the number type to be in scope during
2030  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2031  return data[0];
2032  }
2033 
2034 
2035 
2036  template <int dim, typename Number>
2037  constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2038  symmetric_tensor_access(const TableIndices<2> &indices,
2039  const typename SymmetricTensorAccessors::
2040  StorageType<2, dim, Number>::base_tensor_type &data)
2041  {
2042  // 1d is very simple and done first
2043  if (dim == 1)
2044  return data[0];
2045 
2046  // first treat the main diagonal elements, which are stored consecutively
2047  // at the beginning
2048  if (indices[0] == indices[1])
2049  return data[indices[0]];
2050 
2051  // the rest is messier and requires a few switches.
2052  switch (dim)
2053  {
2054  case 2:
2055  // at least for the 2x2 case it is reasonably simple
2056  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
2057  ((indices[0] == 0) && (indices[1] == 1)),
2058  ExcInternalError());
2059  return data[2];
2060 
2061  default:
2062  // to do the rest, sort our indices before comparing
2063  {
2064  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
2065  std::max(indices[0], indices[1]));
2066  for (unsigned int d = 0, c = 0; d < dim; ++d)
2067  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2068  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2069  return data[dim + c];
2070  Assert(false, ExcInternalError());
2071  }
2072  }
2073 
2074  // The code should never reach here.
2075  // We cannot return a static variable, as this class must support number
2076  // types that require no instances of the number type to be in scope during
2077  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2078  return data[0];
2079  }
2080 
2081 
2082 
2083  template <int dim, typename Number>
2084  constexpr inline Number &
2085  symmetric_tensor_access(const TableIndices<4> &indices,
2086  typename SymmetricTensorAccessors::
2087  StorageType<4, dim, Number>::base_tensor_type &data)
2088  {
2089  switch (dim)
2090  {
2091  case 1:
2092  return data[0][0];
2093 
2094  case 2:
2095  // each entry of the tensor can be thought of as an entry in a
2096  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2097  // rank-2 tensors. this is the format in which we store rank-4
2098  // tensors. determine which position the present entry is
2099  // stored in
2100  {
2101  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2102  return data[base_index[indices[0]][indices[1]]]
2103  [base_index[indices[2]][indices[3]]];
2104  }
2105  case 3:
2106  // each entry of the tensor can be thought of as an entry in a
2107  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2108  // rank-2 tensors. this is the format in which we store rank-4
2109  // tensors. determine which position the present entry is
2110  // stored in
2111  {
2112  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2113  {3, 1, 5},
2114  {4, 5, 2}};
2115  return data[base_index[indices[0]][indices[1]]]
2116  [base_index[indices[2]][indices[3]]];
2117  }
2118 
2119  default:
2120  Assert(false, ExcNotImplemented());
2121  }
2122 
2123  // The code should never reach here.
2124  // We cannot return a static variable, as this class must support number
2125  // types that require no instances of the number type to be in scope during
2126  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2127  return data[0][0];
2128  }
2129 
2130 
2131  template <int dim, typename Number>
2132  constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2133  symmetric_tensor_access(const TableIndices<4> &indices,
2134  const typename SymmetricTensorAccessors::
2135  StorageType<4, dim, Number>::base_tensor_type &data)
2136  {
2137  switch (dim)
2138  {
2139  case 1:
2140  return data[0][0];
2141 
2142  case 2:
2143  // each entry of the tensor can be thought of as an entry in a
2144  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2145  // rank-2 tensors. this is the format in which we store rank-4
2146  // tensors. determine which position the present entry is
2147  // stored in
2148  {
2149  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
2150  return data[base_index[indices[0]][indices[1]]]
2151  [base_index[indices[2]][indices[3]]];
2152  }
2153  case 3:
2154  // each entry of the tensor can be thought of as an entry in a
2155  // matrix that maps the rolled-out rank-2 tensors into rolled-out
2156  // rank-2 tensors. this is the format in which we store rank-4
2157  // tensors. determine which position the present entry is
2158  // stored in
2159  {
2160  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
2161  {3, 1, 5},
2162  {4, 5, 2}};
2163  return data[base_index[indices[0]][indices[1]]]
2164  [base_index[indices[2]][indices[3]]];
2165  }
2166 
2167  default:
2168  Assert(false, ExcNotImplemented());
2169  }
2170 
2171  // The code should never reach here.
2172  // We cannot return a static variable, as this class must support number
2173  // types that require no instances of the number type to be in scope during
2174  // a reinitialization procedure (e.g. ADOL-C adtl::adouble).
2175  return data[0][0];
2176  }
2177 
2178 } // end of namespace internal
2179 
2180 
2181 
2182 template <int rank_, int dim, typename Number>
2183 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2185  const TableIndices<rank_> &indices)
2186 {
2187  for (unsigned int r = 0; r < rank; ++r)
2188  AssertIndexRange(indices[r], dimension);
2189  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2190 }
2191 
2192 
2193 
2194 template <int rank_, int dim, typename Number>
2195 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE const Number &
2197  const TableIndices<rank_> &indices) const
2198 {
2199  for (unsigned int r = 0; r < rank; ++r)
2200  AssertIndexRange(indices[r], dimension);
2201  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2202 }
2203 
2204 
2205 
2206 namespace internal
2207 {
2208  namespace SymmetricTensorImplementation
2209  {
2210  template <int rank_>
2211  constexpr TableIndices<rank_>
2212  get_partially_filled_indices(const unsigned int row,
2213  const std::integral_constant<int, 2> &)
2214  {
2216  }
2217 
2218 
2219  template <int rank_>
2220  constexpr TableIndices<rank_>
2221  get_partially_filled_indices(const unsigned int row,
2222  const std::integral_constant<int, 4> &)
2223  {
2224  return TableIndices<rank_>(row,
2228  }
2229  } // namespace SymmetricTensorImplementation
2230 } // namespace internal
2231 
2232 
2233 template <int rank_, int dim, typename Number>
2234 DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE internal::
2235  SymmetricTensorAccessors::Accessor<rank_, dim, true, rank_ - 1, Number>
2236  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row) const
2237 {
2238  return internal::SymmetricTensorAccessors::
2239  Accessor<rank_, dim, true, rank_ - 1, Number>(
2240  *this,
2241  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2242  rank_>(row, std::integral_constant<int, rank_>()));
2243 }
2244 
2245 
2246 
2247 template <int rank_, int dim, typename Number>
2248 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE internal::
2249  SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2250  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2251 {
2252  return internal::SymmetricTensorAccessors::
2253  Accessor<rank_, dim, false, rank_ - 1, Number>(
2254  *this,
2255  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2256  rank_>(row, std::integral_constant<int, rank_>()));
2257 }
2258 
2259 
2260 
2261 template <int rank_, int dim, typename Number>
2262 DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE const Number &
2264  const TableIndices<rank_> &indices) const
2265 {
2266  return operator()(indices);
2267 }
2268 
2269 
2270 
2271 template <int rank_, int dim, typename Number>
2272 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number &
2274  const TableIndices<rank_> &indices)
2275 {
2276  return operator()(indices);
2277 }
2278 
2279 
2280 
2281 template <int rank_, int dim, typename Number>
2282 inline Number *
2284 {
2285  return std::addressof(this->access_raw_entry(0));
2286 }
2287 
2288 
2289 
2290 template <int rank_, int dim, typename Number>
2291 inline const Number *
2293 {
2294  return std::addressof(this->access_raw_entry(0));
2295 }
2296 
2297 
2298 
2299 template <int rank_, int dim, typename Number>
2300 inline Number *
2302 {
2303  return begin_raw() + n_independent_components;
2304 }
2305 
2306 
2307 
2308 template <int rank_, int dim, typename Number>
2309 inline const Number *
2311 {
2312  return begin_raw() + n_independent_components;
2313 }
2314 
2315 
2316 
2317 namespace internal
2318 {
2319  namespace SymmetricTensorImplementation
2320  {
2321  template <int dim, typename Number>
2322  constexpr unsigned int
2323  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2324  const unsigned int index)
2325  {
2326  return index;
2327  }
2328 
2329 
2330  template <int dim, typename Number>
2331  constexpr ::TableIndices<2>
2332  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2333  const unsigned int index)
2334  {
2337  }
2338 
2339  } // namespace SymmetricTensorImplementation
2340 } // namespace internal
2341 
2342 
2343 
2344 template <int rank_, int dim, typename Number>
2345 DEAL_II_HOST constexpr inline const Number &
2347  const unsigned int index) const
2348 {
2349  AssertIndexRange(index, n_independent_components);
2350  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2351  index)];
2352 }
2353 
2354 
2355 
2356 template <int rank_, int dim, typename Number>
2357 DEAL_II_HOST constexpr inline Number &
2359 {
2360  AssertIndexRange(index, n_independent_components);
2361  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2362  index)];
2363 }
2364 
2365 
2366 
2367 namespace internal
2368 {
2369  template <int dim, typename Number>
2370  constexpr inline typename numbers::NumberTraits<Number>::real_type
2371  compute_norm(const typename SymmetricTensorAccessors::
2372  StorageType<2, dim, Number>::base_tensor_type &data)
2373  {
2374  // Make things work with AD types
2375  using std::sqrt;
2376  switch (dim)
2377  {
2378  case 1:
2379  return numbers::NumberTraits<Number>::abs(data[0]);
2380 
2381  case 2:
2385 
2386  case 3:
2393 
2394  default:
2395  {
2396  typename numbers::NumberTraits<Number>::real_type return_value =
2398 
2399  for (unsigned int d = 0; d < dim; ++d)
2400  return_value +=
2402  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2403  return_value +=
2405 
2406  return sqrt(return_value);
2407  }
2408  }
2409  }
2410 
2411 
2412 
2413  template <int dim, typename Number>
2414  constexpr inline typename numbers::NumberTraits<Number>::real_type
2415  compute_norm(const typename SymmetricTensorAccessors::
2416  StorageType<4, dim, Number>::base_tensor_type &data)
2417  {
2418  // Make things work with AD types
2419  using std::sqrt;
2420  switch (dim)
2421  {
2422  case 1:
2423  return numbers::NumberTraits<Number>::abs(data[0][0]);
2424 
2425  default:
2426  {
2427  typename numbers::NumberTraits<Number>::real_type return_value =
2429 
2430  const unsigned int n_independent_components = data.dimension;
2431 
2432  for (unsigned int i = 0; i < dim; ++i)
2433  for (unsigned int j = 0; j < dim; ++j)
2434  return_value +=
2436  for (unsigned int i = 0; i < dim; ++i)
2437  for (unsigned int j = dim; j < n_independent_components; ++j)
2438  return_value +=
2440  for (unsigned int i = dim; i < n_independent_components; ++i)
2441  for (unsigned int j = 0; j < dim; ++j)
2442  return_value +=
2444  for (unsigned int i = dim; i < n_independent_components; ++i)
2445  for (unsigned int j = dim; j < n_independent_components; ++j)
2446  return_value +=
2448 
2449  return sqrt(return_value);
2450  }
2451  }
2452  }
2453 
2454 } // end of namespace internal
2455 
2456 
2457 
2458 template <int rank_, int dim, typename Number>
2461 {
2462  return internal::compute_norm<dim, Number>(data);
2463 }
2464 
2465 
2466 
2467 namespace internal
2468 {
2469  namespace SymmetricTensorImplementation
2470  {
2471  // a function to do the unrolling from a set of indices to a
2472  // scalar index into the array in which we store the elements of
2473  // a symmetric tensor
2474  //
2475  // this function is for rank-2 tensors
2476  template <int dim>
2477  constexpr inline DEAL_II_ALWAYS_INLINE unsigned int
2479  {
2480  AssertIndexRange(indices[0], dim);
2481  AssertIndexRange(indices[1], dim);
2482 
2483  switch (dim)
2484  {
2485  case 1:
2486  {
2487  return 0;
2488  }
2489 
2490  case 2:
2491  {
2492  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2493  return table[indices[0]][indices[1]];
2494  }
2495 
2496  case 3:
2497  {
2498  constexpr unsigned int table[3][3] = {{0, 3, 4},
2499  {3, 1, 5},
2500  {4, 5, 2}};
2501  return table[indices[0]][indices[1]];
2502  }
2503 
2504  case 4:
2505  {
2506  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2507  {4, 1, 7, 8},
2508  {5, 7, 2, 9},
2509  {6, 8, 9, 3}};
2510  return table[indices[0]][indices[1]];
2511  }
2512 
2513  default:
2514  // for the remainder, manually figure out the numbering
2515  {
2516  if (indices[0] == indices[1])
2517  return indices[0];
2518 
2519  TableIndices<2> sorted_indices(indices);
2520  sorted_indices.sort();
2521 
2522  for (unsigned int d = 0, c = 0; d < dim; ++d)
2523  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2524  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2525  return dim + c;
2526 
2527  // should never get here:
2528  Assert(false, ExcInternalError());
2529  return 0;
2530  }
2531  }
2532  }
2533 
2534  // a function to do the unrolling from a set of indices to a
2535  // scalar index into the array in which we store the elements of
2536  // a symmetric tensor
2537  //
2538  // this function is for tensors of ranks not already handled
2539  // above
2540  template <int dim, int rank_>
2541  constexpr inline unsigned int
2543  {
2544  (void)indices;
2545  Assert(false, ExcNotImplemented());
2547  }
2548  } // namespace SymmetricTensorImplementation
2549 } // namespace internal
2550 
2551 
2552 template <int rank_, int dim, typename Number>
2553 DEAL_II_HOST constexpr unsigned int
2555  const TableIndices<rank_> &indices)
2556 {
2557  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2558  dim>(indices);
2559 }
2560 
2561 
2562 
2563 namespace internal
2564 {
2565  namespace SymmetricTensorImplementation
2566  {
2567  // a function to do the inverse of the unrolling from a set of
2568  // indices to a scalar index into the array in which we store
2569  // the elements of a symmetric tensor. in other words, it goes
2570  // from the scalar index into the array to a set of indices of
2571  // the tensor
2572  //
2573  // this function is for rank-2 tensors
2574  template <int dim>
2575  constexpr inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2576  unrolled_to_component_indices(const unsigned int i,
2577  const std::integral_constant<int, 2> &)
2578  {
2579  Assert(
2581  ExcIndexRange(
2582  i,
2583  0,
2585  switch (dim)
2586  {
2587  case 1:
2588  {
2589  return {0, 0};
2590  }
2591 
2592  case 2:
2593  {
2594  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2595  TableIndices<2>(1, 1),
2596  TableIndices<2>(0, 1)};
2597  return table[i];
2598  }
2599 
2600  case 3:
2601  {
2602  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2603  TableIndices<2>(1, 1),
2604  TableIndices<2>(2, 2),
2605  TableIndices<2>(0, 1),
2606  TableIndices<2>(0, 2),
2607  TableIndices<2>(1, 2)};
2608  return table[i];
2609  }
2610 
2611  default:
2612  if (i < dim)
2613  return {i, i};
2614 
2615  for (unsigned int d = 0, c = dim; d < dim; ++d)
2616  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2617  if (c == i)
2618  return {d, e};
2619 
2620  // should never get here:
2621  Assert(false, ExcInternalError());
2622  return {0, 0};
2623  }
2624  }
2625 
2626  // a function to do the inverse of the unrolling from a set of
2627  // indices to a scalar index into the array in which we store
2628  // the elements of a symmetric tensor. in other words, it goes
2629  // from the scalar index into the array to a set of indices of
2630  // the tensor
2631  //
2632  // this function is for tensors of a rank not already handled
2633  // above
2634  template <int dim, int rank_>
2635  constexpr inline std::enable_if_t<rank_ != 2, TableIndices<rank_>>
2636  unrolled_to_component_indices(const unsigned int i,
2637  const std::integral_constant<int, rank_> &)
2638  {
2639  (void)i;
2640  Assert(
2641  (i <
2643  ExcIndexRange(i,
2644  0,
2646  n_independent_components));
2647  Assert(false, ExcNotImplemented());
2648  return TableIndices<rank_>();
2649  }
2650 
2651  } // namespace SymmetricTensorImplementation
2652 } // namespace internal
2653 
2654 template <int rank_, int dim, typename Number>
2657  const unsigned int i)
2658 {
2659  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2660  dim>(i, std::integral_constant<int, rank_>());
2661 }
2662 
2663 
2664 
2665 template <int rank_, int dim, typename Number>
2666 template <class Archive>
2667 inline void
2668 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2669 {
2670  ar &data;
2671 }
2672 
2673 
2674 #endif // DOXYGEN
2675 
2676 /* ----------------- Non-member functions operating on tensors. ------------ */
2677 
2678 
2691 template <int rank_, int dim, typename Number, typename OtherNumber>
2692 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2696 {
2698  tmp = left;
2699  tmp += right;
2700  return tmp;
2701 }
2702 
2703 
2716 template <int rank_, int dim, typename Number, typename OtherNumber>
2717 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
2721 {
2723  tmp = left;
2724  tmp -= right;
2725  return tmp;
2726 }
2727 
2728 
2736 template <int rank_, int dim, typename Number, typename OtherNumber>
2740  const Tensor<rank_, dim, OtherNumber> & right)
2741 {
2742  return Tensor<rank_, dim, Number>(left) + right;
2743 }
2744 
2745 
2753 template <int rank_, int dim, typename Number, typename OtherNumber>
2758 {
2759  return left + Tensor<rank_, dim, OtherNumber>(right);
2760 }
2761 
2762 
2770 template <int rank_, int dim, typename Number, typename OtherNumber>
2774  const Tensor<rank_, dim, OtherNumber> & right)
2775 {
2776  return Tensor<rank_, dim, Number>(left) - right;
2777 }
2778 
2779 
2787 template <int rank_, int dim, typename Number, typename OtherNumber>
2792 {
2793  return left - Tensor<rank_, dim, OtherNumber>(right);
2794 }
2795 
2796 
2797 
2798 template <int dim, typename Number>
2801 {
2802  switch (dim)
2803  {
2804  case 1:
2805  return t.data[0];
2806  case 2:
2807  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2808  case 3:
2809  {
2810  // in analogy to general tensors, but
2811  // there's something to be simplified for
2812  // the present case
2813  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2814  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2815  t.data[0] * t.data[5] * t.data[5] -
2816  t.data[1] * t.data[4] * t.data[4] -
2817  t.data[2] * t.data[3] * t.data[3]);
2818  }
2819  default:
2820  Assert(false, ExcNotImplemented());
2822  }
2823 }
2824 
2825 
2826 
2838 template <int dim, typename Number>
2841 {
2842  return determinant(t);
2843 }
2844 
2845 
2846 
2847 template <int dim, typename Number>
2848 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE Number
2850 {
2851  Number t = d.data[0];
2852  for (unsigned int i = 1; i < dim; ++i)
2853  t += d.data[i];
2854  return t;
2855 }
2856 
2857 
2869 template <int dim, typename Number>
2870 DEAL_II_HOST constexpr Number
2872 {
2873  return trace(t);
2874 }
2875 
2876 
2888 template <typename Number>
2889 DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2891 {
2893 }
2894 
2895 
2896 
2915 template <typename Number>
2916 DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2918 {
2919  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2920 }
2921 
2922 
2923 
2932 template <typename Number>
2933 DEAL_II_HOST constexpr DEAL_II_ALWAYS_INLINE Number
2935 {
2936  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2937  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2938 }
2939 
2940 
2941 
2949 template <typename Number>
2950 std::array<Number, 1>
2952 
2953 
2954 
2977 template <typename Number>
2978 std::array<Number, 2>
2980 
2981 
2982 
3005 template <typename Number>
3006 std::array<Number, 3>
3008 
3009 
3010 
3011 namespace internal
3012 {
3013  namespace SymmetricTensorImplementation
3014  {
3052  template <int dim, typename Number>
3053  void
3054  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
3055  ::Tensor<2, dim, Number> & Q,
3056  std::array<Number, dim> & d,
3057  std::array<Number, dim - 1> & e);
3058 
3059 
3060 
3100  template <int dim, typename Number>
3101  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3102  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
3103 
3104 
3105 
3145  template <int dim, typename Number>
3146  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3148 
3149 
3150 
3164  template <typename Number>
3165  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3166  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3167 
3168 
3169 
3202  template <typename Number>
3203  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3204  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3205 
3210  template <int dim, typename Number>
3212  {
3213  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3214  bool
3215  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3216  {
3217  return lhs.first > rhs.first;
3218  }
3219  };
3220 
3221  } // namespace SymmetricTensorImplementation
3222 
3223 } // namespace internal
3224 
3225 
3226 
3227 // The line below is to ensure that doxygen puts the full description
3228 // of this global enumeration into the documentation
3229 // See https://stackoverflow.com/a/1717984
3259 {
3269  hybrid,
3287  jacobi
3288 };
3289 
3290 
3291 
3320 template <int dim, typename Number>
3321 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3322  std::integral_constant<int, dim>::value>
3324  const SymmetricTensorEigenvectorMethod method =
3326 
3327 
3328 
3337 template <int rank_, int dim, typename Number>
3340 {
3341  return t;
3342 }
3343 
3344 
3345 
3346 template <int dim, typename Number>
3347 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3350 {
3352 
3353  // subtract scaled trace from the diagonal
3354  const Number tr = trace(t) * internal::NumberType<Number>::value(1.0 / dim);
3355  for (unsigned int i = 0; i < dim; ++i)
3356  tmp.data[i] -= tr;
3357 
3358  return tmp;
3359 }
3360 
3361 
3362 
3363 template <int dim, typename Number>
3367 {
3368  // create a default constructed matrix filled with
3369  // zeros, then set the diagonal elements to one
3371  switch (dim)
3372  {
3373  case 1:
3375  break;
3376  case 2:
3377  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3378  break;
3379  case 3:
3380  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3382  break;
3383  default:
3384  for (unsigned int d = 0; d < dim; ++d)
3386  }
3387  return tmp;
3388 }
3389 
3390 
3391 
3392 template <int dim, typename Number>
3395 {
3397 
3398  // fill the elements treating the diagonal
3399  for (unsigned int i = 0; i < dim; ++i)
3400  for (unsigned int j = 0; j < dim; ++j)
3401  tmp.data[i][j] =
3402  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3403 
3404  // then fill the ones that copy over the
3405  // non-diagonal elements. note that during
3406  // the double-contraction, we handle the
3407  // off-diagonal elements twice, so simply
3408  // copying requires a weight of 1/2
3409  for (unsigned int i = dim;
3410  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3411  n_rank2_components;
3412  ++i)
3413  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3414 
3415  return tmp;
3416 }
3417 
3418 
3419 
3420 template <int dim, typename Number>
3424 {
3426 
3427  // fill the elements treating the diagonal
3428  for (unsigned int i = 0; i < dim; ++i)
3429  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3430 
3431  // then fill the ones that copy over the
3432  // non-diagonal elements. note that during
3433  // the double-contraction, we handle the
3434  // off-diagonal elements twice, so simply
3435  // copying requires a weight of 1/2
3436  for (unsigned int i = dim;
3437  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3438  n_rank2_components;
3439  ++i)
3440  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3441 
3442  return tmp;
3443 }
3444 
3445 
3446 
3456 template <int dim, typename Number>
3459 {
3461  value(t);
3462 }
3463 
3464 
3465 
3476 template <int dim, typename Number>
3479 {
3481  value(t);
3482 }
3483 
3484 
3485 
3507 template <int dim, typename Number>
3511 {
3513 
3514  // fill only the elements really needed
3515  for (unsigned int i = 0; i < dim; ++i)
3516  for (unsigned int j = i; j < dim; ++j)
3517  for (unsigned int k = 0; k < dim; ++k)
3518  for (unsigned int l = k; l < dim; ++l)
3519  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3520 
3521  return tmp;
3522 }
3523 
3524 
3525 
3533 template <int dim, typename Number>
3534 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3537 {
3539  for (unsigned int d = 0; d < dim; ++d)
3540  result[d][d] = t[d][d];
3541 
3542  const Number half = internal::NumberType<Number>::value(0.5);
3543  for (unsigned int d = 0; d < dim; ++d)
3544  for (unsigned int e = d + 1; e < dim; ++e)
3545  result[d][e] = (t[d][e] + t[e][d]) * half;
3546  return result;
3547 }
3548 
3549 
3550 
3558 template <int rank_, int dim, typename Number>
3559 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3561  operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3562 {
3564  tt *= factor;
3565  return tt;
3566 }
3567 
3568 
3569 
3577 template <int rank_, int dim, typename Number>
3579 operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3580 {
3581  // simply forward to the other operator
3582  return t * factor;
3583 }
3584 
3585 
3586 
3611 template <int rank_, int dim, typename Number, typename OtherNumber>
3613  rank_,
3614  dim,
3615  typename ProductType<Number,
3616  typename EnableIfScalar<OtherNumber>::type>::type>
3618  const OtherNumber & factor)
3619 {
3620  // form the product. we have to convert the two factors into the final
3621  // type via explicit casts because, for awkward reasons, the C++
3622  // standard committee saw it fit to not define an
3623  // operator*(float,std::complex<double>)
3624  // (as well as with switched arguments and double<->float).
3625  using product_type = typename ProductType<Number, OtherNumber>::type;
3628  return tt;
3629 }
3630 
3631 
3632 
3640 template <int rank_, int dim, typename Number, typename OtherNumber>
3642  rank_,
3643  dim,
3644  typename ProductType<OtherNumber,
3645  typename EnableIfScalar<Number>::type>::type>
3646 operator*(const Number & factor,
3648 {
3649  // simply forward to the other operator with switched arguments
3650  return (t * factor);
3651 }
3652 
3653 
3654 
3660 template <int rank_, int dim, typename Number, typename OtherNumber>
3661 DEAL_II_HOST constexpr inline SymmetricTensor<
3662  rank_,
3663  dim,
3664  typename ProductType<Number,
3665  typename EnableIfScalar<OtherNumber>::type>::type>
3667  const OtherNumber & factor)
3668 {
3669  using product_type = typename ProductType<Number, OtherNumber>::type;
3672  return tt;
3673 }
3674 
3675 
3676 
3683 template <int rank_, int dim>
3685 operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3686 {
3688  tt *= factor;
3689  return tt;
3690 }
3691 
3692 
3693 
3700 template <int rank_, int dim>
3702 operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3703 {
3705  tt *= factor;
3706  return tt;
3707 }
3708 
3709 
3710 
3716 template <int rank_, int dim>
3718 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3719 {
3721  tt /= factor;
3722  return tt;
3723 }
3724 
3734 template <int dim, typename Number, typename OtherNumber>
3739 {
3740  return (t1 * t2);
3741 }
3742 
3743 
3756 template <int dim, typename Number, typename OtherNumber>
3757 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE
3760  const Tensor<2, dim, OtherNumber> & t2)
3761 {
3763  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3764  for (unsigned int i = 0; i < dim; ++i)
3765  for (unsigned int j = 0; j < dim; ++j)
3766  s += t1[i][j] * t2[i][j];
3767  return s;
3768 }
3769 
3770 
3783 template <int dim, typename Number, typename OtherNumber>
3788 {
3789  return scalar_product(t2, t1);
3790 }
3791 
3792 
3807 template <typename Number, typename OtherNumber>
3808 DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE void
3813 {
3814  tmp[0][0] = t[0][0][0][0] * s[0][0];
3815 }
3816 
3817 
3818 
3833 template <typename Number, typename OtherNumber>
3834 DEAL_II_HOST constexpr inline void
3839 {
3840  tmp[0][0] = t[0][0][0][0] * s[0][0];
3841 }
3842 
3843 
3844 
3859 template <typename Number, typename OtherNumber>
3860 DEAL_II_HOST constexpr inline void
3865 {
3866  const unsigned int dim = 2;
3867 
3868  for (unsigned int i = 0; i < dim; ++i)
3869  for (unsigned int j = i; j < dim; ++j)
3870  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3871  2 * t[i][j][0][1] * s[0][1];
3872 }
3873 
3874 
3875 
3890 template <typename Number, typename OtherNumber>
3891 DEAL_II_HOST constexpr inline void
3896 {
3897  const unsigned int dim = 2;
3898 
3899  for (unsigned int i = 0; i < dim; ++i)
3900  for (unsigned int j = i; j < dim; ++j)
3901  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3902  2 * s[0][1] * t[0][1][i][j];
3903 }
3904 
3905 
3906 
3921 template <typename Number, typename OtherNumber>
3922 DEAL_II_HOST constexpr inline void
3927 {
3928  const unsigned int dim = 3;
3929 
3930  for (unsigned int i = 0; i < dim; ++i)
3931  for (unsigned int j = i; j < dim; ++j)
3932  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3933  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3934  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3935 }
3936 
3937 
3938 
3953 template <typename Number, typename OtherNumber>
3954 DEAL_II_HOST constexpr inline void
3959 {
3960  const unsigned int dim = 3;
3961 
3962  for (unsigned int i = 0; i < dim; ++i)
3963  for (unsigned int j = i; j < dim; ++j)
3964  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3965  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3966  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3967 }
3968 
3969 
3970 
3977 template <int dim, typename Number, typename OtherNumber>
3978 DEAL_II_HOST constexpr Tensor<1,
3979  dim,
3982  const Tensor<1, dim, OtherNumber> & src2)
3983 {
3985  for (unsigned int i = 0; i < dim; ++i)
3986  for (unsigned int j = 0; j < dim; ++j)
3987  dest[i] += src1[i][j] * src2[j];
3988  return dest;
3989 }
3990 
3991 
3998 template <int dim, typename Number, typename OtherNumber>
3999 DEAL_II_HOST constexpr Tensor<1,
4000  dim,
4004 {
4005  // this is easy for symmetric tensors:
4006  return src2 * src1;
4007 }
4008 
4009 
4010 
4030 template <int rank_1,
4031  int rank_2,
4032  int dim,
4033  typename Number,
4034  typename OtherNumber>
4036  typename Tensor<rank_1 + rank_2 - 2,
4037  dim,
4038  typename ProductType<Number, OtherNumber>::type>::tensor_type
4041 {
4042  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4043 }
4044 
4045 
4046 
4066 template <int rank_1,
4067  int rank_2,
4068  int dim,
4069  typename Number,
4070  typename OtherNumber>
4072  typename Tensor<rank_1 + rank_2 - 2,
4073  dim,
4074  typename ProductType<Number, OtherNumber>::type>::tensor_type
4076  const Tensor<rank_2, dim, OtherNumber> & src2)
4077 {
4078  return Tensor<rank_1, dim, Number>(src1) * src2;
4079 }
4080 
4081 
4082 
4092 template <int dim, typename Number>
4093 inline std::ostream &
4094 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4095 {
4096  // make our lives a bit simpler by outputting
4097  // the tensor through the operator for the
4098  // general Tensor class
4100 
4101  for (unsigned int i = 0; i < dim; ++i)
4102  for (unsigned int j = 0; j < dim; ++j)
4103  tt[i][j] = t[i][j];
4104 
4105  return out << tt;
4106 }
4107 
4108 
4109 
4119 template <int dim, typename Number>
4120 inline std::ostream &
4121 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4122 {
4123  // make our lives a bit simpler by outputting
4124  // the tensor through the operator for the
4125  // general Tensor class
4127 
4128  for (unsigned int i = 0; i < dim; ++i)
4129  for (unsigned int j = 0; j < dim; ++j)
4130  for (unsigned int k = 0; k < dim; ++k)
4131  for (unsigned int l = 0; l < dim; ++l)
4132  tt[i][j][k][l] = t[i][j][k][l];
4133 
4134  return out << tt;
4135 }
4136 
4137 
4139 
4140 #endif
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
constexpr DEAL_II_HOST const Number & operator[](const TableIndices< rank_ > &indices) const
constexpr DEAL_II_HOST Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const Tensor< 1, dim, Number > &src1, const SymmetricTensor< 2, dim, OtherNumber > &src2)
constexpr DEAL_II_HOST 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 DEAL_II_HOST internal::SymmetricTensorAccessors::Accessor< rank_, dim, false, rank_ - 1, Number > operator[](const unsigned int row)
constexpr DEAL_II_HOST 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 DEAL_II_HOST Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
constexpr DEAL_II_HOST Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
constexpr DEAL_II_HOST Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr DEAL_II_HOST Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Number * begin_raw()
constexpr DEAL_II_HOST Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Number * end_raw()
constexpr DEAL_II_HOST Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr DEAL_II_HOST SymmetricTensor(const Number(&array)[n_independent_components])
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
constexpr DEAL_II_HOST internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
constexpr DEAL_II_HOST Number & operator()(const TableIndices< rank_ > &indices)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr DEAL_II_HOST ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, OtherNumber > &t2)
std::array< Number, 2 > eigenvalues(const SymmetricTensor< 2, 2, Number > &T)
void serialize(Archive &ar, const unsigned int version)
constexpr DEAL_II_HOST SymmetricTensor()=default
constexpr DEAL_II_HOST SymmetricTensor(const SymmetricTensor< rank_, dim, OtherNumber > &initializer)
constexpr DEAL_II_HOST SymmetricTensor operator-() const
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr DEAL_II_HOST bool operator!=(const SymmetricTensor &) const
static constexpr DEAL_II_HOST std::size_t memory_consumption()
constexpr DEAL_II_HOST Number trace(const SymmetricTensor< 2, dim2, Number > &)
typename base_tensor_descriptor::base_tensor_type base_tensor_type
constexpr DEAL_II_HOST SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr DEAL_II_HOST internal::SymmetricTensorAccessors::double_contraction_result< rank_, 4, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 4, dim, OtherNumber > &s) const
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr DEAL_II_HOST void clear()
const Number * end_raw() const
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr DEAL_II_HOST numbers::NumberTraits< Number >::real_type norm() const
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, Number > operator*(const Number &factor, const SymmetricTensor< rank_, dim, Number > &t)
constexpr DEAL_II_HOST SymmetricTensor & operator=(const Number &d)
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > identity_tensor()
constexpr DEAL_II_HOST void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, OtherNumber > &s)
constexpr DEAL_II_HOST 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 DEAL_II_HOST const Number & operator()(const TableIndices< rank_ > &indices) const
constexpr DEAL_II_HOST Tensor< 1, dim, typename ProductType< Number, OtherNumber >::type > operator*(const SymmetricTensor< 2, dim, Number > &src1, const Tensor< 1, dim, OtherNumber > &src2)
base_tensor_type data
constexpr DEAL_II_HOST void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
const Number * begin_raw() const
SymmetricTensor(const Tensor< 2, dim, OtherNumber > &t)
constexpr DEAL_II_HOST ProductType< Number, OtherNumber >::type scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr DEAL_II_HOST Number & operator[](const TableIndices< rank_ > &indices)
constexpr DEAL_II_HOST Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, typename ProductType< OtherNumber, typename EnableIfScalar< Number >::type >::type > operator*(const Number &factor, const SymmetricTensor< rank_, dim, OtherNumber > &t)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
constexpr DEAL_II_HOST SymmetricTensor & operator/=(const OtherNumber &factor)
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
constexpr DEAL_II_HOST void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, OtherNumber > &t)
constexpr DEAL_II_HOST const Number & access_raw_entry(const unsigned int unrolled_index) const
constexpr DEAL_II_HOST bool operator==(const SymmetricTensor &) const
constexpr DEAL_II_HOST Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
static constexpr DEAL_II_HOST unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
constexpr DEAL_II_HOST 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)
static constexpr DEAL_II_HOST TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim > operator/(const SymmetricTensor< rank_, dim > &t, const double factor)
constexpr DEAL_II_HOST ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr DEAL_II_HOST SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim > operator*(const double factor, const SymmetricTensor< rank_, dim > &t)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
std::array< Number, 3 > eigenvalues(const SymmetricTensor< 2, 3, Number > &T)
constexpr DEAL_II_HOST SymmetricTensor< rank_, dim > operator*(const SymmetricTensor< rank_, dim > &t, const double factor)
constexpr DEAL_II_HOST Number & access_raw_entry(const unsigned int unrolled_index)
constexpr DEAL_II_HOST Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr DEAL_II_HOST 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 DEAL_II_HOST SymmetricTensor & operator*=(const OtherNumber &factor)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: tensor.h:516
constexpr DEAL_II_HOST reference operator[](const unsigned int)
typename AccessorTypes< rank, dim, constness, Number >::reference reference
constexpr DEAL_II_HOST Accessor(tensor_type &tensor, const TableIndices< rank > &previous_indices)
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
constexpr DEAL_II_HOST reference operator[](const unsigned int) const
constexpr DEAL_II_HOST Accessor(tensor_type &tensor, const TableIndices< rank > &previous_indices)
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
constexpr DEAL_II_HOST Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i) const
constexpr DEAL_II_HOST Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i)
constexpr DEAL_II_HOST 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:174
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_CONSTEXPR
Definition: config.h:187
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
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 DEAL_II_HOST TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
constexpr DEAL_II_HOST 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:939
static const unsigned int invalid_unsigned_int
Definition: types.h:213
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
#define DEAL_II_HOST
Definition: numbers.h:47
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition: numbers.h:702
typename ProductType< Number, OtherNumber >::type value_type
::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:593
constexpr DEAL_II_HOST Number trace(const SymmetricTensor< 2, dim2, Number > &)
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > identity_tensor()
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
constexpr DEAL_II_HOST Number determinant(const SymmetricTensor< 2, dim, Number > &)
SymmetricTensorEigenvectorMethod
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()