Reference documentation for deal.II version Git 1dafe5d 2017-11-23 14:01:03 +0100
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 
550  static const unsigned int rank = rank_;
551 
557  static const unsigned int n_independent_components
560 
564  SymmetricTensor ();
565 
576  template <typename OtherNumber>
577  explicit
579 
595  SymmetricTensor (const Number (&array) [n_independent_components]);
596 
602  template <typename OtherNumber>
603  explicit
605 
611  template <typename OtherNumber>
613 
620  SymmetricTensor &operator = (const Number &d);
621 
626  operator Tensor<rank_,dim,Number> () const;
627 
631  bool operator == (const SymmetricTensor &) const;
632 
636  bool operator != (const SymmetricTensor &) const;
637 
641  template <typename OtherNumber>
643 
647  template <typename OtherNumber>
649 
654  template <typename OtherNumber>
655  SymmetricTensor &operator *= (const OtherNumber &factor);
656 
660  template <typename OtherNumber>
661  SymmetricTensor &operator /= (const OtherNumber &factor);
662 
667 
692  template <typename OtherNumber>
695 
700  template <typename OtherNumber>
703 
707  Number &operator() (const TableIndices<rank_> &indices);
708 
712  Number operator() (const TableIndices<rank_> &indices) const;
713 
718  internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
719  operator [] (const unsigned int row) const;
720 
725  internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
726  operator [] (const unsigned int row);
727 
733  Number
734  operator [] (const TableIndices<rank_> &indices) const;
735 
741  Number &
742  operator [] (const TableIndices<rank_> &indices);
743 
749  Number
750  access_raw_entry (const unsigned int unrolled_index) const;
751 
757  Number &
758  access_raw_entry (const unsigned int unrolled_index);
759 
770  norm () const;
771 
779  static
780  unsigned int
782 
788  static
790  unrolled_to_component_indices (const unsigned int i);
791 
804  void clear ();
805 
810  static std::size_t memory_consumption ();
811 
816  template <class Archive>
817  void serialize(Archive &ar, const unsigned int version);
818 
819 private:
823  typedef
826 
830  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
831 
835  base_tensor_type data;
836 
840  template <int, int, typename> friend class SymmetricTensor;
841 
845  template <int dim2, typename Number2>
846  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
847 
848  template <int dim2, typename Number2>
849  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
850 
851  template <int dim2, typename Number2>
853  deviator (const SymmetricTensor<2,dim2,Number2> &t);
854 
855  template <int dim2, typename Number2>
857 
858  template <int dim2, typename Number2>
860 
861  template <int dim2, typename Number2>
863 
864 
868  friend struct internal::SymmetricTensor::Inverse<2,dim,Number>;
869 
870  friend struct internal::SymmetricTensor::Inverse<4,dim,Number>;
871 };
872 
873 
874 
875 // ------------------------- inline functions ------------------------
876 
877 #ifndef DOXYGEN
878 
879 namespace internal
880 {
881  namespace SymmetricTensorAccessors
882  {
883  template <int rank_, int dim, bool constness, int P, typename Number>
884  Accessor<rank_,dim,constness,P,Number>::
885  Accessor (tensor_type &tensor,
886  const TableIndices<rank_> &previous_indices)
887  :
888  tensor (tensor),
889  previous_indices (previous_indices)
890  {}
891 
892 
893 
894  template <int rank_, int dim, bool constness, int P, typename Number>
895  Accessor<rank_,dim,constness,P-1,Number>
896  Accessor<rank_,dim,constness,P,Number>::operator[] (const unsigned int i)
897  {
898  return Accessor<rank_,dim,constness,P-1,Number> (tensor,
899  merge (previous_indices, i, rank_-P));
900  }
901 
902 
903 
904  template <int rank_, int dim, bool constness, int P, typename Number>
905  Accessor<rank_,dim,constness,P-1,Number>
906  Accessor<rank_,dim,constness,P,Number>::operator[] (const unsigned int i) const
907  {
908  return Accessor<rank_,dim,constness,P-1,Number> (tensor,
909  merge (previous_indices, i, rank_-P));
910  }
911 
912 
913 
914  template <int rank_, int dim, bool constness, typename Number>
915  Accessor<rank_,dim,constness,1,Number>::
916  Accessor (tensor_type &tensor,
917  const TableIndices<rank_> &previous_indices)
918  :
919  tensor (tensor),
920  previous_indices (previous_indices)
921  {}
922 
923 
924 
925  template <int rank_, int dim, bool constness, typename Number>
926  typename Accessor<rank_,dim,constness,1,Number>::reference
927  Accessor<rank_,dim,constness,1,Number>::operator[] (const unsigned int i)
928  {
929  return tensor(merge (previous_indices, i, rank_-1));
930  }
931 
932 
933  template <int rank_, int dim, bool constness, typename Number>
934  typename Accessor<rank_,dim,constness,1,Number>::reference
935  Accessor<rank_,dim,constness,1,Number>::operator[] (const unsigned int i) const
936  {
937  return tensor(merge (previous_indices, i, rank_-1));
938  }
939  }
940 }
941 
942 
943 
944 template <int rank_, int dim, typename Number>
945 inline
947 {
948  // Some auto-differentiable numbers need explicit
949  // zero initialization.
950  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
952 }
953 
954 
955 template <int rank_, int dim, typename Number>
956 template <typename OtherNumber>
957 inline
959 {
960  Assert (rank == 2, ExcNotImplemented());
961  switch (dim)
962  {
963  case 2:
964  Assert (t[0][1] == t[1][0], ExcInternalError());
965 
966  data[0] = t[0][0];
967  data[1] = t[1][1];
968  data[2] = t[0][1];
969 
970  break;
971  case 3:
972  Assert (t[0][1] == t[1][0], ExcInternalError());
973  Assert (t[0][2] == t[2][0], ExcInternalError());
974  Assert (t[1][2] == t[2][1], ExcInternalError());
975 
976  data[0] = t[0][0];
977  data[1] = t[1][1];
978  data[2] = t[2][2];
979  data[3] = t[0][1];
980  data[4] = t[0][2];
981  data[5] = t[1][2];
982 
983  break;
984  default:
985  for (unsigned int d=0; d<dim; ++d)
986  for (unsigned int e=0; e<d; ++e)
987  Assert(t[d][e] == t[e][d], ExcInternalError());
988 
989  for (unsigned int d=0; d<dim; ++d)
990  data[d] = t[d][d];
991 
992  for (unsigned int d=0, c=0; d<dim; ++d)
993  for (unsigned int e=d+1; e<dim; ++e, ++c)
994  data[dim+c] = t[d][e];
995  }
996 }
997 
998 
999 
1000 template <int rank_, int dim, typename Number>
1001 template <typename OtherNumber>
1002 inline
1005 {
1006  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1008 }
1009 
1010 
1011 
1012 
1013 template <int rank_, int dim, typename Number>
1014 inline
1015 SymmetricTensor<rank_,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
1016  :
1017  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
1018 {
1019  // ensure that the reinterpret_cast above actually works
1020  Assert (sizeof(typename base_tensor_type::array_type)
1021  == sizeof(array),
1022  ExcInternalError());
1023 }
1024 
1025 
1026 
1027 template <int rank_, int dim, typename Number>
1028 template <typename OtherNumber>
1029 inline
1032 {
1033  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
1034  data[i] = t.data[i];
1035  return *this;
1036 }
1037 
1038 
1039 
1040 template <int rank_, int dim, typename Number>
1041 inline
1044 {
1045  Assert (d==internal::NumberType<Number>::value(0.0), ExcMessage ("Only assignment with zero is allowed"));
1046  (void) d;
1047 
1049 
1050  return *this;
1051 }
1052 
1053 
1054 namespace internal
1055 {
1056  namespace SymmetricTensor
1057  {
1058  template <int dim, typename Number>
1060  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1061  {
1063 
1064  // diagonal entries are stored first
1065  for (unsigned int d=0; d<dim; ++d)
1066  t[d][d] = s.access_raw_entry(d);
1067 
1068  // off-diagonal entries come next, row by row
1069  for (unsigned int d=0, c=0; d<dim; ++d)
1070  for (unsigned int e=d+1; e<dim; ++e, ++c)
1071  {
1072  t[d][e] = s.access_raw_entry(dim+c);
1073  t[e][d] = s.access_raw_entry(dim+c);
1074  }
1075  return t;
1076  }
1077 
1078 
1079  template <int dim, typename Number>
1081  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1082  {
1083  // utilize the symmetry properties of SymmetricTensor<4,dim>
1084  // discussed in the class documentation to avoid accessing all
1085  // independent elements of the input tensor more than once
1087 
1088  for (unsigned int i=0; i<dim; ++i)
1089  for (unsigned int j=i; j<dim; ++j)
1090  for (unsigned int k=0; k<dim; ++k)
1091  for (unsigned int l=k; l<dim; ++l)
1092  t[TableIndices<4>(i,j,k,l)]
1093  = t[TableIndices<4>(i,j,l,k)]
1094  = t[TableIndices<4>(j,i,k,l)]
1095  = t[TableIndices<4>(j,i,l,k)]
1096  = st[TableIndices<4>(i,j,k,l)];
1097 
1098  return t;
1099  }
1100 
1101 
1102  template <typename Number>
1103  struct Inverse<2,1,Number>
1104  {
1105  static inline ::SymmetricTensor<2,1,Number>
1106  value (const ::SymmetricTensor<2,1,Number> &t)
1107  {
1109 
1110  tmp[0][0] = 1.0/t[0][0];
1111 
1112  return tmp;
1113  }
1114  };
1115 
1116 
1117  template <typename Number>
1118  struct Inverse<2,2,Number>
1119  {
1120  static inline ::SymmetricTensor<2,2,Number>
1121  value (const ::SymmetricTensor<2,2,Number> &t)
1122  {
1124 
1125  // Sympy result: ([
1126  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
1127  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
1128  const TableIndices<2> idx_00 (0,0);
1129  const TableIndices<2> idx_01 (0,1);
1130  const TableIndices<2> idx_11 (1,1);
1131  const Number inv_det_t
1132  = 1.0/(t[idx_00]*t[idx_11]
1133  - t[idx_01]*t[idx_01]);
1134  tmp[idx_00] = t[idx_11];
1135  tmp[idx_01] = -t[idx_01];
1136  tmp[idx_11] = t[idx_00];
1137  tmp *= inv_det_t;
1138 
1139  return tmp;
1140  }
1141  };
1142 
1143 
1144  template <typename Number>
1145  struct Inverse<2,3,Number>
1146  {
1147  static ::SymmetricTensor<2,3,Number>
1148  value (const ::SymmetricTensor<2,3,Number> &t)
1149  {
1151 
1152  // Sympy result: ([
1153  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1154  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1155  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
1156  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1157  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1158  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
1159  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
1160  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
1161  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
1162  const TableIndices<2> idx_00 (0,0);
1163  const TableIndices<2> idx_01 (0,1);
1164  const TableIndices<2> idx_02 (0,2);
1165  const TableIndices<2> idx_11 (1,1);
1166  const TableIndices<2> idx_12 (1,2);
1167  const TableIndices<2> idx_22 (2,2);
1168  const Number inv_det_t
1169  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
1170  - t[idx_00]*t[idx_12]*t[idx_12]
1171  - t[idx_01]*t[idx_01]*t[idx_22]
1172  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
1173  - t[idx_02]*t[idx_02]*t[idx_11]);
1174  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
1175  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
1176  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
1177  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
1178  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
1179  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
1180  tmp *= inv_det_t;
1181 
1182  return tmp;
1183  }
1184  };
1185 
1186 
1187  template <typename Number>
1188  struct Inverse<4,1,Number>
1189  {
1190  static inline ::SymmetricTensor<4,1,Number>
1191  value (const ::SymmetricTensor<4,1,Number> &t)
1192  {
1194  tmp.data[0][0] = 1.0/t.data[0][0];
1195  return tmp;
1196  }
1197  };
1198 
1199 
1200  template <typename Number>
1201  struct Inverse<4,2,Number>
1202  {
1203  static inline ::SymmetricTensor<4,2,Number>
1204  value (const ::SymmetricTensor<4,2,Number> &t)
1205  {
1207 
1208  // Inverting this tensor is a little more complicated than necessary,
1209  // since we store the data of 't' as a 3x3 matrix t.data, but the
1210  // product between a rank-4 and a rank-2 tensor is really not the
1211  // product between this matrix and the 3-vector of a rhs, but rather
1212  //
1213  // B.vec = t.data * mult * A.vec
1214  //
1215  // where mult is a 3x3 matrix with entries [[1,0,0],[0,1,0],[0,0,2]] to
1216  // capture the fact that we need to add up both the c_ij12*a_12 and the
1217  // c_ij21*a_21 terms.
1218  //
1219  // In addition, in this scheme, the identity tensor has the matrix
1220  // representation mult^-1.
1221  //
1222  // The inverse of 't' therefore has the matrix representation
1223  //
1224  // inv.data = mult^-1 * t.data^-1 * mult^-1
1225  //
1226  // in order to compute it, let's first compute the inverse of t.data and
1227  // put it into tmp.data; at the end of the function we then scale the
1228  // last row and column of the inverse by 1/2, corresponding to the left
1229  // and right multiplication with mult^-1.
1230  const Number t4 = t.data[0][0]*t.data[1][1],
1231  t6 = t.data[0][0]*t.data[1][2],
1232  t8 = t.data[0][1]*t.data[1][0],
1233  t00 = t.data[0][2]*t.data[1][0],
1234  t01 = t.data[0][1]*t.data[2][0],
1235  t04 = t.data[0][2]*t.data[2][0],
1236  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
1237  t8*t.data[2][2]+t00*t.data[2][1]+
1238  t01*t.data[1][2]-t04*t.data[1][1]);
1239  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
1240  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
1241  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
1242  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
1243  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
1244  tmp.data[1][2] = -(t6-t00)*t07;
1245  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
1246  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
1247  tmp.data[2][2] = (t4-t8)*t07;
1248 
1249  // scale last row and column as mentioned
1250  // above
1251  tmp.data[2][0] /= 2;
1252  tmp.data[2][1] /= 2;
1253  tmp.data[0][2] /= 2;
1254  tmp.data[1][2] /= 2;
1255  tmp.data[2][2] /= 4;
1256 
1257  return tmp;
1258  }
1259  };
1260 
1261 
1262  template <typename Number>
1263  struct Inverse<4,3,Number>
1264  {
1265  static ::SymmetricTensor<4,3,Number>
1266  value (const ::SymmetricTensor<4,3,Number> &t)
1267  {
1269 
1270  // This function follows the exact same scheme as the 2d case, except
1271  // that hardcoding the inverse of a 6x6 matrix is pretty wasteful.
1272  // Instead, we use the Gauss-Jordan algorithm implemented for
1273  // FullMatrix. For historical reasons the following code is copied from
1274  // there, with the tangential benefit that we do not need to copy the
1275  // tensor entries to and from the FullMatrix.
1276  const unsigned int N = 6;
1277 
1278  // First get an estimate of the size of the elements of this matrix,
1279  // for later checks whether the pivot element is large enough, or
1280  // whether we have to fear that the matrix is not regular.
1281  Number diagonal_sum = internal::NumberType<Number>::value(0.0);
1282  for (unsigned int i=0; i<N; ++i)
1283  diagonal_sum += std::fabs(tmp.data[i][i]);
1284  const Number typical_diagonal_element = diagonal_sum/N;
1285  (void)typical_diagonal_element;
1286 
1287  unsigned int p[N];
1288  for (unsigned int i=0; i<N; ++i)
1289  p[i] = i;
1290 
1291  for (unsigned int j=0; j<N; ++j)
1292  {
1293  // Pivot search: search that part of the line on and right of the
1294  // diagonal for the largest element.
1295  Number max = std::fabs(tmp.data[j][j]);
1296  unsigned int r = j;
1297  for (unsigned int i=j+1; i<N; ++i)
1298  if (std::fabs(tmp.data[i][j]) > max)
1299  {
1300  max = std::fabs(tmp.data[i][j]);
1301  r = i;
1302  }
1303 
1304  // Check whether the pivot is too small
1305  Assert(max > 1.e-16*typical_diagonal_element,
1306  ExcMessage("This tensor seems to be noninvertible"));
1307 
1308  // Row interchange
1309  if (r>j)
1310  {
1311  for (unsigned int k=0; k<N; ++k)
1312  std::swap (tmp.data[j][k], tmp.data[r][k]);
1313 
1314  std::swap (p[j], p[r]);
1315  }
1316 
1317  // Transformation
1318  const Number hr = 1./tmp.data[j][j];
1319  tmp.data[j][j] = hr;
1320  for (unsigned int k=0; k<N; ++k)
1321  {
1322  if (k==j) continue;
1323  for (unsigned int i=0; i<N; ++i)
1324  {
1325  if (i==j) continue;
1326  tmp.data[i][k] -= tmp.data[i][j]*tmp.data[j][k]*hr;
1327  }
1328  }
1329  for (unsigned int i=0; i<N; ++i)
1330  {
1331  tmp.data[i][j] *= hr;
1332  tmp.data[j][i] *= -hr;
1333  }
1334  tmp.data[j][j] = hr;
1335  }
1336 
1337  // Column interchange
1338  Number hv[N];
1339  for (unsigned int i=0; i<N; ++i)
1340  {
1341  for (unsigned int k=0; k<N; ++k)
1342  hv[p[k]] = tmp.data[i][k];
1343  for (unsigned int k=0; k<N; ++k)
1344  tmp.data[i][k] = hv[k];
1345  }
1346 
1347  // Scale rows and columns. The mult matrix
1348  // here is diag[1, 1, 1, 1/2, 1/2, 1/2].
1349  for (unsigned int i=3; i<6; ++i)
1350  for (unsigned int j=0; j<3; ++j)
1351  tmp.data[i][j] /= 2;
1352 
1353  for (unsigned int i=0; i<3; ++i)
1354  for (unsigned int j=3; j<6; ++j)
1355  tmp.data[i][j] /= 2;
1356 
1357  for (unsigned int i=3; i<6; ++i)
1358  for (unsigned int j=3; j<6; ++j)
1359  tmp.data[i][j] /= 4;
1360 
1361  return tmp;
1362  }
1363  };
1364 
1365  }
1366 }
1367 
1368 
1369 
1370 template <int rank_, int dim, typename Number>
1371 inline
1373 operator Tensor<rank_,dim,Number> () const
1374 {
1375  return internal::SymmetricTensor::convert_to_tensor (*this);
1376 }
1377 
1378 
1379 
1380 template <int rank_, int dim, typename Number>
1381 inline
1382 bool
1384 (const SymmetricTensor<rank_,dim,Number> &t) const
1385 {
1386  return data == t.data;
1387 }
1388 
1389 
1390 
1391 template <int rank_, int dim, typename Number>
1392 inline
1393 bool
1394 SymmetricTensor<rank_,dim,Number>::operator !=
1395 (const SymmetricTensor<rank_,dim,Number> &t) const
1396 {
1397  return data != t.data;
1398 }
1399 
1400 
1401 
1402 template <int rank_, int dim, typename Number>
1403 template <typename OtherNumber>
1404 inline
1406 SymmetricTensor<rank_,dim,Number>::operator +=
1408 {
1409  data += t.data;
1410  return *this;
1411 }
1412 
1413 
1414 
1415 template <int rank_, int dim, typename Number>
1416 template <typename OtherNumber>
1417 inline
1419 SymmetricTensor<rank_,dim,Number>::operator -=
1421 {
1422  data -= t.data;
1423  return *this;
1424 }
1425 
1426 
1427 
1428 template <int rank_, int dim, typename Number>
1429 template <typename OtherNumber>
1430 inline
1433 {
1434  data *= d;
1435  return *this;
1436 }
1437 
1438 
1439 
1440 template <int rank_, int dim, typename Number>
1441 template <typename OtherNumber>
1442 inline
1445 {
1446  data /= d;
1447  return *this;
1448 }
1449 
1450 
1451 
1452 template <int rank_, int dim, typename Number>
1453 inline
1456 {
1457  SymmetricTensor tmp = *this;
1458  tmp.data = -tmp.data;
1459  return tmp;
1460 }
1461 
1462 
1463 
1464 template <int rank_, int dim, typename Number>
1465 inline
1466 void
1468 {
1469  data.clear ();
1470 }
1471 
1472 
1473 
1474 template <int rank_, int dim, typename Number>
1475 inline
1476 std::size_t
1478 {
1479  // all memory consists of statically allocated memory of the current
1480  // object, no pointers
1481  return sizeof(SymmetricTensor<rank_,dim,Number>);
1482 }
1483 
1484 
1485 
1486 namespace internal
1487 {
1488 
1489  template <int dim, typename Number, typename OtherNumber = Number>
1490  inline
1491  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type
1492  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1493  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1494  {
1495  typedef typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number,OtherNumber>::type result_type;
1496 
1497  switch (dim)
1498  {
1499  case 1:
1500  return data[0] * sdata[0];
1501  default:
1502  // Start with the non-diagonal part to avoid some multiplications by
1503  // 2.
1504 
1505  result_type sum = data[dim] * sdata[dim];
1506  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1507  sum += data[d] * sdata[d];
1508  sum += sum; // sum = sum * 2.;
1509 
1510  // Now add the contributions from the diagonal
1511  for (unsigned int d=0; d<dim; ++d)
1512  sum += data[d] * sdata[d];
1513  return sum;
1514  }
1515  }
1516 
1517 
1518 
1519  template <int dim, typename Number, typename OtherNumber = Number>
1520  inline
1521  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type
1522  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1523  const typename SymmetricTensorAccessors::StorageType<2,dim,OtherNumber>::base_tensor_type &sdata)
1524  {
1525  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::type result_type;
1526  typedef typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number,OtherNumber>::value_type value_type;
1527 
1528  const unsigned int data_dim =
1529  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1530  value_type tmp [data_dim];
1531  for (unsigned int i=0; i<data_dim; ++i)
1532  tmp[i] = perform_double_contraction<dim,Number,OtherNumber>(data[i], sdata);
1533  return result_type(tmp);
1534  }
1535 
1536 
1537 
1538  template <int dim, typename Number, typename OtherNumber = Number>
1539  inline
1540  typename SymmetricTensorAccessors::StorageType<2,dim,
1541  typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type
1542  >::base_tensor_type
1543  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1544  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1545  {
1546  typedef typename SymmetricTensorAccessors::double_contraction_result<2,4,dim,Number,OtherNumber>::value_type value_type;
1547  typedef typename SymmetricTensorAccessors::StorageType<2,dim,value_type>::base_tensor_type base_tensor_type;
1548 
1549  base_tensor_type tmp;
1550  for (unsigned int i=0; i<tmp.dimension; ++i)
1551  {
1552  // Start with the non-diagonal part
1553  value_type sum = data[dim] * sdata[dim][i];
1554  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1555  sum += data[d] * sdata[d][i];
1556  sum += sum; // sum = sum * 2.;
1557 
1558  // Now add the contributions from the diagonal
1559  for (unsigned int d=0; d<dim; ++d)
1560  sum += data[d] * sdata[d][i];
1561  tmp[i] = sum;
1562  }
1563  return tmp;
1564  }
1565 
1566 
1567 
1568  template <int dim, typename Number, typename OtherNumber = Number>
1569  inline
1570  typename SymmetricTensorAccessors::StorageType<4,dim,
1571  typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type
1572  >::base_tensor_type
1573  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1574  const typename SymmetricTensorAccessors::StorageType<4,dim,OtherNumber>::base_tensor_type &sdata)
1575  {
1576  typedef typename SymmetricTensorAccessors::double_contraction_result<4,4,dim,Number,OtherNumber>::value_type value_type;
1577  typedef typename SymmetricTensorAccessors::StorageType<4,dim,value_type>::base_tensor_type base_tensor_type;
1578 
1579  const unsigned int data_dim =
1580  SymmetricTensorAccessors::StorageType<2,dim,value_type>::n_independent_components;
1581  base_tensor_type tmp;
1582  for (unsigned int i=0; i<data_dim; ++i)
1583  for (unsigned int j=0; j<data_dim; ++j)
1584  {
1585  // Start with the non-diagonal part
1586  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1587  tmp[i][j] += data[i][d] * sdata[d][j];
1588  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1589 
1590  // Now add the contributions from the diagonal
1591  for (unsigned int d=0; d<dim; ++d)
1592  tmp[i][j] += data[i][d] * sdata[d][j];
1593  }
1594  return tmp;
1595  }
1596 
1597 } // end of namespace internal
1598 
1599 
1600 
1601 template <int rank_, int dim, typename Number>
1602 template <typename OtherNumber>
1603 inline
1606 {
1607  // need to have two different function calls
1608  // because a scalar and rank-2 tensor are not
1609  // the same data type (see internal function
1610  // above)
1611  return internal::perform_double_contraction<dim,Number,OtherNumber> (data, s.data);
1612 }
1613 
1614 
1615 
1616 template <int rank_, int dim, typename Number>
1617 template <typename OtherNumber>
1618 inline
1621 {
1624  tmp.data = internal::perform_double_contraction<dim,Number,OtherNumber> (data,s.data);
1625  return tmp;
1626 }
1627 
1628 
1629 
1630 // internal namespace to switch between the
1631 // access of different tensors. There used to
1632 // be explicit instantiations before for
1633 // different ranks and dimensions, but since
1634 // we now allow for templates on the data
1635 // type, and since we cannot partially
1636 // specialize the implementation, this got
1637 // into a separate namespace
1638 namespace internal
1639 {
1640  template <int dim, typename Number>
1641  inline
1642  Number &
1643  symmetric_tensor_access (const TableIndices<2> &indices,
1645  {
1646  // 1d is very simple and done first
1647  if (dim == 1)
1648  return data[0];
1649 
1650  // first treat the main diagonal elements, which are stored consecutively
1651  // at the beginning
1652  if (indices[0] == indices[1])
1653  return data[indices[0]];
1654 
1655  // the rest is messier and requires a few switches.
1656  switch (dim)
1657  {
1658  case 2:
1659  // at least for the 2x2 case it is reasonably simple
1660  Assert (((indices[0]==1) && (indices[1]==0)) ||
1661  ((indices[0]==0) && (indices[1]==1)),
1662  ExcInternalError());
1663  return data[2];
1664 
1665  default:
1666  // to do the rest, sort our indices before comparing
1667  {
1668  TableIndices<2> sorted_indices (indices);
1669  sorted_indices.sort ();
1670 
1671  for (unsigned int d=0, c=0; d<dim; ++d)
1672  for (unsigned int e=d+1; e<dim; ++e, ++c)
1673  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1674  return data[dim+c];
1675  Assert (false, ExcInternalError());
1676  }
1677  }
1678 
1679  static Number dummy_but_referenceable = Number();
1680  return dummy_but_referenceable;
1681  }
1682 
1683 
1684 
1685  template <int dim, typename Number>
1686  inline
1687  Number
1688  symmetric_tensor_access (const TableIndices<2> &indices,
1690  {
1691  // 1d is very simple and done first
1692  if (dim == 1)
1693  return data[0];
1694 
1695  // first treat the main diagonal elements, which are stored consecutively
1696  // at the beginning
1697  if (indices[0] == indices[1])
1698  return data[indices[0]];
1699 
1700  // the rest is messier and requires a few switches.
1701  switch (dim)
1702  {
1703  case 2:
1704  // at least for the 2x2 case it is reasonably simple
1705  Assert (((indices[0]==1) && (indices[1]==0)) ||
1706  ((indices[0]==0) && (indices[1]==1)),
1707  ExcInternalError());
1708  return data[2];
1709 
1710  default:
1711  // to do the rest, sort our indices before comparing
1712  {
1713  TableIndices<2> sorted_indices (indices);
1714  sorted_indices.sort ();
1715 
1716  for (unsigned int d=0, c=0; d<dim; ++d)
1717  for (unsigned int e=d+1; e<dim; ++e, ++c)
1718  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1719  return data[dim+c];
1720  Assert (false, ExcInternalError());
1721  }
1722  }
1723 
1724  static Number dummy_but_referenceable = Number();
1725  return dummy_but_referenceable;
1726  }
1727 
1728 
1729 
1730  template <int dim, typename Number>
1731  inline
1732  Number &
1733  symmetric_tensor_access (const TableIndices<4> &indices,
1735  {
1736  switch (dim)
1737  {
1738  case 1:
1739  return data[0][0];
1740 
1741  case 2:
1742  // each entry of the tensor can be
1743  // thought of as an entry in a
1744  // matrix that maps the rolled-out
1745  // rank-2 tensors into rolled-out
1746  // rank-2 tensors. this is the
1747  // format in which we store rank-4
1748  // tensors. determine which
1749  // position the present entry is
1750  // stored in
1751  {
1752  unsigned int base_index[2] ;
1753  if ((indices[0] == 0) && (indices[1] == 0))
1754  base_index[0] = 0;
1755  else if ((indices[0] == 1) && (indices[1] == 1))
1756  base_index[0] = 1;
1757  else
1758  base_index[0] = 2;
1759 
1760  if ((indices[2] == 0) && (indices[3] == 0))
1761  base_index[1] = 0;
1762  else if ((indices[2] == 1) && (indices[3] == 1))
1763  base_index[1] = 1;
1764  else
1765  base_index[1] = 2;
1766 
1767  return data[base_index[0]][base_index[1]];
1768  }
1769 
1770  case 3:
1771  // each entry of the tensor can be
1772  // thought of as an entry in a
1773  // matrix that maps the rolled-out
1774  // rank-2 tensors into rolled-out
1775  // rank-2 tensors. this is the
1776  // format in which we store rank-4
1777  // tensors. determine which
1778  // position the present entry is
1779  // stored in
1780  {
1781  unsigned int base_index[2] ;
1782  if ((indices[0] == 0) && (indices[1] == 0))
1783  base_index[0] = 0;
1784  else if ((indices[0] == 1) && (indices[1] == 1))
1785  base_index[0] = 1;
1786  else if ((indices[0] == 2) && (indices[1] == 2))
1787  base_index[0] = 2;
1788  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1789  ((indices[0] == 1) && (indices[1] == 0)))
1790  base_index[0] = 3;
1791  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1792  ((indices[0] == 2) && (indices[1] == 0)))
1793  base_index[0] = 4;
1794  else
1795  {
1796  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1797  ((indices[0] == 2) && (indices[1] == 1)),
1798  ExcInternalError());
1799  base_index[0] = 5;
1800  }
1801 
1802  if ((indices[2] == 0) && (indices[3] == 0))
1803  base_index[1] = 0;
1804  else if ((indices[2] == 1) && (indices[3] == 1))
1805  base_index[1] = 1;
1806  else if ((indices[2] == 2) && (indices[3] == 2))
1807  base_index[1] = 2;
1808  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1809  ((indices[2] == 1) && (indices[3] == 0)))
1810  base_index[1] = 3;
1811  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1812  ((indices[2] == 2) && (indices[3] == 0)))
1813  base_index[1] = 4;
1814  else
1815  {
1816  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1817  ((indices[2] == 2) && (indices[3] == 1)),
1818  ExcInternalError());
1819  base_index[1] = 5;
1820  }
1821 
1822  return data[base_index[0]][base_index[1]];
1823  }
1824 
1825  default:
1826  Assert (false, ExcNotImplemented());
1827  }
1828 
1829  static Number dummy;
1830  return dummy;
1831  }
1832 
1833 
1834  template <int dim, typename Number>
1835  inline
1836  Number
1837  symmetric_tensor_access (const TableIndices<4> &indices,
1839  {
1840  switch (dim)
1841  {
1842  case 1:
1843  return data[0][0];
1844 
1845  case 2:
1846  // each entry of the tensor can be
1847  // thought of as an entry in a
1848  // matrix that maps the rolled-out
1849  // rank-2 tensors into rolled-out
1850  // rank-2 tensors. this is the
1851  // format in which we store rank-4
1852  // tensors. determine which
1853  // position the present entry is
1854  // stored in
1855  {
1856  unsigned int base_index[2] ;
1857  if ((indices[0] == 0) && (indices[1] == 0))
1858  base_index[0] = 0;
1859  else if ((indices[0] == 1) && (indices[1] == 1))
1860  base_index[0] = 1;
1861  else
1862  base_index[0] = 2;
1863 
1864  if ((indices[2] == 0) && (indices[3] == 0))
1865  base_index[1] = 0;
1866  else if ((indices[2] == 1) && (indices[3] == 1))
1867  base_index[1] = 1;
1868  else
1869  base_index[1] = 2;
1870 
1871  return data[base_index[0]][base_index[1]];
1872  }
1873 
1874  case 3:
1875  // each entry of the tensor can be
1876  // thought of as an entry in a
1877  // matrix that maps the rolled-out
1878  // rank-2 tensors into rolled-out
1879  // rank-2 tensors. this is the
1880  // format in which we store rank-4
1881  // tensors. determine which
1882  // position the present entry is
1883  // stored in
1884  {
1885  unsigned int base_index[2] ;
1886  if ((indices[0] == 0) && (indices[1] == 0))
1887  base_index[0] = 0;
1888  else if ((indices[0] == 1) && (indices[1] == 1))
1889  base_index[0] = 1;
1890  else if ((indices[0] == 2) && (indices[1] == 2))
1891  base_index[0] = 2;
1892  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1893  ((indices[0] == 1) && (indices[1] == 0)))
1894  base_index[0] = 3;
1895  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1896  ((indices[0] == 2) && (indices[1] == 0)))
1897  base_index[0] = 4;
1898  else
1899  {
1900  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1901  ((indices[0] == 2) && (indices[1] == 1)),
1902  ExcInternalError());
1903  base_index[0] = 5;
1904  }
1905 
1906  if ((indices[2] == 0) && (indices[3] == 0))
1907  base_index[1] = 0;
1908  else if ((indices[2] == 1) && (indices[3] == 1))
1909  base_index[1] = 1;
1910  else if ((indices[2] == 2) && (indices[3] == 2))
1911  base_index[1] = 2;
1912  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1913  ((indices[2] == 1) && (indices[3] == 0)))
1914  base_index[1] = 3;
1915  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1916  ((indices[2] == 2) && (indices[3] == 0)))
1917  base_index[1] = 4;
1918  else
1919  {
1920  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1921  ((indices[2] == 2) && (indices[3] == 1)),
1922  ExcInternalError());
1923  base_index[1] = 5;
1924  }
1925 
1926  return data[base_index[0]][base_index[1]];
1927  }
1928 
1929  default:
1930  Assert (false, ExcNotImplemented());
1931  }
1932 
1933  static Number dummy;
1934  return dummy;
1935  }
1936 
1937 } // end of namespace internal
1938 
1939 
1940 
1941 template <int rank_, int dim, typename Number>
1942 inline
1943 Number &
1945 {
1946  for (unsigned int r=0; r<rank; ++r)
1947  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1948  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1949 }
1950 
1951 
1952 
1953 template <int rank_, int dim, typename Number>
1954 inline
1955 Number
1957 (const TableIndices<rank_> &indices) const
1958 {
1959  for (unsigned int r=0; r<rank; ++r)
1960  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1961  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1962 }
1963 
1964 
1965 
1966 namespace internal
1967 {
1968  namespace SymmetricTensor
1969  {
1970  template <int rank_>
1972  get_partially_filled_indices (const unsigned int row,
1973  const std::integral_constant<int, 2> &)
1974  {
1975  return TableIndices<rank_> (row,
1977 
1978  }
1979 
1980 
1981  template <int rank_>
1983  get_partially_filled_indices (const unsigned int row,
1984  const std::integral_constant<int, 4> &)
1985  {
1986  return TableIndices<rank_> (row,
1990 
1991  }
1992  }
1993 }
1994 
1995 
1996 template <int rank_, int dim, typename Number>
1997 internal::SymmetricTensorAccessors::Accessor<rank_,dim,true,rank_-1,Number>
1998 SymmetricTensor<rank_,dim,Number>::operator [] (const unsigned int row) const
1999 {
2000  return
2001  internal::SymmetricTensorAccessors::
2002  Accessor<rank_,dim,true,rank_-1,Number> (*this,
2003  internal::SymmetricTensor::get_partially_filled_indices<rank_> (row,
2004  std::integral_constant<int, rank_>()));
2005 }
2006 
2007 
2008 
2009 template <int rank_, int dim, typename Number>
2010 internal::SymmetricTensorAccessors::Accessor<rank_,dim,false,rank_-1,Number>
2012 {
2013  return
2014  internal::SymmetricTensorAccessors::
2015  Accessor<rank_,dim,false,rank_-1,Number> (*this,
2016  internal::SymmetricTensor::get_partially_filled_indices<rank_> (row,
2017  std::integral_constant<int, rank_>()));
2018 }
2019 
2020 
2021 
2022 template <int rank_, int dim, typename Number>
2023 inline
2024 Number
2026 {
2027  return operator()(indices);
2028 }
2029 
2030 
2031 
2032 template <int rank_, int dim, typename Number>
2033 inline
2034 Number &
2036 {
2037  return operator()(indices);
2038 }
2039 
2040 
2041 
2042 
2043 namespace internal
2044 {
2045  namespace SymmetricTensor
2046  {
2047  template <int dim, typename Number>
2048  unsigned int
2049  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
2050  const unsigned int index)
2051  {
2052  return index;
2053  }
2054 
2055 
2056  template <int dim, typename Number>
2058  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
2059  const unsigned int index)
2060  {
2061  return
2064  }
2065 
2066  }
2067 }
2068 
2069 
2070 
2071 template <int rank_, int dim, typename Number>
2072 inline
2073 Number
2074 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index) const
2075 {
2076  AssertIndexRange (index, n_independent_components);
2077  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2078 }
2079 
2080 
2081 
2082 template <int rank_, int dim, typename Number>
2083 inline
2084 Number &
2085 SymmetricTensor<rank_,dim,Number>::access_raw_entry (const unsigned int index)
2086 {
2087  AssertIndexRange (index, n_independent_components);
2088  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
2089 }
2090 
2091 
2092 
2093 namespace internal
2094 {
2095  template <int dim, typename Number>
2096  inline
2099  {
2100  switch (dim)
2101  {
2102  case 1:
2103  return numbers::NumberTraits<Number>::abs(data[0]);
2104 
2105  case 2:
2106  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2109 
2110  case 3:
2111  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
2117 
2118  default:
2119  {
2120  typename numbers::NumberTraits<Number>::real_type return_value
2122 
2123  for (unsigned int d=0; d<dim; ++d)
2124  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
2125  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
2126  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
2127 
2128  return std::sqrt(return_value);
2129  }
2130  }
2131  }
2132 
2133 
2134 
2135  template <int dim, typename Number>
2136  inline
2139  {
2140  switch (dim)
2141  {
2142  case 1:
2143  return numbers::NumberTraits<Number>::abs (data[0][0]);
2144 
2145  default:
2146  {
2147  typename numbers::NumberTraits<Number>::real_type return_value
2149 
2150  const unsigned int n_independent_components = data.dimension;
2151 
2152  for (unsigned int i=0; i<dim; ++i)
2153  for (unsigned int j=0; j<dim; ++j)
2154  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
2155  for (unsigned int i=0; i<dim; ++i)
2156  for (unsigned int j=dim; j<n_independent_components; ++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=0; j<dim; ++j)
2160  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2161  for (unsigned int i=dim; i<n_independent_components; ++i)
2162  for (unsigned int j=dim; j<n_independent_components; ++j)
2163  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
2164 
2165  return std::sqrt(return_value);
2166  }
2167  }
2168  }
2169 
2170 } // end of namespace internal
2171 
2172 
2173 
2174 template <int rank_, int dim, typename Number>
2175 inline
2178 {
2179  return internal::compute_norm<dim,Number> (data);
2180 }
2181 
2182 
2183 
2184 namespace internal
2185 {
2186  namespace SymmetricTensor
2187  {
2188  namespace
2189  {
2190  // a function to do the unrolling from a set of indices to a
2191  // scalar index into the array in which we store the elements of
2192  // a symmetric tensor
2193  //
2194  // this function is for rank-2 tensors
2195  template <int dim>
2196  inline
2197  unsigned int
2198  component_to_unrolled_index
2199  (const TableIndices<2> &indices)
2200  {
2201  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2202  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2203 
2204  switch (dim)
2205  {
2206  case 1:
2207  {
2208  return 0;
2209  }
2210 
2211  case 2:
2212  {
2213  static const unsigned int table[2][2] = {{0, 2},
2214  {2, 1}
2215  };
2216  return table[indices[0]][indices[1]];
2217  }
2218 
2219  case 3:
2220  {
2221  static const unsigned int table[3][3] = {{0, 3, 4},
2222  {3, 1, 5},
2223  {4, 5, 2}
2224  };
2225  return table[indices[0]][indices[1]];
2226  }
2227 
2228  case 4:
2229  {
2230  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2231  {4, 1, 7, 8},
2232  {5, 7, 2, 9},
2233  {6, 8, 9, 3}
2234  };
2235  return table[indices[0]][indices[1]];
2236  }
2237 
2238  default:
2239  // for the remainder, manually figure out the numbering
2240  {
2241  if (indices[0] == indices[1])
2242  return indices[0];
2243 
2244  TableIndices<2> sorted_indices (indices);
2245  sorted_indices.sort ();
2246 
2247  for (unsigned int d=0, c=0; d<dim; ++d)
2248  for (unsigned int e=d+1; e<dim; ++e, ++c)
2249  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
2250  return dim+c;
2251 
2252  // should never get here:
2253  Assert(false, ExcInternalError());
2254  return 0;
2255  }
2256  }
2257  }
2258 
2259  // a function to do the unrolling from a set of indices to a
2260  // scalar index into the array in which we store the elements of
2261  // a symmetric tensor
2262  //
2263  // this function is for tensors of ranks not already handled
2264  // above
2265  template <int dim, int rank_>
2266  inline
2267  unsigned int
2268  component_to_unrolled_index
2269  (const TableIndices<rank_> &indices)
2270  {
2271  (void)indices;
2272  Assert (false, ExcNotImplemented());
2274  }
2275  }
2276  }
2277 }
2278 
2279 
2280 template <int rank_, int dim, typename Number>
2281 inline
2282 unsigned int
2284 (const TableIndices<rank_> &indices)
2285 {
2286  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
2287 }
2288 
2289 
2290 
2291 namespace internal
2292 {
2293  namespace SymmetricTensor
2294  {
2295  namespace
2296  {
2297  // a function to do the inverse of the unrolling from a set of
2298  // indices to a scalar index into the array in which we store
2299  // the elements of a symmetric tensor. in other words, it goes
2300  // from the scalar index into the array to a set of indices of
2301  // the tensor
2302  //
2303  // this function is for rank-2 tensors
2304  template <int dim>
2305  inline
2307  unrolled_to_component_indices
2308  (const unsigned int i,
2309  const std::integral_constant<int, 2> &)
2310  {
2313  switch (dim)
2314  {
2315  case 1:
2316  {
2317  return TableIndices<2>(0,0);
2318  }
2319 
2320  case 2:
2321  {
2322  const TableIndices<2> table[3] =
2323  {
2324  TableIndices<2> (0,0),
2325  TableIndices<2> (1,1),
2326  TableIndices<2> (0,1)
2327  };
2328  return table[i];
2329  }
2330 
2331  case 3:
2332  {
2333  const TableIndices<2> table[6] =
2334  {
2335  TableIndices<2> (0,0),
2336  TableIndices<2> (1,1),
2337  TableIndices<2> (2,2),
2338  TableIndices<2> (0,1),
2339  TableIndices<2> (0,2),
2340  TableIndices<2> (1,2)
2341  };
2342  return table[i];
2343  }
2344 
2345  default:
2346  if (i<dim)
2347  return TableIndices<2> (i,i);
2348 
2349  for (unsigned int d=0, c=0; d<dim; ++d)
2350  for (unsigned int e=d+1; e<dim; ++e, ++c)
2351  if (c==i)
2352  return TableIndices<2>(d,e);
2353 
2354  // should never get here:
2355  Assert(false, ExcInternalError());
2356  return TableIndices<2>(0, 0);
2357  }
2358  }
2359 
2360  // a function to do the inverse of the unrolling from a set of
2361  // indices to a scalar index into the array in which we store
2362  // the elements of a symmetric tensor. in other words, it goes
2363  // from the scalar index into the array to a set of indices of
2364  // the tensor
2365  //
2366  // this function is for tensors of a rank not already handled
2367  // above
2368  template <int dim, int rank_>
2369  inline
2371  unrolled_to_component_indices
2372  (const unsigned int i,
2373  const std::integral_constant<int, rank_> &)
2374  {
2375  (void)i;
2378  Assert (false, ExcNotImplemented());
2379  return TableIndices<rank_>();
2380  }
2381 
2382  }
2383  }
2384 }
2385 
2386 template <int rank_, int dim, typename Number>
2387 inline
2390 (const unsigned int i)
2391 {
2392  return
2393  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2394  std::integral_constant<int, rank_>());
2395 }
2396 
2397 
2398 
2399 template <int rank_, int dim, typename Number>
2400 template <class Archive>
2401 inline
2402 void
2403 SymmetricTensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
2404 {
2405  ar &data;
2406 }
2407 
2408 
2409 #endif // DOXYGEN
2410 
2411 /* ----------------- Non-member functions operating on tensors. ------------ */
2412 
2413 
2426 template <int rank_, int dim, typename Number, typename OtherNumber>
2427 inline
2431 {
2433  tmp += right;
2434  return tmp;
2435 }
2436 
2437 
2450 template <int rank_, int dim, typename Number, typename OtherNumber>
2451 inline
2455 {
2457  tmp -= right;
2458  return tmp;
2459 }
2460 
2461 
2469 template <int rank_, int dim, typename Number, typename OtherNumber>
2470 inline
2473  const Tensor<rank_, dim, OtherNumber> &right)
2474 {
2475  return Tensor<rank_, dim, Number>(left) + right;
2476 }
2477 
2478 
2486 template <int rank_, int dim, typename Number, typename OtherNumber>
2487 inline
2491 {
2492  return left + Tensor<rank_, dim, OtherNumber>(right);
2493 }
2494 
2495 
2503 template <int rank_, int dim, typename Number, typename OtherNumber>
2504 inline
2507  const Tensor<rank_, dim, OtherNumber> &right)
2508 {
2509  return Tensor<rank_, dim, Number>(left) - right;
2510 }
2511 
2512 
2520 template <int rank_, int dim, typename Number, typename OtherNumber>
2521 inline
2525 {
2526  return left - Tensor<rank_, dim, OtherNumber>(right);
2527 }
2528 
2529 
2530 
2544 template <int dim, typename Number>
2545 inline
2547 {
2548  switch (dim)
2549  {
2550  case 1:
2551  return t.data[0];
2552  case 2:
2553  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2554  case 3:
2555  {
2556  // in analogy to general tensors, but
2557  // there's something to be simplified for
2558  // the present case
2559  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2560  return ( tmp + tmp
2561  +t.data[0]*t.data[1]*t.data[2]
2562  -t.data[0]*t.data[5]*t.data[5]
2563  -t.data[1]*t.data[4]*t.data[4]
2564  -t.data[2]*t.data[3]*t.data[3]);
2565  }
2566  default:
2567  Assert (false, ExcNotImplemented());
2569  }
2570 }
2571 
2572 
2573 
2583 template <int dim, typename Number>
2584 inline
2586 {
2587  return determinant (t);
2588 }
2589 
2590 
2591 
2599 template <int dim, typename Number>
2601 {
2602  Number t = d.data[0];
2603  for (unsigned int i=1; i<dim; ++i)
2604  t += d.data[i];
2605  return t;
2606 }
2607 
2608 
2618 template <int dim, typename Number>
2619 inline
2621 {
2622  return trace (t);
2623 }
2624 
2625 
2638 template <typename Number>
2639 inline
2641 {
2643 }
2644 
2645 
2646 
2667 template <typename Number>
2668 inline
2670 {
2671  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2672 }
2673 
2674 
2675 
2685 template <typename Number>
2686 inline
2688 {
2689  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2690  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2691 }
2692 
2693 
2694 
2703 template <typename Number>
2704 std::array<Number,1>
2705 eigenvalues (const SymmetricTensor<2,1,Number> &T);
2706 
2707 
2708 
2730 template <typename Number>
2731 std::array<Number,2>
2732 eigenvalues (const SymmetricTensor<2,2,Number> &T);
2733 
2734 
2735 
2755 template <typename Number>
2756 std::array<Number,3>
2757 eigenvalues (const SymmetricTensor<2,3,Number> &T);
2758 
2759 
2760 
2761 namespace internal
2762 {
2763  namespace SymmetricTensor
2764  {
2801  template <int dim, typename Number>
2802  void
2803  tridiagonalize (const ::SymmetricTensor<2,dim,Number> &A,
2804  ::Tensor<2,dim,Number> &Q,
2805  std::array<Number,dim> &d,
2806  std::array<Number,dim-1> &e);
2807 
2808 
2809 
2848  template <int dim, typename Number>
2849  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2850  ql_implicit_shifts (const ::SymmetricTensor<2,dim,Number> &A);
2851 
2852 
2853 
2892  template <int dim, typename Number>
2893  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
2895 
2896 
2897 
2913  template <typename Number>
2914  std::array<std::pair<Number, Tensor<1,2,Number> >,2>
2915  hybrid (const ::SymmetricTensor<2,2,Number> &A);
2916 
2917 
2918 
2951  template <typename Number>
2952  std::array<std::pair<Number, Tensor<1,3,Number> >,3>
2953  hybrid (const ::SymmetricTensor<2,3,Number> &A);
2954 
2955  namespace
2956  {
2957 
2962  template <int dim, typename Number>
2963  struct SortEigenValuesVectors
2964  {
2965  typedef std::pair<Number, Tensor<1,dim,Number> > EigValsVecs;
2966  bool operator() (const EigValsVecs &lhs,
2967  const EigValsVecs &rhs)
2968  {
2969  return lhs.first > rhs.first;
2970  }
2971  };
2972 
2973  }
2974 
2975  } // namespace SymmetricTensor
2976 
2977 } // namespace internal
2978 
2979 
2980 
2981 // The line below is to ensure that doxygen puts the full description
2982 // of this global enumeration into the documentation
2983 // See https://stackoverflow.com/a/1717984
3011 {
3021  hybrid,
3031  ql_implicit_shifts,
3039  jacobi
3040 };
3041 
3042 
3043 
3050 template <typename Number>
3051 std::array<std::pair<Number, Tensor<1,1,Number> >,1>
3054 {
3055  return { {std::make_pair(T[0][0], Tensor<1,1,Number>({1.0}))} };
3056 }
3057 
3058 
3059 
3087 template <int dim, typename Number>
3088 std::array<std::pair<Number, Tensor<1,dim,Number> >,dim>
3091 {
3092  std::array<std::pair<Number, Tensor<1,dim,Number> >,dim> eig_vals_vecs;
3093 
3094  switch (method)
3095  {
3097  eig_vals_vecs = internal::SymmetricTensor::hybrid(T);
3098  break;
3101  break;
3103  eig_vals_vecs = internal::SymmetricTensor::jacobi(T);
3104  break;
3105  default:
3106  AssertThrow(false, ExcNotImplemented());
3107  }
3108 
3109  // Sort in descending order before output.
3110  std::sort(eig_vals_vecs.begin(), eig_vals_vecs.end(),
3111  internal::SymmetricTensor::SortEigenValuesVectors<dim,Number>());
3112  return eig_vals_vecs;
3113 }
3114 
3115 
3116 
3126 template <int rank_, int dim, typename Number>
3127 inline
3130 {
3131  return t;
3132 }
3133 
3134 
3135 
3145 template <int dim, typename Number>
3146 inline
3149 {
3151 
3152  // subtract scaled trace from the diagonal
3153  const Number tr = trace(t) / dim;
3154  for (unsigned int i=0; i<dim; ++i)
3155  tmp.data[i] -= tr;
3156 
3157  return tmp;
3158 }
3159 
3160 
3161 
3169 template <int dim, typename Number>
3170 inline
3172 unit_symmetric_tensor ()
3173 {
3174  // create a default constructed matrix filled with
3175  // zeros, then set the diagonal elements to one
3177  switch (dim)
3178  {
3179  case 1:
3180  tmp.data[0] = 1;
3181  break;
3182  case 2:
3183  tmp.data[0] = tmp.data[1] = 1;
3184  break;
3185  case 3:
3186  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
3187  break;
3188  default:
3189  for (unsigned int d=0; d<dim; ++d)
3190  tmp.data[d] = 1;
3191  }
3192  return tmp;
3193 }
3194 
3195 
3196 
3205 template <int dim>
3206 inline
3208 unit_symmetric_tensor ()
3209 {
3210  return unit_symmetric_tensor<dim,double>();
3211 }
3212 
3213 
3214 
3229 template <int dim, typename Number>
3230 inline
3232 deviator_tensor ()
3233 {
3235 
3236  // fill the elements treating the diagonal
3237  for (unsigned int i=0; i<dim; ++i)
3238  for (unsigned int j=0; j<dim; ++j)
3239  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
3240 
3241  // then fill the ones that copy over the
3242  // non-diagonal elements. note that during
3243  // the double-contraction, we handle the
3244  // off-diagonal elements twice, so simply
3245  // copying requires a weight of 1/2
3246  for (unsigned int i=dim;
3247  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3248  ++i)
3249  tmp.data[i][i] = 0.5;
3250 
3251  return tmp;
3252 }
3253 
3254 
3255 
3270 template <int dim>
3271 inline
3273 deviator_tensor ()
3274 {
3275  return deviator_tensor<dim,double>();
3276 }
3277 
3278 
3279 
3302 template <int dim, typename Number>
3303 inline
3305 identity_tensor ()
3306 {
3308 
3309  // fill the elements treating the diagonal
3310  for (unsigned int i=0; i<dim; ++i)
3311  tmp.data[i][i] = 1;
3312 
3313  // then fill the ones that copy over the
3314  // non-diagonal elements. note that during
3315  // the double-contraction, we handle the
3316  // off-diagonal elements twice, so simply
3317  // copying requires a weight of 1/2
3318  for (unsigned int i=dim;
3319  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
3320  ++i)
3321  tmp.data[i][i] = 0.5;
3322 
3323  return tmp;
3324 }
3325 
3326 
3327 
3349 template <int dim>
3350 inline
3352 identity_tensor ()
3353 {
3354  return identity_tensor<dim,double>();
3355 }
3356 
3357 
3358 
3369 template <int dim, typename Number>
3370 inline
3373 {
3375 }
3376 
3377 
3378 
3390 template <int dim, typename Number>
3391 inline
3394 {
3396 }
3397 
3398 
3399 
3414 template <int dim, typename Number>
3415 inline
3419 {
3421 
3422  // fill only the elements really needed
3423  for (unsigned int i=0; i<dim; ++i)
3424  for (unsigned int j=i; j<dim; ++j)
3425  for (unsigned int k=0; k<dim; ++k)
3426  for (unsigned int l=k; l<dim; ++l)
3427  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
3428 
3429  return tmp;
3430 }
3431 
3432 
3433 
3442 template <int dim,typename Number>
3443 inline
3446 {
3447  Number array[(dim*dim+dim)/2];
3448  for (unsigned int d=0; d<dim; ++d)
3449  array[d] = t[d][d];
3450  for (unsigned int d=0, c=0; d<dim; ++d)
3451  for (unsigned int e=d+1; e<dim; ++e, ++c)
3452  array[dim+c] = (t[d][e]+t[e][d])*0.5;
3453  return SymmetricTensor<2,dim,Number>(array);
3454 }
3455 
3456 
3457 
3465 template <int rank_, int dim, typename Number>
3466 inline
3469  const Number &factor)
3470 {
3472  tt *= factor;
3473  return tt;
3474 }
3475 
3476 
3477 
3485 template <int rank_, int dim, typename Number>
3486 inline
3488 operator * (const Number &factor,
3490 {
3491  // simply forward to the other operator
3492  return t*factor;
3493 }
3494 
3495 
3496 
3522 template <int rank_, int dim, typename Number, typename OtherNumber>
3523 inline
3526  const OtherNumber &factor)
3527 {
3528  // form the product. we have to convert the two factors into the final
3529  // type via explicit casts because, for awkward reasons, the C++
3530  // standard committee saw it fit to not define an
3531  // operator*(float,std::complex<double>)
3532  // (as well as with switched arguments and double<->float).
3533  typedef typename ProductType<Number,OtherNumber>::type product_type;
3535  // we used to shorten the following by 'tt *= product_type(factor);'
3536  // which requires that a converting constructor
3537  // 'product_type::product_type(const OtherNumber) is defined.
3538  // however, a user-defined constructor is not allowed for aggregates,
3539  // e.g. VectorizedArray. therefore, we work around this issue using a
3540  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
3541  // which we assume to be defined.
3542  product_type new_factor;
3543  new_factor = factor;
3544  tt *= new_factor;
3545  return tt;
3546 }
3547 
3548 
3549 
3558 template <int rank_, int dim, typename Number, typename OtherNumber>
3559 inline
3561 operator * (const Number &factor,
3563 {
3564  // simply forward to the other operator with switched arguments
3565  return (t*factor);
3566 }
3567 
3568 
3569 
3575 template <int rank_, int dim, typename Number, typename OtherNumber>
3576 inline
3579  const OtherNumber &factor)
3580 {
3582  tt /= factor;
3583  return tt;
3584 }
3585 
3586 
3587 
3594 template <int rank_, int dim>
3595 inline
3597 operator * (const SymmetricTensor<rank_,dim> &t,
3598  const double factor)
3599 {
3601  tt *= factor;
3602  return tt;
3603 }
3604 
3605 
3606 
3613 template <int rank_, int dim>
3614 inline
3616 operator * (const double factor,
3617  const SymmetricTensor<rank_,dim> &t)
3618 {
3620  tt *= factor;
3621  return tt;
3622 }
3623 
3624 
3625 
3631 template <int rank_, int dim>
3632 inline
3634 operator / (const SymmetricTensor<rank_,dim> &t,
3635  const double factor)
3636 {
3638  tt /= factor;
3639  return tt;
3640 }
3641 
3651 template <int dim, typename Number, typename OtherNumber>
3652 inline
3653 typename ProductType<Number, OtherNumber>::type
3656 {
3657  return (t1*t2);
3658 }
3659 
3660 
3670 template <int dim, typename Number, typename OtherNumber>
3671 inline
3672 typename ProductType<Number, OtherNumber>::type
3674  const Tensor<2,dim,OtherNumber> &t2)
3675 {
3676  typename ProductType<Number, OtherNumber>::type s = internal::NumberType<typename ProductType<Number, OtherNumber>::type>::value(0.0);
3677  for (unsigned int i=0; i<dim; ++i)
3678  for (unsigned int j=0; j<dim; ++j)
3679  s += t1[i][j] * t2[i][j];
3680  return s;
3681 }
3682 
3683 
3693 template <int dim, typename Number, typename OtherNumber>
3694 inline
3695 typename ProductType<Number, OtherNumber>::type
3698 {
3699  return scalar_product(t2, t1);
3700 }
3701 
3702 
3718 template <typename Number, typename OtherNumber>
3719 inline
3720 void
3721 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3722  const SymmetricTensor<4,1,Number> &t,
3724 {
3725  tmp[0][0] = t[0][0][0][0] * s[0][0];
3726 }
3727 
3728 
3729 
3745 template <typename Number, typename OtherNumber>
3746 inline
3747 void
3748 double_contract (SymmetricTensor<2,1,typename ProductType<Number, OtherNumber>::type> &tmp,
3749  const SymmetricTensor<2,1,Number> &s,
3751 {
3752  tmp[0][0] = t[0][0][0][0] * s[0][0];
3753 }
3754 
3755 
3756 
3771 template <typename Number, typename OtherNumber>
3772 inline
3773 void
3774 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3775  const SymmetricTensor<4,2,Number> &t,
3777 {
3778  const unsigned int dim = 2;
3779 
3780  for (unsigned int i=0; i<dim; ++i)
3781  for (unsigned int j=i; j<dim; ++j)
3782  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3783  t[i][j][1][1] * s[1][1] +
3784  2 * t[i][j][0][1] * s[0][1];
3785 }
3786 
3787 
3788 
3804 template <typename Number, typename OtherNumber>
3805 inline
3806 void
3807 double_contract (SymmetricTensor<2,2,typename ProductType<Number, OtherNumber>::type> &tmp,
3808  const SymmetricTensor<2,2,Number> &s,
3810 {
3811  const unsigned int dim = 2;
3812 
3813  for (unsigned int i=0; i<dim; ++i)
3814  for (unsigned int j=i; j<dim; ++j)
3815  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3816  s[1][1] * t[1][1][i][j] +
3817  2 * s[0][1] * t[0][1][i][j];
3818 }
3819 
3820 
3821 
3837 template <typename Number, typename OtherNumber>
3838 inline
3839 void
3840 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3841  const SymmetricTensor<4,3,Number> &t,
3843 {
3844  const unsigned int dim = 3;
3845 
3846  for (unsigned int i=0; i<dim; ++i)
3847  for (unsigned int j=i; j<dim; ++j)
3848  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3849  t[i][j][1][1] * s[1][1] +
3850  t[i][j][2][2] * s[2][2] +
3851  2 * t[i][j][0][1] * s[0][1] +
3852  2 * t[i][j][0][2] * s[0][2] +
3853  2 * t[i][j][1][2] * s[1][2];
3854 }
3855 
3856 
3857 
3873 template <typename Number, typename OtherNumber>
3874 inline
3875 void
3876 double_contract (SymmetricTensor<2,3,typename ProductType<Number, OtherNumber>::type> &tmp,
3877  const SymmetricTensor<2,3,Number> &s,
3879 {
3880  const unsigned int dim = 3;
3881 
3882  for (unsigned int i=0; i<dim; ++i)
3883  for (unsigned int j=i; j<dim; ++j)
3884  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3885  s[1][1] * t[1][1][i][j] +
3886  s[2][2] * t[2][2][i][j] +
3887  2 * s[0][1] * t[0][1][i][j] +
3888  2 * s[0][2] * t[0][2][i][j] +
3889  2 * s[1][2] * t[1][2][i][j];
3890 }
3891 
3892 
3893 
3901 template <int dim, typename Number, typename OtherNumber>
3903 operator * (const SymmetricTensor<2,dim,Number> &src1,
3904  const Tensor<1,dim,OtherNumber> &src2)
3905 {
3907  for (unsigned int i=0; i<dim; ++i)
3908  for (unsigned int j=0; j<dim; ++j)
3909  dest[i] += src1[i][j] * src2[j];
3910  return dest;
3911 }
3912 
3913 
3921 template <int dim, typename Number, typename OtherNumber>
3923 operator * (const Tensor<1,dim,Number> &src1,
3925 {
3926  // this is easy for symmetric tensors:
3927  return src2 * src1;
3928 }
3929 
3930 
3931 
3952 template <int rank_1, int rank_2, int dim,
3953  typename Number, typename OtherNumber>
3954 inline DEAL_II_ALWAYS_INLINE
3955 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3956 operator * (const Tensor<rank_1, dim, Number> &src1,
3958 {
3959  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3960  const Tensor<rank_2, dim, OtherNumber> src2 (src2s);
3961  return src1*src2;
3962 }
3963 
3964 
3965 
3986 template <int rank_1, int rank_2, int dim,
3987  typename Number, typename OtherNumber>
3988 inline DEAL_II_ALWAYS_INLINE
3989 typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type
3990 operator * (const SymmetricTensor<rank_1, dim, Number> &src1s,
3992 {
3993  typename Tensor<rank_1 + rank_2 - 2, dim, typename ProductType<Number, OtherNumber>::type>::tensor_type result;
3994  const Tensor<rank_2, dim, OtherNumber> src1 (src1s);
3995  return src1*src2;
3996 }
3997 
3998 
3999 
4009 template <int dim, typename Number>
4010 inline
4011 std::ostream &operator << (std::ostream &out,
4013 {
4014  //make out lives a bit simpler by outputing
4015  //the tensor through the operator for the
4016  //general Tensor class
4018 
4019  for (unsigned int i=0; i<dim; ++i)
4020  for (unsigned int j=0; j<dim; ++j)
4021  tt[i][j] = t[i][j];
4022 
4023  return out << tt;
4024 }
4025 
4026 
4027 
4037 template <int dim, typename Number>
4038 inline
4039 std::ostream &operator << (std::ostream &out,
4041 {
4042  //make out lives a bit simpler by outputing
4043  //the tensor through the operator for the
4044  //general Tensor class
4046 
4047  for (unsigned int i=0; i<dim; ++i)
4048  for (unsigned int j=0; j<dim; ++j)
4049  for (unsigned int k=0; k<dim; ++k)
4050  for (unsigned int l=0; l<dim; ++l)
4051  tt[i][j][k][l] = t[i][j][k][l];
4052 
4053  return out << tt;
4054 }
4055 
4056 
4057 DEAL_II_NAMESPACE_CLOSE
4058 
4059 #endif
friend SymmetricTensor< 4, dim2, Number2 > identity_tensor()
void tridiagonalize(const ::SymmetricTensor< 2, dim, Number > &A,::Tensor< 2, dim, Number > &Q, std::array< Number, dim > &d, std::array< Number, dim-1 > &e)
static const unsigned int invalid_unsigned_int
Definition: types.h:173
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static const unsigned int n_independent_components
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
std::array< std::pair< Number, Tensor< 1, 2, Number > >, 2 > hybrid(const ::SymmetricTensor< 2, 2, Number > &A)
void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 2, Number > &s, const SymmetricTensor< 4, 2, OtherNumber > &t)
void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, OtherNumber > &s)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1207
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, OtherNumber > &t2)
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
bool operator==(const SymmetricTensor &) const
SymmetricTensor & operator=(const SymmetricTensor< rank_, dim, OtherNumber > &rhs)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
static std::size_t memory_consumption()
numbers::NumberTraits< Number >::real_type norm() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:398
static real_type abs(const number &x)
Definition: numbers.h:352
SymmetricTensorEigenvectorMethod
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
internal::SymmetricTensorAccessors::StorageType< rank_, dim, Number > base_tensor_descriptor
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor & operator/=(const OtherNumber &factor)
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
void double_contract(SymmetricTensor< 2, 2, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, OtherNumber > &s)
static const unsigned int dimension
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static::ExceptionBase & ExcMessage(std::string arg1)
TableIndices< 4 > merge(const TableIndices< 4 > &previous_indices, const unsigned int new_index, const unsigned int position)
internal::SymmetricTensorAccessors::Accessor< rank_, dim, true, rank_-1, Number > operator[](const unsigned int row) const
ProductType< Number, OtherNumber >::type scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
base_tensor_descriptor::base_tensor_type base_tensor_type
friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
static const unsigned int rank
T sum(const T &t, const MPI_Comm &mpi_communicator)
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
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_type data
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
Number trace(const SymmetricTensor< 2, dim, Number > &d)
void serialize(Archive &ar, const unsigned int version)
friend SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const Tensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Number & operator()(const TableIndices< rank_ > &indices)
Number access_raw_entry(const unsigned int unrolled_index) const
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
internal::SymmetricTensorAccessors::double_contraction_result< rank_, 2, dim, Number, OtherNumber >::type operator*(const SymmetricTensor< 2, dim, OtherNumber > &s) const
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
Definition: mpi.h:53
SymmetricTensor & operator*=(const OtherNumber &factor)
SymmetricTensor & operator+=(const SymmetricTensor< rank_, dim, OtherNumber > &)
StreamType & operator<<(StreamType &s, const UpdateFlags u)
bool operator!=(const SymmetricTensor &) const
Number third_invariant(const SymmetricTensor< 2, dim, Number > &t)
Tensor< 1, n_independent_components, Number > base_tensor_type
static::ExceptionBase & ExcNotImplemented()
SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > ql_implicit_shifts(const ::SymmetricTensor< 2, dim, Number > &A)
Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
friend SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
Tensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const Tensor< rank_, dim, OtherNumber > &right)
void double_contract(SymmetricTensor< 2, 3, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, OtherNumber > &t)
T max(const T &t, const MPI_Comm &mpi_communicator)
void double_contract(SymmetricTensor< 2, 1, typename ProductType< Number, OtherNumber >::type > &tmp, const SymmetricTensor< 2, 1, Number > &s, const SymmetricTensor< 4, 1, OtherNumber > &t)
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
std::array< std::pair< Number, Tensor< 1, 1, Number > >, 1 > eigenvectors(const SymmetricTensor< 2, 1, Number > &T, const SymmetricTensorEigenvectorMethod=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
SymmetricTensor & operator-=(const SymmetricTensor< rank_, dim, OtherNumber > &)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor operator-() const
static::ExceptionBase & ExcInternalError()