Reference documentation for deal.II version 9.2.0
\(\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>
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>
41  unit_symmetric_tensor();
42 
43 template <int dim, typename Number>
44 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
45  deviator_tensor();
46 
47 template <int dim, typename Number>
48 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
49  identity_tensor();
50 
51 template <int dim, typename Number>
52 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
53  invert(const SymmetricTensor<2, dim, Number> &);
54 
55 template <int dim, typename Number>
56 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
57  invert(const SymmetricTensor<4, dim, Number> &);
58 
59 template <int dim2, typename Number>
60 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
61  trace(const SymmetricTensor<2, dim2, Number> &);
62 
63 template <int dim, typename Number>
64 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
65  deviator(const SymmetricTensor<2, dim, Number> &);
66 
67 template <int dim, typename Number>
68 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
69  determinant(const SymmetricTensor<2, dim, 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  // - Jean-Paul Pelteret, Matthias Maier, Daniel Arndt 2020
80  template <int rank, int dim, typename T, typename U>
81  struct ProductTypeImpl<SymmetricTensor<rank, dim, T>, std::complex<U>>
82  {
83  using type =
84  SymmetricTensor<rank,
85  dim,
86  std::complex<typename ProductType<T, U>::type>>;
87  };
88 
89  template <int rank, int dim, typename T, typename U>
90  struct ProductTypeImpl<SymmetricTensor<rank, dim, std::complex<T>>,
91  std::complex<U>>
92  {
93  using type =
94  SymmetricTensor<rank,
95  dim,
96  std::complex<typename ProductType<T, U>::type>>;
97  };
98 
99  template <typename T, int rank, int dim, typename U>
100  struct ProductTypeImpl<std::complex<T>, SymmetricTensor<rank, dim, U>>
101  {
102  using type =
103  SymmetricTensor<rank,
104  dim,
105  std::complex<typename ProductType<T, U>::type>>;
106  };
107 
108  template <int rank, int dim, typename T, typename U>
109  struct ProductTypeImpl<std::complex<T>,
110  SymmetricTensor<rank, dim, std::complex<U>>>
111  {
112  using type =
113  SymmetricTensor<rank,
114  dim,
115  std::complex<typename ProductType<T, U>::type>>;
116  };
117  // end workaround
118 
123  namespace SymmetricTensorImplementation
124  {
129  template <int rank, int dim, typename Number>
130  struct Inverse;
131  } // namespace SymmetricTensorImplementation
132 
137  namespace SymmetricTensorAccessors
138  {
145  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
146  merge(const TableIndices<2> &previous_indices,
147  const unsigned int new_index,
148  const unsigned int position)
149  {
150  AssertIndexRange(position, 2);
151 
152  if (position == 0)
153  return {new_index, numbers::invalid_unsigned_int};
154  else
155  return {previous_indices[0], new_index};
156  }
157 
158 
159 
166  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<4>
167  merge(const TableIndices<4> &previous_indices,
168  const unsigned int new_index,
169  const unsigned int position)
170  {
171  AssertIndexRange(position, 4);
172 
173  switch (position)
174  {
175  case 0:
176  return {new_index,
177  numbers::invalid_unsigned_int,
178  numbers::invalid_unsigned_int,
179  numbers::invalid_unsigned_int};
180  case 1:
181  return {previous_indices[0],
182  new_index,
183  numbers::invalid_unsigned_int,
184  numbers::invalid_unsigned_int};
185  case 2:
186  return {previous_indices[0],
187  previous_indices[1],
188  new_index,
189  numbers::invalid_unsigned_int};
190  case 3:
191  return {previous_indices[0],
192  previous_indices[1],
193  previous_indices[2],
194  new_index};
195  default:
196  Assert(false, ExcInternalError());
197  return {};
198  }
199  }
200 
201 
210  template <int rank1,
211  int rank2,
212  int dim,
213  typename Number,
214  typename OtherNumber = Number>
215  struct double_contraction_result
216  {
217  using value_type = typename ProductType<Number, OtherNumber>::type;
218  using type =
219  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
220  };
221 
222 
231  template <int dim, typename Number, typename OtherNumber>
232  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
233  {
234  using type = typename ProductType<Number, OtherNumber>::type;
235  };
236 
237 
238 
251  template <int rank, int dim, typename Number>
252  struct StorageType;
253 
257  template <int dim, typename Number>
258  struct StorageType<2, dim, Number>
259  {
264  static const unsigned int n_independent_components =
265  (dim * dim + dim) / 2;
266 
270  using base_tensor_type = Tensor<1, n_independent_components, Number>;
271  };
272 
273 
274 
278  template <int dim, typename Number>
279  struct StorageType<4, dim, Number>
280  {
286  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
287 
291  static const unsigned int n_independent_components =
292  (n_rank2_components *
293  StorageType<2, dim, Number>::n_independent_components);
294 
301  using base_tensor_type = Tensor<2, n_rank2_components, Number>;
302  };
303 
304 
305 
310  template <int rank, int dim, bool constness, typename Number>
311  struct AccessorTypes;
312 
319  template <int rank, int dim, typename Number>
320  struct AccessorTypes<rank, dim, true, Number>
321  {
322  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
323 
324  using reference = Number;
325  };
326 
333  template <int rank, int dim, typename Number>
334  struct AccessorTypes<rank, dim, false, Number>
335  {
336  using tensor_type = ::SymmetricTensor<rank, dim, Number>;
337 
338  using reference = Number &;
339  };
340 
341 
376  template <int rank, int dim, bool constness, int P, typename Number>
377  class Accessor
378  {
379  public:
383  using reference =
384  typename AccessorTypes<rank, dim, constness, Number>::reference;
385  using tensor_type =
386  typename AccessorTypes<rank, dim, constness, Number>::tensor_type;
387 
388  private:
407  constexpr Accessor(tensor_type & tensor,
408  const TableIndices<rank> &previous_indices);
409 
413  constexpr DEAL_II_ALWAYS_INLINE
414  Accessor(const Accessor &) = default;
415 
416  public:
420  DEAL_II_CONSTEXPR Accessor<rank, dim, constness, P - 1, Number>
421  operator[](const unsigned int i);
422 
426  constexpr Accessor<rank, dim, constness, P - 1, Number>
427  operator[](const unsigned int i) const;
428 
429  private:
433  tensor_type & tensor;
434  const TableIndices<rank> previous_indices;
435 
436  // Declare some other classes as friends. Make sure to work around bugs
437  // in some compilers:
438  template <int, int, typename>
439  friend class ::SymmetricTensor;
440  template <int, int, bool, int, typename>
441  friend class Accessor;
442 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
443  friend class ::SymmetricTensor<rank, dim, Number>;
444  friend class Accessor<rank, dim, constness, P + 1, Number>;
445 #endif
446  };
447 
448 
449 
459  template <int rank, int dim, bool constness, typename Number>
460  class Accessor<rank, dim, constness, 1, Number>
461  {
462  public:
466  using reference =
467  typename AccessorTypes<rank, dim, constness, Number>::reference;
468  using tensor_type =
469  typename AccessorTypes<rank, dim, constness, Number>::tensor_type;
470 
471  private:
493  constexpr Accessor(tensor_type & tensor,
494  const TableIndices<rank> &previous_indices);
495 
499  constexpr DEAL_II_ALWAYS_INLINE
500  Accessor(const Accessor &) = default;
501 
502  public:
506  DEAL_II_CONSTEXPR reference operator[](const unsigned int);
507 
511  constexpr reference operator[](const unsigned int) const;
512 
513  private:
517  tensor_type & tensor;
518  const TableIndices<rank> previous_indices;
519 
520  // Declare some other classes as friends. Make sure to work around bugs
521  // in some compilers:
522  template <int, int, typename>
523  friend class ::SymmetricTensor;
524  template <int, int, bool, int, typename>
525  friend class SymmetricTensorAccessors::Accessor;
526 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
527  friend class ::SymmetricTensor<rank, dim, Number>;
528  friend class SymmetricTensorAccessors::
529  Accessor<rank, dim, constness, 2, Number>;
530 #endif
531  };
532  } // namespace SymmetricTensorAccessors
533 } // namespace internal
534 
535 
536 
610 template <int rank_, int dim, typename Number>
611 class SymmetricTensor
612 {
613 public:
614  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
615 
624  static const unsigned int dimension = dim;
625 
629  static const unsigned int rank = rank_;
630 
636  static constexpr unsigned int n_independent_components =
637  internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>::
638  n_independent_components;
639 
643  constexpr DEAL_II_ALWAYS_INLINE
644  SymmetricTensor() = default;
645 
659  template <typename OtherNumber>
660  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
661 
677  DEAL_II_CONSTEXPR
678  SymmetricTensor(const Number (&array)[n_independent_components]);
679 
685  template <typename OtherNumber>
686  constexpr explicit SymmetricTensor(
687  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
688 
692  Number *
693  begin_raw();
694 
698  const Number *
699  begin_raw() const;
700 
704  Number *
705  end_raw();
706 
711  const Number *
712  end_raw() const;
713 
720  template <typename OtherNumber>
721  DEAL_II_CONSTEXPR SymmetricTensor &
722  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
723 
730  DEAL_II_CONSTEXPR SymmetricTensor &
731  operator=(const Number &d);
732 
737  constexpr operator Tensor<rank_, dim, Number>() const;
738 
742  constexpr bool
743  operator==(const SymmetricTensor &) const;
744 
748  constexpr bool
749  operator!=(const SymmetricTensor &) const;
750 
754  template <typename OtherNumber>
755  DEAL_II_CONSTEXPR SymmetricTensor &
756  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
757 
761  template <typename OtherNumber>
762  DEAL_II_CONSTEXPR SymmetricTensor &
763  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
764 
769  template <typename OtherNumber>
770  DEAL_II_CONSTEXPR SymmetricTensor &
771  operator*=(const OtherNumber &factor);
772 
776  template <typename OtherNumber>
777  DEAL_II_CONSTEXPR SymmetricTensor &
778  operator/=(const OtherNumber &factor);
779 
783  DEAL_II_CONSTEXPR SymmetricTensor
784  operator-() const;
785 
812  template <typename OtherNumber>
813  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
814  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
815  operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const;
816 
821  template <typename OtherNumber>
822  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
823  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
824  operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const;
825 
829  DEAL_II_CONSTEXPR Number &
830  operator()(const TableIndices<rank_> &indices);
831 
835  DEAL_II_CONSTEXPR const Number &
836  operator()(const TableIndices<rank_> &indices) const;
837 
842  constexpr internal::SymmetricTensorAccessors::
843  Accessor<rank_, dim, true, rank_ - 1, Number>
844  operator[](const unsigned int row) const;
845 
850  DEAL_II_CONSTEXPR internal::SymmetricTensorAccessors::
851  Accessor<rank_, dim, false, rank_ - 1, Number>
852  operator[](const unsigned int row);
853 
859  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
860 
866  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
867 
874  DEAL_II_CONSTEXPR const Number &
875  access_raw_entry(const unsigned int unrolled_index) const;
876 
883  DEAL_II_CONSTEXPR Number &
884  access_raw_entry(const unsigned int unrolled_index);
885 
895  constexpr typename numbers::NumberTraits<Number>::real_type
896  norm() const;
897 
905  static constexpr unsigned int
906  component_to_unrolled_index(const TableIndices<rank_> &indices);
907 
913  static constexpr TableIndices<rank_>
914  unrolled_to_component_indices(const unsigned int i);
915 
928  DEAL_II_CONSTEXPR void
929  clear();
930 
935  static constexpr std::size_t
936  memory_consumption();
937 
942  template <class Archive>
943  void
944  serialize(Archive &ar, const unsigned int version);
945 
946 private:
950  using base_tensor_descriptor =
951  internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>;
952 
956  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
957 
961  base_tensor_type data;
962 
963  // Make all other symmetric tensors friends.
964  template <int, int, typename>
965  friend class SymmetricTensor;
966 
967  // Make a few more functions friends.
968  template <int dim2, typename Number2>
969  friend DEAL_II_CONSTEXPR Number2
970  trace(const SymmetricTensor<2, dim2, Number2> &d);
971 
972  template <int dim2, typename Number2>
973  friend DEAL_II_CONSTEXPR Number2
974  determinant(const SymmetricTensor<2, dim2, Number2> &t);
975 
976  template <int dim2, typename Number2>
977  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
978  deviator(const SymmetricTensor<2, dim2, Number2> &t);
979 
980  template <int dim2, typename Number2>
981  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
982  unit_symmetric_tensor();
983 
984  template <int dim2, typename Number2>
985  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
986  deviator_tensor();
987 
988  template <int dim2, typename Number2>
989  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
990  identity_tensor();
991 
992 
993  // Make a few helper classes friends as well.
994  friend struct internal::SymmetricTensorImplementation::
995  Inverse<2, dim, Number>;
996 
997  friend struct internal::SymmetricTensorImplementation::
998  Inverse<4, dim, Number>;
999 };
1000 
1001 
1002 
1003 // ------------------------- inline functions ------------------------
1004 
1005 #ifndef DOXYGEN
1006 
1007 // provide declarations for static members
1008 template <int rank, int dim, typename Number>
1009 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
1010 
1011 template <int rank_, int dim, typename Number>
1012 constexpr unsigned int
1013  SymmetricTensor<rank_, dim, Number>::n_independent_components;
1014 
1015 namespace internal
1016 {
1017  namespace SymmetricTensorAccessors
1018  {
1019  template <int rank_, int dim, bool constness, int P, typename Number>
1020  constexpr DEAL_II_ALWAYS_INLINE
1021  Accessor<rank_, dim, constness, P, Number>::Accessor(
1022  tensor_type & tensor,
1023  const TableIndices<rank_> &previous_indices)
1024  : tensor(tensor)
1025  , previous_indices(previous_indices)
1026  {}
1027 
1028 
1029 
1030  template <int rank_, int dim, bool constness, int P, typename Number>
1031  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1032  Accessor<rank_, dim, constness, P - 1, Number>
1033  Accessor<rank_, dim, constness, P, Number>::
1034  operator[](const unsigned int i)
1035  {
1036  return Accessor<rank_, dim, constness, P - 1, Number>(
1037  tensor, merge(previous_indices, i, rank_ - P));
1038  }
1039 
1040 
1041 
1042  template <int rank_, int dim, bool constness, int P, typename Number>
1043  constexpr DEAL_II_ALWAYS_INLINE
1044  Accessor<rank_, dim, constness, P - 1, Number>
1045  Accessor<rank_, dim, constness, P, Number>::
1046  operator[](const unsigned int i) const
1047  {
1048  return Accessor<rank_, dim, constness, P - 1, Number>(
1049  tensor, merge(previous_indices, i, rank_ - P));
1050  }
1051 
1052 
1053 
1054  template <int rank_, int dim, bool constness, typename Number>
1055  constexpr DEAL_II_ALWAYS_INLINE
1056  Accessor<rank_, dim, constness, 1, Number>::Accessor(
1057  tensor_type & tensor,
1058  const TableIndices<rank_> &previous_indices)
1059  : tensor(tensor)
1060  , previous_indices(previous_indices)
1061  {}
1062 
1063 
1064 
1065  template <int rank_, int dim, bool constness, typename Number>
1066  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1067  typename Accessor<rank_, dim, constness, 1, Number>::reference
1068  Accessor<rank_, dim, constness, 1, Number>::
1069  operator[](const unsigned int i)
1070  {
1071  return tensor(merge(previous_indices, i, rank_ - 1));
1072  }
1073 
1074 
1075  template <int rank_, int dim, bool constness, typename Number>
1076  constexpr DEAL_II_ALWAYS_INLINE
1077  typename Accessor<rank_, dim, constness, 1, Number>::reference
1078  Accessor<rank_, dim, constness, 1, Number>::
1079  operator[](const unsigned int i) const
1080  {
1081  return tensor(merge(previous_indices, i, rank_ - 1));
1082  }
1083  } // namespace SymmetricTensorAccessors
1084 } // namespace internal
1085 
1086 
1087 
1088 template <int rank_, int dim, typename Number>
1089 template <typename OtherNumber>
1090 inline DEAL_II_ALWAYS_INLINE
1091 SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1092  const Tensor<2, dim, OtherNumber> &t)
1093 {
1094  static_assert(rank == 2, "This function is only implemented for rank==2");
1095  for (unsigned int d = 0; d < dim; ++d)
1096  for (unsigned int e = 0; e < d; ++e)
1097  Assert(t[d][e] == t[e][d],
1098  ExcMessage("The incoming Tensor must be exactly symmetric."));
1099 
1100  for (unsigned int d = 0; d < dim; ++d)
1101  data[d] = t[d][d];
1102 
1103  for (unsigned int d = 0, c = 0; d < dim; ++d)
1104  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1105  data[dim + c] = t[d][e];
1106 }
1107 
1108 
1109 
1110 template <int rank_, int dim, typename Number>
1111 template <typename OtherNumber>
1112 constexpr DEAL_II_ALWAYS_INLINE
1113 SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1114  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1115  : data(initializer.data)
1116 {}
1117 
1118 
1119 
1120 template <int rank_, int dim, typename Number>
1121 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1122 SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1123  const Number (&array)[n_independent_components])
1124  : data(
1125  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1126 {
1127  // ensure that the reinterpret_cast above actually works
1128  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1129  ExcInternalError());
1130 }
1131 
1132 
1133 
1134 template <int rank_, int dim, typename Number>
1135 template <typename OtherNumber>
1136 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1137  SymmetricTensor<rank_, dim, Number> &
1138  SymmetricTensor<rank_, dim, Number>::
1139  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1140 {
1141  data = t.data;
1142  return *this;
1143 }
1144 
1145 
1146 
1147 template <int rank_, int dim, typename Number>
1148 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1149  SymmetricTensor<rank_, dim, Number> &
1150  SymmetricTensor<rank_, dim, Number>::operator=(const Number &d)
1151 {
1152  Assert(numbers::value_is_zero(d),
1153  ExcMessage("Only assignment with zero is allowed"));
1154  (void)d;
1155 
1156  data = internal::NumberType<Number>::value(0.0);
1157 
1158  return *this;
1159 }
1160 
1161 
1162 namespace internal
1163 {
1164  namespace SymmetricTensorImplementation
1165  {
1166  template <int dim, typename Number>
1167  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1168  ::Tensor<2, dim, Number>
1169  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1170  {
1171  ::Tensor<2, dim, Number> t;
1172 
1173  // diagonal entries are stored first
1174  for (unsigned int d = 0; d < dim; ++d)
1175  t[d][d] = s.access_raw_entry(d);
1176 
1177  // off-diagonal entries come next, row by row
1178  for (unsigned int d = 0, c = 0; d < dim; ++d)
1179  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1180  {
1181  t[d][e] = s.access_raw_entry(dim + c);
1182  t[e][d] = s.access_raw_entry(dim + c);
1183  }
1184  return t;
1185  }
1186 
1187 
1188  template <int dim, typename Number>
1189  DEAL_II_CONSTEXPR ::Tensor<4, dim, Number>
1190  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1191  {
1192  // utilize the symmetry properties of SymmetricTensor<4,dim>
1193  // discussed in the class documentation to avoid accessing all
1194  // independent elements of the input tensor more than once
1195  ::Tensor<4, dim, Number> t;
1196 
1197  for (unsigned int i = 0; i < dim; ++i)
1198  for (unsigned int j = i; j < dim; ++j)
1199  for (unsigned int k = 0; k < dim; ++k)
1200  for (unsigned int l = k; l < dim; ++l)
1201  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1202  t[TableIndices<4>(j, i, k, l)] =
1203  t[TableIndices<4>(j, i, l, k)] =
1204  st[TableIndices<4>(i, j, k, l)];
1205 
1206  return t;
1207  }
1208 
1209 
1210  template <typename Number>
1211  struct Inverse<2, 1, Number>
1212  {
1213  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1214  ::SymmetricTensor<2, 1, Number>
1215  value(const ::SymmetricTensor<2, 1, Number> &t)
1216  {
1217  ::SymmetricTensor<2, 1, Number> tmp;
1218 
1219  tmp[0][0] = 1.0 / t[0][0];
1220 
1221  return tmp;
1222  }
1223  };
1224 
1225 
1226  template <typename Number>
1227  struct Inverse<2, 2, Number>
1228  {
1229  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1230  ::SymmetricTensor<2, 2, Number>
1231  value(const ::SymmetricTensor<2, 2, Number> &t)
1232  {
1233  ::SymmetricTensor<2, 2, Number> tmp;
1234 
1235  // Sympy result: ([
1236  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1237  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1238  const TableIndices<2> idx_00(0, 0);
1239  const TableIndices<2> idx_01(0, 1);
1240  const TableIndices<2> idx_11(1, 1);
1241  const Number inv_det_t =
1242  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1243  tmp[idx_00] = t[idx_11];
1244  tmp[idx_01] = -t[idx_01];
1245  tmp[idx_11] = t[idx_00];
1246  tmp *= inv_det_t;
1247 
1248  return tmp;
1249  }
1250  };
1251 
1252 
1253  template <typename Number>
1254  struct Inverse<2, 3, Number>
1255  {
1256  DEAL_II_CONSTEXPR static ::SymmetricTensor<2, 3, Number>
1257  value(const ::SymmetricTensor<2, 3, Number> &t)
1258  {
1259  ::SymmetricTensor<2, 3, Number> tmp;
1260 
1261  // Sympy result: ([
1262  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1263  // 2*t01*t02*t12 - t02**2*t11),
1264  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1265  // 2*t01*t02*t12 - t02**2*t11),
1266  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1267  // 2*t01*t02*t12 - t02**2*t11)],
1268  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1269  // 2*t01*t02*t12 - t02**2*t11),
1270  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1271  // 2*t01*t02*t12 - t02**2*t11),
1272  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1273  // 2*t01*t02*t12 + t02**2*t11)],
1274  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1275  // 2*t01*t02*t12 - t02**2*t11),
1276  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1277  // 2*t01*t02*t12 + t02**2*t11),
1278  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1279  // 2*t01*t02*t12 + t02**2*t11)] ])
1280  //
1281  // =
1282  //
1283  // [ (t11*t22 - t12**2)/det_t,
1284  // (-t01*t22 + t02*t12)/det_t,
1285  // (t01*t12 - t02*t11)/det_t],
1286  // [ (-t01*t22 + t02*t12)/det_t,
1287  // (t00*t22 - t02**2)/det_t,
1288  // (-t00*t12 + t01*t02)/det_t],
1289  // [ (t01*t12 - t02*t11)/det_t,
1290  // (-t00*t12 + t01*t02)/det_t,
1291  // (t00*t11 - t01**2)/det_t] ])
1292  //
1293  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1294  // 2*t01*t02*t12 - t02**2*t11)
1295  const TableIndices<2> idx_00(0, 0);
1296  const TableIndices<2> idx_01(0, 1);
1297  const TableIndices<2> idx_02(0, 2);
1298  const TableIndices<2> idx_11(1, 1);
1299  const TableIndices<2> idx_12(1, 2);
1300  const TableIndices<2> idx_22(2, 2);
1301  const Number inv_det_t =
1302  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1303  t[idx_00] * t[idx_12] * t[idx_12] -
1304  t[idx_01] * t[idx_01] * t[idx_22] +
1305  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1306  t[idx_02] * t[idx_02] * t[idx_11]);
1307  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1308  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1309  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1310  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1311  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1312  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1313  tmp *= inv_det_t;
1314 
1315  return tmp;
1316  }
1317  };
1318 
1319 
1320  template <typename Number>
1321  struct Inverse<4, 1, Number>
1322  {
1323  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 1, Number>
1324  value(const ::SymmetricTensor<4, 1, Number> &t)
1325  {
1326  ::SymmetricTensor<4, 1, Number> tmp;
1327  tmp.data[0][0] = 1.0 / t.data[0][0];
1328  return tmp;
1329  }
1330  };
1331 
1332 
1333  template <typename Number>
1334  struct Inverse<4, 2, Number>
1335  {
1336  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 2, Number>
1337  value(const ::SymmetricTensor<4, 2, Number> &t)
1338  {
1339  ::SymmetricTensor<4, 2, Number> tmp;
1340 
1341  // Inverting this tensor is a little more complicated than necessary,
1342  // since we store the data of 't' as a 3x3 matrix t.data, but the
1343  // product between a rank-4 and a rank-2 tensor is really not the
1344  // product between this matrix and the 3-vector of a rhs, but rather
1345  //
1346  // B.vec = t.data * mult * A.vec
1347  //
1348  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1349  // capture the fact that we need to add up both the c_ij12*a_12 and the
1350  // c_ij21*a_21 terms.
1351  //
1352  // In addition, in this scheme, the identity tensor has the matrix
1353  // representation mult^-1.
1354  //
1355  // The inverse of 't' therefore has the matrix representation
1356  //
1357  // inv.data = mult^-1 * t.data^-1 * mult^-1
1358  //
1359  // in order to compute it, let's first compute the inverse of t.data and
1360  // put it into tmp.data; at the end of the function we then scale the
1361  // last row and column of the inverse by 1/2, corresponding to the left
1362  // and right multiplication with mult^-1.
1363  const Number t4 = t.data[0][0] * t.data[1][1],
1364  t6 = t.data[0][0] * t.data[1][2],
1365  t8 = t.data[0][1] * t.data[1][0],
1366  t00 = t.data[0][2] * t.data[1][0],
1367  t01 = t.data[0][1] * t.data[2][0],
1368  t04 = t.data[0][2] * t.data[2][0],
1369  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1370  t8 * t.data[2][2] + t00 * t.data[2][1] +
1371  t01 * t.data[1][2] - t04 * t.data[1][1]);
1372  tmp.data[0][0] =
1373  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1374  tmp.data[0][1] =
1375  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1376  tmp.data[0][2] =
1377  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1378  tmp.data[1][0] =
1379  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1380  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1381  tmp.data[1][2] = -(t6 - t00) * t07;
1382  tmp.data[2][0] =
1383  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1384  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1385  tmp.data[2][2] = (t4 - t8) * t07;
1386 
1387  // scale last row and column as mentioned
1388  // above
1389  tmp.data[2][0] /= 2;
1390  tmp.data[2][1] /= 2;
1391  tmp.data[0][2] /= 2;
1392  tmp.data[1][2] /= 2;
1393  tmp.data[2][2] /= 4;
1394 
1395  return tmp;
1396  }
1397  };
1398 
1399 
1400  template <typename Number>
1401  struct Inverse<4, 3, Number>
1402  {
1403  static ::SymmetricTensor<4, 3, Number>
1404  value(const ::SymmetricTensor<4, 3, Number> &t)
1405  {
1406  ::SymmetricTensor<4, 3, Number> tmp = t;
1407 
1408  // This function follows the exact same scheme as the 2d case, except
1409  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1410  // Instead, we use the Gauss-Jordan algorithm implemented for
1411  // FullMatrix. For historical reasons the following code is copied from
1412  // there, with the tangential benefit that we do not need to copy the
1413  // tensor entries to and from the FullMatrix.
1414  const unsigned int N = 6;
1415 
1416  // First get an estimate of the size of the elements of this matrix,
1417  // for later checks whether the pivot element is large enough, or
1418  // whether we have to fear that the matrix is not regular.
1419  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1420  for (unsigned int i = 0; i < N; ++i)
1421  diagonal_sum += std::fabs(tmp.data[i][i]);
1422  const Number typical_diagonal_element =
1423  diagonal_sum / static_cast<double>(N);
1424  (void)typical_diagonal_element;
1425 
1426  unsigned int p[N];
1427  for (unsigned int i = 0; i < N; ++i)
1428  p[i] = i;
1429 
1430  for (unsigned int j = 0; j < N; ++j)
1431  {
1432  // Pivot search: search that part of the line on and right of the
1433  // diagonal for the largest element.
1434  Number max = std::fabs(tmp.data[j][j]);
1435  unsigned int r = j;
1436  for (unsigned int i = j + 1; i < N; ++i)
1437  if (std::fabs(tmp.data[i][j]) > max)
1438  {
1439  max = std::fabs(tmp.data[i][j]);
1440  r = i;
1441  }
1442 
1443  // Check whether the pivot is too small
1444  Assert(max > 1.e-16 * typical_diagonal_element,
1445  ExcMessage("This tensor seems to be noninvertible"));
1446 
1447  // Row interchange
1448  if (r > j)
1449  {
1450  for (unsigned int k = 0; k < N; ++k)
1451  std::swap(tmp.data[j][k], tmp.data[r][k]);
1452 
1453  std::swap(p[j], p[r]);
1454  }
1455 
1456  // Transformation
1457  const Number hr = 1. / tmp.data[j][j];
1458  tmp.data[j][j] = hr;
1459  for (unsigned int k = 0; k < N; ++k)
1460  {
1461  if (k == j)
1462  continue;
1463  for (unsigned int i = 0; i < N; ++i)
1464  {
1465  if (i == j)
1466  continue;
1467  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1468  }
1469  }
1470  for (unsigned int i = 0; i < N; ++i)
1471  {
1472  tmp.data[i][j] *= hr;
1473  tmp.data[j][i] *= -hr;
1474  }
1475  tmp.data[j][j] = hr;
1476  }
1477 
1478  // Column interchange
1479  Number hv[N];
1480  for (unsigned int i = 0; i < N; ++i)
1481  {
1482  for (unsigned int k = 0; k < N; ++k)
1483  hv[p[k]] = tmp.data[i][k];
1484  for (unsigned int k = 0; k < N; ++k)
1485  tmp.data[i][k] = hv[k];
1486  }
1487 
1488  // Scale rows and columns. The mult matrix
1489  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1490  for (unsigned int i = 3; i < 6; ++i)
1491  for (unsigned int j = 0; j < 3; ++j)
1492  tmp.data[i][j] /= 2;
1493 
1494  for (unsigned int i = 0; i < 3; ++i)
1495  for (unsigned int j = 3; j < 6; ++j)
1496  tmp.data[i][j] /= 2;
1497 
1498  for (unsigned int i = 3; i < 6; ++i)
1499  for (unsigned int j = 3; j < 6; ++j)
1500  tmp.data[i][j] /= 4;
1501 
1502  return tmp;
1503  }
1504  };
1505 
1506  } // namespace SymmetricTensorImplementation
1507 } // namespace internal
1508 
1509 
1510 
1511 template <int rank_, int dim, typename Number>
1512 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>::
1513  operator Tensor<rank_, dim, Number>() const
1514 {
1515  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1516 }
1517 
1518 
1519 
1520 template <int rank_, int dim, typename Number>
1521 constexpr bool
1522 SymmetricTensor<rank_, dim, Number>::
1523 operator==(const SymmetricTensor<rank_, dim, Number> &t) const
1524 {
1525  return data == t.data;
1526 }
1527 
1528 
1529 
1530 template <int rank_, int dim, typename Number>
1531 constexpr bool
1532 SymmetricTensor<rank_, dim, Number>::
1533 operator!=(const SymmetricTensor<rank_, dim, Number> &t) const
1534 {
1535  return data != t.data;
1536 }
1537 
1538 
1539 
1540 template <int rank_, int dim, typename Number>
1541 template <typename OtherNumber>
1542 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1543  SymmetricTensor<rank_, dim, Number> &
1544  SymmetricTensor<rank_, dim, Number>::
1545  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1546 {
1547  data += t.data;
1548  return *this;
1549 }
1550 
1551 
1552 
1553 template <int rank_, int dim, typename Number>
1554 template <typename OtherNumber>
1555 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1556  SymmetricTensor<rank_, dim, Number> &
1557  SymmetricTensor<rank_, dim, Number>::
1558  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1559 {
1560  data -= t.data;
1561  return *this;
1562 }
1563 
1564 
1565 
1566 template <int rank_, int dim, typename Number>
1567 template <typename OtherNumber>
1568 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1569  SymmetricTensor<rank_, dim, Number> &
1570  SymmetricTensor<rank_, dim, Number>::operator*=(const OtherNumber &d)
1571 {
1572  data *= d;
1573  return *this;
1574 }
1575 
1576 
1577 
1578 template <int rank_, int dim, typename Number>
1579 template <typename OtherNumber>
1580 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1581  SymmetricTensor<rank_, dim, Number> &
1582  SymmetricTensor<rank_, dim, Number>::operator/=(const OtherNumber &d)
1583 {
1584  data /= d;
1585  return *this;
1586 }
1587 
1588 
1589 
1590 template <int rank_, int dim, typename Number>
1591 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1592  SymmetricTensor<rank_, dim, Number>
1593  SymmetricTensor<rank_, dim, Number>::operator-() const
1594 {
1595  SymmetricTensor tmp = *this;
1596  tmp.data = -tmp.data;
1597  return tmp;
1598 }
1599 
1600 
1601 
1602 template <int rank_, int dim, typename Number>
1603 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void
1604 SymmetricTensor<rank_, dim, Number>::clear()
1605 {
1606  data.clear();
1607 }
1608 
1609 
1610 
1611 template <int rank_, int dim, typename Number>
1612 constexpr std::size_t
1613 SymmetricTensor<rank_, dim, Number>::memory_consumption()
1614 {
1615  // all memory consists of statically allocated memory of the current
1616  // object, no pointers
1617  return sizeof(SymmetricTensor<rank_, dim, Number>);
1618 }
1619 
1620 
1621 
1622 namespace internal
1623 {
1624  template <int dim, typename Number, typename OtherNumber = Number>
1625  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1626  typename SymmetricTensorAccessors::
1627  double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1628  perform_double_contraction(
1629  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1630  base_tensor_type &data,
1631  const typename SymmetricTensorAccessors::
1632  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1633  {
1634  using result_type = typename SymmetricTensorAccessors::
1635  double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1636 
1637  switch (dim)
1638  {
1639  case 1:
1640  return data[0] * sdata[0];
1641  default:
1642  // Start with the non-diagonal part to avoid some multiplications by
1643  // 2.
1644 
1645  result_type sum = data[dim] * sdata[dim];
1646  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1647  sum += data[d] * sdata[d];
1648  sum += sum; // sum = sum * 2.;
1649 
1650  // Now add the contributions from the diagonal
1651  for (unsigned int d = 0; d < dim; ++d)
1652  sum += data[d] * sdata[d];
1653  return sum;
1654  }
1655  }
1656 
1657 
1658 
1659  template <int dim, typename Number, typename OtherNumber = Number>
1660  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1661  typename SymmetricTensorAccessors::
1662  double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1663  perform_double_contraction(
1664  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1665  base_tensor_type &data,
1666  const typename SymmetricTensorAccessors::
1667  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1668  {
1669  using result_type = typename SymmetricTensorAccessors::
1670  double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1671  using value_type = typename SymmetricTensorAccessors::
1672  double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1673 
1674  const unsigned int data_dim = SymmetricTensorAccessors::
1675  StorageType<2, dim, value_type>::n_independent_components;
1676  value_type tmp[data_dim]{};
1677  for (unsigned int i = 0; i < data_dim; ++i)
1678  tmp[i] =
1679  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1680  return result_type(tmp);
1681  }
1682 
1683 
1684 
1685  template <int dim, typename Number, typename OtherNumber = Number>
1686  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1687  typename SymmetricTensorAccessors::StorageType<
1688  2,
1689  dim,
1690  typename SymmetricTensorAccessors::
1691  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1692  base_tensor_type
1693  perform_double_contraction(
1694  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1695  base_tensor_type &data,
1696  const typename SymmetricTensorAccessors::
1697  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1698  {
1699  using value_type = typename SymmetricTensorAccessors::
1700  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1701  using base_tensor_type = typename SymmetricTensorAccessors::
1702  StorageType<2, dim, value_type>::base_tensor_type;
1703 
1704  base_tensor_type tmp;
1705  for (unsigned int i = 0; i < tmp.dimension; ++i)
1706  {
1707  // Start with the non-diagonal part
1708  value_type sum = data[dim] * sdata[dim][i];
1709  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1710  sum += data[d] * sdata[d][i];
1711  sum += sum; // sum = sum * 2.;
1712 
1713  // Now add the contributions from the diagonal
1714  for (unsigned int d = 0; d < dim; ++d)
1715  sum += data[d] * sdata[d][i];
1716  tmp[i] = sum;
1717  }
1718  return tmp;
1719  }
1720 
1721 
1722 
1723  template <int dim, typename Number, typename OtherNumber = Number>
1724  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1725  typename SymmetricTensorAccessors::StorageType<
1726  4,
1727  dim,
1728  typename SymmetricTensorAccessors::
1729  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1730  base_tensor_type
1731  perform_double_contraction(
1732  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1733  base_tensor_type &data,
1734  const typename SymmetricTensorAccessors::
1735  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1736  {
1737  using value_type = typename SymmetricTensorAccessors::
1738  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1739  using base_tensor_type = typename SymmetricTensorAccessors::
1740  StorageType<4, dim, value_type>::base_tensor_type;
1741 
1742  const unsigned int data_dim = SymmetricTensorAccessors::
1743  StorageType<2, dim, value_type>::n_independent_components;
1744  base_tensor_type tmp;
1745  for (unsigned int i = 0; i < data_dim; ++i)
1746  for (unsigned int j = 0; j < data_dim; ++j)
1747  {
1748  // Start with the non-diagonal part
1749  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1750  tmp[i][j] += data[i][d] * sdata[d][j];
1751  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1752 
1753  // Now add the contributions from the diagonal
1754  for (unsigned int d = 0; d < dim; ++d)
1755  tmp[i][j] += data[i][d] * sdata[d][j];
1756  }
1757  return tmp;
1758  }
1759 
1760 } // end of namespace internal
1761 
1762 
1763 
1764 template <int rank_, int dim, typename Number>
1765 template <typename OtherNumber>
1766 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1767  typename internal::SymmetricTensorAccessors::
1768  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1769  SymmetricTensor<rank_, dim, Number>::
1770  operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const
1771 {
1772  // need to have two different function calls
1773  // because a scalar and rank-2 tensor are not
1774  // the same data type (see internal function
1775  // above)
1776  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1777  s.data);
1778 }
1779 
1780 
1781 
1782 template <int rank_, int dim, typename Number>
1783 template <typename OtherNumber>
1784 DEAL_II_CONSTEXPR inline typename internal::SymmetricTensorAccessors::
1785  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1786  SymmetricTensor<rank_, dim, Number>::
1787  operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const
1788 {
1789  typename internal::SymmetricTensorAccessors::
1790  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1791  tmp.data =
1792  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1793  s.data);
1794  return tmp;
1795 }
1796 
1797 
1798 
1799 // internal namespace to switch between the
1800 // access of different tensors. There used to
1801 // be explicit instantiations before for
1802 // different ranks and dimensions, but since
1803 // we now allow for templates on the data
1804 // type, and since we cannot partially
1805 // specialize the implementation, this got
1806 // into a separate namespace
1807 namespace internal
1808 {
1809  // The variables within this struct will be referenced in the next functions.
1810  // It is a workaround that allows returning a reference to a static variable
1811  // while allowing constexpr evaluation of the function.
1812  // It has to be defined outside the function because constexpr functions
1813  // cannot define static variables.
1814  // A similar struct has also been defined in tensor.h
1815  template <typename Type>
1816  struct Uninitialized
1817  {
1818  static Type value;
1819  };
1820 
1821  template <typename Type>
1822  Type Uninitialized<Type>::value;
1823 
1824  template <int dim, typename Number>
1825  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1826  symmetric_tensor_access(const TableIndices<2> &indices,
1827  typename SymmetricTensorAccessors::
1828  StorageType<2, dim, Number>::base_tensor_type &data)
1829  {
1830  // 1d is very simple and done first
1831  if (dim == 1)
1832  return data[0];
1833 
1834  // first treat the main diagonal elements, which are stored consecutively
1835  // at the beginning
1836  if (indices[0] == indices[1])
1837  return data[indices[0]];
1838 
1839  // the rest is messier and requires a few switches.
1840  switch (dim)
1841  {
1842  case 2:
1843  // at least for the 2x2 case it is reasonably simple
1844  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1845  ((indices[0] == 0) && (indices[1] == 1)),
1846  ExcInternalError());
1847  return data[2];
1848 
1849  default:
1850  // to do the rest, sort our indices before comparing
1851  {
1852  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1853  std::max(indices[0], indices[1]));
1854  for (unsigned int d = 0, c = 0; d < dim; ++d)
1855  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1856  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1857  return data[dim + c];
1858  Assert(false, ExcInternalError());
1859  }
1860  }
1861 
1862  // The code should never reach there.
1863  // Returns a dummy reference to a dummy variable just to make the
1864  // compiler happy.
1865  return Uninitialized<Number>::value;
1866  }
1867 
1868 
1869 
1870  template <int dim, typename Number>
1871  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1872  symmetric_tensor_access(const TableIndices<2> &indices,
1873  const typename SymmetricTensorAccessors::
1874  StorageType<2, dim, Number>::base_tensor_type &data)
1875  {
1876  // 1d is very simple and done first
1877  if (dim == 1)
1878  return data[0];
1879 
1880  // first treat the main diagonal elements, which are stored consecutively
1881  // at the beginning
1882  if (indices[0] == indices[1])
1883  return data[indices[0]];
1884 
1885  // the rest is messier and requires a few switches.
1886  switch (dim)
1887  {
1888  case 2:
1889  // at least for the 2x2 case it is reasonably simple
1890  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1891  ((indices[0] == 0) && (indices[1] == 1)),
1892  ExcInternalError());
1893  return data[2];
1894 
1895  default:
1896  // to do the rest, sort our indices before comparing
1897  {
1898  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1899  std::max(indices[0], indices[1]));
1900  for (unsigned int d = 0, c = 0; d < dim; ++d)
1901  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1902  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1903  return data[dim + c];
1904  Assert(false, ExcInternalError());
1905  }
1906  }
1907 
1908  // The code should never reach there.
1909  // Returns a dummy reference to a dummy variable just to make the
1910  // compiler happy.
1911  return Uninitialized<Number>::value;
1912  }
1913 
1914 
1915 
1916  template <int dim, typename Number>
1917  DEAL_II_CONSTEXPR inline Number &
1918  symmetric_tensor_access(const TableIndices<4> &indices,
1919  typename SymmetricTensorAccessors::
1920  StorageType<4, dim, Number>::base_tensor_type &data)
1921  {
1922  switch (dim)
1923  {
1924  case 1:
1925  return data[0][0];
1926 
1927  case 2:
1928  // each entry of the tensor can be thought of as an entry in a
1929  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1930  // rank-2 tensors. this is the format in which we store rank-4
1931  // tensors. determine which position the present entry is
1932  // stored in
1933  {
1934  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1935  return data[base_index[indices[0]][indices[1]]]
1936  [base_index[indices[2]][indices[3]]];
1937  }
1938  case 3:
1939  // each entry of the tensor can be thought of as an entry in a
1940  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1941  // rank-2 tensors. this is the format in which we store rank-4
1942  // tensors. determine which position the present entry is
1943  // stored in
1944  {
1945  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1946  {3, 1, 5},
1947  {4, 5, 2}};
1948  return data[base_index[indices[0]][indices[1]]]
1949  [base_index[indices[2]][indices[3]]];
1950  }
1951 
1952  default:
1953  Assert(false, ExcNotImplemented());
1954  }
1955 
1956  // The code should never reach there.
1957  // Returns a dummy reference to a dummy variable just to make the
1958  // compiler happy.
1959  return Uninitialized<Number>::value;
1960  }
1961 
1962 
1963  template <int dim, typename Number>
1964  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1965  symmetric_tensor_access(const TableIndices<4> &indices,
1966  const typename SymmetricTensorAccessors::
1967  StorageType<4, dim, Number>::base_tensor_type &data)
1968  {
1969  switch (dim)
1970  {
1971  case 1:
1972  return data[0][0];
1973 
1974  case 2:
1975  // each entry of the tensor can be thought of as an entry in a
1976  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1977  // rank-2 tensors. this is the format in which we store rank-4
1978  // tensors. determine which position the present entry is
1979  // stored in
1980  {
1981  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1982  return data[base_index[indices[0]][indices[1]]]
1983  [base_index[indices[2]][indices[3]]];
1984  }
1985  case 3:
1986  // each entry of the tensor can be thought of as an entry in a
1987  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1988  // rank-2 tensors. this is the format in which we store rank-4
1989  // tensors. determine which position the present entry is
1990  // stored in
1991  {
1992  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1993  {3, 1, 5},
1994  {4, 5, 2}};
1995  return data[base_index[indices[0]][indices[1]]]
1996  [base_index[indices[2]][indices[3]]];
1997  }
1998 
1999  default:
2000  Assert(false, ExcNotImplemented());
2001  }
2002 
2003  // The code should never reach there.
2004  // Returns a dummy reference to a dummy variable just to make the
2005  // compiler happy.
2006  return Uninitialized<Number>::value;
2007  }
2008 
2009 } // end of namespace internal
2010 
2011 
2012 
2013 template <int rank_, int dim, typename Number>
2014 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
2015  SymmetricTensor<rank_, dim, Number>::
2016  operator()(const TableIndices<rank_> &indices)
2017 {
2018  for (unsigned int r = 0; r < rank; ++r)
2019  AssertIndexRange(indices[r], dimension);
2020  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2021 }
2022 
2023 
2024 
2025 template <int rank_, int dim, typename Number>
2026 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
2027  SymmetricTensor<rank_, dim, Number>::
2028  operator()(const TableIndices<rank_> &indices) const
2029 {
2030  for (unsigned int r = 0; r < rank; ++r)
2031  AssertIndexRange(indices[r], dimension);
2032  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2033 }
2034 
2035 
2036 
2037 namespace internal
2038 {
2039  namespace SymmetricTensorImplementation
2040  {
2041  template <int rank_>
2042  constexpr TableIndices<rank_>
2043  get_partially_filled_indices(const unsigned int row,
2044  const std::integral_constant<int, 2> &)
2045  {
2046  return TableIndices<rank_>(row, numbers::invalid_unsigned_int);
2047  }
2048 
2049 
2050  template <int rank_>
2051  constexpr TableIndices<rank_>
2052  get_partially_filled_indices(const unsigned int row,
2053  const std::integral_constant<int, 4> &)
2054  {
2055  return TableIndices<rank_>(row,
2056  numbers::invalid_unsigned_int,
2057  numbers::invalid_unsigned_int,
2058  numbers::invalid_unsigned_int);
2059  }
2060  } // namespace SymmetricTensorImplementation
2061 } // namespace internal
2062 
2063 
2064 template <int rank_, int dim, typename Number>
2065 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2066  Accessor<rank_, dim, true, rank_ - 1, Number>
2067  SymmetricTensor<rank_, dim, Number>::
2068  operator[](const unsigned int row) const
2069 {
2070  return internal::SymmetricTensorAccessors::
2071  Accessor<rank_, dim, true, 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 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE internal::
2081  SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2082  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2083 {
2084  return internal::SymmetricTensorAccessors::
2085  Accessor<rank_, dim, false, rank_ - 1, Number>(
2086  *this,
2087  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2088  rank_>(row, std::integral_constant<int, rank_>()));
2089 }
2090 
2091 
2092 
2093 template <int rank_, int dim, typename Number>
2094 constexpr DEAL_II_ALWAYS_INLINE const Number &
2095  SymmetricTensor<rank_, dim, Number>::
2096  operator[](const TableIndices<rank_> &indices) const
2097 {
2098  return operator()(indices);
2099 }
2100 
2101 
2102 
2103 template <int rank_, int dim, typename Number>
2104 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
2105  SymmetricTensor<rank_, dim, Number>::
2106  operator[](const TableIndices<rank_> &indices)
2107 {
2108  return operator()(indices);
2109 }
2110 
2111 
2112 
2113 template <int rank_, int dim, typename Number>
2114 inline Number *
2115 SymmetricTensor<rank_, dim, Number>::begin_raw()
2116 {
2117  return std::addressof(this->access_raw_entry(0));
2118 }
2119 
2120 
2121 
2122 template <int rank_, int dim, typename Number>
2123 inline const Number *
2124 SymmetricTensor<rank_, dim, Number>::begin_raw() const
2125 {
2126  return std::addressof(this->access_raw_entry(0));
2127 }
2128 
2129 
2130 
2131 template <int rank_, int dim, typename Number>
2132 inline Number *
2133 SymmetricTensor<rank_, dim, Number>::end_raw()
2134 {
2135  return begin_raw() + n_independent_components;
2136 }
2137 
2138 
2139 
2140 template <int rank_, int dim, typename Number>
2141 inline const Number *
2142 SymmetricTensor<rank_, dim, Number>::end_raw() const
2143 {
2144  return begin_raw() + n_independent_components;
2145 }
2146 
2147 
2148 
2149 namespace internal
2150 {
2151  namespace SymmetricTensorImplementation
2152  {
2153  template <int dim, typename Number>
2154  constexpr unsigned int
2155  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2156  const unsigned int index)
2157  {
2158  return index;
2159  }
2160 
2161 
2162  template <int dim, typename Number>
2163  constexpr ::TableIndices<2>
2164  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2165  const unsigned int index)
2166  {
2167  return internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
2168  base_tensor_type::unrolled_to_component_indices(index);
2169  }
2170 
2171  } // namespace SymmetricTensorImplementation
2172 } // namespace internal
2173 
2174 
2175 
2176 template <int rank_, int dim, typename Number>
2177 DEAL_II_CONSTEXPR inline const Number &
2178 SymmetricTensor<rank_, dim, Number>::access_raw_entry(
2179  const unsigned int index) const
2180 {
2181  AssertIndexRange(index, n_independent_components);
2182  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2183  index)];
2184 }
2185 
2186 
2187 
2188 template <int rank_, int dim, typename Number>
2189 DEAL_II_CONSTEXPR inline Number &
2190 SymmetricTensor<rank_, dim, Number>::access_raw_entry(const unsigned int index)
2191 {
2192  AssertIndexRange(index, n_independent_components);
2193  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2194  index)];
2195 }
2196 
2197 
2198 
2199 namespace internal
2200 {
2201  template <int dim, typename Number>
2202  DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
2203  compute_norm(const typename SymmetricTensorAccessors::
2204  StorageType<2, dim, Number>::base_tensor_type &data)
2205  {
2206  switch (dim)
2207  {
2208  case 1:
2209  return numbers::NumberTraits<Number>::abs(data[0]);
2210 
2211  case 2:
2212  return std::sqrt(
2213  numbers::NumberTraits<Number>::abs_square(data[0]) +
2214  numbers::NumberTraits<Number>::abs_square(data[1]) +
2215  2. * numbers::NumberTraits<Number>::abs_square(data[2]));
2216 
2217  case 3:
2218  return std::sqrt(
2219  numbers::NumberTraits<Number>::abs_square(data[0]) +
2220  numbers::NumberTraits<Number>::abs_square(data[1]) +
2221  numbers::NumberTraits<Number>::abs_square(data[2]) +
2222  2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
2223  2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
2224  2. * numbers::NumberTraits<Number>::abs_square(data[5]));
2225 
2226  default:
2227  {
2228  typename numbers::NumberTraits<Number>::real_type return_value =
2229  typename numbers::NumberTraits<Number>::real_type();
2230 
2231  for (unsigned int d = 0; d < dim; ++d)
2232  return_value +=
2233  numbers::NumberTraits<Number>::abs_square(data[d]);
2234  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2235  return_value +=
2236  2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2237 
2238  return std::sqrt(return_value);
2239  }
2240  }
2241  }
2242 
2243 
2244 
2245  template <int dim, typename Number>
2246  DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
2247  compute_norm(const typename SymmetricTensorAccessors::
2248  StorageType<4, dim, Number>::base_tensor_type &data)
2249  {
2250  switch (dim)
2251  {
2252  case 1:
2253  return numbers::NumberTraits<Number>::abs(data[0][0]);
2254 
2255  default:
2256  {
2257  typename numbers::NumberTraits<Number>::real_type return_value =
2258  typename numbers::NumberTraits<Number>::real_type();
2259 
2260  const unsigned int n_independent_components = data.dimension;
2261 
2262  for (unsigned int i = 0; i < dim; ++i)
2263  for (unsigned int j = 0; j < dim; ++j)
2264  return_value +=
2265  numbers::NumberTraits<Number>::abs_square(data[i][j]);
2266  for (unsigned int i = 0; i < dim; ++i)
2267  for (unsigned int j = dim; j < n_independent_components; ++j)
2268  return_value +=
2269  2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2270  for (unsigned int i = dim; i < n_independent_components; ++i)
2271  for (unsigned int j = 0; j < dim; ++j)
2272  return_value +=
2273  2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2274  for (unsigned int i = dim; i < n_independent_components; ++i)
2275  for (unsigned int j = dim; j < n_independent_components; ++j)
2276  return_value +=
2277  4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2278 
2279  return std::sqrt(return_value);
2280  }
2281  }
2282  }
2283 
2284 } // end of namespace internal
2285 
2286 
2287 
2288 template <int rank_, int dim, typename Number>
2289 constexpr typename numbers::NumberTraits<Number>::real_type
2290 SymmetricTensor<rank_, dim, Number>::norm() const
2291 {
2292  return internal::compute_norm<dim, Number>(data);
2293 }
2294 
2295 
2296 
2297 namespace internal
2298 {
2299  namespace SymmetricTensorImplementation
2300  {
2301  // a function to do the unrolling from a set of indices to a
2302  // scalar index into the array in which we store the elements of
2303  // a symmetric tensor
2304  //
2305  // this function is for rank-2 tensors
2306  template <int dim>
2307  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE unsigned int
2308  component_to_unrolled_index(const TableIndices<2> &indices)
2309  {
2310  AssertIndexRange(indices[0], dim);
2311  AssertIndexRange(indices[1], dim);
2312 
2313  switch (dim)
2314  {
2315  case 1:
2316  {
2317  return 0;
2318  }
2319 
2320  case 2:
2321  {
2322  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2323  return table[indices[0]][indices[1]];
2324  }
2325 
2326  case 3:
2327  {
2328  constexpr unsigned int table[3][3] = {{0, 3, 4},
2329  {3, 1, 5},
2330  {4, 5, 2}};
2331  return table[indices[0]][indices[1]];
2332  }
2333 
2334  case 4:
2335  {
2336  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2337  {4, 1, 7, 8},
2338  {5, 7, 2, 9},
2339  {6, 8, 9, 3}};
2340  return table[indices[0]][indices[1]];
2341  }
2342 
2343  default:
2344  // for the remainder, manually figure out the numbering
2345  {
2346  if (indices[0] == indices[1])
2347  return indices[0];
2348 
2349  TableIndices<2> sorted_indices(indices);
2350  sorted_indices.sort();
2351 
2352  for (unsigned int d = 0, c = 0; d < dim; ++d)
2353  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2354  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2355  return dim + c;
2356 
2357  // should never get here:
2358  Assert(false, ExcInternalError());
2359  return 0;
2360  }
2361  }
2362  }
2363 
2364  // a function to do the unrolling from a set of indices to a
2365  // scalar index into the array in which we store the elements of
2366  // a symmetric tensor
2367  //
2368  // this function is for tensors of ranks not already handled
2369  // above
2370  template <int dim, int rank_>
2371  DEAL_II_CONSTEXPR inline unsigned int
2372  component_to_unrolled_index(const TableIndices<rank_> &indices)
2373  {
2374  (void)indices;
2375  Assert(false, ExcNotImplemented());
2376  return numbers::invalid_unsigned_int;
2377  }
2378  } // namespace SymmetricTensorImplementation
2379 } // namespace internal
2380 
2381 
2382 template <int rank_, int dim, typename Number>
2383 constexpr unsigned int
2384 SymmetricTensor<rank_, dim, Number>::component_to_unrolled_index(
2385  const TableIndices<rank_> &indices)
2386 {
2387  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2388  dim>(indices);
2389 }
2390 
2391 
2392 
2393 namespace internal
2394 {
2395  namespace SymmetricTensorImplementation
2396  {
2397  // a function to do the inverse of the unrolling from a set of
2398  // indices to a scalar index into the array in which we store
2399  // the elements of a symmetric tensor. in other words, it goes
2400  // from the scalar index into the array to a set of indices of
2401  // the tensor
2402  //
2403  // this function is for rank-2 tensors
2404  template <int dim>
2405  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2406  unrolled_to_component_indices(const unsigned int i,
2407  const std::integral_constant<int, 2> &)
2408  {
2409  Assert(
2410  (i < ::SymmetricTensor<2, dim, double>::n_independent_components),
2411  ExcIndexRange(
2412  i,
2413  0,
2414  ::SymmetricTensor<2, dim, double>::n_independent_components));
2415  switch (dim)
2416  {
2417  case 1:
2418  {
2419  return {0, 0};
2420  }
2421 
2422  case 2:
2423  {
2424  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2425  TableIndices<2>(1, 1),
2426  TableIndices<2>(0, 1)};
2427  return table[i];
2428  }
2429 
2430  case 3:
2431  {
2432  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2433  TableIndices<2>(1, 1),
2434  TableIndices<2>(2, 2),
2435  TableIndices<2>(0, 1),
2436  TableIndices<2>(0, 2),
2437  TableIndices<2>(1, 2)};
2438  return table[i];
2439  }
2440 
2441  default:
2442  if (i < dim)
2443  return {i, i};
2444 
2445  for (unsigned int d = 0, c = 0; d < dim; ++d)
2446  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2447  if (c == i)
2448  return {d, e};
2449 
2450  // should never get here:
2451  Assert(false, ExcInternalError());
2452  return {0, 0};
2453  }
2454  }
2455 
2456  // a function to do the inverse of the unrolling from a set of
2457  // indices to a scalar index into the array in which we store
2458  // the elements of a symmetric tensor. in other words, it goes
2459  // from the scalar index into the array to a set of indices of
2460  // the tensor
2461  //
2462  // this function is for tensors of a rank not already handled
2463  // above
2464  template <int dim, int rank_>
2465  DEAL_II_CONSTEXPR inline
2466  typename std::enable_if<rank_ != 2, TableIndices<rank_>>::type
2467  unrolled_to_component_indices(const unsigned int i,
2468  const std::integral_constant<int, rank_> &)
2469  {
2470  (void)i;
2471  Assert(
2472  (i <
2473  ::SymmetricTensor<rank_, dim, double>::n_independent_components),
2474  ExcIndexRange(i,
2475  0,
2476  ::SymmetricTensor<rank_, dim, double>::
2477  n_independent_components));
2478  Assert(false, ExcNotImplemented());
2479  return TableIndices<rank_>();
2480  }
2481 
2482  } // namespace SymmetricTensorImplementation
2483 } // namespace internal
2484 
2485 template <int rank_, int dim, typename Number>
2486 constexpr DEAL_II_ALWAYS_INLINE TableIndices<rank_>
2487  SymmetricTensor<rank_, dim, Number>::unrolled_to_component_indices(
2488  const unsigned int i)
2489 {
2490  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2491  dim>(i, std::integral_constant<int, rank_>());
2492 }
2493 
2494 
2495 
2496 template <int rank_, int dim, typename Number>
2497 template <class Archive>
2498 inline void
2499 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2500 {
2501  ar &data;
2502 }
2503 
2504 
2505 #endif // DOXYGEN
2506 
2507 /* ----------------- Non-member functions operating on tensors. ------------ */
2508 
2509 
2522 template <int rank_, int dim, typename Number, typename OtherNumber>
2523 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2524  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2525  operator+(const SymmetricTensor<rank_, dim, Number> & left,
2526  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2527 {
2528  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2529  tmp = left;
2530  tmp += right;
2531  return tmp;
2532 }
2533 
2534 
2547 template <int rank_, int dim, typename Number, typename OtherNumber>
2548 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2549  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2550  operator-(const SymmetricTensor<rank_, dim, Number> & left,
2551  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2552 {
2553  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2554  tmp = left;
2555  tmp -= right;
2556  return tmp;
2557 }
2558 
2559 
2567 template <int rank_, int dim, typename Number, typename OtherNumber>
2568 constexpr DEAL_II_ALWAYS_INLINE
2569  Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2570  operator+(const SymmetricTensor<rank_, dim, Number> &left,
2571  const Tensor<rank_, dim, OtherNumber> & right)
2572 {
2573  return Tensor<rank_, dim, Number>(left) + right;
2574 }
2575 
2576 
2584 template <int rank_, int dim, typename Number, typename OtherNumber>
2585 constexpr DEAL_II_ALWAYS_INLINE
2586  Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2587  operator+(const Tensor<rank_, dim, Number> & left,
2588  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2589 {
2590  return left + Tensor<rank_, dim, OtherNumber>(right);
2591 }
2592 
2593 
2601 template <int rank_, int dim, typename Number, typename OtherNumber>
2602 constexpr DEAL_II_ALWAYS_INLINE
2603  Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2604  operator-(const SymmetricTensor<rank_, dim, Number> &left,
2605  const Tensor<rank_, dim, OtherNumber> & right)
2606 {
2607  return Tensor<rank_, dim, Number>(left) - right;
2608 }
2609 
2610 
2618 template <int rank_, int dim, typename Number, typename OtherNumber>
2619 constexpr DEAL_II_ALWAYS_INLINE
2620  Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2621  operator-(const Tensor<rank_, dim, Number> & left,
2622  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2623 {
2624  return left - Tensor<rank_, dim, OtherNumber>(right);
2625 }
2626 
2627 
2628 
2643 template <int dim, typename Number>
2644 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2645  determinant(const SymmetricTensor<2, dim, Number> &t)
2646 {
2647  switch (dim)
2648  {
2649  case 1:
2650  return t.data[0];
2651  case 2:
2652  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2653  case 3:
2654  {
2655  // in analogy to general tensors, but
2656  // there's something to be simplified for
2657  // the present case
2658  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2659  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2660  t.data[0] * t.data[5] * t.data[5] -
2661  t.data[1] * t.data[4] * t.data[4] -
2662  t.data[2] * t.data[3] * t.data[3]);
2663  }
2664  default:
2665  Assert(false, ExcNotImplemented());
2666  return internal::NumberType<Number>::value(0.0);
2667  }
2668 }
2669 
2670 
2671 
2684 template <int dim, typename Number>
2685 constexpr DEAL_II_ALWAYS_INLINE Number
2686  third_invariant(const SymmetricTensor<2, dim, Number> &t)
2687 {
2688  return determinant(t);
2689 }
2690 
2691 
2692 
2703 template <int dim, typename Number>
2704 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2705  trace(const SymmetricTensor<2, dim, Number> &d)
2706 {
2707  Number t = d.data[0];
2708  for (unsigned int i = 1; i < dim; ++i)
2709  t += d.data[i];
2710  return t;
2711 }
2712 
2713 
2726 template <int dim, typename Number>
2727 constexpr Number
2728 first_invariant(const SymmetricTensor<2, dim, Number> &t)
2729 {
2730  return trace(t);
2731 }
2732 
2733 
2746 template <typename Number>
2747 constexpr DEAL_II_ALWAYS_INLINE Number
2748  second_invariant(const SymmetricTensor<2, 1, Number> &)
2749 {
2750  return internal::NumberType<Number>::value(0.0);
2751 }
2752 
2753 
2754 
2774 template <typename Number>
2775 constexpr DEAL_II_ALWAYS_INLINE Number
2776  second_invariant(const SymmetricTensor<2, 2, Number> &t)
2777 {
2778  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2779 }
2780 
2781 
2782 
2792 template <typename Number>
2793 constexpr DEAL_II_ALWAYS_INLINE Number
2794  second_invariant(const SymmetricTensor<2, 3, Number> &t)
2795 {
2796  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2797  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2798 }
2799 
2800 
2801 
2810 template <typename Number>
2811 std::array<Number, 1>
2812 eigenvalues(const SymmetricTensor<2, 1, Number> &T);
2813 
2814 
2815 
2839 template <typename Number>
2840 std::array<Number, 2>
2841 eigenvalues(const SymmetricTensor<2, 2, Number> &T);
2842 
2843 
2844 
2868 template <typename Number>
2869 std::array<Number, 3>
2870 eigenvalues(const SymmetricTensor<2, 3, Number> &T);
2871 
2872 
2873 
2874 namespace internal
2875 {
2876  namespace SymmetricTensorImplementation
2877  {
2917  template <int dim, typename Number>
2918  void
2919  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2920  ::Tensor<2, dim, Number> & Q,
2921  std::array<Number, dim> & d,
2922  std::array<Number, dim - 1> & e);
2923 
2924 
2925 
2967  template <int dim, typename Number>
2968  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2969  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2970 
2971 
2972 
3014  template <int dim, typename Number>
3015  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3016  jacobi(::SymmetricTensor<2, dim, Number> A);
3017 
3018 
3019 
3035  template <typename Number>
3036  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3037  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3038 
3039 
3040 
3075  template <typename Number>
3076  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3077  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3078 
3083  template <int dim, typename Number>
3084  struct SortEigenValuesVectors
3085  {
3086  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3087  bool
3088  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3089  {
3090  return lhs.first > rhs.first;
3091  }
3092  };
3093 
3094  } // namespace SymmetricTensorImplementation
3095 
3096 } // namespace internal
3097 
3098 
3099 
3100 // The line below is to ensure that doxygen puts the full description
3101 // of this global enumeration into the documentation
3102 // See https://stackoverflow.com/a/1717984
3131 enum struct SymmetricTensorEigenvectorMethod
3132 {
3142  hybrid,
3152  ql_implicit_shifts,
3160  jacobi
3161 };
3162 
3163 
3164 
3194 template <int dim, typename Number>
3195 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3196  std::integral_constant<int, dim>::value>
3197 eigenvectors(const SymmetricTensor<2, dim, Number> &T,
3198  const SymmetricTensorEigenvectorMethod method =
3199  SymmetricTensorEigenvectorMethod::ql_implicit_shifts);
3200 
3201 
3202 
3212 template <int rank_, int dim, typename Number>
3213 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
3214  transpose(const SymmetricTensor<rank_, dim, Number> &t)
3215 {
3216  return t;
3217 }
3218 
3219 
3220 
3232 template <int dim, typename Number>
3233 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3234  deviator(const SymmetricTensor<2, dim, Number> &t)
3235 {
3236  SymmetricTensor<2, dim, Number> tmp = t;
3237 
3238  // subtract scaled trace from the diagonal
3239  const Number tr = trace(t) / dim;
3240  for (unsigned int i = 0; i < dim; ++i)
3241  tmp.data[i] -= tr;
3242 
3243  return tmp;
3244 }
3245 
3246 
3247 
3255 template <int dim, typename Number>
3256 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3257  unit_symmetric_tensor()
3258 {
3259  // create a default constructed matrix filled with
3260  // zeros, then set the diagonal elements to one
3261  SymmetricTensor<2, dim, Number> tmp;
3262  switch (dim)
3263  {
3264  case 1:
3265  tmp.data[0] = internal::NumberType<Number>::value(1.);
3266  break;
3267  case 2:
3268  tmp.data[0] = tmp.data[1] = internal::NumberType<Number>::value(1.);
3269  break;
3270  case 3:
3271  tmp.data[0] = tmp.data[1] = tmp.data[2] =
3272  internal::NumberType<Number>::value(1.);
3273  break;
3274  default:
3275  for (unsigned int d = 0; d < dim; ++d)
3276  tmp.data[d] = internal::NumberType<Number>::value(1.);
3277  }
3278  return tmp;
3279 }
3280 
3281 
3282 
3291 template <int dim>
3292 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim>
3293  unit_symmetric_tensor()
3294 {
3295  return unit_symmetric_tensor<dim, double>();
3296 }
3297 
3298 
3299 
3329 template <int dim, typename Number>
3330 DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
3331 deviator_tensor()
3332 {
3333  SymmetricTensor<4, dim, Number> tmp;
3334 
3335  // fill the elements treating the diagonal
3336  for (unsigned int i = 0; i < dim; ++i)
3337  for (unsigned int j = 0; j < dim; ++j)
3338  tmp.data[i][j] =
3339  internal::NumberType<Number>::value((i == j ? 1. : 0.) - 1. / dim);
3340 
3341  // then fill the ones that copy over the
3342  // non-diagonal elements. note that during
3343  // the double-contraction, we handle the
3344  // off-diagonal elements twice, so simply
3345  // copying requires a weight of 1/2
3346  for (unsigned int i = dim;
3347  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3348  n_rank2_components;
3349  ++i)
3350  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3351 
3352  return tmp;
3353 }
3354 
3355 
3356 
3365 template <int dim>
3366 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
3367  deviator_tensor()
3368 {
3369  return deviator_tensor<dim, double>();
3370 }
3371 
3372 
3373 
3412 template <int dim, typename Number>
3413 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
3414  identity_tensor()
3415 {
3416  SymmetricTensor<4, dim, Number> tmp;
3417 
3418  // fill the elements treating the diagonal
3419  for (unsigned int i = 0; i < dim; ++i)
3420  tmp.data[i][i] = internal::NumberType<Number>::value(1.);
3421 
3422  // then fill the ones that copy over the
3423  // non-diagonal elements. note that during
3424  // the double-contraction, we handle the
3425  // off-diagonal elements twice, so simply
3426  // copying requires a weight of 1/2
3427  for (unsigned int i = dim;
3428  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3429  n_rank2_components;
3430  ++i)
3431  tmp.data[i][i] = internal::NumberType<Number>::value(0.5);
3432 
3433  return tmp;
3434 }
3435 
3436 
3437 
3446 template <int dim>
3447 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
3448  identity_tensor()
3449 {
3450  return identity_tensor<dim, double>();
3451 }
3452 
3453 
3454 
3465 template <int dim, typename Number>
3466 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3467  invert(const SymmetricTensor<2, dim, Number> &t)
3468 {
3469  return internal::SymmetricTensorImplementation::Inverse<2, dim, Number>::
3470  value(t);
3471 }
3472 
3473 
3474 
3486 template <int dim, typename Number>
3487 constexpr SymmetricTensor<4, dim, Number>
3488 invert(const SymmetricTensor<4, dim, Number> &t)
3489 {
3490  return internal::SymmetricTensorImplementation::Inverse<4, dim, Number>::
3491  value(t);
3492 }
3493 
3494 
3495 
3518 template <int dim, typename Number>
3519 DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
3520 outer_product(const SymmetricTensor<2, dim, Number> &t1,
3521  const SymmetricTensor<2, dim, Number> &t2)
3522 {
3523  SymmetricTensor<4, dim, Number> tmp;
3524 
3525  // fill only the elements really needed
3526  for (unsigned int i = 0; i < dim; ++i)
3527  for (unsigned int j = i; j < dim; ++j)
3528  for (unsigned int k = 0; k < dim; ++k)
3529  for (unsigned int l = k; l < dim; ++l)
3530  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3531 
3532  return tmp;
3533 }
3534 
3535 
3536 
3545 template <int dim, typename Number>
3546 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3547  symmetrize(const Tensor<2, dim, Number> &t)
3548 {
3549  SymmetricTensor<2, dim, Number> result;
3550  for (unsigned int d = 0; d < dim; ++d)
3551  result[d][d] = t[d][d];
3552 
3553  const Number half = internal::NumberType<Number>::value(0.5);
3554  for (unsigned int d = 0; d < dim; ++d)
3555  for (unsigned int e = d + 1; e < dim; ++e)
3556  result[d][e] = (t[d][e] + t[e][d]) * half;
3557  return result;
3558 }
3559 
3560 
3561 
3569 template <int rank_, int dim, typename Number>
3570 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
3571  SymmetricTensor<rank_, dim, Number>
3572  operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3573 {
3574  SymmetricTensor<rank_, dim, Number> tt = t;
3575  tt *= factor;
3576  return tt;
3577 }
3578 
3579 
3580 
3588 template <int rank_, int dim, typename Number>
3589 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
3590  operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3591 {
3592  // simply forward to the other operator
3593  return t * factor;
3594 }
3595 
3596 
3597 
3623 template <int rank_, int dim, typename Number, typename OtherNumber>
3624 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3625  rank_,
3626  dim,
3627  typename ProductType<Number,
3628  typename EnableIfScalar<OtherNumber>::type>::type>
3629 operator*(const SymmetricTensor<rank_, dim, Number> &t,
3630  const OtherNumber & factor)
3631 {
3632  // form the product. we have to convert the two factors into the final
3633  // type via explicit casts because, for awkward reasons, the C++
3634  // standard committee saw it fit to not define an
3635  // operator*(float,std::complex<double>)
3636  // (as well as with switched arguments and double<->float).
3637  using product_type = typename ProductType<Number, OtherNumber>::type;
3638  SymmetricTensor<rank_, dim, product_type> tt(t);
3639  tt *= internal::NumberType<product_type>::value(factor);
3640  return tt;
3641 }
3642 
3643 
3644 
3653 template <int rank_, int dim, typename Number, typename OtherNumber>
3654 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3655  rank_,
3656  dim,
3657  typename ProductType<OtherNumber,
3658  typename EnableIfScalar<Number>::type>::type>
3659 operator*(const Number & factor,
3660  const SymmetricTensor<rank_, dim, OtherNumber> &t)
3661 {
3662  // simply forward to the other operator with switched arguments
3663  return (t * factor);
3664 }
3665 
3666 
3667 
3673 template <int rank_, int dim, typename Number, typename OtherNumber>
3674 DEAL_II_CONSTEXPR inline SymmetricTensor<
3675  rank_,
3676  dim,
3677  typename ProductType<Number,
3678  typename EnableIfScalar<OtherNumber>::type>::type>
3679 operator/(const SymmetricTensor<rank_, dim, Number> &t,
3680  const OtherNumber & factor)
3681 {
3682  using product_type = typename ProductType<Number, OtherNumber>::type;
3683  SymmetricTensor<rank_, dim, product_type> tt(t);
3684  tt /= internal::NumberType<product_type>::value(factor);
3685  return tt;
3686 }
3687 
3688 
3689 
3696 template <int rank_, int dim>
3697 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim>
3698  operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3699 {
3700  SymmetricTensor<rank_, dim> tt(t);
3701  tt *= factor;
3702  return tt;
3703 }
3704 
3705 
3706 
3713 template <int rank_, int dim>
3714 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim>
3715  operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3716 {
3717  SymmetricTensor<rank_, dim> tt(t);
3718  tt *= factor;
3719  return tt;
3720 }
3721 
3722 
3723 
3729 template <int rank_, int dim>
3730 DEAL_II_CONSTEXPR inline SymmetricTensor<rank_, dim>
3731 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3732 {
3733  SymmetricTensor<rank_, dim> tt(t);
3734  tt /= factor;
3735  return tt;
3736 }
3737 
3747 template <int dim, typename Number, typename OtherNumber>
3748 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
3749 scalar_product(const SymmetricTensor<2, dim, Number> & t1,
3750  const SymmetricTensor<2, dim, OtherNumber> &t2)
3751 {
3752  return (t1 * t2);
3753 }
3754 
3755 
3769 template <int dim, typename Number, typename OtherNumber>
3770 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
3771  typename ProductType<Number, OtherNumber>::type
3772  scalar_product(const SymmetricTensor<2, dim, Number> &t1,
3773  const Tensor<2, dim, OtherNumber> & t2)
3774 {
3775  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<
3776  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3777  for (unsigned int i = 0; i < dim; ++i)
3778  for (unsigned int j = 0; j < dim; ++j)
3779  s += t1[i][j] * t2[i][j];
3780  return s;
3781 }
3782 
3783 
3797 template <int dim, typename Number, typename OtherNumber>
3798 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
3799 scalar_product(const Tensor<2, dim, Number> & t1,
3800  const SymmetricTensor<2, dim, OtherNumber> &t2)
3801 {
3802  return scalar_product(t2, t1);
3803 }
3804 
3805 
3821 template <typename Number, typename OtherNumber>
3822 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void double_contract(
3823  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3824  const SymmetricTensor<4, 1, Number> & t,
3825  const SymmetricTensor<2, 1, OtherNumber> & s)
3826 {
3827  tmp[0][0] = t[0][0][0][0] * s[0][0];
3828 }
3829 
3830 
3831 
3847 template <typename Number, typename OtherNumber>
3848 DEAL_II_CONSTEXPR inline void double_contract(
3849  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3850  const SymmetricTensor<2, 1, Number> & s,
3851  const SymmetricTensor<4, 1, OtherNumber> & t)
3852 {
3853  tmp[0][0] = t[0][0][0][0] * s[0][0];
3854 }
3855 
3856 
3857 
3873 template <typename Number, typename OtherNumber>
3874 DEAL_II_CONSTEXPR inline void double_contract(
3875  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3876  const SymmetricTensor<4, 2, Number> & t,
3877  const SymmetricTensor<2, 2, OtherNumber> & s)
3878 {
3879  const unsigned int dim = 2;
3880 
3881  for (unsigned int i = 0; i < dim; ++i)
3882  for (unsigned int j = i; j < dim; ++j)
3883  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3884  2 * t[i][j][0][1] * s[0][1];
3885 }
3886 
3887 
3888 
3904 template <typename Number, typename OtherNumber>
3905 DEAL_II_CONSTEXPR inline void double_contract(
3906  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3907  const SymmetricTensor<2, 2, Number> & s,
3908  const SymmetricTensor<4, 2, OtherNumber> & t)
3909 {
3910  const unsigned int dim = 2;
3911 
3912  for (unsigned int i = 0; i < dim; ++i)
3913  for (unsigned int j = i; j < dim; ++j)
3914  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3915  2 * s[0][1] * t[0][1][i][j];
3916 }
3917 
3918 
3919 
3935 template <typename Number, typename OtherNumber>
3936 DEAL_II_CONSTEXPR inline void double_contract(
3937  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3938  const SymmetricTensor<4, 3, Number> & t,
3939  const SymmetricTensor<2, 3, OtherNumber> & s)
3940 {
3941  const unsigned int dim = 3;
3942 
3943  for (unsigned int i = 0; i < dim; ++i)
3944  for (unsigned int j = i; j < dim; ++j)
3945  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3946  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3947  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3948 }
3949 
3950 
3951 
3967 template <typename Number, typename OtherNumber>
3968 DEAL_II_CONSTEXPR inline void double_contract(
3969  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3970  const SymmetricTensor<2, 3, Number> & s,
3971  const SymmetricTensor<4, 3, OtherNumber> & t)
3972 {
3973  const unsigned int dim = 3;
3974 
3975  for (unsigned int i = 0; i < dim; ++i)
3976  for (unsigned int j = i; j < dim; ++j)
3977  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3978  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3979  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3980 }
3981 
3982 
3983 
3991 template <int dim, typename Number, typename OtherNumber>
3992 DEAL_II_CONSTEXPR
3993  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
3994  operator*(const SymmetricTensor<2, dim, Number> &src1,
3995  const Tensor<1, dim, OtherNumber> & src2)
3996 {
3997  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type> dest;
3998  for (unsigned int i = 0; i < dim; ++i)
3999  for (unsigned int j = 0; j < dim; ++j)
4000  dest[i] += src1[i][j] * src2[j];
4001  return dest;
4002 }
4003 
4004 
4012 template <int dim, typename Number, typename OtherNumber>
4013 constexpr Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
4014 operator*(const Tensor<1, dim, Number> & src1,
4015  const SymmetricTensor<2, dim, OtherNumber> &src2)
4016 {
4017  // this is easy for symmetric tensors:
4018  return src2 * src1;
4019 }
4020 
4021 
4022 
4043 template <int rank_1,
4044  int rank_2,
4045  int dim,
4046  typename Number,
4047  typename OtherNumber>
4048 constexpr DEAL_II_ALWAYS_INLINE
4049  typename Tensor<rank_1 + rank_2 - 2,
4050  dim,
4051  typename ProductType<Number, OtherNumber>::type>::tensor_type
4052  operator*(const Tensor<rank_1, dim, Number> & src1,
4053  const SymmetricTensor<rank_2, dim, OtherNumber> &src2)
4054 {
4055  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
4056 }
4057 
4058 
4059 
4080 template <int rank_1,
4081  int rank_2,
4082  int dim,
4083  typename Number,
4084  typename OtherNumber>
4085 constexpr DEAL_II_ALWAYS_INLINE
4086  typename Tensor<rank_1 + rank_2 - 2,
4087  dim,
4088  typename ProductType<Number, OtherNumber>::type>::tensor_type
4089  operator*(const SymmetricTensor<rank_1, dim, Number> &src1,
4090  const Tensor<rank_2, dim, OtherNumber> & src2)
4091 {
4092  return Tensor<rank_2, dim, OtherNumber>(src1) * src2;
4093 }
4094 
4095 
4096 
4106 template <int dim, typename Number>
4107 inline std::ostream &
4108 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4109 {
4110  // make our lives a bit simpler by outputting
4111  // the tensor through the operator for the
4112  // general Tensor class
4113  Tensor<2, dim, Number> tt;
4114 
4115  for (unsigned int i = 0; i < dim; ++i)
4116  for (unsigned int j = 0; j < dim; ++j)
4117  tt[i][j] = t[i][j];
4118 
4119  return out << tt;
4120 }
4121 
4122 
4123 
4133 template <int dim, typename Number>
4134 inline std::ostream &
4135 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4136 {
4137  // make our lives a bit simpler by outputting
4138  // the tensor through the operator for the
4139  // general Tensor class
4140  Tensor<4, dim, Number> tt;
4141 
4142  for (unsigned int i = 0; i < dim; ++i)
4143  for (unsigned int j = 0; j < dim; ++j)
4144  for (unsigned int k = 0; k < dim; ++k)
4145  for (unsigned int l = 0; l < dim; ++l)
4146  tt[i][j][k][l] = t[i][j][k][l];
4147 
4148  return out << tt;
4149 }
4150 
4151 
4152 DEAL_II_NAMESPACE_CLOSE
4153 
4154 #endif
SymmetricTensor::trace
constexpr Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2705
SymmetricTensor::access_raw_entry
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
operator!=
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
Definition: aligned_vector.h:1192
SymmetricTensor::data
base_tensor_type data
Definition: symmetric_tensor.h:961
TableIndices
Definition: table_indices.h:45
SymmetricTensor
Definition: symmetric_tensor.h:611
internal::SymmetricTensorAccessors::Accessor::operator[]
constexpr Accessor< rank, dim, constness, P - 1, Number > operator[](const unsigned int i)
internal::SymmetricTensorAccessors::merge
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
Definition: symmetric_tensor.h:146
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
operator==
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
Definition: aligned_vector.h:1170
determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2645
SymmetricTensor::operator+
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Definition: symmetric_tensor.h:2525
TransposeTableIterators::Accessor
MatrixTableIterators::Accessor< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Accessor
Definition: table.h:1907
SymmetricTensor::operator==
constexpr bool operator==(const SymmetricTensor &) const
internal::NumberType::value
static constexpr const T & value(const T &t)
Definition: numbers.h:703
SymmetricTensor::scalar_product
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Definition: symmetric_tensor.h:3749
SymmetricTensor::outer_product
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition: symmetric_tensor.h:3520
SymmetricTensor::deviator
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3234
internal::SymmetricTensorAccessors::Accessor::previous_indices
const TableIndices< rank > previous_indices
Definition: symmetric_tensor.h:434
internal::SymmetricTensorImplementation::SortEigenValuesVectors::EigValsVecs
std::pair< Number, Tensor< 1, dim, Number > > EigValsVecs
Definition: symmetric_tensor.h:3086
SymmetricTensor::third_invariant
constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2686
SymmetricTensor::eigenvectors
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)
SymmetricTensorEigenvectorMethod::jacobi
@ jacobi
SymmetricTensor::operator-
constexpr SymmetricTensor operator-() const
SymmetricTensor::norm
constexpr numbers::NumberTraits< Number >::real_type norm() const
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
ProductType::type
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
Definition: template_constraints.h:426
internal::SymmetricTensorAccessors::Accessor< rank, dim, constness, 1, Number >::previous_indices
const TableIndices< rank > previous_indices
Definition: symmetric_tensor.h:518
SymmetricTensor::invert
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3467
internal::SymmetricTensorAccessors::AccessorTypes
Definition: symmetric_tensor.h:311
deviator
constexpr SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3234
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor::unrolled_to_component_indices
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
internal::SymmetricTensorAccessors::Accessor< rank, dim, constness, 1, Number >::tensor_type
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
Definition: symmetric_tensor.h:469
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:99
internal::SymmetricTensorImplementation
Definition: symmetric_tensor.h:123
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
DEAL_II_CONSTEXPR
#define DEAL_II_CONSTEXPR
Definition: config.h:102
internal::SymmetricTensorAccessors::double_contraction_result::type
::SymmetricTensor< rank1+rank2 - 4, dim, value_type > type
Definition: symmetric_tensor.h:219
ProductType
Definition: template_constraints.h:422
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition: symengine_math.cc:273
identity_tensor
constexpr SymmetricTensor< 4, dim, Number > identity_tensor()
Definition: symmetric_tensor.h:3414
LAPACKSupport::T
static const char T
Definition: lapack_support.h:163
tensor.h
SymmetricTensor::operator*=
constexpr SymmetricTensor & operator*=(const OtherNumber &factor)
Point::operator*
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
unit_symmetric_tensor
constexpr SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: symmetric_tensor.h:3257
SymmetricTensor::operator/=
constexpr SymmetricTensor & operator/=(const OtherNumber &factor)
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
internal::SymmetricTensorImplementation::Inverse
Definition: symmetric_tensor.h:130
SymmetricTensorEigenvectorMethod::ql_implicit_shifts
@ ql_implicit_shifts
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor::serialize
void serialize(Archive &ar, const unsigned int version)
internal::SymmetricTensorAccessors::Accessor::Accessor
friend class Accessor
Definition: symmetric_tensor.h:441
internal::SymmetricTensorAccessors::Accessor::tensor_type
typename AccessorTypes< rank, dim, constness, Number >::tensor_type tensor_type
Definition: symmetric_tensor.h:386
SymmetricTensor::operator!=
constexpr bool operator!=(const SymmetricTensor &) const
Tensor
Definition: tensor.h:450
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
trace
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:268
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor::first_invariant
constexpr Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2728
internal::SymmetricTensorImplementation::SortEigenValuesVectors
Definition: symmetric_tensor.h:3084
internal::SymmetricTensorAccessors::AccessorTypes< rank, dim, true, Number >::reference
Number reference
Definition: symmetric_tensor.h:324
table_indices.h
SymmetricTensor::operator*
constexpr internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
internal::SymmetricTensorAccessors::Accessor< rank, dim, constness, 1, Number >::tensor
tensor_type & tensor
Definition: symmetric_tensor.h:517
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
numbers::NumberTraits::real_type
number real_type
Definition: numbers.h:437
MemorySpace::swap
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
Definition: memory_space.h:103
SymmetricTensorEigenvectorMethod
SymmetricTensorEigenvectorMethod
Definition: symmetric_tensor.h:3131
operator-
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Definition: symmetric_tensor.h:2550
SymmetricTensor::operator-
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Definition: symmetric_tensor.h:2550
internal::SymmetricTensorAccessors::AccessorTypes< rank, dim, false, Number >::reference
Number & reference
Definition: symmetric_tensor.h:338
SymmetricTensor::operator-=
constexpr SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
internal::SymmetricTensorAccessors::double_contraction_result< 2, 2, dim, Number, OtherNumber >::type
typename ProductType< Number, OtherNumber >::type type
Definition: symmetric_tensor.h:234
SymmetricTensor::clear
constexpr void clear()
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
operator*
constexpr SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
Definition: symmetric_tensor.h:3572
SymmetricTensor::SymmetricTensor
friend class SymmetricTensor
Definition: symmetric_tensor.h:965
value
static const bool value
Definition: dof_tools_constraints.cc:433
SymmetricTensor::end_raw
Number * end_raw()
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
internal::SymmetricTensorImplementation::SortEigenValuesVectors::operator()
bool operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
Definition: symmetric_tensor.h:3088
std::sqrt
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
SymmetricTensor::symmetrize
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3547
internal::SymmetricTensorAccessors::double_contraction_result
Definition: symmetric_tensor.h:215
SymmetricTensor::determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2645
Algorithms::OutputOperator::operator<<
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:173
internal::NumberType
Definition: numbers.h:700
internal::SymmetricTensorAccessors::Accessor::reference
typename AccessorTypes< rank, dim, constness, Number >::reference reference
Definition: symmetric_tensor.h:384
ProductType::operator/
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator/(const std::complex< T > &left, const std::complex< U > &right)
Definition: complex_overloads.h:61
numbers::NumberTraits::abs
static real_type abs(const number &x)
Definition: numbers.h:609
SymmetricTensor< 2, spacedim >::base_tensor_type
typename base_tensor_descriptor::base_tensor_type base_tensor_type
Definition: symmetric_tensor.h:956
deviator_tensor
constexpr SymmetricTensor< 4, dim, Number > deviator_tensor()
Definition: symmetric_tensor.h:3331
internal::SymmetricTensorAccessors::double_contraction_result::value_type
typename ProductType< Number, OtherNumber >::type value_type
Definition: symmetric_tensor.h:217
SymmetricTensor::double_contract
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)
Definition: symmetric_tensor.h:3822
internal::SymmetricTensorAccessors::Accessor< rank, dim, constness, 1, Number >::reference
typename AccessorTypes< rank, dim, constness, Number >::reference reference
Definition: symmetric_tensor.h:467
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
SymmetricTensor::eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
template_constraints.h
config.h
EnableIfScalar
Definition: template_constraints.h:534
internal::SymmetricTensorAccessors::StorageType
Definition: symmetric_tensor.h:252
internal
Definition: aligned_vector.h:369
LAPACKSupport::N
static const char N
Definition: lapack_support.h:159
SymmetricTensor::operator+=
constexpr SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
SymmetricTensorEigenvectorMethod::hybrid
@ hybrid
value_type
SymmetricTensor::operator[]
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
internal::SymmetricTensorAccessors::Accessor
Definition: symmetric_tensor.h:377
SymmetricTensor::operator=
constexpr SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
internal::ProductTypeImpl
Definition: template_constraints.h:364
SymmetricTensor::second_invariant
constexpr Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Definition: symmetric_tensor.h:2748
numbers.h
numbers::value_is_zero
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:941
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
SymmetricTensor::begin_raw
Number * begin_raw()
SymmetricTensor::operator()
constexpr Number & operator()(const TableIndices< rank_ > &indices)
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
int
internal::SymmetricTensorImplementation::tridiagonalize
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A, ::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim - 1 > &e)
internal::SymmetricTensorAccessors::Accessor::tensor
tensor_type & tensor
Definition: symmetric_tensor.h:433
numbers::NumberTraits
Definition: numbers.h:422
SymmetricTensor::component_to_unrolled_index
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensor::memory_consumption
static constexpr std::size_t memory_consumption()
Generated by   doxygen 1.8.17