Reference documentation for deal.II version 9.1.0-pre
symmetric_tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2018 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 TableIndices<2>(new_index, numbers::invalid_unsigned_int);
105  else
106  return TableIndices<2>(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 TableIndices<4>(new_index,
128  numbers::invalid_unsigned_int,
129  numbers::invalid_unsigned_int,
130  numbers::invalid_unsigned_int);
131  case 1:
132  return TableIndices<4>(previous_indices[0],
133  new_index,
134  numbers::invalid_unsigned_int,
135  numbers::invalid_unsigned_int);
136  case 2:
137  return TableIndices<4>(previous_indices[0],
138  previous_indices[1],
139  new_index,
140  numbers::invalid_unsigned_int);
141  case 3:
142  return TableIndices<4>(previous_indices[0],
143  previous_indices[1],
144  previous_indices[2],
145  new_index);
146  }
147  Assert(false, ExcInternalError());
148  return TableIndices<4>();
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() = delete;
449 
453  Accessor(const Accessor &) = default;
454 
455  public:
459  reference operator[](const unsigned int);
460 
464  reference operator[](const unsigned int) const;
465 
466  private:
470  tensor_type & tensor;
471  const TableIndices<rank> previous_indices;
472 
473  // declare some other classes
474  // as friends. make sure to
475  // work around bugs in some
476  // compilers
477  template <int, int, typename>
478  friend class ::SymmetricTensor;
479  template <int, int, bool, int, typename>
480  friend class SymmetricTensorAccessors::Accessor;
481 #ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
482  friend class ::SymmetricTensor<rank, dim, Number>;
483  friend class SymmetricTensorAccessors::
484  Accessor<rank, dim, constness, 2, Number>;
485 #endif
486  };
487  } // namespace SymmetricTensorAccessors
488 } // namespace internal
489 
490 
491 
555 template <int rank_, int dim, typename Number>
556 class SymmetricTensor
557 {
558 public:
559  static_assert(rank_ % 2 == 0, "A SymmetricTensor must have even rank!");
560 
569  static const unsigned int dimension = dim;
570 
574  static const unsigned int rank = rank_;
575 
581  static constexpr unsigned int n_independent_components =
582  internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>::
583  n_independent_components;
584 
588  SymmetricTensor();
589 
600  template <typename OtherNumber>
601  explicit SymmetricTensor(const Tensor<2, dim, OtherNumber> &t);
602 
618  SymmetricTensor(const Number (&array)[n_independent_components]);
619 
625  template <typename OtherNumber>
626  explicit SymmetricTensor(
627  const SymmetricTensor<rank_, dim, OtherNumber> &initializer);
628 
632  Number *
633  begin_raw();
634 
638  const Number *
639  begin_raw() const;
640 
644  Number *
645  end_raw();
646 
651  const Number *
652  end_raw() const;
653 
660  template <typename OtherNumber>
661  SymmetricTensor &
662  operator=(const SymmetricTensor<rank_, dim, OtherNumber> &rhs);
663 
670  SymmetricTensor &
671  operator=(const Number &d);
672 
677  operator Tensor<rank_, dim, Number>() const;
678 
682  bool
683  operator==(const SymmetricTensor &) const;
684 
688  bool
689  operator!=(const SymmetricTensor &) const;
690 
694  template <typename OtherNumber>
695  SymmetricTensor &
696  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &);
697 
701  template <typename OtherNumber>
702  SymmetricTensor &
703  operator-=(const SymmetricTensor<rank_, dim, OtherNumber> &);
704 
709  template <typename OtherNumber>
710  SymmetricTensor &
711  operator*=(const OtherNumber &factor);
712 
716  template <typename OtherNumber>
717  SymmetricTensor &
718  operator/=(const OtherNumber &factor);
719 
723  SymmetricTensor
724  operator-() const;
725 
750  template <typename OtherNumber>
751  typename internal::SymmetricTensorAccessors::
752  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
753  operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const;
754 
759  template <typename OtherNumber>
760  typename internal::SymmetricTensorAccessors::
761  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
762  operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const;
763 
767  Number &
768  operator()(const TableIndices<rank_> &indices);
769 
773  const Number &
774  operator()(const TableIndices<rank_> &indices) const;
775 
780  internal::SymmetricTensorAccessors::
781  Accessor<rank_, dim, true, rank_ - 1, Number>
782  operator[](const unsigned int row) const;
783 
788  internal::SymmetricTensorAccessors::
789  Accessor<rank_, dim, false, rank_ - 1, Number>
790  operator[](const unsigned int row);
791 
797  const Number &operator[](const TableIndices<rank_> &indices) const;
798 
804  Number &operator[](const TableIndices<rank_> &indices);
805 
811  const Number &
812  access_raw_entry(const unsigned int unrolled_index) const;
813 
819  Number &
820  access_raw_entry(const unsigned int unrolled_index);
821 
831  typename numbers::NumberTraits<Number>::real_type
832  norm() const;
833 
841  static unsigned int
842  component_to_unrolled_index(const TableIndices<rank_> &indices);
843 
849  static TableIndices<rank_>
850  unrolled_to_component_indices(const unsigned int i);
851 
864  void
865  clear();
866 
871  static std::size_t
872  memory_consumption();
873 
878  template <class Archive>
879  void
880  serialize(Archive &ar, const unsigned int version);
881 
882 private:
886  using base_tensor_descriptor =
887  internal::SymmetricTensorAccessors::StorageType<rank_, dim, Number>;
888 
892  using base_tensor_type = typename base_tensor_descriptor::base_tensor_type;
893 
897  base_tensor_type data;
898 
902  template <int, int, typename>
903  friend class SymmetricTensor;
904 
908  template <int dim2, typename Number2>
909  friend Number2
910  trace(const SymmetricTensor<2, dim2, Number2> &d);
911 
912  template <int dim2, typename Number2>
913  friend Number2
914  determinant(const SymmetricTensor<2, dim2, Number2> &t);
915 
916  template <int dim2, typename Number2>
917  friend SymmetricTensor<2, dim2, Number2>
918  deviator(const SymmetricTensor<2, dim2, Number2> &t);
919 
920  template <int dim2, typename Number2>
921  friend SymmetricTensor<2, dim2, Number2>
922  unit_symmetric_tensor();
923 
924  template <int dim2, typename Number2>
925  friend SymmetricTensor<4, dim2, Number2>
926  deviator_tensor();
927 
928  template <int dim2, typename Number2>
929  friend SymmetricTensor<4, dim2, Number2>
930  identity_tensor();
931 
932 
936  friend struct internal::SymmetricTensorImplementation::
937  Inverse<2, dim, Number>;
938 
939  friend struct internal::SymmetricTensorImplementation::
940  Inverse<4, dim, Number>;
941 };
942 
943 
944 
945 // ------------------------- inline functions ------------------------
946 
947 #ifndef DOXYGEN
948 
949 // provide declarations for static members
950 template <int rank, int dim, typename Number>
951 const unsigned int SymmetricTensor<rank, dim, Number>::dimension;
952 
953 template <int rank_, int dim, typename Number>
954 constexpr unsigned int
955  SymmetricTensor<rank_, dim, Number>::n_independent_components;
956 
957 namespace internal
958 {
959  namespace SymmetricTensorAccessors
960  {
961  template <int rank_, int dim, bool constness, int P, typename Number>
962  Accessor<rank_, dim, constness, P, Number>::Accessor(
963  tensor_type & tensor,
964  const TableIndices<rank_> &previous_indices)
965  : tensor(tensor)
966  , previous_indices(previous_indices)
967  {}
968 
969 
970 
971  template <int rank_, int dim, bool constness, int P, typename Number>
972  Accessor<rank_, dim, constness, P - 1, Number>
973  Accessor<rank_, dim, constness, P, Number>::
974  operator[](const unsigned int i)
975  {
976  return Accessor<rank_, dim, constness, P - 1, Number>(
977  tensor, merge(previous_indices, i, rank_ - P));
978  }
979 
980 
981 
982  template <int rank_, int dim, bool constness, int P, typename Number>
983  Accessor<rank_, dim, constness, P - 1, Number>
984  Accessor<rank_, dim, constness, P, Number>::
985  operator[](const unsigned int i) const
986  {
987  return Accessor<rank_, dim, constness, P - 1, Number>(
988  tensor, merge(previous_indices, i, rank_ - P));
989  }
990 
991 
992 
993  template <int rank_, int dim, bool constness, typename Number>
994  Accessor<rank_, dim, constness, 1, Number>::Accessor(
995  tensor_type & tensor,
996  const TableIndices<rank_> &previous_indices)
997  : tensor(tensor)
998  , previous_indices(previous_indices)
999  {}
1000 
1001 
1002 
1003  template <int rank_, int dim, bool constness, typename Number>
1004  typename Accessor<rank_, dim, constness, 1, Number>::reference
1005  Accessor<rank_, dim, constness, 1, Number>::
1006  operator[](const unsigned int i)
1007  {
1008  return tensor(merge(previous_indices, i, rank_ - 1));
1009  }
1010 
1011 
1012  template <int rank_, int dim, bool constness, typename Number>
1013  typename Accessor<rank_, dim, constness, 1, Number>::reference
1014  Accessor<rank_, dim, constness, 1, Number>::
1015  operator[](const unsigned int i) const
1016  {
1017  return tensor(merge(previous_indices, i, rank_ - 1));
1018  }
1019  } // namespace SymmetricTensorAccessors
1020 } // namespace internal
1021 
1022 
1023 
1024 template <int rank_, int dim, typename Number>
1025 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor()
1026 {
1027  // Some auto-differentiable numbers need explicit
1028  // zero initialization.
1029  for (unsigned int i = 0; i < base_tensor_type::dimension; ++i)
1030  data[i] = internal::NumberType<Number>::value(0.0);
1031 }
1032 
1033 
1034 template <int rank_, int dim, typename Number>
1035 template <typename OtherNumber>
1036 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1037  const Tensor<2, dim, OtherNumber> &t)
1038 {
1039  Assert(rank == 2, ExcNotImplemented());
1040  switch (dim)
1041  {
1042  case 2:
1043  Assert(t[0][1] == t[1][0], ExcInternalError());
1044 
1045  data[0] = t[0][0];
1046  data[1] = t[1][1];
1047  data[2] = t[0][1];
1048 
1049  break;
1050  case 3:
1051  Assert(t[0][1] == t[1][0], ExcInternalError());
1052  Assert(t[0][2] == t[2][0], ExcInternalError());
1053  Assert(t[1][2] == t[2][1], ExcInternalError());
1054 
1055  data[0] = t[0][0];
1056  data[1] = t[1][1];
1057  data[2] = t[2][2];
1058  data[3] = t[0][1];
1059  data[4] = t[0][2];
1060  data[5] = t[1][2];
1061 
1062  break;
1063  default:
1064  for (unsigned int d = 0; d < dim; ++d)
1065  for (unsigned int e = 0; e < d; ++e)
1066  Assert(t[d][e] == t[e][d], ExcInternalError());
1067 
1068  for (unsigned int d = 0; d < dim; ++d)
1069  data[d] = t[d][d];
1070 
1071  for (unsigned int d = 0, c = 0; d < dim; ++d)
1072  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1073  data[dim + c] = t[d][e];
1074  }
1075 }
1076 
1077 
1078 
1079 template <int rank_, int dim, typename Number>
1080 template <typename OtherNumber>
1081 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1082  const SymmetricTensor<rank_, dim, OtherNumber> &initializer)
1083 {
1084  for (unsigned int i = 0; i < base_tensor_type::dimension; ++i)
1085  data[i] =
1086  internal::NumberType<typename base_tensor_type::value_type>::value(
1087  initializer.data[i]);
1088 }
1089 
1090 
1091 
1092 template <int rank_, int dim, typename Number>
1093 inline SymmetricTensor<rank_, dim, Number>::SymmetricTensor(
1094  const Number (&array)[n_independent_components])
1095  : data(
1096  *reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1097 {
1098  // ensure that the reinterpret_cast above actually works
1099  Assert(sizeof(typename base_tensor_type::array_type) == sizeof(array),
1100  ExcInternalError());
1101 }
1102 
1103 
1104 
1105 template <int rank_, int dim, typename Number>
1106 template <typename OtherNumber>
1107 inline SymmetricTensor<rank_, dim, Number> &
1108 SymmetricTensor<rank_, dim, Number>::
1109 operator=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1110 {
1111  for (unsigned int i = 0; i < base_tensor_type::dimension; ++i)
1112  data[i] = t.data[i];
1113  return *this;
1114 }
1115 
1116 
1117 
1118 template <int rank_, int dim, typename Number>
1119 inline SymmetricTensor<rank_, dim, Number> &
1120 SymmetricTensor<rank_, dim, Number>::operator=(const Number &d)
1121 {
1122  Assert(numbers::value_is_zero(d),
1123  ExcMessage("Only assignment with zero is allowed"));
1124  (void)d;
1125 
1126  data = internal::NumberType<Number>::value(0.0);
1127 
1128  return *this;
1129 }
1130 
1131 
1132 namespace internal
1133 {
1134  namespace SymmetricTensorImplementation
1135  {
1136  template <int dim, typename Number>
1137  inline DEAL_II_ALWAYS_INLINE ::Tensor<2, dim, Number>
1138  convert_to_tensor(const ::SymmetricTensor<2, dim, Number> &s)
1139  {
1140  ::Tensor<2, dim, Number> t;
1141 
1142  // diagonal entries are stored first
1143  for (unsigned int d = 0; d < dim; ++d)
1144  t[d][d] = s.access_raw_entry(d);
1145 
1146  // off-diagonal entries come next, row by row
1147  for (unsigned int d = 0, c = 0; d < dim; ++d)
1148  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1149  {
1150  t[d][e] = s.access_raw_entry(dim + c);
1151  t[e][d] = s.access_raw_entry(dim + c);
1152  }
1153  return t;
1154  }
1155 
1156 
1157  template <int dim, typename Number>
1158  ::Tensor<4, dim, Number>
1159  convert_to_tensor(const ::SymmetricTensor<4, dim, Number> &st)
1160  {
1161  // utilize the symmetry properties of SymmetricTensor<4,dim>
1162  // discussed in the class documentation to avoid accessing all
1163  // independent elements of the input tensor more than once
1164  ::Tensor<4, dim, Number> t;
1165 
1166  for (unsigned int i = 0; i < dim; ++i)
1167  for (unsigned int j = i; j < dim; ++j)
1168  for (unsigned int k = 0; k < dim; ++k)
1169  for (unsigned int l = k; l < dim; ++l)
1170  t[TableIndices<4>(i, j, k, l)] = t[TableIndices<4>(i, j, l, k)] =
1171  t[TableIndices<4>(j, i, k, l)] =
1172  t[TableIndices<4>(j, i, l, k)] =
1173  st[TableIndices<4>(i, j, k, l)];
1174 
1175  return t;
1176  }
1177 
1178 
1179  template <typename Number>
1180  struct Inverse<2, 1, Number>
1181  {
1182  static inline ::SymmetricTensor<2, 1, Number>
1183  value(const ::SymmetricTensor<2, 1, Number> &t)
1184  {
1185  ::SymmetricTensor<2, 1, Number> tmp;
1186 
1187  tmp[0][0] = 1.0 / t[0][0];
1188 
1189  return tmp;
1190  }
1191  };
1192 
1193 
1194  template <typename Number>
1195  struct Inverse<2, 2, Number>
1196  {
1197  static inline ::SymmetricTensor<2, 2, Number>
1198  value(const ::SymmetricTensor<2, 2, Number> &t)
1199  {
1200  ::SymmetricTensor<2, 2, Number> tmp;
1201 
1202  // Sympy result: ([
1203  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1204  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1205  const TableIndices<2> idx_00(0, 0);
1206  const TableIndices<2> idx_01(0, 1);
1207  const TableIndices<2> idx_11(1, 1);
1208  const Number inv_det_t =
1209  1.0 / (t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01]);
1210  tmp[idx_00] = t[idx_11];
1211  tmp[idx_01] = -t[idx_01];
1212  tmp[idx_11] = t[idx_00];
1213  tmp *= inv_det_t;
1214 
1215  return tmp;
1216  }
1217  };
1218 
1219 
1220  template <typename Number>
1221  struct Inverse<2, 3, Number>
1222  {
1223  static ::SymmetricTensor<2, 3, Number>
1224  value(const ::SymmetricTensor<2, 3, Number> &t)
1225  {
1226  ::SymmetricTensor<2, 3, Number> tmp;
1227 
1228  // Sympy result: ([
1229  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1230  // 2*t01*t02*t12 - t02**2*t11),
1231  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1232  // 2*t01*t02*t12 - t02**2*t11), (t01*t12 - t02*t11)/(t00*t11*t22 -
1233  // t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
1234  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1235  // 2*t01*t02*t12 - t02**2*t11),
1236  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1237  // 2*t01*t02*t12 - t02**2*t11), (t00*t12 - t01*t02)/(-t00*t11*t22 +
1238  // t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
1239  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 +
1240  // 2*t01*t02*t12 - t02**2*t11),
1241  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1242  // 2*t01*t02*t12 + t02**2*t11),
1243  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 -
1244  // 2*t01*t02*t12 + t02**2*t11)] ])
1245  const TableIndices<2> idx_00(0, 0);
1246  const TableIndices<2> idx_01(0, 1);
1247  const TableIndices<2> idx_02(0, 2);
1248  const TableIndices<2> idx_11(1, 1);
1249  const TableIndices<2> idx_12(1, 2);
1250  const TableIndices<2> idx_22(2, 2);
1251  const Number inv_det_t =
1252  1.0 / (t[idx_00] * t[idx_11] * t[idx_22] -
1253  t[idx_00] * t[idx_12] * t[idx_12] -
1254  t[idx_01] * t[idx_01] * t[idx_22] +
1255  2.0 * t[idx_01] * t[idx_02] * t[idx_12] -
1256  t[idx_02] * t[idx_02] * t[idx_11]);
1257  tmp[idx_00] = t[idx_11] * t[idx_22] - t[idx_12] * t[idx_12];
1258  tmp[idx_01] = -t[idx_01] * t[idx_22] + t[idx_02] * t[idx_12];
1259  tmp[idx_02] = t[idx_01] * t[idx_12] - t[idx_02] * t[idx_11];
1260  tmp[idx_11] = t[idx_00] * t[idx_22] - t[idx_02] * t[idx_02];
1261  tmp[idx_12] = -t[idx_00] * t[idx_12] + t[idx_01] * t[idx_02];
1262  tmp[idx_22] = t[idx_00] * t[idx_11] - t[idx_01] * t[idx_01];
1263  tmp *= inv_det_t;
1264 
1265  return tmp;
1266  }
1267  };
1268 
1269 
1270  template <typename Number>
1271  struct Inverse<4, 1, Number>
1272  {
1273  static inline ::SymmetricTensor<4, 1, Number>
1274  value(const ::SymmetricTensor<4, 1, Number> &t)
1275  {
1276  ::SymmetricTensor<4, 1, Number> tmp;
1277  tmp.data[0][0] = 1.0 / t.data[0][0];
1278  return tmp;
1279  }
1280  };
1281 
1282 
1283  template <typename Number>
1284  struct Inverse<4, 2, Number>
1285  {
1286  static inline ::SymmetricTensor<4, 2, Number>
1287  value(const ::SymmetricTensor<4, 2, Number> &t)
1288  {
1289  ::SymmetricTensor<4, 2, Number> tmp;
1290 
1291  // Inverting this tensor is a little more complicated than necessary,
1292  // since we store the data of 't' as a 3x3 matrix t.data, but the
1293  // product between a rank-4 and a rank-2 tensor is really not the
1294  // product between this matrix and the 3-vector of a rhs, but rather
1295  //
1296  // B.vec = t.data * mult * A.vec
1297  //
1298  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1299  // capture the fact that we need to add up both the c_ij12*a_12 and the
1300  // c_ij21*a_21 terms.
1301  //
1302  // In addition, in this scheme, the identity tensor has the matrix
1303  // representation mult^-1.
1304  //
1305  // The inverse of 't' therefore has the matrix representation
1306  //
1307  // inv.data = mult^-1 * t.data^-1 * mult^-1
1308  //
1309  // in order to compute it, let's first compute the inverse of t.data and
1310  // put it into tmp.data; at the end of the function we then scale the
1311  // last row and column of the inverse by 1/2, corresponding to the left
1312  // and right multiplication with mult^-1.
1313  const Number t4 = t.data[0][0] * t.data[1][1],
1314  t6 = t.data[0][0] * t.data[1][2],
1315  t8 = t.data[0][1] * t.data[1][0],
1316  t00 = t.data[0][2] * t.data[1][0],
1317  t01 = t.data[0][1] * t.data[2][0],
1318  t04 = t.data[0][2] * t.data[2][0],
1319  t07 = 1.0 / (t4 * t.data[2][2] - t6 * t.data[2][1] -
1320  t8 * t.data[2][2] + t00 * t.data[2][1] +
1321  t01 * t.data[1][2] - t04 * t.data[1][1]);
1322  tmp.data[0][0] =
1323  (t.data[1][1] * t.data[2][2] - t.data[1][2] * t.data[2][1]) * t07;
1324  tmp.data[0][1] =
1325  -(t.data[0][1] * t.data[2][2] - t.data[0][2] * t.data[2][1]) * t07;
1326  tmp.data[0][2] =
1327  -(-t.data[0][1] * t.data[1][2] + t.data[0][2] * t.data[1][1]) * t07;
1328  tmp.data[1][0] =
1329  -(t.data[1][0] * t.data[2][2] - t.data[1][2] * t.data[2][0]) * t07;
1330  tmp.data[1][1] = (t.data[0][0] * t.data[2][2] - t04) * t07;
1331  tmp.data[1][2] = -(t6 - t00) * t07;
1332  tmp.data[2][0] =
1333  -(-t.data[1][0] * t.data[2][1] + t.data[1][1] * t.data[2][0]) * t07;
1334  tmp.data[2][1] = -(t.data[0][0] * t.data[2][1] - t01) * t07;
1335  tmp.data[2][2] = (t4 - t8) * t07;
1336 
1337  // scale last row and column as mentioned
1338  // above
1339  tmp.data[2][0] /= 2;
1340  tmp.data[2][1] /= 2;
1341  tmp.data[0][2] /= 2;
1342  tmp.data[1][2] /= 2;
1343  tmp.data[2][2] /= 4;
1344 
1345  return tmp;
1346  }
1347  };
1348 
1349 
1350  template <typename Number>
1351  struct Inverse<4, 3, Number>
1352  {
1353  static ::SymmetricTensor<4, 3, Number>
1354  value(const ::SymmetricTensor<4, 3, Number> &t)
1355  {
1356  ::SymmetricTensor<4, 3, Number> tmp = t;
1357 
1358  // This function follows the exact same scheme as the 2d case, except
1359  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1360  // Instead, we use the Gauss-Jordan algorithm implemented for
1361  // FullMatrix. For historical reasons the following code is copied from
1362  // there, with the tangential benefit that we do not need to copy the
1363  // tensor entries to and from the FullMatrix.
1364  const unsigned int N = 6;
1365 
1366  // First get an estimate of the size of the elements of this matrix,
1367  // for later checks whether the pivot element is large enough, or
1368  // whether we have to fear that the matrix is not regular.
1369  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1370  for (unsigned int i = 0; i < N; ++i)
1371  diagonal_sum += std::fabs(tmp.data[i][i]);
1372  const Number typical_diagonal_element =
1373  diagonal_sum / static_cast<double>(N);
1374  (void)typical_diagonal_element;
1375 
1376  unsigned int p[N];
1377  for (unsigned int i = 0; i < N; ++i)
1378  p[i] = i;
1379 
1380  for (unsigned int j = 0; j < N; ++j)
1381  {
1382  // Pivot search: search that part of the line on and right of the
1383  // diagonal for the largest element.
1384  Number max = std::fabs(tmp.data[j][j]);
1385  unsigned int r = j;
1386  for (unsigned int i = j + 1; i < N; ++i)
1387  if (std::fabs(tmp.data[i][j]) > max)
1388  {
1389  max = std::fabs(tmp.data[i][j]);
1390  r = i;
1391  }
1392 
1393  // Check whether the pivot is too small
1394  Assert(max > 1.e-16 * typical_diagonal_element,
1395  ExcMessage("This tensor seems to be noninvertible"));
1396 
1397  // Row interchange
1398  if (r > j)
1399  {
1400  for (unsigned int k = 0; k < N; ++k)
1401  std::swap(tmp.data[j][k], tmp.data[r][k]);
1402 
1403  std::swap(p[j], p[r]);
1404  }
1405 
1406  // Transformation
1407  const Number hr = 1. / tmp.data[j][j];
1408  tmp.data[j][j] = hr;
1409  for (unsigned int k = 0; k < N; ++k)
1410  {
1411  if (k == j)
1412  continue;
1413  for (unsigned int i = 0; i < N; ++i)
1414  {
1415  if (i == j)
1416  continue;
1417  tmp.data[i][k] -= tmp.data[i][j] * tmp.data[j][k] * hr;
1418  }
1419  }
1420  for (unsigned int i = 0; i < N; ++i)
1421  {
1422  tmp.data[i][j] *= hr;
1423  tmp.data[j][i] *= -hr;
1424  }
1425  tmp.data[j][j] = hr;
1426  }
1427 
1428  // Column interchange
1429  Number hv[N];
1430  for (unsigned int i = 0; i < N; ++i)
1431  {
1432  for (unsigned int k = 0; k < N; ++k)
1433  hv[p[k]] = tmp.data[i][k];
1434  for (unsigned int k = 0; k < N; ++k)
1435  tmp.data[i][k] = hv[k];
1436  }
1437 
1438  // Scale rows and columns. The mult matrix
1439  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1440  for (unsigned int i = 3; i < 6; ++i)
1441  for (unsigned int j = 0; j < 3; ++j)
1442  tmp.data[i][j] /= 2;
1443 
1444  for (unsigned int i = 0; i < 3; ++i)
1445  for (unsigned int j = 3; j < 6; ++j)
1446  tmp.data[i][j] /= 2;
1447 
1448  for (unsigned int i = 3; i < 6; ++i)
1449  for (unsigned int j = 3; j < 6; ++j)
1450  tmp.data[i][j] /= 4;
1451 
1452  return tmp;
1453  }
1454  };
1455 
1456  } // namespace SymmetricTensorImplementation
1457 } // namespace internal
1458 
1459 
1460 
1461 template <int rank_, int dim, typename Number>
1462 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>::
1463  operator Tensor<rank_, dim, Number>() const
1464 {
1465  return internal::SymmetricTensorImplementation::convert_to_tensor(*this);
1466 }
1467 
1468 
1469 
1470 template <int rank_, int dim, typename Number>
1471 inline bool
1472 SymmetricTensor<rank_, dim, Number>::
1473 operator==(const SymmetricTensor<rank_, dim, Number> &t) const
1474 {
1475  return data == t.data;
1476 }
1477 
1478 
1479 
1480 template <int rank_, int dim, typename Number>
1481 inline bool
1482 SymmetricTensor<rank_, dim, Number>::
1483 operator!=(const SymmetricTensor<rank_, dim, Number> &t) const
1484 {
1485  return data != t.data;
1486 }
1487 
1488 
1489 
1490 template <int rank_, int dim, typename Number>
1491 template <typename OtherNumber>
1492 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
1493  SymmetricTensor<rank_, dim, Number>::
1494  operator+=(const SymmetricTensor<rank_, dim, OtherNumber> &t)
1495 {
1496  data += t.data;
1497  return *this;
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>::operator*=(const OtherNumber &d)
1518 {
1519  data *= d;
1520  return *this;
1521 }
1522 
1523 
1524 
1525 template <int rank_, int dim, typename Number>
1526 template <typename OtherNumber>
1527 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number> &
1528 SymmetricTensor<rank_, dim, Number>::operator/=(const OtherNumber &d)
1529 {
1530  data /= d;
1531  return *this;
1532 }
1533 
1534 
1535 
1536 template <int rank_, int dim, typename Number>
1537 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<rank_, dim, Number>
1538 SymmetricTensor<rank_, dim, Number>::operator-() const
1539 {
1540  SymmetricTensor tmp = *this;
1541  tmp.data = -tmp.data;
1542  return tmp;
1543 }
1544 
1545 
1546 
1547 template <int rank_, int dim, typename Number>
1548 inline DEAL_II_ALWAYS_INLINE void
1549 SymmetricTensor<rank_, dim, Number>::clear()
1550 {
1551  data.clear();
1552 }
1553 
1554 
1555 
1556 template <int rank_, int dim, typename Number>
1557 inline std::size_t
1558 SymmetricTensor<rank_, dim, Number>::memory_consumption()
1559 {
1560  // all memory consists of statically allocated memory of the current
1561  // object, no pointers
1562  return sizeof(SymmetricTensor<rank_, dim, Number>);
1563 }
1564 
1565 
1566 
1567 namespace internal
1568 {
1569  template <int dim, typename Number, typename OtherNumber = Number>
1570  inline DEAL_II_ALWAYS_INLINE typename SymmetricTensorAccessors::
1571  double_contraction_result<2, 2, dim, Number, OtherNumber>::type
1572  perform_double_contraction(
1573  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1574  base_tensor_type &data,
1575  const typename SymmetricTensorAccessors::
1576  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1577  {
1578  using result_type = typename SymmetricTensorAccessors::
1579  double_contraction_result<2, 2, dim, Number, OtherNumber>::type;
1580 
1581  switch (dim)
1582  {
1583  case 1:
1584  return data[0] * sdata[0];
1585  default:
1586  // Start with the non-diagonal part to avoid some multiplications by
1587  // 2.
1588 
1589  result_type sum = data[dim] * sdata[dim];
1590  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1591  sum += data[d] * sdata[d];
1592  sum += sum; // sum = sum * 2.;
1593 
1594  // Now add the contributions from the diagonal
1595  for (unsigned int d = 0; d < dim; ++d)
1596  sum += data[d] * sdata[d];
1597  return sum;
1598  }
1599  }
1600 
1601 
1602 
1603  template <int dim, typename Number, typename OtherNumber = Number>
1604  inline typename SymmetricTensorAccessors::
1605  double_contraction_result<4, 2, dim, Number, OtherNumber>::type
1606  perform_double_contraction(
1607  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1608  base_tensor_type &data,
1609  const typename SymmetricTensorAccessors::
1610  StorageType<2, dim, OtherNumber>::base_tensor_type &sdata)
1611  {
1612  using result_type = typename SymmetricTensorAccessors::
1613  double_contraction_result<4, 2, dim, Number, OtherNumber>::type;
1614  using value_type = typename SymmetricTensorAccessors::
1615  double_contraction_result<4, 2, dim, Number, OtherNumber>::value_type;
1616 
1617  const unsigned int data_dim = SymmetricTensorAccessors::
1618  StorageType<2, dim, value_type>::n_independent_components;
1619  value_type tmp[data_dim];
1620  for (unsigned int i = 0; i < data_dim; ++i)
1621  tmp[i] =
1622  perform_double_contraction<dim, Number, OtherNumber>(data[i], sdata);
1623  return result_type(tmp);
1624  }
1625 
1626 
1627 
1628  template <int dim, typename Number, typename OtherNumber = Number>
1629  inline typename SymmetricTensorAccessors::StorageType<
1630  2,
1631  dim,
1632  typename SymmetricTensorAccessors::
1633  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type>::
1634  base_tensor_type
1635  perform_double_contraction(
1636  const typename SymmetricTensorAccessors::StorageType<2, dim, Number>::
1637  base_tensor_type &data,
1638  const typename SymmetricTensorAccessors::
1639  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1640  {
1641  using value_type = typename SymmetricTensorAccessors::
1642  double_contraction_result<2, 4, dim, Number, OtherNumber>::value_type;
1643  using base_tensor_type = typename SymmetricTensorAccessors::
1644  StorageType<2, dim, value_type>::base_tensor_type;
1645 
1646  base_tensor_type tmp;
1647  for (unsigned int i = 0; i < tmp.dimension; ++i)
1648  {
1649  // Start with the non-diagonal part
1650  value_type sum = data[dim] * sdata[dim][i];
1651  for (unsigned int d = dim + 1; d < (dim * (dim + 1) / 2); ++d)
1652  sum += data[d] * sdata[d][i];
1653  sum += sum; // sum = sum * 2.;
1654 
1655  // Now add the contributions from the diagonal
1656  for (unsigned int d = 0; d < dim; ++d)
1657  sum += data[d] * sdata[d][i];
1658  tmp[i] = sum;
1659  }
1660  return tmp;
1661  }
1662 
1663 
1664 
1665  template <int dim, typename Number, typename OtherNumber = Number>
1666  inline typename SymmetricTensorAccessors::StorageType<
1667  4,
1668  dim,
1669  typename SymmetricTensorAccessors::
1670  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type>::
1671  base_tensor_type
1672  perform_double_contraction(
1673  const typename SymmetricTensorAccessors::StorageType<4, dim, Number>::
1674  base_tensor_type &data,
1675  const typename SymmetricTensorAccessors::
1676  StorageType<4, dim, OtherNumber>::base_tensor_type &sdata)
1677  {
1678  using value_type = typename SymmetricTensorAccessors::
1679  double_contraction_result<4, 4, dim, Number, OtherNumber>::value_type;
1680  using base_tensor_type = typename SymmetricTensorAccessors::
1681  StorageType<4, dim, value_type>::base_tensor_type;
1682 
1683  const unsigned int data_dim = SymmetricTensorAccessors::
1684  StorageType<2, dim, value_type>::n_independent_components;
1685  base_tensor_type tmp;
1686  for (unsigned int i = 0; i < data_dim; ++i)
1687  for (unsigned int j = 0; j < data_dim; ++j)
1688  {
1689  // Start with the non-diagonal part
1690  for (unsigned int d = dim; d < (dim * (dim + 1) / 2); ++d)
1691  tmp[i][j] += data[i][d] * sdata[d][j];
1692  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1693 
1694  // Now add the contributions from the diagonal
1695  for (unsigned int d = 0; d < dim; ++d)
1696  tmp[i][j] += data[i][d] * sdata[d][j];
1697  }
1698  return tmp;
1699  }
1700 
1701 } // end of namespace internal
1702 
1703 
1704 
1705 template <int rank_, int dim, typename Number>
1706 template <typename OtherNumber>
1707 inline DEAL_II_ALWAYS_INLINE typename internal::SymmetricTensorAccessors::
1708  double_contraction_result<rank_, 2, dim, Number, OtherNumber>::type
1709  SymmetricTensor<rank_, dim, Number>::
1710  operator*(const SymmetricTensor<2, dim, OtherNumber> &s) const
1711 {
1712  // need to have two different function calls
1713  // because a scalar and rank-2 tensor are not
1714  // the same data type (see internal function
1715  // above)
1716  return internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1717  s.data);
1718 }
1719 
1720 
1721 
1722 template <int rank_, int dim, typename Number>
1723 template <typename OtherNumber>
1724 inline typename internal::SymmetricTensorAccessors::
1725  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type
1726  SymmetricTensor<rank_, dim, Number>::
1727  operator*(const SymmetricTensor<4, dim, OtherNumber> &s) const
1728 {
1729  typename internal::SymmetricTensorAccessors::
1730  double_contraction_result<rank_, 4, dim, Number, OtherNumber>::type tmp;
1731  tmp.data =
1732  internal::perform_double_contraction<dim, Number, OtherNumber>(data,
1733  s.data);
1734  return tmp;
1735 }
1736 
1737 
1738 
1739 // internal namespace to switch between the
1740 // access of different tensors. There used to
1741 // be explicit instantiations before for
1742 // different ranks and dimensions, but since
1743 // we now allow for templates on the data
1744 // type, and since we cannot partially
1745 // specialize the implementation, this got
1746 // into a separate namespace
1747 namespace internal
1748 {
1749  template <int dim, typename Number>
1750  inline Number &
1751  symmetric_tensor_access(const TableIndices<2> &indices,
1752  typename SymmetricTensorAccessors::
1753  StorageType<2, dim, Number>::base_tensor_type &data)
1754  {
1755  // 1d is very simple and done first
1756  if (dim == 1)
1757  return data[0];
1758 
1759  // first treat the main diagonal elements, which are stored consecutively
1760  // at the beginning
1761  if (indices[0] == indices[1])
1762  return data[indices[0]];
1763 
1764  // the rest is messier and requires a few switches.
1765  switch (dim)
1766  {
1767  case 2:
1768  // at least for the 2x2 case it is reasonably simple
1769  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1770  ((indices[0] == 0) && (indices[1] == 1)),
1771  ExcInternalError());
1772  return data[2];
1773 
1774  default:
1775  // to do the rest, sort our indices before comparing
1776  {
1777  TableIndices<2> sorted_indices(indices);
1778  sorted_indices.sort();
1779 
1780  for (unsigned int d = 0, c = 0; d < dim; ++d)
1781  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1782  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1783  return data[dim + c];
1784  Assert(false, ExcInternalError());
1785  }
1786  }
1787 
1788  static Number dummy_but_referenceable = Number();
1789  return dummy_but_referenceable;
1790  }
1791 
1792 
1793 
1794  template <int dim, typename Number>
1795  inline const Number &
1796  symmetric_tensor_access(const TableIndices<2> &indices,
1797  const typename SymmetricTensorAccessors::
1798  StorageType<2, dim, Number>::base_tensor_type &data)
1799  {
1800  // 1d is very simple and done first
1801  if (dim == 1)
1802  return data[0];
1803 
1804  // first treat the main diagonal elements, which are stored consecutively
1805  // at the beginning
1806  if (indices[0] == indices[1])
1807  return data[indices[0]];
1808 
1809  // the rest is messier and requires a few switches.
1810  switch (dim)
1811  {
1812  case 2:
1813  // at least for the 2x2 case it is reasonably simple
1814  Assert(((indices[0] == 1) && (indices[1] == 0)) ||
1815  ((indices[0] == 0) && (indices[1] == 1)),
1816  ExcInternalError());
1817  return data[2];
1818 
1819  default:
1820  // to do the rest, sort our indices before comparing
1821  {
1822  TableIndices<2> sorted_indices(indices);
1823  sorted_indices.sort();
1824 
1825  for (unsigned int d = 0, c = 0; d < dim; ++d)
1826  for (unsigned int e = d + 1; e < dim; ++e, ++c)
1827  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
1828  return data[dim + c];
1829  Assert(false, ExcInternalError());
1830  }
1831  }
1832 
1833  static Number dummy_but_referenceable = Number();
1834  return dummy_but_referenceable;
1835  }
1836 
1837 
1838 
1839  template <int dim, typename Number>
1840  inline Number &
1841  symmetric_tensor_access(const TableIndices<4> &indices,
1842  typename SymmetricTensorAccessors::
1843  StorageType<4, dim, Number>::base_tensor_type &data)
1844  {
1845  switch (dim)
1846  {
1847  case 1:
1848  return data[0][0];
1849 
1850  case 2:
1851  // each entry of the tensor can be
1852  // thought of as an entry in a
1853  // matrix that maps the rolled-out
1854  // rank-2 tensors into rolled-out
1855  // rank-2 tensors. this is the
1856  // format in which we store rank-4
1857  // tensors. determine which
1858  // position the present entry is
1859  // stored in
1860  {
1861  unsigned int base_index[2];
1862  if ((indices[0] == 0) && (indices[1] == 0))
1863  base_index[0] = 0;
1864  else if ((indices[0] == 1) && (indices[1] == 1))
1865  base_index[0] = 1;
1866  else
1867  base_index[0] = 2;
1868 
1869  if ((indices[2] == 0) && (indices[3] == 0))
1870  base_index[1] = 0;
1871  else if ((indices[2] == 1) && (indices[3] == 1))
1872  base_index[1] = 1;
1873  else
1874  base_index[1] = 2;
1875 
1876  return data[base_index[0]][base_index[1]];
1877  }
1878 
1879  case 3:
1880  // each entry of the tensor can be
1881  // thought of as an entry in a
1882  // matrix that maps the rolled-out
1883  // rank-2 tensors into rolled-out
1884  // rank-2 tensors. this is the
1885  // format in which we store rank-4
1886  // tensors. determine which
1887  // position the present entry is
1888  // stored in
1889  {
1890  unsigned int base_index[2];
1891  if ((indices[0] == 0) && (indices[1] == 0))
1892  base_index[0] = 0;
1893  else if ((indices[0] == 1) && (indices[1] == 1))
1894  base_index[0] = 1;
1895  else if ((indices[0] == 2) && (indices[1] == 2))
1896  base_index[0] = 2;
1897  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1898  ((indices[0] == 1) && (indices[1] == 0)))
1899  base_index[0] = 3;
1900  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1901  ((indices[0] == 2) && (indices[1] == 0)))
1902  base_index[0] = 4;
1903  else
1904  {
1905  Assert(((indices[0] == 1) && (indices[1] == 2)) ||
1906  ((indices[0] == 2) && (indices[1] == 1)),
1907  ExcInternalError());
1908  base_index[0] = 5;
1909  }
1910 
1911  if ((indices[2] == 0) && (indices[3] == 0))
1912  base_index[1] = 0;
1913  else if ((indices[2] == 1) && (indices[3] == 1))
1914  base_index[1] = 1;
1915  else if ((indices[2] == 2) && (indices[3] == 2))
1916  base_index[1] = 2;
1917  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1918  ((indices[2] == 1) && (indices[3] == 0)))
1919  base_index[1] = 3;
1920  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1921  ((indices[2] == 2) && (indices[3] == 0)))
1922  base_index[1] = 4;
1923  else
1924  {
1925  Assert(((indices[2] == 1) && (indices[3] == 2)) ||
1926  ((indices[2] == 2) && (indices[3] == 1)),
1927  ExcInternalError());
1928  base_index[1] = 5;
1929  }
1930 
1931  return data[base_index[0]][base_index[1]];
1932  }
1933 
1934  default:
1935  Assert(false, ExcNotImplemented());
1936  }
1937 
1938  static Number dummy;
1939  return dummy;
1940  }
1941 
1942 
1943  template <int dim, typename Number>
1944  inline const Number &
1945  symmetric_tensor_access(const TableIndices<4> &indices,
1946  const typename SymmetricTensorAccessors::
1947  StorageType<4, dim, Number>::base_tensor_type &data)
1948  {
1949  switch (dim)
1950  {
1951  case 1:
1952  return data[0][0];
1953 
1954  case 2:
1955  // each entry of the tensor can be
1956  // thought of as an entry in a
1957  // matrix that maps the rolled-out
1958  // rank-2 tensors into rolled-out
1959  // rank-2 tensors. this is the
1960  // format in which we store rank-4
1961  // tensors. determine which
1962  // position the present entry is
1963  // stored in
1964  {
1965  unsigned int base_index[2];
1966  if ((indices[0] == 0) && (indices[1] == 0))
1967  base_index[0] = 0;
1968  else if ((indices[0] == 1) && (indices[1] == 1))
1969  base_index[0] = 1;
1970  else
1971  base_index[0] = 2;
1972 
1973  if ((indices[2] == 0) && (indices[3] == 0))
1974  base_index[1] = 0;
1975  else if ((indices[2] == 1) && (indices[3] == 1))
1976  base_index[1] = 1;
1977  else
1978  base_index[1] = 2;
1979 
1980  return data[base_index[0]][base_index[1]];
1981  }
1982 
1983  case 3:
1984  // each entry of the tensor can be
1985  // thought of as an entry in a
1986  // matrix that maps the rolled-out
1987  // rank-2 tensors into rolled-out
1988  // rank-2 tensors. this is the
1989  // format in which we store rank-4
1990  // tensors. determine which
1991  // position the present entry is
1992  // stored in
1993  {
1994  unsigned int base_index[2];
1995  if ((indices[0] == 0) && (indices[1] == 0))
1996  base_index[0] = 0;
1997  else if ((indices[0] == 1) && (indices[1] == 1))
1998  base_index[0] = 1;
1999  else if ((indices[0] == 2) && (indices[1] == 2))
2000  base_index[0] = 2;
2001  else if (((indices[0] == 0) && (indices[1] == 1)) ||
2002  ((indices[0] == 1) && (indices[1] == 0)))
2003  base_index[0] = 3;
2004  else if (((indices[0] == 0) && (indices[1] == 2)) ||
2005  ((indices[0] == 2) && (indices[1] == 0)))
2006  base_index[0] = 4;
2007  else
2008  {
2009  Assert(((indices[0] == 1) && (indices[1] == 2)) ||
2010  ((indices[0] == 2) && (indices[1] == 1)),
2011  ExcInternalError());
2012  base_index[0] = 5;
2013  }
2014 
2015  if ((indices[2] == 0) && (indices[3] == 0))
2016  base_index[1] = 0;
2017  else if ((indices[2] == 1) && (indices[3] == 1))
2018  base_index[1] = 1;
2019  else if ((indices[2] == 2) && (indices[3] == 2))
2020  base_index[1] = 2;
2021  else if (((indices[2] == 0) && (indices[3] == 1)) ||
2022  ((indices[2] == 1) && (indices[3] == 0)))
2023  base_index[1] = 3;
2024  else if (((indices[2] == 0) && (indices[3] == 2)) ||
2025  ((indices[2] == 2) && (indices[3] == 0)))
2026  base_index[1] = 4;
2027  else
2028  {
2029  Assert(((indices[2] == 1) && (indices[3] == 2)) ||
2030  ((indices[2] == 2) && (indices[3] == 1)),
2031  ExcInternalError());
2032  base_index[1] = 5;
2033  }
2034 
2035  return data[base_index[0]][base_index[1]];
2036  }
2037 
2038  default:
2039  Assert(false, ExcNotImplemented());
2040  }
2041 
2042  static Number dummy;
2043  return dummy;
2044  }
2045 
2046 } // end of namespace internal
2047 
2048 
2049 
2050 template <int rank_, int dim, typename Number>
2051 inline Number &
2052 SymmetricTensor<rank_, dim, Number>::
2053 operator()(const TableIndices<rank_> &indices)
2054 {
2055  for (unsigned int r = 0; r < rank; ++r)
2056  Assert(indices[r] < dimension, ExcIndexRange(indices[r], 0, dimension));
2057  return internal::symmetric_tensor_access<dim, Number>(indices, data);
2058 }
2059 
2060 
2061 
2062 template <int rank_, int dim, typename Number>
2063 inline const Number &
2064 SymmetricTensor<rank_, dim, Number>::
2065 operator()(const TableIndices<rank_> &indices) const
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 namespace internal
2075 {
2076  namespace SymmetricTensorImplementation
2077  {
2078  template <int rank_>
2079  TableIndices<rank_>
2080  get_partially_filled_indices(const unsigned int row,
2081  const std::integral_constant<int, 2> &)
2082  {
2083  return TableIndices<rank_>(row, numbers::invalid_unsigned_int);
2084  }
2085 
2086 
2087  template <int rank_>
2088  TableIndices<rank_>
2089  get_partially_filled_indices(const unsigned int row,
2090  const std::integral_constant<int, 4> &)
2091  {
2092  return TableIndices<rank_>(row,
2093  numbers::invalid_unsigned_int,
2094  numbers::invalid_unsigned_int,
2095  numbers::invalid_unsigned_int);
2096  }
2097  } // namespace SymmetricTensorImplementation
2098 } // namespace internal
2099 
2100 
2101 template <int rank_, int dim, typename Number>
2102 internal::SymmetricTensorAccessors::
2103  Accessor<rank_, dim, true, rank_ - 1, Number>
2104  SymmetricTensor<rank_, dim, Number>::
2105  operator[](const unsigned int row) const
2106 {
2107  return internal::SymmetricTensorAccessors::
2108  Accessor<rank_, dim, true, rank_ - 1, Number>(
2109  *this,
2110  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2111  rank_>(row, std::integral_constant<int, rank_>()));
2112 }
2113 
2114 
2115 
2116 template <int rank_, int dim, typename Number>
2117 internal::SymmetricTensorAccessors::
2118  Accessor<rank_, dim, false, rank_ - 1, Number>
2119  SymmetricTensor<rank_, dim, Number>::operator[](const unsigned int row)
2120 {
2121  return internal::SymmetricTensorAccessors::
2122  Accessor<rank_, dim, false, rank_ - 1, Number>(
2123  *this,
2124  internal::SymmetricTensorImplementation::get_partially_filled_indices<
2125  rank_>(row, std::integral_constant<int, rank_>()));
2126 }
2127 
2128 
2129 
2130 template <int rank_, int dim, typename Number>
2131 inline const Number &SymmetricTensor<rank_, dim, Number>::
2132  operator[](const TableIndices<rank_> &indices) const
2133 {
2134  return operator()(indices);
2135 }
2136 
2137 
2138 
2139 template <int rank_, int dim, typename Number>
2140 inline Number &SymmetricTensor<rank_, dim, Number>::
2141  operator[](const TableIndices<rank_> &indices)
2142 {
2143  return operator()(indices);
2144 }
2145 
2146 
2147 
2148 template <int rank_, int dim, typename Number>
2149 inline Number *
2150 SymmetricTensor<rank_, dim, Number>::begin_raw()
2151 {
2152  return std::addressof(this->access_raw_entry(0));
2153 }
2154 
2155 
2156 
2157 template <int rank_, int dim, typename Number>
2158 inline const Number *
2159 SymmetricTensor<rank_, dim, Number>::begin_raw() const
2160 {
2161  return std::addressof(this->access_raw_entry(0));
2162 }
2163 
2164 
2165 
2166 template <int rank_, int dim, typename Number>
2167 inline Number *
2168 SymmetricTensor<rank_, dim, Number>::end_raw()
2169 {
2170  return begin_raw() + n_independent_components;
2171 }
2172 
2173 
2174 
2175 template <int rank_, int dim, typename Number>
2176 inline const Number *
2177 SymmetricTensor<rank_, dim, Number>::end_raw() const
2178 {
2179  return begin_raw() + n_independent_components;
2180 }
2181 
2182 
2183 
2184 namespace internal
2185 {
2186  namespace SymmetricTensorImplementation
2187  {
2188  template <int dim, typename Number>
2189  unsigned int
2190  entry_to_indices(const ::SymmetricTensor<2, dim, Number> &,
2191  const unsigned int index)
2192  {
2193  return index;
2194  }
2195 
2196 
2197  template <int dim, typename Number>
2198  ::TableIndices<2>
2199  entry_to_indices(const ::SymmetricTensor<4, dim, Number> &,
2200  const unsigned int index)
2201  {
2202  return internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
2203  base_tensor_type::unrolled_to_component_indices(index);
2204  }
2205 
2206  } // namespace SymmetricTensorImplementation
2207 } // namespace internal
2208 
2209 
2210 
2211 template <int rank_, int dim, typename Number>
2212 inline const Number &
2213 SymmetricTensor<rank_, dim, Number>::access_raw_entry(
2214  const unsigned int index) const
2215 {
2216  AssertIndexRange(index, n_independent_components);
2217  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2218  index)];
2219 }
2220 
2221 
2222 
2223 template <int rank_, int dim, typename Number>
2224 inline Number &
2225 SymmetricTensor<rank_, dim, Number>::access_raw_entry(const unsigned int index)
2226 {
2227  AssertIndexRange(index, n_independent_components);
2228  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this,
2229  index)];
2230 }
2231 
2232 
2233 
2234 namespace internal
2235 {
2236  template <int dim, typename Number>
2237  inline typename numbers::NumberTraits<Number>::real_type
2238  compute_norm(const typename SymmetricTensorAccessors::
2239  StorageType<2, dim, Number>::base_tensor_type &data)
2240  {
2241  switch (dim)
2242  {
2243  case 1:
2244  return numbers::NumberTraits<Number>::abs(data[0]);
2245 
2246  case 2:
2247  return std::sqrt(
2248  numbers::NumberTraits<Number>::abs_square(data[0]) +
2249  numbers::NumberTraits<Number>::abs_square(data[1]) +
2250  2. * numbers::NumberTraits<Number>::abs_square(data[2]));
2251 
2252  case 3:
2253  return std::sqrt(
2254  numbers::NumberTraits<Number>::abs_square(data[0]) +
2255  numbers::NumberTraits<Number>::abs_square(data[1]) +
2256  numbers::NumberTraits<Number>::abs_square(data[2]) +
2257  2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
2258  2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
2259  2. * numbers::NumberTraits<Number>::abs_square(data[5]));
2260 
2261  default:
2262  {
2263  typename numbers::NumberTraits<Number>::real_type return_value =
2264  typename numbers::NumberTraits<Number>::real_type();
2265 
2266  for (unsigned int d = 0; d < dim; ++d)
2267  return_value +=
2268  numbers::NumberTraits<Number>::abs_square(data[d]);
2269  for (unsigned int d = dim; d < (dim * dim + dim) / 2; ++d)
2270  return_value +=
2271  2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2272 
2273  return std::sqrt(return_value);
2274  }
2275  }
2276  }
2277 
2278 
2279 
2280  template <int dim, typename Number>
2281  inline typename numbers::NumberTraits<Number>::real_type
2282  compute_norm(const typename SymmetricTensorAccessors::
2283  StorageType<4, dim, Number>::base_tensor_type &data)
2284  {
2285  switch (dim)
2286  {
2287  case 1:
2288  return numbers::NumberTraits<Number>::abs(data[0][0]);
2289 
2290  default:
2291  {
2292  typename numbers::NumberTraits<Number>::real_type return_value =
2293  typename numbers::NumberTraits<Number>::real_type();
2294 
2295  const unsigned int n_independent_components = data.dimension;
2296 
2297  for (unsigned int i = 0; i < dim; ++i)
2298  for (unsigned int j = 0; j < dim; ++j)
2299  return_value +=
2300  numbers::NumberTraits<Number>::abs_square(data[i][j]);
2301  for (unsigned int i = 0; i < dim; ++i)
2302  for (unsigned int j = dim; j < n_independent_components; ++j)
2303  return_value +=
2304  2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2305  for (unsigned int i = dim; i < n_independent_components; ++i)
2306  for (unsigned int j = 0; j < dim; ++j)
2307  return_value +=
2308  2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2309  for (unsigned int i = dim; i < n_independent_components; ++i)
2310  for (unsigned int j = dim; j < n_independent_components; ++j)
2311  return_value +=
2312  4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2313 
2314  return std::sqrt(return_value);
2315  }
2316  }
2317  }
2318 
2319 } // end of namespace internal
2320 
2321 
2322 
2323 template <int rank_, int dim, typename Number>
2324 inline typename numbers::NumberTraits<Number>::real_type
2325 SymmetricTensor<rank_, dim, Number>::norm() const
2326 {
2327  return internal::compute_norm<dim, Number>(data);
2328 }
2329 
2330 
2331 
2332 namespace internal
2333 {
2334  namespace SymmetricTensorImplementation
2335  {
2336  // a function to do the unrolling from a set of indices to a
2337  // scalar index into the array in which we store the elements of
2338  // a symmetric tensor
2339  //
2340  // this function is for rank-2 tensors
2341  template <int dim>
2342  inline unsigned int
2343  component_to_unrolled_index(const TableIndices<2> &indices)
2344  {
2345  Assert(indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2346  Assert(indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2347 
2348  switch (dim)
2349  {
2350  case 1:
2351  {
2352  return 0;
2353  }
2354 
2355  case 2:
2356  {
2357  static const unsigned int table[2][2] = {{0, 2}, {2, 1}};
2358  return table[indices[0]][indices[1]];
2359  }
2360 
2361  case 3:
2362  {
2363  static const unsigned int table[3][3] = {{0, 3, 4},
2364  {3, 1, 5},
2365  {4, 5, 2}};
2366  return table[indices[0]][indices[1]];
2367  }
2368 
2369  case 4:
2370  {
2371  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2372  {4, 1, 7, 8},
2373  {5, 7, 2, 9},
2374  {6, 8, 9, 3}};
2375  return table[indices[0]][indices[1]];
2376  }
2377 
2378  default:
2379  // for the remainder, manually figure out the numbering
2380  {
2381  if (indices[0] == indices[1])
2382  return indices[0];
2383 
2384  TableIndices<2> sorted_indices(indices);
2385  sorted_indices.sort();
2386 
2387  for (unsigned int d = 0, c = 0; d < dim; ++d)
2388  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2389  if ((sorted_indices[0] == d) && (sorted_indices[1] == e))
2390  return dim + c;
2391 
2392  // should never get here:
2393  Assert(false, ExcInternalError());
2394  return 0;
2395  }
2396  }
2397  }
2398 
2399  // a function to do the unrolling from a set of indices to a
2400  // scalar index into the array in which we store the elements of
2401  // a symmetric tensor
2402  //
2403  // this function is for tensors of ranks not already handled
2404  // above
2405  template <int dim, int rank_>
2406  inline unsigned int
2407  component_to_unrolled_index(const TableIndices<rank_> &indices)
2408  {
2409  (void)indices;
2410  Assert(false, ExcNotImplemented());
2411  return numbers::invalid_unsigned_int;
2412  }
2413  } // namespace SymmetricTensorImplementation
2414 } // namespace internal
2415 
2416 
2417 template <int rank_, int dim, typename Number>
2418 inline unsigned int
2419 SymmetricTensor<rank_, dim, Number>::component_to_unrolled_index(
2420  const TableIndices<rank_> &indices)
2421 {
2422  return internal::SymmetricTensorImplementation::component_to_unrolled_index<
2423  dim>(indices);
2424 }
2425 
2426 
2427 
2428 namespace internal
2429 {
2430  namespace SymmetricTensorImplementation
2431  {
2432  // a function to do the inverse of the unrolling from a set of
2433  // indices to a scalar index into the array in which we store
2434  // the elements of a symmetric tensor. in other words, it goes
2435  // from the scalar index into the array to a set of indices of
2436  // the tensor
2437  //
2438  // this function is for rank-2 tensors
2439  template <int dim>
2440  inline TableIndices<2>
2441  unrolled_to_component_indices(const unsigned int i,
2442  const std::integral_constant<int, 2> &)
2443  {
2444  Assert(
2445  (i < ::SymmetricTensor<2, dim, double>::n_independent_components),
2446  ExcIndexRange(
2447  i,
2448  0,
2449  ::SymmetricTensor<2, dim, double>::n_independent_components));
2450  switch (dim)
2451  {
2452  case 1:
2453  {
2454  return TableIndices<2>(0, 0);
2455  }
2456 
2457  case 2:
2458  {
2459  const TableIndices<2> table[3] = {TableIndices<2>(0, 0),
2460  TableIndices<2>(1, 1),
2461  TableIndices<2>(0, 1)};
2462  return table[i];
2463  }
2464 
2465  case 3:
2466  {
2467  const TableIndices<2> table[6] = {TableIndices<2>(0, 0),
2468  TableIndices<2>(1, 1),
2469  TableIndices<2>(2, 2),
2470  TableIndices<2>(0, 1),
2471  TableIndices<2>(0, 2),
2472  TableIndices<2>(1, 2)};
2473  return table[i];
2474  }
2475 
2476  default:
2477  if (i < dim)
2478  return TableIndices<2>(i, i);
2479 
2480  for (unsigned int d = 0, c = 0; d < dim; ++d)
2481  for (unsigned int e = d + 1; e < dim; ++e, ++c)
2482  if (c == i)
2483  return TableIndices<2>(d, e);
2484 
2485  // should never get here:
2486  Assert(false, ExcInternalError());
2487  return TableIndices<2>(0, 0);
2488  }
2489  }
2490 
2491  // a function to do the inverse of the unrolling from a set of
2492  // indices to a scalar index into the array in which we store
2493  // the elements of a symmetric tensor. in other words, it goes
2494  // from the scalar index into the array to a set of indices of
2495  // the tensor
2496  //
2497  // this function is for tensors of a rank not already handled
2498  // above
2499  template <int dim, int rank_>
2500  inline TableIndices<rank_>
2501  unrolled_to_component_indices(const unsigned int i,
2502  const std::integral_constant<int, rank_> &)
2503  {
2504  (void)i;
2505  Assert(
2506  (i <
2507  ::SymmetricTensor<rank_, dim, double>::n_independent_components),
2508  ExcIndexRange(i,
2509  0,
2510  ::SymmetricTensor<rank_, dim, double>::
2511  n_independent_components));
2512  Assert(false, ExcNotImplemented());
2513  return TableIndices<rank_>();
2514  }
2515 
2516  } // namespace SymmetricTensorImplementation
2517 } // namespace internal
2518 
2519 template <int rank_, int dim, typename Number>
2520 inline TableIndices<rank_>
2521 SymmetricTensor<rank_, dim, Number>::unrolled_to_component_indices(
2522  const unsigned int i)
2523 {
2524  return internal::SymmetricTensorImplementation::unrolled_to_component_indices<
2525  dim>(i, std::integral_constant<int, rank_>());
2526 }
2527 
2528 
2529 
2530 template <int rank_, int dim, typename Number>
2531 template <class Archive>
2532 inline void
2533 SymmetricTensor<rank_, dim, Number>::serialize(Archive &ar, const unsigned int)
2534 {
2535  ar &data;
2536 }
2537 
2538 
2539 #endif // DOXYGEN
2540 
2541 /* ----------------- Non-member functions operating on tensors. ------------ */
2542 
2543 
2556 template <int rank_, int dim, typename Number, typename OtherNumber>
2557 inline SymmetricTensor<rank_,
2558  dim,
2559  typename ProductType<Number, OtherNumber>::type>
2560 operator+(const SymmetricTensor<rank_, dim, Number> & left,
2561  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2562 {
2563  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2564  tmp = left;
2565  tmp += right;
2566  return tmp;
2567 }
2568 
2569 
2582 template <int rank_, int dim, typename Number, typename OtherNumber>
2583 inline SymmetricTensor<rank_,
2584  dim,
2585  typename ProductType<Number, OtherNumber>::type>
2586 operator-(const SymmetricTensor<rank_, dim, Number> & left,
2587  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2588 {
2589  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2590  tmp = left;
2591  tmp -= right;
2592  return tmp;
2593 }
2594 
2595 
2603 template <int rank_, int dim, typename Number, typename OtherNumber>
2604 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2605 operator+(const SymmetricTensor<rank_, dim, Number> &left,
2606  const Tensor<rank_, dim, OtherNumber> & right)
2607 {
2608  return Tensor<rank_, dim, Number>(left) + right;
2609 }
2610 
2611 
2619 template <int rank_, int dim, typename Number, typename OtherNumber>
2620 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2621 operator+(const Tensor<rank_, dim, Number> & left,
2622  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2623 {
2624  return left + Tensor<rank_, dim, OtherNumber>(right);
2625 }
2626 
2627 
2635 template <int rank_, int dim, typename Number, typename OtherNumber>
2636 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2637 operator-(const SymmetricTensor<rank_, dim, Number> &left,
2638  const Tensor<rank_, dim, OtherNumber> & right)
2639 {
2640  return Tensor<rank_, dim, Number>(left) - right;
2641 }
2642 
2643 
2651 template <int rank_, int dim, typename Number, typename OtherNumber>
2652 inline Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2653 operator-(const Tensor<rank_, dim, Number> & left,
2654  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2655 {
2656  return left - Tensor<rank_, dim, OtherNumber>(right);
2657 }
2658 
2659 
2660 
2674 template <int dim, typename Number>
2675 inline Number
2676 determinant(const SymmetricTensor<2, dim, Number> &t)
2677 {
2678  switch (dim)
2679  {
2680  case 1:
2681  return t.data[0];
2682  case 2:
2683  return (t.data[0] * t.data[1] - t.data[2] * t.data[2]);
2684  case 3:
2685  {
2686  // in analogy to general tensors, but
2687  // there's something to be simplified for
2688  // the present case
2689  const Number tmp = t.data[3] * t.data[4] * t.data[5];
2690  return (tmp + tmp + t.data[0] * t.data[1] * t.data[2] -
2691  t.data[0] * t.data[5] * t.data[5] -
2692  t.data[1] * t.data[4] * t.data[4] -
2693  t.data[2] * t.data[3] * t.data[3]);
2694  }
2695  default:
2696  Assert(false, ExcNotImplemented());
2697  return internal::NumberType<Number>::value(0.0);
2698  }
2699 }
2700 
2701 
2702 
2712 template <int dim, typename Number>
2713 inline Number
2714 third_invariant(const SymmetricTensor<2, dim, Number> &t)
2715 {
2716  return determinant(t);
2717 }
2718 
2719 
2720 
2728 template <int dim, typename Number>
2729 Number
2730 trace(const SymmetricTensor<2, dim, Number> &d)
2731 {
2732  Number t = d.data[0];
2733  for (unsigned int i = 1; i < dim; ++i)
2734  t += d.data[i];
2735  return t;
2736 }
2737 
2738 
2748 template <int dim, typename Number>
2749 inline Number
2750 first_invariant(const SymmetricTensor<2, dim, Number> &t)
2751 {
2752  return trace(t);
2753 }
2754 
2755 
2768 template <typename Number>
2769 inline Number
2770 second_invariant(const SymmetricTensor<2, 1, Number> &)
2771 {
2772  return internal::NumberType<Number>::value(0.0);
2773 }
2774 
2775 
2776 
2796 template <typename Number>
2797 inline Number
2798 second_invariant(const SymmetricTensor<2, 2, Number> &t)
2799 {
2800  return t[0][0] * t[1][1] - t[0][1] * t[0][1];
2801 }
2802 
2803 
2804 
2814 template <typename Number>
2815 inline Number
2816 second_invariant(const SymmetricTensor<2, 3, Number> &t)
2817 {
2818  return (t[0][0] * t[1][1] + t[1][1] * t[2][2] + t[2][2] * t[0][0] -
2819  t[0][1] * t[0][1] - t[0][2] * t[0][2] - t[1][2] * t[1][2]);
2820 }
2821 
2822 
2823 
2832 template <typename Number>
2833 std::array<Number, 1>
2834 eigenvalues(const SymmetricTensor<2, 1, Number> &T);
2835 
2836 
2837 
2861 template <typename Number>
2862 std::array<Number, 2>
2863 eigenvalues(const SymmetricTensor<2, 2, Number> &T);
2864 
2865 
2866 
2889 template <typename Number>
2890 std::array<Number, 3>
2891 eigenvalues(const SymmetricTensor<2, 3, Number> &T);
2892 
2893 
2894 
2895 namespace internal
2896 {
2897  namespace SymmetricTensorImplementation
2898  {
2938  template <int dim, typename Number>
2939  void
2940  tridiagonalize(const ::SymmetricTensor<2, dim, Number> &A,
2941  ::Tensor<2, dim, Number> & Q,
2942  std::array<Number, dim> & d,
2943  std::array<Number, dim - 1> & e);
2944 
2945 
2946 
2988  template <int dim, typename Number>
2989  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
2990  ql_implicit_shifts(const ::SymmetricTensor<2, dim, Number> &A);
2991 
2992 
2993 
3035  template <int dim, typename Number>
3036  std::array<std::pair<Number, Tensor<1, dim, Number>>, dim>
3037  jacobi(::SymmetricTensor<2, dim, Number> A);
3038 
3039 
3040 
3056  template <typename Number>
3057  std::array<std::pair<Number, Tensor<1, 2, Number>>, 2>
3058  hybrid(const ::SymmetricTensor<2, 2, Number> &A);
3059 
3060 
3061 
3096  template <typename Number>
3097  std::array<std::pair<Number, Tensor<1, 3, Number>>, 3>
3098  hybrid(const ::SymmetricTensor<2, 3, Number> &A);
3099 
3104  template <int dim, typename Number>
3105  struct SortEigenValuesVectors
3106  {
3107  using EigValsVecs = std::pair<Number, Tensor<1, dim, Number>>;
3108  bool
3109  operator()(const EigValsVecs &lhs, const EigValsVecs &rhs)
3110  {
3111  return lhs.first > rhs.first;
3112  }
3113  };
3114 
3115  } // namespace SymmetricTensorImplementation
3116 
3117 } // namespace internal
3118 
3119 
3120 
3121 // The line below is to ensure that doxygen puts the full description
3122 // of this global enumeration into the documentation
3123 // See https://stackoverflow.com/a/1717984
3152 enum struct SymmetricTensorEigenvectorMethod
3153 {
3163  hybrid,
3173  ql_implicit_shifts,
3181  jacobi
3182 };
3183 
3184 
3185 
3215 template <int dim, typename Number>
3216 std::array<std::pair<Number, Tensor<1, dim, Number>>,
3217  std::integral_constant<int, dim>::value>
3218 eigenvectors(const SymmetricTensor<2, dim, Number> &T,
3219  const SymmetricTensorEigenvectorMethod method =
3220  SymmetricTensorEigenvectorMethod::ql_implicit_shifts);
3221 
3222 
3223 
3233 template <int rank_, int dim, typename Number>
3234 inline SymmetricTensor<rank_, dim, Number>
3235 transpose(const SymmetricTensor<rank_, dim, Number> &t)
3236 {
3237  return t;
3238 }
3239 
3240 
3241 
3251 template <int dim, typename Number>
3252 inline SymmetricTensor<2, dim, Number>
3253 deviator(const SymmetricTensor<2, dim, Number> &t)
3254 {
3255  SymmetricTensor<2, dim, Number> tmp = t;
3256 
3257  // subtract scaled trace from the diagonal
3258  const Number tr = trace(t) / dim;
3259  for (unsigned int i = 0; i < dim; ++i)
3260  tmp.data[i] -= tr;
3261 
3262  return tmp;
3263 }
3264 
3265 
3266 
3274 template <int dim, typename Number>
3275 inline SymmetricTensor<2, dim, Number>
3276 unit_symmetric_tensor()
3277 {
3278  // create a default constructed matrix filled with
3279  // zeros, then set the diagonal elements to one
3280  SymmetricTensor<2, dim, Number> tmp;
3281  switch (dim)
3282  {
3283  case 1:
3284  tmp.data[0] = 1;
3285  break;
3286  case 2:
3287  tmp.data[0] = tmp.data[1] = 1;
3288  break;
3289  case 3:
3290  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3291  break;
3292  default:
3293  for (unsigned int d = 0; d < dim; ++d)
3294  tmp.data[d] = 1;
3295  }
3296  return tmp;
3297 }
3298 
3299 
3300 
3309 template <int dim>
3310 inline SymmetricTensor<2, dim>
3311 unit_symmetric_tensor()
3312 {
3313  return unit_symmetric_tensor<dim, double>();
3314 }
3315 
3316 
3317 
3332 template <int dim, typename Number>
3333 inline SymmetricTensor<4, dim, Number>
3334 deviator_tensor()
3335 {
3336  SymmetricTensor<4, dim, Number> tmp;
3337 
3338  // fill the elements treating the diagonal
3339  for (unsigned int i = 0; i < dim; ++i)
3340  for (unsigned int j = 0; j < dim; ++j)
3341  tmp.data[i][j] = (i == j ? 1 : 0) - 1. / dim;
3342 
3343  // then fill the ones that copy over the
3344  // non-diagonal elements. note that during
3345  // the double-contraction, we handle the
3346  // off-diagonal elements twice, so simply
3347  // copying requires a weight of 1/2
3348  for (unsigned int i = dim;
3349  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3350  n_rank2_components;
3351  ++i)
3352  tmp.data[i][i] = 0.5;
3353 
3354  return tmp;
3355 }
3356 
3357 
3358 
3373 template <int dim>
3374 inline SymmetricTensor<4, dim>
3375 deviator_tensor()
3376 {
3377  return deviator_tensor<dim, double>();
3378 }
3379 
3380 
3381 
3404 template <int dim, typename Number>
3405 inline SymmetricTensor<4, dim, Number>
3406 identity_tensor()
3407 {
3408  SymmetricTensor<4, dim, Number> tmp;
3409 
3410  // fill the elements treating the diagonal
3411  for (unsigned int i = 0; i < dim; ++i)
3412  tmp.data[i][i] = 1;
3413 
3414  // then fill the ones that copy over the
3415  // non-diagonal elements. note that during
3416  // the double-contraction, we handle the
3417  // off-diagonal elements twice, so simply
3418  // copying requires a weight of 1/2
3419  for (unsigned int i = dim;
3420  i < internal::SymmetricTensorAccessors::StorageType<4, dim, Number>::
3421  n_rank2_components;
3422  ++i)
3423  tmp.data[i][i] = 0.5;
3424 
3425  return tmp;
3426 }
3427 
3428 
3429 
3451 template <int dim>
3452 inline SymmetricTensor<4, dim>
3453 identity_tensor()
3454 {
3455  return identity_tensor<dim, double>();
3456 }
3457 
3458 
3459 
3470 template <int dim, typename Number>
3471 inline SymmetricTensor<2, dim, Number>
3472 invert(const SymmetricTensor<2, dim, Number> &t)
3473 {
3474  return internal::SymmetricTensorImplementation::Inverse<2, dim, Number>::
3475  value(t);
3476 }
3477 
3478 
3479 
3491 template <int dim, typename Number>
3492 inline SymmetricTensor<4, dim, Number>
3493 invert(const SymmetricTensor<4, dim, Number> &t)
3494 {
3495  return internal::SymmetricTensorImplementation::Inverse<4, dim, Number>::
3496  value(t);
3497 }
3498 
3499 
3500 
3515 template <int dim, typename Number>
3516 inline SymmetricTensor<4, dim, Number>
3517 outer_product(const SymmetricTensor<2, dim, Number> &t1,
3518  const SymmetricTensor<2, dim, Number> &t2)
3519 {
3520  SymmetricTensor<4, dim, Number> tmp;
3521 
3522  // fill only the elements really needed
3523  for (unsigned int i = 0; i < dim; ++i)
3524  for (unsigned int j = i; j < dim; ++j)
3525  for (unsigned int k = 0; k < dim; ++k)
3526  for (unsigned int l = k; l < dim; ++l)
3527  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3528 
3529  return tmp;
3530 }
3531 
3532 
3533 
3542 template <int dim, typename Number>
3543 inline SymmetricTensor<2, dim, Number>
3544 symmetrize(const Tensor<2, dim, Number> &t)
3545 {
3546  Number array[(dim * dim + dim) / 2];
3547  for (unsigned int d = 0; d < dim; ++d)
3548  array[d] = t[d][d];
3549  for (unsigned int d = 0, c = 0; d < dim; ++d)
3550  for (unsigned int e = d + 1; e < dim; ++e, ++c)
3551  array[dim + c] = (t[d][e] + t[e][d]) * 0.5;
3552  return SymmetricTensor<2, dim, Number>(array);
3553 }
3554 
3555 
3556 
3564 template <int rank_, int dim, typename Number>
3565 inline SymmetricTensor<rank_, dim, Number>
3566 operator*(const SymmetricTensor<rank_, dim, Number> &t, const Number &factor)
3567 {
3568  SymmetricTensor<rank_, dim, Number> tt = t;
3569  tt *= factor;
3570  return tt;
3571 }
3572 
3573 
3574 
3582 template <int rank_, int dim, typename Number>
3583 inline SymmetricTensor<rank_, dim, Number>
3584 operator*(const Number &factor, const SymmetricTensor<rank_, dim, Number> &t)
3585 {
3586  // simply forward to the other operator
3587  return t * factor;
3588 }
3589 
3590 
3591 
3617 template <int rank_, int dim, typename Number, typename OtherNumber>
3618 inline SymmetricTensor<
3619  rank_,
3620  dim,
3621  typename ProductType<Number,
3622  typename EnableIfScalar<OtherNumber>::type>::type>
3623 operator*(const SymmetricTensor<rank_, dim, Number> &t,
3624  const OtherNumber & factor)
3625 {
3626  // form the product. we have to convert the two factors into the final
3627  // type via explicit casts because, for awkward reasons, the C++
3628  // standard committee saw it fit to not define an
3629  // operator*(float,std::complex<double>)
3630  // (as well as with switched arguments and double<->float).
3631  using product_type = typename ProductType<Number, OtherNumber>::type;
3632  SymmetricTensor<rank_, dim, product_type> tt(t);
3633  // we used to shorten the following by 'tt *= product_type(factor);'
3634  // which requires that a converting constructor
3635  // 'product_type::product_type(const OtherNumber) is defined.
3636  // however, a user-defined constructor is not allowed for aggregates,
3637  // e.g. VectorizedArray. therefore, we work around this issue using a
3638  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3639  // which we assume to be defined.
3640  product_type new_factor;
3641  new_factor = factor;
3642  tt *= new_factor;
3643  return tt;
3644 }
3645 
3646 
3647 
3656 template <int rank_, int dim, typename Number, typename OtherNumber>
3657 inline SymmetricTensor<
3658  rank_,
3659  dim,
3660  typename ProductType<OtherNumber,
3661  typename EnableIfScalar<Number>::type>::type>
3662 operator*(const Number & factor,
3663  const SymmetricTensor<rank_, dim, OtherNumber> &t)
3664 {
3665  // simply forward to the other operator with switched arguments
3666  return (t * factor);
3667 }
3668 
3669 
3670 
3676 template <int rank_, int dim, typename Number, typename OtherNumber>
3677 inline SymmetricTensor<
3678  rank_,
3679  dim,
3680  typename ProductType<Number,
3681  typename EnableIfScalar<OtherNumber>::type>::type>
3682 operator/(const SymmetricTensor<rank_, dim, Number> &t,
3683  const OtherNumber & factor)
3684 {
3685  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
3686  tt = t;
3687  tt /= factor;
3688  return tt;
3689 }
3690 
3691 
3692 
3699 template <int rank_, int dim>
3700 inline SymmetricTensor<rank_, dim>
3701 operator*(const SymmetricTensor<rank_, dim> &t, const double factor)
3702 {
3703  SymmetricTensor<rank_, dim> tt = t;
3704  tt *= factor;
3705  return tt;
3706 }
3707 
3708 
3709 
3716 template <int rank_, int dim>
3717 inline SymmetricTensor<rank_, dim>
3718 operator*(const double factor, const SymmetricTensor<rank_, dim> &t)
3719 {
3720  SymmetricTensor<rank_, dim> tt = t;
3721  tt *= factor;
3722  return tt;
3723 }
3724 
3725 
3726 
3732 template <int rank_, int dim>
3733 inline SymmetricTensor<rank_, dim>
3734 operator/(const SymmetricTensor<rank_, dim> &t, const double factor)
3735 {
3736  SymmetricTensor<rank_, dim> tt = t;
3737  tt /= factor;
3738  return tt;
3739 }
3740 
3750 template <int dim, typename Number, typename OtherNumber>
3751 inline typename ProductType<Number, OtherNumber>::type
3752 scalar_product(const SymmetricTensor<2, dim, Number> & t1,
3753  const SymmetricTensor<2, dim, OtherNumber> &t2)
3754 {
3755  return (t1 * t2);
3756 }
3757 
3758 
3768 template <int dim, typename Number, typename OtherNumber>
3769 inline typename ProductType<Number, OtherNumber>::type
3770 scalar_product(const SymmetricTensor<2, dim, Number> &t1,
3771  const Tensor<2, dim, OtherNumber> & t2)
3772 {
3773  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<
3774  typename ProductType<Number, OtherNumber>::type>::value(0.0);
3775  for (unsigned int i = 0; i < dim; ++i)
3776  for (unsigned int j = 0; j < dim; ++j)
3777  s += t1[i][j] * t2[i][j];
3778  return s;
3779 }
3780 
3781 
3791 template <int dim, typename Number, typename OtherNumber>
3792 inline typename ProductType<Number, OtherNumber>::type
3793 scalar_product(const Tensor<2, dim, Number> & t1,
3794  const SymmetricTensor<2, dim, OtherNumber> &t2)
3795 {
3796  return scalar_product(t2, t1);
3797 }
3798 
3799 
3815 template <typename Number, typename OtherNumber>
3816 inline void double_contract(
3817  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3818  const SymmetricTensor<4, 1, Number> & t,
3819  const SymmetricTensor<2, 1, OtherNumber> & s)
3820 {
3821  tmp[0][0] = t[0][0][0][0] * s[0][0];
3822 }
3823 
3824 
3825 
3841 template <typename Number, typename OtherNumber>
3842 inline void double_contract(
3843  SymmetricTensor<2, 1, typename ProductType<Number, OtherNumber>::type> &tmp,
3844  const SymmetricTensor<2, 1, Number> & s,
3845  const SymmetricTensor<4, 1, OtherNumber> & t)
3846 {
3847  tmp[0][0] = t[0][0][0][0] * s[0][0];
3848 }
3849 
3850 
3851 
3866 template <typename Number, typename OtherNumber>
3867 inline void double_contract(
3868  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3869  const SymmetricTensor<4, 2, Number> & t,
3870  const SymmetricTensor<2, 2, OtherNumber> & s)
3871 {
3872  const unsigned int dim = 2;
3873 
3874  for (unsigned int i = 0; i < dim; ++i)
3875  for (unsigned int j = i; j < dim; ++j)
3876  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3877  2 * t[i][j][0][1] * s[0][1];
3878 }
3879 
3880 
3881 
3897 template <typename Number, typename OtherNumber>
3898 inline void double_contract(
3899  SymmetricTensor<2, 2, typename ProductType<Number, OtherNumber>::type> &tmp,
3900  const SymmetricTensor<2, 2, Number> & s,
3901  const SymmetricTensor<4, 2, OtherNumber> & t)
3902 {
3903  const unsigned int dim = 2;
3904 
3905  for (unsigned int i = 0; i < dim; ++i)
3906  for (unsigned int j = i; j < dim; ++j)
3907  tmp[i][j] = s[0][0] * t[0][0][i][j] * +s[1][1] * t[1][1][i][j] +
3908  2 * s[0][1] * t[0][1][i][j];
3909 }
3910 
3911 
3912 
3928 template <typename Number, typename OtherNumber>
3929 inline void double_contract(
3930  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3931  const SymmetricTensor<4, 3, Number> & t,
3932  const SymmetricTensor<2, 3, OtherNumber> & s)
3933 {
3934  const unsigned int dim = 3;
3935 
3936  for (unsigned int i = 0; i < dim; ++i)
3937  for (unsigned int j = i; j < dim; ++j)
3938  tmp[i][j] = t[i][j][0][0] * s[0][0] + t[i][j][1][1] * s[1][1] +
3939  t[i][j][2][2] * s[2][2] + 2 * t[i][j][0][1] * s[0][1] +
3940  2 * t[i][j][0][2] * s[0][2] + 2 * t[i][j][1][2] * s[1][2];
3941 }
3942 
3943 
3944 
3960 template <typename Number, typename OtherNumber>
3961 inline void double_contract(
3962  SymmetricTensor<2, 3, typename ProductType<Number, OtherNumber>::type> &tmp,
3963  const SymmetricTensor<2, 3, Number> & s,
3964  const SymmetricTensor<4, 3, OtherNumber> & t)
3965 {
3966  const unsigned int dim = 3;
3967 
3968  for (unsigned int i = 0; i < dim; ++i)
3969  for (unsigned int j = i; j < dim; ++j)
3970  tmp[i][j] = s[0][0] * t[0][0][i][j] + s[1][1] * t[1][1][i][j] +
3971  s[2][2] * t[2][2][i][j] + 2 * s[0][1] * t[0][1][i][j] +
3972  2 * s[0][2] * t[0][2][i][j] + 2 * s[1][2] * t[1][2][i][j];
3973 }
3974 
3975 
3976 
3984 template <int dim, typename Number, typename OtherNumber>
3985 Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
3986 operator*(const SymmetricTensor<2, dim, Number> &src1,
3987  const Tensor<1, dim, OtherNumber> & src2)
3988 {
3989  Tensor<1, dim, typename ProductType<Number, OtherNumber>::type> dest;
3990  for (unsigned int i = 0; i < dim; ++i)
3991  for (unsigned int j = 0; j < dim; ++j)
3992  dest[i] += src1[i][j] * src2[j];
3993  return dest;
3994 }
3995 
3996 
4004 template <int dim, typename Number, typename OtherNumber>
4005 Tensor<1, dim, typename ProductType<Number, OtherNumber>::type>
4006 operator*(const Tensor<1, dim, Number> & src1,
4007  const SymmetricTensor<2, dim, OtherNumber> &src2)
4008 {
4009  // this is easy for symmetric tensors:
4010  return src2 * src1;
4011 }
4012 
4013 
4014 
4035 template <int rank_1,
4036  int rank_2,
4037  int dim,
4038  typename Number,
4039  typename OtherNumber>
4040 inline DEAL_II_ALWAYS_INLINE
4041  typename Tensor<rank_1 + rank_2 - 2,
4042  dim,
4043  typename ProductType<Number, OtherNumber>::type>::tensor_type
4044  operator*(const Tensor<rank_1, dim, Number> & src1,
4045  const SymmetricTensor<rank_2, dim, OtherNumber> &src2s)
4046 {
4047  typename Tensor<rank_1 + rank_2 - 2,
4048  dim,
4049  typename ProductType<Number, OtherNumber>::type>::tensor_type
4050  result;
4051  const Tensor<rank_2, dim, OtherNumber> src2(src2s);
4052  return src1 * src2;
4053 }
4054 
4055 
4056 
4077 template <int rank_1,
4078  int rank_2,
4079  int dim,
4080  typename Number,
4081  typename OtherNumber>
4082 inline DEAL_II_ALWAYS_INLINE
4083  typename Tensor<rank_1 + rank_2 - 2,
4084  dim,
4085  typename ProductType<Number, OtherNumber>::type>::tensor_type
4086  operator*(const SymmetricTensor<rank_1, dim, Number> &src1s,
4087  const Tensor<rank_2, dim, OtherNumber> & src2)
4088 {
4089  typename Tensor<rank_1 + rank_2 - 2,
4090  dim,
4091  typename ProductType<Number, OtherNumber>::type>::tensor_type
4092  result;
4093  const Tensor<rank_2, dim, OtherNumber> src1(src1s);
4094  return src1 * src2;
4095 }
4096 
4097 
4098 
4108 template <int dim, typename Number>
4109 inline std::ostream &
4110 operator<<(std::ostream &out, const SymmetricTensor<2, dim, Number> &t)
4111 {
4112  // make our lives a bit simpler by outputting
4113  // the tensor through the operator for the
4114  // general Tensor class
4115  Tensor<2, dim, Number> tt;
4116 
4117  for (unsigned int i = 0; i < dim; ++i)
4118  for (unsigned int j = 0; j < dim; ++j)
4119  tt[i][j] = t[i][j];
4120 
4121  return out << tt;
4122 }
4123 
4124 
4125 
4135 template <int dim, typename Number>
4136 inline std::ostream &
4137 operator<<(std::ostream &out, const SymmetricTensor<4, dim, Number> &t)
4138 {
4139  // make our lives a bit simpler by outputting
4140  // the tensor through the operator for the
4141  // general Tensor class
4142  Tensor<4, dim, Number> tt;
4143 
4144  for (unsigned int i = 0; i < dim; ++i)
4145  for (unsigned int j = 0; j < dim; ++j)
4146  for (unsigned int k = 0; k < dim; ++k)
4147  for (unsigned int l = 0; l < dim; ++l)
4148  tt[i][j][k][l] = t[i][j][k][l];
4149 
4150  return out << tt;
4151 }
4152 
4153 
4154 DEAL_II_NAMESPACE_CLOSE
4155 
4156 #endif
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:2676
SymmetricTensor::component_to_unrolled_index
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
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)
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:3682
SymmetricTensor::invert
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3472
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
internal::SymmetricTensorAccessors::double_contraction_result
Definition: symmetric_tensor.h:165
trace
Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2730
numbers::value_is_zero
bool value_is_zero(const Number &value)
Definition: numbers.h:802
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:3816
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1407
outer_product
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition: symmetric_tensor.h:3517
deviator_tensor
SymmetricTensor< 4, dim, Number > deviator_tensor()
Definition: symmetric_tensor.h:3334
third_invariant
Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2714
SymmetricTensor::operator[]
internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_-1, Number > operator[](const unsigned int row) const
SymmetricTensor::operator==
bool operator==(const SymmetricTensor &) const
SymmetricTensor::operator=
SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
SymmetricTensor::memory_consumption
static std::size_t memory_consumption()
SymmetricTensor::norm
numbers::NumberTraits< Number >::real_type norm() const
numbers::NumberTraits::abs
static real_type abs(const number &x)
Definition: numbers.h:483
SymmetricTensorEigenvectorMethod
SymmetricTensorEigenvectorMethod
Definition: symmetric_tensor.h:3152
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:2560
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:3566
SymmetricTensor::unrolled_to_component_indices
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
invert
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3472
StandardExceptions::ExcMessage
static::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor::access_raw_entry
const Number & access_raw_entry(const unsigned int unrolled_index) const
internal::SymmetricTensorAccessors::merge
TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
Definition: symmetric_tensor.h:118
SymmetricTensor< 2, spacedim >::base_tensor_type
typename base_tensor_descriptor::base_tensor_type base_tensor_type
Definition: symmetric_tensor.h:892
SymmetricTensor::begin_raw
Number * begin_raw()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1227
SymmetricTensor::data
base_tensor_type data
Definition: symmetric_tensor.h:897
SymmetricTensor::deviator
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3253
SymmetricTensor::trace
Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2730
SymmetricTensor::serialize
void serialize(Archive &ar, const unsigned int version)
SymmetricTensor::end_raw
Number * end_raw()
first_invariant
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2750
internal::SymmetricTensorImplementation::hybrid
std::array< std::pair< Number, Tensor< 1, 2, Number > >, 2 > hybrid(const ::SymmetricTensor< 2, 2, Number > &A)
symmetrize
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3544
transpose
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
Definition: symmetric_tensor.h:3235
SymmetricTensorEigenvectorMethod::ql_implicit_shifts
SymmetricTensor::operator()
Number & operator()(const TableIndices< rank_ > &indices)
swap
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1353
SymmetricTensor::operator*
internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
SymmetricTensor::determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2676
deviator
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &)
Definition: symmetric_tensor.h:3253
internal::SymmetricTensorImplementation::ql_implicit_shifts
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > ql_implicit_shifts(const ::SymmetricTensor< 2, dim, Number > &A)
Tensor
Definition: mpi.h:55
SymmetricTensor::operator*=
SymmetricTensor & operator*=(const OtherNumber &factor)
SymmetricTensor::operator+=
SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
SymmetricTensor::operator!=
bool operator!=(const SymmetricTensor &) const
StandardExceptions::ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
unit_symmetric_tensor
SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: symmetric_tensor.h:3276
scalar_product
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Definition: symmetric_tensor.h:3752
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:2586
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)
identity_tensor
SymmetricTensor< 4, dim, Number > identity_tensor()
Definition: symmetric_tensor.h:3406
second_invariant
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Definition: symmetric_tensor.h:2770
SymmetricTensor::operator-=
SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
internal::SymmetricTensorImplementation::jacobi
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
SymmetricTensor::operator-
SymmetricTensor operator-() const
StandardExceptions::ExcInternalError
static::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.11