Reference documentation for deal.II version Git e8a02dd 2017-08-23 11:09:52 +0200
symmetric_tensor.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2017 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 
552  static const unsigned int n_independent_components
555 
559  SymmetricTensor ();
560 
571  template <typename OtherNumber>
572  explicit
574 
590  SymmetricTensor (const Number (&array) [n_independent_components]);
591 
597  template <typename OtherNumber>
598  explicit
600 
606  template <typename OtherNumber>
608 
615  SymmetricTensor &operator = (const Number &d);
616 
621  operator Tensor<rank,dim,Number> () const;
622 
626  bool operator == (const SymmetricTensor &) const;
627 
631  bool operator != (const SymmetricTensor &) const;
632 
636  template <typename OtherNumber>
638 
642  template <typename OtherNumber>
644 
649  template <typename OtherNumber>
650  SymmetricTensor &operator *= (const OtherNumber &factor);
651 
655  template <typename OtherNumber>
656  SymmetricTensor &operator /= (const OtherNumber &factor);
657 
662 
687  template <typename OtherNumber>
690 
695  template <typename OtherNumber>
698 
702  Number &operator() (const TableIndices<rank> &indices);
703 
707  Number operator() (const TableIndices<rank> &indices) const;
708 
713  internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
714  operator [] (const unsigned int row) const;
715 
720  internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
721  operator [] (const unsigned int row);
722 
728  Number
729  operator [] (const TableIndices<rank> &indices) const;
730 
736  Number &
737  operator [] (const TableIndices<rank> &indices);
738 
744  Number
745  access_raw_entry (const unsigned int unrolled_index) const;
746 
752  Number &
753  access_raw_entry (const unsigned int unrolled_index);
754 
765  norm () const;
766 
774  static
775  unsigned int
777 
783  static
785  unrolled_to_component_indices (const unsigned int i);
786 
799  void clear ();
800 
805  static std::size_t memory_consumption ();
806 
811  template <class Archive>
812  void serialize(Archive &ar, const unsigned int version);
813 
814 private:
818  typedef
821 
825  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
826 
830  base_tensor_type data;
831 
835  template <int, int, typename> friend class SymmetricTensor;
836 
840  template <int dim2, typename Number2>
841  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
842 
843  template <int dim2, typename Number2>
844  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
845 
846  template <int dim2, typename Number2>
848  deviator (const SymmetricTensor<2,dim2,Number2> &t);
849 
850  template <int dim2, typename Number2>
852 
853  template <int dim2, typename Number2>
855 
856  template <int dim2, typename Number2>
858 
859 
863  friend struct internal::SymmetricTensor::Inverse<2,dim,Number>;
864 
865  friend struct internal::SymmetricTensor::Inverse<4,dim,Number>;
866 };
867 
868 
869 
870 // ------------------------- inline functions ------------------------
871 
872 #ifndef DOXYGEN
873 
874 namespace internal
875 {
876  namespace SymmetricTensorAccessors
877  {
878  template <int rank, int dim, bool constness, int P, typename Number>
879  Accessor<rank,dim,constness,P,Number>::
880  Accessor (tensor_type &tensor,
881  const TableIndices<rank> &previous_indices)
882  :
883  tensor (tensor),
884  previous_indices (previous_indices)
885  {}
886 
887 
888 
889  template <int rank, int dim, bool constness, int P, typename Number>
890  Accessor<rank,dim,constness,P-1,Number>
891  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
892  {
893  return Accessor<rank,dim,constness,P-1,Number> (tensor,
894  merge (previous_indices, i, rank-P));
895  }
896 
897 
898 
899  template <int rank, int dim, bool constness, int P, typename Number>
900  Accessor<rank,dim,constness,P-1,Number>
901  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i) const
902  {
903  return Accessor<rank,dim,constness,P-1,Number> (tensor,
904  merge (previous_indices, i, rank-P));
905  }
906 
907 
908 
909  template <int rank, int dim, bool constness, typename Number>
910  Accessor<rank,dim,constness,1,Number>::
911  Accessor (tensor_type &tensor,
912  const TableIndices<rank> &previous_indices)
913  :
914  tensor (tensor),
915  previous_indices (previous_indices)
916  {}
917 
918 
919 
920  template <int rank, int dim, bool constness, typename Number>
921  typename Accessor<rank,dim,constness,1,Number>::reference
922  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
923  {
924  return tensor(merge (previous_indices, i, rank-1));
925  }
926 
927 
928  template <int rank, int dim, bool constness, typename Number>
929  typename Accessor<rank,dim,constness,1,Number>::reference
930  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i) const
931  {
932  return tensor(merge (previous_indices, i, rank-1));
933  }
934  }
935 }
936 
937 
938 
939 template <int rank, int dim, typename Number>
940 inline
942 {
943  // Some auto-differentiable numbers need explicit
944  // zero initialization.
945  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
947 }
948 
949 
950 template <int rank, int dim, typename Number>
951 template <typename OtherNumber>
952 inline
954 {
955  Assert (rank == 2, ExcNotImplemented());
956  switch (dim)
957  {
958  case 2:
959  Assert (t[0][1] == t[1][0], ExcInternalError());
960 
961  data[0] = t[0][0];
962  data[1] = t[1][1];
963  data[2] = t[0][1];
964 
965  break;
966  case 3:
967  Assert (t[0][1] == t[1][0], ExcInternalError());
968  Assert (t[0][2] == t[2][0], ExcInternalError());
969  Assert (t[1][2] == t[2][1], ExcInternalError());
970 
971  data[0] = t[0][0];
972  data[1] = t[1][1];
973  data[2] = t[2][2];
974  data[3] = t[0][1];
975  data[4] = t[0][2];
976  data[5] = t[1][2];
977 
978  break;
979  default:
980  for (unsigned int d=0; d<dim; ++d)
981  for (unsigned int e=0; e<d; ++e)
982  Assert(t[d][e] == t[e][d], ExcInternalError());
983 
984  for (unsigned int d=0; d<dim; ++d)
985  data[d] = t[d][d];
986 
987  for (unsigned int d=0, c=0; d<dim; ++d)
988  for (unsigned int e=d+1; e<dim; ++e, ++c)
989  data[dim+c] = t[d][e];
990  }
991 }
992 
993 
994 
995 template <int rank, int dim, typename Number>
996 template <typename OtherNumber>
997 inline
1000 {
1001  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1002  data[i] = initializer.data[i];
1003 }
1004 
1005 
1006 
1007 
1008 template <int rank, int dim, typename Number>
1009 inline
1010 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
1011  :
1012  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1013 {
1014  // ensure that the reinterpret_cast above actually works
1015  Assert (sizeof(typename base_tensor_type::array_type)
1016  == sizeof(array),
1017  ExcInternalError());
1018 }
1019 
1020 
1021 
1022 template <int rank, int dim, typename Number>
1023 template <typename OtherNumber>
1024 inline
1027 {
1028  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1029  data[i] = t.data[i];
1030  return *this;
1031 }
1032 
1033 
1034 
1035 template <int rank, int dim, typename Number>
1036 inline
1039 {
1040  Assert (d==internal::NumberType<Number>::value(0.0), ExcMessage ("Only assignment with zero is allowed"));
1041  (void) d;
1042 
1044 
1045  return *this;
1046 }
1047 
1048 
1049 namespace internal
1050 {
1051  namespace SymmetricTensor
1052  {
1053  template <int dim, typename Number>
1055  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1056  {
1058 
1059  // diagonal entries are stored first
1060  for (unsigned int d=0; d<dim; ++d)
1061  t[d][d] = s.access_raw_entry(d);
1062 
1063  // off-diagonal entries come next, row by row
1064  for (unsigned int d=0, c=0; d<dim; ++d)
1065  for (unsigned int e=d+1; e<dim; ++e, ++c)
1066  {
1067  t[d][e] = s.access_raw_entry(dim+c);
1068  t[e][d] = s.access_raw_entry(dim+c);
1069  }
1070  return t;
1071  }
1072 
1073 
1074  template <int dim, typename Number>
1076  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1077  {
1078  // utilize the symmetry properties of SymmetricTensor<4,dim>
1079  // discussed in the class documentation to avoid accessing all
1080  // independent elements of the input tensor more than once
1082 
1083  for (unsigned int i=0; i<dim; ++i)
1084  for (unsigned int j=i; j<dim; ++j)
1085  for (unsigned int k=0; k<dim; ++k)
1086  for (unsigned int l=k; l<dim; ++l)
1087  t[TableIndices<4>(i,j,k,l)]
1088  = t[TableIndices<4>(i,j,l,k)]
1089  = t[TableIndices<4>(j,i,k,l)]
1090  = t[TableIndices<4>(j,i,l,k)]
1091  = st[TableIndices<4>(i,j,k,l)];
1092 
1093  return t;
1094  }
1095 
1096 
1097  template <typename Number>
1098  struct Inverse<2,1,Number>
1099  {
1100  static inline ::SymmetricTensor<2,1,Number>
1101  value (const ::SymmetricTensor<2,1,Number> &t)
1102  {
1104 
1105  tmp[0][0] = 1.0/t[0][0];
1106 
1107  return tmp;
1108  }
1109  };
1110 
1111 
1112  template <typename Number>
1113  struct Inverse<2,2,Number>
1114  {
1115  static inline ::SymmetricTensor<2,2,Number>
1116  value (const ::SymmetricTensor<2,2,Number> &t)
1117  {
1119 
1120  // Sympy result: ([
1121  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1122  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1123  const TableIndices<2> idx_00 (0,0);
1124  const TableIndices<2> idx_01 (0,1);
1125  const TableIndices<2> idx_11 (1,1);
1126  const Number inv_det_t
1127  = 1.0/(t[idx_00]*t[idx_11]
1128  - t[idx_01]*t[idx_01]);
1129  tmp[idx_00] = t[idx_11];
1130  tmp[idx_01] = -t[idx_01];
1131  tmp[idx_11] = t[idx_00];
1132  tmp *= inv_det_t;
1133 
1134  return tmp;
1135  }
1136  };
1137 
1138 
1139  template <typename Number>
1140  struct Inverse<2,3,Number>
1141  {
1142  static ::SymmetricTensor<2,3,Number>
1143  value (const ::SymmetricTensor<2,3,Number> &t)
1144  {
1146 
1147  // Sympy result: ([
1148  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1149  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1150  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
1151  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1152  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1153  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
1154  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1155  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
1156  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
1157  const TableIndices<2> idx_00 (0,0);
1158  const TableIndices<2> idx_01 (0,1);
1159  const TableIndices<2> idx_02 (0,2);
1160  const TableIndices<2> idx_11 (1,1);
1161  const TableIndices<2> idx_12 (1,2);
1162  const TableIndices<2> idx_22 (2,2);
1163  const Number inv_det_t
1164  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
1165  - t[idx_00]*t[idx_12]*t[idx_12]
1166  - t[idx_01]*t[idx_01]*t[idx_22]
1167  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
1168  - t[idx_02]*t[idx_02]*t[idx_11]);
1169  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
1170  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
1171  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
1172  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
1173  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
1174  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
1175  tmp *= inv_det_t;
1176 
1177  return tmp;
1178  }
1179  };
1180 
1181 
1182  template <typename Number>
1183  struct Inverse<4,1,Number>
1184  {
1185  static inline ::SymmetricTensor<4,1,Number>
1186  value (const ::SymmetricTensor<4,1,Number> &t)
1187  {
1189  tmp.data[0][0] = 1.0/t.data[0][0];
1190  return tmp;
1191  }
1192  };
1193 
1194 
1195  template <typename Number>
1196  struct Inverse<4,2,Number>
1197  {
1198  static inline ::SymmetricTensor<4,2,Number>
1199  value (const ::SymmetricTensor<4,2,Number> &t)
1200  {
1202 
1203  // Inverting this tensor is a little more complicated than necessary,
1204  // since we store the data of 't' as a 3x3 matrix t.data, but the
1205  // product between a rank-4 and a rank-2 tensor is really not the
1206  // product between this matrix and the 3-vector of a rhs, but rather
1207  //
1208  // B.vec = t.data * mult * A.vec
1209  //
1210  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1211  // capture the fact that we need to add up both the c_ij12*a_12 and the
1212  // c_ij21*a_21 terms.
1213  //
1214  // In addition, in this scheme, the identity tensor has the matrix
1215  // representation mult^-1.
1216  //
1217  // The inverse of 't' therefore has the matrix representation
1218  //
1219  // inv.data = mult^-1 * t.data^-1 * mult^-1
1220  //
1221  // in order to compute it, let's first compute the inverse of t.data and
1222  // put it into tmp.data; at the end of the function we then scale the
1223  // last row and column of the inverse by 1/2, corresponding to the left
1224  // and right multiplication with mult^-1.
1225  const Number t4 = t.data[0][0]*t.data[1][1],
1226  t6 = t.data[0][0]*t.data[1][2],
1227  t8 = t.data[0][1]*t.data[1][0],
1228  t00 = t.data[0][2]*t.data[1][0],
1229  t01 = t.data[0][1]*t.data[2][0],
1230  t04 = t.data[0][2]*t.data[2][0],
1231  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
1232  t8*t.data[2][2]+t00*t.data[2][1]+
1233  t01*t.data[1][2]-t04*t.data[1][1]);
1234  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
1235  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
1236  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
1237  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
1238  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
1239  tmp.data[1][2] = -(t6-t00)*t07;
1240  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
1241  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
1242  tmp.data[2][2] = (t4-t8)*t07;
1243 
1244  // scale last row and column as mentioned
1245  // above
1246  tmp.data[2][0] /= 2;
1247  tmp.data[2][1] /= 2;
1248  tmp.data[0][2] /= 2;
1249  tmp.data[1][2] /= 2;
1250  tmp.data[2][2] /= 4;
1251 
1252  return tmp;
1253  }
1254  };
1255 
1256 
1257  template <typename Number>
1258  struct Inverse<4,3,Number>
1259  {
1260  static ::SymmetricTensor<4,3,Number>
1261  value (const ::SymmetricTensor<4,3,Number> &t)
1262  {
1264 
1265  // This function follows the exact same scheme as the 2d case, except
1266  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1267  // Instead, we use the Gauss-Jordan algorithm implemented for
1268  // FullMatrix. For historical reasons the following code is copied from
1269  // there, with the tangential benefit that we do not need to copy the
1270  // tensor entries to and from the FullMatrix.
1271  const unsigned int N = 6;
1272 
1273  // First get an estimate of the size of the elements of this matrix,
1274  // for later checks whether the pivot element is large enough, or
1275  // whether we have to fear that the matrix is not regular.
1276  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1277  for (unsigned int i=0; i<N; ++i)
1278  diagonal_sum += std::fabs(tmp.data[i][i]);
1279  const Number typical_diagonal_element = diagonal_sum/N;
1280  (void)typical_diagonal_element;
1281 
1282  unsigned int p[N];
1283  for (unsigned int i=0; i<N; ++i)
1284  p[i] = i;
1285 
1286  for (unsigned int j=0; j<N; ++j)
1287  {
1288  // Pivot search: search that part of the line on and right of the
1289  // diagonal for the largest element.
1290  Number max = std::fabs(tmp.data[j][j]);
1291  unsigned int r = j;
1292  for (unsigned int i=j+1; i<N; ++i)
1293  if (std::fabs(tmp.data[i][j]) > max)
1294  {
1295  max = std::fabs(tmp.data[i][j]);
1296  r = i;
1297  }
1298 
1299  // Check whether the pivot is too small
1300  Assert(max > 1.e-16*typical_diagonal_element,
1301  ExcMessage("This tensor seems to be noninvertible"));
1302 
1303  // Row interchange
1304  if (r>j)
1305  {
1306  for (unsigned int k=0; k<N; ++k)
1307  std::swap (tmp.data[j][k], tmp.data[r][k]);
1308 
1309  std::swap (p[j], p[r]);
1310  }
1311 
1312  // Transformation
1313  const Number hr = 1./tmp.data[j][j];
1314  tmp.data[j][j] = hr;
1315  for (unsigned int k=0; k<N; ++k)
1316  {
1317  if (k==j) continue;
1318  for (unsigned int i=0; i<N; ++i)
1319  {
1320  if (i==j) continue;
1321  tmp.data[i][k] -= tmp.data[i][j]*tmp.data[j][k]*hr;
1322  }
1323  }
1324  for (unsigned int i=0; i<N; ++i)
1325  {
1326  tmp.data[i][j] *= hr;
1327  tmp.data[j][i] *= -hr;
1328  }
1329  tmp.data[j][j] = hr;
1330  }
1331 
1332  // Column interchange
1333  Number hv[N];
1334  for (unsigned int i=0; i<N; ++i)
1335  {
1336  for (unsigned int k=0; k<N; ++k)
1337  hv[p[k]] = tmp.data[i][k];
1338  for (unsigned int k=0; k<N; ++k)
1339  tmp.data[i][k] = hv[k];
1340  }
1341 
1342  // Scale rows and columns. The mult matrix
1343  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1344  for (unsigned int i=3; i<6; ++i)
1345  for (unsigned int j=0; j<3; ++j)
1346  tmp.data[i][j] /= 2;
1347 
1348  for (unsigned int i=0; i<3; ++i)
1349  for (unsigned int j=3; j<6; ++j)
1350  tmp.data[i][j] /= 2;
1351 
1352  for (unsigned int i=3; i<6; ++i)
1353  for (unsigned int j=3; j<6; ++j)
1354  tmp.data[i][j] /= 4;
1355 
1356  return tmp;
1357  }
1358  };
1359 
1360  }
1361 }
1362 
1363 
1364 
1365 template <int rank, int dim, typename Number>
1366 inline
1368 operator Tensor<rank,dim,Number> () const
1369 {
1370  return internal::SymmetricTensor::convert_to_tensor (*this);
1371 }
1372 
1373 
1374 
1375 template <int rank, int dim, typename Number>
1376 inline
1377 bool
1379 (const SymmetricTensor<rank,dim,Number> &t) const
1380 {
1381  return data == t.data;
1382 }
1383 
1384 
1385 
1386 template <int rank, int dim, typename Number>
1387 inline
1388 bool
1389 SymmetricTensor<rank,dim,Number>::operator !=
1390 (const SymmetricTensor<rank,dim,Number> &t) const
1391 {
1392  return data != t.data;
1393 }
1394 
1395 
1396 
1397 template <int rank, int dim, typename Number>
1398 template <typename OtherNumber>
1399 inline
1401 SymmetricTensor<rank,dim,Number>::operator +=
1403 {
1404  data += t.data;
1405  return *this;
1406 }
1407 
1408 
1409 
1410 template <int rank, int dim, typename Number>
1411 template <typename OtherNumber>
1412 inline
1414 SymmetricTensor<rank,dim,Number>::operator -=
1416 {
1417  data -= t.data;
1418  return *this;
1419 }
1420 
1421 
1422 
1423 template <int rank, int dim, typename Number>
1424 template <typename OtherNumber>
1425 inline
1428 {
1429  data *= d;
1430  return *this;
1431 }
1432 
1433 
1434 
1435 template <int rank, int dim, typename Number>
1436 template <typename OtherNumber>
1437 inline
1440 {
1441  data /= d;
1442  return *this;
1443 }
1444 
1445 
1446 
1447 template <int rank, int dim, typename Number>
1448 inline
1451 {
1452  SymmetricTensor tmp = *this;
1453  tmp.data = -tmp.data;
1454  return tmp;
1455 }
1456 
1457 
1458 
1459 template <int rank, int dim, typename Number>
1460 inline
1461 void
1463 {
1464  data.clear ();
1465 }
1466 
1467 
1468 
1469 template <int rank, int dim, typename Number>
1470 inline
1471 std::size_t
1473 {
1474  // all memory consists of statically allocated memory of the current
1475  // object, no pointers
1476  return sizeof(SymmetricTensor<rank,dim,Number>);
1477 }
1478 
1479 
1480 
1481 namespace internal
1482 {
1483 
1484  template <int dim, typename Number, typename OtherNumber = Number>
1485  inline
1486  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type
1487  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1488  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1489  {
1490  typedef typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type result_type;
1491 
1492  switch (dim)
1493  {
1494  case 1:
1495  return data[0] * sdata[0];
1496  default:
1497  // Start with the non-diagonal part to avoid some multiplications by
1498  // 2.
1499 
1500  result_type sum = data[dim] * sdata[dim];
1501  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1502  sum += data[d] * sdata[d];
1503  sum += sum; // sum = sum * 2.;
1504 
1505  // Now add the contributions from the diagonal
1506  for (unsigned int d=0; d<dim; ++d)
1507  sum += data[d] * sdata[d];
1508  return sum;
1509  }
1510  }
1511 
1512 
1513 
1514  template <int dim, typename Number, typename OtherNumber = Number>
1515  inline
1516  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type
1517  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,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<4,2,dim,Number,OtherNumber>::type result_type;
1521  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::value_type value_type;
1522 
1523  const unsigned int data_dim =
1524  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1525  value_type tmp [data_dim];
1526  for (unsigned int i=0; i<data_dim; ++i)
1527  tmp[i] = perform_double_contraction<dim,Number,OtherNumber>(data[i], sdata);
1528  return result_type(tmp);
1529  }
1530 
1531 
1532 
1533  template <int dim, typename Number, typename OtherNumber = Number>
1534  inline
1535  typename SymmetricTensorAccessors::StorageType<2,dim,
1536  typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type
1537  >::base_tensor_type
1538  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1539  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1540  {
1541  typedef typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::type result_type;
1542  typedef typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type value_type;
1543  typedef typename SymmetricTensorAccessors::StorageType<2,dim,value_type>::base_tensor_type base_tensor_type;
1544 
1545  base_tensor_type tmp;
1546  for (unsigned int i=0; i<tmp.dimension; ++i)
1547  {
1548  // Start with the non-diagonal part
1549  value_type sum = data[dim] * sdata[dim][i];
1550  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1551  sum += data[d] * sdata[d][i];
1552  sum += sum; // sum = sum * 2.;
1553 
1554  // Now add the contributions from the diagonal
1555  for (unsigned int d=0; d<dim; ++d)
1556  sum += data[d] * sdata[d][i];
1557  tmp[i] = sum;
1558  }
1559  return tmp;
1560  }
1561 
1562 
1563 
1564  template <int dim, typename Number, typename OtherNumber = Number>
1565  inline
1566  typename SymmetricTensorAccessors::StorageType<4,dim,
1567  typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type
1568  >::base_tensor_type
1569  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1570  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1571  {
1572  typedef typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::type result_type;
1573  typedef typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type value_type;
1574  typedef typename SymmetricTensorAccessors::StorageType<4,dim,value_type>::base_tensor_type base_tensor_type;
1575 
1576  const unsigned int data_dim =
1577  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1578  base_tensor_type tmp;
1579  for (unsigned int i=0; i<data_dim; ++i)
1580  for (unsigned int j=0; j<data_dim; ++j)
1581  {
1582  // Start with the non-diagonal part
1583  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1584  tmp[i][j] += data[i][d] * sdata[d][j];
1585  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1586 
1587  // Now add the contributions from the diagonal
1588  for (unsigned int d=0; d<dim; ++d)
1589  tmp[i][j] += data[i][d] * sdata[d][j];
1590  }
1591  return tmp;
1592  }
1593 
1594 } // end of namespace internal
1595 
1596 
1597 
1598 template <int rank, int dim, typename Number>
1599 template <typename OtherNumber>
1600 inline
1603 {
1604  // need to have two different function calls
1605  // because a scalar and rank-2 tensor are not
1606  // the same data type (see internal function
1607  // above)
1608  return internal::perform_double_contraction<dim,Number,OtherNumber> (data, s.data);
1609 }
1610 
1611 
1612 
1613 template <int rank, int dim, typename Number>
1614 template <typename OtherNumber>
1615 inline
1618 {
1621  tmp.data = internal::perform_double_contraction<dim,Number,OtherNumber> (data,s.data);
1622  return tmp;
1623 }
1624 
1625 
1626 
1627 // internal namespace to switch between the
1628 // access of different tensors. There used to
1629 // be explicit instantiations before for
1630 // different ranks and dimensions, but since
1631 // we now allow for templates on the data
1632 // type, and since we cannot partially
1633 // specialize the implementation, this got
1634 // into a separate namespace
1635 namespace internal
1636 {
1637  template <int dim, typename Number>
1638  inline
1639  Number &
1640  symmetric_tensor_access (const TableIndices<2> &indices,
1642  {
1643  // 1d is very simple and done first
1644  if (dim == 1)
1645  return data[0];
1646 
1647  // first treat the main diagonal elements, which are stored consecutively
1648  // at the beginning
1649  if (indices[0] == indices[1])
1650  return data[indices[0]];
1651 
1652  // the rest is messier and requires a few switches.
1653  switch (dim)
1654  {
1655  case 2:
1656  // at least for the 2x2 case it is reasonably simple
1657  Assert (((indices[0]==1) && (indices[1]==0)) ||
1658  ((indices[0]==0) && (indices[1]==1)),
1659  ExcInternalError());
1660  return data[2];
1661 
1662  default:
1663  // to do the rest, sort our indices before comparing
1664  {
1665  TableIndices<2> sorted_indices (indices);
1666  sorted_indices.sort ();
1667 
1668  for (unsigned int d=0, c=0; d<dim; ++d)
1669  for (unsigned int e=d+1; e<dim; ++e, ++c)
1670  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1671  return data[dim+c];
1672  Assert (false, ExcInternalError());
1673  }
1674  }
1675 
1676  static Number dummy_but_referenceable = Number();
1677  return dummy_but_referenceable;
1678  }
1679 
1680 
1681 
1682  template <int dim, typename Number>
1683  inline
1684  Number
1685  symmetric_tensor_access (const TableIndices<2> &indices,
1687  {
1688  // 1d is very simple and done first
1689  if (dim == 1)
1690  return data[0];
1691 
1692  // first treat the main diagonal elements, which are stored consecutively
1693  // at the beginning
1694  if (indices[0] == indices[1])
1695  return data[indices[0]];
1696 
1697  // the rest is messier and requires a few switches.
1698  switch (dim)
1699  {
1700  case 2:
1701  // at least for the 2x2 case it is reasonably simple
1702  Assert (((indices[0]==1) && (indices[1]==0)) ||
1703  ((indices[0]==0) && (indices[1]==1)),
1704  ExcInternalError());
1705  return data[2];
1706 
1707  default:
1708  // to do the rest, sort our indices before comparing
1709  {
1710  TableIndices<2> sorted_indices (indices);
1711  sorted_indices.sort ();
1712 
1713  for (unsigned int d=0, c=0; d<dim; ++d)
1714  for (unsigned int e=d+1; e<dim; ++e, ++c)
1715  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1716  return data[dim+c];
1717  Assert (false, ExcInternalError());
1718  }
1719  }
1720 
1721  static Number dummy_but_referenceable = Number();
1722  return dummy_but_referenceable;
1723  }
1724 
1725 
1726 
1727  template <int dim, typename Number>
1728  inline
1729  Number &
1730  symmetric_tensor_access (const TableIndices<4> &indices,
1732  {
1733  switch (dim)
1734  {
1735  case 1:
1736  return data[0][0];
1737 
1738  case 2:
1739  // each entry of the tensor can be
1740  // thought of as an entry in a
1741  // matrix that maps the rolled-out
1742  // rank-2 tensors into rolled-out
1743  // rank-2 tensors. this is the
1744  // format in which we store rank-4
1745  // tensors. determine which
1746  // position the present entry is
1747  // stored in
1748  {
1749  unsigned int base_index[2] ;
1750  if ((indices[0] == 0) && (indices[1] == 0))
1751  base_index[0] = 0;
1752  else if ((indices[0] == 1) && (indices[1] == 1))
1753  base_index[0] = 1;
1754  else
1755  base_index[0] = 2;
1756 
1757  if ((indices[2] == 0) && (indices[3] == 0))
1758  base_index[1] = 0;
1759  else if ((indices[2] == 1) && (indices[3] == 1))
1760  base_index[1] = 1;
1761  else
1762  base_index[1] = 2;
1763 
1764  return data[base_index[0]][base_index[1]];
1765  }
1766 
1767  case 3:
1768  // each entry of the tensor can be
1769  // thought of as an entry in a
1770  // matrix that maps the rolled-out
1771  // rank-2 tensors into rolled-out
1772  // rank-2 tensors. this is the
1773  // format in which we store rank-4
1774  // tensors. determine which
1775  // position the present entry is
1776  // stored in
1777  {
1778  unsigned int base_index[2] ;
1779  if ((indices[0] == 0) && (indices[1] == 0))
1780  base_index[0] = 0;
1781  else if ((indices[0] == 1) && (indices[1] == 1))
1782  base_index[0] = 1;
1783  else if ((indices[0] == 2) && (indices[1] == 2))
1784  base_index[0] = 2;
1785  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1786  ((indices[0] == 1) && (indices[1] == 0)))
1787  base_index[0] = 3;
1788  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1789  ((indices[0] == 2) && (indices[1] == 0)))
1790  base_index[0] = 4;
1791  else
1792  {
1793  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1794  ((indices[0] == 2) && (indices[1] == 1)),
1795  ExcInternalError());
1796  base_index[0] = 5;
1797  }
1798 
1799  if ((indices[2] == 0) && (indices[3] == 0))
1800  base_index[1] = 0;
1801  else if ((indices[2] == 1) && (indices[3] == 1))
1802  base_index[1] = 1;
1803  else if ((indices[2] == 2) && (indices[3] == 2))
1804  base_index[1] = 2;
1805  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1806  ((indices[2] == 1) && (indices[3] == 0)))
1807  base_index[1] = 3;
1808  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1809  ((indices[2] == 2) && (indices[3] == 0)))
1810  base_index[1] = 4;
1811  else
1812  {
1813  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1814  ((indices[2] == 2) && (indices[3] == 1)),
1815  ExcInternalError());
1816  base_index[1] = 5;
1817  }
1818 
1819  return data[base_index[0]][base_index[1]];
1820  }
1821 
1822  default:
1823  Assert (false, ExcNotImplemented());
1824  }
1825 
1826  static Number dummy;
1827  return dummy;
1828  }
1829 
1830 
1831  template <int dim, typename Number>
1832  inline
1833  Number
1834  symmetric_tensor_access (const TableIndices<4> &indices,
1836  {
1837  switch (dim)
1838  {
1839  case 1:
1840  return data[0][0];
1841 
1842  case 2:
1843  // each entry of the tensor can be
1844  // thought of as an entry in a
1845  // matrix that maps the rolled-out
1846  // rank-2 tensors into rolled-out
1847  // rank-2 tensors. this is the
1848  // format in which we store rank-4
1849  // tensors. determine which
1850  // position the present entry is
1851  // stored in
1852  {
1853  unsigned int base_index[2] ;
1854  if ((indices[0] == 0) && (indices[1] == 0))
1855  base_index[0] = 0;
1856  else if ((indices[0] == 1) && (indices[1] == 1))
1857  base_index[0] = 1;
1858  else
1859  base_index[0] = 2;
1860 
1861  if ((indices[2] == 0) && (indices[3] == 0))
1862  base_index[1] = 0;
1863  else if ((indices[2] == 1) && (indices[3] == 1))
1864  base_index[1] = 1;
1865  else
1866  base_index[1] = 2;
1867 
1868  return data[base_index[0]][base_index[1]];
1869  }
1870 
1871  case 3:
1872  // each entry of the tensor can be
1873  // thought of as an entry in a
1874  // matrix that maps the rolled-out
1875  // rank-2 tensors into rolled-out
1876  // rank-2 tensors. this is the
1877  // format in which we store rank-4
1878  // tensors. determine which
1879  // position the present entry is
1880  // stored in
1881  {
1882  unsigned int base_index[2] ;
1883  if ((indices[0] == 0) && (indices[1] == 0))
1884  base_index[0] = 0;
1885  else if ((indices[0] == 1) && (indices[1] == 1))
1886  base_index[0] = 1;
1887  else if ((indices[0] == 2) && (indices[1] == 2))
1888  base_index[0] = 2;
1889  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1890  ((indices[0] == 1) && (indices[1] == 0)))
1891  base_index[0] = 3;
1892  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1893  ((indices[0] == 2) && (indices[1] == 0)))
1894  base_index[0] = 4;
1895  else
1896  {
1897  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1898  ((indices[0] == 2) && (indices[1] == 1)),
1899  ExcInternalError());
1900  base_index[0] = 5;
1901  }
1902 
1903  if ((indices[2] == 0) && (indices[3] == 0))
1904  base_index[1] = 0;
1905  else if ((indices[2] == 1) && (indices[3] == 1))
1906  base_index[1] = 1;
1907  else if ((indices[2] == 2) && (indices[3] == 2))
1908  base_index[1] = 2;
1909  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1910  ((indices[2] == 1) && (indices[3] == 0)))
1911  base_index[1] = 3;
1912  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1913  ((indices[2] == 2) && (indices[3] == 0)))
1914  base_index[1] = 4;
1915  else
1916  {
1917  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1918  ((indices[2] == 2) && (indices[3] == 1)),
1919  ExcInternalError());
1920  base_index[1] = 5;
1921  }
1922 
1923  return data[base_index[0]][base_index[1]];
1924  }
1925 
1926  default:
1927  Assert (false, ExcNotImplemented());
1928  }
1929 
1930  static Number dummy;
1931  return dummy;
1932  }
1933 
1934 } // end of namespace internal
1935 
1936 
1937 
1938 template <int rank, int dim, typename Number>
1939 inline
1940 Number &
1942 {
1943  for (unsigned int r=0; r<rank; ++r)
1944  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1945  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1946 }
1947 
1948 
1949 
1950 template <int rank, int dim, typename Number>
1951 inline
1952 Number
1954 (const TableIndices<rank> &indices) const
1955 {
1956  for (unsigned int r=0; r<rank; ++r)
1957  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1958  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1959 }
1960 
1961 
1962 
1963 namespace internal
1964 {
1965  namespace SymmetricTensor
1966  {
1967  template <int rank>
1969  get_partially_filled_indices (const unsigned int row,
1970  const std::integral_constant<int, 2> &)
1971  {
1972  return TableIndices<rank> (row,
1974 
1975  }
1976 
1977 
1978  template <int rank>
1980  get_partially_filled_indices (const unsigned int row,
1981  const std::integral_constant<int, 4> &)
1982  {
1983  return TableIndices<rank> (row,
1987 
1988  }
1989  }
1990 }
1991 
1992 
1993 template <int rank, int dim, typename Number>
1994 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1995 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1996 {
1997  return
1998  internal::SymmetricTensorAccessors::
1999  Accessor<rank,dim,true,rank-1,Number> (*this,
2000  internal::SymmetricTensor::get_partially_filled_indices<rank> (row,
2001  std::integral_constant<int, rank>()));
2002 }
2003 
2004 
2005 
2006 template <int rank, int dim, typename Number>
2007 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
2008 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
2009 {
2010  return
2011  internal::SymmetricTensorAccessors::
2012  Accessor<rank,dim,false,rank-1,Number> (*this,
2013  internal::SymmetricTensor::get_partially_filled_indices<rank> (row,
2014  std::integral_constant<int, rank>()));
2015 }
2016 
2017 
2018 
2019 template <int rank, int dim, typename Number>
2020 inline
2021 Number
2023 {
2024  return operator()(indices);
2025 }
2026 
2027 
2028 
2029 template <int rank, int dim, typename Number>
2030 inline
2031 Number &
2033 {
2034  return operator()(indices);
2035 }
2036 
2037 
2038 
2039 
2040 namespace internal
2041 {
2042  namespace SymmetricTensor
2043  {
2044  template <int dim, typename Number>
2045  unsigned int
2046  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
2047  const unsigned int index)
2048  {
2049  return index;
2050  }
2051 
2052 
2053  template <int dim, typename Number>
2055  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
2056  const unsigned int index)
2057  {
2058  return
2061  }
2062 
2063  }
2064 }
2065 
2066 
2067 
2068 template <int rank, int dim, typename Number>
2069 inline
2070 Number
2071 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
2072 {
2073  AssertIndexRange (index, n_independent_components);
2074  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2075 }
2076 
2077 
2078 
2079 template <int rank, int dim, typename Number>
2080 inline
2081 Number &
2082 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
2083 {
2084  AssertIndexRange (index, n_independent_components);
2085  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2086 }
2087 
2088 
2089 
2090 namespace internal
2091 {
2092  template <int dim, typename Number>
2093  inline
2096  {
2097  switch (dim)
2098  {
2099  case 1:
2100  return numbers::NumberTraits<Number>::abs(data[0]);
2101 
2102  case 2:
2103  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2106 
2107  case 3:
2108  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2114 
2115  default:
2116  {
2117  typename numbers::NumberTraits<Number>::real_type return_value
2119 
2120  for (unsigned int d=0; d<dim; ++d)
2121  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
2122  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
2123  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2124 
2125  return std::sqrt(return_value);
2126  }
2127  }
2128  }
2129 
2130 
2131 
2132  template <int dim, typename Number>
2133  inline
2136  {
2137  switch (dim)
2138  {
2139  case 1:
2140  return numbers::NumberTraits<Number>::abs (data[0][0]);
2141 
2142  default:
2143  {
2144  typename numbers::NumberTraits<Number>::real_type return_value
2146 
2147  const unsigned int n_independent_components = data.dimension;
2148 
2149  for (unsigned int i=0; i<dim; ++i)
2150  for (unsigned int j=0; j<dim; ++j)
2151  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
2152  for (unsigned int i=0; i<dim; ++i)
2153  for (unsigned int j=dim; j<n_independent_components; ++j)
2154  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2155  for (unsigned int i=dim; i<n_independent_components; ++i)
2156  for (unsigned int j=0; j<dim; ++j)
2157  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2158  for (unsigned int i=dim; i<n_independent_components; ++i)
2159  for (unsigned int j=dim; j<n_independent_components; ++j)
2160  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2161 
2162  return std::sqrt(return_value);
2163  }
2164  }
2165  }
2166 
2167 } // end of namespace internal
2168 
2169 
2170 
2171 template <int rank, int dim, typename Number>
2172 inline
2175 {
2176  return internal::compute_norm<dim,Number> (data);
2177 }
2178 
2179 
2180 
2181 namespace internal
2182 {
2183  namespace SymmetricTensor
2184  {
2185  namespace
2186  {
2187  // a function to do the unrolling from a set of indices to a
2188  // scalar index into the array in which we store the elements of
2189  // a symmetric tensor
2190  //
2191  // this function is for rank-2 tensors
2192  template <int dim>
2193  inline
2194  unsigned int
2195  component_to_unrolled_index
2196  (const TableIndices<2> &indices)
2197  {
2198  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2199  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2200 
2201  switch (dim)
2202  {
2203  case 1:
2204  {
2205  return 0;
2206  }
2207 
2208  case 2:
2209  {
2210  static const unsigned int table[2][2] = {{0, 2},
2211  {2, 1}
2212  };
2213  return table[indices[0]][indices[1]];
2214  }
2215 
2216  case 3:
2217  {
2218  static const unsigned int table[3][3] = {{0, 3, 4},
2219  {3, 1, 5},
2220  {4, 5, 2}
2221  };
2222  return table[indices[0]][indices[1]];
2223  }
2224 
2225  case 4:
2226  {
2227  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2228  {4, 1, 7, 8},
2229  {5, 7, 2, 9},
2230  {6, 8, 9, 3}
2231  };
2232  return table[indices[0]][indices[1]];
2233  }
2234 
2235  default:
2236  // for the remainder, manually figure out the numbering
2237  {
2238  if (indices[0] == indices[1])
2239  return indices[0];
2240 
2241  TableIndices<2> sorted_indices (indices);
2242  sorted_indices.sort ();
2243 
2244  for (unsigned int d=0, c=0; d<dim; ++d)
2245  for (unsigned int e=d+1; e<dim; ++e, ++c)
2246  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
2247  return dim+c;
2248 
2249  // should never get here:
2250  Assert(false, ExcInternalError());
2251  return 0;
2252  }
2253  }
2254  }
2255 
2256  // a function to do the unrolling from a set of indices to a
2257  // scalar index into the array in which we store the elements of
2258  // a symmetric tensor
2259  //
2260  // this function is for tensors of ranks not already handled
2261  // above
2262  template <int dim, int rank>
2263  inline
2264  unsigned int
2265  component_to_unrolled_index
2266  (const TableIndices<rank> &indices)
2267  {
2268  (void)indices;
2269  Assert (false, ExcNotImplemented());
2271  }
2272  }
2273  }
2274 }
2275 
2276 
2277 template <int rank, int dim, typename Number>
2278 inline
2279 unsigned int
2281 (const TableIndices<rank> &indices)
2282 {
2283  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
2284 }
2285 
2286 
2287 
2288 namespace internal
2289 {
2290  namespace SymmetricTensor
2291  {
2292  namespace
2293  {
2294  // a function to do the inverse of the unrolling from a set of
2295  // indices to a scalar index into the array in which we store
2296  // the elements of a symmetric tensor. in other words, it goes
2297  // from the scalar index into the array to a set of indices of
2298  // the tensor
2299  //
2300  // this function is for rank-2 tensors
2301  template <int dim>
2302  inline
2304  unrolled_to_component_indices
2305  (const unsigned int i,
2306  const std::integral_constant<int, 2> &)
2307  {
2310  switch (dim)
2311  {
2312  case 1:
2313  {
2314  return TableIndices<2>(0,0);
2315  }
2316 
2317  case 2:
2318  {
2319  const TableIndices<2> table[3] =
2320  {
2321  TableIndices<2> (0,0),
2322  TableIndices<2> (1,1),
2323  TableIndices<2> (0,1)
2324  };
2325  return table[i];
2326  }
2327 
2328  case 3:
2329  {
2330  const TableIndices<2> table[6] =
2331  {
2332  TableIndices<2> (0,0),
2333  TableIndices<2> (1,1),
2334  TableIndices<2> (2,2),
2335  TableIndices<2> (0,1),
2336  TableIndices<2> (0,2),
2337  TableIndices<2> (1,2)
2338  };
2339  return table[i];
2340  }
2341 
2342  default:
2343  if (i<dim)
2344  return TableIndices<2> (i,i);
2345 
2346  for (unsigned int d=0, c=0; d<dim; ++d)
2347  for (unsigned int e=d+1; e<dim; ++e, ++c)
2348  if (c==i)
2349  return TableIndices<2>(d,e);
2350 
2351  // should never get here:
2352  Assert(false, ExcInternalError());
2353  return TableIndices<2>(0, 0);
2354  }
2355  }
2356 
2357  // a function to do the inverse of the unrolling from a set of
2358  // indices to a scalar index into the array in which we store
2359  // the elements of a symmetric tensor. in other words, it goes
2360  // from the scalar index into the array to a set of indices of
2361  // the tensor
2362  //
2363  // this function is for tensors of a rank not already handled
2364  // above
2365  template <int dim, int rank>
2366  inline
2368  unrolled_to_component_indices
2369  (const unsigned int i,
2370  const std::integral_constant<int, rank> &)
2371  {
2372  (void)i;
2375  Assert (false, ExcNotImplemented());
2376  return TableIndices<rank>();
2377  }
2378 
2379  }
2380  }
2381 }
2382 
2383 template <int rank, int dim, typename Number>
2384 inline
2387 (const unsigned int i)
2388 {
2389  return
2390  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2391  std::integral_constant<int, rank>());
2392 }
2393 
2394 
2395 
2396 template <int rank, int dim, typename Number>
2397 template <class Archive>
2398 inline
2399 void
2400 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2401 {
2402  ar &data;
2403 }
2404 
2405 
2406 #endif // DOXYGEN
2407 
2408 /* ----------------- Non-member functions operating on tensors. ------------ */
2409 
2410 
2423 template <int rank, int dim, typename Number, typename OtherNumber>
2424 inline
2428 {
2430  tmp += right;
2431  return tmp;
2432 }
2433 
2434 
2447 template <int rank, int dim, typename Number, typename OtherNumber>
2448 inline
2452 {
2454  tmp -= right;
2455  return tmp;
2456 }
2457 
2458 
2466 template <int rank, int dim, typename Number, typename OtherNumber>
2467 inline
2470  const Tensor<rank, dim, OtherNumber> &right)
2471 {
2472  return Tensor<rank, dim, Number>(left) + right;
2473 }
2474 
2475 
2483 template <int rank, int dim, typename Number, typename OtherNumber>
2484 inline
2488 {
2489  return left + Tensor<rank, dim, OtherNumber>(right);
2490 }
2491 
2492 
2500 template <int rank, int dim, typename Number, typename OtherNumber>
2501 inline
2504  const Tensor<rank, dim, OtherNumber> &right)
2505 {
2506  return Tensor<rank, dim, Number>(left) - right;
2507 }
2508 
2509 
2517 template <int rank, int dim, typename Number, typename OtherNumber>
2518 inline
2522 {
2523  return left - Tensor<rank, dim, OtherNumber>(right);
2524 }
2525 
2526 
2527 
2541 template <int dim, typename Number>
2542 inline
2544 {
2545  switch (dim)
2546  {
2547  case 1:
2548  return t.data[0];
2549  case 2:
2550  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2551  case 3:
2552  {
2553  // in analogy to general tensors, but
2554  // there's something to be simplified for
2555  // the present case
2556  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2557  return ( tmp + tmp
2558  +t.data[0]*t.data[1]*t.data[2]
2559  -t.data[0]*t.data[5]*t.data[5]
2560  -t.data[1]*t.data[4]*t.data[4]
2561  -t.data[2]*t.data[3]*t.data[3]);
2562  }
2563  default:
2564  Assert (false, ExcNotImplemented());
2565  return 0;
2566  }
2567 }
2568 
2569 
2570 
2580 template <int dim, typename Number>
2581 inline
2583 {
2584  return determinant (t);
2585 }
2586 
2587 
2588 
2596 template <int dim, typename Number>
2598 {
2599  Number t = d.data[0];
2600  for (unsigned int i=1; i<dim; ++i)
2601  t += d.data[i];
2602  return t;
2603 }
2604 
2605 
2615 template <int dim, typename Number>
2616 inline
2618 {
2619  return trace (t);
2620 }
2621 
2622 
2635 template <typename Number>
2636 inline
2638 {
2639  return 0;
2640 }
2641 
2642 
2643 
2664 template <typename Number>
2665 inline
2667 {
2668  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2669 }
2670 
2671 
2672 
2682 template <typename Number>
2683 inline
2685 {
2686  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2687  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2688 }
2689 
2690 
2691 
2700 template <typename Number>
2701 std::array<Number,1>
2702 eigenvalues (const SymmetricTensor<2,1,Number> &T);
2703 
2704 
2705 
2727 template <typename Number>
2728 std::array<Number,2>
2729 eigenvalues (const SymmetricTensor<2,2,Number> &T);
2730 
2731 
2732 
2752 template <typename Number>
2753 std::array<Number,3>
2754 eigenvalues (const SymmetricTensor<2,3,Number> &T);
2755 
2756 
2757 
2758 namespace internal
2759 {
2760  namespace SymmetricTensor
2761  {
2798  template <int dim, typename Number>
2799  void
2800  tridiagonalize (const ::SymmetricTensor<2,dim,Number> &A,
2801  ::Tensor<2,dim,Number> &Q,
2802  std::array<Number,dim> &d,
2803  std::array<Number,dim-1> &e);
2804 
2805 
2806 
2845  template <int dim, typename Number>
2846  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2847  ql_implicit_shifts (const ::SymmetricTensor<2,dim,Number> &A);
2848 
2849 
2850 
2889  template <int dim, typename Number>
2890  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2892 
2893 
2894 
2910  template <typename Number>
2911  std::array<std::pair<Number, Tensor<1,2,Number> >,2>
2912  hybrid (const ::SymmetricTensor<2,2,Number> &A);
2913 
2914 
2915 
2948  template <typename Number>
2949  std::array<std::pair<Number, Tensor<1,3,Number> >,3>
2950  hybrid (const ::SymmetricTensor<2,3,Number> &A);
2951 
2952  namespace
2953  {
2954 
2959  template <int dim, typename Number>
2960  struct SortEigenValuesVectors
2961  {
2962  typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
2963  bool operator() (const EigValsVecs &lhs,
2964  const EigValsVecs &rhs)
2965  {
2966  return lhs.first > rhs.first;
2967  }
2968  };
2969 
2970  }
2971 
2972  } // namespace SymmetricTensor
2973 
2974 } // namespace internal
2975 
2976 
2977 
2978 // The line below is to ensure that doxygen puts the full description
2979 // of this global enumeration into the documentation
2980 // See https://stackoverflow.com/a/1717984
3008 {
3018  hybrid,
3028  ql_implicit_shifts,
3036  jacobi
3037 };
3038 
3039 
3040 
3047 template <typename Number>
3048 std::array<std::pair<Number, Tensor<1,1,Number> >,1>
3051 {
3052  return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} };
3053 }
3054 
3055 
3056 
3084 template <int dim, typename Number>
3085 std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
3088 {
3089  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
3090 
3091  switch (method)
3092  {
3094  eig_vals_vecs = internal::SymmetricTensor::hybrid(T);
3095  break;
3098  break;
3100  eig_vals_vecs = internal::SymmetricTensor::jacobi(T);
3101  break;
3102  default:
3103  AssertThrow(false, ExcNotImplemented());
3104  }
3105 
3106  // Sort in descending order before output.
3107  std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
3108  internal::SymmetricTensor::SortEigenValuesVectors<dim,Number>());
3109  return eig_vals_vecs;
3110 }
3111 
3112 
3113 
3123 template <int rank, int dim, typename Number>
3124 inline
3127 {
3128  return t;
3129 }
3130 
3131 
3132 
3142 template <int dim, typename Number>
3143 inline
3146 {
3148 
3149  // subtract scaled trace from the diagonal
3150  const Number tr = trace(t) / dim;
3151  for (unsigned int i=0; i<dim; ++i)
3152  tmp.data[i] -= tr;
3153 
3154  return tmp;
3155 }
3156 
3157 
3158 
3166 template <int dim, typename Number>
3167 inline
3169 unit_symmetric_tensor ()
3170 {
3171  // create a default constructed matrix filled with
3172  // zeros, then set the diagonal elements to one
3174  switch (dim)
3175  {
3176  case 1:
3177  tmp.data[0] = 1;
3178  break;
3179  case 2:
3180  tmp.data[0] = tmp.data[1] = 1;
3181  break;
3182  case 3:
3183  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3184  break;
3185  default:
3186  for (unsigned int d=0; d<dim; ++d)
3187  tmp.data[d] = 1;
3188  }
3189  return tmp;
3190 }
3191 
3192 
3193 
3202 template <int dim>
3203 inline
3205 unit_symmetric_tensor ()
3206 {
3207  return unit_symmetric_tensor<dim,double>();
3208 }
3209 
3210 
3211 
3226 template <int dim, typename Number>
3227 inline
3229 deviator_tensor ()
3230 {
3232 
3233  // fill the elements treating the diagonal
3234  for (unsigned int i=0; i<dim; ++i)
3235  for (unsigned int j=0; j<dim; ++j)
3236  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
3237 
3238  // then fill the ones that copy over the
3239  // non-diagonal elements. note that during
3240  // the double-contraction, we handle the
3241  // off-diagonal elements twice, so simply
3242  // copying requires a weight of 1/2
3243  for (unsigned int i=dim;
3244  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3245  ++i)
3246  tmp.data[i][i] = 0.5;
3247 
3248  return tmp;
3249 }
3250 
3251 
3252 
3267 template <int dim>
3268 inline
3270 deviator_tensor ()
3271 {
3272  return deviator_tensor<dim,double>();
3273 }
3274 
3275 
3276 
3299 template <int dim, typename Number>
3300 inline
3302 identity_tensor ()
3303 {
3305 
3306  // fill the elements treating the diagonal
3307  for (unsigned int i=0; i<dim; ++i)
3308  tmp.data[i][i] = 1;
3309 
3310  // then fill the ones that copy over the
3311  // non-diagonal elements. note that during
3312  // the double-contraction, we handle the
3313  // off-diagonal elements twice, so simply
3314  // copying requires a weight of 1/2
3315  for (unsigned int i=dim;
3316  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3317  ++i)
3318  tmp.data[i][i] = 0.5;
3319 
3320  return tmp;
3321 }
3322 
3323 
3324 
3346 template <int dim>
3347 inline
3349 identity_tensor ()
3350 {
3351  return identity_tensor<dim,double>();
3352 }
3353 
3354 
3355 
3366 template <int dim, typename Number>
3367 inline
3370 {
3372 }
3373 
3374 
3375 
3387 template <int dim, typename Number>
3388 inline
3391 {
3393 }
3394 
3395 
3396 
3411 template <int dim, typename Number>
3412 inline
3416 {
3418 
3419  // fill only the elements really needed
3420  for (unsigned int i=0; i<dim; ++i)
3421  for (unsigned int j=i; j<dim; ++j)
3422  for (unsigned int k=0; k<dim; ++k)
3423  for (unsigned int l=k; l<dim; ++l)
3424  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3425 
3426  return tmp;
3427 }
3428 
3429 
3430 
3439 template <int dim,typename Number>
3440 inline
3443 {
3444  Number array[(dim*dim+dim)/2];
3445  for (unsigned int d=0; d<dim; ++d)
3446  array[d] = t[d][d];
3447  for (unsigned int d=0, c=0; d<dim; ++d)
3448  for (unsigned int e=d+1; e<dim; ++e, ++c)
3449  array[dim+c] = (t[d][e]+t[e][d])*0.5;
3450  return SymmetricTensor<2,dim,Number>(array);
3451 }
3452 
3453 
3454 
3462 template <int rank, int dim, typename Number>
3463 inline
3466  const Number &factor)
3467 {
3469  tt *= factor;
3470  return tt;
3471 }
3472 
3473 
3474 
3482 template <int rank, int dim, typename Number>
3483 inline
3485 operator * (const Number &factor,
3487 {
3488  // simply forward to the other operator
3489  return t*factor;
3490 }
3491 
3492 
3493 
3519 template <int rank, int dim, typename Number, typename OtherNumber>
3520 inline
3523  const OtherNumber &factor)
3524 {
3525  // form the product. we have to convert the two factors into the final
3526  // type via explicit casts because, for awkward reasons, the C++
3527  // standard committee saw it fit to not define an
3528  // operator*(float,std::complex<double>)
3529  // (as well as with switched arguments and double<->float).
3530  typedef typename ProductType<Number,OtherNumber>::type product_type;
3532  // we used to shorten the following by 'tt *= product_type(factor);'
3533  // which requires that a converting constructor
3534  // 'product_type::product_type(const OtherNumber) is defined.
3535  // however, a user-defined constructor is not allowed for aggregates,
3536  // e.g. VectorizedArray. therefore, we work around this issue using a
3537  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3538  // which we assume to be defined.
3539  product_type new_factor;
3540  new_factor = factor;
3541  tt *= new_factor;
3542  return tt;
3543 }
3544 
3545 
3546 
3555 template <int rank, int dim, typename Number, typename OtherNumber>
3556 inline
3558 operator * (const Number &factor,
3560 {
3561  // simply forward to the other operator with switched arguments
3562  return (t*factor);
3563 }
3564 
3565 
3566 
3572 template <int rank, int dim, typename Number, typename OtherNumber>
3573 inline
3576  const OtherNumber &factor)
3577 {
3579  tt /= factor;
3580  return tt;
3581 }
3582 
3583 
3584 
3591 template <int rank, int dim>
3592 inline
3594 operator * (const SymmetricTensor<rank,dim> &t,
3595  const double factor)
3596 {
3598  tt *= factor;
3599  return tt;
3600 }
3601 
3602 
3603 
3610 template <int rank, int dim>
3611 inline
3613 operator * (const double factor,
3614  const SymmetricTensor<rank,dim> &t)
3615 {
3617  tt *= factor;
3618  return tt;
3619 }
3620 
3621 
3622 
3628 template <int rank, int dim>
3629 inline
3631 operator / (const SymmetricTensor<rank,dim> &t,
3632  const double factor)
3633 {
3635  tt /= factor;
3636  return tt;
3637 }
3638 
3648 template <int dim, typename Number, typename OtherNumber>
3649 inline
3653 {
3654  return (t1*t2);
3655 }
3656 
3657 
3667 template <int dim, typename Number, typename OtherNumber>
3668 inline
3671  const Tensor<2,dim,OtherNumber> &t2)
3672 {
3673  Number s = 0;
3674  for (unsigned int i=0; i<dim; ++i)
3675  for (unsigned int j=0; j<dim; ++j)
3676  s += t1[i][j] * t2[i][j];
3677  return s;
3678 }
3679 
3680 
3690 template <int dim, typename Number, typename OtherNumber>
3691 inline
3695 {
3696  return scalar_product(t2, t1);
3697 }
3698 
3699 
3715 template <typename Number, typename OtherNumber>
3716 inline
3717 void
3719  const SymmetricTensor<4,1,Number> &t,
3721 {
3722  tmp[0][0] = t[0][0][0][0] * s[0][0];
3723 }
3724 
3725 
3726 
3742 template <typename Number, typename OtherNumber>
3743 inline
3744 void
3746  const SymmetricTensor<2,1,Number> &s,
3748 {
3749  tmp[0][0] = t[0][0][0][0] * s[0][0];
3750 }
3751 
3752 
3753 
3768 template <typename Number, typename OtherNumber>
3769 inline
3770 void
3772  const SymmetricTensor<4,2,Number> &t,
3774 {
3775  const unsigned int dim = 2;
3776 
3777  for (unsigned int i=0; i<dim; ++i)
3778  for (unsigned int j=i; j<dim; ++j)
3779  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3780  t[i][j][1][1] * s[1][1] +
3781  2 * t[i][j][0][1] * s[0][1];
3782 }
3783 
3784 
3785 
3801 template <typename Number, typename OtherNumber>
3802 inline
3803 void
3805  const SymmetricTensor<2,2,Number> &s,
3807 {
3808  const unsigned int dim = 2;
3809 
3810  for (unsigned int i=0; i<dim; ++i)
3811  for (unsigned int j=i; j<dim; ++j)
3812  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3813  s[1][1] * t[1][1][i][j] +
3814  2 * s[0][1] * t[0][1][i][j];
3815 }
3816 
3817 
3818 
3834 template <typename Number, typename OtherNumber>
3835 inline
3836 void
3838  const SymmetricTensor<4,3,Number> &t,
3840 {
3841  const unsigned int dim = 3;
3842 
3843  for (unsigned int i=0; i<dim; ++i)
3844  for (unsigned int j=i; j<dim; ++j)
3845  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3846  t[i][j][1][1] * s[1][1] +
3847  t[i][j][2][2] * s[2][2] +
3848  2 * t[i][j][0][1] * s[0][1] +
3849  2 * t[i][j][0][2] * s[0][2] +
3850  2 * t[i][j][1][2] * s[1][2];
3851 }
3852 
3853 
3854 
3870 template <typename Number, typename OtherNumber>
3871 inline
3872 void
3874  const SymmetricTensor<2,3,Number> &s,
3876 {
3877  const unsigned int dim = 3;
3878 
3879  for (unsigned int i=0; i<dim; ++i)
3880  for (unsigned int j=i; j<dim; ++j)
3881  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3882  s[1][1] * t[1][1][i][j] +
3883  s[2][2] * t[2][2][i][j] +
3884  2 * s[0][1] * t[0][1][i][j] +
3885  2 * s[0][2] * t[0][2][i][j] +
3886  2 * s[1][2] * t[1][2][i][j];
3887 }
3888 
3889 
3890 
3898 template <int dim, typename Number, typename OtherNumber>
3900 operator * (const SymmetricTensor<2,dim,Number> &src1,
3901  const Tensor<1,dim,OtherNumber> &src2)
3902 {
3904  for (unsigned int i=0; i<dim; ++i)
3905  for (unsigned int j=0; j<dim; ++j)
3906  dest[i] += src1[i][j] * src2[j];
3907  return dest;
3908 }
3909 
3910 
3918 template <int dim, typename Number, typename OtherNumber>
3920 operator * (const Tensor<1,dim,Number> &src1,
3922 {
3923  // this is easy for symmetric tensors:
3924  return src2 * src1;
3925 }
3926 
3927 
3928 
3949 template <int rank_1, int rank_2, int dim,
3950  typename Number, typename OtherNumber>
3951 inline DEAL_II_ALWAYS_INLINE
3952 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3953 operator * (const Tensor<rank_1, dim, Number> &src1,
3955 {
3956  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3957  const Tensor<rank_2, dim, OtherNumber> src2 (src2s);
3958  return src1*src2;
3959 }
3960 
3961 
3962 
3983 template <int rank_1, int rank_2, int dim,
3984  typename Number, typename OtherNumber>
3985 inline DEAL_II_ALWAYS_INLINE
3986 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3987 operator * (const SymmetricTensor<rank_1, dim, Number> &src1s,
3989 {
3990  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3991  const Tensor<rank_2, dim, OtherNumber> src1 (src1s);
3992  return src1*src2;
3993 }
3994 
3995 
3996 
4006 template <int dim, typename Number>
4007 inline
4008 std::ostream &operator << (std::ostream &out,
4010 {
4011  //make out lives a bit simpler by outputing
4012  //the tensor through the operator for the
4013  //general Tensor class
4015 
4016  for (unsigned int i=0; i<dim; ++i)
4017  for (unsigned int j=0; j<dim; ++j)
4018  tt[i][j] = t[i][j];
4019 
4020  return out << tt;
4021 }
4022 
4023 
4024 
4034 template <int dim, typename Number>
4035 inline
4036 std::ostream &operator << (std::ostream &out,
4038 {
4039  //make out lives a bit simpler by outputing
4040  //the tensor through the operator for the
4041  //general Tensor class
4043 
4044  for (unsigned int i=0; i<dim; ++i)
4045  for (unsigned int j=0; j<dim; ++j)
4046  for (unsigned int k=0; k<dim; ++k)
4047  for (unsigned int l=0; l<dim; ++l)
4048  tt[i][j][k][l] = t[i][j][k][l];
4049 
4050  return out << tt;
4051 }
4052 
4053 
4054 DEAL_II_NAMESPACE_CLOSE
4055 
4056 #endif
numbers::NumberTraits< Number >::real_type norm() const
internal::SymmetricTensorAccessors::Accessor< rank, dim, true, rank-1, Number > operator[](const unsigned int row) const
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
bool operator!=(const SymmetricTensor &) const
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank, dim, Number > &left, const SymmetricTensor< rank, dim, OtherNumber > &right)
bool operator==(const SymmetricTensor &) const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::array< std::pair< Number, Tensor< 1, 1, Number > >, 1 > eigenvectors(const SymmetricTensor< 2, 1, Number > &T, const enum SymmetricTensorEigenvectorMethod=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
static const unsigned int n_independent_components
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
SymmetricTensor & operator=(const SymmetricTensor< rank, dim, OtherNumber > &rhs)
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:1199
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, OtherNumber > &t2)
static TableIndices< rank > unrolled_to_component_indices(const unsigned int i)
base_tensor_type data
SymmetricTensor & operator/=(const OtherNumber &factor)
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
static real_type abs(const number &x)
Definition: numbers.h:351
SymmetricTensorEigenvectorMethod
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, OtherNumber > &s)
void serialize(Archive &ar, const unsigned int version)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const enum SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
static const unsigned int dimension
static unsigned int component_to_unrolled_index(const TableIndices< rank > &indices)
static::ExceptionBase & ExcMessage(std::string arg1)
static std::size_t memory_consumption()
TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
SymmetricTensor & operator+=(const SymmetricTensor< rank, dim, OtherNumber > &)
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)
friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
T sum(const T &t, const MPI_Comm &mpi_communicator)
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:337
base_tensor_descriptor::base_tensor_type base_tensor_type
SymmetricTensor< rank, dim, Number > transpose(const SymmetricTensor< rank, dim, Number > &t)
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Number trace(const SymmetricTensor< 2, dim, Number > &d)
Number access_raw_entry(const unsigned int unrolled_index) const
internal::SymmetricTensorAccessors::StorageType< rank, dim, Number > base_tensor_descriptor
friend SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank, dim, Number > &left, const SymmetricTensor< rank, dim, OtherNumber > &right)
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
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)
SymmetricTensor operator-() const
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: mpi.h:41
Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Number & operator()(const TableIndices< rank > &indices)
Tensor< 1, n_independent_components, Number > base_tensor_type
static::ExceptionBase & ExcNotImplemented()
SymmetricTensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank, dim, Number > &left, const SymmetricTensor< rank, dim, OtherNumber > &right)
SymmetricTensor & operator*=(const OtherNumber &factor)
SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
SymmetricTensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank, dim, Number > &left, const SymmetricTensor< rank, dim, OtherNumber > &right)
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()
StreamType & operator<<(StreamType &s, UpdateFlags u)
internal::SymmetricTensorAccessors::double_contraction_result< rank, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, OtherNumber > &t)
SymmetricTensor & operator-=(const SymmetricTensor< rank, dim, OtherNumber > &)
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 > &)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static::ExceptionBase & ExcInternalError()