Reference documentation for deal.II version Git 58efed1 2018-02-20 08:43:45 -0500
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>
35 
36 template <int dim, typename Number>
38 
39 template <int dim, typename Number>
41 
42 template <int dim, typename Number>
44 
45 template <int dim, typename Number>
47 
48 template <int dim2, typename Number> Number
50 
51 template <int dim, typename Number>
53 
54 template <int dim, typename Number>
56 
57 
58 
59 namespace internal
60 {
65  namespace SymmetricTensor
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)
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,
122  case 1:
123  return TableIndices<4>(previous_indices[0],
124  new_index,
127  case 2:
128  return TableIndices<4>(previous_indices[0],
129  previous_indices[1],
130  new_index,
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>
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 
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 *
231 
239  };
240 
241 
242 
247  template <int rank, int dim, bool constness, typename Number>
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
566 
570  SymmetricTensor ();
571 
582  template <typename OtherNumber>
583  explicit
585 
601  SymmetricTensor (const Number (&array) [n_independent_components]);
602 
608  template <typename OtherNumber>
609  explicit
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>
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>
674 
678  template <typename OtherNumber>
680 
685  template <typename OtherNumber>
686  SymmetricTensor &operator *= (const OtherNumber &factor);
687 
691  template <typename OtherNumber>
692  SymmetricTensor &operator /= (const OtherNumber &factor);
693 
698 
723  template <typename OtherNumber>
726 
731  template <typename OtherNumber>
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 
801  norm () const;
802 
810  static
811  unsigned int
813 
819  static
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
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>
884  deviator (const SymmetricTensor<2,dim2,Number2> &t);
885 
886  template <int dim2, typename Number2>
888 
889  template <int dim2, typename Number2>
891 
892  template <int dim2, typename Number2>
894 
895 
899  friend struct internal::SymmetricTensor::Inverse<2,dim,Number>;
900 
901  friend struct internal::SymmetricTensor::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
978 {
979  // Some auto-differentiable numbers need explicit
980  // zero initialization.
981  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
983 }
984 
985 
986 template <int rank_, int dim, typename Number>
987 template <typename OtherNumber>
988 inline
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
1036 {
1037  for (unsigned int i=0; i<base_tensor_type::dimension; ++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
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
1075 {
1076  Assert (d==internal::NumberType<Number>::value(0.0), ExcMessage ("Only assignment with zero is allowed"));
1077  (void) d;
1078 
1080 
1081  return *this;
1082 }
1083 
1084 
1085 namespace internal
1086 {
1087  namespace SymmetricTensor
1088  {
1089  template <int dim, typename Number>
1091  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1092  {
1094 
1095  // diagonal entries are stored first
1096  for (unsigned int d=0; d<dim; ++d)
1097  t[d][d] = s.access_raw_entry(d);
1098 
1099  // off-diagonal entries come next, row by row
1100  for (unsigned int d=0, c=0; d<dim; ++d)
1101  for (unsigned int e=d+1; e<dim; ++e, ++c)
1102  {
1103  t[d][e] = s.access_raw_entry(dim+c);
1104  t[e][d] = s.access_raw_entry(dim+c);
1105  }
1106  return t;
1107  }
1108 
1109 
1110  template <int dim, typename Number>
1112  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1113  {
1114  // utilize the symmetry properties of SymmetricTensor<4,dim>
1115  // discussed in the class documentation to avoid accessing all
1116  // independent elements of the input tensor more than once
1118 
1119  for (unsigned int i=0; i<dim; ++i)
1120  for (unsigned int j=i; j<dim; ++j)
1121  for (unsigned int k=0; k<dim; ++k)
1122  for (unsigned int l=k; l<dim; ++l)
1123  t[TableIndices<4>(i,j,k,l)]
1124  = t[TableIndices<4>(i,j,l,k)]
1125  = t[TableIndices<4>(j,i,k,l)]
1126  = t[TableIndices<4>(j,i,l,k)]
1127  = st[TableIndices<4>(i,j,k,l)];
1128 
1129  return t;
1130  }
1131 
1132 
1133  template <typename Number>
1134  struct Inverse<2,1,Number>
1135  {
1136  static inline ::SymmetricTensor<2,1,Number>
1137  value (const ::SymmetricTensor<2,1,Number> &t)
1138  {
1140 
1141  tmp[0][0] = 1.0/t[0][0];
1142 
1143  return tmp;
1144  }
1145  };
1146 
1147 
1148  template <typename Number>
1149  struct Inverse<2,2,Number>
1150  {
1151  static inline ::SymmetricTensor<2,2,Number>
1152  value (const ::SymmetricTensor<2,2,Number> &t)
1153  {
1155 
1156  // Sympy result: ([
1157  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1158  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1159  const TableIndices<2> idx_00 (0,0);
1160  const TableIndices<2> idx_01 (0,1);
1161  const TableIndices<2> idx_11 (1,1);
1162  const Number inv_det_t
1163  = 1.0/(t[idx_00]*t[idx_11]
1164  - t[idx_01]*t[idx_01]);
1165  tmp[idx_00] = t[idx_11];
1166  tmp[idx_01] = -t[idx_01];
1167  tmp[idx_11] = t[idx_00];
1168  tmp *= inv_det_t;
1169 
1170  return tmp;
1171  }
1172  };
1173 
1174 
1175  template <typename Number>
1176  struct Inverse<2,3,Number>
1177  {
1178  static ::SymmetricTensor<2,3,Number>
1179  value (const ::SymmetricTensor<2,3,Number> &t)
1180  {
1182 
1183  // Sympy result: ([
1184  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1185  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1186  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
1187  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1188  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1189  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
1190  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1191  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
1192  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
1193  const TableIndices<2> idx_00 (0,0);
1194  const TableIndices<2> idx_01 (0,1);
1195  const TableIndices<2> idx_02 (0,2);
1196  const TableIndices<2> idx_11 (1,1);
1197  const TableIndices<2> idx_12 (1,2);
1198  const TableIndices<2> idx_22 (2,2);
1199  const Number inv_det_t
1200  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
1201  - t[idx_00]*t[idx_12]*t[idx_12]
1202  - t[idx_01]*t[idx_01]*t[idx_22]
1203  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
1204  - t[idx_02]*t[idx_02]*t[idx_11]);
1205  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
1206  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
1207  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
1208  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
1209  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
1210  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
1211  tmp *= inv_det_t;
1212 
1213  return tmp;
1214  }
1215  };
1216 
1217 
1218  template <typename Number>
1219  struct Inverse<4,1,Number>
1220  {
1221  static inline ::SymmetricTensor<4,1,Number>
1222  value (const ::SymmetricTensor<4,1,Number> &t)
1223  {
1225  tmp.data[0][0] = 1.0/t.data[0][0];
1226  return tmp;
1227  }
1228  };
1229 
1230 
1231  template <typename Number>
1232  struct Inverse<4,2,Number>
1233  {
1234  static inline ::SymmetricTensor<4,2,Number>
1235  value (const ::SymmetricTensor<4,2,Number> &t)
1236  {
1238 
1239  // Inverting this tensor is a little more complicated than necessary,
1240  // since we store the data of 't' as a 3x3 matrix t.data, but the
1241  // product between a rank-4 and a rank-2 tensor is really not the
1242  // product between this matrix and the 3-vector of a rhs, but rather
1243  //
1244  // B.vec = t.data * mult * A.vec
1245  //
1246  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1247  // capture the fact that we need to add up both the c_ij12*a_12 and the
1248  // c_ij21*a_21 terms.
1249  //
1250  // In addition, in this scheme, the identity tensor has the matrix
1251  // representation mult^-1.
1252  //
1253  // The inverse of 't' therefore has the matrix representation
1254  //
1255  // inv.data = mult^-1 * t.data^-1 * mult^-1
1256  //
1257  // in order to compute it, let's first compute the inverse of t.data and
1258  // put it into tmp.data; at the end of the function we then scale the
1259  // last row and column of the inverse by 1/2, corresponding to the left
1260  // and right multiplication with mult^-1.
1261  const Number t4 = t.data[0][0]*t.data[1][1],
1262  t6 = t.data[0][0]*t.data[1][2],
1263  t8 = t.data[0][1]*t.data[1][0],
1264  t00 = t.data[0][2]*t.data[1][0],
1265  t01 = t.data[0][1]*t.data[2][0],
1266  t04 = t.data[0][2]*t.data[2][0],
1267  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
1268  t8*t.data[2][2]+t00*t.data[2][1]+
1269  t01*t.data[1][2]-t04*t.data[1][1]);
1270  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
1271  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
1272  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
1273  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
1274  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
1275  tmp.data[1][2] = -(t6-t00)*t07;
1276  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
1277  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
1278  tmp.data[2][2] = (t4-t8)*t07;
1279 
1280  // scale last row and column as mentioned
1281  // above
1282  tmp.data[2][0] /= 2;
1283  tmp.data[2][1] /= 2;
1284  tmp.data[0][2] /= 2;
1285  tmp.data[1][2] /= 2;
1286  tmp.data[2][2] /= 4;
1287 
1288  return tmp;
1289  }
1290  };
1291 
1292 
1293  template <typename Number>
1294  struct Inverse<4,3,Number>
1295  {
1296  static ::SymmetricTensor<4,3,Number>
1297  value (const ::SymmetricTensor<4,3,Number> &t)
1298  {
1300 
1301  // This function follows the exact same scheme as the 2d case, except
1302  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1303  // Instead, we use the Gauss-Jordan algorithm implemented for
1304  // FullMatrix. For historical reasons the following code is copied from
1305  // there, with the tangential benefit that we do not need to copy the
1306  // tensor entries to and from the FullMatrix.
1307  const unsigned int N = 6;
1308 
1309  // First get an estimate of the size of the elements of this matrix,
1310  // for later checks whether the pivot element is large enough, or
1311  // whether we have to fear that the matrix is not regular.
1312  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1313  for (unsigned int i=0; i<N; ++i)
1314  diagonal_sum += std::fabs(tmp.data[i][i]);
1315  const Number typical_diagonal_element = diagonal_sum/N;
1316  (void)typical_diagonal_element;
1317 
1318  unsigned int p[N];
1319  for (unsigned int i=0; i<N; ++i)
1320  p[i] = i;
1321 
1322  for (unsigned int j=0; j<N; ++j)
1323  {
1324  // Pivot search: search that part of the line on and right of the
1325  // diagonal for the largest element.
1326  Number max = std::fabs(tmp.data[j][j]);
1327  unsigned int r = j;
1328  for (unsigned int i=j+1; i<N; ++i)
1329  if (std::fabs(tmp.data[i][j]) > max)
1330  {
1331  max = std::fabs(tmp.data[i][j]);
1332  r = i;
1333  }
1334 
1335  // Check whether the pivot is too small
1336  Assert(max > 1.e-16*typical_diagonal_element,
1337  ExcMessage("This tensor seems to be noninvertible"));
1338 
1339  // Row interchange
1340  if (r>j)
1341  {
1342  for (unsigned int k=0; k<N; ++k)
1343  std::swap (tmp.data[j][k], tmp.data[r][k]);
1344 
1345  std::swap (p[j], p[r]);
1346  }
1347 
1348  // Transformation
1349  const Number hr = 1./tmp.data[j][j];
1350  tmp.data[j][j] = hr;
1351  for (unsigned int k=0; k<N; ++k)
1352  {
1353  if (k==j) continue;
1354  for (unsigned int i=0; i<N; ++i)
1355  {
1356  if (i==j) continue;
1357  tmp.data[i][k] -= tmp.data[i][j]*tmp.data[j][k]*hr;
1358  }
1359  }
1360  for (unsigned int i=0; i<N; ++i)
1361  {
1362  tmp.data[i][j] *= hr;
1363  tmp.data[j][i] *= -hr;
1364  }
1365  tmp.data[j][j] = hr;
1366  }
1367 
1368  // Column interchange
1369  Number hv[N];
1370  for (unsigned int i=0; i<N; ++i)
1371  {
1372  for (unsigned int k=0; k<N; ++k)
1373  hv[p[k]] = tmp.data[i][k];
1374  for (unsigned int k=0; k<N; ++k)
1375  tmp.data[i][k] = hv[k];
1376  }
1377 
1378  // Scale rows and columns. The mult matrix
1379  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1380  for (unsigned int i=3; i<6; ++i)
1381  for (unsigned int j=0; j<3; ++j)
1382  tmp.data[i][j] /= 2;
1383 
1384  for (unsigned int i=0; i<3; ++i)
1385  for (unsigned int j=3; j<6; ++j)
1386  tmp.data[i][j] /= 2;
1387 
1388  for (unsigned int i=3; i<6; ++i)
1389  for (unsigned int j=3; j<6; ++j)
1390  tmp.data[i][j] /= 4;
1391 
1392  return tmp;
1393  }
1394  };
1395 
1396  }
1397 }
1398 
1399 
1400 
1401 template <int rank_, int dim, typename Number>
1402 inline
1404 operator Tensor<rank_,dim,Number> () const
1405 {
1406  return internal::SymmetricTensor::convert_to_tensor (*this);
1407 }
1408 
1409 
1410 
1411 template <int rank_, int dim, typename Number>
1412 inline
1413 bool
1415 (const SymmetricTensor<rank_,dim,Number> &t) const
1416 {
1417  return data == t.data;
1418 }
1419 
1420 
1421 
1422 template <int rank_, int dim, typename Number>
1423 inline
1424 bool
1425 SymmetricTensor<rank_,dim,Number>::operator !=
1426 (const SymmetricTensor<rank_,dim,Number> &t) const
1427 {
1428  return data != t.data;
1429 }
1430 
1431 
1432 
1433 template <int rank_, int dim, typename Number>
1434 template <typename OtherNumber>
1435 inline
1437 SymmetricTensor<rank_,dim,Number>::operator +=
1439 {
1440  data += t.data;
1441  return *this;
1442 }
1443 
1444 
1445 
1446 template <int rank_, int dim, typename Number>
1447 template <typename OtherNumber>
1448 inline
1450 SymmetricTensor<rank_,dim,Number>::operator -=
1452 {
1453  data -= t.data;
1454  return *this;
1455 }
1456 
1457 
1458 
1459 template <int rank_, int dim, typename Number>
1460 template <typename OtherNumber>
1461 inline
1464 {
1465  data *= d;
1466  return *this;
1467 }
1468 
1469 
1470 
1471 template <int rank_, int dim, typename Number>
1472 template <typename OtherNumber>
1473 inline
1476 {
1477  data /= d;
1478  return *this;
1479 }
1480 
1481 
1482 
1483 template <int rank_, int dim, typename Number>
1484 inline
1487 {
1488  SymmetricTensor tmp = *this;
1489  tmp.data = -tmp.data;
1490  return tmp;
1491 }
1492 
1493 
1494 
1495 template <int rank_, int dim, typename Number>
1496 inline
1497 void
1499 {
1500  data.clear ();
1501 }
1502 
1503 
1504 
1505 template <int rank_, int dim, typename Number>
1506 inline
1507 std::size_t
1509 {
1510  // all memory consists of statically allocated memory of the current
1511  // object, no pointers
1512  return sizeof(SymmetricTensor<rank_,dim,Number>);
1513 }
1514 
1515 
1516 
1517 namespace internal
1518 {
1519 
1520  template <int dim, typename Number, typename OtherNumber = Number>
1521  inline
1522  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type
1523  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1524  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1525  {
1526  typedef typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type result_type;
1527 
1528  switch (dim)
1529  {
1530  case 1:
1531  return data[0] * sdata[0];
1532  default:
1533  // Start with the non-diagonal part to avoid some multiplications by
1534  // 2.
1535 
1536  result_type sum = data[dim] * sdata[dim];
1537  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1538  sum += data[d] * sdata[d];
1539  sum += sum; // sum = sum * 2.;
1540 
1541  // Now add the contributions from the diagonal
1542  for (unsigned int d=0; d<dim; ++d)
1543  sum += data[d] * sdata[d];
1544  return sum;
1545  }
1546  }
1547 
1548 
1549 
1550  template <int dim, typename Number, typename OtherNumber = Number>
1551  inline
1552  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type
1553  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1554  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1555  {
1556  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type result_type;
1557  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::value_type value_type;
1558 
1559  const unsigned int data_dim =
1560  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1561  value_type tmp [data_dim];
1562  for (unsigned int i=0; i<data_dim; ++i)
1563  tmp[i] = perform_double_contraction<dim,Number,OtherNumber>(data[i], sdata);
1564  return result_type(tmp);
1565  }
1566 
1567 
1568 
1569  template <int dim, typename Number, typename OtherNumber = Number>
1570  inline
1571  typename SymmetricTensorAccessors::StorageType<2,dim,
1572  typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type
1573  >::base_tensor_type
1574  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1575  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1576  {
1577  typedef typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type value_type;
1578  typedef typename SymmetricTensorAccessors::StorageType<2,dim,value_type>::base_tensor_type base_tensor_type;
1579 
1580  base_tensor_type tmp;
1581  for (unsigned int i=0; i<tmp.dimension; ++i)
1582  {
1583  // Start with the non-diagonal part
1584  value_type sum = data[dim] * sdata[dim][i];
1585  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1586  sum += data[d] * sdata[d][i];
1587  sum += sum; // sum = sum * 2.;
1588 
1589  // Now add the contributions from the diagonal
1590  for (unsigned int d=0; d<dim; ++d)
1591  sum += data[d] * sdata[d][i];
1592  tmp[i] = sum;
1593  }
1594  return tmp;
1595  }
1596 
1597 
1598 
1599  template <int dim, typename Number, typename OtherNumber = Number>
1600  inline
1601  typename SymmetricTensorAccessors::StorageType<4,dim,
1602  typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type
1603  >::base_tensor_type
1604  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1605  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1606  {
1607  typedef typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type value_type;
1608  typedef typename SymmetricTensorAccessors::StorageType<4,dim,value_type>::base_tensor_type base_tensor_type;
1609 
1610  const unsigned int data_dim =
1611  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1612  base_tensor_type tmp;
1613  for (unsigned int i=0; i<data_dim; ++i)
1614  for (unsigned int j=0; j<data_dim; ++j)
1615  {
1616  // Start with the non-diagonal part
1617  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1618  tmp[i][j] += data[i][d] * sdata[d][j];
1619  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1620 
1621  // Now add the contributions from the diagonal
1622  for (unsigned int d=0; d<dim; ++d)
1623  tmp[i][j] += data[i][d] * sdata[d][j];
1624  }
1625  return tmp;
1626  }
1627 
1628 } // end of namespace internal
1629 
1630 
1631 
1632 template <int rank_, int dim, typename Number>
1633 template <typename OtherNumber>
1634 inline
1637 {
1638  // need to have two different function calls
1639  // because a scalar and rank-2 tensor are not
1640  // the same data type (see internal function
1641  // above)
1642  return internal::perform_double_contraction<dim,Number,OtherNumber> (data, s.data);
1643 }
1644 
1645 
1646 
1647 template <int rank_, int dim, typename Number>
1648 template <typename OtherNumber>
1649 inline
1652 {
1655  tmp.data = internal::perform_double_contraction<dim,Number,OtherNumber> (data,s.data);
1656  return tmp;
1657 }
1658 
1659 
1660 
1661 // internal namespace to switch between the
1662 // access of different tensors. There used to
1663 // be explicit instantiations before for
1664 // different ranks and dimensions, but since
1665 // we now allow for templates on the data
1666 // type, and since we cannot partially
1667 // specialize the implementation, this got
1668 // into a separate namespace
1669 namespace internal
1670 {
1671  template <int dim, typename Number>
1672  inline
1673  Number &
1674  symmetric_tensor_access (const TableIndices<2> &indices,
1676  {
1677  // 1d is very simple and done first
1678  if (dim == 1)
1679  return data[0];
1680 
1681  // first treat the main diagonal elements, which are stored consecutively
1682  // at the beginning
1683  if (indices[0] == indices[1])
1684  return data[indices[0]];
1685 
1686  // the rest is messier and requires a few switches.
1687  switch (dim)
1688  {
1689  case 2:
1690  // at least for the 2x2 case it is reasonably simple
1691  Assert (((indices[0]==1) && (indices[1]==0)) ||
1692  ((indices[0]==0) && (indices[1]==1)),
1693  ExcInternalError());
1694  return data[2];
1695 
1696  default:
1697  // to do the rest, sort our indices before comparing
1698  {
1699  TableIndices<2> sorted_indices (indices);
1700  sorted_indices.sort ();
1701 
1702  for (unsigned int d=0, c=0; d<dim; ++d)
1703  for (unsigned int e=d+1; e<dim; ++e, ++c)
1704  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1705  return data[dim+c];
1706  Assert (false, ExcInternalError());
1707  }
1708  }
1709 
1710  static Number dummy_but_referenceable = Number();
1711  return dummy_but_referenceable;
1712  }
1713 
1714 
1715 
1716  template <int dim, typename Number>
1717  inline
1718  const Number &
1719  symmetric_tensor_access (const TableIndices<2> &indices,
1721  {
1722  // 1d is very simple and done first
1723  if (dim == 1)
1724  return data[0];
1725 
1726  // first treat the main diagonal elements, which are stored consecutively
1727  // at the beginning
1728  if (indices[0] == indices[1])
1729  return data[indices[0]];
1730 
1731  // the rest is messier and requires a few switches.
1732  switch (dim)
1733  {
1734  case 2:
1735  // at least for the 2x2 case it is reasonably simple
1736  Assert (((indices[0]==1) && (indices[1]==0)) ||
1737  ((indices[0]==0) && (indices[1]==1)),
1738  ExcInternalError());
1739  return data[2];
1740 
1741  default:
1742  // to do the rest, sort our indices before comparing
1743  {
1744  TableIndices<2> sorted_indices (indices);
1745  sorted_indices.sort ();
1746 
1747  for (unsigned int d=0, c=0; d<dim; ++d)
1748  for (unsigned int e=d+1; e<dim; ++e, ++c)
1749  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1750  return data[dim+c];
1751  Assert (false, ExcInternalError());
1752  }
1753  }
1754 
1755  static Number dummy_but_referenceable = Number();
1756  return dummy_but_referenceable;
1757  }
1758 
1759 
1760 
1761  template <int dim, typename Number>
1762  inline
1763  Number &
1764  symmetric_tensor_access (const TableIndices<4> &indices,
1766  {
1767  switch (dim)
1768  {
1769  case 1:
1770  return data[0][0];
1771 
1772  case 2:
1773  // each entry of the tensor can be
1774  // thought of as an entry in a
1775  // matrix that maps the rolled-out
1776  // rank-2 tensors into rolled-out
1777  // rank-2 tensors. this is the
1778  // format in which we store rank-4
1779  // tensors. determine which
1780  // position the present entry is
1781  // stored in
1782  {
1783  unsigned int base_index[2] ;
1784  if ((indices[0] == 0) && (indices[1] == 0))
1785  base_index[0] = 0;
1786  else if ((indices[0] == 1) && (indices[1] == 1))
1787  base_index[0] = 1;
1788  else
1789  base_index[0] = 2;
1790 
1791  if ((indices[2] == 0) && (indices[3] == 0))
1792  base_index[1] = 0;
1793  else if ((indices[2] == 1) && (indices[3] == 1))
1794  base_index[1] = 1;
1795  else
1796  base_index[1] = 2;
1797 
1798  return data[base_index[0]][base_index[1]];
1799  }
1800 
1801  case 3:
1802  // each entry of the tensor can be
1803  // thought of as an entry in a
1804  // matrix that maps the rolled-out
1805  // rank-2 tensors into rolled-out
1806  // rank-2 tensors. this is the
1807  // format in which we store rank-4
1808  // tensors. determine which
1809  // position the present entry is
1810  // stored in
1811  {
1812  unsigned int base_index[2] ;
1813  if ((indices[0] == 0) && (indices[1] == 0))
1814  base_index[0] = 0;
1815  else if ((indices[0] == 1) && (indices[1] == 1))
1816  base_index[0] = 1;
1817  else if ((indices[0] == 2) && (indices[1] == 2))
1818  base_index[0] = 2;
1819  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1820  ((indices[0] == 1) && (indices[1] == 0)))
1821  base_index[0] = 3;
1822  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1823  ((indices[0] == 2) && (indices[1] == 0)))
1824  base_index[0] = 4;
1825  else
1826  {
1827  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1828  ((indices[0] == 2) && (indices[1] == 1)),
1829  ExcInternalError());
1830  base_index[0] = 5;
1831  }
1832 
1833  if ((indices[2] == 0) && (indices[3] == 0))
1834  base_index[1] = 0;
1835  else if ((indices[2] == 1) && (indices[3] == 1))
1836  base_index[1] = 1;
1837  else if ((indices[2] == 2) && (indices[3] == 2))
1838  base_index[1] = 2;
1839  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1840  ((indices[2] == 1) && (indices[3] == 0)))
1841  base_index[1] = 3;
1842  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1843  ((indices[2] == 2) && (indices[3] == 0)))
1844  base_index[1] = 4;
1845  else
1846  {
1847  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1848  ((indices[2] == 2) && (indices[3] == 1)),
1849  ExcInternalError());
1850  base_index[1] = 5;
1851  }
1852 
1853  return data[base_index[0]][base_index[1]];
1854  }
1855 
1856  default:
1857  Assert (false, ExcNotImplemented());
1858  }
1859 
1860  static Number dummy;
1861  return dummy;
1862  }
1863 
1864 
1865  template <int dim, typename Number>
1866  inline
1867  const Number &
1868  symmetric_tensor_access (const TableIndices<4> &indices,
1870  {
1871  switch (dim)
1872  {
1873  case 1:
1874  return data[0][0];
1875 
1876  case 2:
1877  // each entry of the tensor can be
1878  // thought of as an entry in a
1879  // matrix that maps the rolled-out
1880  // rank-2 tensors into rolled-out
1881  // rank-2 tensors. this is the
1882  // format in which we store rank-4
1883  // tensors. determine which
1884  // position the present entry is
1885  // stored in
1886  {
1887  unsigned int base_index[2] ;
1888  if ((indices[0] == 0) && (indices[1] == 0))
1889  base_index[0] = 0;
1890  else if ((indices[0] == 1) && (indices[1] == 1))
1891  base_index[0] = 1;
1892  else
1893  base_index[0] = 2;
1894 
1895  if ((indices[2] == 0) && (indices[3] == 0))
1896  base_index[1] = 0;
1897  else if ((indices[2] == 1) && (indices[3] == 1))
1898  base_index[1] = 1;
1899  else
1900  base_index[1] = 2;
1901 
1902  return data[base_index[0]][base_index[1]];
1903  }
1904 
1905  case 3:
1906  // each entry of the tensor can be
1907  // thought of as an entry in a
1908  // matrix that maps the rolled-out
1909  // rank-2 tensors into rolled-out
1910  // rank-2 tensors. this is the
1911  // format in which we store rank-4
1912  // tensors. determine which
1913  // position the present entry is
1914  // stored in
1915  {
1916  unsigned int base_index[2] ;
1917  if ((indices[0] == 0) && (indices[1] == 0))
1918  base_index[0] = 0;
1919  else if ((indices[0] == 1) && (indices[1] == 1))
1920  base_index[0] = 1;
1921  else if ((indices[0] == 2) && (indices[1] == 2))
1922  base_index[0] = 2;
1923  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1924  ((indices[0] == 1) && (indices[1] == 0)))
1925  base_index[0] = 3;
1926  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1927  ((indices[0] == 2) && (indices[1] == 0)))
1928  base_index[0] = 4;
1929  else
1930  {
1931  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1932  ((indices[0] == 2) && (indices[1] == 1)),
1933  ExcInternalError());
1934  base_index[0] = 5;
1935  }
1936 
1937  if ((indices[2] == 0) && (indices[3] == 0))
1938  base_index[1] = 0;
1939  else if ((indices[2] == 1) && (indices[3] == 1))
1940  base_index[1] = 1;
1941  else if ((indices[2] == 2) && (indices[3] == 2))
1942  base_index[1] = 2;
1943  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1944  ((indices[2] == 1) && (indices[3] == 0)))
1945  base_index[1] = 3;
1946  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1947  ((indices[2] == 2) && (indices[3] == 0)))
1948  base_index[1] = 4;
1949  else
1950  {
1951  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1952  ((indices[2] == 2) && (indices[3] == 1)),
1953  ExcInternalError());
1954  base_index[1] = 5;
1955  }
1956 
1957  return data[base_index[0]][base_index[1]];
1958  }
1959 
1960  default:
1961  Assert (false, ExcNotImplemented());
1962  }
1963 
1964  static Number dummy;
1965  return dummy;
1966  }
1967 
1968 } // end of namespace internal
1969 
1970 
1971 
1972 template <int rank_, int dim, typename Number>
1973 inline
1974 Number &
1976 {
1977  for (unsigned int r=0; r<rank; ++r)
1978  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1979  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1980 }
1981 
1982 
1983 
1984 template <int rank_, int dim, typename Number>
1985 inline
1986 const Number &
1988 (const TableIndices<rank_> &indices) const
1989 {
1990  for (unsigned int r=0; r<rank; ++r)
1991  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1992  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1993 }
1994 
1995 
1996 
1997 namespace internal
1998 {
1999  namespace SymmetricTensor
2000  {
2001  template <int rank_>
2003  get_partially_filled_indices (const unsigned int row,
2004  const std::integral_constant<int, 2> &)
2005  {
2006  return TableIndices<rank_> (row,
2008 
2009  }
2010 
2011 
2012  template <int rank_>
2014  get_partially_filled_indices (const unsigned int row,
2015  const std::integral_constant<int, 4> &)
2016  {
2017  return TableIndices<rank_> (row,
2021 
2022  }
2023  }
2024 }
2025 
2026 
2027 template <int rank_, int dim, typename Number>
2028 internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
2029 SymmetricTensor<rank_,dim,Number>::operator [] (const unsigned int row) const
2030 {
2031  return
2032  internal::SymmetricTensorAccessors::
2033  Accessor<rank_,dim,true,rank_-1,Number> (*this,
2034  internal::SymmetricTensor::get_partially_filled_indices<rank_> (row,
2035  std::integral_constant<int, rank_>()));
2036 }
2037 
2038 
2039 
2040 template <int rank_, int dim, typename Number>
2041 internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
2043 {
2044  return
2045  internal::SymmetricTensorAccessors::
2046  Accessor<rank_,dim,false,rank_-1,Number> (*this,
2047  internal::SymmetricTensor::get_partially_filled_indices<rank_> (row,
2048  std::integral_constant<int, rank_>()));
2049 }
2050 
2051 
2052 
2053 template <int rank_, int dim, typename Number>
2054 inline
2055 const Number &
2057 {
2058  return operator()(indices);
2059 }
2060 
2061 
2062 
2063 template <int rank_, int dim, typename Number>
2064 inline
2065 Number &
2067 {
2068  return operator()(indices);
2069 }
2070 
2071 
2072 
2073 template <int rank_, int dim, typename Number>
2074 inline
2075 Number *
2077 {
2078  return std::addressof(this->access_raw_entry(0));
2079 }
2080 
2081 
2082 
2083 template <int rank_, int dim, typename Number>
2084 inline
2085 const Number *
2087 {
2088  return std::addressof(this->access_raw_entry(0));
2089 }
2090 
2091 
2092 
2093 template <int rank_, int dim, typename Number>
2094 inline
2095 Number *
2097 {
2098  return begin_raw()+n_independent_components;
2099 }
2100 
2101 
2102 
2103 template <int rank_, int dim, typename Number>
2104 inline
2105 const Number *
2107 {
2108  return begin_raw()+n_independent_components;
2109 }
2110 
2111 
2112 
2113 namespace internal
2114 {
2115  namespace SymmetricTensor
2116  {
2117  template <int dim, typename Number>
2118  unsigned int
2119  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
2120  const unsigned int index)
2121  {
2122  return index;
2123  }
2124 
2125 
2126  template <int dim, typename Number>
2128  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
2129  const unsigned int index)
2130  {
2131  return
2134  }
2135 
2136  }
2137 }
2138 
2139 
2140 
2141 template <int rank_, int dim, typename Number>
2142 inline
2143 const Number &
2144 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index) const
2145 {
2146  AssertIndexRange (index, n_independent_components);
2147  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2148 }
2149 
2150 
2151 
2152 template <int rank_, int dim, typename Number>
2153 inline
2154 Number &
2155 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index)
2156 {
2157  AssertIndexRange (index, n_independent_components);
2158  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2159 }
2160 
2161 
2162 
2163 namespace internal
2164 {
2165  template <int dim, typename Number>
2166  inline
2169  {
2170  switch (dim)
2171  {
2172  case 1:
2173  return numbers::NumberTraits<Number>::abs(data[0]);
2174 
2175  case 2:
2176  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2179 
2180  case 3:
2181  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2187 
2188  default:
2189  {
2190  typename numbers::NumberTraits<Number>::real_type return_value
2192 
2193  for (unsigned int d=0; d<dim; ++d)
2194  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
2195  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
2196  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2197 
2198  return std::sqrt(return_value);
2199  }
2200  }
2201  }
2202 
2203 
2204 
2205  template <int dim, typename Number>
2206  inline
2209  {
2210  switch (dim)
2211  {
2212  case 1:
2213  return numbers::NumberTraits<Number>::abs (data[0][0]);
2214 
2215  default:
2216  {
2217  typename numbers::NumberTraits<Number>::real_type return_value
2219 
2220  const unsigned int n_independent_components = data.dimension;
2221 
2222  for (unsigned int i=0; i<dim; ++i)
2223  for (unsigned int j=0; j<dim; ++j)
2224  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
2225  for (unsigned int i=0; i<dim; ++i)
2226  for (unsigned int j=dim; j<n_independent_components; ++j)
2227  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2228  for (unsigned int i=dim; i<n_independent_components; ++i)
2229  for (unsigned int j=0; j<dim; ++j)
2230  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2231  for (unsigned int i=dim; i<n_independent_components; ++i)
2232  for (unsigned int j=dim; j<n_independent_components; ++j)
2233  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2234 
2235  return std::sqrt(return_value);
2236  }
2237  }
2238  }
2239 
2240 } // end of namespace internal
2241 
2242 
2243 
2244 template <int rank_, int dim, typename Number>
2245 inline
2248 {
2249  return internal::compute_norm<dim,Number> (data);
2250 }
2251 
2252 
2253 
2254 namespace internal
2255 {
2256  namespace SymmetricTensor
2257  {
2258  namespace
2259  {
2260  // a function to do the unrolling from a set of indices to a
2261  // scalar index into the array in which we store the elements of
2262  // a symmetric tensor
2263  //
2264  // this function is for rank-2 tensors
2265  template <int dim>
2266  inline
2267  unsigned int
2268  component_to_unrolled_index
2269  (const TableIndices<2> &indices)
2270  {
2271  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2272  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2273 
2274  switch (dim)
2275  {
2276  case 1:
2277  {
2278  return 0;
2279  }
2280 
2281  case 2:
2282  {
2283  static const unsigned int table[2][2] = {{0, 2},
2284  {2, 1}
2285  };
2286  return table[indices[0]][indices[1]];
2287  }
2288 
2289  case 3:
2290  {
2291  static const unsigned int table[3][3] = {{0, 3, 4},
2292  {3, 1, 5},
2293  {4, 5, 2}
2294  };
2295  return table[indices[0]][indices[1]];
2296  }
2297 
2298  case 4:
2299  {
2300  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2301  {4, 1, 7, 8},
2302  {5, 7, 2, 9},
2303  {6, 8, 9, 3}
2304  };
2305  return table[indices[0]][indices[1]];
2306  }
2307 
2308  default:
2309  // for the remainder, manually figure out the numbering
2310  {
2311  if (indices[0] == indices[1])
2312  return indices[0];
2313 
2314  TableIndices<2> sorted_indices (indices);
2315  sorted_indices.sort ();
2316 
2317  for (unsigned int d=0, c=0; d<dim; ++d)
2318  for (unsigned int e=d+1; e<dim; ++e, ++c)
2319  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
2320  return dim+c;
2321 
2322  // should never get here:
2323  Assert(false, ExcInternalError());
2324  return 0;
2325  }
2326  }
2327  }
2328 
2329  // a function to do the unrolling from a set of indices to a
2330  // scalar index into the array in which we store the elements of
2331  // a symmetric tensor
2332  //
2333  // this function is for tensors of ranks not already handled
2334  // above
2335  template <int dim, int rank_>
2336  inline
2337  unsigned int
2338  component_to_unrolled_index
2339  (const TableIndices<rank_> &indices)
2340  {
2341  (void)indices;
2342  Assert (false, ExcNotImplemented());
2344  }
2345  }
2346  }
2347 }
2348 
2349 
2350 template <int rank_, int dim, typename Number>
2351 inline
2352 unsigned int
2354 (const TableIndices<rank_> &indices)
2355 {
2356  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
2357 }
2358 
2359 
2360 
2361 namespace internal
2362 {
2363  namespace SymmetricTensor
2364  {
2365  namespace
2366  {
2367  // a function to do the inverse of the unrolling from a set of
2368  // indices to a scalar index into the array in which we store
2369  // the elements of a symmetric tensor. in other words, it goes
2370  // from the scalar index into the array to a set of indices of
2371  // the tensor
2372  //
2373  // this function is for rank-2 tensors
2374  template <int dim>
2375  inline
2377  unrolled_to_component_indices
2378  (const unsigned int i,
2379  const std::integral_constant<int, 2> &)
2380  {
2383  switch (dim)
2384  {
2385  case 1:
2386  {
2387  return TableIndices<2>(0,0);
2388  }
2389 
2390  case 2:
2391  {
2392  const TableIndices<2> table[3] =
2393  {
2394  TableIndices<2> (0,0),
2395  TableIndices<2> (1,1),
2396  TableIndices<2> (0,1)
2397  };
2398  return table[i];
2399  }
2400 
2401  case 3:
2402  {
2403  const TableIndices<2> table[6] =
2404  {
2405  TableIndices<2> (0,0),
2406  TableIndices<2> (1,1),
2407  TableIndices<2> (2,2),
2408  TableIndices<2> (0,1),
2409  TableIndices<2> (0,2),
2410  TableIndices<2> (1,2)
2411  };
2412  return table[i];
2413  }
2414 
2415  default:
2416  if (i<dim)
2417  return TableIndices<2> (i,i);
2418 
2419  for (unsigned int d=0, c=0; d<dim; ++d)
2420  for (unsigned int e=d+1; e<dim; ++e, ++c)
2421  if (c==i)
2422  return TableIndices<2>(d,e);
2423 
2424  // should never get here:
2425  Assert(false, ExcInternalError());
2426  return TableIndices<2>(0, 0);
2427  }
2428  }
2429 
2430  // a function to do the inverse of the unrolling from a set of
2431  // indices to a scalar index into the array in which we store
2432  // the elements of a symmetric tensor. in other words, it goes
2433  // from the scalar index into the array to a set of indices of
2434  // the tensor
2435  //
2436  // this function is for tensors of a rank not already handled
2437  // above
2438  template <int dim, int rank_>
2439  inline
2441  unrolled_to_component_indices
2442  (const unsigned int i,
2443  const std::integral_constant<int, rank_> &)
2444  {
2445  (void)i;
2448  Assert (false, ExcNotImplemented());
2449  return TableIndices<rank_>();
2450  }
2451 
2452  }
2453  }
2454 }
2455 
2456 template <int rank_, int dim, typename Number>
2457 inline
2460 (const unsigned int i)
2461 {
2462  return
2463  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2464  std::integral_constant<int, rank_>());
2465 }
2466 
2467 
2468 
2469 template <int rank_, int dim, typename Number>
2470 template <class Archive>
2471 inline
2472 void
2473 SymmetricTensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
2474 {
2475  ar &data;
2476 }
2477 
2478 
2479 #endif // DOXYGEN
2480 
2481 /* ----------------- Non-member functions operating on tensors. ------------ */
2482 
2483 
2496 template <int rank_, int dim, typename Number, typename OtherNumber>
2497 inline
2501 {
2503  tmp += right;
2504  return tmp;
2505 }
2506 
2507 
2520 template <int rank_, int dim, typename Number, typename OtherNumber>
2521 inline
2525 {
2527  tmp -= right;
2528  return tmp;
2529 }
2530 
2531 
2539 template <int rank_, int dim, typename Number, typename OtherNumber>
2540 inline
2543  const Tensor<rank_, dim, OtherNumber> &right)
2544 {
2545  return Tensor<rank_, dim, Number>(left) + right;
2546 }
2547 
2548 
2556 template <int rank_, int dim, typename Number, typename OtherNumber>
2557 inline
2561 {
2562  return left + Tensor<rank_, dim, OtherNumber>(right);
2563 }
2564 
2565 
2573 template <int rank_, int dim, typename Number, typename OtherNumber>
2574 inline
2577  const Tensor<rank_, dim, OtherNumber> &right)
2578 {
2579  return Tensor<rank_, dim, Number>(left) - right;
2580 }
2581 
2582 
2590 template <int rank_, int dim, typename Number, typename OtherNumber>
2591 inline
2595 {
2596  return left - Tensor<rank_, dim, OtherNumber>(right);
2597 }
2598 
2599 
2600 
2614 template <int dim, typename Number>
2615 inline
2617 {
2618  switch (dim)
2619  {
2620  case 1:
2621  return t.data[0];
2622  case 2:
2623  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2624  case 3:
2625  {
2626  // in analogy to general tensors, but
2627  // there's something to be simplified for
2628  // the present case
2629  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2630  return ( tmp + tmp
2631  +t.data[0]*t.data[1]*t.data[2]
2632  -t.data[0]*t.data[5]*t.data[5]
2633  -t.data[1]*t.data[4]*t.data[4]
2634  -t.data[2]*t.data[3]*t.data[3]);
2635  }
2636  default:
2637  Assert (false, ExcNotImplemented());
2639  }
2640 }
2641 
2642 
2643 
2653 template <int dim, typename Number>
2654 inline
2656 {
2657  return determinant (t);
2658 }
2659 
2660 
2661 
2669 template <int dim, typename Number>
2671 {
2672  Number t = d.data[0];
2673  for (unsigned int i=1; i<dim; ++i)
2674  t += d.data[i];
2675  return t;
2676 }
2677 
2678 
2688 template <int dim, typename Number>
2689 inline
2691 {
2692  return trace (t);
2693 }
2694 
2695 
2708 template <typename Number>
2709 inline
2711 {
2713 }
2714 
2715 
2716 
2737 template <typename Number>
2738 inline
2740 {
2741  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2742 }
2743 
2744 
2745 
2755 template <typename Number>
2756 inline
2758 {
2759  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2760  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2761 }
2762 
2763 
2764 
2773 template <typename Number>
2774 std::array<Number,1>
2775 eigenvalues (const SymmetricTensor<2,1,Number> &T);
2776 
2777 
2778 
2800 template <typename Number>
2801 std::array<Number,2>
2802 eigenvalues (const SymmetricTensor<2,2,Number> &T);
2803 
2804 
2805 
2825 template <typename Number>
2826 std::array<Number,3>
2827 eigenvalues (const SymmetricTensor<2,3,Number> &T);
2828 
2829 
2830 
2831 namespace internal
2832 {
2833  namespace SymmetricTensor
2834  {
2871  template <int dim, typename Number>
2872  void
2873  tridiagonalize (const ::SymmetricTensor<2,dim,Number> &A,
2874  ::Tensor<2,dim,Number> &Q,
2875  std::array<Number,dim> &d,
2876  std::array<Number,dim-1> &e);
2877 
2878 
2879 
2918  template <int dim, typename Number>
2919  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2920  ql_implicit_shifts (const ::SymmetricTensor<2,dim,Number> &A);
2921 
2922 
2923 
2962  template <int dim, typename Number>
2963  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2965 
2966 
2967 
2983  template <typename Number>
2984  std::array<std::pair<Number, Tensor<1,2,Number> >,2>
2985  hybrid (const ::SymmetricTensor<2,2,Number> &A);
2986 
2987 
2988 
3021  template <typename Number>
3022  std::array<std::pair<Number, Tensor<1,3,Number> >,3>
3023  hybrid (const ::SymmetricTensor<2,3,Number> &A);
3024 
3025  namespace
3026  {
3027 
3032  template <int dim, typename Number>
3033  struct SortEigenValuesVectors
3034  {
3035  typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
3036  bool operator() (const EigValsVecs &lhs,
3037  const EigValsVecs &rhs)
3038  {
3039  return lhs.first > rhs.first;
3040  }
3041  };
3042 
3043  }
3044 
3045  } // namespace SymmetricTensor
3046 
3047 } // namespace internal
3048 
3049 
3050 
3051 // The line below is to ensure that doxygen puts the full description
3052 // of this global enumeration into the documentation
3053 // See https://stackoverflow.com/a/1717984
3081 {
3091  hybrid,
3101  ql_implicit_shifts,
3109  jacobi
3110 };
3111 
3112 
3113 
3120 template <typename Number>
3121 std::array<std::pair<Number, Tensor<1,1,Number> >,1>
3124 {
3125  return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} };
3126 }
3127 
3128 
3129 
3157 template <int dim, typename Number>
3158 std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
3161 {
3162  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
3163 
3164  switch (method)
3165  {
3167  eig_vals_vecs = internal::SymmetricTensor::hybrid(T);
3168  break;
3171  break;
3173  eig_vals_vecs = internal::SymmetricTensor::jacobi(T);
3174  break;
3175  default:
3176  AssertThrow(false, ExcNotImplemented());
3177  }
3178 
3179  // Sort in descending order before output.
3180  std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
3181  internal::SymmetricTensor::SortEigenValuesVectors<dim,Number>());
3182  return eig_vals_vecs;
3183 }
3184 
3185 
3186 
3196 template <int rank_, int dim, typename Number>
3197 inline
3200 {
3201  return t;
3202 }
3203 
3204 
3205 
3215 template <int dim, typename Number>
3216 inline
3219 {
3221 
3222  // subtract scaled trace from the diagonal
3223  const Number tr = trace(t) / dim;
3224  for (unsigned int i=0; i<dim; ++i)
3225  tmp.data[i] -= tr;
3226 
3227  return tmp;
3228 }
3229 
3230 
3231 
3239 template <int dim, typename Number>
3240 inline
3242 unit_symmetric_tensor ()
3243 {
3244  // create a default constructed matrix filled with
3245  // zeros, then set the diagonal elements to one
3247  switch (dim)
3248  {
3249  case 1:
3250  tmp.data[0] = 1;
3251  break;
3252  case 2:
3253  tmp.data[0] = tmp.data[1] = 1;
3254  break;
3255  case 3:
3256  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3257  break;
3258  default:
3259  for (unsigned int d=0; d<dim; ++d)
3260  tmp.data[d] = 1;
3261  }
3262  return tmp;
3263 }
3264 
3265 
3266 
3275 template <int dim>
3276 inline
3278 unit_symmetric_tensor ()
3279 {
3280  return unit_symmetric_tensor<dim,double>();
3281 }
3282 
3283 
3284 
3299 template <int dim, typename Number>
3300 inline
3302 deviator_tensor ()
3303 {
3305 
3306  // fill the elements treating the diagonal
3307  for (unsigned int i=0; i<dim; ++i)
3308  for (unsigned int j=0; j<dim; ++j)
3309  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
3310 
3311  // then fill the ones that copy over the
3312  // non-diagonal elements. note that during
3313  // the double-contraction, we handle the
3314  // off-diagonal elements twice, so simply
3315  // copying requires a weight of 1/2
3316  for (unsigned int i=dim;
3317  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3318  ++i)
3319  tmp.data[i][i] = 0.5;
3320 
3321  return tmp;
3322 }
3323 
3324 
3325 
3340 template <int dim>
3341 inline
3343 deviator_tensor ()
3344 {
3345  return deviator_tensor<dim,double>();
3346 }
3347 
3348 
3349 
3372 template <int dim, typename Number>
3373 inline
3375 identity_tensor ()
3376 {
3378 
3379  // fill the elements treating the diagonal
3380  for (unsigned int i=0; i<dim; ++i)
3381  tmp.data[i][i] = 1;
3382 
3383  // then fill the ones that copy over the
3384  // non-diagonal elements. note that during
3385  // the double-contraction, we handle the
3386  // off-diagonal elements twice, so simply
3387  // copying requires a weight of 1/2
3388  for (unsigned int i=dim;
3389  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3390  ++i)
3391  tmp.data[i][i] = 0.5;
3392 
3393  return tmp;
3394 }
3395 
3396 
3397 
3419 template <int dim>
3420 inline
3422 identity_tensor ()
3423 {
3424  return identity_tensor<dim,double>();
3425 }
3426 
3427 
3428 
3439 template <int dim, typename Number>
3440 inline
3443 {
3445 }
3446 
3447 
3448 
3460 template <int dim, typename Number>
3461 inline
3464 {
3466 }
3467 
3468 
3469 
3484 template <int dim, typename Number>
3485 inline
3489 {
3491 
3492  // fill only the elements really needed
3493  for (unsigned int i=0; i<dim; ++i)
3494  for (unsigned int j=i; j<dim; ++j)
3495  for (unsigned int k=0; k<dim; ++k)
3496  for (unsigned int l=k; l<dim; ++l)
3497  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3498 
3499  return tmp;
3500 }
3501 
3502 
3503 
3512 template <int dim,typename Number>
3513 inline
3516 {
3517  Number array[(dim*dim+dim)/2];
3518  for (unsigned int d=0; d<dim; ++d)
3519  array[d] = t[d][d];
3520  for (unsigned int d=0, c=0; d<dim; ++d)
3521  for (unsigned int e=d+1; e<dim; ++e, ++c)
3522  array[dim+c] = (t[d][e]+t[e][d])*0.5;
3523  return SymmetricTensor<2,dim,Number>(array);
3524 }
3525 
3526 
3527 
3535 template <int rank_, int dim, typename Number>
3536 inline
3539  const Number &factor)
3540 {
3542  tt *= factor;
3543  return tt;
3544 }
3545 
3546 
3547 
3555 template <int rank_, int dim, typename Number>
3556 inline
3558 operator * (const Number &factor,
3560 {
3561  // simply forward to the other operator
3562  return t*factor;
3563 }
3564 
3565 
3566 
3592 template <int rank_, int dim, typename Number, typename OtherNumber>
3593 inline
3596  const OtherNumber &factor)
3597 {
3598  // form the product. we have to convert the two factors into the final
3599  // type via explicit casts because, for awkward reasons, the C++
3600  // standard committee saw it fit to not define an
3601  // operator*(float,std::complex<double>)
3602  // (as well as with switched arguments and double<->float).
3603  typedef typename ProductType<Number,OtherNumber>::type product_type;
3605  // we used to shorten the following by 'tt *= product_type(factor);'
3606  // which requires that a converting constructor
3607  // 'product_type::product_type(const OtherNumber) is defined.
3608  // however, a user-defined constructor is not allowed for aggregates,
3609  // e.g. VectorizedArray. therefore, we work around this issue using a
3610  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3611  // which we assume to be defined.
3612  product_type new_factor;
3613  new_factor = factor;
3614  tt *= new_factor;
3615  return tt;
3616 }
3617 
3618 
3619 
3628 template <int rank_, int dim, typename Number, typename OtherNumber>
3629 inline
3631 operator * (const Number &factor,
3633 {
3634  // simply forward to the other operator with switched arguments
3635  return (t*factor);
3636 }
3637 
3638 
3639 
3645 template <int rank_, int dim, typename Number, typename OtherNumber>
3646 inline
3649  const OtherNumber &factor)
3650 {
3652  tt /= factor;
3653  return tt;
3654 }
3655 
3656 
3657 
3664 template <int rank_, int dim>
3665 inline
3667 operator * (const SymmetricTensor<rank_,dim> &t,
3668  const double factor)
3669 {
3671  tt *= factor;
3672  return tt;
3673 }
3674 
3675 
3676 
3683 template <int rank_, int dim>
3684 inline
3686 operator * (const double factor,
3687  const SymmetricTensor<rank_,dim> &t)
3688 {
3690  tt *= factor;
3691  return tt;
3692 }
3693 
3694 
3695 
3701 template <int rank_, int dim>
3702 inline
3704 operator / (const SymmetricTensor<rank_,dim> &t,
3705  const double factor)
3706 {
3708  tt /= factor;
3709  return tt;
3710 }
3711 
3721 template <int dim, typename Number, typename OtherNumber>
3722 inline
3723 typename ProductType<Number, OtherNumber>::type
3726 {
3727  return (t1*t2);
3728 }
3729 
3730 
3740 template <int dim, typename Number, typename OtherNumber>
3741 inline
3742 typename ProductType<Number, OtherNumber>::type
3744  const Tensor<2,dim,OtherNumber> &t2)
3745 {
3746  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<typename ProductType<Number, OtherNumber>::type>::value(0.0);
3747  for (unsigned int i=0; i<dim; ++i)
3748  for (unsigned int j=0; j<dim; ++j)
3749  s += t1[i][j] * t2[i][j];
3750  return s;
3751 }
3752 
3753 
3763 template <int dim, typename Number, typename OtherNumber>
3764 inline
3765 typename ProductType<Number, OtherNumber>::type
3768 {
3769  return scalar_product(t2, t1);
3770 }
3771 
3772 
3788 template <typename Number, typename OtherNumber>
3789 inline
3790 void
3791 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3792  const SymmetricTensor<4,1,Number> &t,
3794 {
3795  tmp[0][0] = t[0][0][0][0] * s[0][0];
3796 }
3797 
3798 
3799 
3815 template <typename Number, typename OtherNumber>
3816 inline
3817 void
3818 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3819  const SymmetricTensor<2,1,Number> &s,
3821 {
3822  tmp[0][0] = t[0][0][0][0] * s[0][0];
3823 }
3824 
3825 
3826 
3841 template <typename Number, typename OtherNumber>
3842 inline
3843 void
3844 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3845  const SymmetricTensor<4,2,Number> &t,
3847 {
3848  const unsigned int dim = 2;
3849 
3850  for (unsigned int i=0; i<dim; ++i)
3851  for (unsigned int j=i; j<dim; ++j)
3852  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3853  t[i][j][1][1] * s[1][1] +
3854  2 * t[i][j][0][1] * s[0][1];
3855 }
3856 
3857 
3858 
3874 template <typename Number, typename OtherNumber>
3875 inline
3876 void
3877 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3878  const SymmetricTensor<2,2,Number> &s,
3880 {
3881  const unsigned int dim = 2;
3882 
3883  for (unsigned int i=0; i<dim; ++i)
3884  for (unsigned int j=i; j<dim; ++j)
3885  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3886  s[1][1] * t[1][1][i][j] +
3887  2 * s[0][1] * t[0][1][i][j];
3888 }
3889 
3890 
3891 
3907 template <typename Number, typename OtherNumber>
3908 inline
3909 void
3910 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3911  const SymmetricTensor<4,3,Number> &t,
3913 {
3914  const unsigned int dim = 3;
3915 
3916  for (unsigned int i=0; i<dim; ++i)
3917  for (unsigned int j=i; j<dim; ++j)
3918  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3919  t[i][j][1][1] * s[1][1] +
3920  t[i][j][2][2] * s[2][2] +
3921  2 * t[i][j][0][1] * s[0][1] +
3922  2 * t[i][j][0][2] * s[0][2] +
3923  2 * t[i][j][1][2] * s[1][2];
3924 }
3925 
3926 
3927 
3943 template <typename Number, typename OtherNumber>
3944 inline
3945 void
3946 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3947  const SymmetricTensor<2,3,Number> &s,
3949 {
3950  const unsigned int dim = 3;
3951 
3952  for (unsigned int i=0; i<dim; ++i)
3953  for (unsigned int j=i; j<dim; ++j)
3954  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3955  s[1][1] * t[1][1][i][j] +
3956  s[2][2] * t[2][2][i][j] +
3957  2 * s[0][1] * t[0][1][i][j] +
3958  2 * s[0][2] * t[0][2][i][j] +
3959  2 * s[1][2] * t[1][2][i][j];
3960 }
3961 
3962 
3963 
3971 template <int dim, typename Number, typename OtherNumber>
3973 operator * (const SymmetricTensor<2,dim,Number> &src1,
3974  const Tensor<1,dim,OtherNumber> &src2)
3975 {
3977  for (unsigned int i=0; i<dim; ++i)
3978  for (unsigned int j=0; j<dim; ++j)
3979  dest[i] += src1[i][j] * src2[j];
3980  return dest;
3981 }
3982 
3983 
3991 template <int dim, typename Number, typename OtherNumber>
3993 operator * (const Tensor<1,dim,Number> &src1,
3995 {
3996  // this is easy for symmetric tensors:
3997  return src2 * src1;
3998 }
3999 
4000 
4001 
4022 template <int rank_1, int rank_2, int dim,
4023  typename Number, typename OtherNumber>
4024 inline DEAL_II_ALWAYS_INLINE
4025 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
4026 operator * (const Tensor<rank_1, dim, Number> &src1,
4028 {
4029  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
4030  const Tensor<rank_2, dim, OtherNumber> src2 (src2s);
4031  return src1*src2;
4032 }
4033 
4034 
4035 
4056 template <int rank_1, int rank_2, int dim,
4057  typename Number, typename OtherNumber>
4058 inline DEAL_II_ALWAYS_INLINE
4059 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
4060 operator * (const SymmetricTensor<rank_1, dim, Number> &src1s,
4062 {
4063  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
4064  const Tensor<rank_2, dim, OtherNumber> src1 (src1s);
4065  return src1*src2;
4066 }
4067 
4068 
4069 
4079 template <int dim, typename Number>
4080 inline
4081 std::ostream &operator << (std::ostream &out,
4083 {
4084  //make out lives a bit simpler by outputing
4085  //the tensor through the operator for the
4086  //general Tensor class
4088 
4089  for (unsigned int i=0; i<dim; ++i)
4090  for (unsigned int j=0; j<dim; ++j)
4091  tt[i][j] = t[i][j];
4092 
4093  return out << tt;
4094 }
4095 
4096 
4097 
4107 template <int dim, typename Number>
4108 inline
4109 std::ostream &operator << (std::ostream &out,
4111 {
4112  //make out lives a bit simpler by outputing
4113  //the tensor through the operator for the
4114  //general Tensor class
4116 
4117  for (unsigned int i=0; i<dim; ++i)
4118  for (unsigned int j=0; j<dim; ++j)
4119  for (unsigned int k=0; k<dim; ++k)
4120  for (unsigned int l=0; l<dim; ++l)
4121  tt[i][j][k][l] = t[i][j][k][l];
4122 
4123  return out << tt;
4124 }
4125 
4126 
4127 DEAL_II_NAMESPACE_CLOSE
4128 
4129 #endif
friend SymmetricTensor< 4, dim2, Number2 > identity_tensor()
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A,::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim-1 > &e)
static const unsigned int invalid_unsigned_int
Definition: types.h:173
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static const unsigned int n_independent_components
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
std::array< std::pair< Number, Tensor< 1, 2, Number > >, 2 > hybrid(const ::SymmetricTensor< 2, 2, Number > &A)
void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 2, Number > &s, const SymmetricTensor< 4, 2, OtherNumber > &t)
void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1237
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, OtherNumber > &t2)
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
bool operator==(const SymmetricTensor &) const
SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
static std::size_t memory_consumption()
numbers::NumberTraits< Number >::real_type norm() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:410
static real_type abs(const number &x)
Definition: numbers.h:354
SymmetricTensorEigenvectorMethod
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
internal::SymmetricTensorAccessors::StorageType< rank_, dim, Number > base_tensor_descriptor
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor & operator/=(const OtherNumber &factor)
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, OtherNumber > &s)
static const unsigned int dimension
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static::ExceptionBase & ExcMessage(std::string arg1)
const Number & access_raw_entry(const unsigned int unrolled_index) const
TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_-1, Number > operator[](const unsigned int row) const
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
base_tensor_descriptor::base_tensor_type base_tensor_type
friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
static const unsigned int rank
T sum(const T &t, const MPI_Comm &mpi_communicator)
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
Number * begin_raw()
void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 3, Number > &t, const SymmetricTensor< 2, 3, OtherNumber > &s)
#define Assert(cond, exc)
Definition: exceptions.h:349
base_tensor_type data
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Number trace(const SymmetricTensor< 2, dim, Number > &d)
void serialize(Archive &ar, const unsigned int version)
Number * end_raw()
friend SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Number & operator()(const TableIndices< rank_ > &indices)
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
Definition: mpi.h:53
SymmetricTensor & operator*=(const OtherNumber &factor)
SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
StreamType & operator<<(StreamType &s, const UpdateFlags u)
bool operator!=(const SymmetricTensor &) const
Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Tensor< 1, n_independent_components, Number > base_tensor_type
static::ExceptionBase & ExcNotImplemented()
SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > ql_implicit_shifts(const ::SymmetricTensor< 2, dim, Number > &A)
Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
friend SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, OtherNumber > &t)
T max(const T &t, const MPI_Comm &mpi_communicator)
void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 1, Number > &s, const SymmetricTensor< 4, 1, OtherNumber > &t)
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
std::array< std::pair< Number, Tensor< 1, 1, Number > >, 1 > eigenvectors(const SymmetricTensor< 2, 1, Number > &T, const SymmetricTensorEigenvectorMethod=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor operator-() const
static::ExceptionBase & ExcInternalError()