Reference documentation for deal.II version Git ee55ce0147 2021-01-21 22:25:19 -0500
\(\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 - 2020 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 <algorithm>
28 #include <array>
29 #include <functional>
30 
32 
33 // Forward declaration
34 #ifndef DOXYGEN
35 template <int rank, int dim, typename Number = double>
36 class SymmetricTensor;
37 #endif
38 
39 template <int dim, typename Number>
42 
43 template <int dim, typename Number>
46 
47 template <int dim, typename Number>
50 
51 template <int dim, typename Number>
54 
55 template <int dim, typename Number>
58 
59 template <int dim2, typename Number>
62 
63 template <int dim, typename Number>
66 
67 template <int dim, typename Number>
70 
71 
72 
73 namespace internal
74 {
75  // Workaround: The following 4 overloads are necessary to be able to
76  // compile the library with Apple Clang 8 and older. We should remove
77  // these overloads again when we bump the minimal required version to
78  // something later than clang-3.6 / Apple Clang 6.3.
79  template <int rank, int dim, typename T, typename U>
80  struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
81  {
82  using type =
83  SymmetricTensor<rank,
84  dim,
85  std::complex<typename ProductType<T, U>::type>>;
86  };
87 
88  template <int rank, int dim, typename T, typename U>
89  struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
90  std::complex<U>>
91  {
92  using type =
93  SymmetricTensor<rank,
94  dim,
95  std::complex<typename ProductType<T, U>::type>>;
96  };
97 
98  template <typename T, int rank, int dim, typename U>
99  struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
100  {
101  using type =
102  SymmetricTensor<rank,
103  dim,
104  std::complex<typename ProductType<T, U>::type>>;
105  };
106 
107  template <int rank, int dim, typename T, typename U>
108  struct ProductTypeImpl<std::complex<T>,
109  SymmetricTensor<rank, dim, std::complex<U>>>
110  {
111  using type =
112  SymmetricTensor<rank,
113  dim,
114  std::complex<typename ProductType<T, U>::type>>;
115  };
116  // end workaround
117 
122  namespace SymmetricTensorImplementation
123  {
128  template <int rank, int dim, typename Number>
129  struct Inverse;
130  } // namespace SymmetricTensorImplementation
131 
136  namespace SymmetricTensorAccessors
137  {
145  merge(const TableIndices<2> &previous_indices,
146  const unsigned int new_index,
147  const unsigned int position)
148  {
149  AssertIndexRange(position, 2);
150 
151  if (position == 0)
152  return {new_index, numbers::invalid_unsigned_int};
153  else
154  return {previous_indices[0], new_index};
155  }
156 
157 
158 
166  merge(const TableIndices<4> &previous_indices,
167  const unsigned int new_index,
168  const unsigned int position)
169  {
170  AssertIndexRange(position, 4);
171 
172  switch (position)
173  {
174  case 0:
175  return {new_index,
178  numbers::invalid_unsigned_int};
179  case 1:
180  return {previous_indices[0],
181  new_index,
183  numbers::invalid_unsigned_int};
184  case 2:
185  return {previous_indices[0],
186  previous_indices[1],
187  new_index,
188  numbers::invalid_unsigned_int};
189  case 3:
190  return {previous_indices[0],
191  previous_indices[1],
192  previous_indices[2],
193  new_index};
194  default:
195  Assert(false, ExcInternalError());
196  return {};
197  }
198  }
199 
200 
207  template <int rank1,
208  int rank2,
209  int dim,
210  typename Number,
211  typename OtherNumber = Number>
213  {
215  using type =
216  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
217  };
218 
219 
226  template <int dim, typename Number, typename OtherNumber>
227  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
228  {
230  };
231 
232 
233 
246  template <int rank, int dim, typename Number>
247  struct StorageType;
248 
252  template <int dim, typename Number>
253  struct StorageType<2, dim, Number>
254  {
259  static const unsigned int n_independent_components =
260  (dim * dim + dim) / 2;
261 
266  };
267 
268 
269 
273  template <int dim, typename Number>
274  struct StorageType<4, dim, Number>
275  {
281  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
282 
286  static const unsigned int n_independent_components =
287  (n_rank2_components *
289 
297  };
298 
299 
300 
305  template <int rank, int dim, bool constness, typename Number>
307 
314  template <int rank, int dim, typename Number>
315  struct AccessorTypes<rank, dim, true, Number>
316  {
317  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
318 
319  using reference = Number;
320  };
321 
328  template <int rank, int dim, typename Number>
329  struct AccessorTypes<rank, dim, false, Number>
330  {
332 
333  using reference = Number &;
334  };
335 
336 
369  template <int rank, int dim, bool constness, int P, typename Number>
370  class Accessor
371  {
372  public:
376  using reference =
378  using tensor_type =
380 
381  private:
400  constexpr Accessor(tensor_type & tensor,
401  const TableIndices<rank> &previous_indices);
402 
406  constexpr DEAL_II_ALWAYS_INLINE
407  Accessor(const Accessor &) = default;
408 
409  public:
413  DEAL_II_CONSTEXPR Accessor<rank, dim, constness, P - 1, Number>
414  operator[](const unsigned int i);
415 
419  constexpr Accessor<rank, dim, constness, P - 1, Number>
420  operator[](const unsigned int i) const;
421 
422  private:
428 
429  // Declare some other classes as friends. Make sure to work around bugs
430  // in some compilers:
431  template <int, int, typename>
432  friend class ::SymmetricTensor;
433  template <int, int, bool, int, typename>
434  friend class Accessor;
435  friend class ::SymmetricTensor<rank, dim, Number>;
436  friend class Accessor<rank, dim, constness, P + 1, Number>;
437  };
438 
439 
440 
448  template <int rank, int dim, bool constness, typename Number>
449  class Accessor<rank, dim, constness, 1, Number>
450  {
451  public:
455  using reference =
457  using tensor_type =
459 
460  private:
482  constexpr Accessor(tensor_type & tensor,
483  const TableIndices<rank> &previous_indices);
484 
488  constexpr DEAL_II_ALWAYS_INLINE
489  Accessor(const Accessor &) = default;
490 
491  public:
495  DEAL_II_CONSTEXPR reference operator[](const unsigned int);
496 
500  constexpr reference operator[](const unsigned int) const;
501 
502  private:
508 
509  // Declare some other classes as friends. Make sure to work around bugs
510  // in some compilers:
511  template <int, int, typename>
512  friend class ::SymmetricTensor;
513  template <int, int, bool, int, typename>
515  friend class ::SymmetricTensor<rank, dim, Number>;
516  friend class SymmetricTensorAccessors::
517  Accessor<rank, dim, constness, 2, Number>;
518  };
519  } // namespace SymmetricTensorAccessors
520 } // namespace internal
521 
522 
523 
596 template <int rank_, int dim, typename Number>
598 {
599 public:
600  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
601 
610  static const unsigned int dimension = dim;
611 
615  static const unsigned int rank = rank_;
616 
622  static constexpr unsigned int n_independent_components =
624  n_independent_components;
625 
629  constexpr DEAL_II_ALWAYS_INLINE
630  SymmetricTensor() = default;
631 
645  template <typename OtherNumber>
646  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
647 
664  SymmetricTensor(const Number (&array)[n_independent_components]);
665 
671  template <typename OtherNumber>
672  constexpr explicit SymmetricTensor(
673  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
674 
678  Number *
679  begin_raw();
680 
684  const Number *
685  begin_raw() const;
686 
690  Number *
691  end_raw();
692 
697  const Number *
698  end_raw() const;
699 
706  template <typename OtherNumber>
708  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
709 
717  operator=(const Number &d);
718 
723  constexpr operator Tensor<rank_, dim, Number>() const;
724 
728  constexpr bool
729  operator==(const SymmetricTensor &) const;
730 
734  constexpr bool
735  operator!=(const SymmetricTensor &) const;
736 
740  template <typename OtherNumber>
742  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
743 
747  template <typename OtherNumber>
749  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
750 
755  template <typename OtherNumber>
757  operator*=(const OtherNumber &factor);
758 
762  template <typename OtherNumber>
764  operator/=(const OtherNumber &factor);
765 
770  operator-() const;
771 
798  template <typename OtherNumber>
802 
807  template <typename OtherNumber>
811 
815  DEAL_II_CONSTEXPR Number &
816  operator()(const TableIndices<rank_> &indices);
817 
821  DEAL_II_CONSTEXPR const Number &
822  operator()(const TableIndices<rank_> &indices) const;
823 
828  constexpr internal::SymmetricTensorAccessors::
829  Accessor<rank_, dim, true, rank_ - 1, Number>
830  operator[](const unsigned int row) const;
831 
836  DEAL_II_CONSTEXPR internal::SymmetricTensorAccessors::
837  Accessor<rank_, dim, false, rank_ - 1, Number>
838  operator[](const unsigned int row);
839 
845  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
846 
852  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
853 
860  DEAL_II_CONSTEXPR const Number &
861  access_raw_entry(const unsigned int unrolled_index) const;
862 
869  DEAL_II_CONSTEXPR Number &
870  access_raw_entry(const unsigned int unrolled_index);
871 
882  norm() const;
883 
891  static constexpr unsigned int
892  component_to_unrolled_index(const TableIndices<rank_> &indices);
893 
899  static constexpr TableIndices<rank_>
900  unrolled_to_component_indices(const unsigned int i);
901 
914  DEAL_II_CONSTEXPR void
915  clear();
916 
921  static constexpr std::size_t
923 
929  template <class Archive>
930  void
931  serialize(Archive &ar, const unsigned int version);
932 
933 private:
937  using base_tensor_descriptor =
939 
943  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
944 
949 
950  // Make all other symmetric tensors friends.
951  template <int, int, typename>
952  friend class SymmetricTensor;
953 
954  // Make a few more functions friends.
955  template <int dim2, typename Number2>
956  friend DEAL_II_CONSTEXPR Number2
958 
959  template <int dim2, typename Number2>
960  friend DEAL_II_CONSTEXPR Number2
962 
963  template <int dim2, typename Number2>
966 
967  template <int dim2, typename Number2>
970 
971  template <int dim2, typename Number2>
973  deviator_tensor();
974 
975  template <int dim2, typename Number2>
977  identity_tensor();
978 
979 
980  // Make a few helper classes friends as well.
982  Inverse<2, dim, Number>;
983 
985  Inverse<4, dim, Number>;
986 };
987 
988 
989 
990 // ------------------------- inline functions ------------------------
991 
992 #ifndef DOXYGEN
993 
994 // provide declarations for static members
995 template <int rank, int dim, typename Number>
996 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
997 
998 template <int rank_, int dim, typename Number>
999 constexpr unsigned int
1000  SymmetricTensor<rank_, dim, Number>::n_independent_components;
1001 
1002 namespace internal
1003 {
1004  namespace SymmetricTensorAccessors
1005  {
1006  template <int rank_, int dim, bool constness, int P, typename Number>
1007  constexpr DEAL_II_ALWAYS_INLINE
1008  Accessor<rank_, dim, constness, P, Number>::Accessor(
1009  tensor_type & tensor,
1010  const TableIndices<rank_> &previous_indices)
1011  : tensor(tensor)
1012  , previous_indices(previous_indices)
1013  {}
1014 
1015 
1016 
1017  template <int rank_, int dim, bool constness, int P, typename Number>
1019  Accessor<rank_, dim, constness, P - 1, Number>
1020  Accessor<rank_, dim, constness, P, Number>::
1021  operator[](const unsigned int i)
1022  {
1023  return Accessor<rank_, dim, constness, P - 1, Number>(
1024  tensor, merge(previous_indices, i, rank_ - P));
1025  }
1026 
1027 
1028 
1029  template <int rank_, int dim, bool constness, int P, typename Number>
1030  constexpr DEAL_II_ALWAYS_INLINE
1031  Accessor<rank_, dim, constness, P - 1, Number>
1032  Accessor<rank_, dim, constness, P, Number>::
1033  operator[](const unsigned int i) const
1034  {
1035  return Accessor<rank_, dim, constness, P - 1, Number>(
1036  tensor, merge(previous_indices, i, rank_ - P));
1037  }
1038 
1039 
1040 
1041  template <int rank_, int dim, bool constness, typename Number>
1042  constexpr DEAL_II_ALWAYS_INLINE
1043  Accessor<rank_, dim, constness, 1, Number>::Accessor(
1044  tensor_type & tensor,
1045  const TableIndices<rank_> &previous_indices)
1046  : tensor(tensor)
1047  , previous_indices(previous_indices)
1048  {}
1049 
1050 
1051 
1052  template <int rank_, int dim, bool constness, typename Number>
1054  typename Accessor<rank_, dim, constness, 1, Number>::reference
1055  Accessor<rank_, dim, constness, 1, Number>::
1056  operator[](const unsigned int i)
1057  {
1058  return tensor(merge(previous_indices, i, rank_ - 1));
1059  }
1060 
1061 
1062  template <int rank_, int dim, bool constness, typename Number>
1063  constexpr DEAL_II_ALWAYS_INLINE
1064  typename Accessor<rank_, dim, constness, 1, Number>::reference
1065  Accessor<rank_, dim, constness, 1, Number>::
1066  operator[](const unsigned int i) const
1067  {
1068  return tensor(merge(previous_indices, i, rank_ - 1));
1069  }
1070  } // namespace SymmetricTensorAccessors
1071 } // namespace internal
1072 
1073 
1074 
1075 template <int rank_, int dim, typename Number>
1076 template <typename OtherNumber>
1077 inline DEAL_II_ALWAYS_INLINE
1079  const Tensor<2, dim, OtherNumber> &t)
1080 {
1081  static_assert(rank == 2, "This function is only implemented for rank==2");
1082  for (unsigned int d = 0; d < dim; ++d)
1083  for (unsigned int e = 0; e < d; ++e)
1084  Assert(t[d][e] == t[e][d],
1085  ExcMessage("The incoming Tensor must be exactly symmetric."));
1086 
1087  for (unsigned int d = 0; d < dim; ++d)
1088  data[d] = t[d][d];
1089 
1090  for (unsigned int d = 0, c = 0; d < dim; ++d)
1091  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1092  data[dim + c] = t[d][e];
1093 }
1094 
1095 
1096 
1097 template <int rank_, int dim, typename Number>
1098 template <typename OtherNumber>
1099 constexpr DEAL_II_ALWAYS_INLINE
1101  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1102  : data(initializer.data)
1103 {}
1104 
1105 
1106 
1107 template <int rank_, int dim, typename Number>
1110  const Number (&array)[n_independent_components])
1111  : data(
1112  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1113 {
1114  // ensure that the reinterpret_cast above actually works
1115  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1116  ExcInternalError());
1117 }
1118 
1119 
1120 
1121 template <int rank_, int dim, typename Number>
1122 template <typename OtherNumber>
1127 {
1128  data = t.data;
1129  return *this;
1130 }
1131 
1132 
1133 
1134 template <int rank_, int dim, typename Number>
1138 {
1140  ExcMessage("Only assignment with zero is allowed"));
1141  (void)d;
1142 
1144 
1145  return *this;
1146 }
1147 
1148 
1149 namespace internal
1150 {
1151  namespace SymmetricTensorImplementation
1152  {
1153  template <int dim, typename Number>
1154  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1155  ::Tensor<2, dim, Number>
1156  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1157  {
1159 
1160  // diagonal entries are stored first
1161  for (unsigned int d = 0; d < dim; ++d)
1162  t[d][d] = s.access_raw_entry(d);
1163 
1164  // off-diagonal entries come next, row by row
1165  for (unsigned int d = 0, c = 0; d < dim; ++d)
1166  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1167  {
1168  t[d][e] = s.access_raw_entry(dim + c);
1169  t[e][d] = s.access_raw_entry(dim + c);
1170  }
1171  return t;
1172  }
1173 
1174 
1175  template <int dim, typename Number>
1176  DEAL_II_CONSTEXPR ::Tensor<4, dim, Number>
1177  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1178  {
1179  // utilize the symmetry properties of SymmetricTensor<4,dim>
1180  // discussed in the class documentation to avoid accessing all
1181  // independent elements of the input tensor more than once
1183 
1184  for (unsigned int i = 0; i < dim; ++i)
1185  for (unsigned int j = i; j < dim; ++j)
1186  for (unsigned int k = 0; k < dim; ++k)
1187  for (unsigned int l = k; l < dim; ++l)
1188  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1189  t[TableIndices<4>(j, i, k, l)] =
1190  t[TableIndices<4>(j, i, l, k)] =
1191  st[TableIndices<4>(i, j, k, l)];
1192 
1193  return t;
1194  }
1195 
1196 
1197  template <typename Number>
1198  struct Inverse<2, 1, Number>
1199  {
1200  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1201  ::SymmetricTensor<2, 1, Number>
1202  value(const ::SymmetricTensor<2, 1, Number> &t)
1203  {
1205 
1206  tmp[0][0] = 1.0 / t[0][0];
1207 
1208  return tmp;
1209  }
1210  };
1211 
1212 
1213  template <typename Number>
1214  struct Inverse<2, 2, Number>
1215  {
1216  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1217  ::SymmetricTensor<2, 2, Number>
1218  value(const ::SymmetricTensor<2, 2, Number> &t)
1219  {
1221 
1222  // Sympy result: ([
1223  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1224  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1225  const TableIndices<2> idx_00(0, 0);
1226  const TableIndices<2> idx_01(0, 1);
1227  const TableIndices<2> idx_11(1, 1);
1228  const Number inv_det_t =
1229  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1230  tmp[idx_00] = t[idx_11];
1231  tmp[idx_01] = -t[idx_01];
1232  tmp[idx_11] = t[idx_00];
1233  tmp *= inv_det_t;
1234 
1235  return tmp;
1236  }
1237  };
1238 
1239 
1240  template <typename Number>
1241  struct Inverse<2, 3, Number>
1242  {
1243  DEAL_II_CONSTEXPR static ::SymmetricTensor<2, 3, Number>
1244  value(const ::SymmetricTensor<2, 3, Number> &t)
1245  {
1247 
1248  // Sympy result: ([
1249  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1250  // 2*t01*t02*t12 - t02**2*t11),
1251  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1252  // 2*t01*t02*t12 - t02**2*t11),
1253  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1254  // 2*t01*t02*t12 - t02**2*t11)],
1255  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1256  // 2*t01*t02*t12 - t02**2*t11),
1257  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1258  // 2*t01*t02*t12 - t02**2*t11),
1259  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1260  // 2*t01*t02*t12 + t02**2*t11)],
1261  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1262  // 2*t01*t02*t12 - t02**2*t11),
1263  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1264  // 2*t01*t02*t12 + t02**2*t11),
1265  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1266  // 2*t01*t02*t12 + t02**2*t11)] ])
1267  //
1268  // =
1269  //
1270  // [ (t11*t22 - t12**2)/det_t,
1271  // (-t01*t22 + t02*t12)/det_t,
1272  // (t01*t12 - t02*t11)/det_t],
1273  // [ (-t01*t22 + t02*t12)/det_t,
1274  // (t00*t22 - t02**2)/det_t,
1275  // (-t00*t12 + t01*t02)/det_t],
1276  // [ (t01*t12 - t02*t11)/det_t,
1277  // (-t00*t12 + t01*t02)/det_t,
1278  // (t00*t11 - t01**2)/det_t] ])
1279  //
1280  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1281  // 2*t01*t02*t12 - t02**2*t11)
1282  const TableIndices<2> idx_00(0, 0);
1283  const TableIndices<2> idx_01(0, 1);
1284  const TableIndices<2> idx_02(0, 2);
1285  const TableIndices<2> idx_11(1, 1);
1286  const TableIndices<2> idx_12(1, 2);
1287  const TableIndices<2> idx_22(2, 2);
1288  const Number inv_det_t =
1289  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1290  t[idx_00] * t[idx_12] * t[idx_12] -
1291  t[idx_01] * t[idx_01] * t[idx_22] +
1292  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1293  t[idx_02] * t[idx_02] * t[idx_11]);
1294  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1295  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1296  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1297  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1298  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1299  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1300  tmp *= inv_det_t;
1301 
1302  return tmp;
1303  }
1304  };
1305 
1306 
1307  template <typename Number>
1308  struct Inverse<4, 1, Number>
1309  {
1310  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 1, Number>
1311  value(const ::SymmetricTensor<4, 1, Number> &t)
1312  {
1314  tmp.data[0][0] = 1.0 / t.data[0][0];
1315  return tmp;
1316  }
1317  };
1318 
1319 
1320  template <typename Number>
1321  struct Inverse<4, 2, Number>
1322  {
1323  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 2, Number>
1324  value(const ::SymmetricTensor<4, 2, Number> &t)
1325  {
1327 
1328  // Inverting this tensor is a little more complicated than necessary,
1329  // since we store the data of 't' as a 3x3 matrix t.data, but the
1330  // product between a rank-4 and a rank-2 tensor is really not the
1331  // product between this matrix and the 3-vector of a rhs, but rather
1332  //
1333  // B.vec = t.data * mult * A.vec
1334  //
1335  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1336  // capture the fact that we need to add up both the c_ij12*a_12 and the
1337  // c_ij21*a_21 terms.
1338  //
1339  // In addition, in this scheme, the identity tensor has the matrix
1340  // representation mult^-1.
1341  //
1342  // The inverse of 't' therefore has the matrix representation
1343  //
1344  // inv.data = mult^-1 * t.data^-1 * mult^-1
1345  //
1346  // in order to compute it, let's first compute the inverse of t.data and
1347  // put it into tmp.data; at the end of the function we then scale the
1348  // last row and column of the inverse by 1/2, corresponding to the left
1349  // and right multiplication with mult^-1.
1350  const Number t4 = t.data[0][0] * t.data[1][1],
1351  t6 = t.data[0][0] * t.data[1][2],
1352  t8 = t.data[0][1] * t.data[1][0],
1353  t00 = t.data[0][2] * t.data[1][0],
1354  t01 = t.data[0][1] * t.data[2][0],
1355  t04 = t.data[0][2] * t.data[2][0],
1356  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1357  t8 * t.data[2][2] + t00 * t.data[2][1] +
1358  t01 * t.data[1][2] - t04 * t.data[1][1]);
1359  tmp.data[0][0] =
1360  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1361  tmp.data[0][1] =
1362  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1363  tmp.data[0][2] =
1364  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1365  tmp.data[1][0] =
1366  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1367  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1368  tmp.data[1][2] = -(t6 - t00) * t07;
1369  tmp.data[2][0] =
1370  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1371  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1372  tmp.data[2][2] = (t4 - t8) * t07;
1373 
1374  // scale last row and column as mentioned
1375  // above
1376  tmp.data[2][0] /= 2;
1377  tmp.data[2][1] /= 2;
1378  tmp.data[0][2] /= 2;
1379  tmp.data[1][2] /= 2;
1380  tmp.data[2][2] /= 4;
1381 
1382  return tmp;
1383  }
1384  };
1385 
1386 
1387  template <typename Number>
1388  struct Inverse<4, 3, Number>
1389  {
1390  static ::SymmetricTensor<4, 3, Number>
1391  value(const ::SymmetricTensor<4, 3, Number> &t)
1392  {
1394 
1395  // This function follows the exact same scheme as the 2d case, except
1396  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1397  // Instead, we use the Gauss-Jordan algorithm implemented for
1398  // FullMatrix. For historical reasons the following code is copied from
1399  // there, with the tangential benefit that we do not need to copy the
1400  // tensor entries to and from the FullMatrix.
1401  const unsigned int N = 6;
1402 
1403  // First get an estimate of the size of the elements of this matrix,
1404  // for later checks whether the pivot element is large enough, or
1405  // whether we have to fear that the matrix is not regular.
1406  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1407  for (unsigned int i = 0; i < N; ++i)
1408  diagonal_sum += std::fabs(tmp.data[i][i]);
1409  const Number typical_diagonal_element =
1410  diagonal_sum / static_cast<double>(N);
1411  (void)typical_diagonal_element;
1412 
1413  unsigned int p[N];
1414  for (unsigned int i = 0; i < N; ++i)
1415  p[i] = i;
1416 
1417  for (unsigned int j = 0; j < N; ++j)
1418  {
1419  // Pivot search: search that part of the line on and right of the
1420  // diagonal for the largest element.
1421  Number max = std::fabs(tmp.data[j][j]);
1422  unsigned int r = j;
1423  for (unsigned int i = j + 1; i < N; ++i)
1424  if (std::fabs(tmp.data[i][j]) > max)
1425  {
1426  max = std::fabs(tmp.data[i][j]);
1427  r = i;
1428  }
1429 
1430  // Check whether the pivot is too small
1431  Assert(max > 1.e-16 * typical_diagonal_element,
1432  ExcMessage("This tensor seems to be noninvertible"));
1433 
1434  // Row interchange
1435  if (r > j)
1436  {
1437  for (unsigned int k = 0; k < N; ++k)
1438  std::swap(tmp.data[j][k], tmp.data[r][k]);
1439 
1440  std::swap(p[j], p[r]);
1441  }
1442 
1443  // Transformation
1444  const Number hr = 1. / tmp.data[j][j];
1445  tmp.data[j][j] = hr;
1446  for (unsigned int k = 0; k < N; ++k)
1447  {
1448  if (k == j)
1449  continue;
1450  for (unsigned int i = 0; i < N; ++i)
1451  {
1452  if (i == j)
1453  continue;
1454  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1455  }
1456  }
1457  for (unsigned int i = 0; i < N; ++i)
1458  {
1459  tmp.data[i][j] *= hr;
1460  tmp.data[j][i] *= -hr;
1461  }
1462  tmp.data[j][j] = hr;
1463  }
1464 
1465  // Column interchange
1466  Number hv[N];
1467  for (unsigned int i = 0; i < N; ++i)
1468  {
1469  for (unsigned int k = 0; k < N; ++k)
1470  hv[p[k]] = tmp.data[i][k];
1471  for (unsigned int k = 0; k < N; ++k)
1472  tmp.data[i][k] = hv[k];
1473  }
1474 
1475  // Scale rows and columns. The mult matrix
1476  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1477  for (unsigned int i = 3; i < 6; ++i)
1478  for (unsigned int j = 0; j < 3; ++j)
1479  tmp.data[i][j] /= 2;
1480 
1481  for (unsigned int i = 0; i < 3; ++i)
1482  for (unsigned int j = 3; j < 6; ++j)
1483  tmp.data[i][j] /= 2;
1484 
1485  for (unsigned int i = 3; i < 6; ++i)
1486  for (unsigned int j = 3; j < 6; ++j)
1487  tmp.data[i][j] /= 4;
1488 
1489  return tmp;
1490  }
1491  };
1492 
1493  } // namespace SymmetricTensorImplementation
1494 } // namespace internal
1495 
1496 
1497 
1498 template <int rank_, int dim, typename Number>
1500  operator Tensor<rank_, dim, Number>() const
1501 {
1502  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1503 }
1504 
1505 
1506 
1507 template <int rank_, int dim, typename Number>
1508 constexpr bool
1511 {
1512  return data == t.data;
1513 }
1514 
1515 
1516 
1517 template <int rank_, int dim, typename Number>
1518 constexpr bool
1521 {
1522  return data != t.data;
1523 }
1524 
1525 
1526 
1527 template <int rank_, int dim, typename Number>
1528 template <typename OtherNumber>
1533 {
1534  data += t.data;
1535  return *this;
1536 }
1537 
1538 
1539 
1540 template <int rank_, int dim, typename Number>
1541 template <typename OtherNumber>
1546 {
1547  data -= t.data;
1548  return *this;
1549 }
1550 
1551 
1552 
1553 template <int rank_, int dim, typename Number>
1554 template <typename OtherNumber>
1558 {
1559  data *= d;
1560  return *this;
1561 }
1562 
1563 
1564 
1565 template <int rank_, int dim, typename Number>
1566 template <typename OtherNumber>
1570 {
1571  data /= d;
1572  return *this;
1573 }
1574 
1575 
1576 
1577 template <int rank_, int dim, typename Number>
1581 {
1582  SymmetricTensor tmp = *this;
1583  tmp.data = -tmp.data;
1584  return tmp;
1585 }
1586 
1587 
1588 
1589 template <int rank_, int dim, typename Number>
1592 {
1593  data.clear();
1594 }
1595 
1596 
1597 
1598 template <int rank_, int dim, typename Number>
1599 constexpr std::size_t
1601 {
1602  // all memory consists of statically allocated memory of the current
1603  // object, no pointers
1604  return sizeof(SymmetricTensor<rank_, dim, Number>);
1605 }
1606 
1607 
1608 
1609 namespace internal
1610 {
1611  template <int dim, typename Number, typename OtherNumber = Number>
1615  perform_double_contraction(
1616  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1617  base_tensor_type &data,
1618  const typename SymmetricTensorAccessors::
1619  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1620  {
1621  using result_type = typename SymmetricTensorAccessors::
1623 
1624  switch (dim)
1625  {
1626  case 1:
1627  return data[0] * sdata[0];
1628  default:
1629  // Start with the non-diagonal part to avoid some multiplications by
1630  // 2.
1631 
1632  result_type sum = data[dim] * sdata[dim];
1633  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1634  sum += data[d] * sdata[d];
1635  sum += sum; // sum = sum * 2.;
1636 
1637  // Now add the contributions from the diagonal
1638  for (unsigned int d = 0; d < dim; ++d)
1639  sum += data[d] * sdata[d];
1640  return sum;
1641  }
1642  }
1643 
1644 
1645 
1646  template <int dim, typename Number, typename OtherNumber = Number>
1650  perform_double_contraction(
1651  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1652  base_tensor_type &data,
1653  const typename SymmetricTensorAccessors::
1654  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1655  {
1656  using result_type = typename SymmetricTensorAccessors::
1658  using value_type = typename SymmetricTensorAccessors::
1660 
1661  const unsigned int data_dim = SymmetricTensorAccessors::
1662  StorageType<2, dim, value_type>::n_independent_components;
1663  value_type tmp[data_dim]{};
1664  for (unsigned int i = 0; i < data_dim; ++i)
1665  tmp[i] =
1666  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1667  return result_type(tmp);
1668  }
1669 
1670 
1671 
1672  template <int dim, typename Number, typename OtherNumber = Number>
1674  typename SymmetricTensorAccessors::StorageType<
1675  2,
1676  dim,
1679  base_tensor_type
1680  perform_double_contraction(
1681  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1682  base_tensor_type &data,
1683  const typename SymmetricTensorAccessors::
1684  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1685  {
1686  using value_type = typename SymmetricTensorAccessors::
1688  using base_tensor_type = typename SymmetricTensorAccessors::
1689  StorageType<2, dim, value_type>::base_tensor_type;
1690 
1691  base_tensor_type tmp;
1692  for (unsigned int i = 0; i < tmp.dimension; ++i)
1693  {
1694  // Start with the non-diagonal part
1695  value_type sum = data[dim] * sdata[dim][i];
1696  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1697  sum += data[d] * sdata[d][i];
1698  sum += sum; // sum = sum * 2.;
1699 
1700  // Now add the contributions from the diagonal
1701  for (unsigned int d = 0; d < dim; ++d)
1702  sum += data[d] * sdata[d][i];
1703  tmp[i] = sum;
1704  }
1705  return tmp;
1706  }
1707 
1708 
1709 
1710  template <int dim, typename Number, typename OtherNumber = Number>
1712  typename SymmetricTensorAccessors::StorageType<
1713  4,
1714  dim,
1717  base_tensor_type
1718  perform_double_contraction(
1719  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1720  base_tensor_type &data,
1721  const typename SymmetricTensorAccessors::
1722  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1723  {
1724  using value_type = typename SymmetricTensorAccessors::
1726  using base_tensor_type = typename SymmetricTensorAccessors::
1727  StorageType<4, dim, value_type>::base_tensor_type;
1728 
1729  const unsigned int data_dim = SymmetricTensorAccessors::
1730  StorageType<2, dim, value_type>::n_independent_components;
1731  base_tensor_type tmp;
1732  for (unsigned int i = 0; i < data_dim; ++i)
1733  for (unsigned int j = 0; j < data_dim; ++j)
1734  {
1735  // Start with the non-diagonal part
1736  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1737  tmp[i][j] += data[i][d] * sdata[d][j];
1738  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1739 
1740  // Now add the contributions from the diagonal
1741  for (unsigned int d = 0; d < dim; ++d)
1742  tmp[i][j] += data[i][d] * sdata[d][j];
1743  }
1744  return tmp;
1745  }
1746 
1747 } // end of namespace internal
1748 
1749 
1750 
1751 template <int rank_, int dim, typename Number>
1752 template <typename OtherNumber>
1758 {
1759  // need to have two different function calls
1760  // because a scalar and rank-2 tensor are not
1761  // the same data type (see internal function
1762  // above)
1763  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1764  s.data);
1765 }
1766 
1767 
1768 
1769 template <int rank_, int dim, typename Number>
1770 template <typename OtherNumber>
1775 {
1778  tmp.data =
1779  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1780  s.data);
1781  return tmp;
1782 }
1783 
1784 
1785 
1786 // internal namespace to switch between the
1787 // access of different tensors. There used to
1788 // be explicit instantiations before for
1789 // different ranks and dimensions, but since
1790 // we now allow for templates on the data
1791 // type, and since we cannot partially
1792 // specialize the implementation, this got
1793 // into a separate namespace
1794 namespace internal
1795 {
1796  // The variables within this struct will be referenced in the next functions.
1797  // It is a workaround that allows returning a reference to a static variable
1798  // while allowing constexpr evaluation of the function.
1799  // It has to be defined outside the function because constexpr functions
1800  // cannot define static variables.
1801  // A similar struct has also been defined in tensor.h
1802  template <typename Type>
1803  struct Uninitialized
1804  {
1805  static Type value;
1806  };
1807 
1808  template <typename Type>
1810 
1811  template <int dim, typename Number>
1813  symmetric_tensor_access(const TableIndices<2> &indices,
1814  typename SymmetricTensorAccessors::
1815  StorageType<2, dim, Number>::base_tensor_type &data)
1816  {
1817  // 1d is very simple and done first
1818  if (dim == 1)
1819  return data[0];
1820 
1821  // first treat the main diagonal elements, which are stored consecutively
1822  // at the beginning
1823  if (indices[0] == indices[1])
1824  return data[indices[0]];
1825 
1826  // the rest is messier and requires a few switches.
1827  switch (dim)
1828  {
1829  case 2:
1830  // at least for the 2x2 case it is reasonably simple
1831  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1832  ((indices[0] == 0) && (indices[1] == 1)),
1833  ExcInternalError());
1834  return data[2];
1835 
1836  default:
1837  // to do the rest, sort our indices before comparing
1838  {
1839  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1840  std::max(indices[0], indices[1]));
1841  for (unsigned int d = 0, c = 0; d < dim; ++d)
1842  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1843  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1844  return data[dim + c];
1845  Assert(false, ExcInternalError());
1846  }
1847  }
1848 
1849  // The code should never reach there.
1850  // Returns a dummy reference to a dummy variable just to make the
1851  // compiler happy.
1853  }
1854 
1855 
1856 
1857  template <int dim, typename Number>
1858  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1859  symmetric_tensor_access(const TableIndices<2> &indices,
1860  const typename SymmetricTensorAccessors::
1861  StorageType<2, dim, Number>::base_tensor_type &data)
1862  {
1863  // 1d is very simple and done first
1864  if (dim == 1)
1865  return data[0];
1866 
1867  // first treat the main diagonal elements, which are stored consecutively
1868  // at the beginning
1869  if (indices[0] == indices[1])
1870  return data[indices[0]];
1871 
1872  // the rest is messier and requires a few switches.
1873  switch (dim)
1874  {
1875  case 2:
1876  // at least for the 2x2 case it is reasonably simple
1877  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1878  ((indices[0] == 0) && (indices[1] == 1)),
1879  ExcInternalError());
1880  return data[2];
1881 
1882  default:
1883  // to do the rest, sort our indices before comparing
1884  {
1885  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1886  std::max(indices[0], indices[1]));
1887  for (unsigned int d = 0, c = 0; d < dim; ++d)
1888  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1889  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1890  return data[dim + c];
1891  Assert(false, ExcInternalError());
1892  }
1893  }
1894 
1895  // The code should never reach there.
1896  // Returns a dummy reference to a dummy variable just to make the
1897  // compiler happy.
1899  }
1900 
1901 
1902 
1903  template <int dim, typename Number>
1904  DEAL_II_CONSTEXPR inline Number &
1905  symmetric_tensor_access(const TableIndices<4> &indices,
1906  typename SymmetricTensorAccessors::
1907  StorageType<4, dim, Number>::base_tensor_type &data)
1908  {
1909  switch (dim)
1910  {
1911  case 1:
1912  return data[0][0];
1913 
1914  case 2:
1915  // each entry of the tensor can be thought of as an entry in a
1916  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1917  // rank-2 tensors. this is the format in which we store rank-4
1918  // tensors. determine which position the present entry is
1919  // stored in
1920  {
1921  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1922  return data[base_index[indices[0]][indices[1]]]
1923  [base_index[indices[2]][indices[3]]];
1924  }
1925  case 3:
1926  // each entry of the tensor can be thought of as an entry in a
1927  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1928  // rank-2 tensors. this is the format in which we store rank-4
1929  // tensors. determine which position the present entry is
1930  // stored in
1931  {
1932  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1933  {3, 1, 5},
1934  {4, 5, 2}};
1935  return data[base_index[indices[0]][indices[1]]]
1936  [base_index[indices[2]][indices[3]]];
1937  }
1938 
1939  default:
1940  Assert(false, ExcNotImplemented());
1941  }
1942 
1943  // The code should never reach there.
1944  // Returns a dummy reference to a dummy variable just to make the
1945  // compiler happy.
1947  }
1948 
1949 
1950  template <int dim, typename Number>
1951  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1952  symmetric_tensor_access(const TableIndices<4> &indices,
1953  const typename SymmetricTensorAccessors::
1954  StorageType<4, dim, Number>::base_tensor_type &data)
1955  {
1956  switch (dim)
1957  {
1958  case 1:
1959  return data[0][0];
1960 
1961  case 2:
1962  // each entry of the tensor can be thought of as an entry in a
1963  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1964  // rank-2 tensors. this is the format in which we store rank-4
1965  // tensors. determine which position the present entry is
1966  // stored in
1967  {
1968  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1969  return data[base_index[indices[0]][indices[1]]]
1970  [base_index[indices[2]][indices[3]]];
1971  }
1972  case 3:
1973  // each entry of the tensor can be thought of as an entry in a
1974  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1975  // rank-2 tensors. this is the format in which we store rank-4
1976  // tensors. determine which position the present entry is
1977  // stored in
1978  {
1979  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1980  {3, 1, 5},
1981  {4, 5, 2}};
1982  return data[base_index[indices[0]][indices[1]]]
1983  [base_index[indices[2]][indices[3]]];
1984  }
1985 
1986  default:
1987  Assert(false, ExcNotImplemented());
1988  }
1989 
1990  // The code should never reach there.
1991  // Returns a dummy reference to a dummy variable just to make the
1992  // compiler happy.
1994  }
1995 
1996 } // end of namespace internal
1997 
1998 
1999 
2000 template <int rank_, int dim, typename Number>
2003  operator()(const TableIndices<rank_> &indices)
2004 {
2005  for (unsigned int r = 0; r < rank; ++r)
2006  AssertIndexRange(indices[r], dimension);
2007  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2008 }
2009 
2010 
2011 
2012 template <int rank_, int dim, typename Number>
2013 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
2015  operator()(const TableIndices<rank_> &indices) const
2016 {
2017  for (unsigned int r = 0; r < rank; ++r)
2018  AssertIndexRange(indices[r], dimension);
2019  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2020 }
2021 
2022 
2023 
2024 namespace internal
2025 {
2026  namespace SymmetricTensorImplementation
2027  {
2028  template <int rank_>
2029  constexpr TableIndices<rank_>
2030  get_partially_filled_indices(const unsigned int row,
2031  const std::integral_constant<int, 2> &)
2032  {
2034  }
2035 
2036 
2037  template <int rank_>
2038  constexpr TableIndices<rank_>
2039  get_partially_filled_indices(const unsigned int row,
2040  const std::integral_constant<int, 4> &)
2041  {
2042  return TableIndices<rank_>(row,
2046  }
2047  } // namespace SymmetricTensorImplementation
2048 } // namespace internal
2049 
2050 
2051 template <int rank_, int dim, typename Number>
2052 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2053  Accessor<rank_, dim, true, rank_ - 1, Number>
2055  operator[](const unsigned int row) const
2056 {
2057  return internal::SymmetricTensorAccessors::
2058  Accessor<rank_, dim, true, rank_ - 1, Number>(
2059  *this,
2060  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2061  rank_>(row, std::integral_constant<int, rank_>()));
2062 }
2063 
2064 
2065 
2066 template <int rank_, int dim, typename Number>
2067 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE internal::
2068  SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2069  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2070 {
2071  return internal::SymmetricTensorAccessors::
2072  Accessor<rank_, dim, false, rank_ - 1, Number>(
2073  *this,
2074  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2075  rank_>(row, std::integral_constant<int, rank_>()));
2076 }
2077 
2078 
2079 
2080 template <int rank_, int dim, typename Number>
2081 constexpr DEAL_II_ALWAYS_INLINE const Number &
2083  operator[](const TableIndices<rank_> &indices) const
2084 {
2085  return operator()(indices);
2086 }
2087 
2088 
2089 
2090 template <int rank_, int dim, typename Number>
2093  operator[](const TableIndices<rank_> &indices)
2094 {
2095  return operator()(indices);
2096 }
2097 
2098 
2099 
2100 template <int rank_, int dim, typename Number>
2101 inline Number *
2103 {
2104  return std::addressof(this->access_raw_entry(0));
2105 }
2106 
2107 
2108 
2109 template <int rank_, int dim, typename Number>
2110 inline const Number *
2112 {
2113  return std::addressof(this->access_raw_entry(0));
2114 }
2115 
2116 
2117 
2118 template <int rank_, int dim, typename Number>
2119 inline Number *
2121 {
2122  return begin_raw() + n_independent_components;
2123 }
2124 
2125 
2126 
2127 template <int rank_, int dim, typename Number>
2128 inline const Number *
2130 {
2131  return begin_raw() + n_independent_components;
2132 }
2133 
2134 
2135 
2136 namespace internal
2137 {
2138  namespace SymmetricTensorImplementation
2139  {
2140  template <int dim, typename Number>
2141  constexpr unsigned int
2142  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2143  const unsigned int index)
2144  {
2145  return index;
2146  }
2147 
2148 
2149  template <int dim, typename Number>
2150  constexpr ::TableIndices<2>
2151  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2152  const unsigned int index)
2153  {
2156  }
2157 
2158  } // namespace SymmetricTensorImplementation
2159 } // namespace internal
2160 
2161 
2162 
2163 template <int rank_, int dim, typename Number>
2164 DEAL_II_CONSTEXPR inline const Number &
2166  const unsigned int index) const
2167 {
2168  AssertIndexRange(index, n_independent_components);
2169  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2170  index)];
2171 }
2172 
2173 
2174 
2175 template <int rank_, int dim, typename Number>
2176 DEAL_II_CONSTEXPR inline Number &
2178 {
2179  AssertIndexRange(index, n_independent_components);
2180  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2181  index)];
2182 }
2183 
2184 
2185 
2186 namespace internal
2187 {
2188  template <int dim, typename Number>
2190  compute_norm(const typename SymmetricTensorAccessors::
2191  StorageType<2, dim, Number>::base_tensor_type &data)
2192  {
2193  switch (dim)
2194  {
2195  case 1:
2196  return numbers::NumberTraits<Number>::abs(data[0]);
2197 
2198  case 2:
2199  return std::sqrt(
2203 
2204  case 3:
2205  return std::sqrt(
2212 
2213  default:
2214  {
2215  typename numbers::NumberTraits<Number>::real_type return_value =
2217 
2218  for (unsigned int d = 0; d < dim; ++d)
2219  return_value +=
2221  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2222  return_value +=
2224 
2225  return std::sqrt(return_value);
2226  }
2227  }
2228  }
2229 
2230 
2231 
2232  template <int dim, typename Number>
2234  compute_norm(const typename SymmetricTensorAccessors::
2235  StorageType<4, dim, Number>::base_tensor_type &data)
2236  {
2237  switch (dim)
2238  {
2239  case 1:
2240  return numbers::NumberTraits<Number>::abs(data[0][0]);
2241 
2242  default:
2243  {
2244  typename numbers::NumberTraits<Number>::real_type return_value =
2246 
2247  const unsigned int n_independent_components = data.dimension;
2248 
2249  for (unsigned int i = 0; i < dim; ++i)
2250  for (unsigned int j = 0; j < dim; ++j)
2251  return_value +=
2253  for (unsigned int i = 0; i < dim; ++i)
2254  for (unsigned int j = dim; j < n_independent_components; ++j)
2255  return_value +=
2257  for (unsigned int i = dim; i < n_independent_components; ++i)
2258  for (unsigned int j = 0; j < dim; ++j)
2259  return_value +=
2261  for (unsigned int i = dim; i < n_independent_components; ++i)
2262  for (unsigned int j = dim; j < n_independent_components; ++j)
2263  return_value +=
2265 
2266  return std::sqrt(return_value);
2267  }
2268  }
2269  }
2270 
2271 } // end of namespace internal
2272 
2273 
2274 
2275 template <int rank_, int dim, typename Number>
2278 {
2279  return internal::compute_norm<dim, Number>(data);
2280 }
2281 
2282 
2283 
2284 namespace internal
2285 {
2286  namespace SymmetricTensorImplementation
2287  {
2288  // a function to do the unrolling from a set of indices to a
2289  // scalar index into the array in which we store the elements of
2290  // a symmetric tensor
2291  //
2292  // this function is for rank-2 tensors
2293  template <int dim>
2294  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE unsigned int
2295  component_to_unrolled_index(const TableIndices<2> &indices)
2296  {
2297  AssertIndexRange(indices[0], dim);
2298  AssertIndexRange(indices[1], dim);
2299 
2300  switch (dim)
2301  {
2302  case 1:
2303  {
2304  return 0;
2305  }
2306 
2307  case 2:
2308  {
2309  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2310  return table[indices[0]][indices[1]];
2311  }
2312 
2313  case 3:
2314  {
2315  constexpr unsigned int table[3][3] = {{0, 3, 4},
2316  {3, 1, 5},
2317  {4, 5, 2}};
2318  return table[indices[0]][indices[1]];
2319  }
2320 
2321  case 4:
2322  {
2323  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2324  {4, 1, 7, 8},
2325  {5, 7, 2, 9},
2326  {6, 8, 9, 3}};
2327  return table[indices[0]][indices[1]];
2328  }
2329 
2330  default:
2331  // for the remainder, manually figure out the numbering
2332  {
2333  if (indices[0] == indices[1])
2334  return indices[0];
2335 
2336  TableIndices<2> sorted_indices(indices);
2337  sorted_indices.sort();
2338 
2339  for (unsigned int d = 0, c = 0; d < dim; ++d)
2340  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2341  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2342  return dim + c;
2343 
2344  // should never get here:
2345  Assert(false, ExcInternalError());
2346  return 0;
2347  }
2348  }
2349  }
2350 
2351  // a function to do the unrolling from a set of indices to a
2352  // scalar index into the array in which we store the elements of
2353  // a symmetric tensor
2354  //
2355  // this function is for tensors of ranks not already handled
2356  // above
2357  template <int dim, int rank_>
2358  DEAL_II_CONSTEXPR inline unsigned int
2359  component_to_unrolled_index(const TableIndices<rank_> &indices)
2360  {
2361  (void)indices;
2362  Assert(false, ExcNotImplemented());
2364  }
2365  } // namespace SymmetricTensorImplementation
2366 } // namespace internal
2367 
2368 
2369 template <int rank_, int dim, typename Number>
2370 constexpr unsigned int
2372  const TableIndices<rank_> &indices)
2373 {
2374  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2375  dim>(indices);
2376 }
2377 
2378 
2379 
2380 namespace internal
2381 {
2382  namespace SymmetricTensorImplementation
2383  {
2384  // a function to do the inverse of the unrolling from a set of
2385  // indices to a scalar index into the array in which we store
2386  // the elements of a symmetric tensor. in other words, it goes
2387  // from the scalar index into the array to a set of indices of
2388  // the tensor
2389  //
2390  // this function is for rank-2 tensors
2391  template <int dim>
2393  unrolled_to_component_indices(const unsigned int i,
2394  const std::integral_constant<int, 2> &)
2395  {
2396  Assert(
2398  ExcIndexRange(
2399  i,
2400  0,
2402  switch (dim)
2403  {
2404  case 1:
2405  {
2406  return {0, 0};
2407  }
2408 
2409  case 2:
2410  {
2411  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2412  TableIndices<2>(1, 1),
2413  TableIndices<2>(0, 1)};
2414  return table[i];
2415  }
2416 
2417  case 3:
2418  {
2419  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2420  TableIndices<2>(1, 1),
2421  TableIndices<2>(2, 2),
2422  TableIndices<2>(0, 1),
2423  TableIndices<2>(0, 2),
2424  TableIndices<2>(1, 2)};
2425  return table[i];
2426  }
2427 
2428  default:
2429  if (i < dim)
2430  return {i, i};
2431 
2432  for (unsigned int d = 0, c = dim; d < dim; ++d)
2433  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2434  if (c == i)
2435  return {d, e};
2436 
2437  // should never get here:
2438  Assert(false, ExcInternalError());
2439  return {0, 0};
2440  }
2441  }
2442 
2443  // a function to do the inverse of the unrolling from a set of
2444  // indices to a scalar index into the array in which we store
2445  // the elements of a symmetric tensor. in other words, it goes
2446  // from the scalar index into the array to a set of indices of
2447  // the tensor
2448  //
2449  // this function is for tensors of a rank not already handled
2450  // above
2451  template <int dim, int rank_>
2452  DEAL_II_CONSTEXPR inline
2453  typename std::enable_if<rank_ != 2, TableIndices<rank_>>::type
2454  unrolled_to_component_indices(const unsigned int i,
2455  const std::integral_constant<int, rank_> &)
2456  {
2457  (void)i;
2458  Assert(
2459  (i <
2461  ExcIndexRange(i,
2462  0,
2464  n_independent_components));
2465  Assert(false, ExcNotImplemented());
2466  return TableIndices<rank_>();
2467  }
2468 
2469  } // namespace SymmetricTensorImplementation
2470 } // namespace internal
2471 
2472 template <int rank_, int dim, typename Number>
2475  const unsigned int i)
2476 {
2477  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2478  dim>(i, std::integral_constant<int, rank_>());
2479 }
2480 
2481 
2482 
2483 template <int rank_, int dim, typename Number>
2484 template <class Archive>
2485 inline void
2486 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2487 {
2488  ar &data;
2489 }
2490 
2491 
2492 #endif // DOXYGEN
2493 
2494 /* ----------------- Non-member functions operating on tensors. ------------ */
2495 
2496 
2509 template <int rank_, int dim, typename Number, typename OtherNumber>
2514 {
2516  tmp = left;
2517  tmp += right;
2518  return tmp;
2519 }
2520 
2521 
2534 template <int rank_, int dim, typename Number, typename OtherNumber>
2539 {
2541  tmp = left;
2542  tmp -= right;
2543  return tmp;
2544 }
2545 
2546 
2554 template <int rank_, int dim, typename Number, typename OtherNumber>
2555 constexpr DEAL_II_ALWAYS_INLINE
2558  const Tensor<rank_, dim, OtherNumber> & right)
2559 {
2560  return Tensor<rank_, dim, Number>(left) + right;
2561 }
2562 
2563 
2571 template <int rank_, int dim, typename Number, typename OtherNumber>
2572 constexpr DEAL_II_ALWAYS_INLINE
2576 {
2577  return left + Tensor<rank_, dim, OtherNumber>(right);
2578 }
2579 
2580 
2588 template <int rank_, int dim, typename Number, typename OtherNumber>
2589 constexpr DEAL_II_ALWAYS_INLINE
2592  const Tensor<rank_, dim, OtherNumber> & right)
2593 {
2594  return Tensor<rank_, dim, Number>(left) - right;
2595 }
2596 
2597 
2605 template <int rank_, int dim, typename Number, typename OtherNumber>
2606 constexpr DEAL_II_ALWAYS_INLINE
2610 {
2611  return left - Tensor<rank_, dim, OtherNumber>(right);
2612 }
2613 
2614 
2615 
2629 template <int dim, typename Number>
2632 {
2633  switch (dim)
2634  {
2635  case 1:
2636  return t.data[0];
2637  case 2:
2638  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2639  case 3:
2640  {
2641  // in analogy to general tensors, but
2642  // there's something to be simplified for
2643  // the present case
2644  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2645  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2646  t.data[0] * t.data[5] * t.data[5] -
2647  t.data[1] * t.data[4] * t.data[4] -
2648  t.data[2] * t.data[3] * t.data[3]);
2649  }
2650  default:
2651  Assert(false, ExcNotImplemented());
2653  }
2654 }
2655 
2656 
2657 
2669 template <int dim, typename Number>
2670 constexpr DEAL_II_ALWAYS_INLINE Number
2672 {
2673  return determinant(t);
2674 }
2675 
2676 
2677 
2687 template <int dim, typename Number>
2690 {
2691  Number t = d.data[0];
2692  for (unsigned int i = 1; i < dim; ++i)
2693  t += d.data[i];
2694  return t;
2695 }
2696 
2697 
2709 template <int dim, typename Number>
2710 constexpr Number
2712 {
2713  return trace(t);
2714 }
2715 
2716 
2728 template <typename Number>
2729 constexpr DEAL_II_ALWAYS_INLINE Number
2731 {
2733 }
2734 
2735 
2736 
2755 template <typename Number>
2756 constexpr DEAL_II_ALWAYS_INLINE Number
2758 {
2759  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2760 }
2761 
2762 
2763 
2772 template <typename Number>
2773 constexpr DEAL_II_ALWAYS_INLINE Number
2775 {
2776  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2777  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2778 }
2779 
2780 
2781 
2789 template <typename Number>
2790 std::array<Number, 1>
2792 
2793 
2794 
2817 template <typename Number>
2818 std::array<Number, 2>
2820 
2821 
2822 
2845 template <typename Number>
2846 std::array<Number, 3>
2848 
2849 
2850 
2851 namespace internal
2852 {
2853  namespace SymmetricTensorImplementation
2854  {
2892  template <int dim, typename Number>
2893  void
2894  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2895  ::Tensor<2, dim, Number> & Q,
2896  std::array<Number, dim> & d,
2897  std::array<Number, dim - 1> & e);
2898 
2899 
2900 
2940  template <int dim, typename Number>
2941  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2942  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2943 
2944 
2945 
2985  template <int dim, typename Number>
2986  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2987  jacobi(::SymmetricTensor<2, dim, Number> A);
2988 
2989 
2990 
3004  template <typename Number>
3005  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3006  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3007 
3008 
3009 
3042  template <typename Number>
3043  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3044  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3045 
3050  template <int dim, typename Number>
3052  {
3053  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3054  bool
3055  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3056  {
3057  return lhs.first > rhs.first;
3058  }
3059  };
3060 
3061  } // namespace SymmetricTensorImplementation
3062 
3063 } // namespace internal
3064 
3065 
3066 
3067 // The line below is to ensure that doxygen puts the full description
3068 // of this global enumeration into the documentation
3069 // See https://stackoverflow.com/a/1717984
3099 {
3109  hybrid,
3119  ql_implicit_shifts,
3127  jacobi
3128 };
3129 
3130 
3131 
3160 template <int dim, typename Number>
3161 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3162  std::integral_constant<int, dim>::value>
3164  const SymmetricTensorEigenvectorMethod method =
3166 
3167 
3168 
3177 template <int rank_, int dim, typename Number>
3180 {
3181  return t;
3182 }
3183 
3184 
3185 
3196 template <int dim, typename Number>
3199 {
3201 
3202  // subtract scaled trace from the diagonal
3203  const Number tr = trace(t) / dim;
3204  for (unsigned int i = 0; i < dim; ++i)
3205  tmp.data[i] -= tr;
3206 
3207  return tmp;
3208 }
3209 
3210 
3211 
3218 template <int dim, typename Number>
3221 {
3222  // create a default constructed matrix filled with
3223  // zeros, then set the diagonal elements to one
3225  switch (dim)
3226  {
3227  case 1:
3229  break;
3230  case 2:
3231  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3232  break;
3233  case 3:
3234  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3236  break;
3237  default:
3238  for (unsigned int d = 0; d < dim; ++d)
3240  }
3241  return tmp;
3242 }
3243 
3244 
3245 
3253 template <int dim>
3256 {
3257  return unit_symmetric_tensor<dim, double>();
3258 }
3259 
3260 
3261 
3290 template <int dim, typename Number>
3293 {
3295 
3296  // fill the elements treating the diagonal
3297  for (unsigned int i = 0; i < dim; ++i)
3298  for (unsigned int j = 0; j < dim; ++j)
3299  tmp.data[i][j] =
3300  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3301 
3302  // then fill the ones that copy over the
3303  // non-diagonal elements. note that during
3304  // the double-contraction, we handle the
3305  // off-diagonal elements twice, so simply
3306  // copying requires a weight of 1/2
3307  for (unsigned int i = dim;
3308  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3309  n_rank2_components;
3310  ++i)
3311  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3312 
3313  return tmp;
3314 }
3315 
3316 
3317 
3325 template <int dim>
3328 {
3329  return deviator_tensor<dim, double>();
3330 }
3331 
3332 
3333 
3371 template <int dim, typename Number>
3374 {
3376 
3377  // fill the elements treating the diagonal
3378  for (unsigned int i = 0; i < dim; ++i)
3379  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3380 
3381  // then fill the ones that copy over the
3382  // non-diagonal elements. note that during
3383  // the double-contraction, we handle the
3384  // off-diagonal elements twice, so simply
3385  // copying requires a weight of 1/2
3386  for (unsigned int i = dim;
3387  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3388  n_rank2_components;
3389  ++i)
3390  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3391 
3392  return tmp;
3393 }
3394 
3395 
3396 
3404 template <int dim>
3407 {
3408  return identity_tensor<dim, double>();
3409 }
3410 
3411 
3412 
3422 template <int dim, typename Number>
3425 {
3427  value(t);
3428 }
3429 
3430 
3431 
3442 template <int dim, typename Number>
3445 {
3447  value(t);
3448 }
3449 
3450 
3451 
3473 template <int dim, typename Number>
3477 {
3479 
3480  // fill only the elements really needed
3481  for (unsigned int i = 0; i < dim; ++i)
3482  for (unsigned int j = i; j < dim; ++j)
3483  for (unsigned int k = 0; k < dim; ++k)
3484  for (unsigned int l = k; l < dim; ++l)
3485  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3486 
3487  return tmp;
3488 }
3489 
3490 
3491 
3499 template <int dim, typename Number>
3502 {
3504  for (unsigned int d = 0; d < dim; ++d)
3505  result[d][d] = t[d][d];
3506 
3507  const Number half = internal::NumberType<Number>::value(0.5);
3508  for (unsigned int d = 0; d < dim; ++d)
3509  for (unsigned int e = d + 1; e < dim; ++e)
3510  result[d][e] = (t[d][e] + t[e][d]) * half;
3511  return result;
3512 }
3513 
3514 
3515 
3523 template <int rank_, int dim, typename Number>
3526  operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3527 {
3529  tt *= factor;
3530  return tt;
3531 }
3532 
3533 
3534 
3542 template <int rank_, int dim, typename Number>
3544  operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3545 {
3546  // simply forward to the other operator
3547  return t * factor;
3548 }
3549 
3550 
3551 
3577 template <int rank_, int dim, typename Number, typename OtherNumber>
3579  rank_,
3580  dim,
3581  typename ProductType<Number,
3584  const OtherNumber & factor)
3585 {
3586  // form the product. we have to convert the two factors into the final
3587  // type via explicit casts because, for awkward reasons, the C++
3588  // standard committee saw it fit to not define an
3589  // operator*(float,std::complex<double>)
3590  // (as well as with switched arguments and double<->float).
3591  using product_type = typename ProductType<Number, OtherNumber>::type;
3594  return tt;
3595 }
3596 
3597 
3598 
3607 template <int rank_, int dim, typename Number, typename OtherNumber>
3609  rank_,
3610  dim,
3611  typename ProductType<OtherNumber,
3613 operator*(const Number & factor,
3615 {
3616  // simply forward to the other operator with switched arguments
3617  return (t * factor);
3618 }
3619 
3620 
3621 
3627 template <int rank_, int dim, typename Number, typename OtherNumber>
3629  rank_,
3630  dim,
3631  typename ProductType<Number,
3634  const OtherNumber & factor)
3635 {
3636  using product_type = typename ProductType<Number, OtherNumber>::type;
3639  return tt;
3640 }
3641 
3642 
3643 
3650 template <int rank_, int dim>
3652  operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3653 {
3655  tt *= factor;
3656  return tt;
3657 }
3658 
3659 
3660 
3667 template <int rank_, int dim>
3669  operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3670 {
3672  tt *= factor;
3673  return tt;
3674 }
3675 
3676 
3677 
3683 template <int rank_, int dim>
3685 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3686 {
3688  tt /= factor;
3689  return tt;
3690 }
3691 
3701 template <int dim, typename Number, typename OtherNumber>
3705 {
3706  return (t1 * t2);
3707 }
3708 
3709 
3723 template <int dim, typename Number, typename OtherNumber>
3727  const Tensor<2, dim, OtherNumber> & t2)
3728 {
3731  for (unsigned int i = 0; i < dim; ++i)
3732  for (unsigned int j = 0; j < dim; ++j)
3733  s += t1[i][j] * t2[i][j];
3734  return s;
3735 }
3736 
3737 
3751 template <int dim, typename Number, typename OtherNumber>
3755 {
3756  return scalar_product(t2, t1);
3757 }
3758 
3759 
3774 template <typename Number, typename OtherNumber>
3779 {
3780  tmp[0][0] = t[0][0][0][0] * s[0][0];
3781 }
3782 
3783 
3784 
3799 template <typename Number, typename OtherNumber>
3804 {
3805  tmp[0][0] = t[0][0][0][0] * s[0][0];
3806 }
3807 
3808 
3809 
3824 template <typename Number, typename OtherNumber>
3829 {
3830  const unsigned int dim = 2;
3831 
3832  for (unsigned int i = 0; i < dim; ++i)
3833  for (unsigned int j = i; j < dim; ++j)
3834  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3835  2 * t[i][j][0][1] * s[0][1];
3836 }
3837 
3838 
3839 
3854 template <typename Number, typename OtherNumber>
3859 {
3860  const unsigned int dim = 2;
3861 
3862  for (unsigned int i = 0; i < dim; ++i)
3863  for (unsigned int j = i; j < dim; ++j)
3864  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3865  2 * s[0][1] * t[0][1][i][j];
3866 }
3867 
3868 
3869 
3884 template <typename Number, typename OtherNumber>
3889 {
3890  const unsigned int dim = 3;
3891 
3892  for (unsigned int i = 0; i < dim; ++i)
3893  for (unsigned int j = i; j < dim; ++j)
3894  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3895  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3896  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3897 }
3898 
3899 
3900 
3915 template <typename Number, typename OtherNumber>
3920 {
3921  const unsigned int dim = 3;
3922 
3923  for (unsigned int i = 0; i < dim; ++i)
3924  for (unsigned int j = i; j < dim; ++j)
3925  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3926  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3927  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3928 }
3929 
3930 
3931 
3938 template <int dim, typename Number, typename OtherNumber>
3942  const Tensor<1, dim, OtherNumber> & src2)
3943 {
3945  for (unsigned int i = 0; i < dim; ++i)
3946  for (unsigned int j = 0; j < dim; ++j)
3947  dest[i] += src1[i][j] * src2[j];
3948  return dest;
3949 }
3950 
3951 
3958 template <int dim, typename Number, typename OtherNumber>
3962 {
3963  // this is easy for symmetric tensors:
3964  return src2 * src1;
3965 }
3966 
3967 
3968 
3988 template <int rank_1,
3989  int rank_2,
3990  int dim,
3991  typename Number,
3992  typename OtherNumber>
3993 constexpr DEAL_II_ALWAYS_INLINE
3994  typename Tensor<rank_1 + rank_2 - 2,
3995  dim,
3996  typename ProductType<Number, OtherNumber>::type>::tensor_type
3999 {
4000  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4001 }
4002 
4003 
4004 
4024 template <int rank_1,
4025  int rank_2,
4026  int dim,
4027  typename Number,
4028  typename OtherNumber>
4029 constexpr DEAL_II_ALWAYS_INLINE
4030  typename Tensor<rank_1 + rank_2 - 2,
4031  dim,
4032  typename ProductType<Number, OtherNumber>::type>::tensor_type
4034  const Tensor<rank_2, dim, OtherNumber> & src2)
4035 {
4036  return Tensor<rank_1, dim, Number>(src1) * src2;
4037 }
4038 
4039 
4040 
4050 template <int dim, typename Number>
4051 inline std::ostream &
4052 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4053 {
4054  // make our lives a bit simpler by outputting
4055  // the tensor through the operator for the
4056  // general Tensor class
4058 
4059  for (unsigned int i = 0; i < dim; ++i)
4060  for (unsigned int j = 0; j < dim; ++j)
4061  tt[i][j] = t[i][j];
4062 
4063  return out << tt;
4064 }
4065 
4066 
4067 
4077 template <int dim, typename Number>
4078 inline std::ostream &
4079 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4080 {
4081  // make our lives a bit simpler by outputting
4082  // the tensor through the operator for the
4083  // general Tensor class
4085 
4086  for (unsigned int i = 0; i < dim; ++i)
4087  for (unsigned int j = 0; j < dim; ++j)
4088  for (unsigned int k = 0; k < dim; ++k)
4089  for (unsigned int l = 0; l < dim; ++l)
4090  tt[i][j][k][l] = t[i][j][k][l];
4091 
4092  return out << tt;
4093 }
4094 
4095 
4097 
4098 #endif
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
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 SymmetricTensor operator-() const
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
static constexpr const T & value(const T &t)
Definition: numbers.h:693
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
constexpr SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr numbers::NumberTraits< Number >::real_type norm() const
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
typename AccessorTypes< rank, dim, constness, Number >::reference reference
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
constexpr bool operator==(const SymmetricTensor &) const
STL namespace.
static real_type abs(const number &x)
Definition: numbers.h:599
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
SymmetricTensorEigenvectorMethod
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
constexpr void sort()
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
constexpr void clear()
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
constexpr SymmetricTensor & operator/=(const OtherNumber &factor)
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr bool operator!=(const SymmetricTensor &) const
constexpr SymmetricTensor()=default
typename base_tensor_descriptor::base_tensor_type base_tensor_type
static const char T
T sum(const T &t, const MPI_Comm &mpi_communicator)
Number * begin_raw()
#define Assert(cond, exc)
Definition: exceptions.h:1466
std::pair< Number, Tensor< 1, dim, Number > > EigValsVecs
bool operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
constexpr TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
base_tensor_type data
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
typename ProductType< Number, OtherNumber >::type value_type
constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
void serialize(Archive &ar, const unsigned int version)
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A, ::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim - 1 > &e)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:94
Number * end_raw()
typename AccessorTypes< rank, dim, constness, Number >::reference reference
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 Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
::SymmetricTensor< rank1+rank2 - 4, dim, value_type > type
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
static const char A
static constexpr std::size_t memory_consumption()
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
constexpr Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: tensor.h:449
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
T min(const T &t, const MPI_Comm &mpi_communicator)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:931
static const char N
decltype(std::declval< T >() *std::declval< U >()) type
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
static ::ExceptionBase & ExcNotImplemented()
constexpr Number & operator()(const TableIndices< rank_ > &indices)
constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_CONSTEXPR
Definition: config.h:155
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)