Reference documentation for deal.II version Git 32ab9f15fd 2020-11-24 23:04:10 -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 
928  template <class Archive>
929  void
930  serialize(Archive &ar, const unsigned int version);
931 
932 private:
936  using base_tensor_descriptor =
938 
942  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
943 
948 
949  // Make all other symmetric tensors friends.
950  template <int, int, typename>
951  friend class SymmetricTensor;
952 
953  // Make a few more functions friends.
954  template <int dim2, typename Number2>
955  friend DEAL_II_CONSTEXPR Number2
957 
958  template <int dim2, typename Number2>
959  friend DEAL_II_CONSTEXPR Number2
961 
962  template <int dim2, typename Number2>
965 
966  template <int dim2, typename Number2>
969 
970  template <int dim2, typename Number2>
972  deviator_tensor();
973 
974  template <int dim2, typename Number2>
976  identity_tensor();
977 
978 
979  // Make a few helper classes friends as well.
981  Inverse<2, dim, Number>;
982 
984  Inverse<4, dim, Number>;
985 };
986 
987 
988 
989 // ------------------------- inline functions ------------------------
990 
991 #ifndef DOXYGEN
992 
993 // provide declarations for static members
994 template <int rank, int dim, typename Number>
995 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
996 
997 template <int rank_, int dim, typename Number>
998 constexpr unsigned int
999  SymmetricTensor<rank_, dim, Number>::n_independent_components;
1000 
1001 namespace internal
1002 {
1003  namespace SymmetricTensorAccessors
1004  {
1005  template <int rank_, int dim, bool constness, int P, typename Number>
1006  constexpr DEAL_II_ALWAYS_INLINE
1007  Accessor<rank_, dim, constness, P, Number>::Accessor(
1008  tensor_type & tensor,
1009  const TableIndices<rank_> &previous_indices)
1010  : tensor(tensor)
1011  , previous_indices(previous_indices)
1012  {}
1013 
1014 
1015 
1016  template <int rank_, int dim, bool constness, int P, typename Number>
1018  Accessor<rank_, dim, constness, P - 1, Number>
1019  Accessor<rank_, dim, constness, P, Number>::
1020  operator[](const unsigned int i)
1021  {
1022  return Accessor<rank_, dim, constness, P - 1, Number>(
1023  tensor, merge(previous_indices, i, rank_ - P));
1024  }
1025 
1026 
1027 
1028  template <int rank_, int dim, bool constness, int P, typename Number>
1029  constexpr DEAL_II_ALWAYS_INLINE
1030  Accessor<rank_, dim, constness, P - 1, Number>
1031  Accessor<rank_, dim, constness, P, Number>::
1032  operator[](const unsigned int i) const
1033  {
1034  return Accessor<rank_, dim, constness, P - 1, Number>(
1035  tensor, merge(previous_indices, i, rank_ - P));
1036  }
1037 
1038 
1039 
1040  template <int rank_, int dim, bool constness, typename Number>
1041  constexpr DEAL_II_ALWAYS_INLINE
1042  Accessor<rank_, dim, constness, 1, Number>::Accessor(
1043  tensor_type & tensor,
1044  const TableIndices<rank_> &previous_indices)
1045  : tensor(tensor)
1046  , previous_indices(previous_indices)
1047  {}
1048 
1049 
1050 
1051  template <int rank_, int dim, bool constness, typename Number>
1053  typename Accessor<rank_, dim, constness, 1, Number>::reference
1054  Accessor<rank_, dim, constness, 1, Number>::
1055  operator[](const unsigned int i)
1056  {
1057  return tensor(merge(previous_indices, i, rank_ - 1));
1058  }
1059 
1060 
1061  template <int rank_, int dim, bool constness, typename Number>
1062  constexpr DEAL_II_ALWAYS_INLINE
1063  typename Accessor<rank_, dim, constness, 1, Number>::reference
1064  Accessor<rank_, dim, constness, 1, Number>::
1065  operator[](const unsigned int i) const
1066  {
1067  return tensor(merge(previous_indices, i, rank_ - 1));
1068  }
1069  } // namespace SymmetricTensorAccessors
1070 } // namespace internal
1071 
1072 
1073 
1074 template <int rank_, int dim, typename Number>
1075 template <typename OtherNumber>
1076 inline DEAL_II_ALWAYS_INLINE
1078  const Tensor<2, dim, OtherNumber> &t)
1079 {
1080  static_assert(rank == 2, "This function is only implemented for rank==2");
1081  for (unsigned int d = 0; d < dim; ++d)
1082  for (unsigned int e = 0; e < d; ++e)
1083  Assert(t[d][e] == t[e][d],
1084  ExcMessage("The incoming Tensor must be exactly symmetric."));
1085 
1086  for (unsigned int d = 0; d < dim; ++d)
1087  data[d] = t[d][d];
1088 
1089  for (unsigned int d = 0, c = 0; d < dim; ++d)
1090  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1091  data[dim + c] = t[d][e];
1092 }
1093 
1094 
1095 
1096 template <int rank_, int dim, typename Number>
1097 template <typename OtherNumber>
1098 constexpr DEAL_II_ALWAYS_INLINE
1100  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1101  : data(initializer.data)
1102 {}
1103 
1104 
1105 
1106 template <int rank_, int dim, typename Number>
1109  const Number (&array)[n_independent_components])
1110  : data(
1111  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1112 {
1113  // ensure that the reinterpret_cast above actually works
1114  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1115  ExcInternalError());
1116 }
1117 
1118 
1119 
1120 template <int rank_, int dim, typename Number>
1121 template <typename OtherNumber>
1126 {
1127  data = t.data;
1128  return *this;
1129 }
1130 
1131 
1132 
1133 template <int rank_, int dim, typename Number>
1137 {
1139  ExcMessage("Only assignment with zero is allowed"));
1140  (void)d;
1141 
1143 
1144  return *this;
1145 }
1146 
1147 
1148 namespace internal
1149 {
1150  namespace SymmetricTensorImplementation
1151  {
1152  template <int dim, typename Number>
1153  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1154  ::Tensor<2, dim, Number>
1155  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1156  {
1158 
1159  // diagonal entries are stored first
1160  for (unsigned int d = 0; d < dim; ++d)
1161  t[d][d] = s.access_raw_entry(d);
1162 
1163  // off-diagonal entries come next, row by row
1164  for (unsigned int d = 0, c = 0; d < dim; ++d)
1165  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1166  {
1167  t[d][e] = s.access_raw_entry(dim + c);
1168  t[e][d] = s.access_raw_entry(dim + c);
1169  }
1170  return t;
1171  }
1172 
1173 
1174  template <int dim, typename Number>
1175  DEAL_II_CONSTEXPR ::Tensor<4, dim, Number>
1176  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1177  {
1178  // utilize the symmetry properties of SymmetricTensor<4,dim>
1179  // discussed in the class documentation to avoid accessing all
1180  // independent elements of the input tensor more than once
1182 
1183  for (unsigned int i = 0; i < dim; ++i)
1184  for (unsigned int j = i; j < dim; ++j)
1185  for (unsigned int k = 0; k < dim; ++k)
1186  for (unsigned int l = k; l < dim; ++l)
1187  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1188  t[TableIndices<4>(j, i, k, l)] =
1189  t[TableIndices<4>(j, i, l, k)] =
1190  st[TableIndices<4>(i, j, k, l)];
1191 
1192  return t;
1193  }
1194 
1195 
1196  template <typename Number>
1197  struct Inverse<2, 1, Number>
1198  {
1199  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1200  ::SymmetricTensor<2, 1, Number>
1201  value(const ::SymmetricTensor<2, 1, Number> &t)
1202  {
1204 
1205  tmp[0][0] = 1.0 / t[0][0];
1206 
1207  return tmp;
1208  }
1209  };
1210 
1211 
1212  template <typename Number>
1213  struct Inverse<2, 2, Number>
1214  {
1215  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1216  ::SymmetricTensor<2, 2, Number>
1217  value(const ::SymmetricTensor<2, 2, Number> &t)
1218  {
1220 
1221  // Sympy result: ([
1222  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1223  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1224  const TableIndices<2> idx_00(0, 0);
1225  const TableIndices<2> idx_01(0, 1);
1226  const TableIndices<2> idx_11(1, 1);
1227  const Number inv_det_t =
1228  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1229  tmp[idx_00] = t[idx_11];
1230  tmp[idx_01] = -t[idx_01];
1231  tmp[idx_11] = t[idx_00];
1232  tmp *= inv_det_t;
1233 
1234  return tmp;
1235  }
1236  };
1237 
1238 
1239  template <typename Number>
1240  struct Inverse<2, 3, Number>
1241  {
1242  DEAL_II_CONSTEXPR static ::SymmetricTensor<2, 3, Number>
1243  value(const ::SymmetricTensor<2, 3, Number> &t)
1244  {
1246 
1247  // Sympy result: ([
1248  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1249  // 2*t01*t02*t12 - t02**2*t11),
1250  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1251  // 2*t01*t02*t12 - t02**2*t11),
1252  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1253  // 2*t01*t02*t12 - t02**2*t11)],
1254  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1255  // 2*t01*t02*t12 - t02**2*t11),
1256  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1257  // 2*t01*t02*t12 - t02**2*t11),
1258  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1259  // 2*t01*t02*t12 + t02**2*t11)],
1260  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1261  // 2*t01*t02*t12 - t02**2*t11),
1262  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1263  // 2*t01*t02*t12 + t02**2*t11),
1264  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1265  // 2*t01*t02*t12 + t02**2*t11)] ])
1266  //
1267  // =
1268  //
1269  // [ (t11*t22 - t12**2)/det_t,
1270  // (-t01*t22 + t02*t12)/det_t,
1271  // (t01*t12 - t02*t11)/det_t],
1272  // [ (-t01*t22 + t02*t12)/det_t,
1273  // (t00*t22 - t02**2)/det_t,
1274  // (-t00*t12 + t01*t02)/det_t],
1275  // [ (t01*t12 - t02*t11)/det_t,
1276  // (-t00*t12 + t01*t02)/det_t,
1277  // (t00*t11 - t01**2)/det_t] ])
1278  //
1279  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1280  // 2*t01*t02*t12 - t02**2*t11)
1281  const TableIndices<2> idx_00(0, 0);
1282  const TableIndices<2> idx_01(0, 1);
1283  const TableIndices<2> idx_02(0, 2);
1284  const TableIndices<2> idx_11(1, 1);
1285  const TableIndices<2> idx_12(1, 2);
1286  const TableIndices<2> idx_22(2, 2);
1287  const Number inv_det_t =
1288  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1289  t[idx_00] * t[idx_12] * t[idx_12] -
1290  t[idx_01] * t[idx_01] * t[idx_22] +
1291  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1292  t[idx_02] * t[idx_02] * t[idx_11]);
1293  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1294  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1295  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1296  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1297  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1298  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1299  tmp *= inv_det_t;
1300 
1301  return tmp;
1302  }
1303  };
1304 
1305 
1306  template <typename Number>
1307  struct Inverse<4, 1, Number>
1308  {
1309  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 1, Number>
1310  value(const ::SymmetricTensor<4, 1, Number> &t)
1311  {
1313  tmp.data[0][0] = 1.0 / t.data[0][0];
1314  return tmp;
1315  }
1316  };
1317 
1318 
1319  template <typename Number>
1320  struct Inverse<4, 2, Number>
1321  {
1322  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 2, Number>
1323  value(const ::SymmetricTensor<4, 2, Number> &t)
1324  {
1326 
1327  // Inverting this tensor is a little more complicated than necessary,
1328  // since we store the data of 't' as a 3x3 matrix t.data, but the
1329  // product between a rank-4 and a rank-2 tensor is really not the
1330  // product between this matrix and the 3-vector of a rhs, but rather
1331  //
1332  // B.vec = t.data * mult * A.vec
1333  //
1334  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1335  // capture the fact that we need to add up both the c_ij12*a_12 and the
1336  // c_ij21*a_21 terms.
1337  //
1338  // In addition, in this scheme, the identity tensor has the matrix
1339  // representation mult^-1.
1340  //
1341  // The inverse of 't' therefore has the matrix representation
1342  //
1343  // inv.data = mult^-1 * t.data^-1 * mult^-1
1344  //
1345  // in order to compute it, let's first compute the inverse of t.data and
1346  // put it into tmp.data; at the end of the function we then scale the
1347  // last row and column of the inverse by 1/2, corresponding to the left
1348  // and right multiplication with mult^-1.
1349  const Number t4 = t.data[0][0] * t.data[1][1],
1350  t6 = t.data[0][0] * t.data[1][2],
1351  t8 = t.data[0][1] * t.data[1][0],
1352  t00 = t.data[0][2] * t.data[1][0],
1353  t01 = t.data[0][1] * t.data[2][0],
1354  t04 = t.data[0][2] * t.data[2][0],
1355  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1356  t8 * t.data[2][2] + t00 * t.data[2][1] +
1357  t01 * t.data[1][2] - t04 * t.data[1][1]);
1358  tmp.data[0][0] =
1359  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1360  tmp.data[0][1] =
1361  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1362  tmp.data[0][2] =
1363  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1364  tmp.data[1][0] =
1365  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1366  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1367  tmp.data[1][2] = -(t6 - t00) * t07;
1368  tmp.data[2][0] =
1369  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1370  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1371  tmp.data[2][2] = (t4 - t8) * t07;
1372 
1373  // scale last row and column as mentioned
1374  // above
1375  tmp.data[2][0] /= 2;
1376  tmp.data[2][1] /= 2;
1377  tmp.data[0][2] /= 2;
1378  tmp.data[1][2] /= 2;
1379  tmp.data[2][2] /= 4;
1380 
1381  return tmp;
1382  }
1383  };
1384 
1385 
1386  template <typename Number>
1387  struct Inverse<4, 3, Number>
1388  {
1389  static ::SymmetricTensor<4, 3, Number>
1390  value(const ::SymmetricTensor<4, 3, Number> &t)
1391  {
1393 
1394  // This function follows the exact same scheme as the 2d case, except
1395  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1396  // Instead, we use the Gauss-Jordan algorithm implemented for
1397  // FullMatrix. For historical reasons the following code is copied from
1398  // there, with the tangential benefit that we do not need to copy the
1399  // tensor entries to and from the FullMatrix.
1400  const unsigned int N = 6;
1401 
1402  // First get an estimate of the size of the elements of this matrix,
1403  // for later checks whether the pivot element is large enough, or
1404  // whether we have to fear that the matrix is not regular.
1405  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1406  for (unsigned int i = 0; i < N; ++i)
1407  diagonal_sum += std::fabs(tmp.data[i][i]);
1408  const Number typical_diagonal_element =
1409  diagonal_sum / static_cast<double>(N);
1410  (void)typical_diagonal_element;
1411 
1412  unsigned int p[N];
1413  for (unsigned int i = 0; i < N; ++i)
1414  p[i] = i;
1415 
1416  for (unsigned int j = 0; j < N; ++j)
1417  {
1418  // Pivot search: search that part of the line on and right of the
1419  // diagonal for the largest element.
1420  Number max = std::fabs(tmp.data[j][j]);
1421  unsigned int r = j;
1422  for (unsigned int i = j + 1; i < N; ++i)
1423  if (std::fabs(tmp.data[i][j]) > max)
1424  {
1425  max = std::fabs(tmp.data[i][j]);
1426  r = i;
1427  }
1428 
1429  // Check whether the pivot is too small
1430  Assert(max > 1.e-16 * typical_diagonal_element,
1431  ExcMessage("This tensor seems to be noninvertible"));
1432 
1433  // Row interchange
1434  if (r > j)
1435  {
1436  for (unsigned int k = 0; k < N; ++k)
1437  std::swap(tmp.data[j][k], tmp.data[r][k]);
1438 
1439  std::swap(p[j], p[r]);
1440  }
1441 
1442  // Transformation
1443  const Number hr = 1. / tmp.data[j][j];
1444  tmp.data[j][j] = hr;
1445  for (unsigned int k = 0; k < N; ++k)
1446  {
1447  if (k == j)
1448  continue;
1449  for (unsigned int i = 0; i < N; ++i)
1450  {
1451  if (i == j)
1452  continue;
1453  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1454  }
1455  }
1456  for (unsigned int i = 0; i < N; ++i)
1457  {
1458  tmp.data[i][j] *= hr;
1459  tmp.data[j][i] *= -hr;
1460  }
1461  tmp.data[j][j] = hr;
1462  }
1463 
1464  // Column interchange
1465  Number hv[N];
1466  for (unsigned int i = 0; i < N; ++i)
1467  {
1468  for (unsigned int k = 0; k < N; ++k)
1469  hv[p[k]] = tmp.data[i][k];
1470  for (unsigned int k = 0; k < N; ++k)
1471  tmp.data[i][k] = hv[k];
1472  }
1473 
1474  // Scale rows and columns. The mult matrix
1475  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1476  for (unsigned int i = 3; i < 6; ++i)
1477  for (unsigned int j = 0; j < 3; ++j)
1478  tmp.data[i][j] /= 2;
1479 
1480  for (unsigned int i = 0; i < 3; ++i)
1481  for (unsigned int j = 3; j < 6; ++j)
1482  tmp.data[i][j] /= 2;
1483 
1484  for (unsigned int i = 3; i < 6; ++i)
1485  for (unsigned int j = 3; j < 6; ++j)
1486  tmp.data[i][j] /= 4;
1487 
1488  return tmp;
1489  }
1490  };
1491 
1492  } // namespace SymmetricTensorImplementation
1493 } // namespace internal
1494 
1495 
1496 
1497 template <int rank_, int dim, typename Number>
1499  operator Tensor<rank_, dim, Number>() const
1500 {
1501  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1502 }
1503 
1504 
1505 
1506 template <int rank_, int dim, typename Number>
1507 constexpr bool
1510 {
1511  return data == t.data;
1512 }
1513 
1514 
1515 
1516 template <int rank_, int dim, typename Number>
1517 constexpr bool
1520 {
1521  return data != t.data;
1522 }
1523 
1524 
1525 
1526 template <int rank_, int dim, typename Number>
1527 template <typename OtherNumber>
1532 {
1533  data += t.data;
1534  return *this;
1535 }
1536 
1537 
1538 
1539 template <int rank_, int dim, typename Number>
1540 template <typename OtherNumber>
1545 {
1546  data -= t.data;
1547  return *this;
1548 }
1549 
1550 
1551 
1552 template <int rank_, int dim, typename Number>
1553 template <typename OtherNumber>
1557 {
1558  data *= d;
1559  return *this;
1560 }
1561 
1562 
1563 
1564 template <int rank_, int dim, typename Number>
1565 template <typename OtherNumber>
1569 {
1570  data /= d;
1571  return *this;
1572 }
1573 
1574 
1575 
1576 template <int rank_, int dim, typename Number>
1580 {
1581  SymmetricTensor tmp = *this;
1582  tmp.data = -tmp.data;
1583  return tmp;
1584 }
1585 
1586 
1587 
1588 template <int rank_, int dim, typename Number>
1591 {
1592  data.clear();
1593 }
1594 
1595 
1596 
1597 template <int rank_, int dim, typename Number>
1598 constexpr std::size_t
1600 {
1601  // all memory consists of statically allocated memory of the current
1602  // object, no pointers
1603  return sizeof(SymmetricTensor<rank_, dim, Number>);
1604 }
1605 
1606 
1607 
1608 namespace internal
1609 {
1610  template <int dim, typename Number, typename OtherNumber = Number>
1614  perform_double_contraction(
1615  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1616  base_tensor_type &data,
1617  const typename SymmetricTensorAccessors::
1618  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1619  {
1620  using result_type = typename SymmetricTensorAccessors::
1622 
1623  switch (dim)
1624  {
1625  case 1:
1626  return data[0] * sdata[0];
1627  default:
1628  // Start with the non-diagonal part to avoid some multiplications by
1629  // 2.
1630 
1631  result_type sum = data[dim] * sdata[dim];
1632  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1633  sum += data[d] * sdata[d];
1634  sum += sum; // sum = sum * 2.;
1635 
1636  // Now add the contributions from the diagonal
1637  for (unsigned int d = 0; d < dim; ++d)
1638  sum += data[d] * sdata[d];
1639  return sum;
1640  }
1641  }
1642 
1643 
1644 
1645  template <int dim, typename Number, typename OtherNumber = Number>
1649  perform_double_contraction(
1650  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1651  base_tensor_type &data,
1652  const typename SymmetricTensorAccessors::
1653  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1654  {
1655  using result_type = typename SymmetricTensorAccessors::
1657  using value_type = typename SymmetricTensorAccessors::
1659 
1660  const unsigned int data_dim = SymmetricTensorAccessors::
1661  StorageType<2, dim, value_type>::n_independent_components;
1662  value_type tmp[data_dim]{};
1663  for (unsigned int i = 0; i < data_dim; ++i)
1664  tmp[i] =
1665  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1666  return result_type(tmp);
1667  }
1668 
1669 
1670 
1671  template <int dim, typename Number, typename OtherNumber = Number>
1673  typename SymmetricTensorAccessors::StorageType<
1674  2,
1675  dim,
1678  base_tensor_type
1679  perform_double_contraction(
1680  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1681  base_tensor_type &data,
1682  const typename SymmetricTensorAccessors::
1683  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1684  {
1685  using value_type = typename SymmetricTensorAccessors::
1687  using base_tensor_type = typename SymmetricTensorAccessors::
1688  StorageType<2, dim, value_type>::base_tensor_type;
1689 
1690  base_tensor_type tmp;
1691  for (unsigned int i = 0; i < tmp.dimension; ++i)
1692  {
1693  // Start with the non-diagonal part
1694  value_type sum = data[dim] * sdata[dim][i];
1695  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1696  sum += data[d] * sdata[d][i];
1697  sum += sum; // sum = sum * 2.;
1698 
1699  // Now add the contributions from the diagonal
1700  for (unsigned int d = 0; d < dim; ++d)
1701  sum += data[d] * sdata[d][i];
1702  tmp[i] = sum;
1703  }
1704  return tmp;
1705  }
1706 
1707 
1708 
1709  template <int dim, typename Number, typename OtherNumber = Number>
1711  typename SymmetricTensorAccessors::StorageType<
1712  4,
1713  dim,
1716  base_tensor_type
1717  perform_double_contraction(
1718  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1719  base_tensor_type &data,
1720  const typename SymmetricTensorAccessors::
1721  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1722  {
1723  using value_type = typename SymmetricTensorAccessors::
1725  using base_tensor_type = typename SymmetricTensorAccessors::
1726  StorageType<4, dim, value_type>::base_tensor_type;
1727 
1728  const unsigned int data_dim = SymmetricTensorAccessors::
1729  StorageType<2, dim, value_type>::n_independent_components;
1730  base_tensor_type tmp;
1731  for (unsigned int i = 0; i < data_dim; ++i)
1732  for (unsigned int j = 0; j < data_dim; ++j)
1733  {
1734  // Start with the non-diagonal part
1735  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1736  tmp[i][j] += data[i][d] * sdata[d][j];
1737  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1738 
1739  // Now add the contributions from the diagonal
1740  for (unsigned int d = 0; d < dim; ++d)
1741  tmp[i][j] += data[i][d] * sdata[d][j];
1742  }
1743  return tmp;
1744  }
1745 
1746 } // end of namespace internal
1747 
1748 
1749 
1750 template <int rank_, int dim, typename Number>
1751 template <typename OtherNumber>
1757 {
1758  // need to have two different function calls
1759  // because a scalar and rank-2 tensor are not
1760  // the same data type (see internal function
1761  // above)
1762  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1763  s.data);
1764 }
1765 
1766 
1767 
1768 template <int rank_, int dim, typename Number>
1769 template <typename OtherNumber>
1774 {
1777  tmp.data =
1778  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1779  s.data);
1780  return tmp;
1781 }
1782 
1783 
1784 
1785 // internal namespace to switch between the
1786 // access of different tensors. There used to
1787 // be explicit instantiations before for
1788 // different ranks and dimensions, but since
1789 // we now allow for templates on the data
1790 // type, and since we cannot partially
1791 // specialize the implementation, this got
1792 // into a separate namespace
1793 namespace internal
1794 {
1795  // The variables within this struct will be referenced in the next functions.
1796  // It is a workaround that allows returning a reference to a static variable
1797  // while allowing constexpr evaluation of the function.
1798  // It has to be defined outside the function because constexpr functions
1799  // cannot define static variables.
1800  // A similar struct has also been defined in tensor.h
1801  template <typename Type>
1802  struct Uninitialized
1803  {
1804  static Type value;
1805  };
1806 
1807  template <typename Type>
1809 
1810  template <int dim, typename Number>
1812  symmetric_tensor_access(const TableIndices<2> &indices,
1813  typename SymmetricTensorAccessors::
1814  StorageType<2, dim, Number>::base_tensor_type &data)
1815  {
1816  // 1d is very simple and done first
1817  if (dim == 1)
1818  return data[0];
1819 
1820  // first treat the main diagonal elements, which are stored consecutively
1821  // at the beginning
1822  if (indices[0] == indices[1])
1823  return data[indices[0]];
1824 
1825  // the rest is messier and requires a few switches.
1826  switch (dim)
1827  {
1828  case 2:
1829  // at least for the 2x2 case it is reasonably simple
1830  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1831  ((indices[0] == 0) && (indices[1] == 1)),
1832  ExcInternalError());
1833  return data[2];
1834 
1835  default:
1836  // to do the rest, sort our indices before comparing
1837  {
1838  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1839  std::max(indices[0], indices[1]));
1840  for (unsigned int d = 0, c = 0; d < dim; ++d)
1841  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1842  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1843  return data[dim + c];
1844  Assert(false, ExcInternalError());
1845  }
1846  }
1847 
1848  // The code should never reach there.
1849  // Returns a dummy reference to a dummy variable just to make the
1850  // compiler happy.
1852  }
1853 
1854 
1855 
1856  template <int dim, typename Number>
1857  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1858  symmetric_tensor_access(const TableIndices<2> &indices,
1859  const typename SymmetricTensorAccessors::
1860  StorageType<2, dim, Number>::base_tensor_type &data)
1861  {
1862  // 1d is very simple and done first
1863  if (dim == 1)
1864  return data[0];
1865 
1866  // first treat the main diagonal elements, which are stored consecutively
1867  // at the beginning
1868  if (indices[0] == indices[1])
1869  return data[indices[0]];
1870 
1871  // the rest is messier and requires a few switches.
1872  switch (dim)
1873  {
1874  case 2:
1875  // at least for the 2x2 case it is reasonably simple
1876  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1877  ((indices[0] == 0) && (indices[1] == 1)),
1878  ExcInternalError());
1879  return data[2];
1880 
1881  default:
1882  // to do the rest, sort our indices before comparing
1883  {
1884  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1885  std::max(indices[0], indices[1]));
1886  for (unsigned int d = 0, c = 0; d < dim; ++d)
1887  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1888  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1889  return data[dim + c];
1890  Assert(false, ExcInternalError());
1891  }
1892  }
1893 
1894  // The code should never reach there.
1895  // Returns a dummy reference to a dummy variable just to make the
1896  // compiler happy.
1898  }
1899 
1900 
1901 
1902  template <int dim, typename Number>
1903  DEAL_II_CONSTEXPR inline Number &
1904  symmetric_tensor_access(const TableIndices<4> &indices,
1905  typename SymmetricTensorAccessors::
1906  StorageType<4, dim, Number>::base_tensor_type &data)
1907  {
1908  switch (dim)
1909  {
1910  case 1:
1911  return data[0][0];
1912 
1913  case 2:
1914  // each entry of the tensor can be thought of as an entry in a
1915  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1916  // rank-2 tensors. this is the format in which we store rank-4
1917  // tensors. determine which position the present entry is
1918  // stored in
1919  {
1920  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1921  return data[base_index[indices[0]][indices[1]]]
1922  [base_index[indices[2]][indices[3]]];
1923  }
1924  case 3:
1925  // each entry of the tensor can be thought of as an entry in a
1926  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1927  // rank-2 tensors. this is the format in which we store rank-4
1928  // tensors. determine which position the present entry is
1929  // stored in
1930  {
1931  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1932  {3, 1, 5},
1933  {4, 5, 2}};
1934  return data[base_index[indices[0]][indices[1]]]
1935  [base_index[indices[2]][indices[3]]];
1936  }
1937 
1938  default:
1939  Assert(false, ExcNotImplemented());
1940  }
1941 
1942  // The code should never reach there.
1943  // Returns a dummy reference to a dummy variable just to make the
1944  // compiler happy.
1946  }
1947 
1948 
1949  template <int dim, typename Number>
1950  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1951  symmetric_tensor_access(const TableIndices<4> &indices,
1952  const typename SymmetricTensorAccessors::
1953  StorageType<4, dim, Number>::base_tensor_type &data)
1954  {
1955  switch (dim)
1956  {
1957  case 1:
1958  return data[0][0];
1959 
1960  case 2:
1961  // each entry of the tensor can be thought of as an entry in a
1962  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1963  // rank-2 tensors. this is the format in which we store rank-4
1964  // tensors. determine which position the present entry is
1965  // stored in
1966  {
1967  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1968  return data[base_index[indices[0]][indices[1]]]
1969  [base_index[indices[2]][indices[3]]];
1970  }
1971  case 3:
1972  // each entry of the tensor can be thought of as an entry in a
1973  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1974  // rank-2 tensors. this is the format in which we store rank-4
1975  // tensors. determine which position the present entry is
1976  // stored in
1977  {
1978  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1979  {3, 1, 5},
1980  {4, 5, 2}};
1981  return data[base_index[indices[0]][indices[1]]]
1982  [base_index[indices[2]][indices[3]]];
1983  }
1984 
1985  default:
1986  Assert(false, ExcNotImplemented());
1987  }
1988 
1989  // The code should never reach there.
1990  // Returns a dummy reference to a dummy variable just to make the
1991  // compiler happy.
1993  }
1994 
1995 } // end of namespace internal
1996 
1997 
1998 
1999 template <int rank_, int dim, typename Number>
2002  operator()(const TableIndices<rank_> &indices)
2003 {
2004  for (unsigned int r = 0; r < rank; ++r)
2005  AssertIndexRange(indices[r], dimension);
2006  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2007 }
2008 
2009 
2010 
2011 template <int rank_, int dim, typename Number>
2012 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
2014  operator()(const TableIndices<rank_> &indices) const
2015 {
2016  for (unsigned int r = 0; r < rank; ++r)
2017  AssertIndexRange(indices[r], dimension);
2018  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2019 }
2020 
2021 
2022 
2023 namespace internal
2024 {
2025  namespace SymmetricTensorImplementation
2026  {
2027  template <int rank_>
2028  constexpr TableIndices<rank_>
2029  get_partially_filled_indices(const unsigned int row,
2030  const std::integral_constant<int, 2> &)
2031  {
2033  }
2034 
2035 
2036  template <int rank_>
2037  constexpr TableIndices<rank_>
2038  get_partially_filled_indices(const unsigned int row,
2039  const std::integral_constant<int, 4> &)
2040  {
2041  return TableIndices<rank_>(row,
2045  }
2046  } // namespace SymmetricTensorImplementation
2047 } // namespace internal
2048 
2049 
2050 template <int rank_, int dim, typename Number>
2051 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2052  Accessor<rank_, dim, true, rank_ - 1, Number>
2054  operator[](const unsigned int row) const
2055 {
2056  return internal::SymmetricTensorAccessors::
2057  Accessor<rank_, dim, true, rank_ - 1, Number>(
2058  *this,
2059  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2060  rank_>(row, std::integral_constant<int, rank_>()));
2061 }
2062 
2063 
2064 
2065 template <int rank_, int dim, typename Number>
2066 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE internal::
2067  SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2068  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2069 {
2070  return internal::SymmetricTensorAccessors::
2071  Accessor<rank_, dim, false, rank_ - 1, Number>(
2072  *this,
2073  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2074  rank_>(row, std::integral_constant<int, rank_>()));
2075 }
2076 
2077 
2078 
2079 template <int rank_, int dim, typename Number>
2080 constexpr DEAL_II_ALWAYS_INLINE const Number &
2082  operator[](const TableIndices<rank_> &indices) const
2083 {
2084  return operator()(indices);
2085 }
2086 
2087 
2088 
2089 template <int rank_, int dim, typename Number>
2092  operator[](const TableIndices<rank_> &indices)
2093 {
2094  return operator()(indices);
2095 }
2096 
2097 
2098 
2099 template <int rank_, int dim, typename Number>
2100 inline Number *
2102 {
2103  return std::addressof(this->access_raw_entry(0));
2104 }
2105 
2106 
2107 
2108 template <int rank_, int dim, typename Number>
2109 inline const Number *
2111 {
2112  return std::addressof(this->access_raw_entry(0));
2113 }
2114 
2115 
2116 
2117 template <int rank_, int dim, typename Number>
2118 inline Number *
2120 {
2121  return begin_raw() + n_independent_components;
2122 }
2123 
2124 
2125 
2126 template <int rank_, int dim, typename Number>
2127 inline const Number *
2129 {
2130  return begin_raw() + n_independent_components;
2131 }
2132 
2133 
2134 
2135 namespace internal
2136 {
2137  namespace SymmetricTensorImplementation
2138  {
2139  template <int dim, typename Number>
2140  constexpr unsigned int
2141  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2142  const unsigned int index)
2143  {
2144  return index;
2145  }
2146 
2147 
2148  template <int dim, typename Number>
2149  constexpr ::TableIndices<2>
2150  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2151  const unsigned int index)
2152  {
2155  }
2156 
2157  } // namespace SymmetricTensorImplementation
2158 } // namespace internal
2159 
2160 
2161 
2162 template <int rank_, int dim, typename Number>
2163 DEAL_II_CONSTEXPR inline const Number &
2165  const unsigned int index) const
2166 {
2167  AssertIndexRange(index, n_independent_components);
2168  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2169  index)];
2170 }
2171 
2172 
2173 
2174 template <int rank_, int dim, typename Number>
2175 DEAL_II_CONSTEXPR inline Number &
2177 {
2178  AssertIndexRange(index, n_independent_components);
2179  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2180  index)];
2181 }
2182 
2183 
2184 
2185 namespace internal
2186 {
2187  template <int dim, typename Number>
2189  compute_norm(const typename SymmetricTensorAccessors::
2190  StorageType<2, dim, Number>::base_tensor_type &data)
2191  {
2192  switch (dim)
2193  {
2194  case 1:
2195  return numbers::NumberTraits<Number>::abs(data[0]);
2196 
2197  case 2:
2198  return std::sqrt(
2202 
2203  case 3:
2204  return std::sqrt(
2211 
2212  default:
2213  {
2214  typename numbers::NumberTraits<Number>::real_type return_value =
2216 
2217  for (unsigned int d = 0; d < dim; ++d)
2218  return_value +=
2220  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2221  return_value +=
2223 
2224  return std::sqrt(return_value);
2225  }
2226  }
2227  }
2228 
2229 
2230 
2231  template <int dim, typename Number>
2233  compute_norm(const typename SymmetricTensorAccessors::
2234  StorageType<4, dim, Number>::base_tensor_type &data)
2235  {
2236  switch (dim)
2237  {
2238  case 1:
2239  return numbers::NumberTraits<Number>::abs(data[0][0]);
2240 
2241  default:
2242  {
2243  typename numbers::NumberTraits<Number>::real_type return_value =
2245 
2246  const unsigned int n_independent_components = data.dimension;
2247 
2248  for (unsigned int i = 0; i < dim; ++i)
2249  for (unsigned int j = 0; j < dim; ++j)
2250  return_value +=
2252  for (unsigned int i = 0; i < dim; ++i)
2253  for (unsigned int j = dim; j < n_independent_components; ++j)
2254  return_value +=
2256  for (unsigned int i = dim; i < n_independent_components; ++i)
2257  for (unsigned int j = 0; j < dim; ++j)
2258  return_value +=
2260  for (unsigned int i = dim; i < n_independent_components; ++i)
2261  for (unsigned int j = dim; j < n_independent_components; ++j)
2262  return_value +=
2264 
2265  return std::sqrt(return_value);
2266  }
2267  }
2268  }
2269 
2270 } // end of namespace internal
2271 
2272 
2273 
2274 template <int rank_, int dim, typename Number>
2277 {
2278  return internal::compute_norm<dim, Number>(data);
2279 }
2280 
2281 
2282 
2283 namespace internal
2284 {
2285  namespace SymmetricTensorImplementation
2286  {
2287  // a function to do the unrolling from a set of indices to a
2288  // scalar index into the array in which we store the elements of
2289  // a symmetric tensor
2290  //
2291  // this function is for rank-2 tensors
2292  template <int dim>
2293  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE unsigned int
2294  component_to_unrolled_index(const TableIndices<2> &indices)
2295  {
2296  AssertIndexRange(indices[0], dim);
2297  AssertIndexRange(indices[1], dim);
2298 
2299  switch (dim)
2300  {
2301  case 1:
2302  {
2303  return 0;
2304  }
2305 
2306  case 2:
2307  {
2308  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2309  return table[indices[0]][indices[1]];
2310  }
2311 
2312  case 3:
2313  {
2314  constexpr unsigned int table[3][3] = {{0, 3, 4},
2315  {3, 1, 5},
2316  {4, 5, 2}};
2317  return table[indices[0]][indices[1]];
2318  }
2319 
2320  case 4:
2321  {
2322  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2323  {4, 1, 7, 8},
2324  {5, 7, 2, 9},
2325  {6, 8, 9, 3}};
2326  return table[indices[0]][indices[1]];
2327  }
2328 
2329  default:
2330  // for the remainder, manually figure out the numbering
2331  {
2332  if (indices[0] == indices[1])
2333  return indices[0];
2334 
2335  TableIndices<2> sorted_indices(indices);
2336  sorted_indices.sort();
2337 
2338  for (unsigned int d = 0, c = 0; d < dim; ++d)
2339  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2340  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2341  return dim + c;
2342 
2343  // should never get here:
2344  Assert(false, ExcInternalError());
2345  return 0;
2346  }
2347  }
2348  }
2349 
2350  // a function to do the unrolling from a set of indices to a
2351  // scalar index into the array in which we store the elements of
2352  // a symmetric tensor
2353  //
2354  // this function is for tensors of ranks not already handled
2355  // above
2356  template <int dim, int rank_>
2357  DEAL_II_CONSTEXPR inline unsigned int
2358  component_to_unrolled_index(const TableIndices<rank_> &indices)
2359  {
2360  (void)indices;
2361  Assert(false, ExcNotImplemented());
2363  }
2364  } // namespace SymmetricTensorImplementation
2365 } // namespace internal
2366 
2367 
2368 template <int rank_, int dim, typename Number>
2369 constexpr unsigned int
2371  const TableIndices<rank_> &indices)
2372 {
2373  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2374  dim>(indices);
2375 }
2376 
2377 
2378 
2379 namespace internal
2380 {
2381  namespace SymmetricTensorImplementation
2382  {
2383  // a function to do the inverse of the unrolling from a set of
2384  // indices to a scalar index into the array in which we store
2385  // the elements of a symmetric tensor. in other words, it goes
2386  // from the scalar index into the array to a set of indices of
2387  // the tensor
2388  //
2389  // this function is for rank-2 tensors
2390  template <int dim>
2392  unrolled_to_component_indices(const unsigned int i,
2393  const std::integral_constant<int, 2> &)
2394  {
2395  Assert(
2397  ExcIndexRange(
2398  i,
2399  0,
2401  switch (dim)
2402  {
2403  case 1:
2404  {
2405  return {0, 0};
2406  }
2407 
2408  case 2:
2409  {
2410  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2411  TableIndices<2>(1, 1),
2412  TableIndices<2>(0, 1)};
2413  return table[i];
2414  }
2415 
2416  case 3:
2417  {
2418  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2419  TableIndices<2>(1, 1),
2420  TableIndices<2>(2, 2),
2421  TableIndices<2>(0, 1),
2422  TableIndices<2>(0, 2),
2423  TableIndices<2>(1, 2)};
2424  return table[i];
2425  }
2426 
2427  default:
2428  if (i < dim)
2429  return {i, i};
2430 
2431  for (unsigned int d = 0, c = dim; d < dim; ++d)
2432  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2433  if (c == i)
2434  return {d, e};
2435 
2436  // should never get here:
2437  Assert(false, ExcInternalError());
2438  return {0, 0};
2439  }
2440  }
2441 
2442  // a function to do the inverse of the unrolling from a set of
2443  // indices to a scalar index into the array in which we store
2444  // the elements of a symmetric tensor. in other words, it goes
2445  // from the scalar index into the array to a set of indices of
2446  // the tensor
2447  //
2448  // this function is for tensors of a rank not already handled
2449  // above
2450  template <int dim, int rank_>
2451  DEAL_II_CONSTEXPR inline
2452  typename std::enable_if<rank_ != 2, TableIndices<rank_>>::type
2453  unrolled_to_component_indices(const unsigned int i,
2454  const std::integral_constant<int, rank_> &)
2455  {
2456  (void)i;
2457  Assert(
2458  (i <
2460  ExcIndexRange(i,
2461  0,
2463  n_independent_components));
2464  Assert(false, ExcNotImplemented());
2465  return TableIndices<rank_>();
2466  }
2467 
2468  } // namespace SymmetricTensorImplementation
2469 } // namespace internal
2470 
2471 template <int rank_, int dim, typename Number>
2474  const unsigned int i)
2475 {
2476  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2477  dim>(i, std::integral_constant<int, rank_>());
2478 }
2479 
2480 
2481 
2482 template <int rank_, int dim, typename Number>
2483 template <class Archive>
2484 inline void
2485 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2486 {
2487  ar &data;
2488 }
2489 
2490 
2491 #endif // DOXYGEN
2492 
2493 /* ----------------- Non-member functions operating on tensors. ------------ */
2494 
2495 
2508 template <int rank_, int dim, typename Number, typename OtherNumber>
2513 {
2515  tmp = left;
2516  tmp += right;
2517  return tmp;
2518 }
2519 
2520 
2533 template <int rank_, int dim, typename Number, typename OtherNumber>
2538 {
2540  tmp = left;
2541  tmp -= right;
2542  return tmp;
2543 }
2544 
2545 
2553 template <int rank_, int dim, typename Number, typename OtherNumber>
2554 constexpr DEAL_II_ALWAYS_INLINE
2557  const Tensor<rank_, dim, OtherNumber> & right)
2558 {
2559  return Tensor<rank_, dim, Number>(left) + right;
2560 }
2561 
2562 
2570 template <int rank_, int dim, typename Number, typename OtherNumber>
2571 constexpr DEAL_II_ALWAYS_INLINE
2575 {
2576  return left + Tensor<rank_, dim, OtherNumber>(right);
2577 }
2578 
2579 
2587 template <int rank_, int dim, typename Number, typename OtherNumber>
2588 constexpr DEAL_II_ALWAYS_INLINE
2591  const Tensor<rank_, dim, OtherNumber> & right)
2592 {
2593  return Tensor<rank_, dim, Number>(left) - right;
2594 }
2595 
2596 
2604 template <int rank_, int dim, typename Number, typename OtherNumber>
2605 constexpr DEAL_II_ALWAYS_INLINE
2609 {
2610  return left - Tensor<rank_, dim, OtherNumber>(right);
2611 }
2612 
2613 
2614 
2628 template <int dim, typename Number>
2631 {
2632  switch (dim)
2633  {
2634  case 1:
2635  return t.data[0];
2636  case 2:
2637  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2638  case 3:
2639  {
2640  // in analogy to general tensors, but
2641  // there's something to be simplified for
2642  // the present case
2643  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2644  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2645  t.data[0] * t.data[5] * t.data[5] -
2646  t.data[1] * t.data[4] * t.data[4] -
2647  t.data[2] * t.data[3] * t.data[3]);
2648  }
2649  default:
2650  Assert(false, ExcNotImplemented());
2652  }
2653 }
2654 
2655 
2656 
2668 template <int dim, typename Number>
2669 constexpr DEAL_II_ALWAYS_INLINE Number
2671 {
2672  return determinant(t);
2673 }
2674 
2675 
2676 
2686 template <int dim, typename Number>
2689 {
2690  Number t = d.data[0];
2691  for (unsigned int i = 1; i < dim; ++i)
2692  t += d.data[i];
2693  return t;
2694 }
2695 
2696 
2708 template <int dim, typename Number>
2709 constexpr Number
2711 {
2712  return trace(t);
2713 }
2714 
2715 
2727 template <typename Number>
2728 constexpr DEAL_II_ALWAYS_INLINE Number
2730 {
2732 }
2733 
2734 
2735 
2754 template <typename Number>
2755 constexpr DEAL_II_ALWAYS_INLINE Number
2757 {
2758  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2759 }
2760 
2761 
2762 
2771 template <typename Number>
2772 constexpr DEAL_II_ALWAYS_INLINE Number
2774 {
2775  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2776  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2777 }
2778 
2779 
2780 
2788 template <typename Number>
2789 std::array<Number, 1>
2791 
2792 
2793 
2816 template <typename Number>
2817 std::array<Number, 2>
2819 
2820 
2821 
2844 template <typename Number>
2845 std::array<Number, 3>
2847 
2848 
2849 
2850 namespace internal
2851 {
2852  namespace SymmetricTensorImplementation
2853  {
2891  template <int dim, typename Number>
2892  void
2893  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2894  ::Tensor<2, dim, Number> & Q,
2895  std::array<Number, dim> & d,
2896  std::array<Number, dim - 1> & e);
2897 
2898 
2899 
2939  template <int dim, typename Number>
2940  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2941  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2942 
2943 
2944 
2984  template <int dim, typename Number>
2985  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2986  jacobi(::SymmetricTensor<2, dim, Number> A);
2987 
2988 
2989 
3003  template <typename Number>
3004  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3005  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3006 
3007 
3008 
3041  template <typename Number>
3042  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3043  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3044 
3049  template <int dim, typename Number>
3051  {
3052  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3053  bool
3054  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3055  {
3056  return lhs.first > rhs.first;
3057  }
3058  };
3059 
3060  } // namespace SymmetricTensorImplementation
3061 
3062 } // namespace internal
3063 
3064 
3065 
3066 // The line below is to ensure that doxygen puts the full description
3067 // of this global enumeration into the documentation
3068 // See https://stackoverflow.com/a/1717984
3098 {
3108  hybrid,
3118  ql_implicit_shifts,
3126  jacobi
3127 };
3128 
3129 
3130 
3159 template <int dim, typename Number>
3160 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3161  std::integral_constant<int, dim>::value>
3163  const SymmetricTensorEigenvectorMethod method =
3165 
3166 
3167 
3176 template <int rank_, int dim, typename Number>
3179 {
3180  return t;
3181 }
3182 
3183 
3184 
3195 template <int dim, typename Number>
3198 {
3200 
3201  // subtract scaled trace from the diagonal
3202  const Number tr = trace(t) / dim;
3203  for (unsigned int i = 0; i < dim; ++i)
3204  tmp.data[i] -= tr;
3205 
3206  return tmp;
3207 }
3208 
3209 
3210 
3217 template <int dim, typename Number>
3220 {
3221  // create a default constructed matrix filled with
3222  // zeros, then set the diagonal elements to one
3224  switch (dim)
3225  {
3226  case 1:
3228  break;
3229  case 2:
3230  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3231  break;
3232  case 3:
3233  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3235  break;
3236  default:
3237  for (unsigned int d = 0; d < dim; ++d)
3239  }
3240  return tmp;
3241 }
3242 
3243 
3244 
3252 template <int dim>
3255 {
3256  return unit_symmetric_tensor<dim, double>();
3257 }
3258 
3259 
3260 
3289 template <int dim, typename Number>
3292 {
3294 
3295  // fill the elements treating the diagonal
3296  for (unsigned int i = 0; i < dim; ++i)
3297  for (unsigned int j = 0; j < dim; ++j)
3298  tmp.data[i][j] =
3299  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3300 
3301  // then fill the ones that copy over the
3302  // non-diagonal elements. note that during
3303  // the double-contraction, we handle the
3304  // off-diagonal elements twice, so simply
3305  // copying requires a weight of 1/2
3306  for (unsigned int i = dim;
3307  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3308  n_rank2_components;
3309  ++i)
3310  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3311 
3312  return tmp;
3313 }
3314 
3315 
3316 
3324 template <int dim>
3327 {
3328  return deviator_tensor<dim, double>();
3329 }
3330 
3331 
3332 
3370 template <int dim, typename Number>
3373 {
3375 
3376  // fill the elements treating the diagonal
3377  for (unsigned int i = 0; i < dim; ++i)
3378  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3379 
3380  // then fill the ones that copy over the
3381  // non-diagonal elements. note that during
3382  // the double-contraction, we handle the
3383  // off-diagonal elements twice, so simply
3384  // copying requires a weight of 1/2
3385  for (unsigned int i = dim;
3386  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3387  n_rank2_components;
3388  ++i)
3389  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3390 
3391  return tmp;
3392 }
3393 
3394 
3395 
3403 template <int dim>
3406 {
3407  return identity_tensor<dim, double>();
3408 }
3409 
3410 
3411 
3421 template <int dim, typename Number>
3424 {
3426  value(t);
3427 }
3428 
3429 
3430 
3441 template <int dim, typename Number>
3444 {
3446  value(t);
3447 }
3448 
3449 
3450 
3472 template <int dim, typename Number>
3476 {
3478 
3479  // fill only the elements really needed
3480  for (unsigned int i = 0; i < dim; ++i)
3481  for (unsigned int j = i; j < dim; ++j)
3482  for (unsigned int k = 0; k < dim; ++k)
3483  for (unsigned int l = k; l < dim; ++l)
3484  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3485 
3486  return tmp;
3487 }
3488 
3489 
3490 
3498 template <int dim, typename Number>
3501 {
3503  for (unsigned int d = 0; d < dim; ++d)
3504  result[d][d] = t[d][d];
3505 
3506  const Number half = internal::NumberType<Number>::value(0.5);
3507  for (unsigned int d = 0; d < dim; ++d)
3508  for (unsigned int e = d + 1; e < dim; ++e)
3509  result[d][e] = (t[d][e] + t[e][d]) * half;
3510  return result;
3511 }
3512 
3513 
3514 
3522 template <int rank_, int dim, typename Number>
3525  operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3526 {
3528  tt *= factor;
3529  return tt;
3530 }
3531 
3532 
3533 
3541 template <int rank_, int dim, typename Number>
3543  operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3544 {
3545  // simply forward to the other operator
3546  return t * factor;
3547 }
3548 
3549 
3550 
3576 template <int rank_, int dim, typename Number, typename OtherNumber>
3578  rank_,
3579  dim,
3580  typename ProductType<Number,
3583  const OtherNumber & factor)
3584 {
3585  // form the product. we have to convert the two factors into the final
3586  // type via explicit casts because, for awkward reasons, the C++
3587  // standard committee saw it fit to not define an
3588  // operator*(float,std::complex<double>)
3589  // (as well as with switched arguments and double<->float).
3590  using product_type = typename ProductType<Number, OtherNumber>::type;
3593  return tt;
3594 }
3595 
3596 
3597 
3606 template <int rank_, int dim, typename Number, typename OtherNumber>
3608  rank_,
3609  dim,
3610  typename ProductType<OtherNumber,
3612 operator*(const Number & factor,
3614 {
3615  // simply forward to the other operator with switched arguments
3616  return (t * factor);
3617 }
3618 
3619 
3620 
3626 template <int rank_, int dim, typename Number, typename OtherNumber>
3628  rank_,
3629  dim,
3630  typename ProductType<Number,
3633  const OtherNumber & factor)
3634 {
3635  using product_type = typename ProductType<Number, OtherNumber>::type;
3638  return tt;
3639 }
3640 
3641 
3642 
3649 template <int rank_, int dim>
3651  operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3652 {
3654  tt *= factor;
3655  return tt;
3656 }
3657 
3658 
3659 
3666 template <int rank_, int dim>
3668  operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3669 {
3671  tt *= factor;
3672  return tt;
3673 }
3674 
3675 
3676 
3682 template <int rank_, int dim>
3684 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3685 {
3687  tt /= factor;
3688  return tt;
3689 }
3690 
3700 template <int dim, typename Number, typename OtherNumber>
3704 {
3705  return (t1 * t2);
3706 }
3707 
3708 
3722 template <int dim, typename Number, typename OtherNumber>
3726  const Tensor<2, dim, OtherNumber> & t2)
3727 {
3730  for (unsigned int i = 0; i < dim; ++i)
3731  for (unsigned int j = 0; j < dim; ++j)
3732  s += t1[i][j] * t2[i][j];
3733  return s;
3734 }
3735 
3736 
3750 template <int dim, typename Number, typename OtherNumber>
3754 {
3755  return scalar_product(t2, t1);
3756 }
3757 
3758 
3773 template <typename Number, typename OtherNumber>
3778 {
3779  tmp[0][0] = t[0][0][0][0] * s[0][0];
3780 }
3781 
3782 
3783 
3798 template <typename Number, typename OtherNumber>
3803 {
3804  tmp[0][0] = t[0][0][0][0] * s[0][0];
3805 }
3806 
3807 
3808 
3823 template <typename Number, typename OtherNumber>
3828 {
3829  const unsigned int dim = 2;
3830 
3831  for (unsigned int i = 0; i < dim; ++i)
3832  for (unsigned int j = i; j < dim; ++j)
3833  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3834  2 * t[i][j][0][1] * s[0][1];
3835 }
3836 
3837 
3838 
3853 template <typename Number, typename OtherNumber>
3858 {
3859  const unsigned int dim = 2;
3860 
3861  for (unsigned int i = 0; i < dim; ++i)
3862  for (unsigned int j = i; j < dim; ++j)
3863  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3864  2 * s[0][1] * t[0][1][i][j];
3865 }
3866 
3867 
3868 
3883 template <typename Number, typename OtherNumber>
3888 {
3889  const unsigned int dim = 3;
3890 
3891  for (unsigned int i = 0; i < dim; ++i)
3892  for (unsigned int j = i; j < dim; ++j)
3893  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3894  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3895  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3896 }
3897 
3898 
3899 
3914 template <typename Number, typename OtherNumber>
3919 {
3920  const unsigned int dim = 3;
3921 
3922  for (unsigned int i = 0; i < dim; ++i)
3923  for (unsigned int j = i; j < dim; ++j)
3924  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3925  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3926  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3927 }
3928 
3929 
3930 
3937 template <int dim, typename Number, typename OtherNumber>
3941  const Tensor<1, dim, OtherNumber> & src2)
3942 {
3944  for (unsigned int i = 0; i < dim; ++i)
3945  for (unsigned int j = 0; j < dim; ++j)
3946  dest[i] += src1[i][j] * src2[j];
3947  return dest;
3948 }
3949 
3950 
3957 template <int dim, typename Number, typename OtherNumber>
3961 {
3962  // this is easy for symmetric tensors:
3963  return src2 * src1;
3964 }
3965 
3966 
3967 
3987 template <int rank_1,
3988  int rank_2,
3989  int dim,
3990  typename Number,
3991  typename OtherNumber>
3992 constexpr DEAL_II_ALWAYS_INLINE
3993  typename Tensor<rank_1 + rank_2 - 2,
3994  dim,
3995  typename ProductType<Number, OtherNumber>::type>::tensor_type
3998 {
3999  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4000 }
4001 
4002 
4003 
4023 template <int rank_1,
4024  int rank_2,
4025  int dim,
4026  typename Number,
4027  typename OtherNumber>
4028 constexpr DEAL_II_ALWAYS_INLINE
4029  typename Tensor<rank_1 + rank_2 - 2,
4030  dim,
4031  typename ProductType<Number, OtherNumber>::type>::tensor_type
4033  const Tensor<rank_2, dim, OtherNumber> & src2)
4034 {
4035  return Tensor<rank_1, dim, Number>(src1) * src2;
4036 }
4037 
4038 
4039 
4049 template <int dim, typename Number>
4050 inline std::ostream &
4051 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4052 {
4053  // make our lives a bit simpler by outputting
4054  // the tensor through the operator for the
4055  // general Tensor class
4057 
4058  for (unsigned int i = 0; i < dim; ++i)
4059  for (unsigned int j = 0; j < dim; ++j)
4060  tt[i][j] = t[i][j];
4061 
4062  return out << tt;
4063 }
4064 
4065 
4066 
4076 template <int dim, typename Number>
4077 inline std::ostream &
4078 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4079 {
4080  // make our lives a bit simpler by outputting
4081  // the tensor through the operator for the
4082  // general Tensor class
4084 
4085  for (unsigned int i = 0; i < dim; ++i)
4086  for (unsigned int j = 0; j < dim; ++j)
4087  for (unsigned int k = 0; k < dim; ++k)
4088  for (unsigned int l = 0; l < dim; ++l)
4089  tt[i][j][k][l] = t[i][j][k][l];
4090 
4091  return out << tt;
4092 }
4093 
4094 
4096 
4097 #endif
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:372
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:448
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
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:154
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)