Reference documentation for deal.II version Git 8a320d4 2018-02-19 08:48:02 +0100
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> SymmetricTensor<2,dim,Number>
35 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
37 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
39 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
41 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
43 template <int dim2, typename Number> Number
45 
46 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
48 template <int dim, typename Number> Number
50 
51 
52 
53 namespace internal
54 {
59  namespace SymmetricTensor
60  {
65  template <int rank, int dim, typename Number>
66  struct Inverse;
67  }
68 
73  namespace SymmetricTensorAccessors
74  {
81  inline
82  TableIndices<2> merge (const TableIndices<2> &previous_indices,
83  const unsigned int new_index,
84  const unsigned int position)
85  {
86  Assert (position < 2, ExcIndexRange (position, 0, 2));
87 
88  if (position == 0)
90  else
91  return TableIndices<2>(previous_indices[0], new_index);
92  }
93 
94 
95 
102  inline
103  TableIndices<4> merge (const TableIndices<4> &previous_indices,
104  const unsigned int new_index,
105  const unsigned int position)
106  {
107  Assert (position < 4, ExcIndexRange (position, 0, 4));
108 
109  switch (position)
110  {
111  case 0:
112  return TableIndices<4>(new_index,
116  case 1:
117  return TableIndices<4>(previous_indices[0],
118  new_index,
121  case 2:
122  return TableIndices<4>(previous_indices[0],
123  previous_indices[1],
124  new_index,
126  case 3:
127  return TableIndices<4>(previous_indices[0],
128  previous_indices[1],
129  previous_indices[2],
130  new_index);
131  }
132  Assert (false, ExcInternalError());
133  return TableIndices<4>();
134  }
135 
136 
145  template <int rank1, int rank2, int dim, typename Number, typename OtherNumber = Number>
147  {
148  typedef typename ProductType<Number,OtherNumber>::type value_type;
149  typedef ::SymmetricTensor<rank1+rank2-4,dim,value_type> type;
150  };
151 
152 
161  template <int dim, typename Number, typename OtherNumber>
162  struct double_contraction_result<2,2,dim,Number,OtherNumber>
163  {
164  typedef typename ProductType<Number,OtherNumber>::type type;
165  };
166 
167 
168 
181  template <int rank, int dim, typename Number>
182  struct StorageType;
183 
187  template <int dim, typename Number>
188  struct StorageType<2,dim,Number>
189  {
194  static const unsigned int
195  n_independent_components = (dim*dim + dim)/2;
196 
201  };
202 
203 
204 
208  template <int dim, typename Number>
209  struct StorageType<4,dim,Number>
210  {
216  static const unsigned int
217  n_rank2_components = (dim*dim + dim)/2;
218 
222  static const unsigned int
223  n_independent_components = (n_rank2_components *
225 
233  };
234 
235 
236 
241  template <int rank, int dim, bool constness, typename Number>
243 
250  template <int rank, int dim, typename Number>
251  struct AccessorTypes<rank,dim,true,Number>
252  {
253  typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
254 
255  typedef Number reference;
256  };
257 
264  template <int rank, int dim, typename Number>
265  struct AccessorTypes<rank,dim,false,Number>
266  {
267  typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
268 
269  typedef Number &reference;
270  };
271 
272 
307  template <int rank, int dim, bool constness, int P, typename Number>
308  class Accessor
309  {
310  public:
314  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
315  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
316 
317  private:
336  Accessor (tensor_type &tensor,
337  const TableIndices<rank> &previous_indices);
338 
342  Accessor (const Accessor &) = default;
343 
344  public:
345 
349  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
350 
354  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;
355 
356  private:
360  tensor_type &tensor;
361  const TableIndices<rank> previous_indices;
362 
363  // declare some other classes
364  // as friends. make sure to
365  // work around bugs in some
366  // compilers
367  template <int,int,typename> friend class ::SymmetricTensor;
368  template <int,int,bool,int,typename>
369  friend class Accessor;
370 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
371  friend class ::SymmetricTensor<rank,dim,Number>;
372  friend class Accessor<rank,dim,constness,P+1,Number>;
373 # endif
374  };
375 
376 
377 
387  template <int rank, int dim, bool constness, typename Number>
388  class Accessor<rank,dim,constness,1,Number>
389  {
390  public:
394  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
395  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
396 
397  private:
419  Accessor (tensor_type &tensor,
420  const TableIndices<rank> &previous_indices);
421 
425  Accessor () = delete;
426 
430  Accessor (const Accessor &) = default;
431 
432  public:
433 
437  reference operator [] (const unsigned int);
438 
442  reference operator [] (const unsigned int) const;
443 
444  private:
448  tensor_type &tensor;
449  const TableIndices<rank> previous_indices;
450 
451  // declare some other classes
452  // as friends. make sure to
453  // work around bugs in some
454  // compilers
455  template <int,int,typename> friend class ::SymmetricTensor;
456  template <int,int,bool,int,typename>
457  friend class SymmetricTensorAccessors::Accessor;
458 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
459  friend class ::SymmetricTensor<rank,dim,Number>;
460  friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
461 # endif
462  };
463  }
464 }
465 
466 
467 
531 template <int rank_, int dim, typename Number>
532 class SymmetricTensor
533 {
534 public:
535  static_assert(rank_%2==0, "A SymmetricTensor must have even rank!");
536 
545  static const unsigned int dimension = dim;
546 
550  static const unsigned int rank = rank_;
551 
557  static const unsigned int n_independent_components
560 
564  SymmetricTensor ();
565 
576  template <typename OtherNumber>
577  explicit
579 
595  SymmetricTensor (const Number (&array) [n_independent_components]);
596 
602  template <typename OtherNumber>
603  explicit
605 
609  Number *
610  begin_raw();
611 
615  const Number *
616  begin_raw() const;
617 
621  Number *
622  end_raw();
623 
628  const Number *
629  end_raw() const;
630 
636  template <typename OtherNumber>
638 
645  SymmetricTensor &operator = (const Number &d);
646 
651  operator Tensor<rank_,dim,Number> () const;
652 
656  bool operator == (const SymmetricTensor &) const;
657 
661  bool operator != (const SymmetricTensor &) const;
662 
666  template <typename OtherNumber>
668 
672  template <typename OtherNumber>
674 
679  template <typename OtherNumber>
680  SymmetricTensor &operator *= (const OtherNumber &factor);
681 
685  template <typename OtherNumber>
686  SymmetricTensor &operator /= (const OtherNumber &factor);
687 
692 
717  template <typename OtherNumber>
720 
725  template <typename OtherNumber>
728 
732  Number &operator() (const TableIndices<rank_> &indices);
733 
737  const Number &operator() (const TableIndices<rank_> &indices) const;
738 
743  internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
744  operator [] (const unsigned int row) const;
745 
750  internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
751  operator [] (const unsigned int row);
752 
758  const Number &
759  operator [] (const TableIndices<rank_> &indices) const;
760 
766  Number &
767  operator [] (const TableIndices<rank_> &indices);
768 
774  const Number &
775  access_raw_entry (const unsigned int unrolled_index) const;
776 
782  Number &
783  access_raw_entry (const unsigned int unrolled_index);
784 
795  norm () const;
796 
804  static
805  unsigned int
807 
813  static
815  unrolled_to_component_indices (const unsigned int i);
816 
829  void clear ();
830 
835  static std::size_t memory_consumption ();
836 
841  template <class Archive>
842  void serialize(Archive &ar, const unsigned int version);
843 
844 private:
848  typedef
851 
855  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
856 
860  base_tensor_type data;
861 
865  template <int, int, typename> friend class SymmetricTensor;
866 
870  template <int dim2, typename Number2>
871  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
872 
873  template <int dim2, typename Number2>
874  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
875 
876  template <int dim2, typename Number2>
878  deviator (const SymmetricTensor<2,dim2,Number2> &t);
879 
880  template <int dim2, typename Number2>
882 
883  template <int dim2, typename Number2>
885 
886  template <int dim2, typename Number2>
888 
889 
893  friend struct internal::SymmetricTensor::Inverse<2,dim,Number>;
894 
895  friend struct internal::SymmetricTensor::Inverse<4,dim,Number>;
896 };
897 
898 
899 
900 // ------------------------- inline functions ------------------------
901 
902 #ifndef DOXYGEN
903 
904 namespace internal
905 {
906  namespace SymmetricTensorAccessors
907  {
908  template <int rank_, int dim, bool constness, int P, typename Number>
909  Accessor<rank_,dim,constness,P,Number>::
910  Accessor (tensor_type &tensor,
911  const TableIndices<rank_> &previous_indices)
912  :
913  tensor (tensor),
914  previous_indices (previous_indices)
915  {}
916 
917 
918 
919  template <int rank_, int dim, bool constness, int P, typename Number>
920  Accessor<rank_,dim,constness,P-1,Number>
921  Accessor<rank_,dim,constness,P,Number>::operator[] (const unsigned int i)
922  {
923  return Accessor<rank_,dim,constness,P-1,Number> (tensor,
924  merge (previous_indices, i, rank_-P));
925  }
926 
927 
928 
929  template <int rank_, int dim, bool constness, int P, typename Number>
930  Accessor<rank_,dim,constness,P-1,Number>
931  Accessor<rank_,dim,constness,P,Number>::operator[] (const unsigned int i) const
932  {
933  return Accessor<rank_,dim,constness,P-1,Number> (tensor,
934  merge (previous_indices, i, rank_-P));
935  }
936 
937 
938 
939  template <int rank_, int dim, bool constness, typename Number>
940  Accessor<rank_,dim,constness,1,Number>::
941  Accessor (tensor_type &tensor,
942  const TableIndices<rank_> &previous_indices)
943  :
944  tensor (tensor),
945  previous_indices (previous_indices)
946  {}
947 
948 
949 
950  template <int rank_, int dim, bool constness, typename Number>
951  typename Accessor<rank_,dim,constness,1,Number>::reference
952  Accessor<rank_,dim,constness,1,Number>::operator[] (const unsigned int i)
953  {
954  return tensor(merge (previous_indices, i, rank_-1));
955  }
956 
957 
958  template <int rank_, int dim, bool constness, typename Number>
959  typename Accessor<rank_,dim,constness,1,Number>::reference
960  Accessor<rank_,dim,constness,1,Number>::operator[] (const unsigned int i) const
961  {
962  return tensor(merge (previous_indices, i, rank_-1));
963  }
964  }
965 }
966 
967 
968 
969 template <int rank_, int dim, typename Number>
970 inline
972 {
973  // Some auto-differentiable numbers need explicit
974  // zero initialization.
975  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
977 }
978 
979 
980 template <int rank_, int dim, typename Number>
981 template <typename OtherNumber>
982 inline
984 {
985  Assert (rank == 2, ExcNotImplemented());
986  switch (dim)
987  {
988  case 2:
989  Assert (t[0][1] == t[1][0], ExcInternalError());
990 
991  data[0] = t[0][0];
992  data[1] = t[1][1];
993  data[2] = t[0][1];
994 
995  break;
996  case 3:
997  Assert (t[0][1] == t[1][0], ExcInternalError());
998  Assert (t[0][2] == t[2][0], ExcInternalError());
999  Assert (t[1][2] == t[2][1], ExcInternalError());
1000 
1001  data[0] = t[0][0];
1002  data[1] = t[1][1];
1003  data[2] = t[2][2];
1004  data[3] = t[0][1];
1005  data[4] = t[0][2];
1006  data[5] = t[1][2];
1007 
1008  break;
1009  default:
1010  for (unsigned int d=0; d<dim; ++d)
1011  for (unsigned int e=0; e<d; ++e)
1012  Assert(t[d][e] == t[e][d], ExcInternalError());
1013 
1014  for (unsigned int d=0; d<dim; ++d)
1015  data[d] = t[d][d];
1016 
1017  for (unsigned int d=0, c=0; d<dim; ++d)
1018  for (unsigned int e=d+1; e<dim; ++e, ++c)
1019  data[dim+c] = t[d][e];
1020  }
1021 }
1022 
1023 
1024 
1025 template <int rank_, int dim, typename Number>
1026 template <typename OtherNumber>
1027 inline
1030 {
1031  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1033 }
1034 
1035 
1036 
1037 
1038 template <int rank_, int dim, typename Number>
1039 inline
1040 SymmetricTensor<rank_,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
1041  :
1042  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1043 {
1044  // ensure that the reinterpret_cast above actually works
1045  Assert (sizeof(typename base_tensor_type::array_type)
1046  == sizeof(array),
1047  ExcInternalError());
1048 }
1049 
1050 
1051 
1052 template <int rank_, int dim, typename Number>
1053 template <typename OtherNumber>
1054 inline
1057 {
1058  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1059  data[i] = t.data[i];
1060  return *this;
1061 }
1062 
1063 
1064 
1065 template <int rank_, int dim, typename Number>
1066 inline
1069 {
1070  Assert (d==internal::NumberType<Number>::value(0.0), ExcMessage ("Only assignment with zero is allowed"));
1071  (void) d;
1072 
1074 
1075  return *this;
1076 }
1077 
1078 
1079 namespace internal
1080 {
1081  namespace SymmetricTensor
1082  {
1083  template <int dim, typename Number>
1085  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1086  {
1088 
1089  // diagonal entries are stored first
1090  for (unsigned int d=0; d<dim; ++d)
1091  t[d][d] = s.access_raw_entry(d);
1092 
1093  // off-diagonal entries come next, row by row
1094  for (unsigned int d=0, c=0; d<dim; ++d)
1095  for (unsigned int e=d+1; e<dim; ++e, ++c)
1096  {
1097  t[d][e] = s.access_raw_entry(dim+c);
1098  t[e][d] = s.access_raw_entry(dim+c);
1099  }
1100  return t;
1101  }
1102 
1103 
1104  template <int dim, typename Number>
1106  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1107  {
1108  // utilize the symmetry properties of SymmetricTensor<4,dim>
1109  // discussed in the class documentation to avoid accessing all
1110  // independent elements of the input tensor more than once
1112 
1113  for (unsigned int i=0; i<dim; ++i)
1114  for (unsigned int j=i; j<dim; ++j)
1115  for (unsigned int k=0; k<dim; ++k)
1116  for (unsigned int l=k; l<dim; ++l)
1117  t[TableIndices<4>(i,j,k,l)]
1118  = t[TableIndices<4>(i,j,l,k)]
1119  = t[TableIndices<4>(j,i,k,l)]
1120  = t[TableIndices<4>(j,i,l,k)]
1121  = st[TableIndices<4>(i,j,k,l)];
1122 
1123  return t;
1124  }
1125 
1126 
1127  template <typename Number>
1128  struct Inverse<2,1,Number>
1129  {
1130  static inline ::SymmetricTensor<2,1,Number>
1131  value (const ::SymmetricTensor<2,1,Number> &t)
1132  {
1134 
1135  tmp[0][0] = 1.0/t[0][0];
1136 
1137  return tmp;
1138  }
1139  };
1140 
1141 
1142  template <typename Number>
1143  struct Inverse<2,2,Number>
1144  {
1145  static inline ::SymmetricTensor<2,2,Number>
1146  value (const ::SymmetricTensor<2,2,Number> &t)
1147  {
1149 
1150  // Sympy result: ([
1151  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1152  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1153  const TableIndices<2> idx_00 (0,0);
1154  const TableIndices<2> idx_01 (0,1);
1155  const TableIndices<2> idx_11 (1,1);
1156  const Number inv_det_t
1157  = 1.0/(t[idx_00]*t[idx_11]
1158  - t[idx_01]*t[idx_01]);
1159  tmp[idx_00] = t[idx_11];
1160  tmp[idx_01] = -t[idx_01];
1161  tmp[idx_11] = t[idx_00];
1162  tmp *= inv_det_t;
1163 
1164  return tmp;
1165  }
1166  };
1167 
1168 
1169  template <typename Number>
1170  struct Inverse<2,3,Number>
1171  {
1172  static ::SymmetricTensor<2,3,Number>
1173  value (const ::SymmetricTensor<2,3,Number> &t)
1174  {
1176 
1177  // Sympy result: ([
1178  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1179  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1180  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
1181  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1182  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1183  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
1184  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1185  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
1186  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
1187  const TableIndices<2> idx_00 (0,0);
1188  const TableIndices<2> idx_01 (0,1);
1189  const TableIndices<2> idx_02 (0,2);
1190  const TableIndices<2> idx_11 (1,1);
1191  const TableIndices<2> idx_12 (1,2);
1192  const TableIndices<2> idx_22 (2,2);
1193  const Number inv_det_t
1194  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
1195  - t[idx_00]*t[idx_12]*t[idx_12]
1196  - t[idx_01]*t[idx_01]*t[idx_22]
1197  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
1198  - t[idx_02]*t[idx_02]*t[idx_11]);
1199  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
1200  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
1201  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
1202  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
1203  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
1204  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
1205  tmp *= inv_det_t;
1206 
1207  return tmp;
1208  }
1209  };
1210 
1211 
1212  template <typename Number>
1213  struct Inverse<4,1,Number>
1214  {
1215  static inline ::SymmetricTensor<4,1,Number>
1216  value (const ::SymmetricTensor<4,1,Number> &t)
1217  {
1219  tmp.data[0][0] = 1.0/t.data[0][0];
1220  return tmp;
1221  }
1222  };
1223 
1224 
1225  template <typename Number>
1226  struct Inverse<4,2,Number>
1227  {
1228  static inline ::SymmetricTensor<4,2,Number>
1229  value (const ::SymmetricTensor<4,2,Number> &t)
1230  {
1232 
1233  // Inverting this tensor is a little more complicated than necessary,
1234  // since we store the data of 't' as a 3x3 matrix t.data, but the
1235  // product between a rank-4 and a rank-2 tensor is really not the
1236  // product between this matrix and the 3-vector of a rhs, but rather
1237  //
1238  // B.vec = t.data * mult * A.vec
1239  //
1240  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1241  // capture the fact that we need to add up both the c_ij12*a_12 and the
1242  // c_ij21*a_21 terms.
1243  //
1244  // In addition, in this scheme, the identity tensor has the matrix
1245  // representation mult^-1.
1246  //
1247  // The inverse of 't' therefore has the matrix representation
1248  //
1249  // inv.data = mult^-1 * t.data^-1 * mult^-1
1250  //
1251  // in order to compute it, let's first compute the inverse of t.data and
1252  // put it into tmp.data; at the end of the function we then scale the
1253  // last row and column of the inverse by 1/2, corresponding to the left
1254  // and right multiplication with mult^-1.
1255  const Number t4 = t.data[0][0]*t.data[1][1],
1256  t6 = t.data[0][0]*t.data[1][2],
1257  t8 = t.data[0][1]*t.data[1][0],
1258  t00 = t.data[0][2]*t.data[1][0],
1259  t01 = t.data[0][1]*t.data[2][0],
1260  t04 = t.data[0][2]*t.data[2][0],
1261  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
1262  t8*t.data[2][2]+t00*t.data[2][1]+
1263  t01*t.data[1][2]-t04*t.data[1][1]);
1264  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
1265  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
1266  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
1267  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
1268  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
1269  tmp.data[1][2] = -(t6-t00)*t07;
1270  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
1271  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
1272  tmp.data[2][2] = (t4-t8)*t07;
1273 
1274  // scale last row and column as mentioned
1275  // above
1276  tmp.data[2][0] /= 2;
1277  tmp.data[2][1] /= 2;
1278  tmp.data[0][2] /= 2;
1279  tmp.data[1][2] /= 2;
1280  tmp.data[2][2] /= 4;
1281 
1282  return tmp;
1283  }
1284  };
1285 
1286 
1287  template <typename Number>
1288  struct Inverse<4,3,Number>
1289  {
1290  static ::SymmetricTensor<4,3,Number>
1291  value (const ::SymmetricTensor<4,3,Number> &t)
1292  {
1294 
1295  // This function follows the exact same scheme as the 2d case, except
1296  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1297  // Instead, we use the Gauss-Jordan algorithm implemented for
1298  // FullMatrix. For historical reasons the following code is copied from
1299  // there, with the tangential benefit that we do not need to copy the
1300  // tensor entries to and from the FullMatrix.
1301  const unsigned int N = 6;
1302 
1303  // First get an estimate of the size of the elements of this matrix,
1304  // for later checks whether the pivot element is large enough, or
1305  // whether we have to fear that the matrix is not regular.
1306  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1307  for (unsigned int i=0; i<N; ++i)
1308  diagonal_sum += std::fabs(tmp.data[i][i]);
1309  const Number typical_diagonal_element = diagonal_sum/N;
1310  (void)typical_diagonal_element;
1311 
1312  unsigned int p[N];
1313  for (unsigned int i=0; i<N; ++i)
1314  p[i] = i;
1315 
1316  for (unsigned int j=0; j<N; ++j)
1317  {
1318  // Pivot search: search that part of the line on and right of the
1319  // diagonal for the largest element.
1320  Number max = std::fabs(tmp.data[j][j]);
1321  unsigned int r = j;
1322  for (unsigned int i=j+1; i<N; ++i)
1323  if (std::fabs(tmp.data[i][j]) > max)
1324  {
1325  max = std::fabs(tmp.data[i][j]);
1326  r = i;
1327  }
1328 
1329  // Check whether the pivot is too small
1330  Assert(max > 1.e-16*typical_diagonal_element,
1331  ExcMessage("This tensor seems to be noninvertible"));
1332 
1333  // Row interchange
1334  if (r>j)
1335  {
1336  for (unsigned int k=0; k<N; ++k)
1337  std::swap (tmp.data[j][k], tmp.data[r][k]);
1338 
1339  std::swap (p[j], p[r]);
1340  }
1341 
1342  // Transformation
1343  const Number hr = 1./tmp.data[j][j];
1344  tmp.data[j][j] = hr;
1345  for (unsigned int k=0; k<N; ++k)
1346  {
1347  if (k==j) continue;
1348  for (unsigned int i=0; i<N; ++i)
1349  {
1350  if (i==j) continue;
1351  tmp.data[i][k] -= tmp.data[i][j]*tmp.data[j][k]*hr;
1352  }
1353  }
1354  for (unsigned int i=0; i<N; ++i)
1355  {
1356  tmp.data[i][j] *= hr;
1357  tmp.data[j][i] *= -hr;
1358  }
1359  tmp.data[j][j] = hr;
1360  }
1361 
1362  // Column interchange
1363  Number hv[N];
1364  for (unsigned int i=0; i<N; ++i)
1365  {
1366  for (unsigned int k=0; k<N; ++k)
1367  hv[p[k]] = tmp.data[i][k];
1368  for (unsigned int k=0; k<N; ++k)
1369  tmp.data[i][k] = hv[k];
1370  }
1371 
1372  // Scale rows and columns. The mult matrix
1373  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1374  for (unsigned int i=3; i<6; ++i)
1375  for (unsigned int j=0; j<3; ++j)
1376  tmp.data[i][j] /= 2;
1377 
1378  for (unsigned int i=0; i<3; ++i)
1379  for (unsigned int j=3; j<6; ++j)
1380  tmp.data[i][j] /= 2;
1381 
1382  for (unsigned int i=3; i<6; ++i)
1383  for (unsigned int j=3; j<6; ++j)
1384  tmp.data[i][j] /= 4;
1385 
1386  return tmp;
1387  }
1388  };
1389 
1390  }
1391 }
1392 
1393 
1394 
1395 template <int rank_, int dim, typename Number>
1396 inline
1398 operator Tensor<rank_,dim,Number> () const
1399 {
1400  return internal::SymmetricTensor::convert_to_tensor (*this);
1401 }
1402 
1403 
1404 
1405 template <int rank_, int dim, typename Number>
1406 inline
1407 bool
1409 (const SymmetricTensor<rank_,dim,Number> &t) const
1410 {
1411  return data == t.data;
1412 }
1413 
1414 
1415 
1416 template <int rank_, int dim, typename Number>
1417 inline
1418 bool
1419 SymmetricTensor<rank_,dim,Number>::operator !=
1420 (const SymmetricTensor<rank_,dim,Number> &t) const
1421 {
1422  return data != t.data;
1423 }
1424 
1425 
1426 
1427 template <int rank_, int dim, typename Number>
1428 template <typename OtherNumber>
1429 inline
1431 SymmetricTensor<rank_,dim,Number>::operator +=
1433 {
1434  data += t.data;
1435  return *this;
1436 }
1437 
1438 
1439 
1440 template <int rank_, int dim, typename Number>
1441 template <typename OtherNumber>
1442 inline
1444 SymmetricTensor<rank_,dim,Number>::operator -=
1446 {
1447  data -= t.data;
1448  return *this;
1449 }
1450 
1451 
1452 
1453 template <int rank_, int dim, typename Number>
1454 template <typename OtherNumber>
1455 inline
1458 {
1459  data *= d;
1460  return *this;
1461 }
1462 
1463 
1464 
1465 template <int rank_, int dim, typename Number>
1466 template <typename OtherNumber>
1467 inline
1470 {
1471  data /= d;
1472  return *this;
1473 }
1474 
1475 
1476 
1477 template <int rank_, int dim, typename Number>
1478 inline
1481 {
1482  SymmetricTensor tmp = *this;
1483  tmp.data = -tmp.data;
1484  return tmp;
1485 }
1486 
1487 
1488 
1489 template <int rank_, int dim, typename Number>
1490 inline
1491 void
1493 {
1494  data.clear ();
1495 }
1496 
1497 
1498 
1499 template <int rank_, int dim, typename Number>
1500 inline
1501 std::size_t
1503 {
1504  // all memory consists of statically allocated memory of the current
1505  // object, no pointers
1506  return sizeof(SymmetricTensor<rank_,dim,Number>);
1507 }
1508 
1509 
1510 
1511 namespace internal
1512 {
1513 
1514  template <int dim, typename Number, typename OtherNumber = Number>
1515  inline
1516  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type
1517  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1518  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1519  {
1520  typedef typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type result_type;
1521 
1522  switch (dim)
1523  {
1524  case 1:
1525  return data[0] * sdata[0];
1526  default:
1527  // Start with the non-diagonal part to avoid some multiplications by
1528  // 2.
1529 
1530  result_type sum = data[dim] * sdata[dim];
1531  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1532  sum += data[d] * sdata[d];
1533  sum += sum; // sum = sum * 2.;
1534 
1535  // Now add the contributions from the diagonal
1536  for (unsigned int d=0; d<dim; ++d)
1537  sum += data[d] * sdata[d];
1538  return sum;
1539  }
1540  }
1541 
1542 
1543 
1544  template <int dim, typename Number, typename OtherNumber = Number>
1545  inline
1546  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type
1547  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1548  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1549  {
1550  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type result_type;
1551  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::value_type value_type;
1552 
1553  const unsigned int data_dim =
1554  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1555  value_type tmp [data_dim];
1556  for (unsigned int i=0; i<data_dim; ++i)
1557  tmp[i] = perform_double_contraction<dim,Number,OtherNumber>(data[i], sdata);
1558  return result_type(tmp);
1559  }
1560 
1561 
1562 
1563  template <int dim, typename Number, typename OtherNumber = Number>
1564  inline
1565  typename SymmetricTensorAccessors::StorageType<2,dim,
1566  typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type
1567  >::base_tensor_type
1568  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1569  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1570  {
1571  typedef typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type value_type;
1572  typedef typename SymmetricTensorAccessors::StorageType<2,dim,value_type>::base_tensor_type base_tensor_type;
1573 
1574  base_tensor_type tmp;
1575  for (unsigned int i=0; i<tmp.dimension; ++i)
1576  {
1577  // Start with the non-diagonal part
1578  value_type sum = data[dim] * sdata[dim][i];
1579  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1580  sum += data[d] * sdata[d][i];
1581  sum += sum; // sum = sum * 2.;
1582 
1583  // Now add the contributions from the diagonal
1584  for (unsigned int d=0; d<dim; ++d)
1585  sum += data[d] * sdata[d][i];
1586  tmp[i] = sum;
1587  }
1588  return tmp;
1589  }
1590 
1591 
1592 
1593  template <int dim, typename Number, typename OtherNumber = Number>
1594  inline
1595  typename SymmetricTensorAccessors::StorageType<4,dim,
1596  typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type
1597  >::base_tensor_type
1598  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1599  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1600  {
1601  typedef typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type value_type;
1602  typedef typename SymmetricTensorAccessors::StorageType<4,dim,value_type>::base_tensor_type base_tensor_type;
1603 
1604  const unsigned int data_dim =
1605  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1606  base_tensor_type tmp;
1607  for (unsigned int i=0; i<data_dim; ++i)
1608  for (unsigned int j=0; j<data_dim; ++j)
1609  {
1610  // Start with the non-diagonal part
1611  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1612  tmp[i][j] += data[i][d] * sdata[d][j];
1613  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1614 
1615  // Now add the contributions from the diagonal
1616  for (unsigned int d=0; d<dim; ++d)
1617  tmp[i][j] += data[i][d] * sdata[d][j];
1618  }
1619  return tmp;
1620  }
1621 
1622 } // end of namespace internal
1623 
1624 
1625 
1626 template <int rank_, int dim, typename Number>
1627 template <typename OtherNumber>
1628 inline
1631 {
1632  // need to have two different function calls
1633  // because a scalar and rank-2 tensor are not
1634  // the same data type (see internal function
1635  // above)
1636  return internal::perform_double_contraction<dim,Number,OtherNumber> (data, s.data);
1637 }
1638 
1639 
1640 
1641 template <int rank_, int dim, typename Number>
1642 template <typename OtherNumber>
1643 inline
1646 {
1649  tmp.data = internal::perform_double_contraction<dim,Number,OtherNumber> (data,s.data);
1650  return tmp;
1651 }
1652 
1653 
1654 
1655 // internal namespace to switch between the
1656 // access of different tensors. There used to
1657 // be explicit instantiations before for
1658 // different ranks and dimensions, but since
1659 // we now allow for templates on the data
1660 // type, and since we cannot partially
1661 // specialize the implementation, this got
1662 // into a separate namespace
1663 namespace internal
1664 {
1665  template <int dim, typename Number>
1666  inline
1667  Number &
1668  symmetric_tensor_access (const TableIndices<2> &indices,
1670  {
1671  // 1d is very simple and done first
1672  if (dim == 1)
1673  return data[0];
1674 
1675  // first treat the main diagonal elements, which are stored consecutively
1676  // at the beginning
1677  if (indices[0] == indices[1])
1678  return data[indices[0]];
1679 
1680  // the rest is messier and requires a few switches.
1681  switch (dim)
1682  {
1683  case 2:
1684  // at least for the 2x2 case it is reasonably simple
1685  Assert (((indices[0]==1) && (indices[1]==0)) ||
1686  ((indices[0]==0) && (indices[1]==1)),
1687  ExcInternalError());
1688  return data[2];
1689 
1690  default:
1691  // to do the rest, sort our indices before comparing
1692  {
1693  TableIndices<2> sorted_indices (indices);
1694  sorted_indices.sort ();
1695 
1696  for (unsigned int d=0, c=0; d<dim; ++d)
1697  for (unsigned int e=d+1; e<dim; ++e, ++c)
1698  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1699  return data[dim+c];
1700  Assert (false, ExcInternalError());
1701  }
1702  }
1703 
1704  static Number dummy_but_referenceable = Number();
1705  return dummy_but_referenceable;
1706  }
1707 
1708 
1709 
1710  template <int dim, typename Number>
1711  inline
1712  const Number &
1713  symmetric_tensor_access (const TableIndices<2> &indices,
1715  {
1716  // 1d is very simple and done first
1717  if (dim == 1)
1718  return data[0];
1719 
1720  // first treat the main diagonal elements, which are stored consecutively
1721  // at the beginning
1722  if (indices[0] == indices[1])
1723  return data[indices[0]];
1724 
1725  // the rest is messier and requires a few switches.
1726  switch (dim)
1727  {
1728  case 2:
1729  // at least for the 2x2 case it is reasonably simple
1730  Assert (((indices[0]==1) && (indices[1]==0)) ||
1731  ((indices[0]==0) && (indices[1]==1)),
1732  ExcInternalError());
1733  return data[2];
1734 
1735  default:
1736  // to do the rest, sort our indices before comparing
1737  {
1738  TableIndices<2> sorted_indices (indices);
1739  sorted_indices.sort ();
1740 
1741  for (unsigned int d=0, c=0; d<dim; ++d)
1742  for (unsigned int e=d+1; e<dim; ++e, ++c)
1743  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1744  return data[dim+c];
1745  Assert (false, ExcInternalError());
1746  }
1747  }
1748 
1749  static Number dummy_but_referenceable = Number();
1750  return dummy_but_referenceable;
1751  }
1752 
1753 
1754 
1755  template <int dim, typename Number>
1756  inline
1757  Number &
1758  symmetric_tensor_access (const TableIndices<4> &indices,
1760  {
1761  switch (dim)
1762  {
1763  case 1:
1764  return data[0][0];
1765 
1766  case 2:
1767  // each entry of the tensor can be
1768  // thought of as an entry in a
1769  // matrix that maps the rolled-out
1770  // rank-2 tensors into rolled-out
1771  // rank-2 tensors. this is the
1772  // format in which we store rank-4
1773  // tensors. determine which
1774  // position the present entry is
1775  // stored in
1776  {
1777  unsigned int base_index[2] ;
1778  if ((indices[0] == 0) && (indices[1] == 0))
1779  base_index[0] = 0;
1780  else if ((indices[0] == 1) && (indices[1] == 1))
1781  base_index[0] = 1;
1782  else
1783  base_index[0] = 2;
1784 
1785  if ((indices[2] == 0) && (indices[3] == 0))
1786  base_index[1] = 0;
1787  else if ((indices[2] == 1) && (indices[3] == 1))
1788  base_index[1] = 1;
1789  else
1790  base_index[1] = 2;
1791 
1792  return data[base_index[0]][base_index[1]];
1793  }
1794 
1795  case 3:
1796  // each entry of the tensor can be
1797  // thought of as an entry in a
1798  // matrix that maps the rolled-out
1799  // rank-2 tensors into rolled-out
1800  // rank-2 tensors. this is the
1801  // format in which we store rank-4
1802  // tensors. determine which
1803  // position the present entry is
1804  // stored in
1805  {
1806  unsigned int base_index[2] ;
1807  if ((indices[0] == 0) && (indices[1] == 0))
1808  base_index[0] = 0;
1809  else if ((indices[0] == 1) && (indices[1] == 1))
1810  base_index[0] = 1;
1811  else if ((indices[0] == 2) && (indices[1] == 2))
1812  base_index[0] = 2;
1813  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1814  ((indices[0] == 1) && (indices[1] == 0)))
1815  base_index[0] = 3;
1816  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1817  ((indices[0] == 2) && (indices[1] == 0)))
1818  base_index[0] = 4;
1819  else
1820  {
1821  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1822  ((indices[0] == 2) && (indices[1] == 1)),
1823  ExcInternalError());
1824  base_index[0] = 5;
1825  }
1826 
1827  if ((indices[2] == 0) && (indices[3] == 0))
1828  base_index[1] = 0;
1829  else if ((indices[2] == 1) && (indices[3] == 1))
1830  base_index[1] = 1;
1831  else if ((indices[2] == 2) && (indices[3] == 2))
1832  base_index[1] = 2;
1833  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1834  ((indices[2] == 1) && (indices[3] == 0)))
1835  base_index[1] = 3;
1836  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1837  ((indices[2] == 2) && (indices[3] == 0)))
1838  base_index[1] = 4;
1839  else
1840  {
1841  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1842  ((indices[2] == 2) && (indices[3] == 1)),
1843  ExcInternalError());
1844  base_index[1] = 5;
1845  }
1846 
1847  return data[base_index[0]][base_index[1]];
1848  }
1849 
1850  default:
1851  Assert (false, ExcNotImplemented());
1852  }
1853 
1854  static Number dummy;
1855  return dummy;
1856  }
1857 
1858 
1859  template <int dim, typename Number>
1860  inline
1861  const Number &
1862  symmetric_tensor_access (const TableIndices<4> &indices,
1864  {
1865  switch (dim)
1866  {
1867  case 1:
1868  return data[0][0];
1869 
1870  case 2:
1871  // each entry of the tensor can be
1872  // thought of as an entry in a
1873  // matrix that maps the rolled-out
1874  // rank-2 tensors into rolled-out
1875  // rank-2 tensors. this is the
1876  // format in which we store rank-4
1877  // tensors. determine which
1878  // position the present entry is
1879  // stored in
1880  {
1881  unsigned int base_index[2] ;
1882  if ((indices[0] == 0) && (indices[1] == 0))
1883  base_index[0] = 0;
1884  else if ((indices[0] == 1) && (indices[1] == 1))
1885  base_index[0] = 1;
1886  else
1887  base_index[0] = 2;
1888 
1889  if ((indices[2] == 0) && (indices[3] == 0))
1890  base_index[1] = 0;
1891  else if ((indices[2] == 1) && (indices[3] == 1))
1892  base_index[1] = 1;
1893  else
1894  base_index[1] = 2;
1895 
1896  return data[base_index[0]][base_index[1]];
1897  }
1898 
1899  case 3:
1900  // each entry of the tensor can be
1901  // thought of as an entry in a
1902  // matrix that maps the rolled-out
1903  // rank-2 tensors into rolled-out
1904  // rank-2 tensors. this is the
1905  // format in which we store rank-4
1906  // tensors. determine which
1907  // position the present entry is
1908  // stored in
1909  {
1910  unsigned int base_index[2] ;
1911  if ((indices[0] == 0) && (indices[1] == 0))
1912  base_index[0] = 0;
1913  else if ((indices[0] == 1) && (indices[1] == 1))
1914  base_index[0] = 1;
1915  else if ((indices[0] == 2) && (indices[1] == 2))
1916  base_index[0] = 2;
1917  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1918  ((indices[0] == 1) && (indices[1] == 0)))
1919  base_index[0] = 3;
1920  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1921  ((indices[0] == 2) && (indices[1] == 0)))
1922  base_index[0] = 4;
1923  else
1924  {
1925  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1926  ((indices[0] == 2) && (indices[1] == 1)),
1927  ExcInternalError());
1928  base_index[0] = 5;
1929  }
1930 
1931  if ((indices[2] == 0) && (indices[3] == 0))
1932  base_index[1] = 0;
1933  else if ((indices[2] == 1) && (indices[3] == 1))
1934  base_index[1] = 1;
1935  else if ((indices[2] == 2) && (indices[3] == 2))
1936  base_index[1] = 2;
1937  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1938  ((indices[2] == 1) && (indices[3] == 0)))
1939  base_index[1] = 3;
1940  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1941  ((indices[2] == 2) && (indices[3] == 0)))
1942  base_index[1] = 4;
1943  else
1944  {
1945  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1946  ((indices[2] == 2) && (indices[3] == 1)),
1947  ExcInternalError());
1948  base_index[1] = 5;
1949  }
1950 
1951  return data[base_index[0]][base_index[1]];
1952  }
1953 
1954  default:
1955  Assert (false, ExcNotImplemented());
1956  }
1957 
1958  static Number dummy;
1959  return dummy;
1960  }
1961 
1962 } // end of namespace internal
1963 
1964 
1965 
1966 template <int rank_, int dim, typename Number>
1967 inline
1968 Number &
1970 {
1971  for (unsigned int r=0; r<rank; ++r)
1972  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1973  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1974 }
1975 
1976 
1977 
1978 template <int rank_, int dim, typename Number>
1979 inline
1980 const Number &
1982 (const TableIndices<rank_> &indices) const
1983 {
1984  for (unsigned int r=0; r<rank; ++r)
1985  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1986  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1987 }
1988 
1989 
1990 
1991 namespace internal
1992 {
1993  namespace SymmetricTensor
1994  {
1995  template <int rank_>
1997  get_partially_filled_indices (const unsigned int row,
1998  const std::integral_constant<int, 2> &)
1999  {
2000  return TableIndices<rank_> (row,
2002 
2003  }
2004 
2005 
2006  template <int rank_>
2008  get_partially_filled_indices (const unsigned int row,
2009  const std::integral_constant<int, 4> &)
2010  {
2011  return TableIndices<rank_> (row,
2015 
2016  }
2017  }
2018 }
2019 
2020 
2021 template <int rank_, int dim, typename Number>
2022 internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
2023 SymmetricTensor<rank_,dim,Number>::operator [] (const unsigned int row) const
2024 {
2025  return
2026  internal::SymmetricTensorAccessors::
2027  Accessor<rank_,dim,true,rank_-1,Number> (*this,
2028  internal::SymmetricTensor::get_partially_filled_indices<rank_> (row,
2029  std::integral_constant<int, rank_>()));
2030 }
2031 
2032 
2033 
2034 template <int rank_, int dim, typename Number>
2035 internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
2037 {
2038  return
2039  internal::SymmetricTensorAccessors::
2040  Accessor<rank_,dim,false,rank_-1,Number> (*this,
2041  internal::SymmetricTensor::get_partially_filled_indices<rank_> (row,
2042  std::integral_constant<int, rank_>()));
2043 }
2044 
2045 
2046 
2047 template <int rank_, int dim, typename Number>
2048 inline
2049 const Number &
2051 {
2052  return operator()(indices);
2053 }
2054 
2055 
2056 
2057 template <int rank_, int dim, typename Number>
2058 inline
2059 Number &
2061 {
2062  return operator()(indices);
2063 }
2064 
2065 
2066 
2067 template <int rank_, int dim, typename Number>
2068 inline
2069 Number *
2071 {
2072  return std::addressof(this->access_raw_entry(0));
2073 }
2074 
2075 
2076 
2077 template <int rank_, int dim, typename Number>
2078 inline
2079 const Number *
2081 {
2082  return std::addressof(this->access_raw_entry(0));
2083 }
2084 
2085 
2086 
2087 template <int rank_, int dim, typename Number>
2088 inline
2089 Number *
2091 {
2092  return begin_raw()+n_independent_components;
2093 }
2094 
2095 
2096 
2097 template <int rank_, int dim, typename Number>
2098 inline
2099 const Number *
2101 {
2102  return begin_raw()+n_independent_components;
2103 }
2104 
2105 
2106 
2107 namespace internal
2108 {
2109  namespace SymmetricTensor
2110  {
2111  template <int dim, typename Number>
2112  unsigned int
2113  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
2114  const unsigned int index)
2115  {
2116  return index;
2117  }
2118 
2119 
2120  template <int dim, typename Number>
2122  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
2123  const unsigned int index)
2124  {
2125  return
2128  }
2129 
2130  }
2131 }
2132 
2133 
2134 
2135 template <int rank_, int dim, typename Number>
2136 inline
2137 const Number &
2138 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index) const
2139 {
2140  AssertIndexRange (index, n_independent_components);
2141  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2142 }
2143 
2144 
2145 
2146 template <int rank_, int dim, typename Number>
2147 inline
2148 Number &
2149 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index)
2150 {
2151  AssertIndexRange (index, n_independent_components);
2152  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2153 }
2154 
2155 
2156 
2157 namespace internal
2158 {
2159  template <int dim, typename Number>
2160  inline
2163  {
2164  switch (dim)
2165  {
2166  case 1:
2167  return numbers::NumberTraits<Number>::abs(data[0]);
2168 
2169  case 2:
2170  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2173 
2174  case 3:
2175  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2181 
2182  default:
2183  {
2184  typename numbers::NumberTraits<Number>::real_type return_value
2186 
2187  for (unsigned int d=0; d<dim; ++d)
2188  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
2189  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
2190  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2191 
2192  return std::sqrt(return_value);
2193  }
2194  }
2195  }
2196 
2197 
2198 
2199  template <int dim, typename Number>
2200  inline
2203  {
2204  switch (dim)
2205  {
2206  case 1:
2207  return numbers::NumberTraits<Number>::abs (data[0][0]);
2208 
2209  default:
2210  {
2211  typename numbers::NumberTraits<Number>::real_type return_value
2213 
2214  const unsigned int n_independent_components = data.dimension;
2215 
2216  for (unsigned int i=0; i<dim; ++i)
2217  for (unsigned int j=0; j<dim; ++j)
2218  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
2219  for (unsigned int i=0; i<dim; ++i)
2220  for (unsigned int j=dim; j<n_independent_components; ++j)
2221  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2222  for (unsigned int i=dim; i<n_independent_components; ++i)
2223  for (unsigned int j=0; j<dim; ++j)
2224  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2225  for (unsigned int i=dim; i<n_independent_components; ++i)
2226  for (unsigned int j=dim; j<n_independent_components; ++j)
2227  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2228 
2229  return std::sqrt(return_value);
2230  }
2231  }
2232  }
2233 
2234 } // end of namespace internal
2235 
2236 
2237 
2238 template <int rank_, int dim, typename Number>
2239 inline
2242 {
2243  return internal::compute_norm<dim,Number> (data);
2244 }
2245 
2246 
2247 
2248 namespace internal
2249 {
2250  namespace SymmetricTensor
2251  {
2252  namespace
2253  {
2254  // a function to do the unrolling from a set of indices to a
2255  // scalar index into the array in which we store the elements of
2256  // a symmetric tensor
2257  //
2258  // this function is for rank-2 tensors
2259  template <int dim>
2260  inline
2261  unsigned int
2262  component_to_unrolled_index
2263  (const TableIndices<2> &indices)
2264  {
2265  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2266  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2267 
2268  switch (dim)
2269  {
2270  case 1:
2271  {
2272  return 0;
2273  }
2274 
2275  case 2:
2276  {
2277  static const unsigned int table[2][2] = {{0, 2},
2278  {2, 1}
2279  };
2280  return table[indices[0]][indices[1]];
2281  }
2282 
2283  case 3:
2284  {
2285  static const unsigned int table[3][3] = {{0, 3, 4},
2286  {3, 1, 5},
2287  {4, 5, 2}
2288  };
2289  return table[indices[0]][indices[1]];
2290  }
2291 
2292  case 4:
2293  {
2294  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2295  {4, 1, 7, 8},
2296  {5, 7, 2, 9},
2297  {6, 8, 9, 3}
2298  };
2299  return table[indices[0]][indices[1]];
2300  }
2301 
2302  default:
2303  // for the remainder, manually figure out the numbering
2304  {
2305  if (indices[0] == indices[1])
2306  return indices[0];
2307 
2308  TableIndices<2> sorted_indices (indices);
2309  sorted_indices.sort ();
2310 
2311  for (unsigned int d=0, c=0; d<dim; ++d)
2312  for (unsigned int e=d+1; e<dim; ++e, ++c)
2313  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
2314  return dim+c;
2315 
2316  // should never get here:
2317  Assert(false, ExcInternalError());
2318  return 0;
2319  }
2320  }
2321  }
2322 
2323  // a function to do the unrolling from a set of indices to a
2324  // scalar index into the array in which we store the elements of
2325  // a symmetric tensor
2326  //
2327  // this function is for tensors of ranks not already handled
2328  // above
2329  template <int dim, int rank_>
2330  inline
2331  unsigned int
2332  component_to_unrolled_index
2333  (const TableIndices<rank_> &indices)
2334  {
2335  (void)indices;
2336  Assert (false, ExcNotImplemented());
2338  }
2339  }
2340  }
2341 }
2342 
2343 
2344 template <int rank_, int dim, typename Number>
2345 inline
2346 unsigned int
2348 (const TableIndices<rank_> &indices)
2349 {
2350  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
2351 }
2352 
2353 
2354 
2355 namespace internal
2356 {
2357  namespace SymmetricTensor
2358  {
2359  namespace
2360  {
2361  // a function to do the inverse of the unrolling from a set of
2362  // indices to a scalar index into the array in which we store
2363  // the elements of a symmetric tensor. in other words, it goes
2364  // from the scalar index into the array to a set of indices of
2365  // the tensor
2366  //
2367  // this function is for rank-2 tensors
2368  template <int dim>
2369  inline
2371  unrolled_to_component_indices
2372  (const unsigned int i,
2373  const std::integral_constant<int, 2> &)
2374  {
2377  switch (dim)
2378  {
2379  case 1:
2380  {
2381  return TableIndices<2>(0,0);
2382  }
2383 
2384  case 2:
2385  {
2386  const TableIndices<2> table[3] =
2387  {
2388  TableIndices<2> (0,0),
2389  TableIndices<2> (1,1),
2390  TableIndices<2> (0,1)
2391  };
2392  return table[i];
2393  }
2394 
2395  case 3:
2396  {
2397  const TableIndices<2> table[6] =
2398  {
2399  TableIndices<2> (0,0),
2400  TableIndices<2> (1,1),
2401  TableIndices<2> (2,2),
2402  TableIndices<2> (0,1),
2403  TableIndices<2> (0,2),
2404  TableIndices<2> (1,2)
2405  };
2406  return table[i];
2407  }
2408 
2409  default:
2410  if (i<dim)
2411  return TableIndices<2> (i,i);
2412 
2413  for (unsigned int d=0, c=0; d<dim; ++d)
2414  for (unsigned int e=d+1; e<dim; ++e, ++c)
2415  if (c==i)
2416  return TableIndices<2>(d,e);
2417 
2418  // should never get here:
2419  Assert(false, ExcInternalError());
2420  return TableIndices<2>(0, 0);
2421  }
2422  }
2423 
2424  // a function to do the inverse of the unrolling from a set of
2425  // indices to a scalar index into the array in which we store
2426  // the elements of a symmetric tensor. in other words, it goes
2427  // from the scalar index into the array to a set of indices of
2428  // the tensor
2429  //
2430  // this function is for tensors of a rank not already handled
2431  // above
2432  template <int dim, int rank_>
2433  inline
2435  unrolled_to_component_indices
2436  (const unsigned int i,
2437  const std::integral_constant<int, rank_> &)
2438  {
2439  (void)i;
2442  Assert (false, ExcNotImplemented());
2443  return TableIndices<rank_>();
2444  }
2445 
2446  }
2447  }
2448 }
2449 
2450 template <int rank_, int dim, typename Number>
2451 inline
2454 (const unsigned int i)
2455 {
2456  return
2457  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2458  std::integral_constant<int, rank_>());
2459 }
2460 
2461 
2462 
2463 template <int rank_, int dim, typename Number>
2464 template <class Archive>
2465 inline
2466 void
2467 SymmetricTensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
2468 {
2469  ar &data;
2470 }
2471 
2472 
2473 #endif // DOXYGEN
2474 
2475 /* ----------------- Non-member functions operating on tensors. ------------ */
2476 
2477 
2490 template <int rank_, int dim, typename Number, typename OtherNumber>
2491 inline
2495 {
2497  tmp += right;
2498  return tmp;
2499 }
2500 
2501 
2514 template <int rank_, int dim, typename Number, typename OtherNumber>
2515 inline
2519 {
2521  tmp -= right;
2522  return tmp;
2523 }
2524 
2525 
2533 template <int rank_, int dim, typename Number, typename OtherNumber>
2534 inline
2537  const Tensor<rank_, dim, OtherNumber> &right)
2538 {
2539  return Tensor<rank_, dim, Number>(left) + right;
2540 }
2541 
2542 
2550 template <int rank_, int dim, typename Number, typename OtherNumber>
2551 inline
2555 {
2556  return left + Tensor<rank_, dim, OtherNumber>(right);
2557 }
2558 
2559 
2567 template <int rank_, int dim, typename Number, typename OtherNumber>
2568 inline
2571  const Tensor<rank_, dim, OtherNumber> &right)
2572 {
2573  return Tensor<rank_, dim, Number>(left) - right;
2574 }
2575 
2576 
2584 template <int rank_, int dim, typename Number, typename OtherNumber>
2585 inline
2589 {
2590  return left - Tensor<rank_, dim, OtherNumber>(right);
2591 }
2592 
2593 
2594 
2608 template <int dim, typename Number>
2609 inline
2611 {
2612  switch (dim)
2613  {
2614  case 1:
2615  return t.data[0];
2616  case 2:
2617  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2618  case 3:
2619  {
2620  // in analogy to general tensors, but
2621  // there's something to be simplified for
2622  // the present case
2623  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2624  return ( tmp + tmp
2625  +t.data[0]*t.data[1]*t.data[2]
2626  -t.data[0]*t.data[5]*t.data[5]
2627  -t.data[1]*t.data[4]*t.data[4]
2628  -t.data[2]*t.data[3]*t.data[3]);
2629  }
2630  default:
2631  Assert (false, ExcNotImplemented());
2633  }
2634 }
2635 
2636 
2637 
2647 template <int dim, typename Number>
2648 inline
2650 {
2651  return determinant (t);
2652 }
2653 
2654 
2655 
2663 template <int dim, typename Number>
2665 {
2666  Number t = d.data[0];
2667  for (unsigned int i=1; i<dim; ++i)
2668  t += d.data[i];
2669  return t;
2670 }
2671 
2672 
2682 template <int dim, typename Number>
2683 inline
2685 {
2686  return trace (t);
2687 }
2688 
2689 
2702 template <typename Number>
2703 inline
2705 {
2707 }
2708 
2709 
2710 
2731 template <typename Number>
2732 inline
2734 {
2735  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2736 }
2737 
2738 
2739 
2749 template <typename Number>
2750 inline
2752 {
2753  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2754  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2755 }
2756 
2757 
2758 
2767 template <typename Number>
2768 std::array<Number,1>
2769 eigenvalues (const SymmetricTensor<2,1,Number> &T);
2770 
2771 
2772 
2794 template <typename Number>
2795 std::array<Number,2>
2796 eigenvalues (const SymmetricTensor<2,2,Number> &T);
2797 
2798 
2799 
2819 template <typename Number>
2820 std::array<Number,3>
2821 eigenvalues (const SymmetricTensor<2,3,Number> &T);
2822 
2823 
2824 
2825 namespace internal
2826 {
2827  namespace SymmetricTensor
2828  {
2865  template <int dim, typename Number>
2866  void
2867  tridiagonalize (const ::SymmetricTensor<2,dim,Number> &A,
2868  ::Tensor<2,dim,Number> &Q,
2869  std::array<Number,dim> &d,
2870  std::array<Number,dim-1> &e);
2871 
2872 
2873 
2912  template <int dim, typename Number>
2913  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2914  ql_implicit_shifts (const ::SymmetricTensor<2,dim,Number> &A);
2915 
2916 
2917 
2956  template <int dim, typename Number>
2957  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2959 
2960 
2961 
2977  template <typename Number>
2978  std::array<std::pair<Number, Tensor<1,2,Number> >,2>
2979  hybrid (const ::SymmetricTensor<2,2,Number> &A);
2980 
2981 
2982 
3015  template <typename Number>
3016  std::array<std::pair<Number, Tensor<1,3,Number> >,3>
3017  hybrid (const ::SymmetricTensor<2,3,Number> &A);
3018 
3019  namespace
3020  {
3021 
3026  template <int dim, typename Number>
3027  struct SortEigenValuesVectors
3028  {
3029  typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
3030  bool operator() (const EigValsVecs &lhs,
3031  const EigValsVecs &rhs)
3032  {
3033  return lhs.first > rhs.first;
3034  }
3035  };
3036 
3037  }
3038 
3039  } // namespace SymmetricTensor
3040 
3041 } // namespace internal
3042 
3043 
3044 
3045 // The line below is to ensure that doxygen puts the full description
3046 // of this global enumeration into the documentation
3047 // See https://stackoverflow.com/a/1717984
3075 {
3085  hybrid,
3095  ql_implicit_shifts,
3103  jacobi
3104 };
3105 
3106 
3107 
3114 template <typename Number>
3115 std::array<std::pair<Number, Tensor<1,1,Number> >,1>
3118 {
3119  return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} };
3120 }
3121 
3122 
3123 
3151 template <int dim, typename Number>
3152 std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
3155 {
3156  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
3157 
3158  switch (method)
3159  {
3161  eig_vals_vecs = internal::SymmetricTensor::hybrid(T);
3162  break;
3165  break;
3167  eig_vals_vecs = internal::SymmetricTensor::jacobi(T);
3168  break;
3169  default:
3170  AssertThrow(false, ExcNotImplemented());
3171  }
3172 
3173  // Sort in descending order before output.
3174  std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
3175  internal::SymmetricTensor::SortEigenValuesVectors<dim,Number>());
3176  return eig_vals_vecs;
3177 }
3178 
3179 
3180 
3190 template <int rank_, int dim, typename Number>
3191 inline
3194 {
3195  return t;
3196 }
3197 
3198 
3199 
3209 template <int dim, typename Number>
3210 inline
3213 {
3215 
3216  // subtract scaled trace from the diagonal
3217  const Number tr = trace(t) / dim;
3218  for (unsigned int i=0; i<dim; ++i)
3219  tmp.data[i] -= tr;
3220 
3221  return tmp;
3222 }
3223 
3224 
3225 
3233 template <int dim, typename Number>
3234 inline
3236 unit_symmetric_tensor ()
3237 {
3238  // create a default constructed matrix filled with
3239  // zeros, then set the diagonal elements to one
3241  switch (dim)
3242  {
3243  case 1:
3244  tmp.data[0] = 1;
3245  break;
3246  case 2:
3247  tmp.data[0] = tmp.data[1] = 1;
3248  break;
3249  case 3:
3250  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3251  break;
3252  default:
3253  for (unsigned int d=0; d<dim; ++d)
3254  tmp.data[d] = 1;
3255  }
3256  return tmp;
3257 }
3258 
3259 
3260 
3269 template <int dim>
3270 inline
3272 unit_symmetric_tensor ()
3273 {
3274  return unit_symmetric_tensor<dim,double>();
3275 }
3276 
3277 
3278 
3293 template <int dim, typename Number>
3294 inline
3296 deviator_tensor ()
3297 {
3299 
3300  // fill the elements treating the diagonal
3301  for (unsigned int i=0; i<dim; ++i)
3302  for (unsigned int j=0; j<dim; ++j)
3303  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
3304 
3305  // then fill the ones that copy over the
3306  // non-diagonal elements. note that during
3307  // the double-contraction, we handle the
3308  // off-diagonal elements twice, so simply
3309  // copying requires a weight of 1/2
3310  for (unsigned int i=dim;
3311  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3312  ++i)
3313  tmp.data[i][i] = 0.5;
3314 
3315  return tmp;
3316 }
3317 
3318 
3319 
3334 template <int dim>
3335 inline
3337 deviator_tensor ()
3338 {
3339  return deviator_tensor<dim,double>();
3340 }
3341 
3342 
3343 
3366 template <int dim, typename Number>
3367 inline
3369 identity_tensor ()
3370 {
3372 
3373  // fill the elements treating the diagonal
3374  for (unsigned int i=0; i<dim; ++i)
3375  tmp.data[i][i] = 1;
3376 
3377  // then fill the ones that copy over the
3378  // non-diagonal elements. note that during
3379  // the double-contraction, we handle the
3380  // off-diagonal elements twice, so simply
3381  // copying requires a weight of 1/2
3382  for (unsigned int i=dim;
3383  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3384  ++i)
3385  tmp.data[i][i] = 0.5;
3386 
3387  return tmp;
3388 }
3389 
3390 
3391 
3413 template <int dim>
3414 inline
3416 identity_tensor ()
3417 {
3418  return identity_tensor<dim,double>();
3419 }
3420 
3421 
3422 
3433 template <int dim, typename Number>
3434 inline
3437 {
3439 }
3440 
3441 
3442 
3454 template <int dim, typename Number>
3455 inline
3458 {
3460 }
3461 
3462 
3463 
3478 template <int dim, typename Number>
3479 inline
3483 {
3485 
3486  // fill only the elements really needed
3487  for (unsigned int i=0; i<dim; ++i)
3488  for (unsigned int j=i; j<dim; ++j)
3489  for (unsigned int k=0; k<dim; ++k)
3490  for (unsigned int l=k; l<dim; ++l)
3491  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3492 
3493  return tmp;
3494 }
3495 
3496 
3497 
3506 template <int dim,typename Number>
3507 inline
3510 {
3511  Number array[(dim*dim+dim)/2];
3512  for (unsigned int d=0; d<dim; ++d)
3513  array[d] = t[d][d];
3514  for (unsigned int d=0, c=0; d<dim; ++d)
3515  for (unsigned int e=d+1; e<dim; ++e, ++c)
3516  array[dim+c] = (t[d][e]+t[e][d])*0.5;
3517  return SymmetricTensor<2,dim,Number>(array);
3518 }
3519 
3520 
3521 
3529 template <int rank_, int dim, typename Number>
3530 inline
3533  const Number &factor)
3534 {
3536  tt *= factor;
3537  return tt;
3538 }
3539 
3540 
3541 
3549 template <int rank_, int dim, typename Number>
3550 inline
3552 operator * (const Number &factor,
3554 {
3555  // simply forward to the other operator
3556  return t*factor;
3557 }
3558 
3559 
3560 
3586 template <int rank_, int dim, typename Number, typename OtherNumber>
3587 inline
3590  const OtherNumber &factor)
3591 {
3592  // form the product. we have to convert the two factors into the final
3593  // type via explicit casts because, for awkward reasons, the C++
3594  // standard committee saw it fit to not define an
3595  // operator*(float,std::complex<double>)
3596  // (as well as with switched arguments and double<->float).
3597  typedef typename ProductType<Number,OtherNumber>::type product_type;
3599  // we used to shorten the following by 'tt *= product_type(factor);'
3600  // which requires that a converting constructor
3601  // 'product_type::product_type(const OtherNumber) is defined.
3602  // however, a user-defined constructor is not allowed for aggregates,
3603  // e.g. VectorizedArray. therefore, we work around this issue using a
3604  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3605  // which we assume to be defined.
3606  product_type new_factor;
3607  new_factor = factor;
3608  tt *= new_factor;
3609  return tt;
3610 }
3611 
3612 
3613 
3622 template <int rank_, int dim, typename Number, typename OtherNumber>
3623 inline
3625 operator * (const Number &factor,
3627 {
3628  // simply forward to the other operator with switched arguments
3629  return (t*factor);
3630 }
3631 
3632 
3633 
3639 template <int rank_, int dim, typename Number, typename OtherNumber>
3640 inline
3643  const OtherNumber &factor)
3644 {
3646  tt /= factor;
3647  return tt;
3648 }
3649 
3650 
3651 
3658 template <int rank_, int dim>
3659 inline
3661 operator * (const SymmetricTensor<rank_,dim> &t,
3662  const double factor)
3663 {
3665  tt *= factor;
3666  return tt;
3667 }
3668 
3669 
3670 
3677 template <int rank_, int dim>
3678 inline
3680 operator * (const double factor,
3681  const SymmetricTensor<rank_,dim> &t)
3682 {
3684  tt *= factor;
3685  return tt;
3686 }
3687 
3688 
3689 
3695 template <int rank_, int dim>
3696 inline
3698 operator / (const SymmetricTensor<rank_,dim> &t,
3699  const double factor)
3700 {
3702  tt /= factor;
3703  return tt;
3704 }
3705 
3715 template <int dim, typename Number, typename OtherNumber>
3716 inline
3717 typename ProductType<Number, OtherNumber>::type
3720 {
3721  return (t1*t2);
3722 }
3723 
3724 
3734 template <int dim, typename Number, typename OtherNumber>
3735 inline
3736 typename ProductType<Number, OtherNumber>::type
3738  const Tensor<2,dim,OtherNumber> &t2)
3739 {
3740  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<typename ProductType<Number, OtherNumber>::type>::value(0.0);
3741  for (unsigned int i=0; i<dim; ++i)
3742  for (unsigned int j=0; j<dim; ++j)
3743  s += t1[i][j] * t2[i][j];
3744  return s;
3745 }
3746 
3747 
3757 template <int dim, typename Number, typename OtherNumber>
3758 inline
3759 typename ProductType<Number, OtherNumber>::type
3762 {
3763  return scalar_product(t2, t1);
3764 }
3765 
3766 
3782 template <typename Number, typename OtherNumber>
3783 inline
3784 void
3785 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3786  const SymmetricTensor<4,1,Number> &t,
3788 {
3789  tmp[0][0] = t[0][0][0][0] * s[0][0];
3790 }
3791 
3792 
3793 
3809 template <typename Number, typename OtherNumber>
3810 inline
3811 void
3812 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3813  const SymmetricTensor<2,1,Number> &s,
3815 {
3816  tmp[0][0] = t[0][0][0][0] * s[0][0];
3817 }
3818 
3819 
3820 
3835 template <typename Number, typename OtherNumber>
3836 inline
3837 void
3838 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3839  const SymmetricTensor<4,2,Number> &t,
3841 {
3842  const unsigned int dim = 2;
3843 
3844  for (unsigned int i=0; i<dim; ++i)
3845  for (unsigned int j=i; j<dim; ++j)
3846  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3847  t[i][j][1][1] * s[1][1] +
3848  2 * t[i][j][0][1] * s[0][1];
3849 }
3850 
3851 
3852 
3868 template <typename Number, typename OtherNumber>
3869 inline
3870 void
3871 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3872  const SymmetricTensor<2,2,Number> &s,
3874 {
3875  const unsigned int dim = 2;
3876 
3877  for (unsigned int i=0; i<dim; ++i)
3878  for (unsigned int j=i; j<dim; ++j)
3879  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3880  s[1][1] * t[1][1][i][j] +
3881  2 * s[0][1] * t[0][1][i][j];
3882 }
3883 
3884 
3885 
3901 template <typename Number, typename OtherNumber>
3902 inline
3903 void
3904 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3905  const SymmetricTensor<4,3,Number> &t,
3907 {
3908  const unsigned int dim = 3;
3909 
3910  for (unsigned int i=0; i<dim; ++i)
3911  for (unsigned int j=i; j<dim; ++j)
3912  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3913  t[i][j][1][1] * s[1][1] +
3914  t[i][j][2][2] * s[2][2] +
3915  2 * t[i][j][0][1] * s[0][1] +
3916  2 * t[i][j][0][2] * s[0][2] +
3917  2 * t[i][j][1][2] * s[1][2];
3918 }
3919 
3920 
3921 
3937 template <typename Number, typename OtherNumber>
3938 inline
3939 void
3940 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3941  const SymmetricTensor<2,3,Number> &s,
3943 {
3944  const unsigned int dim = 3;
3945 
3946  for (unsigned int i=0; i<dim; ++i)
3947  for (unsigned int j=i; j<dim; ++j)
3948  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3949  s[1][1] * t[1][1][i][j] +
3950  s[2][2] * t[2][2][i][j] +
3951  2 * s[0][1] * t[0][1][i][j] +
3952  2 * s[0][2] * t[0][2][i][j] +
3953  2 * s[1][2] * t[1][2][i][j];
3954 }
3955 
3956 
3957 
3965 template <int dim, typename Number, typename OtherNumber>
3967 operator * (const SymmetricTensor<2,dim,Number> &src1,
3968  const Tensor<1,dim,OtherNumber> &src2)
3969 {
3971  for (unsigned int i=0; i<dim; ++i)
3972  for (unsigned int j=0; j<dim; ++j)
3973  dest[i] += src1[i][j] * src2[j];
3974  return dest;
3975 }
3976 
3977 
3985 template <int dim, typename Number, typename OtherNumber>
3987 operator * (const Tensor<1,dim,Number> &src1,
3989 {
3990  // this is easy for symmetric tensors:
3991  return src2 * src1;
3992 }
3993 
3994 
3995 
4016 template <int rank_1, int rank_2, int dim,
4017  typename Number, typename OtherNumber>
4018 inline DEAL_II_ALWAYS_INLINE
4019 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
4020 operator * (const Tensor<rank_1, dim, Number> &src1,
4022 {
4023  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
4024  const Tensor<rank_2, dim, OtherNumber> src2 (src2s);
4025  return src1*src2;
4026 }
4027 
4028 
4029 
4050 template <int rank_1, int rank_2, int dim,
4051  typename Number, typename OtherNumber>
4052 inline DEAL_II_ALWAYS_INLINE
4053 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
4054 operator * (const SymmetricTensor<rank_1, dim, Number> &src1s,
4056 {
4057  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
4058  const Tensor<rank_2, dim, OtherNumber> src1 (src1s);
4059  return src1*src2;
4060 }
4061 
4062 
4063 
4073 template <int dim, typename Number>
4074 inline
4075 std::ostream &operator << (std::ostream &out,
4077 {
4078  //make out lives a bit simpler by outputing
4079  //the tensor through the operator for the
4080  //general Tensor class
4082 
4083  for (unsigned int i=0; i<dim; ++i)
4084  for (unsigned int j=0; j<dim; ++j)
4085  tt[i][j] = t[i][j];
4086 
4087  return out << tt;
4088 }
4089 
4090 
4091 
4101 template <int dim, typename Number>
4102 inline
4103 std::ostream &operator << (std::ostream &out,
4105 {
4106  //make out lives a bit simpler by outputing
4107  //the tensor through the operator for the
4108  //general Tensor class
4110 
4111  for (unsigned int i=0; i<dim; ++i)
4112  for (unsigned int j=0; j<dim; ++j)
4113  for (unsigned int k=0; k<dim; ++k)
4114  for (unsigned int l=0; l<dim; ++l)
4115  tt[i][j][k][l] = t[i][j][k][l];
4116 
4117  return out << tt;
4118 }
4119 
4120 
4121 DEAL_II_NAMESPACE_CLOSE
4122 
4123 #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()