Reference documentation for deal.II version 9.0.0
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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_symmetric_tensor_h
17 #define dealii_symmetric_tensor_h
18 
19 
20 #include <deal.II/base/tensor.h>
21 #include <deal.II/base/numbers.h>
22 #include <deal.II/base/table_indices.h>
23 #include <deal.II/base/template_constraints.h>
24 
25 #include <array>
26 #include <algorithm>
27 #include <functional>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 template <int rank, int dim, typename Number=double> class SymmetricTensor;
32 
33 template <int dim, typename Number>
34 SymmetricTensor<2,dim,Number> unit_symmetric_tensor ();
35 
36 template <int dim, typename Number>
37 SymmetricTensor<4,dim,Number> deviator_tensor ();
38 
39 template <int dim, typename Number>
40 SymmetricTensor<4,dim,Number> identity_tensor ();
41 
42 template <int dim, typename Number>
43 SymmetricTensor<2,dim,Number> invert (const SymmetricTensor<2,dim,Number> &);
44 
45 template <int dim, typename Number>
46 SymmetricTensor<4,dim,Number> invert (const SymmetricTensor<4,dim,Number> &);
47 
48 template <int dim2, typename Number> Number
49 trace (const SymmetricTensor<2,dim2,Number> &);
50 
51 template <int dim, typename Number>
52 SymmetricTensor<2,dim,Number> deviator (const SymmetricTensor<2,dim,Number> &);
53 
54 template <int dim, typename Number>
55 Number determinant (const SymmetricTensor<2,dim,Number> &);
56 
57 
58 
59 namespace internal
60 {
65  namespace SymmetricTensorImplementation
66  {
71  template <int rank, int dim, typename Number>
72  struct Inverse;
73  }
74 
79  namespace SymmetricTensorAccessors
80  {
87  inline
88  TableIndices<2> merge (const TableIndices<2> &previous_indices,
89  const unsigned int new_index,
90  const unsigned int position)
91  {
92  Assert (position < 2, ExcIndexRange (position, 0, 2));
93 
94  if (position == 0)
95  return TableIndices<2>(new_index, numbers::invalid_unsigned_int);
96  else
97  return TableIndices<2>(previous_indices[0], new_index);
98  }
99 
100 
101 
108  inline
109  TableIndices<4> merge (const TableIndices<4> &previous_indices,
110  const unsigned int new_index,
111  const unsigned int position)
112  {
113  Assert (position < 4, ExcIndexRange (position, 0, 4));
114 
115  switch (position)
116  {
117  case 0:
118  return TableIndices<4>(new_index,
119  numbers::invalid_unsigned_int,
120  numbers::invalid_unsigned_int,
121  numbers::invalid_unsigned_int);
122  case 1:
123  return TableIndices<4>(previous_indices[0],
124  new_index,
125  numbers::invalid_unsigned_int,
126  numbers::invalid_unsigned_int);
127  case 2:
128  return TableIndices<4>(previous_indices[0],
129  previous_indices[1],
130  new_index,
131  numbers::invalid_unsigned_int);
132  case 3:
133  return TableIndices<4>(previous_indices[0],
134  previous_indices[1],
135  previous_indices[2],
136  new_index);
137  }
138  Assert (false, ExcInternalError());
139  return TableIndices<4>();
140  }
141 
142 
151  template <int rank1, int rank2, int dim, typename Number, typename OtherNumber = Number>
152  struct double_contraction_result
153  {
154  typedef typename ProductType<Number,OtherNumber>::type value_type;
155  typedef ::SymmetricTensor<rank1+rank2-4,dim,value_type> type;
156  };
157 
158 
167  template <int dim, typename Number, typename OtherNumber>
168  struct double_contraction_result<2,2,dim,Number,OtherNumber>
169  {
170  typedef typename ProductType<Number,OtherNumber>::type type;
171  };
172 
173 
174 
187  template <int rank, int dim, typename Number>
188  struct StorageType;
189 
193  template <int dim, typename Number>
194  struct StorageType<2,dim,Number>
195  {
200  static const unsigned int
201  n_independent_components = (dim*dim + dim)/2;
202 
206  typedef Tensor<1,n_independent_components,Number> base_tensor_type;
207  };
208 
209 
210 
214  template <int dim, typename Number>
215  struct StorageType<4,dim,Number>
216  {
222  static const unsigned int
223  n_rank2_components = (dim*dim + dim)/2;
224 
228  static const unsigned int
229  n_independent_components = (n_rank2_components *
230  StorageType<2,dim,Number>::n_independent_components);
231 
238  typedef Tensor<2,n_rank2_components,Number> base_tensor_type;
239  };
240 
241 
242 
247  template <int rank, int dim, bool constness, typename Number>
248  struct AccessorTypes;
249 
256  template <int rank, int dim, typename Number>
257  struct AccessorTypes<rank,dim,true,Number>
258  {
259  typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
260 
261  typedef Number reference;
262  };
263 
270  template <int rank, int dim, typename Number>
271  struct AccessorTypes<rank,dim,false,Number>
272  {
273  typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
274 
275  typedef Number &reference;
276  };
277 
278 
313  template <int rank, int dim, bool constness, int P, typename Number>
314  class Accessor
315  {
316  public:
320  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
321  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
322 
323  private:
342  Accessor (tensor_type &tensor,
343  const TableIndices<rank> &previous_indices);
344 
348  Accessor (const Accessor &) = default;
349 
350  public:
351 
355  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
356 
360  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;
361 
362  private:
366  tensor_type &tensor;
367  const TableIndices<rank> previous_indices;
368 
369  // declare some other classes
370  // as friends. make sure to
371  // work around bugs in some
372  // compilers
373  template <int,int,typename> friend class ::SymmetricTensor;
374  template <int,int,bool,int,typename>
375  friend class Accessor;
376 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
377  friend class ::SymmetricTensor<rank,dim,Number>;
378  friend class Accessor<rank,dim,constness,P+1,Number>;
379 # endif
380  };
381 
382 
383 
393  template <int rank, int dim, bool constness, typename Number>
394  class Accessor<rank,dim,constness,1,Number>
395  {
396  public:
400  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
401  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
402 
403  private:
425  Accessor (tensor_type &tensor,
426  const TableIndices<rank> &previous_indices);
427 
431  Accessor () = delete;
432 
436  Accessor (const Accessor &) = default;
437 
438  public:
439 
443  reference operator [] (const unsigned int);
444 
448  reference operator [] (const unsigned int) const;
449 
450  private:
454  tensor_type &tensor;
455  const TableIndices<rank> previous_indices;
456 
457  // declare some other classes
458  // as friends. make sure to
459  // work around bugs in some
460  // compilers
461  template <int,int,typename> friend class ::SymmetricTensor;
462  template <int,int,bool,int,typename>
463  friend class SymmetricTensorAccessors::Accessor;
464 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
465  friend class ::SymmetricTensor<rank,dim,Number>;
466  friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
467 # endif
468  };
469  }
470 }
471 
472 
473 
537 template <int rank_, int dim, typename Number>
538 class SymmetricTensor
539 {
540 public:
541  static_assert(rank_%2==0, "A SymmetricTensor must have even rank!");
542 
551  static const unsigned int dimension = dim;
552 
556  static const unsigned int rank = rank_;
557 
563  static const unsigned int n_independent_components
564  = internal::SymmetricTensorAccessors::StorageType<rank_,dim,Number>::
565  n_independent_components;
566 
570  SymmetricTensor ();
571 
582  template <typename OtherNumber>
583  explicit
584  SymmetricTensor (const Tensor<2,dim,OtherNumber> &t);
585 
601  SymmetricTensor (const Number (&array) [n_independent_components]);
602 
608  template <typename OtherNumber>
609  explicit
610  SymmetricTensor (const SymmetricTensor<rank_,dim,OtherNumber> &initializer);
611 
615  Number *
616  begin_raw();
617 
621  const Number *
622  begin_raw() const;
623 
627  Number *
628  end_raw();
629 
634  const Number *
635  end_raw() const;
636 
642  template <typename OtherNumber>
643  SymmetricTensor &operator = (const SymmetricTensor<rank_,dim,OtherNumber> &rhs);
644 
651  SymmetricTensor &operator = (const Number &d);
652 
657  operator Tensor<rank_,dim,Number> () const;
658 
662  bool operator == (const SymmetricTensor &) const;
663 
667  bool operator != (const SymmetricTensor &) const;
668 
672  template <typename OtherNumber>
673  SymmetricTensor &operator += (const SymmetricTensor<rank_,dim,OtherNumber> &);
674 
678  template <typename OtherNumber>
679  SymmetricTensor &operator -= (const SymmetricTensor<rank_,dim,OtherNumber> &);
680 
685  template <typename OtherNumber>
686  SymmetricTensor &operator *= (const OtherNumber &factor);
687 
691  template <typename OtherNumber>
692  SymmetricTensor &operator /= (const OtherNumber &factor);
693 
697  SymmetricTensor operator - () const;
698 
723  template <typename OtherNumber>
724  typename internal::SymmetricTensorAccessors::double_contraction_result<rank_,2,dim,Number,OtherNumber>::type
725  operator * (const SymmetricTensor<2,dim,OtherNumber> &s) const;
726 
731  template <typename OtherNumber>
732  typename internal::SymmetricTensorAccessors::double_contraction_result<rank_,4,dim,Number,OtherNumber>::type
733  operator * (const SymmetricTensor<4,dim,OtherNumber> &s) const;
734 
738  Number &operator() (const TableIndices<rank_> &indices);
739 
743  const Number &operator() (const TableIndices<rank_> &indices) const;
744 
749  internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
750  operator [] (const unsigned int row) const;
751 
756  internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
757  operator [] (const unsigned int row);
758 
764  const Number &
765  operator [] (const TableIndices<rank_> &indices) const;
766 
772  Number &
773  operator [] (const TableIndices<rank_> &indices);
774 
780  const Number &
781  access_raw_entry (const unsigned int unrolled_index) const;
782 
788  Number &
789  access_raw_entry (const unsigned int unrolled_index);
790 
800  typename numbers::NumberTraits<Number>::real_type
801  norm () const;
802 
810  static
811  unsigned int
812  component_to_unrolled_index (const TableIndices<rank_> &indices);
813 
819  static
820  TableIndices<rank_>
821  unrolled_to_component_indices (const unsigned int i);
822 
835  void clear ();
836 
841  static std::size_t memory_consumption ();
842 
847  template <class Archive>
848  void serialize(Archive &ar, const unsigned int version);
849 
850 private:
854  typedef
855  internal::SymmetricTensorAccessors::StorageType<rank_,dim,Number>
856  base_tensor_descriptor;
857 
861  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
862 
866  base_tensor_type data;
867 
871  template <int, int, typename> friend class SymmetricTensor;
872 
876  template <int dim2, typename Number2>
877  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
878 
879  template <int dim2, typename Number2>
880  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
881 
882  template <int dim2, typename Number2>
883  friend SymmetricTensor<2,dim2,Number2>
884  deviator (const SymmetricTensor<2,dim2,Number2> &t);
885 
886  template <int dim2, typename Number2>
887  friend SymmetricTensor<2,dim2,Number2> unit_symmetric_tensor ();
888 
889  template <int dim2, typename Number2>
890  friend SymmetricTensor<4,dim2,Number2> deviator_tensor ();
891 
892  template <int dim2, typename Number2>
893  friend SymmetricTensor<4,dim2,Number2> identity_tensor ();
894 
895 
899  friend struct internal::SymmetricTensorImplementation::Inverse<2,dim,Number>;
900 
901  friend struct internal::SymmetricTensorImplementation::Inverse<4,dim,Number>;
902 };
903 
904 
905 
906 // ------------------------- inline functions ------------------------
907 
908 #ifndef DOXYGEN
909 
910 namespace internal
911 {
912  namespace SymmetricTensorAccessors
913  {
914  template <int rank_, int dim, bool constness, int P, typename Number>
915  Accessor<rank_,dim,constness,P,Number>::
916  Accessor (tensor_type &tensor,
917  const TableIndices<rank_> &previous_indices)
918  :
919  tensor (tensor),
920  previous_indices (previous_indices)
921  {}
922 
923 
924 
925  template <int rank_, int dim, bool constness, int P, typename Number>
926  Accessor<rank_,dim,constness,P-1,Number>
927  Accessor<rank_,dim,constness,P,Number>::operator[] (const unsigned int i)
928  {
929  return Accessor<rank_,dim,constness,P-1,Number> (tensor,
930  merge (previous_indices, i, rank_-P));
931  }
932 
933 
934 
935  template <int rank_, int dim, bool constness, int P, typename Number>
936  Accessor<rank_,dim,constness,P-1,Number>
937  Accessor<rank_,dim,constness,P,Number>::operator[] (const unsigned int i) const
938  {
939  return Accessor<rank_,dim,constness,P-1,Number> (tensor,
940  merge (previous_indices, i, rank_-P));
941  }
942 
943 
944 
945  template <int rank_, int dim, bool constness, typename Number>
946  Accessor<rank_,dim,constness,1,Number>::
947  Accessor (tensor_type &tensor,
948  const TableIndices<rank_> &previous_indices)
949  :
950  tensor (tensor),
951  previous_indices (previous_indices)
952  {}
953 
954 
955 
956  template <int rank_, int dim, bool constness, typename Number>
957  typename Accessor<rank_,dim,constness,1,Number>::reference
958  Accessor<rank_,dim,constness,1,Number>::operator[] (const unsigned int i)
959  {
960  return tensor(merge (previous_indices, i, rank_-1));
961  }
962 
963 
964  template <int rank_, int dim, bool constness, typename Number>
965  typename Accessor<rank_,dim,constness,1,Number>::reference
966  Accessor<rank_,dim,constness,1,Number>::operator[] (const unsigned int i) const
967  {
968  return tensor(merge (previous_indices, i, rank_-1));
969  }
970  }
971 }
972 
973 
974 
975 template <int rank_, int dim, typename Number>
976 inline
977 SymmetricTensor<rank_,dim,Number>::SymmetricTensor ()
978 {
979  // Some auto-differentiable numbers need explicit
980  // zero initialization.
981  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
982  data[i] = internal::NumberType<Number>::value(0.0);
983 }
984 
985 
986 template <int rank_, int dim, typename Number>
987 template <typename OtherNumber>
988 inline
989 SymmetricTensor<rank_,dim,Number>::SymmetricTensor (const Tensor<2,dim,OtherNumber> &t)
990 {
991  Assert (rank == 2, ExcNotImplemented());
992  switch (dim)
993  {
994  case 2:
995  Assert (t[0][1] == t[1][0], ExcInternalError());
996 
997  data[0] = t[0][0];
998  data[1] = t[1][1];
999  data[2] = t[0][1];
1000 
1001  break;
1002  case 3:
1003  Assert (t[0][1] == t[1][0], ExcInternalError());
1004  Assert (t[0][2] == t[2][0], ExcInternalError());
1005  Assert (t[1][2] == t[2][1], ExcInternalError());
1006 
1007  data[0] = t[0][0];
1008  data[1] = t[1][1];
1009  data[2] = t[2][2];
1010  data[3] = t[0][1];
1011  data[4] = t[0][2];
1012  data[5] = t[1][2];
1013 
1014  break;
1015  default:
1016  for (unsigned int d=0; d<dim; ++d)
1017  for (unsigned int e=0; e<d; ++e)
1018  Assert(t[d][e] == t[e][d], ExcInternalError());
1019 
1020  for (unsigned int d=0; d<dim; ++d)
1021  data[d] = t[d][d];
1022 
1023  for (unsigned int d=0, c=0; d<dim; ++d)
1024  for (unsigned int e=d+1; e<dim; ++e, ++c)
1025  data[dim+c] = t[d][e];
1026  }
1027 }
1028 
1029 
1030 
1031 template <int rank_, int dim, typename Number>
1032 template <typename OtherNumber>
1033 inline
1034 SymmetricTensor<rank_,dim,Number>::
1035 SymmetricTensor (const SymmetricTensor<rank_,dim,OtherNumber> &initializer)
1036 {
1037  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1038  data[i] = internal::NumberType<typename base_tensor_type::value_type>::value(initializer.data[i]);
1039 }
1040 
1041 
1042 
1043 
1044 template <int rank_, int dim, typename Number>
1045 inline
1046 SymmetricTensor<rank_,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
1047  :
1048  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1049 {
1050  // ensure that the reinterpret_cast above actually works
1051  Assert (sizeof(typename base_tensor_type::array_type)
1052  == sizeof(array),
1053  ExcInternalError());
1054 }
1055 
1056 
1057 
1058 template <int rank_, int dim, typename Number>
1059 template <typename OtherNumber>
1060 inline
1061 SymmetricTensor<rank_,dim,Number> &
1062 SymmetricTensor<rank_,dim,Number>::operator = (const SymmetricTensor<rank_,dim,OtherNumber> &t)
1063 {
1064  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1065  data[i] = t.data[i];
1066  return *this;
1067 }
1068 
1069 
1070 
1071 template <int rank_, int dim, typename Number>
1072 inline
1073 SymmetricTensor<rank_,dim,Number> &
1074 SymmetricTensor<rank_,dim,Number>::operator = (const Number &d)
1075 {
1076  Assert (numbers::value_is_zero(d), ExcMessage ("Only assignment with zero is allowed"));
1077  (void) d;
1078 
1079  data = internal::NumberType<Number>::value(0.0);
1080 
1081  return *this;
1082 }
1083 
1084 
1085 namespace internal
1086 {
1087  namespace SymmetricTensorImplementation
1088  {
1089  template <int dim, typename Number>
1090  inline DEAL_II_ALWAYS_INLINE
1091  ::Tensor<2,dim,Number>
1092  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1093  {
1094  ::Tensor<2,dim,Number> t;
1095 
1096  // diagonal entries are stored first
1097  for (unsigned int d=0; d<dim; ++d)
1098  t[d][d] = s.access_raw_entry(d);
1099 
1100  // off-diagonal entries come next, row by row
1101  for (unsigned int d=0, c=0; d<dim; ++d)
1102  for (unsigned int e=d+1; e<dim; ++e, ++c)
1103  {
1104  t[d][e] = s.access_raw_entry(dim+c);
1105  t[e][d] = s.access_raw_entry(dim+c);
1106  }
1107  return t;
1108  }
1109 
1110 
1111  template <int dim, typename Number>
1112  ::Tensor<4,dim,Number>
1113  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1114  {
1115  // utilize the symmetry properties of SymmetricTensor<4,dim>
1116  // discussed in the class documentation to avoid accessing all
1117  // independent elements of the input tensor more than once
1118  ::Tensor<4,dim,Number> t;
1119 
1120  for (unsigned int i=0; i<dim; ++i)
1121  for (unsigned int j=i; j<dim; ++j)
1122  for (unsigned int k=0; k<dim; ++k)
1123  for (unsigned int l=k; l<dim; ++l)
1124  t[TableIndices<4>(i,j,k,l)]
1125  = t[TableIndices<4>(i,j,l,k)]
1126  = t[TableIndices<4>(j,i,k,l)]
1127  = t[TableIndices<4>(j,i,l,k)]
1128  = st[TableIndices<4>(i,j,k,l)];
1129 
1130  return t;
1131  }
1132 
1133 
1134  template <typename Number>
1135  struct Inverse<2,1,Number>
1136  {
1137  static inline ::SymmetricTensor<2,1,Number>
1138  value (const ::SymmetricTensor<2,1,Number> &t)
1139  {
1140  ::SymmetricTensor<2,1,Number> tmp;
1141 
1142  tmp[0][0] = 1.0/t[0][0];
1143 
1144  return tmp;
1145  }
1146  };
1147 
1148 
1149  template <typename Number>
1150  struct Inverse<2,2,Number>
1151  {
1152  static inline ::SymmetricTensor<2,2,Number>
1153  value (const ::SymmetricTensor<2,2,Number> &t)
1154  {
1155  ::SymmetricTensor<2,2,Number> tmp;
1156 
1157  // Sympy result: ([
1158  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1159  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1160  const TableIndices<2> idx_00 (0,0);
1161  const TableIndices<2> idx_01 (0,1);
1162  const TableIndices<2> idx_11 (1,1);
1163  const Number inv_det_t
1164  = 1.0/(t[idx_00]*t[idx_11]
1165  - t[idx_01]*t[idx_01]);
1166  tmp[idx_00] = t[idx_11];
1167  tmp[idx_01] = -t[idx_01];
1168  tmp[idx_11] = t[idx_00];
1169  tmp *= inv_det_t;
1170 
1171  return tmp;
1172  }
1173  };
1174 
1175 
1176  template <typename Number>
1177  struct Inverse<2,3,Number>
1178  {
1179  static ::SymmetricTensor<2,3,Number>
1180  value (const ::SymmetricTensor<2,3,Number> &t)
1181  {
1182  ::SymmetricTensor<2,3,Number> tmp;
1183 
1184  // Sympy result: ([
1185  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1186  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1187  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
1188  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1189  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1190  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
1191  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1192  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
1193  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
1194  const TableIndices<2> idx_00 (0,0);
1195  const TableIndices<2> idx_01 (0,1);
1196  const TableIndices<2> idx_02 (0,2);
1197  const TableIndices<2> idx_11 (1,1);
1198  const TableIndices<2> idx_12 (1,2);
1199  const TableIndices<2> idx_22 (2,2);
1200  const Number inv_det_t
1201  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
1202  - t[idx_00]*t[idx_12]*t[idx_12]
1203  - t[idx_01]*t[idx_01]*t[idx_22]
1204  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
1205  - t[idx_02]*t[idx_02]*t[idx_11]);
1206  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
1207  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
1208  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
1209  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
1210  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
1211  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
1212  tmp *= inv_det_t;
1213 
1214  return tmp;
1215  }
1216  };
1217 
1218 
1219  template <typename Number>
1220  struct Inverse<4,1,Number>
1221  {
1222  static inline ::SymmetricTensor<4,1,Number>
1223  value (const ::SymmetricTensor<4,1,Number> &t)
1224  {
1225  ::SymmetricTensor<4,1,Number> tmp;
1226  tmp.data[0][0] = 1.0/t.data[0][0];
1227  return tmp;
1228  }
1229  };
1230 
1231 
1232  template <typename Number>
1233  struct Inverse<4,2,Number>
1234  {
1235  static inline ::SymmetricTensor<4,2,Number>
1236  value (const ::SymmetricTensor<4,2,Number> &t)
1237  {
1238  ::SymmetricTensor<4,2,Number> tmp;
1239 
1240  // Inverting this tensor is a little more complicated than necessary,
1241  // since we store the data of 't' as a 3x3 matrix t.data, but the
1242  // product between a rank-4 and a rank-2 tensor is really not the
1243  // product between this matrix and the 3-vector of a rhs, but rather
1244  //
1245  // B.vec = t.data * mult * A.vec
1246  //
1247  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1248  // capture the fact that we need to add up both the c_ij12*a_12 and the
1249  // c_ij21*a_21 terms.
1250  //
1251  // In addition, in this scheme, the identity tensor has the matrix
1252  // representation mult^-1.
1253  //
1254  // The inverse of 't' therefore has the matrix representation
1255  //
1256  // inv.data = mult^-1 * t.data^-1 * mult^-1
1257  //
1258  // in order to compute it, let's first compute the inverse of t.data and
1259  // put it into tmp.data; at the end of the function we then scale the
1260  // last row and column of the inverse by 1/2, corresponding to the left
1261  // and right multiplication with mult^-1.
1262  const Number t4 = t.data[0][0]*t.data[1][1],
1263  t6 = t.data[0][0]*t.data[1][2],
1264  t8 = t.data[0][1]*t.data[1][0],
1265  t00 = t.data[0][2]*t.data[1][0],
1266  t01 = t.data[0][1]*t.data[2][0],
1267  t04 = t.data[0][2]*t.data[2][0],
1268  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
1269  t8*t.data[2][2]+t00*t.data[2][1]+
1270  t01*t.data[1][2]-t04*t.data[1][1]);
1271  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
1272  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
1273  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
1274  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
1275  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
1276  tmp.data[1][2] = -(t6-t00)*t07;
1277  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
1278  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
1279  tmp.data[2][2] = (t4-t8)*t07;
1280 
1281  // scale last row and column as mentioned
1282  // above
1283  tmp.data[2][0] /= 2;
1284  tmp.data[2][1] /= 2;
1285  tmp.data[0][2] /= 2;
1286  tmp.data[1][2] /= 2;
1287  tmp.data[2][2] /= 4;
1288 
1289  return tmp;
1290  }
1291  };
1292 
1293 
1294  template <typename Number>
1295  struct Inverse<4,3,Number>
1296  {
1297  static ::SymmetricTensor<4,3,Number>
1298  value (const ::SymmetricTensor<4,3,Number> &t)
1299  {
1300  ::SymmetricTensor<4,3,Number> tmp = t;
1301 
1302  // This function follows the exact same scheme as the 2d case, except
1303  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1304  // Instead, we use the Gauss-Jordan algorithm implemented for
1305  // FullMatrix. For historical reasons the following code is copied from
1306  // there, with the tangential benefit that we do not need to copy the
1307  // tensor entries to and from the FullMatrix.
1308  const unsigned int N = 6;
1309 
1310  // First get an estimate of the size of the elements of this matrix,
1311  // for later checks whether the pivot element is large enough, or
1312  // whether we have to fear that the matrix is not regular.
1313  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1314  for (unsigned int i=0; i<N; ++i)
1315  diagonal_sum += std::fabs(tmp.data[i][i]);
1316  const Number typical_diagonal_element = diagonal_sum/static_cast<double>(N);
1317  (void)typical_diagonal_element;
1318 
1319  unsigned int p[N];
1320  for (unsigned int i=0; i<N; ++i)
1321  p[i] = i;
1322 
1323  for (unsigned int j=0; j<N; ++j)
1324  {
1325  // Pivot search: search that part of the line on and right of the
1326  // diagonal for the largest element.
1327  Number max = std::fabs(tmp.data[j][j]);
1328  unsigned int r = j;
1329  for (unsigned int i=j+1; i<N; ++i)
1330  if (std::fabs(tmp.data[i][j]) > max)
1331  {
1332  max = std::fabs(tmp.data[i][j]);
1333  r = i;
1334  }
1335 
1336  // Check whether the pivot is too small
1337  Assert(max > 1.e-16*typical_diagonal_element,
1338  ExcMessage("This tensor seems to be noninvertible"));
1339 
1340  // Row interchange
1341  if (r>j)
1342  {
1343  for (unsigned int k=0; k<N; ++k)
1344  std::swap (tmp.data[j][k], tmp.data[r][k]);
1345 
1346  std::swap (p[j], p[r]);
1347  }
1348 
1349  // Transformation
1350  const Number hr = 1./tmp.data[j][j];
1351  tmp.data[j][j] = hr;
1352  for (unsigned int k=0; k<N; ++k)
1353  {
1354  if (k==j) continue;
1355  for (unsigned int i=0; i<N; ++i)
1356  {
1357  if (i==j) continue;
1358  tmp.data[i][k] -= tmp.data[i][j]*tmp.data[j][k]*hr;
1359  }
1360  }
1361  for (unsigned int i=0; i<N; ++i)
1362  {
1363  tmp.data[i][j] *= hr;
1364  tmp.data[j][i] *= -hr;
1365  }
1366  tmp.data[j][j] = hr;
1367  }
1368 
1369  // Column interchange
1370  Number hv[N];
1371  for (unsigned int i=0; i<N; ++i)
1372  {
1373  for (unsigned int k=0; k<N; ++k)
1374  hv[p[k]] = tmp.data[i][k];
1375  for (unsigned int k=0; k<N; ++k)
1376  tmp.data[i][k] = hv[k];
1377  }
1378 
1379  // Scale rows and columns. The mult matrix
1380  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1381  for (unsigned int i=3; i<6; ++i)
1382  for (unsigned int j=0; j<3; ++j)
1383  tmp.data[i][j] /= 2;
1384 
1385  for (unsigned int i=0; i<3; ++i)
1386  for (unsigned int j=3; j<6; ++j)
1387  tmp.data[i][j] /= 2;
1388 
1389  for (unsigned int i=3; i<6; ++i)
1390  for (unsigned int j=3; j<6; ++j)
1391  tmp.data[i][j] /= 4;
1392 
1393  return tmp;
1394  }
1395  };
1396 
1397  }
1398 }
1399 
1400 
1401 
1402 template <int rank_, int dim, typename Number>
1403 inline DEAL_II_ALWAYS_INLINE
1404 SymmetricTensor<rank_,dim,Number>::
1405 operator Tensor<rank_,dim,Number> () const
1406 {
1407  return internal::SymmetricTensorImplementation::convert_to_tensor (*this);
1408 }
1409 
1410 
1411 
1412 template <int rank_, int dim, typename Number>
1413 inline
1414 bool
1415 SymmetricTensor<rank_,dim,Number>::operator ==
1416 (const SymmetricTensor<rank_,dim,Number> &t) const
1417 {
1418  return data == t.data;
1419 }
1420 
1421 
1422 
1423 template <int rank_, int dim, typename Number>
1424 inline
1425 bool
1426 SymmetricTensor<rank_,dim,Number>::operator !=
1427 (const SymmetricTensor<rank_,dim,Number> &t) const
1428 {
1429  return data != t.data;
1430 }
1431 
1432 
1433 
1434 template <int rank_, int dim, typename Number>
1435 template <typename OtherNumber>
1436 inline DEAL_II_ALWAYS_INLINE
1437 SymmetricTensor<rank_,dim,Number> &
1438 SymmetricTensor<rank_,dim,Number>::operator +=
1439 (const SymmetricTensor<rank_,dim,OtherNumber> &t)
1440 {
1441  data += t.data;
1442  return *this;
1443 }
1444 
1445 
1446 
1447 template <int rank_, int dim, typename Number>
1448 template <typename OtherNumber>
1449 inline DEAL_II_ALWAYS_INLINE
1450 SymmetricTensor<rank_,dim,Number> &
1451 SymmetricTensor<rank_,dim,Number>::operator -=
1452 (const SymmetricTensor<rank_,dim,OtherNumber> &t)
1453 {
1454  data -= t.data;
1455  return *this;
1456 }
1457 
1458 
1459 
1460 template <int rank_, int dim, typename Number>
1461 template <typename OtherNumber>
1462 inline DEAL_II_ALWAYS_INLINE
1463 SymmetricTensor<rank_,dim,Number> &
1464 SymmetricTensor<rank_,dim,Number>::operator *= (const OtherNumber &d)
1465 {
1466  data *= d;
1467  return *this;
1468 }
1469 
1470 
1471 
1472 template <int rank_, int dim, typename Number>
1473 template <typename OtherNumber>
1474 inline DEAL_II_ALWAYS_INLINE
1475 SymmetricTensor<rank_,dim,Number> &
1476 SymmetricTensor<rank_,dim,Number>::operator /= (const OtherNumber &d)
1477 {
1478  data /= d;
1479  return *this;
1480 }
1481 
1482 
1483 
1484 template <int rank_, int dim, typename Number>
1485 inline DEAL_II_ALWAYS_INLINE
1486 SymmetricTensor<rank_,dim,Number>
1487 SymmetricTensor<rank_,dim,Number>::operator - () const
1488 {
1489  SymmetricTensor tmp = *this;
1490  tmp.data = -tmp.data;
1491  return tmp;
1492 }
1493 
1494 
1495 
1496 template <int rank_, int dim, typename Number>
1497 inline DEAL_II_ALWAYS_INLINE
1498 void
1499 SymmetricTensor<rank_,dim,Number>::clear ()
1500 {
1501  data.clear ();
1502 }
1503 
1504 
1505 
1506 template <int rank_, int dim, typename Number>
1507 inline
1508 std::size_t
1509 SymmetricTensor<rank_,dim,Number>::memory_consumption ()
1510 {
1511  // all memory consists of statically allocated memory of the current
1512  // object, no pointers
1513  return sizeof(SymmetricTensor<rank_,dim,Number>);
1514 }
1515 
1516 
1517 
1518 namespace internal
1519 {
1520 
1521  template <int dim, typename Number, typename OtherNumber = Number>
1522  inline DEAL_II_ALWAYS_INLINE
1523  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type
1524  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1525  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1526  {
1527  typedef typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type result_type;
1528 
1529  switch (dim)
1530  {
1531  case 1:
1532  return data[0] * sdata[0];
1533  default:
1534  // Start with the non-diagonal part to avoid some multiplications by
1535  // 2.
1536 
1537  result_type sum = data[dim] * sdata[dim];
1538  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1539  sum += data[d] * sdata[d];
1540  sum += sum; // sum = sum * 2.;
1541 
1542  // Now add the contributions from the diagonal
1543  for (unsigned int d=0; d<dim; ++d)
1544  sum += data[d] * sdata[d];
1545  return sum;
1546  }
1547  }
1548 
1549 
1550 
1551  template <int dim, typename Number, typename OtherNumber = Number>
1552  inline
1553  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type
1554  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1555  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1556  {
1557  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type result_type;
1558  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::value_type value_type;
1559 
1560  const unsigned int data_dim =
1561  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1562  value_type tmp [data_dim];
1563  for (unsigned int i=0; i<data_dim; ++i)
1564  tmp[i] = perform_double_contraction<dim,Number,OtherNumber>(data[i], sdata);
1565  return result_type(tmp);
1566  }
1567 
1568 
1569 
1570  template <int dim, typename Number, typename OtherNumber = Number>
1571  inline
1572  typename SymmetricTensorAccessors::StorageType<2,dim,
1573  typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type
1574  >::base_tensor_type
1575  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1576  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1577  {
1578  typedef typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type value_type;
1579  typedef typename SymmetricTensorAccessors::StorageType<2,dim,value_type>::base_tensor_type base_tensor_type;
1580 
1581  base_tensor_type tmp;
1582  for (unsigned int i=0; i<tmp.dimension; ++i)
1583  {
1584  // Start with the non-diagonal part
1585  value_type sum = data[dim] * sdata[dim][i];
1586  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1587  sum += data[d] * sdata[d][i];
1588  sum += sum; // sum = sum * 2.;
1589 
1590  // Now add the contributions from the diagonal
1591  for (unsigned int d=0; d<dim; ++d)
1592  sum += data[d] * sdata[d][i];
1593  tmp[i] = sum;
1594  }
1595  return tmp;
1596  }
1597 
1598 
1599 
1600  template <int dim, typename Number, typename OtherNumber = Number>
1601  inline
1602  typename SymmetricTensorAccessors::StorageType<4,dim,
1603  typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type
1604  >::base_tensor_type
1605  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1606  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1607  {
1608  typedef typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type value_type;
1609  typedef typename SymmetricTensorAccessors::StorageType<4,dim,value_type>::base_tensor_type base_tensor_type;
1610 
1611  const unsigned int data_dim =
1612  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1613  base_tensor_type tmp;
1614  for (unsigned int i=0; i<data_dim; ++i)
1615  for (unsigned int j=0; j<data_dim; ++j)
1616  {
1617  // Start with the non-diagonal part
1618  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1619  tmp[i][j] += data[i][d] * sdata[d][j];
1620  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1621 
1622  // Now add the contributions from the diagonal
1623  for (unsigned int d=0; d<dim; ++d)
1624  tmp[i][j] += data[i][d] * sdata[d][j];
1625  }
1626  return tmp;
1627  }
1628 
1629 } // end of namespace internal
1630 
1631 
1632 
1633 template <int rank_, int dim, typename Number>
1634 template <typename OtherNumber>
1635 inline DEAL_II_ALWAYS_INLINE
1636 typename internal::SymmetricTensorAccessors::double_contraction_result<rank_,2,dim,Number,OtherNumber>::type
1637 SymmetricTensor<rank_,dim,Number>::operator * (const SymmetricTensor<2,dim,OtherNumber> &s) const
1638 {
1639  // need to have two different function calls
1640  // because a scalar and rank-2 tensor are not
1641  // the same data type (see internal function
1642  // above)
1643  return internal::perform_double_contraction<dim,Number,OtherNumber> (data, s.data);
1644 }
1645 
1646 
1647 
1648 template <int rank_, int dim, typename Number>
1649 template <typename OtherNumber>
1650 inline
1651 typename internal::SymmetricTensorAccessors::double_contraction_result<rank_,4,dim,Number,OtherNumber>::type
1652 SymmetricTensor<rank_,dim,Number>::operator * (const SymmetricTensor<4,dim,OtherNumber> &s) const
1653 {
1654  typename internal::SymmetricTensorAccessors::
1655  double_contraction_result<rank_,4,dim,Number,OtherNumber>::type tmp;
1656  tmp.data = internal::perform_double_contraction<dim,Number,OtherNumber> (data,s.data);
1657  return tmp;
1658 }
1659 
1660 
1661 
1662 // internal namespace to switch between the
1663 // access of different tensors. There used to
1664 // be explicit instantiations before for
1665 // different ranks and dimensions, but since
1666 // we now allow for templates on the data
1667 // type, and since we cannot partially
1668 // specialize the implementation, this got
1669 // into a separate namespace
1670 namespace internal
1671 {
1672  template <int dim, typename Number>
1673  inline
1674  Number &
1675  symmetric_tensor_access (const TableIndices<2> &indices,
1676  typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
1677  {
1678  // 1d is very simple and done first
1679  if (dim == 1)
1680  return data[0];
1681 
1682  // first treat the main diagonal elements, which are stored consecutively
1683  // at the beginning
1684  if (indices[0] == indices[1])
1685  return data[indices[0]];
1686 
1687  // the rest is messier and requires a few switches.
1688  switch (dim)
1689  {
1690  case 2:
1691  // at least for the 2x2 case it is reasonably simple
1692  Assert (((indices[0]==1) && (indices[1]==0)) ||
1693  ((indices[0]==0) && (indices[1]==1)),
1694  ExcInternalError());
1695  return data[2];
1696 
1697  default:
1698  // to do the rest, sort our indices before comparing
1699  {
1700  TableIndices<2> sorted_indices (indices);
1701  sorted_indices.sort ();
1702 
1703  for (unsigned int d=0, c=0; d<dim; ++d)
1704  for (unsigned int e=d+1; e<dim; ++e, ++c)
1705  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1706  return data[dim+c];
1707  Assert (false, ExcInternalError());
1708  }
1709  }
1710 
1711  static Number dummy_but_referenceable = Number();
1712  return dummy_but_referenceable;
1713  }
1714 
1715 
1716 
1717  template <int dim, typename Number>
1718  inline
1719  const Number &
1720  symmetric_tensor_access (const TableIndices<2> &indices,
1721  const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
1722  {
1723  // 1d is very simple and done first
1724  if (dim == 1)
1725  return data[0];
1726 
1727  // first treat the main diagonal elements, which are stored consecutively
1728  // at the beginning
1729  if (indices[0] == indices[1])
1730  return data[indices[0]];
1731 
1732  // the rest is messier and requires a few switches.
1733  switch (dim)
1734  {
1735  case 2:
1736  // at least for the 2x2 case it is reasonably simple
1737  Assert (((indices[0]==1) && (indices[1]==0)) ||
1738  ((indices[0]==0) && (indices[1]==1)),
1739  ExcInternalError());
1740  return data[2];
1741 
1742  default:
1743  // to do the rest, sort our indices before comparing
1744  {
1745  TableIndices<2> sorted_indices (indices);
1746  sorted_indices.sort ();
1747 
1748  for (unsigned int d=0, c=0; d<dim; ++d)
1749  for (unsigned int e=d+1; e<dim; ++e, ++c)
1750  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1751  return data[dim+c];
1752  Assert (false, ExcInternalError());
1753  }
1754  }
1755 
1756  static Number dummy_but_referenceable = Number();
1757  return dummy_but_referenceable;
1758  }
1759 
1760 
1761 
1762  template <int dim, typename Number>
1763  inline
1764  Number &
1765  symmetric_tensor_access (const TableIndices<4> &indices,
1766  typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
1767  {
1768  switch (dim)
1769  {
1770  case 1:
1771  return data[0][0];
1772 
1773  case 2:
1774  // each entry of the tensor can be
1775  // thought of as an entry in a
1776  // matrix that maps the rolled-out
1777  // rank-2 tensors into rolled-out
1778  // rank-2 tensors. this is the
1779  // format in which we store rank-4
1780  // tensors. determine which
1781  // position the present entry is
1782  // stored in
1783  {
1784  unsigned int base_index[2] ;
1785  if ((indices[0] == 0) && (indices[1] == 0))
1786  base_index[0] = 0;
1787  else if ((indices[0] == 1) && (indices[1] == 1))
1788  base_index[0] = 1;
1789  else
1790  base_index[0] = 2;
1791 
1792  if ((indices[2] == 0) && (indices[3] == 0))
1793  base_index[1] = 0;
1794  else if ((indices[2] == 1) && (indices[3] == 1))
1795  base_index[1] = 1;
1796  else
1797  base_index[1] = 2;
1798 
1799  return data[base_index[0]][base_index[1]];
1800  }
1801 
1802  case 3:
1803  // each entry of the tensor can be
1804  // thought of as an entry in a
1805  // matrix that maps the rolled-out
1806  // rank-2 tensors into rolled-out
1807  // rank-2 tensors. this is the
1808  // format in which we store rank-4
1809  // tensors. determine which
1810  // position the present entry is
1811  // stored in
1812  {
1813  unsigned int base_index[2] ;
1814  if ((indices[0] == 0) && (indices[1] == 0))
1815  base_index[0] = 0;
1816  else if ((indices[0] == 1) && (indices[1] == 1))
1817  base_index[0] = 1;
1818  else if ((indices[0] == 2) && (indices[1] == 2))
1819  base_index[0] = 2;
1820  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1821  ((indices[0] == 1) && (indices[1] == 0)))
1822  base_index[0] = 3;
1823  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1824  ((indices[0] == 2) && (indices[1] == 0)))
1825  base_index[0] = 4;
1826  else
1827  {
1828  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1829  ((indices[0] == 2) && (indices[1] == 1)),
1830  ExcInternalError());
1831  base_index[0] = 5;
1832  }
1833 
1834  if ((indices[2] == 0) && (indices[3] == 0))
1835  base_index[1] = 0;
1836  else if ((indices[2] == 1) && (indices[3] == 1))
1837  base_index[1] = 1;
1838  else if ((indices[2] == 2) && (indices[3] == 2))
1839  base_index[1] = 2;
1840  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1841  ((indices[2] == 1) && (indices[3] == 0)))
1842  base_index[1] = 3;
1843  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1844  ((indices[2] == 2) && (indices[3] == 0)))
1845  base_index[1] = 4;
1846  else
1847  {
1848  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1849  ((indices[2] == 2) && (indices[3] == 1)),
1850  ExcInternalError());
1851  base_index[1] = 5;
1852  }
1853 
1854  return data[base_index[0]][base_index[1]];
1855  }
1856 
1857  default:
1858  Assert (false, ExcNotImplemented());
1859  }
1860 
1861  static Number dummy;
1862  return dummy;
1863  }
1864 
1865 
1866  template <int dim, typename Number>
1867  inline
1868  const Number &
1869  symmetric_tensor_access (const TableIndices<4> &indices,
1870  const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
1871  {
1872  switch (dim)
1873  {
1874  case 1:
1875  return data[0][0];
1876 
1877  case 2:
1878  // each entry of the tensor can be
1879  // thought of as an entry in a
1880  // matrix that maps the rolled-out
1881  // rank-2 tensors into rolled-out
1882  // rank-2 tensors. this is the
1883  // format in which we store rank-4
1884  // tensors. determine which
1885  // position the present entry is
1886  // stored in
1887  {
1888  unsigned int base_index[2] ;
1889  if ((indices[0] == 0) && (indices[1] == 0))
1890  base_index[0] = 0;
1891  else if ((indices[0] == 1) && (indices[1] == 1))
1892  base_index[0] = 1;
1893  else
1894  base_index[0] = 2;
1895 
1896  if ((indices[2] == 0) && (indices[3] == 0))
1897  base_index[1] = 0;
1898  else if ((indices[2] == 1) && (indices[3] == 1))
1899  base_index[1] = 1;
1900  else
1901  base_index[1] = 2;
1902 
1903  return data[base_index[0]][base_index[1]];
1904  }
1905 
1906  case 3:
1907  // each entry of the tensor can be
1908  // thought of as an entry in a
1909  // matrix that maps the rolled-out
1910  // rank-2 tensors into rolled-out
1911  // rank-2 tensors. this is the
1912  // format in which we store rank-4
1913  // tensors. determine which
1914  // position the present entry is
1915  // stored in
1916  {
1917  unsigned int base_index[2] ;
1918  if ((indices[0] == 0) && (indices[1] == 0))
1919  base_index[0] = 0;
1920  else if ((indices[0] == 1) && (indices[1] == 1))
1921  base_index[0] = 1;
1922  else if ((indices[0] == 2) && (indices[1] == 2))
1923  base_index[0] = 2;
1924  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1925  ((indices[0] == 1) && (indices[1] == 0)))
1926  base_index[0] = 3;
1927  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1928  ((indices[0] == 2) && (indices[1] == 0)))
1929  base_index[0] = 4;
1930  else
1931  {
1932  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1933  ((indices[0] == 2) && (indices[1] == 1)),
1934  ExcInternalError());
1935  base_index[0] = 5;
1936  }
1937 
1938  if ((indices[2] == 0) && (indices[3] == 0))
1939  base_index[1] = 0;
1940  else if ((indices[2] == 1) && (indices[3] == 1))
1941  base_index[1] = 1;
1942  else if ((indices[2] == 2) && (indices[3] == 2))
1943  base_index[1] = 2;
1944  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1945  ((indices[2] == 1) && (indices[3] == 0)))
1946  base_index[1] = 3;
1947  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1948  ((indices[2] == 2) && (indices[3] == 0)))
1949  base_index[1] = 4;
1950  else
1951  {
1952  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1953  ((indices[2] == 2) && (indices[3] == 1)),
1954  ExcInternalError());
1955  base_index[1] = 5;
1956  }
1957 
1958  return data[base_index[0]][base_index[1]];
1959  }
1960 
1961  default:
1962  Assert (false, ExcNotImplemented());
1963  }
1964 
1965  static Number dummy;
1966  return dummy;
1967  }
1968 
1969 } // end of namespace internal
1970 
1971 
1972 
1973 template <int rank_, int dim, typename Number>
1974 inline
1975 Number &
1976 SymmetricTensor<rank_,dim,Number>::operator () (const TableIndices<rank_> &indices)
1977 {
1978  for (unsigned int r=0; r<rank; ++r)
1979  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1980  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1981 }
1982 
1983 
1984 
1985 template <int rank_, int dim, typename Number>
1986 inline
1987 const Number &
1988 SymmetricTensor<rank_,dim,Number>::operator ()
1989 (const TableIndices<rank_> &indices) const
1990 {
1991  for (unsigned int r=0; r<rank; ++r)
1992  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1993  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1994 }
1995 
1996 
1997 
1998 namespace internal
1999 {
2000  namespace SymmetricTensorImplementation
2001  {
2002  template <int rank_>
2003  TableIndices<rank_>
2004  get_partially_filled_indices (const unsigned int row,
2005  const std::integral_constant<int, 2> &)
2006  {
2007  return TableIndices<rank_> (row,
2008  numbers::invalid_unsigned_int);
2009 
2010  }
2011 
2012 
2013  template <int rank_>
2014  TableIndices<rank_>
2015  get_partially_filled_indices (const unsigned int row,
2016  const std::integral_constant<int, 4> &)
2017  {
2018  return TableIndices<rank_> (row,
2019  numbers::invalid_unsigned_int,
2020  numbers::invalid_unsigned_int,
2021  numbers::invalid_unsigned_int);
2022 
2023  }
2024  }
2025 }
2026 
2027 
2028 template <int rank_, int dim, typename Number>
2029 internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
2030 SymmetricTensor<rank_,dim,Number>::operator [] (const unsigned int row) const
2031 {
2032  return
2033  internal::SymmetricTensorAccessors::
2034  Accessor<rank_,dim,true,rank_-1,Number> (*this,
2035  internal::SymmetricTensorImplementation::get_partially_filled_indices<rank_> (row,
2036  std::integral_constant<int, rank_>()));
2037 }
2038 
2039 
2040 
2041 template <int rank_, int dim, typename Number>
2042 internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
2043 SymmetricTensor<rank_,dim,Number>::operator [] (const unsigned int row)
2044 {
2045  return
2046  internal::SymmetricTensorAccessors::
2047  Accessor<rank_,dim,false,rank_-1,Number> (*this,
2048  internal::SymmetricTensorImplementation::get_partially_filled_indices<rank_> (row,
2049  std::integral_constant<int, rank_>()));
2050 }
2051 
2052 
2053 
2054 template <int rank_, int dim, typename Number>
2055 inline
2056 const Number &
2057 SymmetricTensor<rank_,dim,Number>::operator [] (const TableIndices<rank_> &indices) const
2058 {
2059  return operator()(indices);
2060 }
2061 
2062 
2063 
2064 template <int rank_, int dim, typename Number>
2065 inline
2066 Number &
2067 SymmetricTensor<rank_,dim,Number>::operator [] (const TableIndices<rank_> &indices)
2068 {
2069  return operator()(indices);
2070 }
2071 
2072 
2073 
2074 template <int rank_, int dim, typename Number>
2075 inline
2076 Number *
2077 SymmetricTensor<rank_,dim,Number>::begin_raw()
2078 {
2079  return std::addressof(this->access_raw_entry(0));
2080 }
2081 
2082 
2083 
2084 template <int rank_, int dim, typename Number>
2085 inline
2086 const Number *
2087 SymmetricTensor<rank_,dim,Number>::begin_raw() const
2088 {
2089  return std::addressof(this->access_raw_entry(0));
2090 }
2091 
2092 
2093 
2094 template <int rank_, int dim, typename Number>
2095 inline
2096 Number *
2097 SymmetricTensor<rank_,dim,Number>::end_raw()
2098 {
2099  return begin_raw()+n_independent_components;
2100 }
2101 
2102 
2103 
2104 template <int rank_, int dim, typename Number>
2105 inline
2106 const Number *
2107 SymmetricTensor<rank_,dim,Number>::end_raw() const
2108 {
2109  return begin_raw()+n_independent_components;
2110 }
2111 
2112 
2113 
2114 namespace internal
2115 {
2116  namespace SymmetricTensorImplementation
2117  {
2118  template <int dim, typename Number>
2119  unsigned int
2120  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
2121  const unsigned int index)
2122  {
2123  return index;
2124  }
2125 
2126 
2127  template <int dim, typename Number>
2128  ::TableIndices<2>
2129  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
2130  const unsigned int index)
2131  {
2132  return
2133  internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type::
2134  unrolled_to_component_indices(index);
2135  }
2136 
2137  }
2138 }
2139 
2140 
2141 
2142 template <int rank_, int dim, typename Number>
2143 inline
2144 const Number &
2145 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index) const
2146 {
2147  AssertIndexRange (index, n_independent_components);
2148  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this, index)];
2149 }
2150 
2151 
2152 
2153 template <int rank_, int dim, typename Number>
2154 inline
2155 Number &
2156 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index)
2157 {
2158  AssertIndexRange (index, n_independent_components);
2159  return data[internal::SymmetricTensorImplementation::entry_to_indices(*this, index)];
2160 }
2161 
2162 
2163 
2164 namespace internal
2165 {
2166  template <int dim, typename Number>
2167  inline
2168  typename numbers::NumberTraits<Number>::real_type
2169  compute_norm (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
2170  {
2171  switch (dim)
2172  {
2173  case 1:
2174  return numbers::NumberTraits<Number>::abs(data[0]);
2175 
2176  case 2:
2177  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2178  numbers::NumberTraits<Number>::abs_square(data[1]) +
2179  2. * numbers::NumberTraits<Number>::abs_square(data[2]));
2180 
2181  case 3:
2182  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2183  numbers::NumberTraits<Number>::abs_square(data[1]) +
2184  numbers::NumberTraits<Number>::abs_square(data[2]) +
2185  2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
2186  2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
2187  2. * numbers::NumberTraits<Number>::abs_square(data[5]));
2188 
2189  default:
2190  {
2191  typename numbers::NumberTraits<Number>::real_type return_value
2192  = typename numbers::NumberTraits<Number>::real_type();
2193 
2194  for (unsigned int d=0; d<dim; ++d)
2195  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
2196  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
2197  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2198 
2199  return std::sqrt(return_value);
2200  }
2201  }
2202  }
2203 
2204 
2205 
2206  template <int dim, typename Number>
2207  inline
2208  typename numbers::NumberTraits<Number>::real_type
2209  compute_norm (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
2210  {
2211  switch (dim)
2212  {
2213  case 1:
2214  return numbers::NumberTraits<Number>::abs (data[0][0]);
2215 
2216  default:
2217  {
2218  typename numbers::NumberTraits<Number>::real_type return_value
2219  = typename numbers::NumberTraits<Number>::real_type();
2220 
2221  const unsigned int n_independent_components = data.dimension;
2222 
2223  for (unsigned int i=0; i<dim; ++i)
2224  for (unsigned int j=0; j<dim; ++j)
2225  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
2226  for (unsigned int i=0; i<dim; ++i)
2227  for (unsigned int j=dim; j<n_independent_components; ++j)
2228  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2229  for (unsigned int i=dim; i<n_independent_components; ++i)
2230  for (unsigned int j=0; j<dim; ++j)
2231  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2232  for (unsigned int i=dim; i<n_independent_components; ++i)
2233  for (unsigned int j=dim; j<n_independent_components; ++j)
2234  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2235 
2236  return std::sqrt(return_value);
2237  }
2238  }
2239  }
2240 
2241 } // end of namespace internal
2242 
2243 
2244 
2245 template <int rank_, int dim, typename Number>
2246 inline
2247 typename numbers::NumberTraits<Number>::real_type
2248 SymmetricTensor<rank_,dim,Number>::norm () const
2249 {
2250  return internal::compute_norm<dim,Number> (data);
2251 }
2252 
2253 
2254 
2255 namespace internal
2256 {
2257  namespace SymmetricTensorImplementation
2258  {
2259  namespace
2260  {
2261  // a function to do the unrolling from a set of indices to a
2262  // scalar index into the array in which we store the elements of
2263  // a symmetric tensor
2264  //
2265  // this function is for rank-2 tensors
2266  template <int dim>
2267  inline
2268  unsigned int
2269  component_to_unrolled_index
2270  (const TableIndices<2> &indices)
2271  {
2272  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2273  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2274 
2275  switch (dim)
2276  {
2277  case 1:
2278  {
2279  return 0;
2280  }
2281 
2282  case 2:
2283  {
2284  static const unsigned int table[2][2] = {{0, 2},
2285  {2, 1}
2286  };
2287  return table[indices[0]][indices[1]];
2288  }
2289 
2290  case 3:
2291  {
2292  static const unsigned int table[3][3] = {{0, 3, 4},
2293  {3, 1, 5},
2294  {4, 5, 2}
2295  };
2296  return table[indices[0]][indices[1]];
2297  }
2298 
2299  case 4:
2300  {
2301  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2302  {4, 1, 7, 8},
2303  {5, 7, 2, 9},
2304  {6, 8, 9, 3}
2305  };
2306  return table[indices[0]][indices[1]];
2307  }
2308 
2309  default:
2310  // for the remainder, manually figure out the numbering
2311  {
2312  if (indices[0] == indices[1])
2313  return indices[0];
2314 
2315  TableIndices<2> sorted_indices (indices);
2316  sorted_indices.sort ();
2317 
2318  for (unsigned int d=0, c=0; d<dim; ++d)
2319  for (unsigned int e=d+1; e<dim; ++e, ++c)
2320  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
2321  return dim+c;
2322 
2323  // should never get here:
2324  Assert(false, ExcInternalError());
2325  return 0;
2326  }
2327  }
2328  }
2329 
2330  // a function to do the unrolling from a set of indices to a
2331  // scalar index into the array in which we store the elements of
2332  // a symmetric tensor
2333  //
2334  // this function is for tensors of ranks not already handled
2335  // above
2336  template <int dim, int rank_>
2337  inline
2338  unsigned int
2339  component_to_unrolled_index
2340  (const TableIndices<rank_> &indices)
2341  {
2342  (void)indices;
2343  Assert (false, ExcNotImplemented());
2344  return numbers::invalid_unsigned_int;
2345  }
2346  }
2347  }
2348 }
2349 
2350 
2351 template <int rank_, int dim, typename Number>
2352 inline
2353 unsigned int
2354 SymmetricTensor<rank_,dim,Number>::component_to_unrolled_index
2355 (const TableIndices<rank_> &indices)
2356 {
2357  return internal::SymmetricTensorImplementation::component_to_unrolled_index<dim> (indices);
2358 }
2359 
2360 
2361 
2362 namespace internal
2363 {
2364  namespace SymmetricTensorImplementation
2365  {
2366  namespace
2367  {
2368  // a function to do the inverse of the unrolling from a set of
2369  // indices to a scalar index into the array in which we store
2370  // the elements of a symmetric tensor. in other words, it goes
2371  // from the scalar index into the array to a set of indices of
2372  // the tensor
2373  //
2374  // this function is for rank-2 tensors
2375  template <int dim>
2376  inline
2377  TableIndices<2>
2378  unrolled_to_component_indices
2379  (const unsigned int i,
2380  const std::integral_constant<int, 2> &)
2381  {
2382  Assert ((i < ::SymmetricTensor<2,dim,double>::n_independent_components),
2383  ExcIndexRange(i, 0, ::SymmetricTensor<2,dim,double>::n_independent_components));
2384  switch (dim)
2385  {
2386  case 1:
2387  {
2388  return TableIndices<2>(0,0);
2389  }
2390 
2391  case 2:
2392  {
2393  const TableIndices<2> table[3] =
2394  {
2395  TableIndices<2> (0,0),
2396  TableIndices<2> (1,1),
2397  TableIndices<2> (0,1)
2398  };
2399  return table[i];
2400  }
2401 
2402  case 3:
2403  {
2404  const TableIndices<2> table[6] =
2405  {
2406  TableIndices<2> (0,0),
2407  TableIndices<2> (1,1),
2408  TableIndices<2> (2,2),
2409  TableIndices<2> (0,1),
2410  TableIndices<2> (0,2),
2411  TableIndices<2> (1,2)
2412  };
2413  return table[i];
2414  }
2415 
2416  default:
2417  if (i<dim)
2418  return TableIndices<2> (i,i);
2419 
2420  for (unsigned int d=0, c=0; d<dim; ++d)
2421  for (unsigned int e=d+1; e<dim; ++e, ++c)
2422  if (c==i)
2423  return TableIndices<2>(d,e);
2424 
2425  // should never get here:
2426  Assert(false, ExcInternalError());
2427  return TableIndices<2>(0, 0);
2428  }
2429  }
2430 
2431  // a function to do the inverse of the unrolling from a set of
2432  // indices to a scalar index into the array in which we store
2433  // the elements of a symmetric tensor. in other words, it goes
2434  // from the scalar index into the array to a set of indices of
2435  // the tensor
2436  //
2437  // this function is for tensors of a rank not already handled
2438  // above
2439  template <int dim, int rank_>
2440  inline
2441  TableIndices<rank_>
2442  unrolled_to_component_indices
2443  (const unsigned int i,
2444  const std::integral_constant<int, rank_> &)
2445  {
2446  (void)i;
2447  Assert ((i < ::SymmetricTensor<rank_,dim,double>::n_independent_components),
2448  ExcIndexRange(i, 0, ::SymmetricTensor<rank_,dim,double>::n_independent_components));
2449  Assert (false, ExcNotImplemented());
2450  return TableIndices<rank_>();
2451  }
2452 
2453  }
2454  }
2455 }
2456 
2457 template <int rank_, int dim, typename Number>
2458 inline
2459 TableIndices<rank_>
2460 SymmetricTensor<rank_,dim,Number>::unrolled_to_component_indices
2461 (const unsigned int i)
2462 {
2463  return
2464  internal::SymmetricTensorImplementation::unrolled_to_component_indices<dim> (i,
2465  std::integral_constant<int, rank_>());
2466 }
2467 
2468 
2469 
2470 template <int rank_, int dim, typename Number>
2471 template <class Archive>
2472 inline
2473 void
2474 SymmetricTensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
2475 {
2476  ar &data;
2477 }
2478 
2479 
2480 #endif // DOXYGEN
2481 
2482 /* ----------------- Non-member functions operating on tensors. ------------ */
2483 
2484 
2497 template <int rank_, int dim, typename Number, typename OtherNumber>
2498 inline
2499 SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2500 operator+(const SymmetricTensor<rank_, dim, Number> &left,
2501  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2502 {
2503  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type> tmp = left;
2504  tmp += right;
2505  return tmp;
2506 }
2507 
2508 
2521 template <int rank_, int dim, typename Number, typename OtherNumber>
2522 inline
2523 SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2524 operator-(const SymmetricTensor<rank_, dim, Number> &left,
2525  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2526 {
2527  SymmetricTensor<rank_, dim, typename ProductType<Number, OtherNumber>::type> tmp = left;
2528  tmp -= right;
2529  return tmp;
2530 }
2531 
2532 
2540 template <int rank_, int dim, typename Number, typename OtherNumber>
2541 inline
2542 Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2543 operator+(const SymmetricTensor<rank_, dim, Number> &left,
2544  const Tensor<rank_, dim, OtherNumber> &right)
2545 {
2546  return Tensor<rank_, dim, Number>(left) + right;
2547 }
2548 
2549 
2557 template <int rank_, int dim, typename Number, typename OtherNumber>
2558 inline
2559 Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2560 operator+(const Tensor<rank_, dim, Number> &left,
2561  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2562 {
2563  return left + Tensor<rank_, dim, OtherNumber>(right);
2564 }
2565 
2566 
2574 template <int rank_, int dim, typename Number, typename OtherNumber>
2575 inline
2576 Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2577 operator-(const SymmetricTensor<rank_, dim, Number> &left,
2578  const Tensor<rank_, dim, OtherNumber> &right)
2579 {
2580  return Tensor<rank_, dim, Number>(left) - right;
2581 }
2582 
2583 
2591 template <int rank_, int dim, typename Number, typename OtherNumber>
2592 inline
2593 Tensor<rank_, dim, typename ProductType<Number, OtherNumber>::type>
2594 operator-(const Tensor<rank_, dim, Number> &left,
2595  const SymmetricTensor<rank_, dim, OtherNumber> &right)
2596 {
2597  return left - Tensor<rank_, dim, OtherNumber>(right);
2598 }
2599 
2600 
2601 
2615 template <int dim, typename Number>
2616 inline
2617 Number determinant (const SymmetricTensor<2,dim,Number> &t)
2618 {
2619  switch (dim)
2620  {
2621  case 1:
2622  return t.data[0];
2623  case 2:
2624  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2625  case 3:
2626  {
2627  // in analogy to general tensors, but
2628  // there's something to be simplified for
2629  // the present case
2630  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2631  return ( tmp + tmp
2632  +t.data[0]*t.data[1]*t.data[2]
2633  -t.data[0]*t.data[5]*t.data[5]
2634  -t.data[1]*t.data[4]*t.data[4]
2635  -t.data[2]*t.data[3]*t.data[3]);
2636  }
2637  default:
2638  Assert (false, ExcNotImplemented());
2639  return internal::NumberType<Number>::value(0.0);
2640  }
2641 }
2642 
2643 
2644 
2654 template <int dim, typename Number>
2655 inline
2656 Number third_invariant (const SymmetricTensor<2,dim,Number> &t)
2657 {
2658  return determinant (t);
2659 }
2660 
2661 
2662 
2670 template <int dim, typename Number>
2671 Number trace (const SymmetricTensor<2,dim,Number> &d)
2672 {
2673  Number t = d.data[0];
2674  for (unsigned int i=1; i<dim; ++i)
2675  t += d.data[i];
2676  return t;
2677 }
2678 
2679 
2689 template <int dim, typename Number>
2690 inline
2691 Number first_invariant (const SymmetricTensor<2,dim,Number> &t)
2692 {
2693  return trace (t);
2694 }
2695 
2696 
2709 template <typename Number>
2710 inline
2711 Number second_invariant (const SymmetricTensor<2,1,Number> &)
2712 {
2713  return internal::NumberType<Number>::value(0.0);
2714 }
2715 
2716 
2717 
2738 template <typename Number>
2739 inline
2740 Number second_invariant (const SymmetricTensor<2,2,Number> &t)
2741 {
2742  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2743 }
2744 
2745 
2746 
2756 template <typename Number>
2757 inline
2758 Number second_invariant (const SymmetricTensor<2,3,Number> &t)
2759 {
2760  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2761  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2762 }
2763 
2764 
2765 
2774 template <typename Number>
2775 std::array<Number,1>
2776 eigenvalues (const SymmetricTensor<2,1,Number> &T);
2777 
2778 
2779 
2801 template <typename Number>
2802 std::array<Number,2>
2803 eigenvalues (const SymmetricTensor<2,2,Number> &T);
2804 
2805 
2806 
2826 template <typename Number>
2827 std::array<Number,3>
2828 eigenvalues (const SymmetricTensor<2,3,Number> &T);
2829 
2830 
2831 
2832 namespace internal
2833 {
2834  namespace SymmetricTensorImplementation
2835  {
2872  template <int dim, typename Number>
2873  void
2874  tridiagonalize (const ::SymmetricTensor<2,dim,Number> &A,
2875  ::Tensor<2,dim,Number> &Q,
2876  std::array<Number,dim> &d,
2877  std::array<Number,dim-1> &e);
2878 
2879 
2880 
2919  template <int dim, typename Number>
2920  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2921  ql_implicit_shifts (const ::SymmetricTensor<2,dim,Number> &A);
2922 
2923 
2924 
2963  template <int dim, typename Number>
2964  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2965  jacobi (::SymmetricTensor<2,dim,Number> A);
2966 
2967 
2968 
2984  template <typename Number>
2985  std::array<std::pair<Number, Tensor<1,2,Number> >,2>
2986  hybrid (const ::SymmetricTensor<2,2,Number> &A);
2987 
2988 
2989 
3022  template <typename Number>
3023  std::array<std::pair<Number, Tensor<1,3,Number> >,3>
3024  hybrid (const ::SymmetricTensor<2,3,Number> &A);
3025 
3026  namespace
3027  {
3028 
3033  template <int dim, typename Number>
3034  struct SortEigenValuesVectors
3035  {
3036  typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
3037  bool operator() (const EigValsVecs &lhs,
3038  const EigValsVecs &rhs)
3039  {
3040  return lhs.first > rhs.first;
3041  }
3042  };
3043 
3044  }
3045 
3046  } // namespace SymmetricTensor
3047 
3048 } // namespace internal
3049 
3050 
3051 
3052 // The line below is to ensure that doxygen puts the full description
3053 // of this global enumeration into the documentation
3054 // See https://stackoverflow.com/a/1717984
3081 enum struct SymmetricTensorEigenvectorMethod
3082 {
3092  hybrid,
3102  ql_implicit_shifts,
3110  jacobi
3111 };
3112 
3113 
3114 
3142 template <int dim, typename Number>
3143 std::array<std::pair<Number, Tensor<1,dim,Number> >,std::integral_constant<int, dim>::value>
3144 eigenvectors (const SymmetricTensor<2,dim,Number> &T,
3145  const SymmetricTensorEigenvectorMethod method = SymmetricTensorEigenvectorMethod::ql_implicit_shifts);
3146 
3147 
3148 
3158 template <int rank_, int dim, typename Number>
3159 inline
3160 SymmetricTensor<rank_,dim,Number>
3161 transpose (const SymmetricTensor<rank_,dim,Number> &t)
3162 {
3163  return t;
3164 }
3165 
3166 
3167 
3177 template <int dim, typename Number>
3178 inline
3179 SymmetricTensor<2,dim,Number>
3180 deviator (const SymmetricTensor<2,dim,Number> &t)
3181 {
3182  SymmetricTensor<2,dim,Number> tmp = t;
3183 
3184  // subtract scaled trace from the diagonal
3185  const Number tr = trace(t) / dim;
3186  for (unsigned int i=0; i<dim; ++i)
3187  tmp.data[i] -= tr;
3188 
3189  return tmp;
3190 }
3191 
3192 
3193 
3201 template <int dim, typename Number>
3202 inline
3203 SymmetricTensor<2,dim,Number>
3204 unit_symmetric_tensor ()
3205 {
3206  // create a default constructed matrix filled with
3207  // zeros, then set the diagonal elements to one
3208  SymmetricTensor<2,dim,Number> tmp;
3209  switch (dim)
3210  {
3211  case 1:
3212  tmp.data[0] = 1;
3213  break;
3214  case 2:
3215  tmp.data[0] = tmp.data[1] = 1;
3216  break;
3217  case 3:
3218  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3219  break;
3220  default:
3221  for (unsigned int d=0; d<dim; ++d)
3222  tmp.data[d] = 1;
3223  }
3224  return tmp;
3225 }
3226 
3227 
3228 
3237 template <int dim>
3238 inline
3239 SymmetricTensor<2,dim>
3240 unit_symmetric_tensor ()
3241 {
3242  return unit_symmetric_tensor<dim,double>();
3243 }
3244 
3245 
3246 
3261 template <int dim, typename Number>
3262 inline
3263 SymmetricTensor<4,dim,Number>
3264 deviator_tensor ()
3265 {
3266  SymmetricTensor<4,dim,Number> tmp;
3267 
3268  // fill the elements treating the diagonal
3269  for (unsigned int i=0; i<dim; ++i)
3270  for (unsigned int j=0; j<dim; ++j)
3271  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
3272 
3273  // then fill the ones that copy over the
3274  // non-diagonal elements. note that during
3275  // the double-contraction, we handle the
3276  // off-diagonal elements twice, so simply
3277  // copying requires a weight of 1/2
3278  for (unsigned int i=dim;
3279  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3280  ++i)
3281  tmp.data[i][i] = 0.5;
3282 
3283  return tmp;
3284 }
3285 
3286 
3287 
3302 template <int dim>
3303 inline
3304 SymmetricTensor<4,dim>
3305 deviator_tensor ()
3306 {
3307  return deviator_tensor<dim,double>();
3308 }
3309 
3310 
3311 
3334 template <int dim, typename Number>
3335 inline
3336 SymmetricTensor<4,dim,Number>
3337 identity_tensor ()
3338 {
3339  SymmetricTensor<4,dim,Number> tmp;
3340 
3341  // fill the elements treating the diagonal
3342  for (unsigned int i=0; i<dim; ++i)
3343  tmp.data[i][i] = 1;
3344 
3345  // then fill the ones that copy over the
3346  // non-diagonal elements. note that during
3347  // the double-contraction, we handle the
3348  // off-diagonal elements twice, so simply
3349  // copying requires a weight of 1/2
3350  for (unsigned int i=dim;
3351  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3352  ++i)
3353  tmp.data[i][i] = 0.5;
3354 
3355  return tmp;
3356 }
3357 
3358 
3359 
3381 template <int dim>
3382 inline
3383 SymmetricTensor<4,dim>
3384 identity_tensor ()
3385 {
3386  return identity_tensor<dim,double>();
3387 }
3388 
3389 
3390 
3401 template <int dim, typename Number>
3402 inline
3403 SymmetricTensor<2,dim,Number>
3404 invert (const SymmetricTensor<2,dim,Number> &t)
3405 {
3406  return internal::SymmetricTensorImplementation::Inverse<2,dim,Number>::value(t);
3407 }
3408 
3409 
3410 
3422 template <int dim, typename Number>
3423 inline
3424 SymmetricTensor<4,dim,Number>
3425 invert (const SymmetricTensor<4,dim,Number> &t)
3426 {
3427  return internal::SymmetricTensorImplementation::Inverse<4,dim,Number>::value(t);
3428 }
3429 
3430 
3431 
3446 template <int dim, typename Number>
3447 inline
3448 SymmetricTensor<4,dim,Number>
3449 outer_product (const SymmetricTensor<2,dim,Number> &t1,
3450  const SymmetricTensor<2,dim,Number> &t2)
3451 {
3452  SymmetricTensor<4,dim,Number> tmp;
3453 
3454  // fill only the elements really needed
3455  for (unsigned int i=0; i<dim; ++i)
3456  for (unsigned int j=i; j<dim; ++j)
3457  for (unsigned int k=0; k<dim; ++k)
3458  for (unsigned int l=k; l<dim; ++l)
3459  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3460 
3461  return tmp;
3462 }
3463 
3464 
3465 
3474 template <int dim,typename Number>
3475 inline
3476 SymmetricTensor<2,dim,Number>
3477 symmetrize (const Tensor<2,dim,Number> &t)
3478 {
3479  Number array[(dim*dim+dim)/2];
3480  for (unsigned int d=0; d<dim; ++d)
3481  array[d] = t[d][d];
3482  for (unsigned int d=0, c=0; d<dim; ++d)
3483  for (unsigned int e=d+1; e<dim; ++e, ++c)
3484  array[dim+c] = (t[d][e]+t[e][d])*0.5;
3485  return SymmetricTensor<2,dim,Number>(array);
3486 }
3487 
3488 
3489 
3497 template <int rank_, int dim, typename Number>
3498 inline
3499 SymmetricTensor<rank_,dim,Number>
3500 operator * (const SymmetricTensor<rank_,dim,Number> &t,
3501  const Number &factor)
3502 {
3503  SymmetricTensor<rank_,dim,Number> tt = t;
3504  tt *= factor;
3505  return tt;
3506 }
3507 
3508 
3509 
3517 template <int rank_, int dim, typename Number>
3518 inline
3519 SymmetricTensor<rank_,dim,Number>
3520 operator * (const Number &factor,
3521  const SymmetricTensor<rank_,dim,Number> &t)
3522 {
3523  // simply forward to the other operator
3524  return t*factor;
3525 }
3526 
3527 
3528 
3554 template <int rank_, int dim, typename Number, typename OtherNumber>
3555 inline
3556 SymmetricTensor<rank_,dim,typename ProductType<Number,typename EnableIfScalar<OtherNumber>::type>::type>
3557 operator * (const SymmetricTensor<rank_,dim,Number> &t,
3558  const OtherNumber &factor)
3559 {
3560  // form the product. we have to convert the two factors into the final
3561  // type via explicit casts because, for awkward reasons, the C++
3562  // standard committee saw it fit to not define an
3563  // operator*(float,std::complex<double>)
3564  // (as well as with switched arguments and double<->float).
3565  typedef typename ProductType<Number,OtherNumber>::type product_type;
3566  SymmetricTensor<rank_,dim,product_type> tt(t);
3567  // we used to shorten the following by 'tt *= product_type(factor);'
3568  // which requires that a converting constructor
3569  // 'product_type::product_type(const OtherNumber) is defined.
3570  // however, a user-defined constructor is not allowed for aggregates,
3571  // e.g. VectorizedArray. therefore, we work around this issue using a
3572  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3573  // which we assume to be defined.
3574  product_type new_factor;
3575  new_factor = factor;
3576  tt *= new_factor;
3577  return tt;
3578 }
3579 
3580 
3581 
3590 template <int rank_, int dim, typename Number, typename OtherNumber>
3591 inline
3592 SymmetricTensor<rank_,dim,typename ProductType<OtherNumber,typename EnableIfScalar<Number>::type>::type>
3593 operator * (const Number &factor,
3594  const SymmetricTensor<rank_,dim,OtherNumber> &t)
3595 {
3596  // simply forward to the other operator with switched arguments
3597  return (t*factor);
3598 }
3599 
3600 
3601 
3607 template <int rank_, int dim, typename Number, typename OtherNumber>
3608 inline
3609 SymmetricTensor<rank_,dim,typename ProductType<Number,typename EnableIfScalar<OtherNumber>::type>::type>
3610 operator / (const SymmetricTensor<rank_,dim,Number> &t,
3611  const OtherNumber &factor)
3612 {
3613  SymmetricTensor<rank_,dim,typename ProductType<Number,OtherNumber>::type> tt = t;
3614  tt /= factor;
3615  return tt;
3616 }
3617 
3618 
3619 
3626 template <int rank_, int dim>
3627 inline
3628 SymmetricTensor<rank_,dim>
3629 operator * (const SymmetricTensor<rank_,dim> &t,
3630  const double factor)
3631 {
3632  SymmetricTensor<rank_,dim> tt = t;
3633  tt *= factor;
3634  return tt;
3635 }
3636 
3637 
3638 
3645 template <int rank_, int dim>
3646 inline
3647 SymmetricTensor<rank_,dim>
3648 operator * (const double factor,
3649  const SymmetricTensor<rank_,dim> &t)
3650 {
3651  SymmetricTensor<rank_,dim> tt = t;
3652  tt *= factor;
3653  return tt;
3654 }
3655 
3656 
3657 
3663 template <int rank_, int dim>
3664 inline
3665 SymmetricTensor<rank_,dim>
3666 operator / (const SymmetricTensor<rank_,dim> &t,
3667  const double factor)
3668 {
3669  SymmetricTensor<rank_,dim> tt = t;
3670  tt /= factor;
3671  return tt;
3672 }
3673 
3683 template <int dim, typename Number, typename OtherNumber>
3684 inline
3685 typename ProductType<Number, OtherNumber>::type
3686 scalar_product (const SymmetricTensor<2,dim,Number> &t1,
3687  const SymmetricTensor<2,dim,OtherNumber> &t2)
3688 {
3689  return (t1*t2);
3690 }
3691 
3692 
3702 template <int dim, typename Number, typename OtherNumber>
3703 inline
3704 typename ProductType<Number, OtherNumber>::type
3705 scalar_product (const SymmetricTensor<2,dim,Number> &t1,
3706  const Tensor<2,dim,OtherNumber> &t2)
3707 {
3708  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<typename ProductType<Number, OtherNumber>::type>::value(0.0);
3709  for (unsigned int i=0; i<dim; ++i)
3710  for (unsigned int j=0; j<dim; ++j)
3711  s += t1[i][j] * t2[i][j];
3712  return s;
3713 }
3714 
3715 
3725 template <int dim, typename Number, typename OtherNumber>
3726 inline
3727 typename ProductType<Number, OtherNumber>::type
3728 scalar_product (const Tensor<2,dim,Number> &t1,
3729  const SymmetricTensor<2,dim,OtherNumber> &t2)
3730 {
3731  return scalar_product(t2, t1);
3732 }
3733 
3734 
3750 template <typename Number, typename OtherNumber>
3751 inline
3752 void
3753 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3754  const SymmetricTensor<4,1,Number> &t,
3755  const SymmetricTensor<2,1,OtherNumber> &s)
3756 {
3757  tmp[0][0] = t[0][0][0][0] * s[0][0];
3758 }
3759 
3760 
3761 
3777 template <typename Number, typename OtherNumber>
3778 inline
3779 void
3780 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3781  const SymmetricTensor<2,1,Number> &s,
3782  const SymmetricTensor<4,1,OtherNumber> &t)
3783 {
3784  tmp[0][0] = t[0][0][0][0] * s[0][0];
3785 }
3786 
3787 
3788 
3803 template <typename Number, typename OtherNumber>
3804 inline
3805 void
3806 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3807  const SymmetricTensor<4,2,Number> &t,
3808  const SymmetricTensor<2,2,OtherNumber> &s)
3809 {
3810  const unsigned int dim = 2;
3811 
3812  for (unsigned int i=0; i<dim; ++i)
3813  for (unsigned int j=i; j<dim; ++j)
3814  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3815  t[i][j][1][1] * s[1][1] +
3816  2 * t[i][j][0][1] * s[0][1];
3817 }
3818 
3819 
3820 
3836 template <typename Number, typename OtherNumber>
3837 inline
3838 void
3839 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3840  const SymmetricTensor<2,2,Number> &s,
3841  const SymmetricTensor<4,2,OtherNumber> &t)
3842 {
3843  const unsigned int dim = 2;
3844 
3845  for (unsigned int i=0; i<dim; ++i)
3846  for (unsigned int j=i; j<dim; ++j)
3847  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3848  s[1][1] * t[1][1][i][j] +
3849  2 * s[0][1] * t[0][1][i][j];
3850 }
3851 
3852 
3853 
3869 template <typename Number, typename OtherNumber>
3870 inline
3871 void
3872 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3873  const SymmetricTensor<4,3,Number> &t,
3874  const SymmetricTensor<2,3,OtherNumber> &s)
3875 {
3876  const unsigned int dim = 3;
3877 
3878  for (unsigned int i=0; i<dim; ++i)
3879  for (unsigned int j=i; j<dim; ++j)
3880  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3881  t[i][j][1][1] * s[1][1] +
3882  t[i][j][2][2] * s[2][2] +
3883  2 * t[i][j][0][1] * s[0][1] +
3884  2 * t[i][j][0][2] * s[0][2] +
3885  2 * t[i][j][1][2] * s[1][2];
3886 }
3887 
3888 
3889 
3905 template <typename Number, typename OtherNumber>
3906 inline
3907 void
3908 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3909  const SymmetricTensor<2,3,Number> &s,
3910  const SymmetricTensor<4,3,OtherNumber> &t)
3911 {
3912  const unsigned int dim = 3;
3913 
3914  for (unsigned int i=0; i<dim; ++i)
3915  for (unsigned int j=i; j<dim; ++j)
3916  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3917  s[1][1] * t[1][1][i][j] +
3918  s[2][2] * t[2][2][i][j] +
3919  2 * s[0][1] * t[0][1][i][j] +
3920  2 * s[0][2] * t[0][2][i][j] +
3921  2 * s[1][2] * t[1][2][i][j];
3922 }
3923 
3924 
3925 
3933 template <int dim, typename Number, typename OtherNumber>
3934 Tensor<1,dim,typename ProductType<Number,OtherNumber>::type>
3935 operator * (const SymmetricTensor<2,dim,Number> &src1,
3936  const Tensor<1,dim,OtherNumber> &src2)
3937 {
3938  Tensor<1,dim,typename ProductType<Number,OtherNumber>::type> dest;
3939  for (unsigned int i=0; i<dim; ++i)
3940  for (unsigned int j=0; j<dim; ++j)
3941  dest[i] += src1[i][j] * src2[j];
3942  return dest;
3943 }
3944 
3945 
3953 template <int dim, typename Number, typename OtherNumber>
3954 Tensor<1,dim,typename ProductType<Number,OtherNumber>::type>
3955 operator * (const Tensor<1,dim,Number> &src1,
3956  const SymmetricTensor<2,dim,OtherNumber> &src2)
3957 {
3958  // this is easy for symmetric tensors:
3959  return src2 * src1;
3960 }
3961 
3962 
3963 
3984 template <int rank_1, int rank_2, int dim,
3985  typename Number, typename OtherNumber>
3986 inline DEAL_II_ALWAYS_INLINE
3987 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3988 operator * (const Tensor<rank_1, dim, Number> &src1,
3989  const SymmetricTensor<rank_2, dim, OtherNumber> &src2s)
3990 {
3991  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3992  const Tensor<rank_2, dim, OtherNumber> src2 (src2s);
3993  return src1*src2;
3994 }
3995 
3996 
3997 
4018 template <int rank_1, int rank_2, int dim,
4019  typename Number, typename OtherNumber>
4020 inline DEAL_II_ALWAYS_INLINE
4021 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
4022 operator * (const SymmetricTensor<rank_1, dim, Number> &src1s,
4023  const Tensor<rank_2, dim, OtherNumber> &src2)
4024 {
4025  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
4026  const Tensor<rank_2, dim, OtherNumber> src1 (src1s);
4027  return src1*src2;
4028 }
4029 
4030 
4031 
4041 template <int dim, typename Number>
4042 inline
4043 std::ostream &operator << (std::ostream &out,
4044  const SymmetricTensor<2,dim,Number> &t)
4045 {
4046  //make our lives a bit simpler by outputting
4047  //the tensor through the operator for the
4048  //general Tensor class
4049  Tensor<2,dim,Number> tt;
4050 
4051  for (unsigned int i=0; i<dim; ++i)
4052  for (unsigned int j=0; j<dim; ++j)
4053  tt[i][j] = t[i][j];
4054 
4055  return out << tt;
4056 }
4057 
4058 
4059 
4069 template <int dim, typename Number>
4070 inline
4071 std::ostream &operator << (std::ostream &out,
4072  const SymmetricTensor<4,dim,Number> &t)
4073 {
4074  //make our lives a bit simpler by outputting
4075  //the tensor through the operator for the
4076  //general Tensor class
4077  Tensor<4,dim,Number> tt;
4078 
4079  for (unsigned int i=0; i<dim; ++i)
4080  for (unsigned int j=0; j<dim; ++j)
4081  for (unsigned int k=0; k<dim; ++k)
4082  for (unsigned int l=0; l<dim; ++l)
4083  tt[i][j][k][l] = t[i][j][k][l];
4084 
4085  return out << tt;
4086 }
4087 
4088 
4089 DEAL_II_NAMESPACE_CLOSE
4090 
4091 #endif
SymmetricTensorEigenvectorMethod::jacobi
SymmetricTensor::identity_tensor
friend SymmetricTensor< 4, dim2, Number2 > identity_tensor()
Definition: symmetric_tensor.h:3337
SymmetricTensor::operator/
SymmetricTensor< rank_, dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator/(const SymmetricTensor< rank_, dim, Number > &t, const OtherNumber &factor)
Definition: symmetric_tensor.h:3610
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
SymmetricTensor::component_to_unrolled_index
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensorEigenvectorMethod::hybrid
SymmetricTensor::invert
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3404
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor::eigenvectors
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
SymmetricTensor::n_independent_components
static const unsigned int n_independent_components
Definition: symmetric_tensor.h:564
SymmetricTensor::symmetrize
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3477
internal::SymmetricTensorAccessors::double_contraction_result
Definition: symmetric_tensor.h:152
numbers::value_is_zero
bool value_is_zero(const Number &value)
Definition: numbers.h:784
SymmetricTensor::double_contract
void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
Definition: symmetric_tensor.h:3753
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1284
deviator_tensor
SymmetricTensor< 4, dim, Number > deviator_tensor()
Definition: symmetric_tensor.h:3264
internal::SymmetricTensorAccessors::merge
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
Definition: symmetric_tensor.h:88
SymmetricTensor::operator=
SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
SymmetricTensor::memory_consumption
static std::size_t memory_consumption()
numbers::NumberTraits::abs
static real_type abs(const number &x)
Definition: numbers.h:465
SymmetricTensorEigenvectorMethod
SymmetricTensorEigenvectorMethod
Definition: symmetric_tensor.h:3081
SymmetricTensor::operator+
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Definition: symmetric_tensor.h:2500
SymmetricTensor::base_tensor_descriptor
internal::SymmetricTensorAccessors::StorageType< rank_, dim, Number > base_tensor_descriptor
Definition: symmetric_tensor.h:856
StandardExceptions::ExcIndexRange
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor::operator/=
SymmetricTensor & operator/=(const OtherNumber &factor)
SymmetricTensor::unrolled_to_component_indices
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
SymmetricTensor::operator*
internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
SymmetricTensor::dimension
static const unsigned int dimension
Definition: symmetric_tensor.h:551
SymmetricTensor::operator-
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Definition: symmetric_tensor.h:2524
internal::SymmetricTensorImplementation::ql_implicit_shifts
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > ql_implicit_shifts(const ::SymmetricTensor< 2, dim, Number > &A)
internal::SymmetricTensorImplementation::jacobi
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::SymmetricTensorImplementation::tridiagonalize
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A, ::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim-1 > &e)
internal::SymmetricTensorAccessors::merge
TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
Definition: symmetric_tensor.h:109
SymmetricTensor::first_invariant
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2691
SymmetricTensor::access_raw_entry
const Number & access_raw_entry(const unsigned int unrolled_index) const
SymmetricTensor::base_tensor_type
base_tensor_descriptor::base_tensor_type base_tensor_type
Definition: symmetric_tensor.h:861
SymmetricTensor::trace
friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
SymmetricTensor::rank
static const unsigned int rank
Definition: symmetric_tensor.h:556
SymmetricTensor::operator!=
bool operator!=(const SymmetricTensor &) const
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
SymmetricTensor::begin_raw
Number * begin_raw()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1142
SymmetricTensor::data
base_tensor_type data
Definition: symmetric_tensor.h:866
SymmetricTensor::deviator
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3180
SymmetricTensor::trace
Number trace(const SymmetricTensor< 2, dim, Number > &d)
Definition: symmetric_tensor.h:2671
SymmetricTensor::serialize
void serialize(Archive &ar, const unsigned int version)
SymmetricTensor::end_raw
Number * end_raw()
SymmetricTensor::unit_symmetric_tensor
friend SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
Definition: symmetric_tensor.h:3204
SymmetricTensor::outer_product
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition: symmetric_tensor.h:3449
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor::scalar_product
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Definition: symmetric_tensor.h:3686
SymmetricTensorEigenvectorMethod::ql_implicit_shifts
SymmetricTensor::operator()
Number & operator()(const TableIndices< rank_ > &indices)
internal::SymmetricTensorImplementation::hybrid
std::array< std::pair< Number, Tensor< 1, 2, Number > >, 2 > hybrid(const ::SymmetricTensor< 2, 2, Number > &A)
swap
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1312
Algorithms::OutputOperator::operator<<
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:163
SymmetricTensor::determinant
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2617
SymmetricTensor::operator-
SymmetricTensor operator-() const
Tensor
Definition: mpi.h:53
SymmetricTensor::operator*=
SymmetricTensor & operator*=(const OtherNumber &factor)
SymmetricTensor::operator+=
SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
SymmetricTensor::third_invariant
Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2656
internal::SymmetricTensorAccessors::StorageType< 2, dim, Number >::base_tensor_type
Tensor< 1, n_independent_components, Number > base_tensor_type
Definition: symmetric_tensor.h:206
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
SymmetricTensor::operator==
bool operator==(const SymmetricTensor &) const
SymmetricTensor::eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
unit_symmetric_tensor
SymmetricTensor< 2, dim, Number > unit_symmetric_tensor()
Definition: symmetric_tensor.h:3204
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:463
SymmetricTensor::deviator_tensor
friend SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
Definition: symmetric_tensor.h:3264
Point::operator*
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
SymmetricTensor::operator[]
internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_-1, Number > operator[](const unsigned int row) const
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
identity_tensor
SymmetricTensor< 4, dim, Number > identity_tensor()
Definition: symmetric_tensor.h:3337
SymmetricTensor::second_invariant
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Definition: symmetric_tensor.h:2711
SymmetricTensor::operator-=
SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
SymmetricTensor::norm
numbers::NumberTraits< Number >::real_type norm() const
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Generated by   doxygen 1.8.14