Reference documentation for deal.II version Git a2518c4284 2020-04-02 19:23:00 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\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 - 2019 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>
23 #include <deal.II/base/table_indices.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/base/tensor.h>
26 
27 #include <algorithm>
28 #include <array>
29 #include <functional>
30 
31 DEAL_II_NAMESPACE_OPEN
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>
40 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
42 
43 template <int dim, typename Number>
44 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
46 
47 template <int dim, typename Number>
48 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
50 
51 template <int dim, typename Number>
52 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
54 
55 template <int dim, typename Number>
56 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
58 
59 template <int dim2, typename Number>
60 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
62 
63 template <int dim, typename Number>
64 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
66 
67 template <int dim, typename Number>
68 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
70 
71 
72 
73 namespace internal
74 {
75  template <int rank, int dim, typename T, typename U>
76  struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
77  {
78  using type =
79  SymmetricTensor<rank,
80  dim,
81  std::complex<typename ProductType<T, U>::type>>;
82  };
83 
84  template <typename T, int rank, int dim, typename U>
85  struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
86  {
87  using type =
88  SymmetricTensor<rank,
89  dim,
90  std::complex<typename ProductType<T, U>::type>>;
91  };
92 
97  namespace SymmetricTensorImplementation
98  {
103  template <int rank, int dim, typename Number>
104  struct Inverse;
105  } // namespace SymmetricTensorImplementation
106 
111  namespace SymmetricTensorAccessors
112  {
119  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
120  merge(const TableIndices<2> &previous_indices,
121  const unsigned int new_index,
122  const unsigned int position)
123  {
124  AssertIndexRange(position, 2);
125 
126  if (position == 0)
127  return {new_index, numbers::invalid_unsigned_int};
128  else
129  return {previous_indices[0], new_index};
130  }
131 
132 
133 
140  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<4>
141  merge(const TableIndices<4> &previous_indices,
142  const unsigned int new_index,
143  const unsigned int position)
144  {
145  AssertIndexRange(position, 4);
146 
147  switch (position)
148  {
149  case 0:
150  return {new_index,
153  numbers::invalid_unsigned_int};
154  case 1:
155  return {previous_indices[0],
156  new_index,
158  numbers::invalid_unsigned_int};
159  case 2:
160  return {previous_indices[0],
161  previous_indices[1],
162  new_index,
163  numbers::invalid_unsigned_int};
164  case 3:
165  return {previous_indices[0],
166  previous_indices[1],
167  previous_indices[2],
168  new_index};
169  default:
170  Assert(false, ExcInternalError());
171  return {};
172  }
173  }
174 
175 
184  template <int rank1,
185  int rank2,
186  int dim,
187  typename Number,
188  typename OtherNumber = Number>
190  {
191  using value_type = typename ProductType<Number, OtherNumber>::type;
192  using type =
193  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
194  };
195 
196 
205  template <int dim, typename Number, typename OtherNumber>
206  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
207  {
208  using type = typename ProductType<Number, OtherNumber>::type;
209  };
210 
211 
212 
225  template <int rank, int dim, typename Number>
226  struct StorageType;
227 
231  template <int dim, typename Number>
232  struct StorageType<2, dim, Number>
233  {
238  static const unsigned int n_independent_components =
239  (dim * dim + dim) / 2;
240 
245  };
246 
247 
248 
252  template <int dim, typename Number>
253  struct StorageType<4, dim, Number>
254  {
260  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
261 
265  static const unsigned int n_independent_components =
266  (n_rank2_components *
268 
276  };
277 
278 
279 
284  template <int rank, int dim, bool constness, typename Number>
286 
293  template <int rank, int dim, typename Number>
294  struct AccessorTypes<rank, dim, true, Number>
295  {
296  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
297 
298  using reference = Number;
299  };
300 
307  template <int rank, int dim, typename Number>
308  struct AccessorTypes<rank, dim, false, Number>
309  {
311 
312  using reference = Number &;
313  };
314 
315 
350  template <int rank, int dim, bool constness, int P, typename Number>
351  class Accessor
352  {
353  public:
357  using reference =
359  using tensor_type =
361 
362  private:
381  constexpr Accessor(tensor_type & tensor,
382  const TableIndices<rank> &previous_indices);
383 
387  constexpr DEAL_II_ALWAYS_INLINE
388  Accessor(const Accessor &) = default;
389 
390  public:
394  DEAL_II_CONSTEXPR Accessor<rank, dim, constness, P - 1, Number>
395  operator[](const unsigned int i);
396 
400  constexpr Accessor<rank, dim, constness, P - 1, Number>
401  operator[](const unsigned int i) const;
402 
403  private:
407  tensor_type & tensor;
408  const TableIndices<rank> previous_indices;
409 
410  // Declare some other classes as friends. Make sure to work around bugs
411  // in some compilers:
412  template <int, int, typename>
413  friend class ::SymmetricTensor;
414  template <int, int, bool, int, typename>
415  friend class Accessor;
416 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
417  friend class ::SymmetricTensor<rank, dim, Number>;
418  friend class Accessor<rank, dim, constness, P + 1, Number>;
419 #endif
420  };
421 
422 
423 
433  template <int rank, int dim, bool constness, typename Number>
434  class Accessor<rank, dim, constness, 1, Number>
435  {
436  public:
440  using reference =
442  using tensor_type =
444 
445  private:
467  constexpr Accessor(tensor_type & tensor,
468  const TableIndices<rank> &previous_indices);
469 
473  constexpr DEAL_II_ALWAYS_INLINE
474  Accessor(const Accessor &) = default;
475 
476  public:
480  DEAL_II_CONSTEXPR reference operator[](const unsigned int);
481 
485  constexpr reference operator[](const unsigned int) const;
486 
487  private:
491  tensor_type & tensor;
492  const TableIndices<rank> previous_indices;
493 
494  // Declare some other classes as friends. Make sure to work around bugs
495  // in some compilers:
496  template <int, int, typename>
497  friend class ::SymmetricTensor;
498  template <int, int, bool, int, typename>
499  friend class SymmetricTensorAccessors::Accessor;
500 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
501  friend class ::SymmetricTensor<rank, dim, Number>;
502  friend class SymmetricTensorAccessors::
503  Accessor<rank, dim, constness, 2, Number>;
504 #endif
505  };
506  } // namespace SymmetricTensorAccessors
507 } // namespace internal
508 
509 
510 
574 template <int rank_, int dim, typename Number>
576 {
577 public:
578  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
579 
588  static const unsigned int dimension = dim;
589 
593  static const unsigned int rank = rank_;
594 
600  static constexpr unsigned int n_independent_components =
602  n_independent_components;
603 
607  constexpr DEAL_II_ALWAYS_INLINE
608  SymmetricTensor() = default;
609 
623  template <typename OtherNumber>
624  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
625 
641  DEAL_II_CONSTEXPR
642  SymmetricTensor(const Number (&array)[n_independent_components]);
643 
649  template <typename OtherNumber>
650  constexpr explicit SymmetricTensor(
651  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
652 
656  Number *
657  begin_raw();
658 
662  const Number *
663  begin_raw() const;
664 
668  Number *
669  end_raw();
670 
675  const Number *
676  end_raw() const;
677 
684  template <typename OtherNumber>
685  DEAL_II_CONSTEXPR SymmetricTensor &
686  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
687 
694  DEAL_II_CONSTEXPR SymmetricTensor &
695  operator=(const Number &d);
696 
701  constexpr operator Tensor<rank_, dim, Number>() const;
702 
706  constexpr bool
707  operator==(const SymmetricTensor &) const;
708 
712  constexpr bool
713  operator!=(const SymmetricTensor &) const;
714 
718  template <typename OtherNumber>
719  DEAL_II_CONSTEXPR SymmetricTensor &
720  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
721 
725  template <typename OtherNumber>
726  DEAL_II_CONSTEXPR SymmetricTensor &
727  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
728 
733  template <typename OtherNumber>
734  DEAL_II_CONSTEXPR SymmetricTensor &
735  operator*=(const OtherNumber &factor);
736 
740  template <typename OtherNumber>
741  DEAL_II_CONSTEXPR SymmetricTensor &
742  operator/=(const OtherNumber &factor);
743 
747  DEAL_II_CONSTEXPR SymmetricTensor
748  operator-() const;
749 
774  template <typename OtherNumber>
775  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
776  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
778 
783  template <typename OtherNumber>
784  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
785  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
787 
791  DEAL_II_CONSTEXPR Number &
792  operator()(const TableIndices<rank_> &indices);
793 
797  DEAL_II_CONSTEXPR const Number &
798  operator()(const TableIndices<rank_> &indices) const;
799 
804  constexpr internal::SymmetricTensorAccessors::
805  Accessor<rank_, dim, true, rank_ - 1, Number>
806  operator[](const unsigned int row) const;
807 
812  DEAL_II_CONSTEXPR internal::SymmetricTensorAccessors::
813  Accessor<rank_, dim, false, rank_ - 1, Number>
814  operator[](const unsigned int row);
815 
821  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
822 
828  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
829 
835  DEAL_II_CONSTEXPR const Number &
836  access_raw_entry(const unsigned int unrolled_index) const;
837 
843  DEAL_II_CONSTEXPR Number &
844  access_raw_entry(const unsigned int unrolled_index);
845 
856  norm() const;
857 
865  static constexpr unsigned int
866  component_to_unrolled_index(const TableIndices<rank_> &indices);
867 
873  static constexpr TableIndices<rank_>
874  unrolled_to_component_indices(const unsigned int i);
875 
888  DEAL_II_CONSTEXPR void
889  clear();
890 
895  static constexpr std::size_t
896  memory_consumption();
897 
902  template <class Archive>
903  void
904  serialize(Archive &ar, const unsigned int version);
905 
906 private:
910  using base_tensor_descriptor =
912 
916  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
917 
922 
923  // Make all other symmetric tensors friends.
924  template <int, int, typename>
925  friend class SymmetricTensor;
926 
927  // Make a few more functions friends.
928  template <int dim2, typename Number2>
929  friend DEAL_II_CONSTEXPR Number2
931 
932  template <int dim2, typename Number2>
933  friend DEAL_II_CONSTEXPR Number2
935 
936  template <int dim2, typename Number2>
937  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
939 
940  template <int dim2, typename Number2>
941  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
943 
944  template <int dim2, typename Number2>
945  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
946  deviator_tensor();
947 
948  template <int dim2, typename Number2>
949  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
950  identity_tensor();
951 
952 
953  // Make a few helper classes friends as well.
955  Inverse<2, dim, Number>;
956 
958  Inverse<4, dim, Number>;
959 };
960 
961 
962 
963 // ------------------------- inline functions ------------------------
964 
965 #ifndef DOXYGEN
966 
967 // provide declarations for static members
968 template <int rank, int dim, typename Number>
969 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
970 
971 template <int rank_, int dim, typename Number>
972 constexpr unsigned int
973  SymmetricTensor<rank_, dim, Number>::n_independent_components;
974 
975 namespace internal
976 {
977  namespace SymmetricTensorAccessors
978  {
979  template <int rank_, int dim, bool constness, int P, typename Number>
980  constexpr DEAL_II_ALWAYS_INLINE
981  Accessor<rank_, dim, constness, P, Number>::Accessor(
982  tensor_type & tensor,
983  const TableIndices<rank_> &previous_indices)
984  : tensor(tensor)
985  , previous_indices(previous_indices)
986  {}
987 
988 
989 
990  template <int rank_, int dim, bool constness, int P, typename Number>
991  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
992  Accessor<rank_, dim, constness, P - 1, Number>
993  Accessor<rank_, dim, constness, P, Number>::
994  operator[](const unsigned int i)
995  {
996  return Accessor<rank_, dim, constness, P - 1, Number>(
997  tensor, merge(previous_indices, i, rank_ - P));
998  }
999 
1000 
1001 
1002  template <int rank_, int dim, bool constness, int P, typename Number>
1003  constexpr DEAL_II_ALWAYS_INLINE
1004  Accessor<rank_, dim, constness, P - 1, Number>
1005  Accessor<rank_, dim, constness, P, Number>::
1006  operator[](const unsigned int i) const
1007  {
1008  return Accessor<rank_, dim, constness, P - 1, Number>(
1009  tensor, merge(previous_indices, i, rank_ - P));
1010  }
1011 
1012 
1013 
1014  template <int rank_, int dim, bool constness, typename Number>
1015  constexpr DEAL_II_ALWAYS_INLINE
1016  Accessor<rank_, dim, constness, 1, Number>::Accessor(
1017  tensor_type & tensor,
1018  const TableIndices<rank_> &previous_indices)
1019  : tensor(tensor)
1020  , previous_indices(previous_indices)
1021  {}
1022 
1023 
1024 
1025  template <int rank_, int dim, bool constness, typename Number>
1026  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1027  typename Accessor<rank_, dim, constness, 1, Number>::reference
1028  Accessor<rank_, dim, constness, 1, Number>::
1029  operator[](const unsigned int i)
1030  {
1031  return tensor(merge(previous_indices, i, rank_ - 1));
1032  }
1033 
1034 
1035  template <int rank_, int dim, bool constness, typename Number>
1036  constexpr DEAL_II_ALWAYS_INLINE
1037  typename Accessor<rank_, dim, constness, 1, Number>::reference
1038  Accessor<rank_, dim, constness, 1, Number>::
1039  operator[](const unsigned int i) const
1040  {
1041  return tensor(merge(previous_indices, i, rank_ - 1));
1042  }
1043  } // namespace SymmetricTensorAccessors
1044 } // namespace internal
1045 
1046 
1047 
1048 template <int rank_, int dim, typename Number>
1049 template <typename OtherNumber>
1050 inline DEAL_II_ALWAYS_INLINE
1052  const Tensor<2, dim, OtherNumber> &t)
1053 {
1054  static_assert(rank == 2, "This function is only implemented for rank==2");
1055  for (unsigned int d = 0; d < dim; ++d)
1056  for (unsigned int e = 0; e < d; ++e)
1057  Assert(t[d][e] == t[e][d], ExcInternalError());
1058 
1059  for (unsigned int d = 0; d < dim; ++d)
1060  data[d] = t[d][d];
1061 
1062  for (unsigned int d = 0, c = 0; d < dim; ++d)
1063  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1064  data[dim + c] = t[d][e];
1065 }
1066 
1067 
1068 
1069 template <int rank_, int dim, typename Number>
1070 template <typename OtherNumber>
1071 constexpr DEAL_II_ALWAYS_INLINE
1073  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1074  : data(initializer.data)
1075 {}
1076 
1077 
1078 
1079 template <int rank_, int dim, typename Number>
1080 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1082  const Number (&array)[n_independent_components])
1083  : data(
1084  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1085 {
1086  // ensure that the reinterpret_cast above actually works
1087  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1088  ExcInternalError());
1089 }
1090 
1091 
1092 
1093 template <int rank_, int dim, typename Number>
1094 template <typename OtherNumber>
1095 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1099 {
1100  data = t.data;
1101  return *this;
1102 }
1103 
1104 
1105 
1106 template <int rank_, int dim, typename Number>
1107 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1110 {
1112  ExcMessage("Only assignment with zero is allowed"));
1113  (void)d;
1114 
1115  data = internal::NumberType<Number>::value(0.0);
1116 
1117  return *this;
1118 }
1119 
1120 
1121 namespace internal
1122 {
1123  namespace SymmetricTensorImplementation
1124  {
1125  template <int dim, typename Number>
1126  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1127  ::Tensor<2, dim, Number>
1128  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1129  {
1131 
1132  // diagonal entries are stored first
1133  for (unsigned int d = 0; d < dim; ++d)
1134  t[d][d] = s.access_raw_entry(d);
1135 
1136  // off-diagonal entries come next, row by row
1137  for (unsigned int d = 0, c = 0; d < dim; ++d)
1138  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1139  {
1140  t[d][e] = s.access_raw_entry(dim + c);
1141  t[e][d] = s.access_raw_entry(dim + c);
1142  }
1143  return t;
1144  }
1145 
1146 
1147  template <int dim, typename Number>
1148  DEAL_II_CONSTEXPR ::Tensor<4, dim, Number>
1149  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1150  {
1151  // utilize the symmetry properties of SymmetricTensor<4,dim>
1152  // discussed in the class documentation to avoid accessing all
1153  // independent elements of the input tensor more than once
1155 
1156  for (unsigned int i = 0; i < dim; ++i)
1157  for (unsigned int j = i; j < dim; ++j)
1158  for (unsigned int k = 0; k < dim; ++k)
1159  for (unsigned int l = k; l < dim; ++l)
1160  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1161  t[TableIndices<4>(j, i, k, l)] =
1162  t[TableIndices<4>(j, i, l, k)] =
1163  st[TableIndices<4>(i, j, k, l)];
1164 
1165  return t;
1166  }
1167 
1168 
1169  template <typename Number>
1170  struct Inverse<2, 1, Number>
1171  {
1172  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1173  ::SymmetricTensor<2, 1, Number>
1174  value(const ::SymmetricTensor<2, 1, Number> &t)
1175  {
1177 
1178  tmp[0][0] = 1.0 / t[0][0];
1179 
1180  return tmp;
1181  }
1182  };
1183 
1184 
1185  template <typename Number>
1186  struct Inverse<2, 2, Number>
1187  {
1188  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1189  ::SymmetricTensor<2, 2, Number>
1190  value(const ::SymmetricTensor<2, 2, Number> &t)
1191  {
1193 
1194  // Sympy result: ([
1195  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1196  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1197  const TableIndices<2> idx_00(0, 0);
1198  const TableIndices<2> idx_01(0, 1);
1199  const TableIndices<2> idx_11(1, 1);
1200  const Number inv_det_t =
1201  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1202  tmp[idx_00] = t[idx_11];
1203  tmp[idx_01] = -t[idx_01];
1204  tmp[idx_11] = t[idx_00];
1205  tmp *= inv_det_t;
1206 
1207  return tmp;
1208  }
1209  };
1210 
1211 
1212  template <typename Number>
1213  struct Inverse<2, 3, Number>
1214  {
1215  DEAL_II_CONSTEXPR static ::SymmetricTensor<2, 3, Number>
1216  value(const ::SymmetricTensor<2, 3, Number> &t)
1217  {
1219 
1220  // Sympy result: ([
1221  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1222  // 2*t01*t02*t12 - t02**2*t11),
1223  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1224  // 2*t01*t02*t12 - t02**2*t11),
1225  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1226  // 2*t01*t02*t12 - t02**2*t11)],
1227  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1228  // 2*t01*t02*t12 - t02**2*t11),
1229  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1230  // 2*t01*t02*t12 - t02**2*t11),
1231  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1232  // 2*t01*t02*t12 + t02**2*t11)],
1233  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1234  // 2*t01*t02*t12 - t02**2*t11),
1235  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1236  // 2*t01*t02*t12 + t02**2*t11),
1237  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1238  // 2*t01*t02*t12 + t02**2*t11)] ])
1239  //
1240  // =
1241  //
1242  // [ (t11*t22 - t12**2)/det_t,
1243  // (-t01*t22 + t02*t12)/det_t,
1244  // (t01*t12 - t02*t11)/det_t],
1245  // [ (-t01*t22 + t02*t12)/det_t,
1246  // (t00*t22 - t02**2)/det_t,
1247  // (-t00*t12 + t01*t02)/det_t],
1248  // [ (t01*t12 - t02*t11)/det_t,
1249  // (-t00*t12 + t01*t02)/det_t,
1250  // (t00*t11 - t01**2)/det_t] ])
1251  //
1252  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1253  // 2*t01*t02*t12 - t02**2*t11)
1254  const TableIndices<2> idx_00(0, 0);
1255  const TableIndices<2> idx_01(0, 1);
1256  const TableIndices<2> idx_02(0, 2);
1257  const TableIndices<2> idx_11(1, 1);
1258  const TableIndices<2> idx_12(1, 2);
1259  const TableIndices<2> idx_22(2, 2);
1260  const Number inv_det_t =
1261  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1262  t[idx_00] * t[idx_12] * t[idx_12] -
1263  t[idx_01] * t[idx_01] * t[idx_22] +
1264  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1265  t[idx_02] * t[idx_02] * t[idx_11]);
1266  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1267  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1268  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1269  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1270  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1271  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1272  tmp *= inv_det_t;
1273 
1274  return tmp;
1275  }
1276  };
1277 
1278 
1279  template <typename Number>
1280  struct Inverse<4, 1, Number>
1281  {
1282  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 1, Number>
1283  value(const ::SymmetricTensor<4, 1, Number> &t)
1284  {
1286  tmp.data[0][0] = 1.0 / t.data[0][0];
1287  return tmp;
1288  }
1289  };
1290 
1291 
1292  template <typename Number>
1293  struct Inverse<4, 2, Number>
1294  {
1295  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 2, Number>
1296  value(const ::SymmetricTensor<4, 2, Number> &t)
1297  {
1299 
1300  // Inverting this tensor is a little more complicated than necessary,
1301  // since we store the data of 't' as a 3x3 matrix t.data, but the
1302  // product between a rank-4 and a rank-2 tensor is really not the
1303  // product between this matrix and the 3-vector of a rhs, but rather
1304  //
1305  // B.vec = t.data * mult * A.vec
1306  //
1307  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1308  // capture the fact that we need to add up both the c_ij12*a_12 and the
1309  // c_ij21*a_21 terms.
1310  //
1311  // In addition, in this scheme, the identity tensor has the matrix
1312  // representation mult^-1.
1313  //
1314  // The inverse of 't' therefore has the matrix representation
1315  //
1316  // inv.data = mult^-1 * t.data^-1 * mult^-1
1317  //
1318  // in order to compute it, let's first compute the inverse of t.data and
1319  // put it into tmp.data; at the end of the function we then scale the
1320  // last row and column of the inverse by 1/2, corresponding to the left
1321  // and right multiplication with mult^-1.
1322  const Number t4 = t.data[0][0] * t.data[1][1],
1323  t6 = t.data[0][0] * t.data[1][2],
1324  t8 = t.data[0][1] * t.data[1][0],
1325  t00 = t.data[0][2] * t.data[1][0],
1326  t01 = t.data[0][1] * t.data[2][0],
1327  t04 = t.data[0][2] * t.data[2][0],
1328  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1329  t8 * t.data[2][2] + t00 * t.data[2][1] +
1330  t01 * t.data[1][2] - t04 * t.data[1][1]);
1331  tmp.data[0][0] =
1332  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1333  tmp.data[0][1] =
1334  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1335  tmp.data[0][2] =
1336  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1337  tmp.data[1][0] =
1338  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1339  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1340  tmp.data[1][2] = -(t6 - t00) * t07;
1341  tmp.data[2][0] =
1342  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1343  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1344  tmp.data[2][2] = (t4 - t8) * t07;
1345 
1346  // scale last row and column as mentioned
1347  // above
1348  tmp.data[2][0] /= 2;
1349  tmp.data[2][1] /= 2;
1350  tmp.data[0][2] /= 2;
1351  tmp.data[1][2] /= 2;
1352  tmp.data[2][2] /= 4;
1353 
1354  return tmp;
1355  }
1356  };
1357 
1358 
1359  template <typename Number>
1360  struct Inverse<4, 3, Number>
1361  {
1362  static ::SymmetricTensor<4, 3, Number>
1363  value(const ::SymmetricTensor<4, 3, Number> &t)
1364  {
1366 
1367  // This function follows the exact same scheme as the 2d case, except
1368  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1369  // Instead, we use the Gauss-Jordan algorithm implemented for
1370  // FullMatrix. For historical reasons the following code is copied from
1371  // there, with the tangential benefit that we do not need to copy the
1372  // tensor entries to and from the FullMatrix.
1373  const unsigned int N = 6;
1374 
1375  // First get an estimate of the size of the elements of this matrix,
1376  // for later checks whether the pivot element is large enough, or
1377  // whether we have to fear that the matrix is not regular.
1378  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1379  for (unsigned int i = 0; i < N; ++i)
1380  diagonal_sum += std::fabs(tmp.data[i][i]);
1381  const Number typical_diagonal_element =
1382  diagonal_sum / static_cast<double>(N);
1383  (void)typical_diagonal_element;
1384 
1385  unsigned int p[N];
1386  for (unsigned int i = 0; i < N; ++i)
1387  p[i] = i;
1388 
1389  for (unsigned int j = 0; j < N; ++j)
1390  {
1391  // Pivot search: search that part of the line on and right of the
1392  // diagonal for the largest element.
1393  Number max = std::fabs(tmp.data[j][j]);
1394  unsigned int r = j;
1395  for (unsigned int i = j + 1; i < N; ++i)
1396  if (std::fabs(tmp.data[i][j]) > max)
1397  {
1398  max = std::fabs(tmp.data[i][j]);
1399  r = i;
1400  }
1401 
1402  // Check whether the pivot is too small
1403  Assert(max > 1.e-16 * typical_diagonal_element,
1404  ExcMessage("This tensor seems to be noninvertible"));
1405 
1406  // Row interchange
1407  if (r > j)
1408  {
1409  for (unsigned int k = 0; k < N; ++k)
1410  std::swap(tmp.data[j][k], tmp.data[r][k]);
1411 
1412  std::swap(p[j], p[r]);
1413  }
1414 
1415  // Transformation
1416  const Number hr = 1. / tmp.data[j][j];
1417  tmp.data[j][j] = hr;
1418  for (unsigned int k = 0; k < N; ++k)
1419  {
1420  if (k == j)
1421  continue;
1422  for (unsigned int i = 0; i < N; ++i)
1423  {
1424  if (i == j)
1425  continue;
1426  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1427  }
1428  }
1429  for (unsigned int i = 0; i < N; ++i)
1430  {
1431  tmp.data[i][j] *= hr;
1432  tmp.data[j][i] *= -hr;
1433  }
1434  tmp.data[j][j] = hr;
1435  }
1436 
1437  // Column interchange
1438  Number hv[N];
1439  for (unsigned int i = 0; i < N; ++i)
1440  {
1441  for (unsigned int k = 0; k < N; ++k)
1442  hv[p[k]] = tmp.data[i][k];
1443  for (unsigned int k = 0; k < N; ++k)
1444  tmp.data[i][k] = hv[k];
1445  }
1446 
1447  // Scale rows and columns. The mult matrix
1448  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1449  for (unsigned int i = 3; i < 6; ++i)
1450  for (unsigned int j = 0; j < 3; ++j)
1451  tmp.data[i][j] /= 2;
1452 
1453  for (unsigned int i = 0; i < 3; ++i)
1454  for (unsigned int j = 3; j < 6; ++j)
1455  tmp.data[i][j] /= 2;
1456 
1457  for (unsigned int i = 3; i < 6; ++i)
1458  for (unsigned int j = 3; j < 6; ++j)
1459  tmp.data[i][j] /= 4;
1460 
1461  return tmp;
1462  }
1463  };
1464 
1465  } // namespace SymmetricTensorImplementation
1466 } // namespace internal
1467 
1468 
1469 
1470 template <int rank_, int dim, typename Number>
1471 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>::
1472  operator Tensor<rank_, dim, Number>() const
1473 {
1474  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1475 }
1476 
1477 
1478 
1479 template <int rank_, int dim, typename Number>
1480 constexpr bool
1483 {
1484  return data == t.data;
1485 }
1486 
1487 
1488 
1489 template <int rank_, int dim, typename Number>
1490 constexpr bool
1493 {
1494  return data != t.data;
1495 }
1496 
1497 
1498 
1499 template <int rank_, int dim, typename Number>
1500 template <typename OtherNumber>
1501 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1505 {
1506  data += t.data;
1507  return *this;
1508 }
1509 
1510 
1511 
1512 template <int rank_, int dim, typename Number>
1513 template <typename OtherNumber>
1514 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1518 {
1519  data -= t.data;
1520  return *this;
1521 }
1522 
1523 
1524 
1525 template <int rank_, int dim, typename Number>
1526 template <typename OtherNumber>
1527 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1530 {
1531  data *= d;
1532  return *this;
1533 }
1534 
1535 
1536 
1537 template <int rank_, int dim, typename Number>
1538 template <typename OtherNumber>
1539 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1542 {
1543  data /= d;
1544  return *this;
1545 }
1546 
1547 
1548 
1549 template <int rank_, int dim, typename Number>
1550 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1553 {
1554  SymmetricTensor tmp = *this;
1555  tmp.data = -tmp.data;
1556  return tmp;
1557 }
1558 
1559 
1560 
1561 template <int rank_, int dim, typename Number>
1562 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void
1564 {
1565  data.clear();
1566 }
1567 
1568 
1569 
1570 template <int rank_, int dim, typename Number>
1571 constexpr std::size_t
1573 {
1574  // all memory consists of statically allocated memory of the current
1575  // object, no pointers
1576  return sizeof(SymmetricTensor<rank_, dim, Number>);
1577 }
1578 
1579 
1580 
1581 namespace internal
1582 {
1583  template <int dim, typename Number, typename OtherNumber = Number>
1584  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1585  typename SymmetricTensorAccessors::
1586  double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1587  perform_double_contraction(
1588  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1589  base_tensor_type &data,
1590  const typename SymmetricTensorAccessors::
1591  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1592  {
1593  using result_type = typename SymmetricTensorAccessors::
1594  double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1595 
1596  switch (dim)
1597  {
1598  case 1:
1599  return data[0] * sdata[0];
1600  default:
1601  // Start with the non-diagonal part to avoid some multiplications by
1602  // 2.
1603 
1604  result_type sum = data[dim] * sdata[dim];
1605  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1606  sum += data[d] * sdata[d];
1607  sum += sum; // sum = sum * 2.;
1608 
1609  // Now add the contributions from the diagonal
1610  for (unsigned int d = 0; d < dim; ++d)
1611  sum += data[d] * sdata[d];
1612  return sum;
1613  }
1614  }
1615 
1616 
1617 
1618  template <int dim, typename Number, typename OtherNumber = Number>
1619  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1620  typename SymmetricTensorAccessors::
1621  double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1622  perform_double_contraction(
1623  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1624  base_tensor_type &data,
1625  const typename SymmetricTensorAccessors::
1626  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1627  {
1628  using result_type = typename SymmetricTensorAccessors::
1629  double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1630  using value_type = typename SymmetricTensorAccessors::
1631  double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1632 
1633  const unsigned int data_dim = SymmetricTensorAccessors::
1634  StorageType<2, dim, value_type>::n_independent_components;
1635  value_type tmp[data_dim]{};
1636  for (unsigned int i = 0; i < data_dim; ++i)
1637  tmp[i] =
1638  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1639  return result_type(tmp);
1640  }
1641 
1642 
1643 
1644  template <int dim, typename Number, typename OtherNumber = Number>
1645  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1646  typename SymmetricTensorAccessors::StorageType<
1647  2,
1648  dim,
1649  typename SymmetricTensorAccessors::
1650  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1651  base_tensor_type
1652  perform_double_contraction(
1653  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1654  base_tensor_type &data,
1655  const typename SymmetricTensorAccessors::
1656  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1657  {
1658  using value_type = typename SymmetricTensorAccessors::
1659  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1660  using base_tensor_type = typename SymmetricTensorAccessors::
1661  StorageType<2, dim, value_type>::base_tensor_type;
1662 
1663  base_tensor_type tmp;
1664  for (unsigned int i = 0; i < tmp.dimension; ++i)
1665  {
1666  // Start with the non-diagonal part
1667  value_type sum = data[dim] * sdata[dim][i];
1668  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1669  sum += data[d] * sdata[d][i];
1670  sum += sum; // sum = sum * 2.;
1671 
1672  // Now add the contributions from the diagonal
1673  for (unsigned int d = 0; d < dim; ++d)
1674  sum += data[d] * sdata[d][i];
1675  tmp[i] = sum;
1676  }
1677  return tmp;
1678  }
1679 
1680 
1681 
1682  template <int dim, typename Number, typename OtherNumber = Number>
1683  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1684  typename SymmetricTensorAccessors::StorageType<
1685  4,
1686  dim,
1687  typename SymmetricTensorAccessors::
1688  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1689  base_tensor_type
1690  perform_double_contraction(
1691  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1692  base_tensor_type &data,
1693  const typename SymmetricTensorAccessors::
1694  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1695  {
1696  using value_type = typename SymmetricTensorAccessors::
1697  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1698  using base_tensor_type = typename SymmetricTensorAccessors::
1699  StorageType<4, dim, value_type>::base_tensor_type;
1700 
1701  const unsigned int data_dim = SymmetricTensorAccessors::
1702  StorageType<2, dim, value_type>::n_independent_components;
1703  base_tensor_type tmp;
1704  for (unsigned int i = 0; i < data_dim; ++i)
1705  for (unsigned int j = 0; j < data_dim; ++j)
1706  {
1707  // Start with the non-diagonal part
1708  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1709  tmp[i][j] += data[i][d] * sdata[d][j];
1710  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1711 
1712  // Now add the contributions from the diagonal
1713  for (unsigned int d = 0; d < dim; ++d)
1714  tmp[i][j] += data[i][d] * sdata[d][j];
1715  }
1716  return tmp;
1717  }
1718 
1719 } // end of namespace internal
1720 
1721 
1722 
1723 template <int rank_, int dim, typename Number>
1724 template <typename OtherNumber>
1725 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1726  typename internal::SymmetricTensorAccessors::
1727  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1730 {
1731  // need to have two different function calls
1732  // because a scalar and rank-2 tensor are not
1733  // the same data type (see internal function
1734  // above)
1735  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1736  s.data);
1737 }
1738 
1739 
1740 
1741 template <int rank_, int dim, typename Number>
1742 template <typename OtherNumber>
1743 DEAL_II_CONSTEXPR inline typename internal::SymmetricTensorAccessors::
1744  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1747 {
1748  typename internal::SymmetricTensorAccessors::
1749  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1750  tmp.data =
1751  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1752  s.data);
1753  return tmp;
1754 }
1755 
1756 
1757 
1758 // internal namespace to switch between the
1759 // access of different tensors. There used to
1760 // be explicit instantiations before for
1761 // different ranks and dimensions, but since
1762 // we now allow for templates on the data
1763 // type, and since we cannot partially
1764 // specialize the implementation, this got
1765 // into a separate namespace
1766 namespace internal
1767 {
1768  // The variables within this struct will be referenced in the next functions.
1769  // It is a workaround that allows returning a reference to a static variable
1770  // while allowing constexpr evaluation of the function.
1771  // It has to be defined outside the function because constexpr functions
1772  // cannot define static variables.
1773  // A similar struct has also been defined in tensor.h
1774  template <typename Type>
1775  struct Uninitialized
1776  {
1777  static Type value;
1778  };
1779 
1780  template <typename Type>
1781  Type Uninitialized<Type>::value;
1782 
1783  template <int dim, typename Number>
1784  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1785  symmetric_tensor_access(const TableIndices<2> &indices,
1786  typename SymmetricTensorAccessors::
1787  StorageType<2, dim, Number>::base_tensor_type &data)
1788  {
1789  // 1d is very simple and done first
1790  if (dim == 1)
1791  return data[0];
1792 
1793  // first treat the main diagonal elements, which are stored consecutively
1794  // at the beginning
1795  if (indices[0] == indices[1])
1796  return data[indices[0]];
1797 
1798  // the rest is messier and requires a few switches.
1799  switch (dim)
1800  {
1801  case 2:
1802  // at least for the 2x2 case it is reasonably simple
1803  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1804  ((indices[0] == 0) && (indices[1] == 1)),
1805  ExcInternalError());
1806  return data[2];
1807 
1808  default:
1809  // to do the rest, sort our indices before comparing
1810  {
1811  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1812  std::max(indices[0], indices[1]));
1813  for (unsigned int d = 0, c = 0; d < dim; ++d)
1814  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1815  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1816  return data[dim + c];
1817  Assert(false, ExcInternalError());
1818  }
1819  }
1820 
1821  // The code should never reach there.
1822  // Returns a dummy reference to a dummy variable just to make the
1823  // compiler happy.
1824  return Uninitialized<Number>::value;
1825  }
1826 
1827 
1828 
1829  template <int dim, typename Number>
1830  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1831  symmetric_tensor_access(const TableIndices<2> &indices,
1832  const typename SymmetricTensorAccessors::
1833  StorageType<2, dim, Number>::base_tensor_type &data)
1834  {
1835  // 1d is very simple and done first
1836  if (dim == 1)
1837  return data[0];
1838 
1839  // first treat the main diagonal elements, which are stored consecutively
1840  // at the beginning
1841  if (indices[0] == indices[1])
1842  return data[indices[0]];
1843 
1844  // the rest is messier and requires a few switches.
1845  switch (dim)
1846  {
1847  case 2:
1848  // at least for the 2x2 case it is reasonably simple
1849  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1850  ((indices[0] == 0) && (indices[1] == 1)),
1851  ExcInternalError());
1852  return data[2];
1853 
1854  default:
1855  // to do the rest, sort our indices before comparing
1856  {
1857  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1858  std::max(indices[0], indices[1]));
1859  for (unsigned int d = 0, c = 0; d < dim; ++d)
1860  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1861  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1862  return data[dim + c];
1863  Assert(false, ExcInternalError());
1864  }
1865  }
1866 
1867  // The code should never reach there.
1868  // Returns a dummy reference to a dummy variable just to make the
1869  // compiler happy.
1870  return Uninitialized<Number>::value;
1871  }
1872 
1873 
1874 
1875  template <int dim, typename Number>
1876  DEAL_II_CONSTEXPR inline Number &
1877  symmetric_tensor_access(const TableIndices<4> &indices,
1878  typename SymmetricTensorAccessors::
1879  StorageType<4, dim, Number>::base_tensor_type &data)
1880  {
1881  switch (dim)
1882  {
1883  case 1:
1884  return data[0][0];
1885 
1886  case 2:
1887  // each entry of the tensor can be thought of as an entry in a
1888  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1889  // rank-2 tensors. this is the format in which we store rank-4
1890  // tensors. determine which position the present entry is
1891  // stored in
1892  {
1893  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1894  return data[base_index[indices[0]][indices[1]]]
1895  [base_index[indices[2]][indices[3]]];
1896  }
1897  case 3:
1898  // each entry of the tensor can be thought of as an entry in a
1899  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1900  // rank-2 tensors. this is the format in which we store rank-4
1901  // tensors. determine which position the present entry is
1902  // stored in
1903  {
1904  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1905  {3, 1, 5},
1906  {4, 5, 2}};
1907  return data[base_index[indices[0]][indices[1]]]
1908  [base_index[indices[2]][indices[3]]];
1909  }
1910 
1911  default:
1912  Assert(false, ExcNotImplemented());
1913  }
1914 
1915  // The code should never reach there.
1916  // Returns a dummy reference to a dummy variable just to make the
1917  // compiler happy.
1918  return Uninitialized<Number>::value;
1919  }
1920 
1921 
1922  template <int dim, typename Number>
1923  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1924  symmetric_tensor_access(const TableIndices<4> &indices,
1925  const typename SymmetricTensorAccessors::
1926  StorageType<4, dim, Number>::base_tensor_type &data)
1927  {
1928  switch (dim)
1929  {
1930  case 1:
1931  return data[0][0];
1932 
1933  case 2:
1934  // each entry of the tensor can be thought of as an entry in a
1935  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1936  // rank-2 tensors. this is the format in which we store rank-4
1937  // tensors. determine which position the present entry is
1938  // stored in
1939  {
1940  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1941  return data[base_index[indices[0]][indices[1]]]
1942  [base_index[indices[2]][indices[3]]];
1943  }
1944  case 3:
1945  // each entry of the tensor can be thought of as an entry in a
1946  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1947  // rank-2 tensors. this is the format in which we store rank-4
1948  // tensors. determine which position the present entry is
1949  // stored in
1950  {
1951  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1952  {3, 1, 5},
1953  {4, 5, 2}};
1954  return data[base_index[indices[0]][indices[1]]]
1955  [base_index[indices[2]][indices[3]]];
1956  }
1957 
1958  default:
1959  Assert(false, ExcNotImplemented());
1960  }
1961 
1962  // The code should never reach there.
1963  // Returns a dummy reference to a dummy variable just to make the
1964  // compiler happy.
1965  return Uninitialized<Number>::value;
1966  }
1967 
1968 } // end of namespace internal
1969 
1970 
1971 
1972 template <int rank_, int dim, typename Number>
1973 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1975  operator()(const TableIndices<rank_> &indices)
1976 {
1977  for (unsigned int r = 0; r < rank; ++r)
1978  AssertIndexRange(indices[r], dimension);
1979  return internal::symmetric_tensor_access<dim, Number>(indices, data);
1980 }
1981 
1982 
1983 
1984 template <int rank_, int dim, typename Number>
1985 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1987  operator()(const TableIndices<rank_> &indices) const
1988 {
1989  for (unsigned int r = 0; r < rank; ++r)
1990  AssertIndexRange(indices[r], dimension);
1991  return internal::symmetric_tensor_access<dim, Number>(indices, data);
1992 }
1993 
1994 
1995 
1996 namespace internal
1997 {
1998  namespace SymmetricTensorImplementation
1999  {
2000  template <int rank_>
2001  constexpr TableIndices<rank_>
2002  get_partially_filled_indices(const unsigned int row,
2003  const std::integral_constant<int, 2> &)
2004  {
2006  }
2007 
2008 
2009  template <int rank_>
2010  constexpr TableIndices<rank_>
2011  get_partially_filled_indices(const unsigned int row,
2012  const std::integral_constant<int, 4> &)
2013  {
2014  return TableIndices<rank_>(row,
2018  }
2019  } // namespace SymmetricTensorImplementation
2020 } // namespace internal
2021 
2022 
2023 template <int rank_, int dim, typename Number>
2024 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2025  Accessor<rank_, dim, true, rank_ - 1, Number>
2027  operator[](const unsigned int row) const
2028 {
2029  return internal::SymmetricTensorAccessors::
2030  Accessor<rank_, dim, true, rank_ - 1, Number>(
2031  *this,
2032  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2033  rank_>(row, std::integral_constant<int, rank_>()));
2034 }
2035 
2036 
2037 
2038 template <int rank_, int dim, typename Number>
2039 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE internal::
2040  SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2041  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2042 {
2043  return internal::SymmetricTensorAccessors::
2044  Accessor<rank_, dim, false, rank_ - 1, Number>(
2045  *this,
2046  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2047  rank_>(row, std::integral_constant<int, rank_>()));
2048 }
2049 
2050 
2051 
2052 template <int rank_, int dim, typename Number>
2053 constexpr DEAL_II_ALWAYS_INLINE const Number &
2055  operator[](const TableIndices<rank_> &indices) const
2056 {
2057  return operator()(indices);
2058 }
2059 
2060 
2061 
2062 template <int rank_, int dim, typename Number>
2063 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
2065  operator[](const TableIndices<rank_> &indices)
2066 {
2067  return operator()(indices);
2068 }
2069 
2070 
2071 
2072 template <int rank_, int dim, typename Number>
2073 inline Number *
2075 {
2076  return std::addressof(this->access_raw_entry(0));
2077 }
2078 
2079 
2080 
2081 template <int rank_, int dim, typename Number>
2082 inline const Number *
2084 {
2085  return std::addressof(this->access_raw_entry(0));
2086 }
2087 
2088 
2089 
2090 template <int rank_, int dim, typename Number>
2091 inline Number *
2093 {
2095 }
2096 
2097 
2098 
2099 template <int rank_, int dim, typename Number>
2100 inline const Number *
2102 {
2104 }
2105 
2106 
2107 
2108 namespace internal
2109 {
2110  namespace SymmetricTensorImplementation
2111  {
2112  template <int dim, typename Number>
2113  constexpr unsigned int
2114  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2115  const unsigned int index)
2116  {
2117  return index;
2118  }
2119 
2120 
2121  template <int dim, typename Number>
2122  constexpr ::TableIndices<2>
2123  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2124  const unsigned int index)
2125  {
2128  }
2129 
2130  } // namespace SymmetricTensorImplementation
2131 } // namespace internal
2132 
2133 
2134 
2135 template <int rank_, int dim, typename Number>
2136 DEAL_II_CONSTEXPR inline const Number &
2138  const unsigned int index) const
2139 {
2140  AssertIndexRange(index, n_independent_components);
2141  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2142  index)];
2143 }
2144 
2145 
2146 
2147 template <int rank_, int dim, typename Number>
2148 DEAL_II_CONSTEXPR inline Number &
2150 {
2151  AssertIndexRange(index, n_independent_components);
2152  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2153  index)];
2154 }
2155 
2156 
2157 
2158 namespace internal
2159 {
2160  template <int dim, typename Number>
2161  DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
2162  compute_norm(const typename SymmetricTensorAccessors::
2163  StorageType<2, dim, Number>::base_tensor_type &data)
2164  {
2165  switch (dim)
2166  {
2167  case 1:
2168  return numbers::NumberTraits<Number>::abs(data[0]);
2169 
2170  case 2:
2171  return std::sqrt(
2175 
2176  case 3:
2177  return std::sqrt(
2184 
2185  default:
2186  {
2187  typename numbers::NumberTraits<Number>::real_type return_value =
2189 
2190  for (unsigned int d = 0; d < dim; ++d)
2191  return_value +=
2193  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2194  return_value +=
2196 
2197  return std::sqrt(return_value);
2198  }
2199  }
2200  }
2201 
2202 
2203 
2204  template <int dim, typename Number>
2205  DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
2206  compute_norm(const typename SymmetricTensorAccessors::
2207  StorageType<4, dim, Number>::base_tensor_type &data)
2208  {
2209  switch (dim)
2210  {
2211  case 1:
2212  return numbers::NumberTraits<Number>::abs(data[0][0]);
2213 
2214  default:
2215  {
2216  typename numbers::NumberTraits<Number>::real_type return_value =
2218 
2219  const unsigned int n_independent_components = data.dimension;
2220 
2221  for (unsigned int i = 0; i < dim; ++i)
2222  for (unsigned int j = 0; j < dim; ++j)
2223  return_value +=
2225  for (unsigned int i = 0; i < dim; ++i)
2226  for (unsigned int j = dim; j < n_independent_components; ++j)
2227  return_value +=
2229  for (unsigned int i = dim; i < n_independent_components; ++i)
2230  for (unsigned int j = 0; j < dim; ++j)
2231  return_value +=
2233  for (unsigned int i = dim; i < n_independent_components; ++i)
2234  for (unsigned int j = dim; j < n_independent_components; ++j)
2235  return_value +=
2237 
2238  return std::sqrt(return_value);
2239  }
2240  }
2241  }
2242 
2243 } // end of namespace internal
2244 
2245 
2246 
2247 template <int rank_, int dim, typename Number>
2250 {
2251  return internal::compute_norm<dim, Number>(data);
2252 }
2253 
2254 
2255 
2256 namespace internal
2257 {
2258  namespace SymmetricTensorImplementation
2259  {
2260  // a function to do the unrolling from a set of indices to a
2261  // scalar index into the array in which we store the elements of
2262  // a symmetric tensor
2263  //
2264  // this function is for rank-2 tensors
2265  template <int dim>
2266  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE unsigned int
2268  {
2269  AssertIndexRange(indices[0], dim);
2270  AssertIndexRange(indices[1], dim);
2271 
2272  switch (dim)
2273  {
2274  case 1:
2275  {
2276  return 0;
2277  }
2278 
2279  case 2:
2280  {
2281  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2282  return table[indices[0]][indices[1]];
2283  }
2284 
2285  case 3:
2286  {
2287  constexpr unsigned int table[3][3] = {{0, 3, 4},
2288  {3, 1, 5},
2289  {4, 5, 2}};
2290  return table[indices[0]][indices[1]];
2291  }
2292 
2293  case 4:
2294  {
2295  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2296  {4, 1, 7, 8},
2297  {5, 7, 2, 9},
2298  {6, 8, 9, 3}};
2299  return table[indices[0]][indices[1]];
2300  }
2301 
2302  default:
2303  // for the remainder, manually figure out the numbering
2304  {
2305  if (indices[0] == indices[1])
2306  return indices[0];
2307 
2308  TableIndices<2> sorted_indices(indices);
2309  sorted_indices.sort();
2310 
2311  for (unsigned int d = 0, c = 0; d < dim; ++d)
2312  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2313  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2314  return dim + c;
2315 
2316  // should never get here:
2317  Assert(false, ExcInternalError());
2318  return 0;
2319  }
2320  }
2321  }
2322 
2323  // a function to do the unrolling from a set of indices to a
2324  // scalar index into the array in which we store the elements of
2325  // a symmetric tensor
2326  //
2327  // this function is for tensors of ranks not already handled
2328  // above
2329  template <int dim, int rank_>
2330  DEAL_II_CONSTEXPR inline unsigned int
2332  {
2333  (void)indices;
2334  Assert(false, ExcNotImplemented());
2336  }
2337  } // namespace SymmetricTensorImplementation
2338 } // namespace internal
2339 
2340 
2341 template <int rank_, int dim, typename Number>
2342 constexpr unsigned int
2344  const TableIndices<rank_> &indices)
2345 {
2346  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2347  dim>(indices);
2348 }
2349 
2350 
2351 
2352 namespace internal
2353 {
2354  namespace SymmetricTensorImplementation
2355  {
2356  // a function to do the inverse of the unrolling from a set of
2357  // indices to a scalar index into the array in which we store
2358  // the elements of a symmetric tensor. in other words, it goes
2359  // from the scalar index into the array to a set of indices of
2360  // the tensor
2361  //
2362  // this function is for rank-2 tensors
2363  template <int dim>
2364  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2365  unrolled_to_component_indices(const unsigned int i,
2366  const std::integral_constant<int, 2> &)
2367  {
2368  Assert(
2370  ExcIndexRange(
2371  i,
2372  0,
2374  switch (dim)
2375  {
2376  case 1:
2377  {
2378  return {0, 0};
2379  }
2380 
2381  case 2:
2382  {
2383  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2384  TableIndices<2>(1, 1),
2385  TableIndices<2>(0, 1)};
2386  return table[i];
2387  }
2388 
2389  case 3:
2390  {
2391  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2392  TableIndices<2>(1, 1),
2393  TableIndices<2>(2, 2),
2394  TableIndices<2>(0, 1),
2395  TableIndices<2>(0, 2),
2396  TableIndices<2>(1, 2)};
2397  return table[i];
2398  }
2399 
2400  default:
2401  if (i < dim)
2402  return {i, i};
2403 
2404  for (unsigned int d = 0, c = 0; d < dim; ++d)
2405  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2406  if (c == i)
2407  return {d, e};
2408 
2409  // should never get here:
2410  Assert(false, ExcInternalError());
2411  return {0, 0};
2412  }
2413  }
2414 
2415  // a function to do the inverse of the unrolling from a set of
2416  // indices to a scalar index into the array in which we store
2417  // the elements of a symmetric tensor. in other words, it goes
2418  // from the scalar index into the array to a set of indices of
2419  // the tensor
2420  //
2421  // this function is for tensors of a rank not already handled
2422  // above
2423  template <int dim, int rank_>
2424  DEAL_II_CONSTEXPR inline
2425  typename std::enable_if<rank_ != 2, TableIndices<rank_>>::type
2426  unrolled_to_component_indices(const unsigned int i,
2427  const std::integral_constant<int, rank_> &)
2428  {
2429  (void)i;
2430  Assert(
2431  (i <
2433  ExcIndexRange(i,
2434  0,
2436  n_independent_components));
2437  Assert(false, ExcNotImplemented());
2438  return TableIndices<rank_>();
2439  }
2440 
2441  } // namespace SymmetricTensorImplementation
2442 } // namespace internal
2443 
2444 template <int rank_, int dim, typename Number>
2445 constexpr DEAL_II_ALWAYS_INLINE TableIndices<rank_>
2447  const unsigned int i)
2448 {
2449  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2450  dim>(i, std::integral_constant<int, rank_>());
2451 }
2452 
2453 
2454 
2455 template <int rank_, int dim, typename Number>
2456 template <class Archive>
2457 inline void
2458 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2459 {
2460  ar &data;
2461 }
2462 
2463 
2464 #endif // DOXYGEN
2465 
2466 /* ----------------- Non-member functions operating on tensors. ------------ */
2467 
2468 
2481 template <int rank_, int dim, typename Number, typename OtherNumber>
2482 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2486 {
2488  tmp = left;
2489  tmp += right;
2490  return tmp;
2491 }
2492 
2493 
2506 template <int rank_, int dim, typename Number, typename OtherNumber>
2507 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2511 {
2513  tmp = left;
2514  tmp -= right;
2515  return tmp;
2516 }
2517 
2518 
2526 template <int rank_, int dim, typename Number, typename OtherNumber>
2527 constexpr DEAL_II_ALWAYS_INLINE
2530  const Tensor<rank_, dim, OtherNumber> & right)
2531 {
2532  return Tensor<rank_, dim, Number>(left) + right;
2533 }
2534 
2535 
2543 template <int rank_, int dim, typename Number, typename OtherNumber>
2544 constexpr DEAL_II_ALWAYS_INLINE
2548 {
2549  return left + Tensor<rank_, dim, OtherNumber>(right);
2550 }
2551 
2552 
2560 template <int rank_, int dim, typename Number, typename OtherNumber>
2561 constexpr DEAL_II_ALWAYS_INLINE
2564  const Tensor<rank_, dim, OtherNumber> & right)
2565 {
2566  return Tensor<rank_, dim, Number>(left) - right;
2567 }
2568 
2569 
2577 template <int rank_, int dim, typename Number, typename OtherNumber>
2578 constexpr DEAL_II_ALWAYS_INLINE
2582 {
2583  return left - Tensor<rank_, dim, OtherNumber>(right);
2584 }
2585 
2586 
2587 
2601 template <int dim, typename Number>
2602 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2604 {
2605  switch (dim)
2606  {
2607  case 1:
2608  return t.data[0];
2609  case 2:
2610  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2611  case 3:
2612  {
2613  // in analogy to general tensors, but
2614  // there's something to be simplified for
2615  // the present case
2616  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2617  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2618  t.data[0] * t.data[5] * t.data[5] -
2619  t.data[1] * t.data[4] * t.data[4] -
2620  t.data[2] * t.data[3] * t.data[3]);
2621  }
2622  default:
2623  Assert(false, ExcNotImplemented());
2624  return internal::NumberType<Number>::value(0.0);
2625  }
2626 }
2627 
2628 
2629 
2639 template <int dim, typename Number>
2640 constexpr DEAL_II_ALWAYS_INLINE Number
2642 {
2643  return determinant(t);
2644 }
2645 
2646 
2647 
2655 template <int dim, typename Number>
2656 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2658 {
2659  Number t = d.data[0];
2660  for (unsigned int i = 1; i < dim; ++i)
2661  t += d.data[i];
2662  return t;
2663 }
2664 
2665 
2675 template <int dim, typename Number>
2676 constexpr Number
2678 {
2679  return trace(t);
2680 }
2681 
2682 
2695 template <typename Number>
2696 constexpr DEAL_II_ALWAYS_INLINE Number
2698 {
2699  return internal::NumberType<Number>::value(0.0);
2700 }
2701 
2702 
2703 
2723 template <typename Number>
2724 constexpr DEAL_II_ALWAYS_INLINE Number
2726 {
2727  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2728 }
2729 
2730 
2731 
2741 template <typename Number>
2742 constexpr DEAL_II_ALWAYS_INLINE Number
2744 {
2745  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2746  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2747 }
2748 
2749 
2750 
2759 template <typename Number>
2760 std::array<Number, 1>
2762 
2763 
2764 
2788 template <typename Number>
2789 std::array<Number, 2>
2791 
2792 
2793 
2816 template <typename Number>
2817 std::array<Number, 3>
2819 
2820 
2821 
2822 namespace internal
2823 {
2824  namespace SymmetricTensorImplementation
2825  {
2865  template <int dim, typename Number>
2866  void
2867  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2868  ::Tensor<2, dim, Number> & Q,
2869  std::array<Number, dim> & d,
2870  std::array<Number, dim - 1> & e);
2871 
2872 
2873 
2915  template <int dim, typename Number>
2916  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2917  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2918 
2919 
2920 
2962  template <int dim, typename Number>
2963  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2965 
2966 
2967 
2983  template <typename Number>
2984  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
2985  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
2986 
2987 
2988 
3023  template <typename Number>
3024  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3025  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3026 
3031  template <int dim, typename Number>
3033  {
3034  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3035  bool
3036  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3037  {
3038  return lhs.first > rhs.first;
3039  }
3040  };
3041 
3042  } // namespace SymmetricTensorImplementation
3043 
3044 } // namespace internal
3045 
3046 
3047 
3048 // The line below is to ensure that doxygen puts the full description
3049 // of this global enumeration into the documentation
3050 // See https://stackoverflow.com/a/1717984
3080 {
3090  hybrid,
3108  jacobi
3109 };
3110 
3111 
3112 
3142 template <int dim, typename Number>
3143 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3144  std::integral_constant<int, dim>::value>
3146  const SymmetricTensorEigenvectorMethod method =
3148 
3149 
3150 
3160 template <int rank_, int dim, typename Number>
3161 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
3163 {
3164  return t;
3165 }
3166 
3167 
3168 
3178 template <int dim, typename Number>
3179 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3181 {
3183 
3184  // subtract scaled trace from the diagonal
3185  const Number tr = trace(t) / dim;
3186  for (unsigned int i = 0; i < dim; ++i)
3187  tmp.data[i] -= tr;
3188 
3189  return tmp;
3190 }
3191 
3192 
3193 
3201 template <int dim, typename Number>
3202 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3204 {
3205  // create a default constructed matrix filled with
3206  // zeros, then set the diagonal elements to one
3208  switch (dim)
3209  {
3210  case 1:
3211  tmp.data[0] = internal::NumberType<Number>::value(1.);
3212  break;
3213  case 2:
3214  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3215  break;
3216  case 3:
3217  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3218  internal::NumberType<Number>::value(1.);
3219  break;
3220  default:
3221  for (unsigned int d = 0; d < dim; ++d)
3222  tmp.data[d] = internal::NumberType<Number>::value(1.);
3223  }
3224  return tmp;
3225 }
3226 
3227 
3228 
3237 template <int dim>
3238 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim>
3240 {
3241  return unit_symmetric_tensor<dim, double>();
3242 }
3243 
3244 
3245 
3260 template <int dim, typename Number>
3261 DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
3263 {
3265 
3266  // fill the elements treating the diagonal
3267  for (unsigned int i = 0; i < dim; ++i)
3268  for (unsigned int j = 0; j < dim; ++j)
3269  tmp.data[i][j] =
3270  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3271 
3272  // then fill the ones that copy over the
3273  // non-diagonal elements. note that during
3274  // the double-contraction, we handle the
3275  // off-diagonal elements twice, so simply
3276  // copying requires a weight of 1/2
3277  for (unsigned int i = dim;
3278  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3279  n_rank2_components;
3280  ++i)
3281  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3282 
3283  return tmp;
3284 }
3285 
3286 
3287 
3302 template <int dim>
3303 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
3305 {
3306  return deviator_tensor<dim, double>();
3307 }
3308 
3309 
3310 
3333 template <int dim, typename Number>
3334 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
3336 {
3338 
3339  // fill the elements treating the diagonal
3340  for (unsigned int i = 0; i < dim; ++i)
3341  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3342 
3343  // then fill the ones that copy over the
3344  // non-diagonal elements. note that during
3345  // the double-contraction, we handle the
3346  // off-diagonal elements twice, so simply
3347  // copying requires a weight of 1/2
3348  for (unsigned int i = dim;
3349  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3350  n_rank2_components;
3351  ++i)
3352  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3353 
3354  return tmp;
3355 }
3356 
3357 
3358 
3380 template <int dim>
3381 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
3383 {
3384  return identity_tensor<dim, double>();
3385 }
3386 
3387 
3388 
3399 template <int dim, typename Number>
3400 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3402 {
3404  value(t);
3405 }
3406 
3407 
3408 
3420 template <int dim, typename Number>
3423 {
3425  value(t);
3426 }
3427 
3428 
3429 
3444 template <int dim, typename Number>
3445 DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
3448 {
3450 
3451  // fill only the elements really needed
3452  for (unsigned int i = 0; i < dim; ++i)
3453  for (unsigned int j = i; j < dim; ++j)
3454  for (unsigned int k = 0; k < dim; ++k)
3455  for (unsigned int l = k; l < dim; ++l)
3456  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3457 
3458  return tmp;
3459 }
3460 
3461 
3462 
3471 template <int dim, typename Number>
3472 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3474 {
3476  for (unsigned int d = 0; d < dim; ++d)
3477  result[d][d] = t[d][d];
3478  Number half = 0.5;
3479  for (unsigned int d = 0; d < dim; ++d)
3480  for (unsigned int e = d + 1; e < dim; ++e)
3481  result[d][e] = (t[d][e] + t[e][d]) * half;
3482  return result;
3483 }
3484 
3485 
3486 
3494 template <int rank_, int dim, typename Number>
3495 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
3497  operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3498 {
3500  tt *= factor;
3501  return tt;
3502 }
3503 
3504 
3505 
3513 template <int rank_, int dim, typename Number>
3514 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
3515  operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3516 {
3517  // simply forward to the other operator
3518  return t * factor;
3519 }
3520 
3521 
3522 
3548 template <int rank_, int dim, typename Number, typename OtherNumber>
3549 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3550  rank_,
3551  dim,
3552  typename ProductType<Number,
3553  typename EnableIfScalar<OtherNumber>::type>::type>
3555  const OtherNumber & factor)
3556 {
3557  // form the product. we have to convert the two factors into the final
3558  // type via explicit casts because, for awkward reasons, the C++
3559  // standard committee saw it fit to not define an
3560  // operator*(float,std::complex<double>)
3561  // (as well as with switched arguments and double<->float).
3562  using product_type = typename ProductType<Number, OtherNumber>::type;
3564  tt *= internal::NumberType<product_type>::value(factor);
3565  return tt;
3566 }
3567 
3568 
3569 
3578 template <int rank_, int dim, typename Number, typename OtherNumber>
3579 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3580  rank_,
3581  dim,
3582  typename ProductType<OtherNumber,
3583  typename EnableIfScalar<Number>::type>::type>
3584 operator*(const Number & factor,
3586 {
3587  // simply forward to the other operator with switched arguments
3588  return (t * factor);
3589 }
3590 
3591 
3592 
3598 template <int rank_, int dim, typename Number, typename OtherNumber>
3599 DEAL_II_CONSTEXPR inline SymmetricTensor<
3600  rank_,
3601  dim,
3602  typename ProductType<Number,
3603  typename EnableIfScalar<OtherNumber>::type>::type>
3605  const OtherNumber & factor)
3606 {
3607  using product_type = typename ProductType<Number, OtherNumber>::type;
3609  tt /= internal::NumberType<product_type>::value(factor);
3610  return tt;
3611 }
3612 
3613 
3614 
3621 template <int rank_, int dim>
3622 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim>
3623  operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3624 {
3626  tt *= factor;
3627  return tt;
3628 }
3629 
3630 
3631 
3638 template <int rank_, int dim>
3639 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim>
3640  operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3641 {
3643  tt *= factor;
3644  return tt;
3645 }
3646 
3647 
3648 
3654 template <int rank_, int dim>
3655 DEAL_II_CONSTEXPR inline SymmetricTensor<rank_, dim>
3656 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3657 {
3659  tt /= factor;
3660  return tt;
3661 }
3662 
3672 template <int dim, typename Number, typename OtherNumber>
3673 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
3676 {
3677  return (t1 * t2);
3678 }
3679 
3680 
3690 template <int dim, typename Number, typename OtherNumber>
3691 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
3692  typename ProductType<Number, OtherNumber>::type
3694  const Tensor<2, dim, OtherNumber> & t2)
3695 {
3696  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<
3697  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3698  for (unsigned int i = 0; i < dim; ++i)
3699  for (unsigned int j = 0; j < dim; ++j)
3700  s += t1[i][j] * t2[i][j];
3701  return s;
3702 }
3703 
3704 
3714 template <int dim, typename Number, typename OtherNumber>
3715 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
3718 {
3719  return scalar_product(t2, t1);
3720 }
3721 
3722 
3738 template <typename Number, typename OtherNumber>
3739 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void double_contract(
3740  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3743 {
3744  tmp[0][0] = t[0][0][0][0] * s[0][0];
3745 }
3746 
3747 
3748 
3764 template <typename Number, typename OtherNumber>
3765 DEAL_II_CONSTEXPR inline void double_contract(
3766  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3769 {
3770  tmp[0][0] = t[0][0][0][0] * s[0][0];
3771 }
3772 
3773 
3774 
3789 template <typename Number, typename OtherNumber>
3790 DEAL_II_CONSTEXPR inline void double_contract(
3791  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3794 {
3795  const unsigned int dim = 2;
3796 
3797  for (unsigned int i = 0; i < dim; ++i)
3798  for (unsigned int j = i; j < dim; ++j)
3799  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3800  2 * t[i][j][0][1] * s[0][1];
3801 }
3802 
3803 
3804 
3820 template <typename Number, typename OtherNumber>
3821 DEAL_II_CONSTEXPR inline void double_contract(
3822  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3825 {
3826  const unsigned int dim = 2;
3827 
3828  for (unsigned int i = 0; i < dim; ++i)
3829  for (unsigned int j = i; j < dim; ++j)
3830  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3831  2 * s[0][1] * t[0][1][i][j];
3832 }
3833 
3834 
3835 
3851 template <typename Number, typename OtherNumber>
3852 DEAL_II_CONSTEXPR inline void double_contract(
3853  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3856 {
3857  const unsigned int dim = 3;
3858 
3859  for (unsigned int i = 0; i < dim; ++i)
3860  for (unsigned int j = i; j < dim; ++j)
3861  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3862  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3863  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3864 }
3865 
3866 
3867 
3883 template <typename Number, typename OtherNumber>
3884 DEAL_II_CONSTEXPR inline void double_contract(
3885  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
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] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3894  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3895  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3896 }
3897 
3898 
3899 
3907 template <int dim, typename Number, typename OtherNumber>
3908 DEAL_II_CONSTEXPR
3911  const Tensor<1, dim, OtherNumber> & src2)
3912 {
3914  for (unsigned int i = 0; i < dim; ++i)
3915  for (unsigned int j = 0; j < dim; ++j)
3916  dest[i] += src1[i][j] * src2[j];
3917  return dest;
3918 }
3919 
3920 
3928 template <int dim, typename Number, typename OtherNumber>
3932 {
3933  // this is easy for symmetric tensors:
3934  return src2 * src1;
3935 }
3936 
3937 
3938 
3959 template <int rank_1,
3960  int rank_2,
3961  int dim,
3962  typename Number,
3963  typename OtherNumber>
3964 constexpr DEAL_II_ALWAYS_INLINE
3965  typename Tensor<rank_1 + rank_2 - 2,
3966  dim,
3967  typename ProductType<Number, OtherNumber>::type>::tensor_type
3970 {
3971  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
3972 }
3973 
3974 
3975 
3996 template <int rank_1,
3997  int rank_2,
3998  int dim,
3999  typename Number,
4000  typename OtherNumber>
4001 constexpr DEAL_II_ALWAYS_INLINE
4002  typename Tensor<rank_1 + rank_2 - 2,
4003  dim,
4004  typename ProductType<Number, OtherNumber>::type>::tensor_type
4006  const Tensor<rank_2, dim, OtherNumber> & src2)
4007 {
4008  return Tensor<rank_2, dim, OtherNumber>(src1) * src2;
4009 }
4010 
4011 
4012 
4022 template <int dim, typename Number>
4023 inline std::ostream &
4024 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4025 {
4026  // make our lives a bit simpler by outputting
4027  // the tensor through the operator for the
4028  // general Tensor class
4030 
4031  for (unsigned int i = 0; i < dim; ++i)
4032  for (unsigned int j = 0; j < dim; ++j)
4033  tt[i][j] = t[i][j];
4034 
4035  return out << tt;
4036 }
4037 
4038 
4039 
4049 template <int dim, typename Number>
4050 inline std::ostream &
4051 operator<<(std::ostream &out, const SymmetricTensor<4, 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  for (unsigned int k = 0; k < dim; ++k)
4061  for (unsigned int l = 0; l < dim; ++l)
4062  tt[i][j][k][l] = t[i][j][k][l];
4063 
4064  return out << tt;
4065 }
4066 
4067 
4068 DEAL_II_NAMESPACE_CLOSE
4069 
4070 #endif
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static const unsigned int invalid_unsigned_int
Definition: types.h:190
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
static constexpr unsigned int n_independent_components
friend constexpr SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
constexpr SymmetricTensor operator-() const
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
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< 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 numbers::NumberTraits< Number >::real_type norm() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr bool operator==(const SymmetricTensor &) const
STL namespace.
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
friend constexpr SymmetricTensor< 4, dim2, Number2 > identity_tensor()
static real_type abs(const number &x)
Definition: numbers.h:625
SymmetricTensorEigenvectorMethod
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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)
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
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 unsigned int rank
constexpr Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Number * begin_raw()
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 first_invariant(const SymmetricTensor< 2, dim, Number > &t)
#define Assert(cond, exc)
Definition: exceptions.h:1419
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
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)
Number * end_raw()
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
static constexpr std::size_t memory_consumption()
constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: tensor.h:417
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
friend constexpr SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:957
static ::ExceptionBase & ExcNotImplemented()
constexpr Number & operator()(const TableIndices< rank_ > &indices)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
static ::ExceptionBase & ExcInternalError()
constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)