Reference documentation for deal.II version Git 3e8fd02 2017-09-22 11:37:14 +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)
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>::value_type value_type;
1542  typedef typename SymmetricTensorAccessors::StorageType<2,dim,value_type>::base_tensor_type base_tensor_type;
1543 
1544  base_tensor_type tmp;
1545  for (unsigned int i=0; i<tmp.dimension; ++i)
1546  {
1547  // Start with the non-diagonal part
1548  value_type sum = data[dim] * sdata[dim][i];
1549  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1550  sum += data[d] * sdata[d][i];
1551  sum += sum; // sum = sum * 2.;
1552 
1553  // Now add the contributions from the diagonal
1554  for (unsigned int d=0; d<dim; ++d)
1555  sum += data[d] * sdata[d][i];
1556  tmp[i] = sum;
1557  }
1558  return tmp;
1559  }
1560 
1561 
1562 
1563  template <int dim, typename Number, typename OtherNumber = Number>
1564  inline
1565  typename SymmetricTensorAccessors::StorageType<4,dim,
1566  typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type
1567  >::base_tensor_type
1568  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,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<4,4,dim,Number,OtherNumber>::value_type value_type;
1572  typedef typename SymmetricTensorAccessors::StorageType<4,dim,value_type>::base_tensor_type base_tensor_type;
1573 
1574  const unsigned int data_dim =
1575  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1576  base_tensor_type tmp;
1577  for (unsigned int i=0; i<data_dim; ++i)
1578  for (unsigned int j=0; j<data_dim; ++j)
1579  {
1580  // Start with the non-diagonal part
1581  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1582  tmp[i][j] += data[i][d] * sdata[d][j];
1583  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1584 
1585  // Now add the contributions from the diagonal
1586  for (unsigned int d=0; d<dim; ++d)
1587  tmp[i][j] += data[i][d] * sdata[d][j];
1588  }
1589  return tmp;
1590  }
1591 
1592 } // end of namespace internal
1593 
1594 
1595 
1596 template <int rank, int dim, typename Number>
1597 template <typename OtherNumber>
1598 inline
1601 {
1602  // need to have two different function calls
1603  // because a scalar and rank-2 tensor are not
1604  // the same data type (see internal function
1605  // above)
1606  return internal::perform_double_contraction<dim,Number,OtherNumber> (data, s.data);
1607 }
1608 
1609 
1610 
1611 template <int rank, int dim, typename Number>
1612 template <typename OtherNumber>
1613 inline
1616 {
1619  tmp.data = internal::perform_double_contraction<dim,Number,OtherNumber> (data,s.data);
1620  return tmp;
1621 }
1622 
1623 
1624 
1625 // internal namespace to switch between the
1626 // access of different tensors. There used to
1627 // be explicit instantiations before for
1628 // different ranks and dimensions, but since
1629 // we now allow for templates on the data
1630 // type, and since we cannot partially
1631 // specialize the implementation, this got
1632 // into a separate namespace
1633 namespace internal
1634 {
1635  template <int dim, typename Number>
1636  inline
1637  Number &
1638  symmetric_tensor_access (const TableIndices<2> &indices,
1640  {
1641  // 1d is very simple and done first
1642  if (dim == 1)
1643  return data[0];
1644 
1645  // first treat the main diagonal elements, which are stored consecutively
1646  // at the beginning
1647  if (indices[0] == indices[1])
1648  return data[indices[0]];
1649 
1650  // the rest is messier and requires a few switches.
1651  switch (dim)
1652  {
1653  case 2:
1654  // at least for the 2x2 case it is reasonably simple
1655  Assert (((indices[0]==1) && (indices[1]==0)) ||
1656  ((indices[0]==0) && (indices[1]==1)),
1657  ExcInternalError());
1658  return data[2];
1659 
1660  default:
1661  // to do the rest, sort our indices before comparing
1662  {
1663  TableIndices<2> sorted_indices (indices);
1664  sorted_indices.sort ();
1665 
1666  for (unsigned int d=0, c=0; d<dim; ++d)
1667  for (unsigned int e=d+1; e<dim; ++e, ++c)
1668  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1669  return data[dim+c];
1670  Assert (false, ExcInternalError());
1671  }
1672  }
1673 
1674  static Number dummy_but_referenceable = Number();
1675  return dummy_but_referenceable;
1676  }
1677 
1678 
1679 
1680  template <int dim, typename Number>
1681  inline
1682  Number
1683  symmetric_tensor_access (const TableIndices<2> &indices,
1685  {
1686  // 1d is very simple and done first
1687  if (dim == 1)
1688  return data[0];
1689 
1690  // first treat the main diagonal elements, which are stored consecutively
1691  // at the beginning
1692  if (indices[0] == indices[1])
1693  return data[indices[0]];
1694 
1695  // the rest is messier and requires a few switches.
1696  switch (dim)
1697  {
1698  case 2:
1699  // at least for the 2x2 case it is reasonably simple
1700  Assert (((indices[0]==1) && (indices[1]==0)) ||
1701  ((indices[0]==0) && (indices[1]==1)),
1702  ExcInternalError());
1703  return data[2];
1704 
1705  default:
1706  // to do the rest, sort our indices before comparing
1707  {
1708  TableIndices<2> sorted_indices (indices);
1709  sorted_indices.sort ();
1710 
1711  for (unsigned int d=0, c=0; d<dim; ++d)
1712  for (unsigned int e=d+1; e<dim; ++e, ++c)
1713  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1714  return data[dim+c];
1715  Assert (false, ExcInternalError());
1716  }
1717  }
1718 
1719  static Number dummy_but_referenceable = Number();
1720  return dummy_but_referenceable;
1721  }
1722 
1723 
1724 
1725  template <int dim, typename Number>
1726  inline
1727  Number &
1728  symmetric_tensor_access (const TableIndices<4> &indices,
1730  {
1731  switch (dim)
1732  {
1733  case 1:
1734  return data[0][0];
1735 
1736  case 2:
1737  // each entry of the tensor can be
1738  // thought of as an entry in a
1739  // matrix that maps the rolled-out
1740  // rank-2 tensors into rolled-out
1741  // rank-2 tensors. this is the
1742  // format in which we store rank-4
1743  // tensors. determine which
1744  // position the present entry is
1745  // stored in
1746  {
1747  unsigned int base_index[2] ;
1748  if ((indices[0] == 0) && (indices[1] == 0))
1749  base_index[0] = 0;
1750  else if ((indices[0] == 1) && (indices[1] == 1))
1751  base_index[0] = 1;
1752  else
1753  base_index[0] = 2;
1754 
1755  if ((indices[2] == 0) && (indices[3] == 0))
1756  base_index[1] = 0;
1757  else if ((indices[2] == 1) && (indices[3] == 1))
1758  base_index[1] = 1;
1759  else
1760  base_index[1] = 2;
1761 
1762  return data[base_index[0]][base_index[1]];
1763  }
1764 
1765  case 3:
1766  // each entry of the tensor can be
1767  // thought of as an entry in a
1768  // matrix that maps the rolled-out
1769  // rank-2 tensors into rolled-out
1770  // rank-2 tensors. this is the
1771  // format in which we store rank-4
1772  // tensors. determine which
1773  // position the present entry is
1774  // stored in
1775  {
1776  unsigned int base_index[2] ;
1777  if ((indices[0] == 0) && (indices[1] == 0))
1778  base_index[0] = 0;
1779  else if ((indices[0] == 1) && (indices[1] == 1))
1780  base_index[0] = 1;
1781  else if ((indices[0] == 2) && (indices[1] == 2))
1782  base_index[0] = 2;
1783  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1784  ((indices[0] == 1) && (indices[1] == 0)))
1785  base_index[0] = 3;
1786  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1787  ((indices[0] == 2) && (indices[1] == 0)))
1788  base_index[0] = 4;
1789  else
1790  {
1791  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1792  ((indices[0] == 2) && (indices[1] == 1)),
1793  ExcInternalError());
1794  base_index[0] = 5;
1795  }
1796 
1797  if ((indices[2] == 0) && (indices[3] == 0))
1798  base_index[1] = 0;
1799  else if ((indices[2] == 1) && (indices[3] == 1))
1800  base_index[1] = 1;
1801  else if ((indices[2] == 2) && (indices[3] == 2))
1802  base_index[1] = 2;
1803  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1804  ((indices[2] == 1) && (indices[3] == 0)))
1805  base_index[1] = 3;
1806  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1807  ((indices[2] == 2) && (indices[3] == 0)))
1808  base_index[1] = 4;
1809  else
1810  {
1811  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1812  ((indices[2] == 2) && (indices[3] == 1)),
1813  ExcInternalError());
1814  base_index[1] = 5;
1815  }
1816 
1817  return data[base_index[0]][base_index[1]];
1818  }
1819 
1820  default:
1821  Assert (false, ExcNotImplemented());
1822  }
1823 
1824  static Number dummy;
1825  return dummy;
1826  }
1827 
1828 
1829  template <int dim, typename Number>
1830  inline
1831  Number
1832  symmetric_tensor_access (const TableIndices<4> &indices,
1834  {
1835  switch (dim)
1836  {
1837  case 1:
1838  return data[0][0];
1839 
1840  case 2:
1841  // each entry of the tensor can be
1842  // thought of as an entry in a
1843  // matrix that maps the rolled-out
1844  // rank-2 tensors into rolled-out
1845  // rank-2 tensors. this is the
1846  // format in which we store rank-4
1847  // tensors. determine which
1848  // position the present entry is
1849  // stored in
1850  {
1851  unsigned int base_index[2] ;
1852  if ((indices[0] == 0) && (indices[1] == 0))
1853  base_index[0] = 0;
1854  else if ((indices[0] == 1) && (indices[1] == 1))
1855  base_index[0] = 1;
1856  else
1857  base_index[0] = 2;
1858 
1859  if ((indices[2] == 0) && (indices[3] == 0))
1860  base_index[1] = 0;
1861  else if ((indices[2] == 1) && (indices[3] == 1))
1862  base_index[1] = 1;
1863  else
1864  base_index[1] = 2;
1865 
1866  return data[base_index[0]][base_index[1]];
1867  }
1868 
1869  case 3:
1870  // each entry of the tensor can be
1871  // thought of as an entry in a
1872  // matrix that maps the rolled-out
1873  // rank-2 tensors into rolled-out
1874  // rank-2 tensors. this is the
1875  // format in which we store rank-4
1876  // tensors. determine which
1877  // position the present entry is
1878  // stored in
1879  {
1880  unsigned int base_index[2] ;
1881  if ((indices[0] == 0) && (indices[1] == 0))
1882  base_index[0] = 0;
1883  else if ((indices[0] == 1) && (indices[1] == 1))
1884  base_index[0] = 1;
1885  else if ((indices[0] == 2) && (indices[1] == 2))
1886  base_index[0] = 2;
1887  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1888  ((indices[0] == 1) && (indices[1] == 0)))
1889  base_index[0] = 3;
1890  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1891  ((indices[0] == 2) && (indices[1] == 0)))
1892  base_index[0] = 4;
1893  else
1894  {
1895  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1896  ((indices[0] == 2) && (indices[1] == 1)),
1897  ExcInternalError());
1898  base_index[0] = 5;
1899  }
1900 
1901  if ((indices[2] == 0) && (indices[3] == 0))
1902  base_index[1] = 0;
1903  else if ((indices[2] == 1) && (indices[3] == 1))
1904  base_index[1] = 1;
1905  else if ((indices[2] == 2) && (indices[3] == 2))
1906  base_index[1] = 2;
1907  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1908  ((indices[2] == 1) && (indices[3] == 0)))
1909  base_index[1] = 3;
1910  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1911  ((indices[2] == 2) && (indices[3] == 0)))
1912  base_index[1] = 4;
1913  else
1914  {
1915  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1916  ((indices[2] == 2) && (indices[3] == 1)),
1917  ExcInternalError());
1918  base_index[1] = 5;
1919  }
1920 
1921  return data[base_index[0]][base_index[1]];
1922  }
1923 
1924  default:
1925  Assert (false, ExcNotImplemented());
1926  }
1927 
1928  static Number dummy;
1929  return dummy;
1930  }
1931 
1932 } // end of namespace internal
1933 
1934 
1935 
1936 template <int rank, int dim, typename Number>
1937 inline
1938 Number &
1940 {
1941  for (unsigned int r=0; r<rank; ++r)
1942  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1943  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1944 }
1945 
1946 
1947 
1948 template <int rank, int dim, typename Number>
1949 inline
1950 Number
1952 (const TableIndices<rank> &indices) const
1953 {
1954  for (unsigned int r=0; r<rank; ++r)
1955  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1956  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1957 }
1958 
1959 
1960 
1961 namespace internal
1962 {
1963  namespace SymmetricTensor
1964  {
1965  template <int rank>
1967  get_partially_filled_indices (const unsigned int row,
1968  const std::integral_constant<int, 2> &)
1969  {
1970  return TableIndices<rank> (row,
1972 
1973  }
1974 
1975 
1976  template <int rank>
1978  get_partially_filled_indices (const unsigned int row,
1979  const std::integral_constant<int, 4> &)
1980  {
1981  return TableIndices<rank> (row,
1985 
1986  }
1987  }
1988 }
1989 
1990 
1991 template <int rank, int dim, typename Number>
1992 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1993 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1994 {
1995  return
1996  internal::SymmetricTensorAccessors::
1997  Accessor<rank,dim,true,rank-1,Number> (*this,
1998  internal::SymmetricTensor::get_partially_filled_indices<rank> (row,
1999  std::integral_constant<int, rank>()));
2000 }
2001 
2002 
2003 
2004 template <int rank, int dim, typename Number>
2005 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
2006 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
2007 {
2008  return
2009  internal::SymmetricTensorAccessors::
2010  Accessor<rank,dim,false,rank-1,Number> (*this,
2011  internal::SymmetricTensor::get_partially_filled_indices<rank> (row,
2012  std::integral_constant<int, rank>()));
2013 }
2014 
2015 
2016 
2017 template <int rank, int dim, typename Number>
2018 inline
2019 Number
2021 {
2022  return operator()(indices);
2023 }
2024 
2025 
2026 
2027 template <int rank, int dim, typename Number>
2028 inline
2029 Number &
2031 {
2032  return operator()(indices);
2033 }
2034 
2035 
2036 
2037 
2038 namespace internal
2039 {
2040  namespace SymmetricTensor
2041  {
2042  template <int dim, typename Number>
2043  unsigned int
2044  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
2045  const unsigned int index)
2046  {
2047  return index;
2048  }
2049 
2050 
2051  template <int dim, typename Number>
2053  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
2054  const unsigned int index)
2055  {
2056  return
2059  }
2060 
2061  }
2062 }
2063 
2064 
2065 
2066 template <int rank, int dim, typename Number>
2067 inline
2068 Number
2069 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
2070 {
2071  AssertIndexRange (index, n_independent_components);
2072  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2073 }
2074 
2075 
2076 
2077 template <int rank, int dim, typename Number>
2078 inline
2079 Number &
2080 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
2081 {
2082  AssertIndexRange (index, n_independent_components);
2083  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2084 }
2085 
2086 
2087 
2088 namespace internal
2089 {
2090  template <int dim, typename Number>
2091  inline
2094  {
2095  switch (dim)
2096  {
2097  case 1:
2098  return numbers::NumberTraits<Number>::abs(data[0]);
2099 
2100  case 2:
2101  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2104 
2105  case 3:
2106  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2112 
2113  default:
2114  {
2115  typename numbers::NumberTraits<Number>::real_type return_value
2117 
2118  for (unsigned int d=0; d<dim; ++d)
2119  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
2120  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
2121  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2122 
2123  return std::sqrt(return_value);
2124  }
2125  }
2126  }
2127 
2128 
2129 
2130  template <int dim, typename Number>
2131  inline
2134  {
2135  switch (dim)
2136  {
2137  case 1:
2138  return numbers::NumberTraits<Number>::abs (data[0][0]);
2139 
2140  default:
2141  {
2142  typename numbers::NumberTraits<Number>::real_type return_value
2144 
2145  const unsigned int n_independent_components = data.dimension;
2146 
2147  for (unsigned int i=0; i<dim; ++i)
2148  for (unsigned int j=0; j<dim; ++j)
2149  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
2150  for (unsigned int i=0; i<dim; ++i)
2151  for (unsigned int j=dim; j<n_independent_components; ++j)
2152  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2153  for (unsigned int i=dim; i<n_independent_components; ++i)
2154  for (unsigned int j=0; j<dim; ++j)
2155  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2156  for (unsigned int i=dim; i<n_independent_components; ++i)
2157  for (unsigned int j=dim; j<n_independent_components; ++j)
2158  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2159 
2160  return std::sqrt(return_value);
2161  }
2162  }
2163  }
2164 
2165 } // end of namespace internal
2166 
2167 
2168 
2169 template <int rank, int dim, typename Number>
2170 inline
2173 {
2174  return internal::compute_norm<dim,Number> (data);
2175 }
2176 
2177 
2178 
2179 namespace internal
2180 {
2181  namespace SymmetricTensor
2182  {
2183  namespace
2184  {
2185  // a function to do the unrolling from a set of indices to a
2186  // scalar index into the array in which we store the elements of
2187  // a symmetric tensor
2188  //
2189  // this function is for rank-2 tensors
2190  template <int dim>
2191  inline
2192  unsigned int
2193  component_to_unrolled_index
2194  (const TableIndices<2> &indices)
2195  {
2196  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2197  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2198 
2199  switch (dim)
2200  {
2201  case 1:
2202  {
2203  return 0;
2204  }
2205 
2206  case 2:
2207  {
2208  static const unsigned int table[2][2] = {{0, 2},
2209  {2, 1}
2210  };
2211  return table[indices[0]][indices[1]];
2212  }
2213 
2214  case 3:
2215  {
2216  static const unsigned int table[3][3] = {{0, 3, 4},
2217  {3, 1, 5},
2218  {4, 5, 2}
2219  };
2220  return table[indices[0]][indices[1]];
2221  }
2222 
2223  case 4:
2224  {
2225  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2226  {4, 1, 7, 8},
2227  {5, 7, 2, 9},
2228  {6, 8, 9, 3}
2229  };
2230  return table[indices[0]][indices[1]];
2231  }
2232 
2233  default:
2234  // for the remainder, manually figure out the numbering
2235  {
2236  if (indices[0] == indices[1])
2237  return indices[0];
2238 
2239  TableIndices<2> sorted_indices (indices);
2240  sorted_indices.sort ();
2241 
2242  for (unsigned int d=0, c=0; d<dim; ++d)
2243  for (unsigned int e=d+1; e<dim; ++e, ++c)
2244  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
2245  return dim+c;
2246 
2247  // should never get here:
2248  Assert(false, ExcInternalError());
2249  return 0;
2250  }
2251  }
2252  }
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 tensors of ranks not already handled
2259  // above
2260  template <int dim, int rank>
2261  inline
2262  unsigned int
2263  component_to_unrolled_index
2264  (const TableIndices<rank> &indices)
2265  {
2266  (void)indices;
2267  Assert (false, ExcNotImplemented());
2269  }
2270  }
2271  }
2272 }
2273 
2274 
2275 template <int rank, int dim, typename Number>
2276 inline
2277 unsigned int
2279 (const TableIndices<rank> &indices)
2280 {
2281  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
2282 }
2283 
2284 
2285 
2286 namespace internal
2287 {
2288  namespace SymmetricTensor
2289  {
2290  namespace
2291  {
2292  // a function to do the inverse of the unrolling from a set of
2293  // indices to a scalar index into the array in which we store
2294  // the elements of a symmetric tensor. in other words, it goes
2295  // from the scalar index into the array to a set of indices of
2296  // the tensor
2297  //
2298  // this function is for rank-2 tensors
2299  template <int dim>
2300  inline
2302  unrolled_to_component_indices
2303  (const unsigned int i,
2304  const std::integral_constant<int, 2> &)
2305  {
2308  switch (dim)
2309  {
2310  case 1:
2311  {
2312  return TableIndices<2>(0,0);
2313  }
2314 
2315  case 2:
2316  {
2317  const TableIndices<2> table[3] =
2318  {
2319  TableIndices<2> (0,0),
2320  TableIndices<2> (1,1),
2321  TableIndices<2> (0,1)
2322  };
2323  return table[i];
2324  }
2325 
2326  case 3:
2327  {
2328  const TableIndices<2> table[6] =
2329  {
2330  TableIndices<2> (0,0),
2331  TableIndices<2> (1,1),
2332  TableIndices<2> (2,2),
2333  TableIndices<2> (0,1),
2334  TableIndices<2> (0,2),
2335  TableIndices<2> (1,2)
2336  };
2337  return table[i];
2338  }
2339 
2340  default:
2341  if (i<dim)
2342  return TableIndices<2> (i,i);
2343 
2344  for (unsigned int d=0, c=0; d<dim; ++d)
2345  for (unsigned int e=d+1; e<dim; ++e, ++c)
2346  if (c==i)
2347  return TableIndices<2>(d,e);
2348 
2349  // should never get here:
2350  Assert(false, ExcInternalError());
2351  return TableIndices<2>(0, 0);
2352  }
2353  }
2354 
2355  // a function to do the inverse of the unrolling from a set of
2356  // indices to a scalar index into the array in which we store
2357  // the elements of a symmetric tensor. in other words, it goes
2358  // from the scalar index into the array to a set of indices of
2359  // the tensor
2360  //
2361  // this function is for tensors of a rank not already handled
2362  // above
2363  template <int dim, int rank>
2364  inline
2366  unrolled_to_component_indices
2367  (const unsigned int i,
2368  const std::integral_constant<int, rank> &)
2369  {
2370  (void)i;
2373  Assert (false, ExcNotImplemented());
2374  return TableIndices<rank>();
2375  }
2376 
2377  }
2378  }
2379 }
2380 
2381 template <int rank, int dim, typename Number>
2382 inline
2385 (const unsigned int i)
2386 {
2387  return
2388  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2389  std::integral_constant<int, rank>());
2390 }
2391 
2392 
2393 
2394 template <int rank, int dim, typename Number>
2395 template <class Archive>
2396 inline
2397 void
2398 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2399 {
2400  ar &data;
2401 }
2402 
2403 
2404 #endif // DOXYGEN
2405 
2406 /* ----------------- Non-member functions operating on tensors. ------------ */
2407 
2408 
2421 template <int rank, int dim, typename Number, typename OtherNumber>
2422 inline
2426 {
2428  tmp += right;
2429  return tmp;
2430 }
2431 
2432 
2445 template <int rank, int dim, typename Number, typename OtherNumber>
2446 inline
2450 {
2452  tmp -= right;
2453  return tmp;
2454 }
2455 
2456 
2464 template <int rank, int dim, typename Number, typename OtherNumber>
2465 inline
2468  const Tensor<rank, dim, OtherNumber> &right)
2469 {
2470  return Tensor<rank, dim, Number>(left) + right;
2471 }
2472 
2473 
2481 template <int rank, int dim, typename Number, typename OtherNumber>
2482 inline
2486 {
2487  return left + Tensor<rank, dim, OtherNumber>(right);
2488 }
2489 
2490 
2498 template <int rank, int dim, typename Number, typename OtherNumber>
2499 inline
2502  const Tensor<rank, dim, OtherNumber> &right)
2503 {
2504  return Tensor<rank, dim, Number>(left) - right;
2505 }
2506 
2507 
2515 template <int rank, int dim, typename Number, typename OtherNumber>
2516 inline
2520 {
2521  return left - Tensor<rank, dim, OtherNumber>(right);
2522 }
2523 
2524 
2525 
2539 template <int dim, typename Number>
2540 inline
2542 {
2543  switch (dim)
2544  {
2545  case 1:
2546  return t.data[0];
2547  case 2:
2548  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2549  case 3:
2550  {
2551  // in analogy to general tensors, but
2552  // there's something to be simplified for
2553  // the present case
2554  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2555  return ( tmp + tmp
2556  +t.data[0]*t.data[1]*t.data[2]
2557  -t.data[0]*t.data[5]*t.data[5]
2558  -t.data[1]*t.data[4]*t.data[4]
2559  -t.data[2]*t.data[3]*t.data[3]);
2560  }
2561  default:
2562  Assert (false, ExcNotImplemented());
2564  }
2565 }
2566 
2567 
2568 
2578 template <int dim, typename Number>
2579 inline
2581 {
2582  return determinant (t);
2583 }
2584 
2585 
2586 
2594 template <int dim, typename Number>
2596 {
2597  Number t = d.data[0];
2598  for (unsigned int i=1; i<dim; ++i)
2599  t += d.data[i];
2600  return t;
2601 }
2602 
2603 
2613 template <int dim, typename Number>
2614 inline
2616 {
2617  return trace (t);
2618 }
2619 
2620 
2633 template <typename Number>
2634 inline
2636 {
2638 }
2639 
2640 
2641 
2662 template <typename Number>
2663 inline
2665 {
2666  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2667 }
2668 
2669 
2670 
2680 template <typename Number>
2681 inline
2683 {
2684  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2685  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2686 }
2687 
2688 
2689 
2698 template <typename Number>
2699 std::array<Number,1>
2700 eigenvalues (const SymmetricTensor<2,1,Number> &T);
2701 
2702 
2703 
2725 template <typename Number>
2726 std::array<Number,2>
2727 eigenvalues (const SymmetricTensor<2,2,Number> &T);
2728 
2729 
2730 
2750 template <typename Number>
2751 std::array<Number,3>
2752 eigenvalues (const SymmetricTensor<2,3,Number> &T);
2753 
2754 
2755 
2756 namespace internal
2757 {
2758  namespace SymmetricTensor
2759  {
2796  template <int dim, typename Number>
2797  void
2798  tridiagonalize (const ::SymmetricTensor<2,dim,Number> &A,
2799  ::Tensor<2,dim,Number> &Q,
2800  std::array<Number,dim> &d,
2801  std::array<Number,dim-1> &e);
2802 
2803 
2804 
2843  template <int dim, typename Number>
2844  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2845  ql_implicit_shifts (const ::SymmetricTensor<2,dim,Number> &A);
2846 
2847 
2848 
2887  template <int dim, typename Number>
2888  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2890 
2891 
2892 
2908  template <typename Number>
2909  std::array<std::pair<Number, Tensor<1,2,Number> >,2>
2910  hybrid (const ::SymmetricTensor<2,2,Number> &A);
2911 
2912 
2913 
2946  template <typename Number>
2947  std::array<std::pair<Number, Tensor<1,3,Number> >,3>
2948  hybrid (const ::SymmetricTensor<2,3,Number> &A);
2949 
2950  namespace
2951  {
2952 
2957  template <int dim, typename Number>
2958  struct SortEigenValuesVectors
2959  {
2960  typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
2961  bool operator() (const EigValsVecs &lhs,
2962  const EigValsVecs &rhs)
2963  {
2964  return lhs.first > rhs.first;
2965  }
2966  };
2967 
2968  }
2969 
2970  } // namespace SymmetricTensor
2971 
2972 } // namespace internal
2973 
2974 
2975 
2976 // The line below is to ensure that doxygen puts the full description
2977 // of this global enumeration into the documentation
2978 // See https://stackoverflow.com/a/1717984
3006 {
3016  hybrid,
3026  ql_implicit_shifts,
3034  jacobi
3035 };
3036 
3037 
3038 
3045 template <typename Number>
3046 std::array<std::pair<Number, Tensor<1,1,Number> >,1>
3049 {
3050  return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} };
3051 }
3052 
3053 
3054 
3082 template <int dim, typename Number>
3083 std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
3086 {
3087  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
3088 
3089  switch (method)
3090  {
3092  eig_vals_vecs = internal::SymmetricTensor::hybrid(T);
3093  break;
3096  break;
3098  eig_vals_vecs = internal::SymmetricTensor::jacobi(T);
3099  break;
3100  default:
3101  AssertThrow(false, ExcNotImplemented());
3102  }
3103 
3104  // Sort in descending order before output.
3105  std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
3106  internal::SymmetricTensor::SortEigenValuesVectors<dim,Number>());
3107  return eig_vals_vecs;
3108 }
3109 
3110 
3111 
3121 template <int rank, int dim, typename Number>
3122 inline
3125 {
3126  return t;
3127 }
3128 
3129 
3130 
3140 template <int dim, typename Number>
3141 inline
3144 {
3146 
3147  // subtract scaled trace from the diagonal
3148  const Number tr = trace(t) / dim;
3149  for (unsigned int i=0; i<dim; ++i)
3150  tmp.data[i] -= tr;
3151 
3152  return tmp;
3153 }
3154 
3155 
3156 
3164 template <int dim, typename Number>
3165 inline
3167 unit_symmetric_tensor ()
3168 {
3169  // create a default constructed matrix filled with
3170  // zeros, then set the diagonal elements to one
3172  switch (dim)
3173  {
3174  case 1:
3175  tmp.data[0] = 1;
3176  break;
3177  case 2:
3178  tmp.data[0] = tmp.data[1] = 1;
3179  break;
3180  case 3:
3181  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3182  break;
3183  default:
3184  for (unsigned int d=0; d<dim; ++d)
3185  tmp.data[d] = 1;
3186  }
3187  return tmp;
3188 }
3189 
3190 
3191 
3200 template <int dim>
3201 inline
3203 unit_symmetric_tensor ()
3204 {
3205  return unit_symmetric_tensor<dim,double>();
3206 }
3207 
3208 
3209 
3224 template <int dim, typename Number>
3225 inline
3227 deviator_tensor ()
3228 {
3230 
3231  // fill the elements treating the diagonal
3232  for (unsigned int i=0; i<dim; ++i)
3233  for (unsigned int j=0; j<dim; ++j)
3234  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
3235 
3236  // then fill the ones that copy over the
3237  // non-diagonal elements. note that during
3238  // the double-contraction, we handle the
3239  // off-diagonal elements twice, so simply
3240  // copying requires a weight of 1/2
3241  for (unsigned int i=dim;
3242  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3243  ++i)
3244  tmp.data[i][i] = 0.5;
3245 
3246  return tmp;
3247 }
3248 
3249 
3250 
3265 template <int dim>
3266 inline
3268 deviator_tensor ()
3269 {
3270  return deviator_tensor<dim,double>();
3271 }
3272 
3273 
3274 
3297 template <int dim, typename Number>
3298 inline
3300 identity_tensor ()
3301 {
3303 
3304  // fill the elements treating the diagonal
3305  for (unsigned int i=0; i<dim; ++i)
3306  tmp.data[i][i] = 1;
3307 
3308  // then fill the ones that copy over the
3309  // non-diagonal elements. note that during
3310  // the double-contraction, we handle the
3311  // off-diagonal elements twice, so simply
3312  // copying requires a weight of 1/2
3313  for (unsigned int i=dim;
3314  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3315  ++i)
3316  tmp.data[i][i] = 0.5;
3317 
3318  return tmp;
3319 }
3320 
3321 
3322 
3344 template <int dim>
3345 inline
3347 identity_tensor ()
3348 {
3349  return identity_tensor<dim,double>();
3350 }
3351 
3352 
3353 
3364 template <int dim, typename Number>
3365 inline
3368 {
3370 }
3371 
3372 
3373 
3385 template <int dim, typename Number>
3386 inline
3389 {
3391 }
3392 
3393 
3394 
3409 template <int dim, typename Number>
3410 inline
3414 {
3416 
3417  // fill only the elements really needed
3418  for (unsigned int i=0; i<dim; ++i)
3419  for (unsigned int j=i; j<dim; ++j)
3420  for (unsigned int k=0; k<dim; ++k)
3421  for (unsigned int l=k; l<dim; ++l)
3422  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3423 
3424  return tmp;
3425 }
3426 
3427 
3428 
3437 template <int dim,typename Number>
3438 inline
3441 {
3442  Number array[(dim*dim+dim)/2];
3443  for (unsigned int d=0; d<dim; ++d)
3444  array[d] = t[d][d];
3445  for (unsigned int d=0, c=0; d<dim; ++d)
3446  for (unsigned int e=d+1; e<dim; ++e, ++c)
3447  array[dim+c] = (t[d][e]+t[e][d])*0.5;
3448  return SymmetricTensor<2,dim,Number>(array);
3449 }
3450 
3451 
3452 
3460 template <int rank, int dim, typename Number>
3461 inline
3464  const Number &factor)
3465 {
3467  tt *= factor;
3468  return tt;
3469 }
3470 
3471 
3472 
3480 template <int rank, int dim, typename Number>
3481 inline
3483 operator * (const Number &factor,
3485 {
3486  // simply forward to the other operator
3487  return t*factor;
3488 }
3489 
3490 
3491 
3517 template <int rank, int dim, typename Number, typename OtherNumber>
3518 inline
3521  const OtherNumber &factor)
3522 {
3523  // form the product. we have to convert the two factors into the final
3524  // type via explicit casts because, for awkward reasons, the C++
3525  // standard committee saw it fit to not define an
3526  // operator*(float,std::complex<double>)
3527  // (as well as with switched arguments and double<->float).
3528  typedef typename ProductType<Number,OtherNumber>::type product_type;
3530  // we used to shorten the following by 'tt *= product_type(factor);'
3531  // which requires that a converting constructor
3532  // 'product_type::product_type(const OtherNumber) is defined.
3533  // however, a user-defined constructor is not allowed for aggregates,
3534  // e.g. VectorizedArray. therefore, we work around this issue using a
3535  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3536  // which we assume to be defined.
3537  product_type new_factor;
3538  new_factor = factor;
3539  tt *= new_factor;
3540  return tt;
3541 }
3542 
3543 
3544 
3553 template <int rank, int dim, typename Number, typename OtherNumber>
3554 inline
3556 operator * (const Number &factor,
3558 {
3559  // simply forward to the other operator with switched arguments
3560  return (t*factor);
3561 }
3562 
3563 
3564 
3570 template <int rank, int dim, typename Number, typename OtherNumber>
3571 inline
3574  const OtherNumber &factor)
3575 {
3577  tt /= factor;
3578  return tt;
3579 }
3580 
3581 
3582 
3589 template <int rank, int dim>
3590 inline
3592 operator * (const SymmetricTensor<rank,dim> &t,
3593  const double factor)
3594 {
3596  tt *= factor;
3597  return tt;
3598 }
3599 
3600 
3601 
3608 template <int rank, int dim>
3609 inline
3611 operator * (const double factor,
3612  const SymmetricTensor<rank,dim> &t)
3613 {
3615  tt *= factor;
3616  return tt;
3617 }
3618 
3619 
3620 
3626 template <int rank, int dim>
3627 inline
3629 operator / (const SymmetricTensor<rank,dim> &t,
3630  const double factor)
3631 {
3633  tt /= factor;
3634  return tt;
3635 }
3636 
3646 template <int dim, typename Number, typename OtherNumber>
3647 inline
3648 typename ProductType<Number, OtherNumber>::type
3651 {
3652  return (t1*t2);
3653 }
3654 
3655 
3665 template <int dim, typename Number, typename OtherNumber>
3666 inline
3667 typename ProductType<Number, OtherNumber>::type
3669  const Tensor<2,dim,OtherNumber> &t2)
3670 {
3671  Number s = 0;
3672  for (unsigned int i=0; i<dim; ++i)
3673  for (unsigned int j=0; j<dim; ++j)
3674  s += t1[i][j] * t2[i][j];
3675  return s;
3676 }
3677 
3678 
3688 template <int dim, typename Number, typename OtherNumber>
3689 inline
3690 typename ProductType<Number, OtherNumber>::type
3693 {
3694  return scalar_product(t2, t1);
3695 }
3696 
3697 
3713 template <typename Number, typename OtherNumber>
3714 inline
3715 void
3716 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3717  const SymmetricTensor<4,1,Number> &t,
3719 {
3720  tmp[0][0] = t[0][0][0][0] * s[0][0];
3721 }
3722 
3723 
3724 
3740 template <typename Number, typename OtherNumber>
3741 inline
3742 void
3743 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3744  const SymmetricTensor<2,1,Number> &s,
3746 {
3747  tmp[0][0] = t[0][0][0][0] * s[0][0];
3748 }
3749 
3750 
3751 
3766 template <typename Number, typename OtherNumber>
3767 inline
3768 void
3769 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3770  const SymmetricTensor<4,2,Number> &t,
3772 {
3773  const unsigned int dim = 2;
3774 
3775  for (unsigned int i=0; i<dim; ++i)
3776  for (unsigned int j=i; j<dim; ++j)
3777  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3778  t[i][j][1][1] * s[1][1] +
3779  2 * t[i][j][0][1] * s[0][1];
3780 }
3781 
3782 
3783 
3799 template <typename Number, typename OtherNumber>
3800 inline
3801 void
3802 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3803  const SymmetricTensor<2,2,Number> &s,
3805 {
3806  const unsigned int dim = 2;
3807 
3808  for (unsigned int i=0; i<dim; ++i)
3809  for (unsigned int j=i; j<dim; ++j)
3810  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3811  s[1][1] * t[1][1][i][j] +
3812  2 * s[0][1] * t[0][1][i][j];
3813 }
3814 
3815 
3816 
3832 template <typename Number, typename OtherNumber>
3833 inline
3834 void
3835 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3836  const SymmetricTensor<4,3,Number> &t,
3838 {
3839  const unsigned int dim = 3;
3840 
3841  for (unsigned int i=0; i<dim; ++i)
3842  for (unsigned int j=i; j<dim; ++j)
3843  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3844  t[i][j][1][1] * s[1][1] +
3845  t[i][j][2][2] * s[2][2] +
3846  2 * t[i][j][0][1] * s[0][1] +
3847  2 * t[i][j][0][2] * s[0][2] +
3848  2 * t[i][j][1][2] * s[1][2];
3849 }
3850 
3851 
3852 
3868 template <typename Number, typename OtherNumber>
3869 inline
3870 void
3871 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3872  const SymmetricTensor<2,3,Number> &s,
3874 {
3875  const unsigned int dim = 3;
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  s[2][2] * t[2][2][i][j] +
3882  2 * s[0][1] * t[0][1][i][j] +
3883  2 * s[0][2] * t[0][2][i][j] +
3884  2 * s[1][2] * t[1][2][i][j];
3885 }
3886 
3887 
3888 
3896 template <int dim, typename Number, typename OtherNumber>
3898 operator * (const SymmetricTensor<2,dim,Number> &src1,
3899  const Tensor<1,dim,OtherNumber> &src2)
3900 {
3902  for (unsigned int i=0; i<dim; ++i)
3903  for (unsigned int j=0; j<dim; ++j)
3904  dest[i] += src1[i][j] * src2[j];
3905  return dest;
3906 }
3907 
3908 
3916 template <int dim, typename Number, typename OtherNumber>
3918 operator * (const Tensor<1,dim,Number> &src1,
3920 {
3921  // this is easy for symmetric tensors:
3922  return src2 * src1;
3923 }
3924 
3925 
3926 
3947 template <int rank_1, int rank_2, int dim,
3948  typename Number, typename OtherNumber>
3949 inline DEAL_II_ALWAYS_INLINE
3950 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3951 operator * (const Tensor<rank_1, dim, Number> &src1,
3953 {
3954  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3955  const Tensor<rank_2, dim, OtherNumber> src2 (src2s);
3956  return src1*src2;
3957 }
3958 
3959 
3960 
3981 template <int rank_1, int rank_2, int dim,
3982  typename Number, typename OtherNumber>
3983 inline DEAL_II_ALWAYS_INLINE
3984 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3985 operator * (const SymmetricTensor<rank_1, dim, Number> &src1s,
3987 {
3988  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3989  const Tensor<rank_2, dim, OtherNumber> src1 (src1s);
3990  return src1*src2;
3991 }
3992 
3993 
3994 
4004 template <int dim, typename Number>
4005 inline
4006 std::ostream &operator << (std::ostream &out,
4008 {
4009  //make out lives a bit simpler by outputing
4010  //the tensor through the operator for the
4011  //general Tensor class
4013 
4014  for (unsigned int i=0; i<dim; ++i)
4015  for (unsigned int j=0; j<dim; ++j)
4016  tt[i][j] = t[i][j];
4017 
4018  return out << tt;
4019 }
4020 
4021 
4022 
4032 template <int dim, typename Number>
4033 inline
4034 std::ostream &operator << (std::ostream &out,
4036 {
4037  //make out lives a bit simpler by outputing
4038  //the tensor through the operator for the
4039  //general Tensor class
4041 
4042  for (unsigned int i=0; i<dim; ++i)
4043  for (unsigned int j=0; j<dim; ++j)
4044  for (unsigned int k=0; k<dim; ++k)
4045  for (unsigned int l=0; l<dim; ++l)
4046  tt[i][j][k][l] = t[i][j][k][l];
4047 
4048  return out << tt;
4049 }
4050 
4051 
4052 DEAL_II_NAMESPACE_CLOSE
4053 
4054 #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)
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:1207
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)
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)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
Definition: mpi.h:51
StreamType & operator<<(StreamType &s, const UpdateFlags u)
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()
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 > &)
std::array< std::pair< Number, Tensor< 1, 1, Number > >, 1 > eigenvectors(const SymmetricTensor< 2, 1, Number > &T, const SymmetricTensorEigenvectorMethod=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static::ExceptionBase & ExcInternalError()