Reference documentation for deal.II version Git c679197b1a 2019-11-19 16:09:58 -0500
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
symmetric_tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_symmetric_tensor_h
17 #define dealii_symmetric_tensor_h
18 
19 
20 #include <deal.II/base/numbers.h>
21 #include <deal.II/base/table_indices.h>
22 #include <deal.II/base/template_constraints.h>
23 #include <deal.II/base/tensor.h>
24 
25 #include <algorithm>
26 #include <array>
27 #include <functional>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 // Forward declaration
32 #ifndef DOXYGEN
33 template <int rank, int dim, typename Number = double>
34 class SymmetricTensor;
35 #endif
36 
37 template <int dim, typename Number>
38 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
40 
41 template <int dim, typename Number>
42 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
44 
45 template <int dim, typename Number>
46 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
48 
49 template <int dim, typename Number>
50 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
52 
53 template <int dim, typename Number>
54 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
56 
57 template <int dim2, typename Number>
58 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
60 
61 template <int dim, typename Number>
62 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
64 
65 template <int dim, typename Number>
66 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
68 
69 
70 
71 namespace internal
72 {
77  namespace SymmetricTensorImplementation
78  {
83  template <int rank, int dim, typename Number>
84  struct Inverse;
85  } // namespace SymmetricTensorImplementation
86 
91  namespace SymmetricTensorAccessors
92  {
99  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
100  merge(const TableIndices<2> &previous_indices,
101  const unsigned int new_index,
102  const unsigned int position)
103  {
104  Assert(position < 2, ExcIndexRange(position, 0, 2));
105 
106  if (position == 0)
107  return {new_index, numbers::invalid_unsigned_int};
108  else
109  return {previous_indices[0], new_index};
110  }
111 
112 
113 
120  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<4>
121  merge(const TableIndices<4> &previous_indices,
122  const unsigned int new_index,
123  const unsigned int position)
124  {
125  Assert(position < 4, ExcIndexRange(position, 0, 4));
126 
127  switch (position)
128  {
129  case 0:
130  return {new_index,
133  numbers::invalid_unsigned_int};
134  case 1:
135  return {previous_indices[0],
136  new_index,
138  numbers::invalid_unsigned_int};
139  case 2:
140  return {previous_indices[0],
141  previous_indices[1],
142  new_index,
143  numbers::invalid_unsigned_int};
144  case 3:
145  return {previous_indices[0],
146  previous_indices[1],
147  previous_indices[2],
148  new_index};
149  default:
150  Assert(false, ExcInternalError());
151  return {};
152  }
153  }
154 
155 
164  template <int rank1,
165  int rank2,
166  int dim,
167  typename Number,
168  typename OtherNumber = Number>
170  {
171  using value_type = typename ProductType<Number, OtherNumber>::type;
172  using type =
173  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
174  };
175 
176 
185  template <int dim, typename Number, typename OtherNumber>
186  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
187  {
188  using type = typename ProductType<Number, OtherNumber>::type;
189  };
190 
191 
192 
205  template <int rank, int dim, typename Number>
206  struct StorageType;
207 
211  template <int dim, typename Number>
212  struct StorageType<2, dim, Number>
213  {
218  static const unsigned int n_independent_components =
219  (dim * dim + dim) / 2;
220 
225  };
226 
227 
228 
232  template <int dim, typename Number>
233  struct StorageType<4, dim, Number>
234  {
240  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
241 
245  static const unsigned int n_independent_components =
246  (n_rank2_components *
248 
256  };
257 
258 
259 
264  template <int rank, int dim, bool constness, typename Number>
266 
273  template <int rank, int dim, typename Number>
274  struct AccessorTypes<rank, dim, true, Number>
275  {
276  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
277 
278  using reference = Number;
279  };
280 
287  template <int rank, int dim, typename Number>
288  struct AccessorTypes<rank, dim, false, Number>
289  {
291 
292  using reference = Number &;
293  };
294 
295 
330  template <int rank, int dim, bool constness, int P, typename Number>
331  class Accessor
332  {
333  public:
337  using reference =
339  using tensor_type =
341 
342  private:
361  constexpr Accessor(tensor_type & tensor,
362  const TableIndices<rank> &previous_indices);
363 
367  constexpr DEAL_II_ALWAYS_INLINE
368  Accessor(const Accessor &) = default;
369 
370  public:
374  DEAL_II_CONSTEXPR Accessor<rank, dim, constness, P - 1, Number>
375  operator[](const unsigned int i);
376 
380  constexpr Accessor<rank, dim, constness, P - 1, Number>
381  operator[](const unsigned int i) const;
382 
383  private:
387  tensor_type & tensor;
388  const TableIndices<rank> previous_indices;
389 
390  // Declare some other classes as friends. Make sure to work around bugs
391  // in some compilers:
392  template <int, int, typename>
393  friend class ::SymmetricTensor;
394  template <int, int, bool, int, typename>
395  friend class Accessor;
396 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
397  friend class ::SymmetricTensor<rank, dim, Number>;
398  friend class Accessor<rank, dim, constness, P + 1, Number>;
399 #endif
400  };
401 
402 
403 
413  template <int rank, int dim, bool constness, typename Number>
414  class Accessor<rank, dim, constness, 1, Number>
415  {
416  public:
420  using reference =
422  using tensor_type =
424 
425  private:
447  constexpr Accessor(tensor_type & tensor,
448  const TableIndices<rank> &previous_indices);
449 
453  constexpr DEAL_II_ALWAYS_INLINE
454  Accessor(const Accessor &) = default;
455 
456  public:
460  DEAL_II_CONSTEXPR reference operator[](const unsigned int);
461 
465  constexpr reference operator[](const unsigned int) const;
466 
467  private:
471  tensor_type & tensor;
472  const TableIndices<rank> previous_indices;
473 
474  // Declare some other classes as friends. Make sure to work around bugs
475  // in some compilers:
476  template <int, int, typename>
477  friend class ::SymmetricTensor;
478  template <int, int, bool, int, typename>
479  friend class SymmetricTensorAccessors::Accessor;
480 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
481  friend class ::SymmetricTensor<rank, dim, Number>;
482  friend class SymmetricTensorAccessors::
483  Accessor<rank, dim, constness, 2, Number>;
484 #endif
485  };
486  } // namespace SymmetricTensorAccessors
487 } // namespace internal
488 
489 
490 
554 template <int rank_, int dim, typename Number>
556 {
557 public:
558  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
559 
568  static const unsigned int dimension = dim;
569 
573  static const unsigned int rank = rank_;
574 
580  static constexpr unsigned int n_independent_components =
582  n_independent_components;
583 
587  constexpr DEAL_II_ALWAYS_INLINE
588  SymmetricTensor() = default;
589 
603  template <typename OtherNumber>
604  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
605 
621  DEAL_II_CONSTEXPR
622  SymmetricTensor(const Number (&array)[n_independent_components]);
623 
629  template <typename OtherNumber>
630  constexpr explicit SymmetricTensor(
631  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
632 
636  Number *
637  begin_raw();
638 
642  const Number *
643  begin_raw() const;
644 
648  Number *
649  end_raw();
650 
655  const Number *
656  end_raw() const;
657 
664  template <typename OtherNumber>
665  DEAL_II_CONSTEXPR SymmetricTensor &
666  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
667 
674  DEAL_II_CONSTEXPR SymmetricTensor &
675  operator=(const Number &d);
676 
681  constexpr operator Tensor<rank_, dim, Number>() const;
682 
686  constexpr bool
687  operator==(const SymmetricTensor &) const;
688 
692  constexpr bool
693  operator!=(const SymmetricTensor &) const;
694 
698  template <typename OtherNumber>
699  DEAL_II_CONSTEXPR SymmetricTensor &
700  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
701 
705  template <typename OtherNumber>
706  DEAL_II_CONSTEXPR SymmetricTensor &
707  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
708 
713  template <typename OtherNumber>
714  DEAL_II_CONSTEXPR SymmetricTensor &
715  operator*=(const OtherNumber &factor);
716 
720  template <typename OtherNumber>
721  DEAL_II_CONSTEXPR SymmetricTensor &
722  operator/=(const OtherNumber &factor);
723 
727  DEAL_II_CONSTEXPR SymmetricTensor
728  operator-() const;
729 
754  template <typename OtherNumber>
755  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
756  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
758 
763  template <typename OtherNumber>
764  DEAL_II_CONSTEXPR typename internal::SymmetricTensorAccessors::
765  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
767 
771  DEAL_II_CONSTEXPR Number &
772  operator()(const TableIndices<rank_> &indices);
773 
777  DEAL_II_CONSTEXPR const Number &
778  operator()(const TableIndices<rank_> &indices) const;
779 
784  constexpr internal::SymmetricTensorAccessors::
785  Accessor<rank_, dim, true, rank_ - 1, Number>
786  operator[](const unsigned int row) const;
787 
792  DEAL_II_CONSTEXPR internal::SymmetricTensorAccessors::
793  Accessor<rank_, dim, false, rank_ - 1, Number>
794  operator[](const unsigned int row);
795 
801  constexpr const Number &operator[](const TableIndices<rank_> &indices) const;
802 
808  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
809 
815  DEAL_II_CONSTEXPR const Number &
816  access_raw_entry(const unsigned int unrolled_index) const;
817 
823  DEAL_II_CONSTEXPR Number &
824  access_raw_entry(const unsigned int unrolled_index);
825 
836  norm() const;
837 
845  static constexpr unsigned int
846  component_to_unrolled_index(const TableIndices<rank_> &indices);
847 
853  static constexpr TableIndices<rank_>
854  unrolled_to_component_indices(const unsigned int i);
855 
868  DEAL_II_CONSTEXPR void
869  clear();
870 
875  static constexpr std::size_t
876  memory_consumption();
877 
882  template <class Archive>
883  void
884  serialize(Archive &ar, const unsigned int version);
885 
886 private:
890  using base_tensor_descriptor =
892 
896  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
897 
902 
903  // Make all other symmetric tensors friends.
904  template <int, int, typename>
905  friend class SymmetricTensor;
906 
907  // Make a few more functions friends.
908  template <int dim2, typename Number2>
909  friend DEAL_II_CONSTEXPR Number2
911 
912  template <int dim2, typename Number2>
913  friend DEAL_II_CONSTEXPR Number2
915 
916  template <int dim2, typename Number2>
917  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
919 
920  template <int dim2, typename Number2>
921  friend DEAL_II_CONSTEXPR SymmetricTensor<2, dim2, Number2>
923 
924  template <int dim2, typename Number2>
925  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
926  deviator_tensor();
927 
928  template <int dim2, typename Number2>
929  friend DEAL_II_CONSTEXPR SymmetricTensor<4, dim2, Number2>
930  identity_tensor();
931 
932 
933  // Make a few helper classes friends as well.
935  Inverse<2, dim, Number>;
936 
938  Inverse<4, dim, Number>;
939 };
940 
941 
942 
943 // ------------------------- inline functions ------------------------
944 
945 #ifndef DOXYGEN
946 
947 // provide declarations for static members
948 template <int rank, int dim, typename Number>
949 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
950 
951 template <int rank_, int dim, typename Number>
952 constexpr unsigned int
953  SymmetricTensor<rank_, dim, Number>::n_independent_components;
954 
955 namespace internal
956 {
957  namespace SymmetricTensorAccessors
958  {
959  template <int rank_, int dim, bool constness, int P, typename Number>
960  constexpr DEAL_II_ALWAYS_INLINE
961  Accessor<rank_, dim, constness, P, Number>::Accessor(
962  tensor_type & tensor,
963  const TableIndices<rank_> &previous_indices)
964  : tensor(tensor)
965  , previous_indices(previous_indices)
966  {}
967 
968 
969 
970  template <int rank_, int dim, bool constness, int P, typename Number>
971  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
972  Accessor<rank_, dim, constness, P - 1, Number>
973  Accessor<rank_, dim, constness, P, Number>::
974  operator[](const unsigned int i)
975  {
976  return Accessor<rank_, dim, constness, P - 1, Number>(
977  tensor, merge(previous_indices, i, rank_ - P));
978  }
979 
980 
981 
982  template <int rank_, int dim, bool constness, int P, typename Number>
983  constexpr DEAL_II_ALWAYS_INLINE
984  Accessor<rank_, dim, constness, P - 1, Number>
985  Accessor<rank_, dim, constness, P, Number>::
986  operator[](const unsigned int i) const
987  {
988  return Accessor<rank_, dim, constness, P - 1, Number>(
989  tensor, merge(previous_indices, i, rank_ - P));
990  }
991 
992 
993 
994  template <int rank_, int dim, bool constness, typename Number>
995  constexpr DEAL_II_ALWAYS_INLINE
996  Accessor<rank_, dim, constness, 1, Number>::Accessor(
997  tensor_type & tensor,
998  const TableIndices<rank_> &previous_indices)
999  : tensor(tensor)
1000  , previous_indices(previous_indices)
1001  {}
1002 
1003 
1004 
1005  template <int rank_, int dim, bool constness, typename Number>
1006  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1007  typename Accessor<rank_, dim, constness, 1, Number>::reference
1008  Accessor<rank_, dim, constness, 1, Number>::
1009  operator[](const unsigned int i)
1010  {
1011  return tensor(merge(previous_indices, i, rank_ - 1));
1012  }
1013 
1014 
1015  template <int rank_, int dim, bool constness, typename Number>
1016  constexpr DEAL_II_ALWAYS_INLINE
1017  typename Accessor<rank_, dim, constness, 1, Number>::reference
1018  Accessor<rank_, dim, constness, 1, Number>::
1019  operator[](const unsigned int i) const
1020  {
1021  return tensor(merge(previous_indices, i, rank_ - 1));
1022  }
1023  } // namespace SymmetricTensorAccessors
1024 } // namespace internal
1025 
1026 
1027 
1028 template <int rank_, int dim, typename Number>
1029 template <typename OtherNumber>
1030 inline DEAL_II_ALWAYS_INLINE
1032  const Tensor<2, dim, OtherNumber> &t)
1033 {
1034  static_assert(rank == 2, "This function is only implemented for rank==2");
1035  for (unsigned int d = 0; d < dim; ++d)
1036  for (unsigned int e = 0; e < d; ++e)
1037  Assert(t[d][e] == t[e][d], ExcInternalError());
1038 
1039  for (unsigned int d = 0; d < dim; ++d)
1040  data[d] = t[d][d];
1041 
1042  for (unsigned int d = 0, c = 0; d < dim; ++d)
1043  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1044  data[dim + c] = t[d][e];
1045 }
1046 
1047 
1048 
1049 template <int rank_, int dim, typename Number>
1050 template <typename OtherNumber>
1051 constexpr DEAL_II_ALWAYS_INLINE
1053  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1054  : data(initializer.data)
1055 {}
1056 
1057 
1058 
1059 template <int rank_, int dim, typename Number>
1060 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1062  const Number (&array)[n_independent_components])
1063  : data(
1064  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1065 {
1066  // ensure that the reinterpret_cast above actually works
1067  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1068  ExcInternalError());
1069 }
1070 
1071 
1072 
1073 template <int rank_, int dim, typename Number>
1074 template <typename OtherNumber>
1075 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1079 {
1080  data = t.data;
1081  return *this;
1082 }
1083 
1084 
1085 
1086 template <int rank_, int dim, typename Number>
1087 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1090 {
1092  ExcMessage("Only assignment with zero is allowed"));
1093  (void)d;
1094 
1095  data = internal::NumberType<Number>::value(0.0);
1096 
1097  return *this;
1098 }
1099 
1100 
1101 namespace internal
1102 {
1103  namespace SymmetricTensorImplementation
1104  {
1105  template <int dim, typename Number>
1106  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1107  ::Tensor<2, dim, Number>
1108  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1109  {
1111 
1112  // diagonal entries are stored first
1113  for (unsigned int d = 0; d < dim; ++d)
1114  t[d][d] = s.access_raw_entry(d);
1115 
1116  // off-diagonal entries come next, row by row
1117  for (unsigned int d = 0, c = 0; d < dim; ++d)
1118  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1119  {
1120  t[d][e] = s.access_raw_entry(dim + c);
1121  t[e][d] = s.access_raw_entry(dim + c);
1122  }
1123  return t;
1124  }
1125 
1126 
1127  template <int dim, typename Number>
1128  DEAL_II_CONSTEXPR ::Tensor<4, dim, Number>
1129  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1130  {
1131  // utilize the symmetry properties of SymmetricTensor<4,dim>
1132  // discussed in the class documentation to avoid accessing all
1133  // independent elements of the input tensor more than once
1135 
1136  for (unsigned int i = 0; i < dim; ++i)
1137  for (unsigned int j = i; j < dim; ++j)
1138  for (unsigned int k = 0; k < dim; ++k)
1139  for (unsigned int l = k; l < dim; ++l)
1140  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1141  t[TableIndices<4>(j, i, k, l)] =
1142  t[TableIndices<4>(j, i, l, k)] =
1143  st[TableIndices<4>(i, j, k, l)];
1144 
1145  return t;
1146  }
1147 
1148 
1149  template <typename Number>
1150  struct Inverse<2, 1, Number>
1151  {
1152  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1153  ::SymmetricTensor<2, 1, Number>
1154  value(const ::SymmetricTensor<2, 1, Number> &t)
1155  {
1157 
1158  tmp[0][0] = 1.0 / t[0][0];
1159 
1160  return tmp;
1161  }
1162  };
1163 
1164 
1165  template <typename Number>
1166  struct Inverse<2, 2, Number>
1167  {
1168  DEAL_II_CONSTEXPR static inline DEAL_II_ALWAYS_INLINE
1169  ::SymmetricTensor<2, 2, Number>
1170  value(const ::SymmetricTensor<2, 2, Number> &t)
1171  {
1173 
1174  // Sympy result: ([
1175  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1176  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1177  const TableIndices<2> idx_00(0, 0);
1178  const TableIndices<2> idx_01(0, 1);
1179  const TableIndices<2> idx_11(1, 1);
1180  const Number inv_det_t =
1181  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1182  tmp[idx_00] = t[idx_11];
1183  tmp[idx_01] = -t[idx_01];
1184  tmp[idx_11] = t[idx_00];
1185  tmp *= inv_det_t;
1186 
1187  return tmp;
1188  }
1189  };
1190 
1191 
1192  template <typename Number>
1193  struct Inverse<2, 3, Number>
1194  {
1195  DEAL_II_CONSTEXPR static ::SymmetricTensor<2, 3, Number>
1196  value(const ::SymmetricTensor<2, 3, Number> &t)
1197  {
1199 
1200  // Sympy result: ([
1201  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1202  // 2*t01*t02*t12 - t02**2*t11),
1203  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1204  // 2*t01*t02*t12 - t02**2*t11),
1205  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1206  // 2*t01*t02*t12 - t02**2*t11)],
1207  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1208  // 2*t01*t02*t12 - t02**2*t11),
1209  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1210  // 2*t01*t02*t12 - t02**2*t11),
1211  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1212  // 2*t01*t02*t12 + t02**2*t11)],
1213  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1214  // 2*t01*t02*t12 - t02**2*t11),
1215  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1216  // 2*t01*t02*t12 + t02**2*t11),
1217  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1218  // 2*t01*t02*t12 + t02**2*t11)] ])
1219  //
1220  // =
1221  //
1222  // [ (t11*t22 - t12**2)/det_t,
1223  // (-t01*t22 + t02*t12)/det_t,
1224  // (t01*t12 - t02*t11)/det_t],
1225  // [ (-t01*t22 + t02*t12)/det_t,
1226  // (t00*t22 - t02**2)/det_t,
1227  // (-t00*t12 + t01*t02)/det_t],
1228  // [ (t01*t12 - t02*t11)/det_t,
1229  // (-t00*t12 + t01*t02)/det_t,
1230  // (t00*t11 - t01**2)/det_t] ])
1231  //
1232  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1233  // 2*t01*t02*t12 - t02**2*t11)
1234  const TableIndices<2> idx_00(0, 0);
1235  const TableIndices<2> idx_01(0, 1);
1236  const TableIndices<2> idx_02(0, 2);
1237  const TableIndices<2> idx_11(1, 1);
1238  const TableIndices<2> idx_12(1, 2);
1239  const TableIndices<2> idx_22(2, 2);
1240  const Number inv_det_t =
1241  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1242  t[idx_00] * t[idx_12] * t[idx_12] -
1243  t[idx_01] * t[idx_01] * t[idx_22] +
1244  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1245  t[idx_02] * t[idx_02] * t[idx_11]);
1246  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1247  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1248  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1249  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1250  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1251  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1252  tmp *= inv_det_t;
1253 
1254  return tmp;
1255  }
1256  };
1257 
1258 
1259  template <typename Number>
1260  struct Inverse<4, 1, Number>
1261  {
1262  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 1, Number>
1263  value(const ::SymmetricTensor<4, 1, Number> &t)
1264  {
1266  tmp.data[0][0] = 1.0 / t.data[0][0];
1267  return tmp;
1268  }
1269  };
1270 
1271 
1272  template <typename Number>
1273  struct Inverse<4, 2, Number>
1274  {
1275  DEAL_II_CONSTEXPR static inline ::SymmetricTensor<4, 2, Number>
1276  value(const ::SymmetricTensor<4, 2, Number> &t)
1277  {
1279 
1280  // Inverting this tensor is a little more complicated than necessary,
1281  // since we store the data of 't' as a 3x3 matrix t.data, but the
1282  // product between a rank-4 and a rank-2 tensor is really not the
1283  // product between this matrix and the 3-vector of a rhs, but rather
1284  //
1285  // B.vec = t.data * mult * A.vec
1286  //
1287  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1288  // capture the fact that we need to add up both the c_ij12*a_12 and the
1289  // c_ij21*a_21 terms.
1290  //
1291  // In addition, in this scheme, the identity tensor has the matrix
1292  // representation mult^-1.
1293  //
1294  // The inverse of 't' therefore has the matrix representation
1295  //
1296  // inv.data = mult^-1 * t.data^-1 * mult^-1
1297  //
1298  // in order to compute it, let's first compute the inverse of t.data and
1299  // put it into tmp.data; at the end of the function we then scale the
1300  // last row and column of the inverse by 1/2, corresponding to the left
1301  // and right multiplication with mult^-1.
1302  const Number t4 = t.data[0][0] * t.data[1][1],
1303  t6 = t.data[0][0] * t.data[1][2],
1304  t8 = t.data[0][1] * t.data[1][0],
1305  t00 = t.data[0][2] * t.data[1][0],
1306  t01 = t.data[0][1] * t.data[2][0],
1307  t04 = t.data[0][2] * t.data[2][0],
1308  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1309  t8 * t.data[2][2] + t00 * t.data[2][1] +
1310  t01 * t.data[1][2] - t04 * t.data[1][1]);
1311  tmp.data[0][0] =
1312  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1313  tmp.data[0][1] =
1314  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1315  tmp.data[0][2] =
1316  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1317  tmp.data[1][0] =
1318  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1319  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1320  tmp.data[1][2] = -(t6 - t00) * t07;
1321  tmp.data[2][0] =
1322  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1323  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1324  tmp.data[2][2] = (t4 - t8) * t07;
1325 
1326  // scale last row and column as mentioned
1327  // above
1328  tmp.data[2][0] /= 2;
1329  tmp.data[2][1] /= 2;
1330  tmp.data[0][2] /= 2;
1331  tmp.data[1][2] /= 2;
1332  tmp.data[2][2] /= 4;
1333 
1334  return tmp;
1335  }
1336  };
1337 
1338 
1339  template <typename Number>
1340  struct Inverse<4, 3, Number>
1341  {
1342  static ::SymmetricTensor<4, 3, Number>
1343  value(const ::SymmetricTensor<4, 3, Number> &t)
1344  {
1346 
1347  // This function follows the exact same scheme as the 2d case, except
1348  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1349  // Instead, we use the Gauss-Jordan algorithm implemented for
1350  // FullMatrix. For historical reasons the following code is copied from
1351  // there, with the tangential benefit that we do not need to copy the
1352  // tensor entries to and from the FullMatrix.
1353  const unsigned int N = 6;
1354 
1355  // First get an estimate of the size of the elements of this matrix,
1356  // for later checks whether the pivot element is large enough, or
1357  // whether we have to fear that the matrix is not regular.
1358  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1359  for (unsigned int i = 0; i < N; ++i)
1360  diagonal_sum += std::fabs(tmp.data[i][i]);
1361  const Number typical_diagonal_element =
1362  diagonal_sum / static_cast<double>(N);
1363  (void)typical_diagonal_element;
1364 
1365  unsigned int p[N];
1366  for (unsigned int i = 0; i < N; ++i)
1367  p[i] = i;
1368 
1369  for (unsigned int j = 0; j < N; ++j)
1370  {
1371  // Pivot search: search that part of the line on and right of the
1372  // diagonal for the largest element.
1373  Number max = std::fabs(tmp.data[j][j]);
1374  unsigned int r = j;
1375  for (unsigned int i = j + 1; i < N; ++i)
1376  if (std::fabs(tmp.data[i][j]) > max)
1377  {
1378  max = std::fabs(tmp.data[i][j]);
1379  r = i;
1380  }
1381 
1382  // Check whether the pivot is too small
1383  Assert(max > 1.e-16 * typical_diagonal_element,
1384  ExcMessage("This tensor seems to be noninvertible"));
1385 
1386  // Row interchange
1387  if (r > j)
1388  {
1389  for (unsigned int k = 0; k < N; ++k)
1390  std::swap(tmp.data[j][k], tmp.data[r][k]);
1391 
1392  std::swap(p[j], p[r]);
1393  }
1394 
1395  // Transformation
1396  const Number hr = 1. / tmp.data[j][j];
1397  tmp.data[j][j] = hr;
1398  for (unsigned int k = 0; k < N; ++k)
1399  {
1400  if (k == j)
1401  continue;
1402  for (unsigned int i = 0; i < N; ++i)
1403  {
1404  if (i == j)
1405  continue;
1406  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1407  }
1408  }
1409  for (unsigned int i = 0; i < N; ++i)
1410  {
1411  tmp.data[i][j] *= hr;
1412  tmp.data[j][i] *= -hr;
1413  }
1414  tmp.data[j][j] = hr;
1415  }
1416 
1417  // Column interchange
1418  Number hv[N];
1419  for (unsigned int i = 0; i < N; ++i)
1420  {
1421  for (unsigned int k = 0; k < N; ++k)
1422  hv[p[k]] = tmp.data[i][k];
1423  for (unsigned int k = 0; k < N; ++k)
1424  tmp.data[i][k] = hv[k];
1425  }
1426 
1427  // Scale rows and columns. The mult matrix
1428  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1429  for (unsigned int i = 3; i < 6; ++i)
1430  for (unsigned int j = 0; j < 3; ++j)
1431  tmp.data[i][j] /= 2;
1432 
1433  for (unsigned int i = 0; i < 3; ++i)
1434  for (unsigned int j = 3; j < 6; ++j)
1435  tmp.data[i][j] /= 2;
1436 
1437  for (unsigned int i = 3; i < 6; ++i)
1438  for (unsigned int j = 3; j < 6; ++j)
1439  tmp.data[i][j] /= 4;
1440 
1441  return tmp;
1442  }
1443  };
1444 
1445  } // namespace SymmetricTensorImplementation
1446 } // namespace internal
1447 
1448 
1449 
1450 template <int rank_, int dim, typename Number>
1451 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>::
1452  operator Tensor<rank_, dim, Number>() const
1453 {
1454  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1455 }
1456 
1457 
1458 
1459 template <int rank_, int dim, typename Number>
1460 constexpr bool
1463 {
1464  return data == t.data;
1465 }
1466 
1467 
1468 
1469 template <int rank_, int dim, typename Number>
1470 constexpr bool
1473 {
1474  return data != t.data;
1475 }
1476 
1477 
1478 
1479 template <int rank_, int dim, typename Number>
1480 template <typename OtherNumber>
1481 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1485 {
1486  data += t.data;
1487  return *this;
1488 }
1489 
1490 
1491 
1492 template <int rank_, int dim, typename Number>
1493 template <typename OtherNumber>
1494 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1498 {
1499  data -= t.data;
1500  return *this;
1501 }
1502 
1503 
1504 
1505 template <int rank_, int dim, typename Number>
1506 template <typename OtherNumber>
1507 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1510 {
1511  data *= d;
1512  return *this;
1513 }
1514 
1515 
1516 
1517 template <int rank_, int dim, typename Number>
1518 template <typename OtherNumber>
1519 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1522 {
1523  data /= d;
1524  return *this;
1525 }
1526 
1527 
1528 
1529 template <int rank_, int dim, typename Number>
1530 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1533 {
1534  SymmetricTensor tmp = *this;
1535  tmp.data = -tmp.data;
1536  return tmp;
1537 }
1538 
1539 
1540 
1541 template <int rank_, int dim, typename Number>
1542 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void
1544 {
1545  data.clear();
1546 }
1547 
1548 
1549 
1550 template <int rank_, int dim, typename Number>
1551 constexpr std::size_t
1553 {
1554  // all memory consists of statically allocated memory of the current
1555  // object, no pointers
1556  return sizeof(SymmetricTensor<rank_, dim, Number>);
1557 }
1558 
1559 
1560 
1561 namespace internal
1562 {
1563  template <int dim, typename Number, typename OtherNumber = Number>
1564  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1565  typename SymmetricTensorAccessors::
1566  double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1567  perform_double_contraction(
1568  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1569  base_tensor_type &data,
1570  const typename SymmetricTensorAccessors::
1571  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1572  {
1573  using result_type = typename SymmetricTensorAccessors::
1574  double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1575 
1576  switch (dim)
1577  {
1578  case 1:
1579  return data[0] * sdata[0];
1580  default:
1581  // Start with the non-diagonal part to avoid some multiplications by
1582  // 2.
1583 
1584  result_type sum = data[dim] * sdata[dim];
1585  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1586  sum += data[d] * sdata[d];
1587  sum += sum; // sum = sum * 2.;
1588 
1589  // Now add the contributions from the diagonal
1590  for (unsigned int d = 0; d < dim; ++d)
1591  sum += data[d] * sdata[d];
1592  return sum;
1593  }
1594  }
1595 
1596 
1597 
1598  template <int dim, typename Number, typename OtherNumber = Number>
1599  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1600  typename SymmetricTensorAccessors::
1601  double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1602  perform_double_contraction(
1603  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1604  base_tensor_type &data,
1605  const typename SymmetricTensorAccessors::
1606  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1607  {
1608  using result_type = typename SymmetricTensorAccessors::
1609  double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1610  using value_type = typename SymmetricTensorAccessors::
1611  double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1612 
1613  const unsigned int data_dim = SymmetricTensorAccessors::
1614  StorageType<2, dim, value_type>::n_independent_components;
1615  value_type tmp[data_dim]{};
1616  for (unsigned int i = 0; i < data_dim; ++i)
1617  tmp[i] =
1618  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1619  return result_type(tmp);
1620  }
1621 
1622 
1623 
1624  template <int dim, typename Number, typename OtherNumber = Number>
1625  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1626  typename SymmetricTensorAccessors::StorageType<
1627  2,
1628  dim,
1629  typename SymmetricTensorAccessors::
1630  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1631  base_tensor_type
1632  perform_double_contraction(
1633  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1634  base_tensor_type &data,
1635  const typename SymmetricTensorAccessors::
1636  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1637  {
1638  using value_type = typename SymmetricTensorAccessors::
1639  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1640  using base_tensor_type = typename SymmetricTensorAccessors::
1641  StorageType<2, dim, value_type>::base_tensor_type;
1642 
1643  base_tensor_type tmp;
1644  for (unsigned int i = 0; i < tmp.dimension; ++i)
1645  {
1646  // Start with the non-diagonal part
1647  value_type sum = data[dim] * sdata[dim][i];
1648  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1649  sum += data[d] * sdata[d][i];
1650  sum += sum; // sum = sum * 2.;
1651 
1652  // Now add the contributions from the diagonal
1653  for (unsigned int d = 0; d < dim; ++d)
1654  sum += data[d] * sdata[d][i];
1655  tmp[i] = sum;
1656  }
1657  return tmp;
1658  }
1659 
1660 
1661 
1662  template <int dim, typename Number, typename OtherNumber = Number>
1663  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1664  typename SymmetricTensorAccessors::StorageType<
1665  4,
1666  dim,
1667  typename SymmetricTensorAccessors::
1668  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1669  base_tensor_type
1670  perform_double_contraction(
1671  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1672  base_tensor_type &data,
1673  const typename SymmetricTensorAccessors::
1674  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1675  {
1676  using value_type = typename SymmetricTensorAccessors::
1677  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1678  using base_tensor_type = typename SymmetricTensorAccessors::
1679  StorageType<4, dim, value_type>::base_tensor_type;
1680 
1681  const unsigned int data_dim = SymmetricTensorAccessors::
1682  StorageType<2, dim, value_type>::n_independent_components;
1683  base_tensor_type tmp;
1684  for (unsigned int i = 0; i < data_dim; ++i)
1685  for (unsigned int j = 0; j < data_dim; ++j)
1686  {
1687  // Start with the non-diagonal part
1688  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1689  tmp[i][j] += data[i][d] * sdata[d][j];
1690  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1691 
1692  // Now add the contributions from the diagonal
1693  for (unsigned int d = 0; d < dim; ++d)
1694  tmp[i][j] += data[i][d] * sdata[d][j];
1695  }
1696  return tmp;
1697  }
1698 
1699 } // end of namespace internal
1700 
1701 
1702 
1703 template <int rank_, int dim, typename Number>
1704 template <typename OtherNumber>
1705 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
1706  typename internal::SymmetricTensorAccessors::
1707  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1710 {
1711  // need to have two different function calls
1712  // because a scalar and rank-2 tensor are not
1713  // the same data type (see internal function
1714  // above)
1715  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1716  s.data);
1717 }
1718 
1719 
1720 
1721 template <int rank_, int dim, typename Number>
1722 template <typename OtherNumber>
1723 DEAL_II_CONSTEXPR inline typename internal::SymmetricTensorAccessors::
1724  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1727 {
1728  typename internal::SymmetricTensorAccessors::
1729  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1730  tmp.data =
1731  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1732  s.data);
1733  return tmp;
1734 }
1735 
1736 
1737 
1738 // internal namespace to switch between the
1739 // access of different tensors. There used to
1740 // be explicit instantiations before for
1741 // different ranks and dimensions, but since
1742 // we now allow for templates on the data
1743 // type, and since we cannot partially
1744 // specialize the implementation, this got
1745 // into a separate namespace
1746 namespace internal
1747 {
1748  // The variables within this struct will be referenced in the next functions.
1749  // It is a workaround that allows returning a reference to a static variable
1750  // while allowing constexpr evaluation of the function.
1751  // It has to be defined outside the function because constexpr functions
1752  // cannot define static variables.
1753  // A similar struct has also been defined in tensor.h
1754  template <typename Type>
1755  struct Uninitialized
1756  {
1757  static Type value;
1758  };
1759 
1760  template <typename Type>
1761  Type Uninitialized<Type>::value;
1762 
1763  template <int dim, typename Number>
1764  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1765  symmetric_tensor_access(const TableIndices<2> &indices,
1766  typename SymmetricTensorAccessors::
1767  StorageType<2, dim, Number>::base_tensor_type &data)
1768  {
1769  // 1d is very simple and done first
1770  if (dim == 1)
1771  return data[0];
1772 
1773  // first treat the main diagonal elements, which are stored consecutively
1774  // at the beginning
1775  if (indices[0] == indices[1])
1776  return data[indices[0]];
1777 
1778  // the rest is messier and requires a few switches.
1779  switch (dim)
1780  {
1781  case 2:
1782  // at least for the 2x2 case it is reasonably simple
1783  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1784  ((indices[0] == 0) && (indices[1] == 1)),
1785  ExcInternalError());
1786  return data[2];
1787 
1788  default:
1789  // to do the rest, sort our indices before comparing
1790  {
1791  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1792  std::max(indices[0], indices[1]));
1793  for (unsigned int d = 0, c = 0; d < dim; ++d)
1794  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1795  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1796  return data[dim + c];
1797  Assert(false, ExcInternalError());
1798  }
1799  }
1800 
1801  // The code should never reach there.
1802  // Returns a dummy reference to a dummy variable just to make the
1803  // compiler happy.
1804  return Uninitialized<Number>::value;
1805  }
1806 
1807 
1808 
1809  template <int dim, typename Number>
1810  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1811  symmetric_tensor_access(const TableIndices<2> &indices,
1812  const typename SymmetricTensorAccessors::
1813  StorageType<2, dim, Number>::base_tensor_type &data)
1814  {
1815  // 1d is very simple and done first
1816  if (dim == 1)
1817  return data[0];
1818 
1819  // first treat the main diagonal elements, which are stored consecutively
1820  // at the beginning
1821  if (indices[0] == indices[1])
1822  return data[indices[0]];
1823 
1824  // the rest is messier and requires a few switches.
1825  switch (dim)
1826  {
1827  case 2:
1828  // at least for the 2x2 case it is reasonably simple
1829  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1830  ((indices[0] == 0) && (indices[1] == 1)),
1831  ExcInternalError());
1832  return data[2];
1833 
1834  default:
1835  // to do the rest, sort our indices before comparing
1836  {
1837  TableIndices<2> sorted_indices(std::min(indices[0], indices[1]),
1838  std::max(indices[0], indices[1]));
1839  for (unsigned int d = 0, c = 0; d < dim; ++d)
1840  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1841  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1842  return data[dim + c];
1843  Assert(false, ExcInternalError());
1844  }
1845  }
1846 
1847  // The code should never reach there.
1848  // Returns a dummy reference to a dummy variable just to make the
1849  // compiler happy.
1850  return Uninitialized<Number>::value;
1851  }
1852 
1853 
1854 
1855  template <int dim, typename Number>
1856  DEAL_II_CONSTEXPR inline Number &
1857  symmetric_tensor_access(const TableIndices<4> &indices,
1858  typename SymmetricTensorAccessors::
1859  StorageType<4, dim, Number>::base_tensor_type &data)
1860  {
1861  switch (dim)
1862  {
1863  case 1:
1864  return data[0][0];
1865 
1866  case 2:
1867  // each entry of the tensor can be thought of as an entry in a
1868  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1869  // rank-2 tensors. this is the format in which we store rank-4
1870  // tensors. determine which position the present entry is
1871  // stored in
1872  {
1873  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1874  return data[base_index[indices[0]][indices[1]]]
1875  [base_index[indices[2]][indices[3]]];
1876  }
1877  case 3:
1878  // each entry of the tensor can be thought of as an entry in a
1879  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1880  // rank-2 tensors. this is the format in which we store rank-4
1881  // tensors. determine which position the present entry is
1882  // stored in
1883  {
1884  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1885  {3, 1, 5},
1886  {4, 5, 2}};
1887  return data[base_index[indices[0]][indices[1]]]
1888  [base_index[indices[2]][indices[3]]];
1889  }
1890 
1891  default:
1892  Assert(false, ExcNotImplemented());
1893  }
1894 
1895  // The code should never reach there.
1896  // Returns a dummy reference to a dummy variable just to make the
1897  // compiler happy.
1898  return Uninitialized<Number>::value;
1899  }
1900 
1901 
1902  template <int dim, typename Number>
1903  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1904  symmetric_tensor_access(const TableIndices<4> &indices,
1905  const typename SymmetricTensorAccessors::
1906  StorageType<4, dim, Number>::base_tensor_type &data)
1907  {
1908  switch (dim)
1909  {
1910  case 1:
1911  return data[0][0];
1912 
1913  case 2:
1914  // each entry of the tensor can be thought of as an entry in a
1915  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1916  // rank-2 tensors. this is the format in which we store rank-4
1917  // tensors. determine which position the present entry is
1918  // stored in
1919  {
1920  constexpr std::size_t base_index[2][2] = {{0, 2}, {2, 1}};
1921  return data[base_index[indices[0]][indices[1]]]
1922  [base_index[indices[2]][indices[3]]];
1923  }
1924  case 3:
1925  // each entry of the tensor can be thought of as an entry in a
1926  // matrix that maps the rolled-out rank-2 tensors into rolled-out
1927  // rank-2 tensors. this is the format in which we store rank-4
1928  // tensors. determine which position the present entry is
1929  // stored in
1930  {
1931  constexpr std::size_t base_index[3][3] = {{0, 3, 4},
1932  {3, 1, 5},
1933  {4, 5, 2}};
1934  return data[base_index[indices[0]][indices[1]]]
1935  [base_index[indices[2]][indices[3]]];
1936  }
1937 
1938  default:
1939  Assert(false, ExcNotImplemented());
1940  }
1941 
1942  // The code should never reach there.
1943  // Returns a dummy reference to a dummy variable just to make the
1944  // compiler happy.
1945  return Uninitialized<Number>::value;
1946  }
1947 
1948 } // end of namespace internal
1949 
1950 
1951 
1952 template <int rank_, int dim, typename Number>
1953 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
1955  operator()(const TableIndices<rank_> &indices)
1956 {
1957  for (unsigned int r = 0; r < rank; ++r)
1958  Assert(indices[r] < dimension, ExcIndexRange(indices[r], 0, dimension));
1959  return internal::symmetric_tensor_access<dim, Number>(indices, data);
1960 }
1961 
1962 
1963 
1964 template <int rank_, int dim, typename Number>
1965 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE const Number &
1967  operator()(const TableIndices<rank_> &indices) const
1968 {
1969  for (unsigned int r = 0; r < rank; ++r)
1970  Assert(indices[r] < dimension, ExcIndexRange(indices[r], 0, dimension));
1971  return internal::symmetric_tensor_access<dim, Number>(indices, data);
1972 }
1973 
1974 
1975 
1976 namespace internal
1977 {
1978  namespace SymmetricTensorImplementation
1979  {
1980  template <int rank_>
1981  constexpr TableIndices<rank_>
1982  get_partially_filled_indices(const unsigned int row,
1983  const std::integral_constant<int, 2> &)
1984  {
1986  }
1987 
1988 
1989  template <int rank_>
1990  constexpr TableIndices<rank_>
1991  get_partially_filled_indices(const unsigned int row,
1992  const std::integral_constant<int, 4> &)
1993  {
1994  return TableIndices<rank_>(row,
1998  }
1999  } // namespace SymmetricTensorImplementation
2000 } // namespace internal
2001 
2002 
2003 template <int rank_, int dim, typename Number>
2004 constexpr DEAL_II_ALWAYS_INLINE internal::SymmetricTensorAccessors::
2005  Accessor<rank_, dim, true, rank_ - 1, Number>
2007  operator[](const unsigned int row) const
2008 {
2009  return internal::SymmetricTensorAccessors::
2010  Accessor<rank_, dim, true, rank_ - 1, Number>(
2011  *this,
2012  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2013  rank_>(row, std::integral_constant<int, rank_>()));
2014 }
2015 
2016 
2017 
2018 template <int rank_, int dim, typename Number>
2019 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE internal::
2020  SymmetricTensorAccessors::Accessor<rank_, dim, false, rank_ - 1, Number>
2021  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2022 {
2023  return internal::SymmetricTensorAccessors::
2024  Accessor<rank_, dim, false, rank_ - 1, Number>(
2025  *this,
2026  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2027  rank_>(row, std::integral_constant<int, rank_>()));
2028 }
2029 
2030 
2031 
2032 template <int rank_, int dim, typename Number>
2033 constexpr DEAL_II_ALWAYS_INLINE const Number &
2035  operator[](const TableIndices<rank_> &indices) const
2036 {
2037  return operator()(indices);
2038 }
2039 
2040 
2041 
2042 template <int rank_, int dim, typename Number>
2043 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number &
2045  operator[](const TableIndices<rank_> &indices)
2046 {
2047  return operator()(indices);
2048 }
2049 
2050 
2051 
2052 template <int rank_, int dim, typename Number>
2053 inline Number *
2055 {
2056  return std::addressof(this->access_raw_entry(0));
2057 }
2058 
2059 
2060 
2061 template <int rank_, int dim, typename Number>
2062 inline const Number *
2064 {
2065  return std::addressof(this->access_raw_entry(0));
2066 }
2067 
2068 
2069 
2070 template <int rank_, int dim, typename Number>
2071 inline Number *
2073 {
2075 }
2076 
2077 
2078 
2079 template <int rank_, int dim, typename Number>
2080 inline const Number *
2082 {
2084 }
2085 
2086 
2087 
2088 namespace internal
2089 {
2090  namespace SymmetricTensorImplementation
2091  {
2092  template <int dim, typename Number>
2093  constexpr unsigned int
2094  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2095  const unsigned int index)
2096  {
2097  return index;
2098  }
2099 
2100 
2101  template <int dim, typename Number>
2102  constexpr ::TableIndices<2>
2103  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2104  const unsigned int index)
2105  {
2108  }
2109 
2110  } // namespace SymmetricTensorImplementation
2111 } // namespace internal
2112 
2113 
2114 
2115 template <int rank_, int dim, typename Number>
2116 DEAL_II_CONSTEXPR inline const Number &
2118  const unsigned int index) const
2119 {
2120  AssertIndexRange(index, n_independent_components);
2121  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2122  index)];
2123 }
2124 
2125 
2126 
2127 template <int rank_, int dim, typename Number>
2128 DEAL_II_CONSTEXPR inline Number &
2130 {
2131  AssertIndexRange(index, n_independent_components);
2132  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2133  index)];
2134 }
2135 
2136 
2137 
2138 namespace internal
2139 {
2140  template <int dim, typename Number>
2141  DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
2142  compute_norm(const typename SymmetricTensorAccessors::
2143  StorageType<2, dim, Number>::base_tensor_type &data)
2144  {
2145  switch (dim)
2146  {
2147  case 1:
2148  return numbers::NumberTraits<Number>::abs(data[0]);
2149 
2150  case 2:
2151  return std::sqrt(
2155 
2156  case 3:
2157  return std::sqrt(
2164 
2165  default:
2166  {
2167  typename numbers::NumberTraits<Number>::real_type return_value =
2169 
2170  for (unsigned int d = 0; d < dim; ++d)
2171  return_value +=
2173  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2174  return_value +=
2176 
2177  return std::sqrt(return_value);
2178  }
2179  }
2180  }
2181 
2182 
2183 
2184  template <int dim, typename Number>
2185  DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
2186  compute_norm(const typename SymmetricTensorAccessors::
2187  StorageType<4, dim, Number>::base_tensor_type &data)
2188  {
2189  switch (dim)
2190  {
2191  case 1:
2192  return numbers::NumberTraits<Number>::abs(data[0][0]);
2193 
2194  default:
2195  {
2196  typename numbers::NumberTraits<Number>::real_type return_value =
2198 
2199  const unsigned int n_independent_components = data.dimension;
2200 
2201  for (unsigned int i = 0; i < dim; ++i)
2202  for (unsigned int j = 0; j < dim; ++j)
2203  return_value +=
2205  for (unsigned int i = 0; i < dim; ++i)
2206  for (unsigned int j = dim; j < n_independent_components; ++j)
2207  return_value +=
2209  for (unsigned int i = dim; i < n_independent_components; ++i)
2210  for (unsigned int j = 0; j < dim; ++j)
2211  return_value +=
2213  for (unsigned int i = dim; i < n_independent_components; ++i)
2214  for (unsigned int j = dim; j < n_independent_components; ++j)
2215  return_value +=
2217 
2218  return std::sqrt(return_value);
2219  }
2220  }
2221  }
2222 
2223 } // end of namespace internal
2224 
2225 
2226 
2227 template <int rank_, int dim, typename Number>
2230 {
2231  return internal::compute_norm<dim, Number>(data);
2232 }
2233 
2234 
2235 
2236 namespace internal
2237 {
2238  namespace SymmetricTensorImplementation
2239  {
2240  // a function to do the unrolling from a set of indices to a
2241  // scalar index into the array in which we store the elements of
2242  // a symmetric tensor
2243  //
2244  // this function is for rank-2 tensors
2245  template <int dim>
2246  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE unsigned int
2248  {
2249  Assert(indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2250  Assert(indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2251 
2252  switch (dim)
2253  {
2254  case 1:
2255  {
2256  return 0;
2257  }
2258 
2259  case 2:
2260  {
2261  constexpr unsigned int table[2][2] = {{0, 2}, {2, 1}};
2262  return table[indices[0]][indices[1]];
2263  }
2264 
2265  case 3:
2266  {
2267  constexpr unsigned int table[3][3] = {{0, 3, 4},
2268  {3, 1, 5},
2269  {4, 5, 2}};
2270  return table[indices[0]][indices[1]];
2271  }
2272 
2273  case 4:
2274  {
2275  constexpr unsigned int table[4][4] = {{0, 4, 5, 6},
2276  {4, 1, 7, 8},
2277  {5, 7, 2, 9},
2278  {6, 8, 9, 3}};
2279  return table[indices[0]][indices[1]];
2280  }
2281 
2282  default:
2283  // for the remainder, manually figure out the numbering
2284  {
2285  if (indices[0] == indices[1])
2286  return indices[0];
2287 
2288  TableIndices<2> sorted_indices(indices);
2289  sorted_indices.sort();
2290 
2291  for (unsigned int d = 0, c = 0; d < dim; ++d)
2292  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2293  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2294  return dim + c;
2295 
2296  // should never get here:
2297  Assert(false, ExcInternalError());
2298  return 0;
2299  }
2300  }
2301  }
2302 
2303  // a function to do the unrolling from a set of indices to a
2304  // scalar index into the array in which we store the elements of
2305  // a symmetric tensor
2306  //
2307  // this function is for tensors of ranks not already handled
2308  // above
2309  template <int dim, int rank_>
2310  DEAL_II_CONSTEXPR inline unsigned int
2312  {
2313  (void)indices;
2314  Assert(false, ExcNotImplemented());
2316  }
2317  } // namespace SymmetricTensorImplementation
2318 } // namespace internal
2319 
2320 
2321 template <int rank_, int dim, typename Number>
2322 constexpr unsigned int
2324  const TableIndices<rank_> &indices)
2325 {
2326  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2327  dim>(indices);
2328 }
2329 
2330 
2331 
2332 namespace internal
2333 {
2334  namespace SymmetricTensorImplementation
2335  {
2336  // a function to do the inverse of the unrolling from a set of
2337  // indices to a scalar index into the array in which we store
2338  // the elements of a symmetric tensor. in other words, it goes
2339  // from the scalar index into the array to a set of indices of
2340  // the tensor
2341  //
2342  // this function is for rank-2 tensors
2343  template <int dim>
2344  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE TableIndices<2>
2345  unrolled_to_component_indices(const unsigned int i,
2346  const std::integral_constant<int, 2> &)
2347  {
2348  Assert(
2350  ExcIndexRange(
2351  i,
2352  0,
2354  switch (dim)
2355  {
2356  case 1:
2357  {
2358  return {0, 0};
2359  }
2360 
2361  case 2:
2362  {
2363  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2364  TableIndices<2>(1, 1),
2365  TableIndices<2>(0, 1)};
2366  return table[i];
2367  }
2368 
2369  case 3:
2370  {
2371  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2372  TableIndices<2>(1, 1),
2373  TableIndices<2>(2, 2),
2374  TableIndices<2>(0, 1),
2375  TableIndices<2>(0, 2),
2376  TableIndices<2>(1, 2)};
2377  return table[i];
2378  }
2379 
2380  default:
2381  if (i < dim)
2382  return {i, i};
2383 
2384  for (unsigned int d = 0, c = 0; d < dim; ++d)
2385  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2386  if (c == i)
2387  return {d, e};
2388 
2389  // should never get here:
2390  Assert(false, ExcInternalError());
2391  return {0, 0};
2392  }
2393  }
2394 
2395  // a function to do the inverse of the unrolling from a set of
2396  // indices to a scalar index into the array in which we store
2397  // the elements of a symmetric tensor. in other words, it goes
2398  // from the scalar index into the array to a set of indices of
2399  // the tensor
2400  //
2401  // this function is for tensors of a rank not already handled
2402  // above
2403  template <int dim, int rank_>
2404  DEAL_II_CONSTEXPR inline TableIndices<rank_>
2405  unrolled_to_component_indices(const unsigned int i,
2406  const std::integral_constant<int, rank_> &)
2407  {
2408  (void)i;
2409  Assert(
2410  (i <
2412  ExcIndexRange(i,
2413  0,
2415  n_independent_components));
2416  Assert(false, ExcNotImplemented());
2417  return TableIndices<rank_>();
2418  }
2419 
2420  } // namespace SymmetricTensorImplementation
2421 } // namespace internal
2422 
2423 template <int rank_, int dim, typename Number>
2424 constexpr DEAL_II_ALWAYS_INLINE TableIndices<rank_>
2426  const unsigned int i)
2427 {
2428  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2429  dim>(i, std::integral_constant<int, rank_>());
2430 }
2431 
2432 
2433 
2434 template <int rank_, int dim, typename Number>
2435 template <class Archive>
2436 inline void
2437 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2438 {
2439  ar &data;
2440 }
2441 
2442 
2443 #endif // DOXYGEN
2444 
2445 /* ----------------- Non-member functions operating on tensors. ------------ */
2446 
2447 
2460 template <int rank_, int dim, typename Number, typename OtherNumber>
2461 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2465 {
2467  tmp = left;
2468  tmp += right;
2469  return tmp;
2470 }
2471 
2472 
2485 template <int rank_, int dim, typename Number, typename OtherNumber>
2486 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
2490 {
2492  tmp = left;
2493  tmp -= right;
2494  return tmp;
2495 }
2496 
2497 
2505 template <int rank_, int dim, typename Number, typename OtherNumber>
2506 constexpr DEAL_II_ALWAYS_INLINE
2509  const Tensor<rank_, dim, OtherNumber> & right)
2510 {
2511  return Tensor<rank_, dim, Number>(left) + right;
2512 }
2513 
2514 
2522 template <int rank_, int dim, typename Number, typename OtherNumber>
2523 constexpr DEAL_II_ALWAYS_INLINE
2527 {
2528  return left + Tensor<rank_, dim, OtherNumber>(right);
2529 }
2530 
2531 
2539 template <int rank_, int dim, typename Number, typename OtherNumber>
2540 constexpr DEAL_II_ALWAYS_INLINE
2543  const Tensor<rank_, dim, OtherNumber> & right)
2544 {
2545  return Tensor<rank_, dim, Number>(left) - right;
2546 }
2547 
2548 
2556 template <int rank_, int dim, typename Number, typename OtherNumber>
2557 constexpr DEAL_II_ALWAYS_INLINE
2561 {
2562  return left - Tensor<rank_, dim, OtherNumber>(right);
2563 }
2564 
2565 
2566 
2580 template <int dim, typename Number>
2581 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2583 {
2584  switch (dim)
2585  {
2586  case 1:
2587  return t.data[0];
2588  case 2:
2589  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2590  case 3:
2591  {
2592  // in analogy to general tensors, but
2593  // there's something to be simplified for
2594  // the present case
2595  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2596  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2597  t.data[0] * t.data[5] * t.data[5] -
2598  t.data[1] * t.data[4] * t.data[4] -
2599  t.data[2] * t.data[3] * t.data[3]);
2600  }
2601  default:
2602  Assert(false, ExcNotImplemented());
2603  return internal::NumberType<Number>::value(0.0);
2604  }
2605 }
2606 
2607 
2608 
2618 template <int dim, typename Number>
2619 constexpr DEAL_II_ALWAYS_INLINE Number
2621 {
2622  return determinant(t);
2623 }
2624 
2625 
2626 
2634 template <int dim, typename Number>
2635 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
2637 {
2638  Number t = d.data[0];
2639  for (unsigned int i = 1; i < dim; ++i)
2640  t += d.data[i];
2641  return t;
2642 }
2643 
2644 
2654 template <int dim, typename Number>
2655 constexpr Number
2657 {
2658  return trace(t);
2659 }
2660 
2661 
2674 template <typename Number>
2675 constexpr DEAL_II_ALWAYS_INLINE Number
2677 {
2678  return internal::NumberType<Number>::value(0.0);
2679 }
2680 
2681 
2682 
2702 template <typename Number>
2703 constexpr DEAL_II_ALWAYS_INLINE Number
2705 {
2706  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2707 }
2708 
2709 
2710 
2720 template <typename Number>
2721 constexpr DEAL_II_ALWAYS_INLINE Number
2723 {
2724  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2725  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2726 }
2727 
2728 
2729 
2738 template <typename Number>
2739 std::array<Number, 1>
2741 
2742 
2743 
2767 template <typename Number>
2768 std::array<Number, 2>
2770 
2771 
2772 
2795 template <typename Number>
2796 std::array<Number, 3>
2798 
2799 
2800 
2801 namespace internal
2802 {
2803  namespace SymmetricTensorImplementation
2804  {
2844  template <int dim, typename Number>
2845  void
2846  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2847  ::Tensor<2, dim, Number> & Q,
2848  std::array<Number, dim> & d,
2849  std::array<Number, dim - 1> & e);
2850 
2851 
2852 
2894  template <int dim, typename Number>
2895  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2896  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2897 
2898 
2899 
2941  template <int dim, typename Number>
2942  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2944 
2945 
2946 
2962  template <typename Number>
2963  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
2964  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
2965 
2966 
2967 
3002  template <typename Number>
3003  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3004  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3005 
3010  template <int dim, typename Number>
3012  {
3013  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3014  bool
3015  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3016  {
3017  return lhs.first > rhs.first;
3018  }
3019  };
3020 
3021  } // namespace SymmetricTensorImplementation
3022 
3023 } // namespace internal
3024 
3025 
3026 
3027 // The line below is to ensure that doxygen puts the full description
3028 // of this global enumeration into the documentation
3029 // See https://stackoverflow.com/a/1717984
3059 {
3069  hybrid,
3087  jacobi
3088 };
3089 
3090 
3091 
3121 template <int dim, typename Number>
3122 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3123  std::integral_constant<int, dim>::value>
3125  const SymmetricTensorEigenvectorMethod method =
3127 
3128 
3129 
3139 template <int rank_, int dim, typename Number>
3140 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
3142 {
3143  return t;
3144 }
3145 
3146 
3147 
3157 template <int dim, typename Number>
3158 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3160 {
3162 
3163  // subtract scaled trace from the diagonal
3164  const Number tr = trace(t) / dim;
3165  for (unsigned int i = 0; i < dim; ++i)
3166  tmp.data[i] -= tr;
3167 
3168  return tmp;
3169 }
3170 
3171 
3172 
3180 template <int dim, typename Number>
3181 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3183 {
3184  // create a default constructed matrix filled with
3185  // zeros, then set the diagonal elements to one
3187  switch (dim)
3188  {
3189  case 1:
3190  tmp.data[0] = Number(1);
3191  break;
3192  case 2:
3193  tmp.data[0] = tmp.data[1] = Number(1);
3194  break;
3195  case 3:
3196  tmp.data[0] = tmp.data[1] = tmp.data[2] = Number(1);
3197  break;
3198  default:
3199  for (unsigned int d = 0; d < dim; ++d)
3200  tmp.data[d] = Number(1);
3201  }
3202  return tmp;
3203 }
3204 
3205 
3206 
3215 template <int dim>
3216 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim>
3218 {
3219  return unit_symmetric_tensor<dim, double>();
3220 }
3221 
3222 
3223 
3238 template <int dim, typename Number>
3239 DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
3241 {
3243 
3244  // fill the elements treating the diagonal
3245  for (unsigned int i = 0; i < dim; ++i)
3246  for (unsigned int j = 0; j < dim; ++j)
3247  tmp.data[i][j] = Number((i == j ? 1 : 0) - 1. / dim);
3248 
3249  // then fill the ones that copy over the
3250  // non-diagonal elements. note that during
3251  // the double-contraction, we handle the
3252  // off-diagonal elements twice, so simply
3253  // copying requires a weight of 1/2
3254  for (unsigned int i = dim;
3255  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3256  n_rank2_components;
3257  ++i)
3258  tmp.data[i][i] = Number(0.5);
3259 
3260  return tmp;
3261 }
3262 
3263 
3264 
3279 template <int dim>
3280 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
3282 {
3283  return deviator_tensor<dim, double>();
3284 }
3285 
3286 
3287 
3310 template <int dim, typename Number>
3311 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim, Number>
3313 {
3315 
3316  // fill the elements treating the diagonal
3317  for (unsigned int i = 0; i < dim; ++i)
3318  tmp.data[i][i] = Number(1);
3319 
3320  // then fill the ones that copy over the
3321  // non-diagonal elements. note that during
3322  // the double-contraction, we handle the
3323  // off-diagonal elements twice, so simply
3324  // copying requires a weight of 1/2
3325  for (unsigned int i = dim;
3326  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3327  n_rank2_components;
3328  ++i)
3329  tmp.data[i][i] = Number(0.5);
3330 
3331  return tmp;
3332 }
3333 
3334 
3335 
3357 template <int dim>
3358 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<4, dim>
3360 {
3361  return identity_tensor<dim, double>();
3362 }
3363 
3364 
3365 
3376 template <int dim, typename Number>
3377 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3379 {
3381  value(t);
3382 }
3383 
3384 
3385 
3397 template <int dim, typename Number>
3400 {
3402  value(t);
3403 }
3404 
3405 
3406 
3421 template <int dim, typename Number>
3422 DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
3425 {
3427 
3428  // fill only the elements really needed
3429  for (unsigned int i = 0; i < dim; ++i)
3430  for (unsigned int j = i; j < dim; ++j)
3431  for (unsigned int k = 0; k < dim; ++k)
3432  for (unsigned int l = k; l < dim; ++l)
3433  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3434 
3435  return tmp;
3436 }
3437 
3438 
3439 
3448 template <int dim, typename Number>
3449 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, Number>
3451 {
3453  for (unsigned int d = 0; d < dim; ++d)
3454  result[d][d] = t[d][d];
3455  Number half = 0.5;
3456  for (unsigned int d = 0; d < dim; ++d)
3457  for (unsigned int e = d + 1; e < dim; ++e)
3458  result[d][e] = (t[d][e] + t[e][d]) * half;
3459  return result;
3460 }
3461 
3462 
3463 
3471 template <int rank_, int dim, typename Number>
3472 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
3474  operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3475 {
3477  tt *= factor;
3478  return tt;
3479 }
3480 
3481 
3482 
3490 template <int rank_, int dim, typename Number>
3491 constexpr DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
3492  operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3493 {
3494  // simply forward to the other operator
3495  return t * factor;
3496 }
3497 
3498 
3499 
3525 template <int rank_, int dim, typename Number, typename OtherNumber>
3526 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3527  rank_,
3528  dim,
3529  typename ProductType<Number,
3530  typename EnableIfScalar<OtherNumber>::type>::type>
3532  const OtherNumber & factor)
3533 {
3534  // form the product. we have to convert the two factors into the final
3535  // type via explicit casts because, for awkward reasons, the C++
3536  // standard committee saw it fit to not define an
3537  // operator*(float,std::complex<double>)
3538  // (as well as with switched arguments and double<->float).
3539  using product_type = typename ProductType<Number, OtherNumber>::type;
3541  tt *= product_type(factor);
3542  return tt;
3543 }
3544 
3545 
3546 
3555 template <int rank_, int dim, typename Number, typename OtherNumber>
3556 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<
3557  rank_,
3558  dim,
3559  typename ProductType<OtherNumber,
3560  typename EnableIfScalar<Number>::type>::type>
3561 operator*(const Number & factor,
3563 {
3564  // simply forward to the other operator with switched arguments
3565  return (t * factor);
3566 }
3567 
3568 
3569 
3575 template <int rank_, int dim, typename Number, typename OtherNumber>
3576 DEAL_II_CONSTEXPR inline SymmetricTensor<
3577  rank_,
3578  dim,
3579  typename ProductType<Number,
3580  typename EnableIfScalar<OtherNumber>::type>::type>
3582  const OtherNumber & factor)
3583 {
3585  tt = t;
3586  tt /= factor;
3587  return tt;
3588 }
3589 
3590 
3591 
3598 template <int rank_, int dim>
3599 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim>
3600  operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3601 {
3603  tt *= factor;
3604  return tt;
3605 }
3606 
3607 
3608 
3615 template <int rank_, int dim>
3616 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim>
3617  operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3618 {
3620  tt *= factor;
3621  return tt;
3622 }
3623 
3624 
3625 
3631 template <int rank_, int dim>
3632 DEAL_II_CONSTEXPR inline SymmetricTensor<rank_, dim>
3633 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3634 {
3636  tt /= factor;
3637  return tt;
3638 }
3639 
3649 template <int dim, typename Number, typename OtherNumber>
3650 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
3653 {
3654  return (t1 * t2);
3655 }
3656 
3657 
3667 template <int dim, typename Number, typename OtherNumber>
3668 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
3669  typename ProductType<Number, OtherNumber>::type
3671  const Tensor<2, dim, OtherNumber> & t2)
3672 {
3673  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<
3674  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3675  for (unsigned int i = 0; i < dim; ++i)
3676  for (unsigned int j = 0; j < dim; ++j)
3677  s += t1[i][j] * t2[i][j];
3678  return s;
3679 }
3680 
3681 
3691 template <int dim, typename Number, typename OtherNumber>
3692 constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
3695 {
3696  return scalar_product(t2, t1);
3697 }
3698 
3699 
3715 template <typename Number, typename OtherNumber>
3716 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void double_contract(
3717  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3720 {
3721  tmp[0][0] = t[0][0][0][0] * s[0][0];
3722 }
3723 
3724 
3725 
3741 template <typename Number, typename OtherNumber>
3742 DEAL_II_CONSTEXPR inline void double_contract(
3743  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3746 {
3747  tmp[0][0] = t[0][0][0][0] * s[0][0];
3748 }
3749 
3750 
3751 
3766 template <typename Number, typename OtherNumber>
3767 DEAL_II_CONSTEXPR inline void double_contract(
3768  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3771 {
3772  const unsigned int dim = 2;
3773 
3774  for (unsigned int i = 0; i < dim; ++i)
3775  for (unsigned int j = i; j < dim; ++j)
3776  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3777  2 * t[i][j][0][1] * s[0][1];
3778 }
3779 
3780 
3781 
3797 template <typename Number, typename OtherNumber>
3798 DEAL_II_CONSTEXPR inline void double_contract(
3799  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3802 {
3803  const unsigned int dim = 2;
3804 
3805  for (unsigned int i = 0; i < dim; ++i)
3806  for (unsigned int j = i; j < dim; ++j)
3807  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3808  2 * s[0][1] * t[0][1][i][j];
3809 }
3810 
3811 
3812 
3828 template <typename Number, typename OtherNumber>
3829 DEAL_II_CONSTEXPR inline void double_contract(
3830  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3833 {
3834  const unsigned int dim = 3;
3835 
3836  for (unsigned int i = 0; i < dim; ++i)
3837  for (unsigned int j = i; j < dim; ++j)
3838  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3839  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3840  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3841 }
3842 
3843 
3844 
3860 template <typename Number, typename OtherNumber>
3861 DEAL_II_CONSTEXPR inline void double_contract(
3862  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3865 {
3866  const unsigned int dim = 3;
3867 
3868  for (unsigned int i = 0; i < dim; ++i)
3869  for (unsigned int j = i; j < dim; ++j)
3870  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3871  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3872  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3873 }
3874 
3875 
3876 
3884 template <int dim, typename Number, typename OtherNumber>
3885 DEAL_II_CONSTEXPR
3888  const Tensor<1, dim, OtherNumber> & src2)
3889 {
3891  for (unsigned int i = 0; i < dim; ++i)
3892  for (unsigned int j = 0; j < dim; ++j)
3893  dest[i] += src1[i][j] * src2[j];
3894  return dest;
3895 }
3896 
3897 
3905 template <int dim, typename Number, typename OtherNumber>
3909 {
3910  // this is easy for symmetric tensors:
3911  return src2 * src1;
3912 }
3913 
3914 
3915 
3936 template <int rank_1,
3937  int rank_2,
3938  int dim,
3939  typename Number,
3940  typename OtherNumber>
3941 constexpr DEAL_II_ALWAYS_INLINE
3942  typename Tensor<rank_1 + rank_2 - 2,
3943  dim,
3944  typename ProductType<Number, OtherNumber>::type>::tensor_type
3947 {
3948  return src1 * Tensor<rank_2, dim, OtherNumber>(src2);
3949 }
3950 
3951 
3952 
3973 template <int rank_1,
3974  int rank_2,
3975  int dim,
3976  typename Number,
3977  typename OtherNumber>
3978 constexpr DEAL_II_ALWAYS_INLINE
3979  typename Tensor<rank_1 + rank_2 - 2,
3980  dim,
3981  typename ProductType<Number, OtherNumber>::type>::tensor_type
3983  const Tensor<rank_2, dim, OtherNumber> & src2)
3984 {
3985  return Tensor<rank_2, dim, OtherNumber>(src1) * src2;
3986 }
3987 
3988 
3989 
3999 template <int dim, typename Number>
4000 inline std::ostream &
4001 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4002 {
4003  // make our lives a bit simpler by outputting
4004  // the tensor through the operator for the
4005  // general Tensor class
4007 
4008  for (unsigned int i = 0; i < dim; ++i)
4009  for (unsigned int j = 0; j < dim; ++j)
4010  tt[i][j] = t[i][j];
4011 
4012  return out << tt;
4013 }
4014 
4015 
4016 
4026 template <int dim, typename Number>
4027 inline std::ostream &
4028 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4029 {
4030  // make our lives a bit simpler by outputting
4031  // the tensor through the operator for the
4032  // general Tensor class
4034 
4035  for (unsigned int i = 0; i < dim; ++i)
4036  for (unsigned int j = 0; j < dim; ++j)
4037  for (unsigned int k = 0; k < dim; ++k)
4038  for (unsigned int l = 0; l < dim; ++l)
4039  tt[i][j][k][l] = t[i][j][k][l];
4040 
4041  return out << tt;
4042 }
4043 
4044 
4045 DEAL_II_NAMESPACE_CLOSE
4046 
4047 #endif
DEAL_II_CONSTEXPR internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
DEAL_II_CONSTEXPR SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
static const unsigned int invalid_unsigned_int
Definition: types.h:187
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
DEAL_II_CONSTEXPR Number & operator()(const TableIndices< rank_ > &indices)
static constexpr unsigned int n_independent_components
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
DEAL_II_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)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
constexpr numbers::NumberTraits< Number >::real_type norm() const
DEAL_II_CONSTEXPR SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
constexpr bool operator==(const SymmetricTensor &) const
friend DEAL_II_CONSTEXPR SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
DEAL_II_CONSTEXPR void clear()
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
DEAL_II_CONSTEXPR SymmetricTensor & operator/=(const OtherNumber &factor)
DEAL_II_CONSTEXPR SymmetricTensor & operator*=(const OtherNumber &factor)
static real_type abs(const number &x)
Definition: numbers.h:627
SymmetricTensorEigenvectorMethod
DEAL_II_CONSTEXPR const Number & access_raw_entry(const unsigned int unrolled_index) const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
DEAL_II_CONSTEXPR SymmetricTensor operator-() const
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
constexpr internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
static ::ExceptionBase & ExcMessage(std::string arg1)
constexpr bool operator!=(const SymmetricTensor &) const
constexpr SymmetricTensor()=default
typename base_tensor_descriptor::base_tensor_type base_tensor_type
static const unsigned int rank
constexpr Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Number * begin_raw()
constexpr Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
#define Assert(cond, exc)
Definition: exceptions.h:1411
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
DEAL_II_CONSTEXPR SymmetricTensor< 4, dim, Number > identity_tensor()
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &t)
base_tensor_type data
DEAL_II_CONSTEXPR Number trace(const SymmetricTensor< 2, dim, Number > &d)
void serialize(Archive &ar, const unsigned int version)
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A, ::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim - 1 > &e)
Number * end_raw()
DEAL_II_CONSTEXPR SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
DEAL_II_CONSTEXPR TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
static constexpr std::size_t memory_consumption()
constexpr Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: tensor.h:422
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:959
DEAL_II_CONSTEXPR Number trace(const SymmetricTensor< 2, dim, Number > &d)
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
friend DEAL_II_CONSTEXPR SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
DEAL_II_CONSTEXPR SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static ::ExceptionBase & ExcNotImplemented()
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
friend DEAL_II_CONSTEXPR SymmetricTensor< 4, dim2, Number2 > identity_tensor()
static ::ExceptionBase & ExcInternalError()
DEAL_II_CONSTEXPR SymmetricTensor< 4, dim, Number > deviator_tensor()