Reference documentation for deal.II version 9.1.1
\(\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 template <int rank, int dim, typename Number = double>
32 class SymmetricTensor;
33 
34 template <int dim, typename Number>
35 SymmetricTensor<2, dim, Number>
36 unit_symmetric_tensor();
37 
38 template <int dim, typename Number>
39 SymmetricTensor<4, dim, Number>
40 deviator_tensor();
41 
42 template <int dim, typename Number>
43 SymmetricTensor<4, dim, Number>
44 identity_tensor();
45 
46 template <int dim, typename Number>
47 SymmetricTensor<2, dim, Number>
48 invert(const SymmetricTensor<2, dim, Number> &);
49 
50 template <int dim, typename Number>
51 SymmetricTensor<4, dim, Number>
52 invert(const SymmetricTensor<4, dim, Number> &);
53 
54 template <int dim2, typename Number>
55 Number
56 trace(const SymmetricTensor<2, dim2, Number> &);
57 
58 template <int dim, typename Number>
59 SymmetricTensor<2, dim, Number>
60 deviator(const SymmetricTensor<2, dim, Number> &);
61 
62 template <int dim, typename Number>
63 Number
64 determinant(const SymmetricTensor<2, dim, Number> &);
65 
66 
67 
68 namespace internal
69 {
74  namespace SymmetricTensorImplementation
75  {
80  template <int rank, int dim, typename Number>
81  struct Inverse;
82  } // namespace SymmetricTensorImplementation
83 
88  namespace SymmetricTensorAccessors
89  {
96  inline TableIndices<2>
97  merge(const TableIndices<2> &previous_indices,
98  const unsigned int new_index,
99  const unsigned int position)
100  {
101  Assert(position < 2, ExcIndexRange(position, 0, 2));
102 
103  if (position == 0)
104  return {new_index, numbers::invalid_unsigned_int};
105  else
106  return {previous_indices[0], new_index};
107  }
108 
109 
110 
117  inline TableIndices<4>
118  merge(const TableIndices<4> &previous_indices,
119  const unsigned int new_index,
120  const unsigned int position)
121  {
122  Assert(position < 4, ExcIndexRange(position, 0, 4));
123 
124  switch (position)
125  {
126  case 0:
127  return {new_index,
128  numbers::invalid_unsigned_int,
129  numbers::invalid_unsigned_int,
130  numbers::invalid_unsigned_int};
131  case 1:
132  return {previous_indices[0],
133  new_index,
134  numbers::invalid_unsigned_int,
135  numbers::invalid_unsigned_int};
136  case 2:
137  return {previous_indices[0],
138  previous_indices[1],
139  new_index,
140  numbers::invalid_unsigned_int};
141  case 3:
142  return {previous_indices[0],
143  previous_indices[1],
144  previous_indices[2],
145  new_index};
146  }
147  Assert(false, ExcInternalError());
148  return {};
149  }
150 
151 
160  template <int rank1,
161  int rank2,
162  int dim,
163  typename Number,
164  typename OtherNumber = Number>
165  struct double_contraction_result
166  {
167  using value_type = typename ProductType<Number, OtherNumber>::type;
168  using type =
169  ::SymmetricTensor<rank1 + rank2 - 4, dim, value_type>;
170  };
171 
172 
181  template <int dim, typename Number, typename OtherNumber>
182  struct double_contraction_result<2, 2, dim, Number, OtherNumber>
183  {
184  using type = typename ProductType<Number, OtherNumber>::type;
185  };
186 
187 
188 
201  template <int rank, int dim, typename Number>
202  struct StorageType;
203 
207  template <int dim, typename Number>
208  struct StorageType<2, dim, Number>
209  {
214  static const unsigned int n_independent_components =
215  (dim * dim + dim) / 2;
216 
220  using base_tensor_type = Tensor<1, n_independent_components, Number>;
221  };
222 
223 
224 
228  template <int dim, typename Number>
229  struct StorageType<4, dim, Number>
230  {
236  static const unsigned int n_rank2_components = (dim * dim + dim) / 2;
237 
241  static const unsigned int n_independent_components =
242  (n_rank2_components *
243  StorageType<2, dim, Number>::n_independent_components);
244 
251  using base_tensor_type = Tensor<2, n_rank2_components, Number>;
252  };
253 
254 
255 
260  template <int rank, int dim, bool constness, typename Number>
261  struct AccessorTypes;
262 
269  template <int rank, int dim, typename Number>
270  struct AccessorTypes<rank, dim, true, Number>
271  {
272  using tensor_type = const ::SymmetricTensor<rank, dim, Number>;
273 
274  using reference = Number;
275  };
276 
283  template <int rank, int dim, typename Number>
284  struct AccessorTypes<rank, dim, false, Number>
285  {
286  using tensor_type = ::SymmetricTensor<rank, dim, Number>;
287 
288  using reference = Number &;
289  };
290 
291 
326  template <int rank, int dim, bool constness, int P, typename Number>
327  class Accessor
328  {
329  public:
333  using reference =
334  typename AccessorTypes<rank, dim, constness, Number>::reference;
335  using tensor_type =
336  typename AccessorTypes<rank, dim, constness, Number>::tensor_type;
337 
338  private:
357  Accessor(tensor_type &tensor, const TableIndices<rank> &previous_indices);
358 
362  Accessor(const Accessor &) = default;
363 
364  public:
368  Accessor<rank, dim, constness, P - 1, Number>
369  operator[](const unsigned int i);
370 
374  Accessor<rank, dim, constness, P - 1, Number>
375  operator[](const unsigned int i) const;
376 
377  private:
381  tensor_type & tensor;
382  const TableIndices<rank> previous_indices;
383 
384  // declare some other classes
385  // as friends. make sure to
386  // work around bugs in some
387  // compilers
388  template <int, int, typename>
389  friend class ::SymmetricTensor;
390  template <int, int, bool, int, typename>
391  friend class Accessor;
392 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
393  friend class ::SymmetricTensor<rank, dim, Number>;
394  friend class Accessor<rank, dim, constness, P + 1, Number>;
395 #endif
396  };
397 
398 
399 
409  template <int rank, int dim, bool constness, typename Number>
410  class Accessor<rank, dim, constness, 1, Number>
411  {
412  public:
416  using reference =
417  typename AccessorTypes<rank, dim, constness, Number>::reference;
418  using tensor_type =
419  typename AccessorTypes<rank, dim, constness, Number>::tensor_type;
420 
421  private:
443  Accessor(tensor_type &tensor, const TableIndices<rank> &previous_indices);
444 
448  Accessor(const Accessor &) = default;
449 
450  public:
454  reference operator[](const unsigned int);
455 
459  reference operator[](const unsigned int) const;
460 
461  private:
465  tensor_type & tensor;
466  const TableIndices<rank> previous_indices;
467 
468  // declare some other classes
469  // as friends. make sure to
470  // work around bugs in some
471  // compilers
472  template <int, int, typename>
473  friend class ::SymmetricTensor;
474  template <int, int, bool, int, typename>
475  friend class SymmetricTensorAccessors::Accessor;
476 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
477  friend class ::SymmetricTensor<rank, dim, Number>;
478  friend class SymmetricTensorAccessors::
479  Accessor<rank, dim, constness, 2, Number>;
480 #endif
481  };
482  } // namespace SymmetricTensorAccessors
483 } // namespace internal
484 
485 
486 
550 template <int rank_, int dim, typename Number>
551 class SymmetricTensor
552 {
553 public:
554  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
555 
564  static const unsigned int dimension = dim;
565 
569  static const unsigned int rank = rank_;
570 
576  static constexpr unsigned int n_independent_components =
577  internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>::
578  n_independent_components;
579 
583  SymmetricTensor();
584 
595  template <typename OtherNumber>
596  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
597 
613  SymmetricTensor(const Number (&array)[n_independent_components]);
614 
620  template <typename OtherNumber>
621  explicit SymmetricTensor(
622  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
623 
627  Number *
628  begin_raw();
629 
633  const Number *
634  begin_raw() const;
635 
639  Number *
640  end_raw();
641 
646  const Number *
647  end_raw() const;
648 
655  template <typename OtherNumber>
656  SymmetricTensor &
657  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
658 
665  SymmetricTensor &
666  operator=(const Number &d);
667 
672  operator Tensor<rank_, dim, Number>() const;
673 
677  bool
678  operator==(const SymmetricTensor &) const;
679 
683  bool
684  operator!=(const SymmetricTensor &) const;
685 
689  template <typename OtherNumber>
690  SymmetricTensor &
691  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
692 
696  template <typename OtherNumber>
697  SymmetricTensor &
698  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
699 
704  template <typename OtherNumber>
705  SymmetricTensor &
706  operator*=(const OtherNumber &factor);
707 
711  template <typename OtherNumber>
712  SymmetricTensor &
713  operator/=(const OtherNumber &factor);
714 
718  SymmetricTensor
719  operator-() const;
720 
745  template <typename OtherNumber>
746  typename internal::SymmetricTensorAccessors::
747  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
748  operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const;
749 
754  template <typename OtherNumber>
755  typename internal::SymmetricTensorAccessors::
756  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
757  operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const;
758 
762  Number &
763  operator()(const TableIndices<rank_> &indices);
764 
768  const Number &
769  operator()(const TableIndices<rank_> &indices) const;
770 
775  internal::SymmetricTensorAccessors::
776  Accessor<rank_, dim, true, rank_ - 1, Number>
777  operator[](const unsigned int row) const;
778 
783  internal::SymmetricTensorAccessors::
784  Accessor<rank_, dim, false, rank_ - 1, Number>
785  operator[](const unsigned int row);
786 
792  const Number &operator[](const TableIndices<rank_> &indices) const;
793 
799  Number &operator[](const TableIndices<rank_> &indices);
800 
806  const Number &
807  access_raw_entry(const unsigned int unrolled_index) const;
808 
814  Number &
815  access_raw_entry(const unsigned int unrolled_index);
816 
826  typename numbers::NumberTraits<Number>::real_type
827  norm() const;
828 
836  static unsigned int
837  component_to_unrolled_index(const TableIndices<rank_> &indices);
838 
844  static TableIndices<rank_>
845  unrolled_to_component_indices(const unsigned int i);
846 
859  void
860  clear();
861 
866  static std::size_t
867  memory_consumption();
868 
873  template <class Archive>
874  void
875  serialize(Archive &ar, const unsigned int version);
876 
877 private:
881  using base_tensor_descriptor =
882  internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>;
883 
887  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
888 
892  base_tensor_type data;
893 
897  template <int, int, typename>
898  friend class SymmetricTensor;
899 
903  template <int dim2, typename Number2>
904  friend Number2
905  trace(const SymmetricTensor<2, dim2, Number2> &d);
906 
907  template <int dim2, typename Number2>
908  friend Number2
909  determinant(const SymmetricTensor<2, dim2, Number2> &t);
910 
911  template <int dim2, typename Number2>
912  friend SymmetricTensor<2, dim2, Number2>
913  deviator(const SymmetricTensor<2, dim2, Number2> &t);
914 
915  template <int dim2, typename Number2>
916  friend SymmetricTensor<2, dim2, Number2>
917  unit_symmetric_tensor();
918 
919  template <int dim2, typename Number2>
920  friend SymmetricTensor<4, dim2, Number2>
921  deviator_tensor();
922 
923  template <int dim2, typename Number2>
924  friend SymmetricTensor<4, dim2, Number2>
925  identity_tensor();
926 
927 
931  friend struct internal::SymmetricTensorImplementation::
932  Inverse<2, dim, Number>;
933 
934  friend struct internal::SymmetricTensorImplementation::
935  Inverse<4, dim, Number>;
936 };
937 
938 
939 
940 // ------------------------- inline functions ------------------------
941 
942 #ifndef DOXYGEN
943 
944 // provide declarations for static members
945 template <int rank, int dim, typename Number>
946 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
947 
948 template <int rank_, int dim, typename Number>
949 constexpr unsigned int
950  SymmetricTensor<rank_, dim, Number>::n_independent_components;
951 
952 namespace internal
953 {
954  namespace SymmetricTensorAccessors
955  {
956  template <int rank_, int dim, bool constness, int P, typename Number>
957  Accessor<rank_, dim, constness, P, Number>::Accessor(
958  tensor_type & tensor,
959  const TableIndices<rank_> &previous_indices)
960  : tensor(tensor)
961  , previous_indices(previous_indices)
962  {}
963 
964 
965 
966  template <int rank_, int dim, bool constness, int P, typename Number>
967  Accessor<rank_, dim, constness, P - 1, Number>
968  Accessor<rank_, dim, constness, P, Number>::
969  operator[](const unsigned int i)
970  {
971  return Accessor<rank_, dim, constness, P - 1, Number>(
972  tensor, merge(previous_indices, i, rank_ - P));
973  }
974 
975 
976 
977  template <int rank_, int dim, bool constness, int P, typename Number>
978  Accessor<rank_, dim, constness, P - 1, Number>
979  Accessor<rank_, dim, constness, P, Number>::
980  operator[](const unsigned int i) const
981  {
982  return Accessor<rank_, dim, constness, P - 1, Number>(
983  tensor, merge(previous_indices, i, rank_ - P));
984  }
985 
986 
987 
988  template <int rank_, int dim, bool constness, typename Number>
989  Accessor<rank_, dim, constness, 1, Number>::Accessor(
990  tensor_type & tensor,
991  const TableIndices<rank_> &previous_indices)
992  : tensor(tensor)
993  , previous_indices(previous_indices)
994  {}
995 
996 
997 
998  template <int rank_, int dim, bool constness, typename Number>
999  typename Accessor<rank_, dim, constness, 1, Number>::reference
1000  Accessor<rank_, dim, constness, 1, Number>::
1001  operator[](const unsigned int i)
1002  {
1003  return tensor(merge(previous_indices, i, rank_ - 1));
1004  }
1005 
1006 
1007  template <int rank_, int dim, bool constness, typename Number>
1008  typename Accessor<rank_, dim, constness, 1, Number>::reference
1009  Accessor<rank_, dim, constness, 1, Number>::
1010  operator[](const unsigned int i) const
1011  {
1012  return tensor(merge(previous_indices, i, rank_ - 1));
1013  }
1014  } // namespace SymmetricTensorAccessors
1015 } // namespace internal
1016 
1017 
1018 
1019 template <int rank_, int dim, typename Number>
1020 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor()
1021 {
1022  // Some auto-differentiable numbers need explicit
1023  // zero initialization.
1024  for (unsigned int i = 0; i < base_tensor_type::dimension; ++i)
1025  data[i] = internal::NumberType<Number>::value(0.0);
1026 }
1027 
1028 
1029 template <int rank_, int dim, typename Number>
1030 template <typename OtherNumber>
1031 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1032  const Tensor<2, dim, OtherNumber> &t)
1033 {
1034  Assert(rank == 2, ExcNotImplemented());
1035  switch (dim)
1036  {
1037  case 2:
1038  Assert(t[0][1] == t[1][0], ExcInternalError());
1039 
1040  data[0] = t[0][0];
1041  data[1] = t[1][1];
1042  data[2] = t[0][1];
1043 
1044  break;
1045  case 3:
1046  Assert(t[0][1] == t[1][0], ExcInternalError());
1047  Assert(t[0][2] == t[2][0], ExcInternalError());
1048  Assert(t[1][2] == t[2][1], ExcInternalError());
1049 
1050  data[0] = t[0][0];
1051  data[1] = t[1][1];
1052  data[2] = t[2][2];
1053  data[3] = t[0][1];
1054  data[4] = t[0][2];
1055  data[5] = t[1][2];
1056 
1057  break;
1058  default:
1059  for (unsigned int d = 0; d < dim; ++d)
1060  for (unsigned int e = 0; e < d; ++e)
1061  Assert(t[d][e] == t[e][d], ExcInternalError());
1062 
1063  for (unsigned int d = 0; d < dim; ++d)
1064  data[d] = t[d][d];
1065 
1066  for (unsigned int d = 0, c = 0; d < dim; ++d)
1067  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1068  data[dim + c] = t[d][e];
1069  }
1070 }
1071 
1072 
1073 
1074 template <int rank_, int dim, typename Number>
1075 template <typename OtherNumber>
1076 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1077  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1078 {
1079  for (unsigned int i = 0; i < base_tensor_type::dimension; ++i)
1080  data[i] =
1081  internal::NumberType<typename base_tensor_type::value_type>::value(
1082  initializer.data[i]);
1083 }
1084 
1085 
1086 
1087 template <int rank_, int dim, typename Number>
1088 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1089  const Number (&array)[n_independent_components])
1090  : data(
1091  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1092 {
1093  // ensure that the reinterpret_cast above actually works
1094  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1095  ExcInternalError());
1096 }
1097 
1098 
1099 
1100 template <int rank_, int dim, typename Number>
1101 template <typename OtherNumber>
1102 inline SymmetricTensor<rank_, dim, Number> &
1103 SymmetricTensor<rank_, dim, Number>::
1104 operator=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1105 {
1106  for (unsigned int i = 0; i < base_tensor_type::dimension; ++i)
1107  data[i] = t.data[i];
1108  return *this;
1109 }
1110 
1111 
1112 
1113 template <int rank_, int dim, typename Number>
1114 inline SymmetricTensor<rank_, dim, Number> &
1115 SymmetricTensor<rank_, dim, Number>::operator=(const Number &d)
1116 {
1117  Assert(numbers::value_is_zero(d),
1118  ExcMessage("Only assignment with zero is allowed"));
1119  (void)d;
1120 
1121  data = internal::NumberType<Number>::value(0.0);
1122 
1123  return *this;
1124 }
1125 
1126 
1127 namespace internal
1128 {
1129  namespace SymmetricTensorImplementation
1130  {
1131  template <int dim, typename Number>
1132  inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1133  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1134  {
1135  ::Tensor<2, dim, Number> t;
1136 
1137  // diagonal entries are stored first
1138  for (unsigned int d = 0; d < dim; ++d)
1139  t[d][d] = s.access_raw_entry(d);
1140 
1141  // off-diagonal entries come next, row by row
1142  for (unsigned int d = 0, c = 0; d < dim; ++d)
1143  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1144  {
1145  t[d][e] = s.access_raw_entry(dim + c);
1146  t[e][d] = s.access_raw_entry(dim + c);
1147  }
1148  return t;
1149  }
1150 
1151 
1152  template <int dim, typename Number>
1153  ::Tensor<4, dim, Number>
1154  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1155  {
1156  // utilize the symmetry properties of SymmetricTensor<4,dim>
1157  // discussed in the class documentation to avoid accessing all
1158  // independent elements of the input tensor more than once
1159  ::Tensor<4, dim, Number> t;
1160 
1161  for (unsigned int i = 0; i < dim; ++i)
1162  for (unsigned int j = i; j < dim; ++j)
1163  for (unsigned int k = 0; k < dim; ++k)
1164  for (unsigned int l = k; l < dim; ++l)
1165  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1166  t[TableIndices<4>(j, i, k, l)] =
1167  t[TableIndices<4>(j, i, l, k)] =
1168  st[TableIndices<4>(i, j, k, l)];
1169 
1170  return t;
1171  }
1172 
1173 
1174  template <typename Number>
1175  struct Inverse<2, 1, Number>
1176  {
1177  static inline ::SymmetricTensor<2, 1, Number>
1178  value(const ::SymmetricTensor<2, 1, Number> &t)
1179  {
1180  ::SymmetricTensor<2, 1, Number> tmp;
1181 
1182  tmp[0][0] = 1.0 / t[0][0];
1183 
1184  return tmp;
1185  }
1186  };
1187 
1188 
1189  template <typename Number>
1190  struct Inverse<2, 2, Number>
1191  {
1192  static inline ::SymmetricTensor<2, 2, Number>
1193  value(const ::SymmetricTensor<2, 2, Number> &t)
1194  {
1195  ::SymmetricTensor<2, 2, Number> tmp;
1196 
1197  // Sympy result: ([
1198  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1199  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1200  const TableIndices<2> idx_00(0, 0);
1201  const TableIndices<2> idx_01(0, 1);
1202  const TableIndices<2> idx_11(1, 1);
1203  const Number inv_det_t =
1204  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1205  tmp[idx_00] = t[idx_11];
1206  tmp[idx_01] = -t[idx_01];
1207  tmp[idx_11] = t[idx_00];
1208  tmp *= inv_det_t;
1209 
1210  return tmp;
1211  }
1212  };
1213 
1214 
1215  template <typename Number>
1216  struct Inverse<2, 3, Number>
1217  {
1218  static ::SymmetricTensor<2, 3, Number>
1219  value(const ::SymmetricTensor<2, 3, Number> &t)
1220  {
1221  ::SymmetricTensor<2, 3, Number> tmp;
1222 
1223  // Sympy result: ([
1224  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1225  // 2*t01*t02*t12 - t02**2*t11),
1226  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1227  // 2*t01*t02*t12 - t02**2*t11),
1228  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1229  // 2*t01*t02*t12 - t02**2*t11)],
1230  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1231  // 2*t01*t02*t12 - t02**2*t11),
1232  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1233  // 2*t01*t02*t12 - t02**2*t11),
1234  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1235  // 2*t01*t02*t12 + t02**2*t11)],
1236  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1237  // 2*t01*t02*t12 - t02**2*t11),
1238  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1239  // 2*t01*t02*t12 + t02**2*t11),
1240  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1241  // 2*t01*t02*t12 + t02**2*t11)] ])
1242  //
1243  // =
1244  //
1245  // [ (t11*t22 - t12**2)/det_t,
1246  // (-t01*t22 + t02*t12)/det_t,
1247  // (t01*t12 - t02*t11)/det_t],
1248  // [ (-t01*t22 + t02*t12)/det_t,
1249  // (t00*t22 - t02**2)/det_t,
1250  // (-t00*t12 + t01*t02)/det_t],
1251  // [ (t01*t12 - t02*t11)/det_t,
1252  // (-t00*t12 + t01*t02)/det_t,
1253  // (t00*t11 - t01**2)/det_t] ])
1254  //
1255  // with det_t = (t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1256  // 2*t01*t02*t12 - t02**2*t11)
1257  const TableIndices<2> idx_00(0, 0);
1258  const TableIndices<2> idx_01(0, 1);
1259  const TableIndices<2> idx_02(0, 2);
1260  const TableIndices<2> idx_11(1, 1);
1261  const TableIndices<2> idx_12(1, 2);
1262  const TableIndices<2> idx_22(2, 2);
1263  const Number inv_det_t =
1264  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1265  t[idx_00] * t[idx_12] * t[idx_12] -
1266  t[idx_01] * t[idx_01] * t[idx_22] +
1267  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1268  t[idx_02] * t[idx_02] * t[idx_11]);
1269  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1270  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1271  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1272  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1273  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1274  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1275  tmp *= inv_det_t;
1276 
1277  return tmp;
1278  }
1279  };
1280 
1281 
1282  template <typename Number>
1283  struct Inverse<4, 1, Number>
1284  {
1285  static inline ::SymmetricTensor<4, 1, Number>
1286  value(const ::SymmetricTensor<4, 1, Number> &t)
1287  {
1288  ::SymmetricTensor<4, 1, Number> tmp;
1289  tmp.data[0][0] = 1.0 / t.data[0][0];
1290  return tmp;
1291  }
1292  };
1293 
1294 
1295  template <typename Number>
1296  struct Inverse<4, 2, Number>
1297  {
1298  static inline ::SymmetricTensor<4, 2, Number>
1299  value(const ::SymmetricTensor<4, 2, Number> &t)
1300  {
1301  ::SymmetricTensor<4, 2, Number> tmp;
1302 
1303  // Inverting this tensor is a little more complicated than necessary,
1304  // since we store the data of 't' as a 3x3 matrix t.data, but the
1305  // product between a rank-4 and a rank-2 tensor is really not the
1306  // product between this matrix and the 3-vector of a rhs, but rather
1307  //
1308  // B.vec = t.data * mult * A.vec
1309  //
1310  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1311  // capture the fact that we need to add up both the c_ij12*a_12 and the
1312  // c_ij21*a_21 terms.
1313  //
1314  // In addition, in this scheme, the identity tensor has the matrix
1315  // representation mult^-1.
1316  //
1317  // The inverse of 't' therefore has the matrix representation
1318  //
1319  // inv.data = mult^-1 * t.data^-1 * mult^-1
1320  //
1321  // in order to compute it, let's first compute the inverse of t.data and
1322  // put it into tmp.data; at the end of the function we then scale the
1323  // last row and column of the inverse by 1/2, corresponding to the left
1324  // and right multiplication with mult^-1.
1325  const Number t4 = t.data[0][0] * t.data[1][1],
1326  t6 = t.data[0][0] * t.data[1][2],
1327  t8 = t.data[0][1] * t.data[1][0],
1328  t00 = t.data[0][2] * t.data[1][0],
1329  t01 = t.data[0][1] * t.data[2][0],
1330  t04 = t.data[0][2] * t.data[2][0],
1331  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1332  t8 * t.data[2][2] + t00 * t.data[2][1] +
1333  t01 * t.data[1][2] - t04 * t.data[1][1]);
1334  tmp.data[0][0] =
1335  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1336  tmp.data[0][1] =
1337  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1338  tmp.data[0][2] =
1339  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1340  tmp.data[1][0] =
1341  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1342  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1343  tmp.data[1][2] = -(t6 - t00) * t07;
1344  tmp.data[2][0] =
1345  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1346  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1347  tmp.data[2][2] = (t4 - t8) * t07;
1348 
1349  // scale last row and column as mentioned
1350  // above
1351  tmp.data[2][0] /= 2;
1352  tmp.data[2][1] /= 2;
1353  tmp.data[0][2] /= 2;
1354  tmp.data[1][2] /= 2;
1355  tmp.data[2][2] /= 4;
1356 
1357  return tmp;
1358  }
1359  };
1360 
1361 
1362  template <typename Number>
1363  struct Inverse<4, 3, Number>
1364  {
1365  static ::SymmetricTensor<4, 3, Number>
1366  value(const ::SymmetricTensor<4, 3, Number> &t)
1367  {
1368  ::SymmetricTensor<4, 3, Number> tmp = t;
1369 
1370  // This function follows the exact same scheme as the 2d case, except
1371  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1372  // Instead, we use the Gauss-Jordan algorithm implemented for
1373  // FullMatrix. For historical reasons the following code is copied from
1374  // there, with the tangential benefit that we do not need to copy the
1375  // tensor entries to and from the FullMatrix.
1376  const unsigned int N = 6;
1377 
1378  // First get an estimate of the size of the elements of this matrix,
1379  // for later checks whether the pivot element is large enough, or
1380  // whether we have to fear that the matrix is not regular.
1381  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1382  for (unsigned int i = 0; i < N; ++i)
1383  diagonal_sum += std::fabs(tmp.data[i][i]);
1384  const Number typical_diagonal_element =
1385  diagonal_sum / static_cast<double>(N);
1386  (void)typical_diagonal_element;
1387 
1388  unsigned int p[N];
1389  for (unsigned int i = 0; i < N; ++i)
1390  p[i] = i;
1391 
1392  for (unsigned int j = 0; j < N; ++j)
1393  {
1394  // Pivot search: search that part of the line on and right of the
1395  // diagonal for the largest element.
1396  Number max = std::fabs(tmp.data[j][j]);
1397  unsigned int r = j;
1398  for (unsigned int i = j + 1; i < N; ++i)
1399  if (std::fabs(tmp.data[i][j]) > max)
1400  {
1401  max = std::fabs(tmp.data[i][j]);
1402  r = i;
1403  }
1404 
1405  // Check whether the pivot is too small
1406  Assert(max > 1.e-16 * typical_diagonal_element,
1407  ExcMessage("This tensor seems to be noninvertible"));
1408 
1409  // Row interchange
1410  if (r > j)
1411  {
1412  for (unsigned int k = 0; k < N; ++k)
1413  std::swap(tmp.data[j][k], tmp.data[r][k]);
1414 
1415  std::swap(p[j], p[r]);
1416  }
1417 
1418  // Transformation
1419  const Number hr = 1. / tmp.data[j][j];
1420  tmp.data[j][j] = hr;
1421  for (unsigned int k = 0; k < N; ++k)
1422  {
1423  if (k == j)
1424  continue;
1425  for (unsigned int i = 0; i < N; ++i)
1426  {
1427  if (i == j)
1428  continue;
1429  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1430  }
1431  }
1432  for (unsigned int i = 0; i < N; ++i)
1433  {
1434  tmp.data[i][j] *= hr;
1435  tmp.data[j][i] *= -hr;
1436  }
1437  tmp.data[j][j] = hr;
1438  }
1439 
1440  // Column interchange
1441  Number hv[N];
1442  for (unsigned int i = 0; i < N; ++i)
1443  {
1444  for (unsigned int k = 0; k < N; ++k)
1445  hv[p[k]] = tmp.data[i][k];
1446  for (unsigned int k = 0; k < N; ++k)
1447  tmp.data[i][k] = hv[k];
1448  }
1449 
1450  // Scale rows and columns. The mult matrix
1451  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1452  for (unsigned int i = 3; i < 6; ++i)
1453  for (unsigned int j = 0; j < 3; ++j)
1454  tmp.data[i][j] /= 2;
1455 
1456  for (unsigned int i = 0; i < 3; ++i)
1457  for (unsigned int j = 3; j < 6; ++j)
1458  tmp.data[i][j] /= 2;
1459 
1460  for (unsigned int i = 3; i < 6; ++i)
1461  for (unsigned int j = 3; j < 6; ++j)
1462  tmp.data[i][j] /= 4;
1463 
1464  return tmp;
1465  }
1466  };
1467 
1468  } // namespace SymmetricTensorImplementation
1469 } // namespace internal
1470 
1471 
1472 
1473 template <int rank_, int dim, typename Number>
1474 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>::
1475  operator Tensor<rank_, dim, Number>() const
1476 {
1477  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1478 }
1479 
1480 
1481 
1482 template <int rank_, int dim, typename Number>
1483 inline bool
1484 SymmetricTensor<rank_, dim, Number>::
1485 operator==(const SymmetricTensor<rank_, dim, Number> &t) const
1486 {
1487  return data == t.data;
1488 }
1489 
1490 
1491 
1492 template <int rank_, int dim, typename Number>
1493 inline bool
1494 SymmetricTensor<rank_, dim, Number>::
1495 operator!=(const SymmetricTensor<rank_, dim, Number> &t) const
1496 {
1497  return data != t.data;
1498 }
1499 
1500 
1501 
1502 template <int rank_, int dim, typename Number>
1503 template <typename OtherNumber>
1504 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
1505  SymmetricTensor<rank_, dim, Number>::
1506  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1507 {
1508  data += t.data;
1509  return *this;
1510 }
1511 
1512 
1513 
1514 template <int rank_, int dim, typename Number>
1515 template <typename OtherNumber>
1516 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
1517  SymmetricTensor<rank_, dim, Number>::
1518  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1519 {
1520  data -= t.data;
1521  return *this;
1522 }
1523 
1524 
1525 
1526 template <int rank_, int dim, typename Number>
1527 template <typename OtherNumber>
1528 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
1529 SymmetricTensor<rank_, dim, Number>::operator*=(const OtherNumber &d)
1530 {
1531  data *= d;
1532  return *this;
1533 }
1534 
1535 
1536 
1537 template <int rank_, int dim, typename Number>
1538 template <typename OtherNumber>
1539 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
1540 SymmetricTensor<rank_, dim, Number>::operator/=(const OtherNumber &d)
1541 {
1542  data /= d;
1543  return *this;
1544 }
1545 
1546 
1547 
1548 template <int rank_, int dim, typename Number>
1549 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
1550 SymmetricTensor<rank_, dim, Number>::operator-() const
1551 {
1552  SymmetricTensor tmp = *this;
1553  tmp.data = -tmp.data;
1554  return tmp;
1555 }
1556 
1557 
1558 
1559 template <int rank_, int dim, typename Number>
1560 inline DEAL_II_ALWAYS_INLINE void
1561 SymmetricTensor<rank_, dim, Number>::clear()
1562 {
1563  data.clear();
1564 }
1565 
1566 
1567 
1568 template <int rank_, int dim, typename Number>
1569 inline std::size_t
1570 SymmetricTensor<rank_, dim, Number>::memory_consumption()
1571 {
1572  // all memory consists of statically allocated memory of the current
1573  // object, no pointers
1574  return sizeof(SymmetricTensor<rank_, dim, Number>);
1575 }
1576 
1577 
1578 
1579 namespace internal
1580 {
1581  template <int dim, typename Number, typename OtherNumber = Number>
1582  inline DEAL_II_ALWAYS_INLINE typename SymmetricTensorAccessors::
1583  double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1584  perform_double_contraction(
1585  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1586  base_tensor_type &data,
1587  const typename SymmetricTensorAccessors::
1588  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1589  {
1590  using result_type = typename SymmetricTensorAccessors::
1591  double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1592 
1593  switch (dim)
1594  {
1595  case 1:
1596  return data[0] * sdata[0];
1597  default:
1598  // Start with the non-diagonal part to avoid some multiplications by
1599  // 2.
1600 
1601  result_type sum = data[dim] * sdata[dim];
1602  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1603  sum += data[d] * sdata[d];
1604  sum += sum; // sum = sum * 2.;
1605 
1606  // Now add the contributions from the diagonal
1607  for (unsigned int d = 0; d < dim; ++d)
1608  sum += data[d] * sdata[d];
1609  return sum;
1610  }
1611  }
1612 
1613 
1614 
1615  template <int dim, typename Number, typename OtherNumber = Number>
1616  inline typename SymmetricTensorAccessors::
1617  double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1618  perform_double_contraction(
1619  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1620  base_tensor_type &data,
1621  const typename SymmetricTensorAccessors::
1622  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1623  {
1624  using result_type = typename SymmetricTensorAccessors::
1625  double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1626  using value_type = typename SymmetricTensorAccessors::
1627  double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1628 
1629  const unsigned int data_dim = SymmetricTensorAccessors::
1630  StorageType<2, dim, value_type>::n_independent_components;
1631  value_type tmp[data_dim];
1632  for (unsigned int i = 0; i < data_dim; ++i)
1633  tmp[i] =
1634  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1635  return result_type(tmp);
1636  }
1637 
1638 
1639 
1640  template <int dim, typename Number, typename OtherNumber = Number>
1641  inline typename SymmetricTensorAccessors::StorageType<
1642  2,
1643  dim,
1644  typename SymmetricTensorAccessors::
1645  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1646  base_tensor_type
1647  perform_double_contraction(
1648  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1649  base_tensor_type &data,
1650  const typename SymmetricTensorAccessors::
1651  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1652  {
1653  using value_type = typename SymmetricTensorAccessors::
1654  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1655  using base_tensor_type = typename SymmetricTensorAccessors::
1656  StorageType<2, dim, value_type>::base_tensor_type;
1657 
1658  base_tensor_type tmp;
1659  for (unsigned int i = 0; i < tmp.dimension; ++i)
1660  {
1661  // Start with the non-diagonal part
1662  value_type sum = data[dim] * sdata[dim][i];
1663  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1664  sum += data[d] * sdata[d][i];
1665  sum += sum; // sum = sum * 2.;
1666 
1667  // Now add the contributions from the diagonal
1668  for (unsigned int d = 0; d < dim; ++d)
1669  sum += data[d] * sdata[d][i];
1670  tmp[i] = sum;
1671  }
1672  return tmp;
1673  }
1674 
1675 
1676 
1677  template <int dim, typename Number, typename OtherNumber = Number>
1678  inline typename SymmetricTensorAccessors::StorageType<
1679  4,
1680  dim,
1681  typename SymmetricTensorAccessors::
1682  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1683  base_tensor_type
1684  perform_double_contraction(
1685  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1686  base_tensor_type &data,
1687  const typename SymmetricTensorAccessors::
1688  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1689  {
1690  using value_type = typename SymmetricTensorAccessors::
1691  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1692  using base_tensor_type = typename SymmetricTensorAccessors::
1693  StorageType<4, dim, value_type>::base_tensor_type;
1694 
1695  const unsigned int data_dim = SymmetricTensorAccessors::
1696  StorageType<2, dim, value_type>::n_independent_components;
1697  base_tensor_type tmp;
1698  for (unsigned int i = 0; i < data_dim; ++i)
1699  for (unsigned int j = 0; j < data_dim; ++j)
1700  {
1701  // Start with the non-diagonal part
1702  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1703  tmp[i][j] += data[i][d] * sdata[d][j];
1704  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1705 
1706  // Now add the contributions from the diagonal
1707  for (unsigned int d = 0; d < dim; ++d)
1708  tmp[i][j] += data[i][d] * sdata[d][j];
1709  }
1710  return tmp;
1711  }
1712 
1713 } // end of namespace internal
1714 
1715 
1716 
1717 template <int rank_, int dim, typename Number>
1718 template <typename OtherNumber>
1719 inline DEAL_II_ALWAYS_INLINE typename internal::SymmetricTensorAccessors::
1720  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1721  SymmetricTensor<rank_, dim, Number>::
1722  operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const
1723 {
1724  // need to have two different function calls
1725  // because a scalar and rank-2 tensor are not
1726  // the same data type (see internal function
1727  // above)
1728  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1729  s.data);
1730 }
1731 
1732 
1733 
1734 template <int rank_, int dim, typename Number>
1735 template <typename OtherNumber>
1736 inline typename internal::SymmetricTensorAccessors::
1737  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1738  SymmetricTensor<rank_, dim, Number>::
1739  operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const
1740 {
1741  typename internal::SymmetricTensorAccessors::
1742  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1743  tmp.data =
1744  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1745  s.data);
1746  return tmp;
1747 }
1748 
1749 
1750 
1751 // internal namespace to switch between the
1752 // access of different tensors. There used to
1753 // be explicit instantiations before for
1754 // different ranks and dimensions, but since
1755 // we now allow for templates on the data
1756 // type, and since we cannot partially
1757 // specialize the implementation, this got
1758 // into a separate namespace
1759 namespace internal
1760 {
1761  template <int dim, typename Number>
1762  inline Number &
1763  symmetric_tensor_access(const TableIndices<2> &indices,
1764  typename SymmetricTensorAccessors::
1765  StorageType<2, dim, Number>::base_tensor_type &data)
1766  {
1767  // 1d is very simple and done first
1768  if (dim == 1)
1769  return data[0];
1770 
1771  // first treat the main diagonal elements, which are stored consecutively
1772  // at the beginning
1773  if (indices[0] == indices[1])
1774  return data[indices[0]];
1775 
1776  // the rest is messier and requires a few switches.
1777  switch (dim)
1778  {
1779  case 2:
1780  // at least for the 2x2 case it is reasonably simple
1781  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1782  ((indices[0] == 0) && (indices[1] == 1)),
1783  ExcInternalError());
1784  return data[2];
1785 
1786  default:
1787  // to do the rest, sort our indices before comparing
1788  {
1789  TableIndices<2> sorted_indices(indices);
1790  sorted_indices.sort();
1791 
1792  for (unsigned int d = 0, c = 0; d < dim; ++d)
1793  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1794  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1795  return data[dim + c];
1796  Assert(false, ExcInternalError());
1797  }
1798  }
1799 
1800  static Number dummy_but_referenceable = Number();
1801  return dummy_but_referenceable;
1802  }
1803 
1804 
1805 
1806  template <int dim, typename Number>
1807  inline const Number &
1808  symmetric_tensor_access(const TableIndices<2> &indices,
1809  const typename SymmetricTensorAccessors::
1810  StorageType<2, dim, Number>::base_tensor_type &data)
1811  {
1812  // 1d is very simple and done first
1813  if (dim == 1)
1814  return data[0];
1815 
1816  // first treat the main diagonal elements, which are stored consecutively
1817  // at the beginning
1818  if (indices[0] == indices[1])
1819  return data[indices[0]];
1820 
1821  // the rest is messier and requires a few switches.
1822  switch (dim)
1823  {
1824  case 2:
1825  // at least for the 2x2 case it is reasonably simple
1826  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1827  ((indices[0] == 0) && (indices[1] == 1)),
1828  ExcInternalError());
1829  return data[2];
1830 
1831  default:
1832  // to do the rest, sort our indices before comparing
1833  {
1834  TableIndices<2> sorted_indices(indices);
1835  sorted_indices.sort();
1836 
1837  for (unsigned int d = 0, c = 0; d < dim; ++d)
1838  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1839  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1840  return data[dim + c];
1841  Assert(false, ExcInternalError());
1842  }
1843  }
1844 
1845  static Number dummy_but_referenceable = Number();
1846  return dummy_but_referenceable;
1847  }
1848 
1849 
1850 
1851  template <int dim, typename Number>
1852  inline Number &
1853  symmetric_tensor_access(const TableIndices<4> &indices,
1854  typename SymmetricTensorAccessors::
1855  StorageType<4, dim, Number>::base_tensor_type &data)
1856  {
1857  switch (dim)
1858  {
1859  case 1:
1860  return data[0][0];
1861 
1862  case 2:
1863  // each entry of the tensor can be
1864  // thought of as an entry in a
1865  // matrix that maps the rolled-out
1866  // rank-2 tensors into rolled-out
1867  // rank-2 tensors. this is the
1868  // format in which we store rank-4
1869  // tensors. determine which
1870  // position the present entry is
1871  // stored in
1872  {
1873  unsigned int base_index[2];
1874  if ((indices[0] == 0) && (indices[1] == 0))
1875  base_index[0] = 0;
1876  else if ((indices[0] == 1) && (indices[1] == 1))
1877  base_index[0] = 1;
1878  else
1879  base_index[0] = 2;
1880 
1881  if ((indices[2] == 0) && (indices[3] == 0))
1882  base_index[1] = 0;
1883  else if ((indices[2] == 1) && (indices[3] == 1))
1884  base_index[1] = 1;
1885  else
1886  base_index[1] = 2;
1887 
1888  return data[base_index[0]][base_index[1]];
1889  }
1890 
1891  case 3:
1892  // each entry of the tensor can be
1893  // thought of as an entry in a
1894  // matrix that maps the rolled-out
1895  // rank-2 tensors into rolled-out
1896  // rank-2 tensors. this is the
1897  // format in which we store rank-4
1898  // tensors. determine which
1899  // position the present entry is
1900  // stored in
1901  {
1902  unsigned int base_index[2];
1903  if ((indices[0] == 0) && (indices[1] == 0))
1904  base_index[0] = 0;
1905  else if ((indices[0] == 1) && (indices[1] == 1))
1906  base_index[0] = 1;
1907  else if ((indices[0] == 2) && (indices[1] == 2))
1908  base_index[0] = 2;
1909  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1910  ((indices[0] == 1) && (indices[1] == 0)))
1911  base_index[0] = 3;
1912  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1913  ((indices[0] == 2) && (indices[1] == 0)))
1914  base_index[0] = 4;
1915  else
1916  {
1917  Assert(((indices[0] == 1) && (indices[1] == 2)) ||
1918  ((indices[0] == 2) && (indices[1] == 1)),
1919  ExcInternalError());
1920  base_index[0] = 5;
1921  }
1922 
1923  if ((indices[2] == 0) && (indices[3] == 0))
1924  base_index[1] = 0;
1925  else if ((indices[2] == 1) && (indices[3] == 1))
1926  base_index[1] = 1;
1927  else if ((indices[2] == 2) && (indices[3] == 2))
1928  base_index[1] = 2;
1929  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1930  ((indices[2] == 1) && (indices[3] == 0)))
1931  base_index[1] = 3;
1932  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1933  ((indices[2] == 2) && (indices[3] == 0)))
1934  base_index[1] = 4;
1935  else
1936  {
1937  Assert(((indices[2] == 1) && (indices[3] == 2)) ||
1938  ((indices[2] == 2) && (indices[3] == 1)),
1939  ExcInternalError());
1940  base_index[1] = 5;
1941  }
1942 
1943  return data[base_index[0]][base_index[1]];
1944  }
1945 
1946  default:
1947  Assert(false, ExcNotImplemented());
1948  }
1949 
1950  static Number dummy;
1951  return dummy;
1952  }
1953 
1954 
1955  template <int dim, typename Number>
1956  inline const Number &
1957  symmetric_tensor_access(const TableIndices<4> &indices,
1958  const typename SymmetricTensorAccessors::
1959  StorageType<4, dim, Number>::base_tensor_type &data)
1960  {
1961  switch (dim)
1962  {
1963  case 1:
1964  return data[0][0];
1965 
1966  case 2:
1967  // each entry of the tensor can be
1968  // thought of as an entry in a
1969  // matrix that maps the rolled-out
1970  // rank-2 tensors into rolled-out
1971  // rank-2 tensors. this is the
1972  // format in which we store rank-4
1973  // tensors. determine which
1974  // position the present entry is
1975  // stored in
1976  {
1977  unsigned int base_index[2];
1978  if ((indices[0] == 0) && (indices[1] == 0))
1979  base_index[0] = 0;
1980  else if ((indices[0] == 1) && (indices[1] == 1))
1981  base_index[0] = 1;
1982  else
1983  base_index[0] = 2;
1984 
1985  if ((indices[2] == 0) && (indices[3] == 0))
1986  base_index[1] = 0;
1987  else if ((indices[2] == 1) && (indices[3] == 1))
1988  base_index[1] = 1;
1989  else
1990  base_index[1] = 2;
1991 
1992  return data[base_index[0]][base_index[1]];
1993  }
1994 
1995  case 3:
1996  // each entry of the tensor can be
1997  // thought of as an entry in a
1998  // matrix that maps the rolled-out
1999  // rank-2 tensors into rolled-out
2000  // rank-2 tensors. this is the
2001  // format in which we store rank-4
2002  // tensors. determine which
2003  // position the present entry is
2004  // stored in
2005  {
2006  unsigned int base_index[2];
2007  if ((indices[0] == 0) && (indices[1] == 0))
2008  base_index[0] = 0;
2009  else if ((indices[0] == 1) && (indices[1] == 1))
2010  base_index[0] = 1;
2011  else if ((indices[0] == 2) && (indices[1] == 2))
2012  base_index[0] = 2;
2013  else if (((indices[0] == 0) && (indices[1] == 1)) ||
2014  ((indices[0] == 1) && (indices[1] == 0)))
2015  base_index[0] = 3;
2016  else if (((indices[0] == 0) && (indices[1] == 2)) ||
2017  ((indices[0] == 2) && (indices[1] == 0)))
2018  base_index[0] = 4;
2019  else
2020  {
2021  Assert(((indices[0] == 1) && (indices[1] == 2)) ||
2022  ((indices[0] == 2) && (indices[1] == 1)),
2023  ExcInternalError());
2024  base_index[0] = 5;
2025  }
2026 
2027  if ((indices[2] == 0) && (indices[3] == 0))
2028  base_index[1] = 0;
2029  else if ((indices[2] == 1) && (indices[3] == 1))
2030  base_index[1] = 1;
2031  else if ((indices[2] == 2) && (indices[3] == 2))
2032  base_index[1] = 2;
2033  else if (((indices[2] == 0) && (indices[3] == 1)) ||
2034  ((indices[2] == 1) && (indices[3] == 0)))
2035  base_index[1] = 3;
2036  else if (((indices[2] == 0) && (indices[3] == 2)) ||
2037  ((indices[2] == 2) && (indices[3] == 0)))
2038  base_index[1] = 4;
2039  else
2040  {
2041  Assert(((indices[2] == 1) && (indices[3] == 2)) ||
2042  ((indices[2] == 2) && (indices[3] == 1)),
2043  ExcInternalError());
2044  base_index[1] = 5;
2045  }
2046 
2047  return data[base_index[0]][base_index[1]];
2048  }
2049 
2050  default:
2051  Assert(false, ExcNotImplemented());
2052  }
2053 
2054  static Number dummy;
2055  return dummy;
2056  }
2057 
2058 } // end of namespace internal
2059 
2060 
2061 
2062 template <int rank_, int dim, typename Number>
2063 inline Number &
2064 SymmetricTensor<rank_, dim, Number>::
2065 operator()(const TableIndices<rank_> &indices)
2066 {
2067  for (unsigned int r = 0; r < rank; ++r)
2068  Assert(indices[r] < dimension, ExcIndexRange(indices[r], 0, dimension));
2069  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2070 }
2071 
2072 
2073 
2074 template <int rank_, int dim, typename Number>
2075 inline const Number &
2076 SymmetricTensor<rank_, dim, Number>::
2077 operator()(const TableIndices<rank_> &indices) const
2078 {
2079  for (unsigned int r = 0; r < rank; ++r)
2080  Assert(indices[r] < dimension, ExcIndexRange(indices[r], 0, dimension));
2081  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2082 }
2083 
2084 
2085 
2086 namespace internal
2087 {
2088  namespace SymmetricTensorImplementation
2089  {
2090  template <int rank_>
2091  TableIndices<rank_>
2092  get_partially_filled_indices(const unsigned int row,
2093  const std::integral_constant<int, 2> &)
2094  {
2095  return TableIndices<rank_>(row, numbers::invalid_unsigned_int);
2096  }
2097 
2098 
2099  template <int rank_>
2100  TableIndices<rank_>
2101  get_partially_filled_indices(const unsigned int row,
2102  const std::integral_constant<int, 4> &)
2103  {
2104  return TableIndices<rank_>(row,
2105  numbers::invalid_unsigned_int,
2106  numbers::invalid_unsigned_int,
2107  numbers::invalid_unsigned_int);
2108  }
2109  } // namespace SymmetricTensorImplementation
2110 } // namespace internal
2111 
2112 
2113 template <int rank_, int dim, typename Number>
2114 internal::SymmetricTensorAccessors::
2115  Accessor<rank_, dim, true, rank_ - 1, Number>
2116  SymmetricTensor<rank_, dim, Number>::
2117  operator[](const unsigned int row) const
2118 {
2119  return internal::SymmetricTensorAccessors::
2120  Accessor<rank_, dim, true, rank_ - 1, Number>(
2121  *this,
2122  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2123  rank_>(row, std::integral_constant<int, rank_>()));
2124 }
2125 
2126 
2127 
2128 template <int rank_, int dim, typename Number>
2129 internal::SymmetricTensorAccessors::
2130  Accessor<rank_, dim, false, rank_ - 1, Number>
2131  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2132 {
2133  return internal::SymmetricTensorAccessors::
2134  Accessor<rank_, dim, false, rank_ - 1, Number>(
2135  *this,
2136  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2137  rank_>(row, std::integral_constant<int, rank_>()));
2138 }
2139 
2140 
2141 
2142 template <int rank_, int dim, typename Number>
2143 inline const Number &SymmetricTensor<rank_, dim, Number>::
2144  operator[](const TableIndices<rank_> &indices) const
2145 {
2146  return operator()(indices);
2147 }
2148 
2149 
2150 
2151 template <int rank_, int dim, typename Number>
2152 inline Number &SymmetricTensor<rank_, dim, Number>::
2153  operator[](const TableIndices<rank_> &indices)
2154 {
2155  return operator()(indices);
2156 }
2157 
2158 
2159 
2160 template <int rank_, int dim, typename Number>
2161 inline Number *
2162 SymmetricTensor<rank_, dim, Number>::begin_raw()
2163 {
2164  return std::addressof(this->access_raw_entry(0));
2165 }
2166 
2167 
2168 
2169 template <int rank_, int dim, typename Number>
2170 inline const Number *
2171 SymmetricTensor<rank_, dim, Number>::begin_raw() const
2172 {
2173  return std::addressof(this->access_raw_entry(0));
2174 }
2175 
2176 
2177 
2178 template <int rank_, int dim, typename Number>
2179 inline Number *
2180 SymmetricTensor<rank_, dim, Number>::end_raw()
2181 {
2182  return begin_raw() + n_independent_components;
2183 }
2184 
2185 
2186 
2187 template <int rank_, int dim, typename Number>
2188 inline const Number *
2189 SymmetricTensor<rank_, dim, Number>::end_raw() const
2190 {
2191  return begin_raw() + n_independent_components;
2192 }
2193 
2194 
2195 
2196 namespace internal
2197 {
2198  namespace SymmetricTensorImplementation
2199  {
2200  template <int dim, typename Number>
2201  unsigned int
2202  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2203  const unsigned int index)
2204  {
2205  return index;
2206  }
2207 
2208 
2209  template <int dim, typename Number>
2210  ::TableIndices<2>
2211  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2212  const unsigned int index)
2213  {
2214  return internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
2215  base_tensor_type::unrolled_to_component_indices(index);
2216  }
2217 
2218  } // namespace SymmetricTensorImplementation
2219 } // namespace internal
2220 
2221 
2222 
2223 template <int rank_, int dim, typename Number>
2224 inline const Number &
2225 SymmetricTensor<rank_, dim, Number>::access_raw_entry(
2226  const unsigned int index) const
2227 {
2228  AssertIndexRange(index, n_independent_components);
2229  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2230  index)];
2231 }
2232 
2233 
2234 
2235 template <int rank_, int dim, typename Number>
2236 inline Number &
2237 SymmetricTensor<rank_, dim, Number>::access_raw_entry(const unsigned int index)
2238 {
2239  AssertIndexRange(index, n_independent_components);
2240  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2241  index)];
2242 }
2243 
2244 
2245 
2246 namespace internal
2247 {
2248  template <int dim, typename Number>
2249  inline typename numbers::NumberTraits<Number>::real_type
2250  compute_norm(const typename SymmetricTensorAccessors::
2251  StorageType<2, dim, Number>::base_tensor_type &data)
2252  {
2253  switch (dim)
2254  {
2255  case 1:
2256  return numbers::NumberTraits<Number>::abs(data[0]);
2257 
2258  case 2:
2259  return std::sqrt(
2260  numbers::NumberTraits<Number>::abs_square(data[0]) +
2261  numbers::NumberTraits<Number>::abs_square(data[1]) +
2262  2. * numbers::NumberTraits<Number>::abs_square(data[2]));
2263 
2264  case 3:
2265  return std::sqrt(
2266  numbers::NumberTraits<Number>::abs_square(data[0]) +
2267  numbers::NumberTraits<Number>::abs_square(data[1]) +
2268  numbers::NumberTraits<Number>::abs_square(data[2]) +
2269  2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
2270  2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
2271  2. * numbers::NumberTraits<Number>::abs_square(data[5]));
2272 
2273  default:
2274  {
2275  typename numbers::NumberTraits<Number>::real_type return_value =
2276  typename numbers::NumberTraits<Number>::real_type();
2277 
2278  for (unsigned int d = 0; d < dim; ++d)
2279  return_value +=
2280  numbers::NumberTraits<Number>::abs_square(data[d]);
2281  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2282  return_value +=
2283  2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2284 
2285  return std::sqrt(return_value);
2286  }
2287  }
2288  }
2289 
2290 
2291 
2292  template <int dim, typename Number>
2293  inline typename numbers::NumberTraits<Number>::real_type
2294  compute_norm(const typename SymmetricTensorAccessors::
2295  StorageType<4, dim, Number>::base_tensor_type &data)
2296  {
2297  switch (dim)
2298  {
2299  case 1:
2300  return numbers::NumberTraits<Number>::abs(data[0][0]);
2301 
2302  default:
2303  {
2304  typename numbers::NumberTraits<Number>::real_type return_value =
2305  typename numbers::NumberTraits<Number>::real_type();
2306 
2307  const unsigned int n_independent_components = data.dimension;
2308 
2309  for (unsigned int i = 0; i < dim; ++i)
2310  for (unsigned int j = 0; j < dim; ++j)
2311  return_value +=
2312  numbers::NumberTraits<Number>::abs_square(data[i][j]);
2313  for (unsigned int i = 0; i < dim; ++i)
2314  for (unsigned int j = dim; j < n_independent_components; ++j)
2315  return_value +=
2316  2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2317  for (unsigned int i = dim; i < n_independent_components; ++i)
2318  for (unsigned int j = 0; j < dim; ++j)
2319  return_value +=
2320  2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2321  for (unsigned int i = dim; i < n_independent_components; ++i)
2322  for (unsigned int j = dim; j < n_independent_components; ++j)
2323  return_value +=
2324  4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2325 
2326  return std::sqrt(return_value);
2327  }
2328  }
2329  }
2330 
2331 } // end of namespace internal
2332 
2333 
2334 
2335 template <int rank_, int dim, typename Number>
2336 inline typename numbers::NumberTraits<Number>::real_type
2337 SymmetricTensor<rank_, dim, Number>::norm() const
2338 {
2339  return internal::compute_norm<dim, Number>(data);
2340 }
2341 
2342 
2343 
2344 namespace internal
2345 {
2346  namespace SymmetricTensorImplementation
2347  {
2348  // a function to do the unrolling from a set of indices to a
2349  // scalar index into the array in which we store the elements of
2350  // a symmetric tensor
2351  //
2352  // this function is for rank-2 tensors
2353  template <int dim>
2354  inline unsigned int
2355  component_to_unrolled_index(const TableIndices<2> &indices)
2356  {
2357  Assert(indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2358  Assert(indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2359 
2360  switch (dim)
2361  {
2362  case 1:
2363  {
2364  return 0;
2365  }
2366 
2367  case 2:
2368  {
2369  static const unsigned int table[2][2] = {{0, 2}, {2, 1}};
2370  return table[indices[0]][indices[1]];
2371  }
2372 
2373  case 3:
2374  {
2375  static const unsigned int table[3][3] = {{0, 3, 4},
2376  {3, 1, 5},
2377  {4, 5, 2}};
2378  return table[indices[0]][indices[1]];
2379  }
2380 
2381  case 4:
2382  {
2383  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2384  {4, 1, 7, 8},
2385  {5, 7, 2, 9},
2386  {6, 8, 9, 3}};
2387  return table[indices[0]][indices[1]];
2388  }
2389 
2390  default:
2391  // for the remainder, manually figure out the numbering
2392  {
2393  if (indices[0] == indices[1])
2394  return indices[0];
2395 
2396  TableIndices<2> sorted_indices(indices);
2397  sorted_indices.sort();
2398 
2399  for (unsigned int d = 0, c = 0; d < dim; ++d)
2400  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2401  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2402  return dim + c;
2403 
2404  // should never get here:
2405  Assert(false, ExcInternalError());
2406  return 0;
2407  }
2408  }
2409  }
2410 
2411  // a function to do the unrolling from a set of indices to a
2412  // scalar index into the array in which we store the elements of
2413  // a symmetric tensor
2414  //
2415  // this function is for tensors of ranks not already handled
2416  // above
2417  template <int dim, int rank_>
2418  inline unsigned int
2419  component_to_unrolled_index(const TableIndices<rank_> &indices)
2420  {
2421  (void)indices;
2422  Assert(false, ExcNotImplemented());
2423  return numbers::invalid_unsigned_int;
2424  }
2425  } // namespace SymmetricTensorImplementation
2426 } // namespace internal
2427 
2428 
2429 template <int rank_, int dim, typename Number>
2430 inline unsigned int
2431 SymmetricTensor<rank_, dim, Number>::component_to_unrolled_index(
2432  const TableIndices<rank_> &indices)
2433 {
2434  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2435  dim>(indices);
2436 }
2437 
2438 
2439 
2440 namespace internal
2441 {
2442  namespace SymmetricTensorImplementation
2443  {
2444  // a function to do the inverse of the unrolling from a set of
2445  // indices to a scalar index into the array in which we store
2446  // the elements of a symmetric tensor. in other words, it goes
2447  // from the scalar index into the array to a set of indices of
2448  // the tensor
2449  //
2450  // this function is for rank-2 tensors
2451  template <int dim>
2452  inline TableIndices<2>
2453  unrolled_to_component_indices(const unsigned int i,
2454  const std::integral_constant<int, 2> &)
2455  {
2456  Assert(
2457  (i < ::SymmetricTensor<2, dim, double>::n_independent_components),
2458  ExcIndexRange(
2459  i,
2460  0,
2461  ::SymmetricTensor<2, dim, double>::n_independent_components));
2462  switch (dim)
2463  {
2464  case 1:
2465  {
2466  return {0, 0};
2467  }
2468 
2469  case 2:
2470  {
2471  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2472  TableIndices<2>(1, 1),
2473  TableIndices<2>(0, 1)};
2474  return table[i];
2475  }
2476 
2477  case 3:
2478  {
2479  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2480  TableIndices<2>(1, 1),
2481  TableIndices<2>(2, 2),
2482  TableIndices<2>(0, 1),
2483  TableIndices<2>(0, 2),
2484  TableIndices<2>(1, 2)};
2485  return table[i];
2486  }
2487 
2488  default:
2489  if (i < dim)
2490  return {i, i};
2491 
2492  for (unsigned int d = 0, c = 0; d < dim; ++d)
2493  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2494  if (c == i)
2495  return {d, e};
2496 
2497  // should never get here:
2498  Assert(false, ExcInternalError());
2499  return {0, 0};
2500  }
2501  }
2502 
2503  // a function to do the inverse of the unrolling from a set of
2504  // indices to a scalar index into the array in which we store
2505  // the elements of a symmetric tensor. in other words, it goes
2506  // from the scalar index into the array to a set of indices of
2507  // the tensor
2508  //
2509  // this function is for tensors of a rank not already handled
2510  // above
2511  template <int dim, int rank_>
2512  inline TableIndices<rank_>
2513  unrolled_to_component_indices(const unsigned int i,
2514  const std::integral_constant<int, rank_> &)
2515  {
2516  (void)i;
2517  Assert(
2518  (i <
2519  ::SymmetricTensor<rank_, dim, double>::n_independent_components),
2520  ExcIndexRange(i,
2521  0,
2522  ::SymmetricTensor<rank_, dim, double>::
2523  n_independent_components));
2524  Assert(false, ExcNotImplemented());
2525  return TableIndices<rank_>();
2526  }
2527 
2528  } // namespace SymmetricTensorImplementation
2529 } // namespace internal
2530 
2531 template <int rank_, int dim, typename Number>
2532 inline TableIndices<rank_>
2533 SymmetricTensor<rank_, dim, Number>::unrolled_to_component_indices(
2534  const unsigned int i)
2535 {
2536  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2537  dim>(i, std::integral_constant<int, rank_>());
2538 }
2539 
2540 
2541 
2542 template <int rank_, int dim, typename Number>
2543 template <class Archive>
2544 inline void
2545 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2546 {
2547  ar &data;
2548 }
2549 
2550 
2551 #endif // DOXYGEN
2552 
2553 /* ----------------- Non-member functions operating on tensors. ------------ */
2554 
2555 
2568 template <int rank_, int dim, typename Number, typename OtherNumber>
2569 inline SymmetricTensor<rank_,
2570  dim,
2571  typename ProductType<Number, OtherNumber>::type>
2572 operator+(const SymmetricTensor<rank_, dim, Number> & left,
2573  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2574 {
2575  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2576  tmp = left;
2577  tmp += right;
2578  return tmp;
2579 }
2580 
2581 
2594 template <int rank_, int dim, typename Number, typename OtherNumber>
2595 inline SymmetricTensor<rank_,
2596  dim,
2597  typename ProductType<Number, OtherNumber>::type>
2598 operator-(const SymmetricTensor<rank_, dim, Number> & left,
2599  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2600 {
2601  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2602  tmp = left;
2603  tmp -= right;
2604  return tmp;
2605 }
2606 
2607 
2615 template <int rank_, int dim, typename Number, typename OtherNumber>
2616 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2617 operator+(const SymmetricTensor<rank_, dim, Number> &left,
2618  const Tensor<rank_, dim, OtherNumber> & right)
2619 {
2620  return Tensor<rank_, dim, Number>(left) + right;
2621 }
2622 
2623 
2631 template <int rank_, int dim, typename Number, typename OtherNumber>
2632 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2633 operator+(const Tensor<rank_, dim, Number> & left,
2634  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2635 {
2636  return left + Tensor<rank_, dim, OtherNumber>(right);
2637 }
2638 
2639 
2647 template <int rank_, int dim, typename Number, typename OtherNumber>
2648 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2649 operator-(const SymmetricTensor<rank_, dim, Number> &left,
2650  const Tensor<rank_, dim, OtherNumber> & right)
2651 {
2652  return Tensor<rank_, dim, Number>(left) - right;
2653 }
2654 
2655 
2663 template <int rank_, int dim, typename Number, typename OtherNumber>
2664 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2665 operator-(const Tensor<rank_, dim, Number> & left,
2666  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2667 {
2668  return left - Tensor<rank_, dim, OtherNumber>(right);
2669 }
2670 
2671 
2672 
2686 template <int dim, typename Number>
2687 inline Number
2688 determinant(const SymmetricTensor<2, dim, Number> &t)
2689 {
2690  switch (dim)
2691  {
2692  case 1:
2693  return t.data[0];
2694  case 2:
2695  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2696  case 3:
2697  {
2698  // in analogy to general tensors, but
2699  // there's something to be simplified for
2700  // the present case
2701  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2702  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2703  t.data[0] * t.data[5] * t.data[5] -
2704  t.data[1] * t.data[4] * t.data[4] -
2705  t.data[2] * t.data[3] * t.data[3]);
2706  }
2707  default:
2708  Assert(false, ExcNotImplemented());
2709  return internal::NumberType<Number>::value(0.0);
2710  }
2711 }
2712 
2713 
2714 
2724 template <int dim, typename Number>
2725 inline Number
2726 third_invariant(const SymmetricTensor<2, dim, Number> &t)
2727 {
2728  return determinant(t);
2729 }
2730 
2731 
2732 
2740 template <int dim, typename Number>
2741 Number
2742 trace(const SymmetricTensor<2, dim, Number> &d)
2743 {
2744  Number t = d.data[0];
2745  for (unsigned int i = 1; i < dim; ++i)
2746  t += d.data[i];
2747  return t;
2748 }
2749 
2750 
2760 template <int dim, typename Number>
2761 inline Number
2762 first_invariant(const SymmetricTensor<2, dim, Number> &t)
2763 {
2764  return trace(t);
2765 }
2766 
2767 
2780 template <typename Number>
2781 inline Number
2782 second_invariant(const SymmetricTensor<2, 1, Number> &)
2783 {
2784  return internal::NumberType<Number>::value(0.0);
2785 }
2786 
2787 
2788 
2808 template <typename Number>
2809 inline Number
2810 second_invariant(const SymmetricTensor<2, 2, Number> &t)
2811 {
2812  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2813 }
2814 
2815 
2816 
2826 template <typename Number>
2827 inline Number
2828 second_invariant(const SymmetricTensor<2, 3, Number> &t)
2829 {
2830  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2831  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2832 }
2833 
2834 
2835 
2844 template <typename Number>
2845 std::array<Number, 1>
2846 eigenvalues(const SymmetricTensor<2, 1, Number> &T);
2847 
2848 
2849 
2873 template <typename Number>
2874 std::array<Number, 2>
2875 eigenvalues(const SymmetricTensor<2, 2, Number> &T);
2876 
2877 
2878 
2901 template <typename Number>
2902 std::array<Number, 3>
2903 eigenvalues(const SymmetricTensor<2, 3, Number> &T);
2904 
2905 
2906 
2907 namespace internal
2908 {
2909  namespace SymmetricTensorImplementation
2910  {
2950  template <int dim, typename Number>
2951  void
2952  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2953  ::Tensor<2, dim, Number> & Q,
2954  std::array<Number, dim> & d,
2955  std::array<Number, dim - 1> & e);
2956 
2957 
2958 
3000  template <int dim, typename Number>
3001  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3002  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
3003 
3004 
3005 
3047  template <int dim, typename Number>
3048  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3049  jacobi(::SymmetricTensor<2, dim, Number> A);
3050 
3051 
3052 
3068  template <typename Number>
3069  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3070  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3071 
3072 
3073 
3108  template <typename Number>
3109  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3110  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3111 
3116  template <int dim, typename Number>
3117  struct SortEigenValuesVectors
3118  {
3119  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3120  bool
3121  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3122  {
3123  return lhs.first > rhs.first;
3124  }
3125  };
3126 
3127  } // namespace SymmetricTensorImplementation
3128 
3129 } // namespace internal
3130 
3131 
3132 
3133 // The line below is to ensure that doxygen puts the full description
3134 // of this global enumeration into the documentation
3135 // See https://stackoverflow.com/a/1717984
3164 enum struct SymmetricTensorEigenvectorMethod
3165 {
3175  hybrid,
3185  ql_implicit_shifts,
3193  jacobi
3194 };
3195 
3196 
3197 
3227 template <int dim, typename Number>
3228 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3229  std::integral_constant<int, dim>::value>
3230 eigenvectors(const SymmetricTensor<2, dim, Number> &T,
3231  const SymmetricTensorEigenvectorMethod method =
3232  SymmetricTensorEigenvectorMethod::ql_implicit_shifts);
3233 
3234 
3235 
3245 template <int rank_, int dim, typename Number>
3246 inline SymmetricTensor<rank_, dim, Number>
3247 transpose(const SymmetricTensor<rank_, dim, Number> &t)
3248 {
3249  return t;
3250 }
3251 
3252 
3253 
3263 template <int dim, typename Number>
3264 inline SymmetricTensor<2, dim, Number>
3265 deviator(const SymmetricTensor<2, dim, Number> &t)
3266 {
3267  SymmetricTensor<2, dim, Number> tmp = t;
3268 
3269  // subtract scaled trace from the diagonal
3270  const Number tr = trace(t) / dim;
3271  for (unsigned int i = 0; i < dim; ++i)
3272  tmp.data[i] -= tr;
3273 
3274  return tmp;
3275 }
3276 
3277 
3278 
3286 template <int dim, typename Number>
3287 inline SymmetricTensor<2, dim, Number>
3288 unit_symmetric_tensor()
3289 {
3290  // create a default constructed matrix filled with
3291  // zeros, then set the diagonal elements to one
3292  SymmetricTensor<2, dim, Number> tmp;
3293  switch (dim)
3294  {
3295  case 1:
3296  tmp.data[0] = 1;
3297  break;
3298  case 2:
3299  tmp.data[0] = tmp.data[1] = 1;
3300  break;
3301  case 3:
3302  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3303  break;
3304  default:
3305  for (unsigned int d = 0; d < dim; ++d)
3306  tmp.data[d] = 1;
3307  }
3308  return tmp;
3309 }
3310 
3311 
3312 
3321 template <int dim>
3322 inline SymmetricTensor<2, dim>
3323 unit_symmetric_tensor()
3324 {
3325  return unit_symmetric_tensor<dim, double>();
3326 }
3327 
3328 
3329 
3344 template <int dim, typename Number>
3345 inline SymmetricTensor<4, dim, Number>
3346 deviator_tensor()
3347 {
3348  SymmetricTensor<4, dim, Number> tmp;
3349 
3350  // fill the elements treating the diagonal
3351  for (unsigned int i = 0; i < dim; ++i)
3352  for (unsigned int j = 0; j < dim; ++j)
3353  tmp.data[i][j] = (i == j ? 1 : 0) - 1. / dim;
3354 
3355  // then fill the ones that copy over the
3356  // non-diagonal elements. note that during
3357  // the double-contraction, we handle the
3358  // off-diagonal elements twice, so simply
3359  // copying requires a weight of 1/2
3360  for (unsigned int i = dim;
3361  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3362  n_rank2_components;
3363  ++i)
3364  tmp.data[i][i] = 0.5;
3365 
3366  return tmp;
3367 }
3368 
3369 
3370 
3385 template <int dim>
3386 inline SymmetricTensor<4, dim>
3387 deviator_tensor()
3388 {
3389  return deviator_tensor<dim, double>();
3390 }
3391 
3392 
3393 
3416 template <int dim, typename Number>
3417 inline SymmetricTensor<4, dim, Number>
3418 identity_tensor()
3419 {
3420  SymmetricTensor<4, dim, Number> tmp;
3421 
3422  // fill the elements treating the diagonal
3423  for (unsigned int i = 0; i < dim; ++i)
3424  tmp.data[i][i] = 1;
3425 
3426  // then fill the ones that copy over the
3427  // non-diagonal elements. note that during
3428  // the double-contraction, we handle the
3429  // off-diagonal elements twice, so simply
3430  // copying requires a weight of 1/2
3431  for (unsigned int i = dim;
3432  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3433  n_rank2_components;
3434  ++i)
3435  tmp.data[i][i] = 0.5;
3436 
3437  return tmp;
3438 }
3439 
3440 
3441 
3463 template <int dim>
3464 inline SymmetricTensor<4, dim>
3465 identity_tensor()
3466 {
3467  return identity_tensor<dim, double>();
3468 }
3469 
3470 
3471 
3482 template <int dim, typename Number>
3483 inline SymmetricTensor<2, dim, Number>
3484 invert(const SymmetricTensor<2, dim, Number> &t)
3485 {
3486  return internal::SymmetricTensorImplementation::Inverse<2, dim, Number>::
3487  value(t);
3488 }
3489 
3490 
3491 
3503 template <int dim, typename Number>
3504 inline SymmetricTensor<4, dim, Number>
3505 invert(const SymmetricTensor<4, dim, Number> &t)
3506 {
3507  return internal::SymmetricTensorImplementation::Inverse<4, dim, Number>::
3508  value(t);
3509 }
3510 
3511 
3512 
3527 template <int dim, typename Number>
3528 inline SymmetricTensor<4, dim, Number>
3529 outer_product(const SymmetricTensor<2, dim, Number> &t1,
3530  const SymmetricTensor<2, dim, Number> &t2)
3531 {
3532  SymmetricTensor<4, dim, Number> tmp;
3533 
3534  // fill only the elements really needed
3535  for (unsigned int i = 0; i < dim; ++i)
3536  for (unsigned int j = i; j < dim; ++j)
3537  for (unsigned int k = 0; k < dim; ++k)
3538  for (unsigned int l = k; l < dim; ++l)
3539  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3540 
3541  return tmp;
3542 }
3543 
3544 
3545 
3554 template <int dim, typename Number>
3555 inline SymmetricTensor<2, dim, Number>
3556 symmetrize(const Tensor<2, dim, Number> &t)
3557 {
3558  Number array[(dim * dim + dim) / 2];
3559  for (unsigned int d = 0; d < dim; ++d)
3560  array[d] = t[d][d];
3561  for (unsigned int d = 0, c = 0; d < dim; ++d)
3562  for (unsigned int e = d + 1; e < dim; ++e, ++c)
3563  array[dim + c] = (t[d][e] + t[e][d]) * 0.5;
3564  return SymmetricTensor<2, dim, Number>(array);
3565 }
3566 
3567 
3568 
3576 template <int rank_, int dim, typename Number>
3577 inline SymmetricTensor<rank_, dim, Number>
3578 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3579 {
3580  SymmetricTensor<rank_, dim, Number> tt = t;
3581  tt *= factor;
3582  return tt;
3583 }
3584 
3585 
3586 
3594 template <int rank_, int dim, typename Number>
3595 inline SymmetricTensor<rank_, dim, Number>
3596 operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3597 {
3598  // simply forward to the other operator
3599  return t * factor;
3600 }
3601 
3602 
3603 
3629 template <int rank_, int dim, typename Number, typename OtherNumber>
3630 inline SymmetricTensor<
3631  rank_,
3632  dim,
3633  typename ProductType<Number,
3634  typename EnableIfScalar<OtherNumber>::type>::type>
3635 operator*(const SymmetricTensor<rank_, dim, Number> &t,
3636  const OtherNumber & factor)
3637 {
3638  // form the product. we have to convert the two factors into the final
3639  // type via explicit casts because, for awkward reasons, the C++
3640  // standard committee saw it fit to not define an
3641  // operator*(float,std::complex<double>)
3642  // (as well as with switched arguments and double<->float).
3643  using product_type = typename ProductType<Number, OtherNumber>::type;
3644  SymmetricTensor<rank_, dim, product_type> tt(t);
3645  // we used to shorten the following by 'tt *= product_type(factor);'
3646  // which requires that a converting constructor
3647  // 'product_type::product_type(const OtherNumber) is defined.
3648  // however, a user-defined constructor is not allowed for aggregates,
3649  // e.g. VectorizedArray. therefore, we work around this issue using a
3650  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3651  // which we assume to be defined.
3652  product_type new_factor;
3653  new_factor = factor;
3654  tt *= new_factor;
3655  return tt;
3656 }
3657 
3658 
3659 
3668 template <int rank_, int dim, typename Number, typename OtherNumber>
3669 inline SymmetricTensor<
3670  rank_,
3671  dim,
3672  typename ProductType<OtherNumber,
3673  typename EnableIfScalar<Number>::type>::type>
3674 operator*(const Number & factor,
3675  const SymmetricTensor<rank_, dim, OtherNumber> &t)
3676 {
3677  // simply forward to the other operator with switched arguments
3678  return (t * factor);
3679 }
3680 
3681 
3682 
3688 template <int rank_, int dim, typename Number, typename OtherNumber>
3689 inline SymmetricTensor<
3690  rank_,
3691  dim,
3692  typename ProductType<Number,
3693  typename EnableIfScalar<OtherNumber>::type>::type>
3694 operator/(const SymmetricTensor<rank_, dim, Number> &t,
3695  const OtherNumber & factor)
3696 {
3697  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
3698  tt = t;
3699  tt /= factor;
3700  return tt;
3701 }
3702 
3703 
3704 
3711 template <int rank_, int dim>
3712 inline SymmetricTensor<rank_, dim>
3713 operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3714 {
3715  SymmetricTensor<rank_, dim> tt = t;
3716  tt *= factor;
3717  return tt;
3718 }
3719 
3720 
3721 
3728 template <int rank_, int dim>
3729 inline SymmetricTensor<rank_, dim>
3730 operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3731 {
3732  SymmetricTensor<rank_, dim> tt = t;
3733  tt *= factor;
3734  return tt;
3735 }
3736 
3737 
3738 
3744 template <int rank_, int dim>
3745 inline SymmetricTensor<rank_, dim>
3746 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3747 {
3748  SymmetricTensor<rank_, dim> tt = t;
3749  tt /= factor;
3750  return tt;
3751 }
3752 
3762 template <int dim, typename Number, typename OtherNumber>
3763 inline typename ProductType<Number, OtherNumber>::type
3764 scalar_product(const SymmetricTensor<2, dim, Number> & t1,
3765  const SymmetricTensor<2, dim, OtherNumber> &t2)
3766 {
3767  return (t1 * t2);
3768 }
3769 
3770 
3780 template <int dim, typename Number, typename OtherNumber>
3781 inline typename ProductType<Number, OtherNumber>::type
3782 scalar_product(const SymmetricTensor<2, dim, Number> &t1,
3783  const Tensor<2, dim, OtherNumber> & t2)
3784 {
3785  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<
3786  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3787  for (unsigned int i = 0; i < dim; ++i)
3788  for (unsigned int j = 0; j < dim; ++j)
3789  s += t1[i][j] * t2[i][j];
3790  return s;
3791 }
3792 
3793 
3803 template <int dim, typename Number, typename OtherNumber>
3804 inline typename ProductType<Number, OtherNumber>::type
3805 scalar_product(const Tensor<2, dim, Number> & t1,
3806  const SymmetricTensor<2, dim, OtherNumber> &t2)
3807 {
3808  return scalar_product(t2, t1);
3809 }
3810 
3811 
3827 template <typename Number, typename OtherNumber>
3828 inline void double_contract(
3829  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3830  const SymmetricTensor<4, 1, Number> & t,
3831  const SymmetricTensor<2, 1, OtherNumber> & s)
3832 {
3833  tmp[0][0] = t[0][0][0][0] * s[0][0];
3834 }
3835 
3836 
3837 
3853 template <typename Number, typename OtherNumber>
3854 inline void double_contract(
3855  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3856  const SymmetricTensor<2, 1, Number> & s,
3857  const SymmetricTensor<4, 1, OtherNumber> & t)
3858 {
3859  tmp[0][0] = t[0][0][0][0] * s[0][0];
3860 }
3861 
3862 
3863 
3878 template <typename Number, typename OtherNumber>
3879 inline void double_contract(
3880  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3881  const SymmetricTensor<4, 2, Number> & t,
3882  const SymmetricTensor<2, 2, OtherNumber> & s)
3883 {
3884  const unsigned int dim = 2;
3885 
3886  for (unsigned int i = 0; i < dim; ++i)
3887  for (unsigned int j = i; j < dim; ++j)
3888  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3889  2 * t[i][j][0][1] * s[0][1];
3890 }
3891 
3892 
3893 
3909 template <typename Number, typename OtherNumber>
3910 inline void double_contract(
3911  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3912  const SymmetricTensor<2, 2, Number> & s,
3913  const SymmetricTensor<4, 2, OtherNumber> & t)
3914 {
3915  const unsigned int dim = 2;
3916 
3917  for (unsigned int i = 0; i < dim; ++i)
3918  for (unsigned int j = i; j < dim; ++j)
3919  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3920  2 * s[0][1] * t[0][1][i][j];
3921 }
3922 
3923 
3924 
3940 template <typename Number, typename OtherNumber>
3941 inline void double_contract(
3942  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3943  const SymmetricTensor<4, 3, Number> & t,
3944  const SymmetricTensor<2, 3, OtherNumber> & s)
3945 {
3946  const unsigned int dim = 3;
3947 
3948  for (unsigned int i = 0; i < dim; ++i)
3949  for (unsigned int j = i; j < dim; ++j)
3950  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3951  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3952  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3953 }
3954 
3955 
3956 
3972 template <typename Number, typename OtherNumber>
3973 inline void double_contract(
3974  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3975  const SymmetricTensor<2, 3, Number> & s,
3976  const SymmetricTensor<4, 3, OtherNumber> & t)
3977 {
3978  const unsigned int dim = 3;
3979 
3980  for (unsigned int i = 0; i < dim; ++i)
3981  for (unsigned int j = i; j < dim; ++j)
3982  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3983  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3984  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3985 }
3986 
3987 
3988 
3996 template <int dim, typename Number, typename OtherNumber>
3997 Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
3998 operator*(const SymmetricTensor<2, dim, Number> &src1,
3999  const Tensor<1, dim, OtherNumber> & src2)
4000 {
4001  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type> dest;
4002  for (unsigned int i = 0; i < dim; ++i)
4003  for (unsigned int j = 0; j < dim; ++j)
4004  dest[i] += src1[i][j] * src2[j];
4005  return dest;
4006 }
4007 
4008 
4016 template <int dim, typename Number, typename OtherNumber>
4017 Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
4018 operator*(const Tensor<1, dim, Number> & src1,
4019  const SymmetricTensor<2, dim, OtherNumber> &src2)
4020 {
4021  // this is easy for symmetric tensors:
4022  return src2 * src1;
4023 }
4024 
4025 
4026 
4047 template <int rank_1,
4048  int rank_2,
4049  int dim,
4050  typename Number,
4051  typename OtherNumber>
4052 inline DEAL_II_ALWAYS_INLINE
4053  typename Tensor<rank_1 + rank_2 - 2,
4054  dim,
4055  typename ProductType<Number, OtherNumber>::type>::tensor_type
4056  operator*(const Tensor<rank_1, dim, Number> & src1,
4057  const SymmetricTensor<rank_2, dim, OtherNumber> &src2s)
4058 {
4059  typename Tensor<rank_1 + rank_2 - 2,
4060  dim,
4061  typename ProductType<Number, OtherNumber>::type>::tensor_type
4062  result;
4063  const Tensor<rank_2, dim, OtherNumber> src2(src2s);
4064  return src1 * src2;
4065 }
4066 
4067 
4068 
4089 template <int rank_1,
4090  int rank_2,
4091  int dim,
4092  typename Number,
4093  typename OtherNumber>
4094 inline DEAL_II_ALWAYS_INLINE
4095  typename Tensor<rank_1 + rank_2 - 2,
4096  dim,
4097  typename ProductType<Number, OtherNumber>::type>::tensor_type
4098  operator*(const SymmetricTensor<rank_1, dim, Number> &src1s,
4099  const Tensor<rank_2, dim, OtherNumber> & src2)
4100 {
4101  typename Tensor<rank_1 + rank_2 - 2,
4102  dim,
4103  typename ProductType<Number, OtherNumber>::type>::tensor_type
4104  result;
4105  const Tensor<rank_2, dim, OtherNumber> src1(src1s);
4106  return src1 * src2;
4107 }
4108 
4109 
4110 
4120 template <int dim, typename Number>
4121 inline std::ostream &
4122 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4123 {
4124  // make our lives a bit simpler by outputting
4125  // the tensor through the operator for the
4126  // general Tensor class
4127  Tensor<2, dim, Number> tt;
4128 
4129  for (unsigned int i = 0; i < dim; ++i)
4130  for (unsigned int j = 0; j < dim; ++j)
4131  tt[i][j] = t[i][j];
4132 
4133  return out << tt;
4134 }
4135 
4136 
4137 
4147 template <int dim, typename Number>
4148 inline std::ostream &
4149 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4150 {
4151  // make our lives a bit simpler by outputting
4152  // the tensor through the operator for the
4153  // general Tensor class
4154  Tensor<4, dim, Number> tt;
4155 
4156  for (unsigned int i = 0; i < dim; ++i)
4157  for (unsigned int j = 0; j < dim; ++j)
4158  for (unsigned int k = 0; k < dim; ++k)
4159  for (unsigned int l = 0; l < dim; ++l)
4160  tt[i][j][k][l] = t[i][j][k][l];
4161 
4162  return out << tt;
4163 }
4164 
4165 
4166 DEAL_II_NAMESPACE_CLOSE
4167 
4168 #endif
SymmetricTensorEigenvectorMethod::jacobi
SymmetricTensor::operator/
SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
Definition: symmetric_tensor.h:3694
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
SymmetricTensor::operator[]
internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_ - 1, Number > operator[](const unsigned int row) const
determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2688
SymmetricTensor::component_to_unrolled_index
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensorEigenvectorMethod::hybrid
SymmetricTensor::invert
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3484
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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)
SymmetricTensor::symmetrize
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3556
internal::SymmetricTensorAccessors::double_contraction_result
Definition: symmetric_tensor.h:165
trace
Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2742
numbers::value_is_zero
bool value_is_zero(const Number &value)
Definition: numbers.h:876
SymmetricTensor::double_contract
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:3828
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
deviator_tensor
SymmetricTensor< 4, dim, Number > deviator_tensor()
Definition: symmetric_tensor.h:3346
internal::SymmetricTensorAccessors::merge
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
Definition: symmetric_tensor.h:97
SymmetricTensor::operator=
SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
SymmetricTensor::memory_consumption
static std::size_t memory_consumption()
numbers::NumberTraits::abs
static real_type abs(const number &x)
Definition: numbers.h:535
SymmetricTensorEigenvectorMethod
SymmetricTensorEigenvectorMethod
Definition: symmetric_tensor.h:3164
SymmetricTensor::operator+
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:2572
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor::operator/=
SymmetricTensor & operator/=(const OtherNumber &factor)
operator*
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
Definition: symmetric_tensor.h:3578
SymmetricTensor::unrolled_to_component_indices
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
SymmetricTensor::operator-
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:2598
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor::first_invariant
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2762
SymmetricTensor::access_raw_entry
const Number & access_raw_entry(const unsigned int unrolled_index) const
SymmetricTensor< 2, spacedim >::base_tensor_type
typename base_tensor_descriptor::base_tensor_type base_tensor_type
Definition: symmetric_tensor.h:887
SymmetricTensor::operator!=
bool operator!=(const SymmetricTensor &) const
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
SymmetricTensor::begin_raw
Number * begin_raw()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1407
SymmetricTensor::data
base_tensor_type data
Definition: symmetric_tensor.h:892
SymmetricTensor::deviator
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3265
SymmetricTensor::trace
Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2742
SymmetricTensor::serialize
void serialize(Archive &ar, const unsigned int version)
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)
SymmetricTensor::end_raw
Number * end_raw()
SymmetricTensor::outer_product
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition: symmetric_tensor.h:3529
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor::scalar_product
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Definition: symmetric_tensor.h:3764
SymmetricTensorEigenvectorMethod::ql_implicit_shifts
SymmetricTensor::operator()
Number & operator()(const TableIndices< rank_ > &indices)
swap
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
SymmetricTensor::determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2688
deviator
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3265
SymmetricTensor::operator-
SymmetricTensor operator-() const
Tensor
Definition: mpi.h:90
SymmetricTensor::operator*=
SymmetricTensor & operator*=(const OtherNumber &factor)
SymmetricTensor::operator*
internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
SymmetricTensor::operator+=
SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
SymmetricTensor::third_invariant
Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2726
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
SymmetricTensor::operator==
bool operator==(const SymmetricTensor &) const
SymmetricTensor::eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
unit_symmetric_tensor
SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: symmetric_tensor.h:3288
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:468
operator-
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:2598
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
identity_tensor
SymmetricTensor< 4, dim, Number > identity_tensor()
Definition: symmetric_tensor.h:3418
SymmetricTensor::second_invariant
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Definition: symmetric_tensor.h:2782
SymmetricTensor::operator-=
SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
SymmetricTensor::norm
numbers::NumberTraits< Number >::real_type norm() const
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Point::operator*
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
Generated by   doxygen 1.8.14