Reference documentation for deal.II version Git c025e4d 2017-04-26 14:50:03 -0500
symmetric_tensor.h
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 DEAL_II_NAMESPACE_OPEN
26 
27 template <int rank, int dim, typename Number=double> class SymmetricTensor;
28 
29 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
31 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
33 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
35 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
37 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
39 template <int dim2, typename Number> Number
41 
42 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
44 template <int dim, typename Number> Number
46 
47 
48 
49 namespace internal
50 {
55  namespace SymmetricTensorAccessors
56  {
63  inline
64  TableIndices<2> merge (const TableIndices<2> &previous_indices,
65  const unsigned int new_index,
66  const unsigned int position)
67  {
68  Assert (position < 2, ExcIndexRange (position, 0, 2));
69 
70  if (position == 0)
71  return TableIndices<2>(new_index);
72  else
73  return TableIndices<2>(previous_indices[0], new_index);
74  }
75 
76 
77 
84  inline
85  TableIndices<4> merge (const TableIndices<4> &previous_indices,
86  const unsigned int new_index,
87  const unsigned int position)
88  {
89  Assert (position < 4, ExcIndexRange (position, 0, 4));
90 
91  switch (position)
92  {
93  case 0:
94  return TableIndices<4>(new_index);
95  case 1:
96  return TableIndices<4>(previous_indices[0],
97  new_index);
98  case 2:
99  return TableIndices<4>(previous_indices[0],
100  previous_indices[1],
101  new_index);
102  case 3:
103  return TableIndices<4>(previous_indices[0],
104  previous_indices[1],
105  previous_indices[2],
106  new_index);
107  }
108  Assert (false, ExcInternalError());
109  return TableIndices<4>();
110  }
111 
112 
121  template <int rank1, int rank2, int dim, typename Number>
123  {
124  typedef ::SymmetricTensor<rank1+rank2-4,dim,Number> type;
125  };
126 
127 
136  template <int dim, typename Number>
137  struct double_contraction_result<2,2,dim,Number>
138  {
139  typedef Number type;
140  };
141 
142 
143 
156  template <int rank, int dim, typename Number>
157  struct StorageType;
158 
162  template <int dim, typename Number>
163  struct StorageType<2,dim,Number>
164  {
169  static const unsigned int
170  n_independent_components = (dim*dim + dim)/2;
171 
176  };
177 
178 
179 
183  template <int dim, typename Number>
184  struct StorageType<4,dim,Number>
185  {
191  static const unsigned int
192  n_rank2_components = (dim*dim + dim)/2;
193 
197  static const unsigned int
198  n_independent_components = (n_rank2_components *
200 
208  };
209 
210 
211 
216  template <int rank, int dim, bool constness, typename Number>
218 
225  template <int rank, int dim, typename Number>
226  struct AccessorTypes<rank,dim,true,Number>
227  {
228  typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
229 
230  typedef Number reference;
231  };
232 
239  template <int rank, int dim, typename Number>
240  struct AccessorTypes<rank,dim,false,Number>
241  {
242  typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
243 
244  typedef Number &reference;
245  };
246 
247 
282  template <int rank, int dim, bool constness, int P, typename Number>
283  class Accessor
284  {
285  public:
289  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
290  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
291 
292  private:
311  Accessor (tensor_type &tensor,
312  const TableIndices<rank> &previous_indices);
313 
317  Accessor (const Accessor &a);
318 
319  public:
320 
324  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
325 
329  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;
330 
331  private:
335  tensor_type &tensor;
336  const TableIndices<rank> previous_indices;
337 
338  // declare some other classes
339  // as friends. make sure to
340  // work around bugs in some
341  // compilers
342  template <int,int,typename> friend class ::SymmetricTensor;
343  template <int,int,bool,int,typename>
344  friend class Accessor;
345 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
346  friend class ::SymmetricTensor<rank,dim,Number>;
347  friend class Accessor<rank,dim,constness,P+1,Number>;
348 # endif
349  };
350 
351 
352 
362  template <int rank, int dim, bool constness, typename Number>
363  class Accessor<rank,dim,constness,1,Number>
364  {
365  public:
369  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
370  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
371 
372  private:
394  Accessor (tensor_type &tensor,
395  const TableIndices<rank> &previous_indices);
396 
400  Accessor ();
401 
405  Accessor (const Accessor &a);
406 
407  public:
408 
412  reference operator [] (const unsigned int);
413 
417  reference operator [] (const unsigned int) const;
418 
419  private:
423  tensor_type &tensor;
424  const TableIndices<rank> previous_indices;
425 
426  // declare some other classes
427  // as friends. make sure to
428  // work around bugs in some
429  // compilers
430  template <int,int,typename> friend class ::SymmetricTensor;
431  template <int,int,bool,int,typename>
432  friend class SymmetricTensorAccessors::Accessor;
433 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
434  friend class ::SymmetricTensor<rank,dim,Number>;
435  friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
436 # endif
437  };
438  }
439 }
440 
441 
442 
506 template <int rank, int dim, typename Number>
507 class SymmetricTensor
508 {
509 public:
518  static const unsigned int dimension = dim;
519 
525  static const unsigned int n_independent_components
528 
532  SymmetricTensor ();
533 
545 
561  SymmetricTensor (const Number (&array) [n_independent_components]);
562 
568  template <typename OtherNumber>
569  explicit
571 
578  SymmetricTensor &operator = (const Number d);
579 
584  operator Tensor<rank,dim,Number> () const;
585 
589  bool operator == (const SymmetricTensor &) const;
590 
594  bool operator != (const SymmetricTensor &) const;
595 
600 
605 
610  SymmetricTensor &operator *= (const Number factor);
611 
615  SymmetricTensor &operator /= (const Number factor);
616 
622 
628 
633 
660 
667 
671  Number &operator() (const TableIndices<rank> &indices);
672 
676  Number operator() (const TableIndices<rank> &indices) const;
677 
682  internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
683  operator [] (const unsigned int row) const;
684 
689  internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
690  operator [] (const unsigned int row);
691 
697  Number
698  operator [] (const TableIndices<rank> &indices) const;
699 
705  Number &
706  operator [] (const TableIndices<rank> &indices);
707 
713  Number
714  access_raw_entry (const unsigned int unrolled_index) const;
715 
721  Number &
722  access_raw_entry (const unsigned int unrolled_index);
723 
734  norm () const;
735 
743  static
744  unsigned int
746 
752  static
754  unrolled_to_component_indices (const unsigned int i);
755 
768  void clear ();
769 
774  static std::size_t memory_consumption ();
775 
780  template <class Archive>
781  void serialize(Archive &ar, const unsigned int version);
782 
783 private:
787  typedef
790 
794  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
795 
799  base_tensor_type data;
800 
804  template <int, int, typename> friend class SymmetricTensor;
805 
809  template <int dim2, typename Number2>
810  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
811 
812  template <int dim2, typename Number2>
813  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
814 
815  template <int dim2, typename Number2>
817  deviator (const SymmetricTensor<2,dim2,Number2> &t);
818 
819  template <int dim2, typename Number2>
821 
822  template <int dim2, typename Number2>
824 
825  template <int dim2, typename Number2>
827 
828  template <int dim2, typename Number2>
830 
831  template <int dim2, typename Number2>
833 };
834 
835 
836 
837 // ------------------------- inline functions ------------------------
838 
839 #ifndef DOXYGEN
840 
841 namespace internal
842 {
843  namespace SymmetricTensorAccessors
844  {
845  template <int rank, int dim, bool constness, int P, typename Number>
846  Accessor<rank,dim,constness,P,Number>::
847  Accessor (tensor_type &tensor,
848  const TableIndices<rank> &previous_indices)
849  :
850  tensor (tensor),
851  previous_indices (previous_indices)
852  {}
853 
854 
855  template <int rank, int dim, bool constness, int P, typename Number>
856  Accessor<rank,dim,constness,P,Number>::
857  Accessor (const Accessor &a)
858  :
859  tensor (a.tensor),
860  previous_indices (a.previous_indices)
861  {}
862 
863 
864 
865  template <int rank, int dim, bool constness, int P, typename Number>
866  Accessor<rank,dim,constness,P-1,Number>
867  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
868  {
869  return Accessor<rank,dim,constness,P-1,Number> (tensor,
870  merge (previous_indices, i, rank-P));
871  }
872 
873 
874 
875  template <int rank, int dim, bool constness, int P, typename Number>
876  Accessor<rank,dim,constness,P-1,Number>
877  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i) const
878  {
879  return Accessor<rank,dim,constness,P-1,Number> (tensor,
880  merge (previous_indices, i, rank-P));
881  }
882 
883 
884 
885  template <int rank, int dim, bool constness, typename Number>
886  Accessor<rank,dim,constness,1,Number>::
887  Accessor (tensor_type &tensor,
888  const TableIndices<rank> &previous_indices)
889  :
890  tensor (tensor),
891  previous_indices (previous_indices)
892  {}
893 
894 
895 
896  template <int rank, int dim, bool constness, typename Number>
897  Accessor<rank,dim,constness,1,Number>::
898  Accessor (const Accessor &a)
899  :
900  tensor (a.tensor),
901  previous_indices (a.previous_indices)
902  {}
903 
904 
905 
906  template <int rank, int dim, bool constness, typename Number>
907  typename Accessor<rank,dim,constness,1,Number>::reference
908  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
909  {
910  return tensor(merge (previous_indices, i, rank-1));
911  }
912 
913 
914  template <int rank, int dim, bool constness, typename Number>
915  typename Accessor<rank,dim,constness,1,Number>::reference
916  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i) const
917  {
918  return tensor(merge (previous_indices, i, rank-1));
919  }
920  }
921 }
922 
923 
924 
925 template <int rank, int dim, typename Number>
926 inline
928 {}
929 
930 
931 
932 template <int rank, int dim, typename Number>
933 inline
935 {
936  Assert (rank == 2, ExcNotImplemented());
937  switch (dim)
938  {
939  case 2:
940  Assert (t[0][1] == t[1][0], ExcInternalError());
941 
942  data[0] = t[0][0];
943  data[1] = t[1][1];
944  data[2] = t[0][1];
945 
946  break;
947  case 3:
948  Assert (t[0][1] == t[1][0], ExcInternalError());
949  Assert (t[0][2] == t[2][0], ExcInternalError());
950  Assert (t[1][2] == t[2][1], ExcInternalError());
951 
952  data[0] = t[0][0];
953  data[1] = t[1][1];
954  data[2] = t[2][2];
955  data[3] = t[0][1];
956  data[4] = t[0][2];
957  data[5] = t[1][2];
958 
959  break;
960  default:
961  for (unsigned int d=0; d<dim; ++d)
962  for (unsigned int e=0; e<d; ++e)
963  Assert(t[d][e] == t[e][d], ExcInternalError());
964 
965  for (unsigned int d=0; d<dim; ++d)
966  data[d] = t[d][d];
967 
968  for (unsigned int d=0, c=0; d<dim; ++d)
969  for (unsigned int e=d+1; e<dim; ++e, ++c)
970  data[dim+c] = t[d][e];
971  }
972 }
973 
974 
975 
976 template <int rank, int dim, typename Number>
977 template <typename OtherNumber>
978 inline
981 {
982  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
983  data[i] = initializer.data[i];
984 }
985 
986 
987 
988 
989 template <int rank, int dim, typename Number>
990 inline
991 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
992  :
993  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
994 {
995  // ensure that the reinterpret_cast above actually works
996  Assert (sizeof(typename base_tensor_type::array_type)
997  == sizeof(array),
998  ExcInternalError());
999 }
1000 
1001 
1002 
1003 template <int rank, int dim, typename Number>
1004 inline
1007 {
1008  Assert (d==Number(), ExcMessage ("Only assignment with zero is allowed"));
1009  (void) d;
1010 
1011  data = 0;
1012 
1013  return *this;
1014 }
1015 
1016 
1017 namespace internal
1018 {
1019  namespace SymmetricTensor
1020  {
1021  template <int dim, typename Number>
1023  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1024  {
1025  Number t[dim][dim];
1026 
1027  // diagonal entries are stored first
1028  for (unsigned int d=0; d<dim; ++d)
1029  t[d][d] = s.access_raw_entry(d);
1030 
1031  // off-diagonal entries come next, row by row
1032  for (unsigned int d=0, c=0; d<dim; ++d)
1033  for (unsigned int e=d+1; e<dim; ++e, ++c)
1034  {
1035  t[d][e] = s.access_raw_entry(dim+c);
1036  t[e][d] = s.access_raw_entry(dim+c);
1037  }
1038  return ::Tensor<2,dim,Number>(t);
1039  }
1040 
1041 
1042  template <int dim, typename Number>
1044  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1045  {
1046  // utilize the symmetry properties of SymmetricTensor<4,dim>
1047  // discussed in the class documentation to avoid accessing all
1048  // independent elements of the input tensor more than once
1050 
1051  for (unsigned int i=0; i<dim; ++i)
1052  for (unsigned int j=i; j<dim; ++j)
1053  for (unsigned int k=0; k<dim; ++k)
1054  for (unsigned int l=k; l<dim; ++l)
1055  t[TableIndices<4>(i,j,k,l)]
1056  = t[TableIndices<4>(i,j,l,k)]
1057  = t[TableIndices<4>(j,i,k,l)]
1058  = t[TableIndices<4>(j,i,l,k)]
1059  = st[TableIndices<4>(i,j,k,l)];
1060 
1061  return t;
1062  }
1063  }
1064 }
1065 
1066 
1067 
1068 template <int rank, int dim, typename Number>
1069 inline
1071 operator Tensor<rank,dim,Number> () const
1072 {
1073  return internal::SymmetricTensor::convert_to_tensor (*this);
1074 }
1075 
1076 
1077 
1078 template <int rank, int dim, typename Number>
1079 inline
1080 bool
1082 (const SymmetricTensor<rank,dim,Number> &t) const
1083 {
1084  return data == t.data;
1085 }
1086 
1087 
1088 
1089 template <int rank, int dim, typename Number>
1090 inline
1091 bool
1092 SymmetricTensor<rank,dim,Number>::operator !=
1093 (const SymmetricTensor<rank,dim,Number> &t) const
1094 {
1095  return data != t.data;
1096 }
1097 
1098 
1099 
1100 template <int rank, int dim, typename Number>
1101 inline
1103 SymmetricTensor<rank,dim,Number>::operator +=
1105 {
1106  data += t.data;
1107  return *this;
1108 }
1109 
1110 
1111 
1112 template <int rank, int dim, typename Number>
1113 inline
1115 SymmetricTensor<rank,dim,Number>::operator -=
1117 {
1118  data -= t.data;
1119  return *this;
1120 }
1121 
1122 
1123 
1124 template <int rank, int dim, typename Number>
1125 inline
1128 {
1129  data *= d;
1130  return *this;
1131 }
1132 
1133 
1134 
1135 template <int rank, int dim, typename Number>
1136 inline
1139 {
1140  data /= d;
1141  return *this;
1142 }
1143 
1144 
1145 
1146 template <int rank, int dim, typename Number>
1147 inline
1150 {
1151  SymmetricTensor tmp = *this;
1152  tmp.data += t.data;
1153  return tmp;
1154 }
1155 
1156 
1157 
1158 template <int rank, int dim, typename Number>
1159 inline
1162 {
1163  SymmetricTensor tmp = *this;
1164  tmp.data -= t.data;
1165  return tmp;
1166 }
1167 
1168 
1169 
1170 template <int rank, int dim, typename Number>
1171 inline
1174 {
1175  SymmetricTensor tmp = *this;
1176  tmp.data = -tmp.data;
1177  return tmp;
1178 }
1179 
1180 
1181 
1182 template <int rank, int dim, typename Number>
1183 inline
1184 void
1186 {
1187  data.clear ();
1188 }
1189 
1190 
1191 
1192 template <int rank, int dim, typename Number>
1193 inline
1194 std::size_t
1196 {
1197  // all memory consists of statically allocated memory of the current
1198  // object, no pointers
1199  return sizeof(SymmetricTensor<rank,dim,Number>);
1200 }
1201 
1202 
1203 
1204 namespace internal
1205 {
1206 
1207  template <int dim, typename Number>
1208  inline
1209  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
1210  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1212  {
1213  switch (dim)
1214  {
1215  case 1:
1216  return data[0] * sdata[0];
1217  default:
1218  // Start with the non-diagonal part to avoid some multiplications by
1219  // 2.
1220  Number sum = data[dim] * sdata[dim];
1221  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1222  sum += data[d] * sdata[d];
1223  sum += sum; // sum = sum * 2.;
1224 
1225  // Now add the contributions from the diagonal
1226  for (unsigned int d=0; d<dim; ++d)
1227  sum += data[d] * sdata[d];
1228  return sum;
1229  }
1230  }
1231 
1232 
1233 
1234  template <int dim, typename Number>
1235  inline
1236  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
1237  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1239  {
1240  const unsigned int data_dim =
1242  Number tmp [data_dim];
1243  for (unsigned int i=0; i<data_dim; ++i)
1244  tmp[i] = perform_double_contraction<dim,Number>(data[i], sdata);
1245  return ::SymmetricTensor<2,dim,Number>(tmp);
1246  }
1247 
1248 
1249 
1250  template <int dim, typename Number>
1251  inline
1253  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1255  {
1257  for (unsigned int i=0; i<tmp.dimension; ++i)
1258  {
1259  // Start with the non-diagonal part
1260  Number sum = data[dim] * sdata[dim][i];
1261  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1262  sum += data[d] * sdata[d][i];
1263  sum += sum; // sum = sum * 2.;
1264 
1265  // Now add the contributions from the diagonal
1266  for (unsigned int d=0; d<dim; ++d)
1267  sum += data[d] * sdata[d][i];
1268  tmp[i] = sum;
1269  }
1270  return tmp;
1271  }
1272 
1273 
1274 
1275  template <int dim, typename Number>
1276  inline
1278  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1280  {
1281  const unsigned int data_dim =
1284  for (unsigned int i=0; i<data_dim; ++i)
1285  for (unsigned int j=0; j<data_dim; ++j)
1286  {
1287  // Start with the non-diagonal part
1288  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1289  tmp[i][j] += data[i][d] * sdata[d][j];
1290  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1291 
1292  // Now add the contributions from the diagonal
1293  for (unsigned int d=0; d<dim; ++d)
1294  tmp[i][j] += data[i][d] * sdata[d][j];
1295  }
1296  return tmp;
1297  }
1298 
1299 } // end of namespace internal
1300 
1301 
1302 
1303 template <int rank, int dim, typename Number>
1304 inline
1307 {
1308  // need to have two different function calls
1309  // because a scalar and rank-2 tensor are not
1310  // the same data type (see internal function
1311  // above)
1312  return internal::perform_double_contraction<dim,Number> (data, s.data);
1313 }
1314 
1315 
1316 
1317 template <int rank, int dim, typename Number>
1318 inline
1321 {
1324  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
1325  return tmp;
1326 }
1327 
1328 
1329 
1330 // internal namespace to switch between the
1331 // access of different tensors. There used to
1332 // be explicit instantiations before for
1333 // different ranks and dimensions, but since
1334 // we now allow for templates on the data
1335 // type, and since we cannot partially
1336 // specialize the implementation, this got
1337 // into a separate namespace
1338 namespace internal
1339 {
1340  template <int dim, typename Number>
1341  inline
1342  Number &
1343  symmetric_tensor_access (const TableIndices<2> &indices,
1345  {
1346  // 1d is very simple and done first
1347  if (dim == 1)
1348  return data[0];
1349 
1350  // first treat the main diagonal elements, which are stored consecutively
1351  // at the beginning
1352  if (indices[0] == indices[1])
1353  return data[indices[0]];
1354 
1355  // the rest is messier and requires a few switches.
1356  switch (dim)
1357  {
1358  case 2:
1359  // at least for the 2x2 case it is reasonably simple
1360  Assert (((indices[0]==1) && (indices[1]==0)) ||
1361  ((indices[0]==0) && (indices[1]==1)),
1362  ExcInternalError());
1363  return data[2];
1364 
1365  default:
1366  // to do the rest, sort our indices before comparing
1367  {
1368  TableIndices<2> sorted_indices (indices);
1369  sorted_indices.sort ();
1370 
1371  for (unsigned int d=0, c=0; d<dim; ++d)
1372  for (unsigned int e=d+1; e<dim; ++e, ++c)
1373  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1374  return data[dim+c];
1375  Assert (false, ExcInternalError());
1376  }
1377  }
1378 
1379  static Number dummy_but_referenceable = Number();
1380  return dummy_but_referenceable;
1381  }
1382 
1383 
1384 
1385  template <int dim, typename Number>
1386  inline
1387  Number
1388  symmetric_tensor_access (const TableIndices<2> &indices,
1390  {
1391  // 1d is very simple and done first
1392  if (dim == 1)
1393  return data[0];
1394 
1395  // first treat the main diagonal elements, which are stored consecutively
1396  // at the beginning
1397  if (indices[0] == indices[1])
1398  return data[indices[0]];
1399 
1400  // the rest is messier and requires a few switches.
1401  switch (dim)
1402  {
1403  case 2:
1404  // at least for the 2x2 case it is reasonably simple
1405  Assert (((indices[0]==1) && (indices[1]==0)) ||
1406  ((indices[0]==0) && (indices[1]==1)),
1407  ExcInternalError());
1408  return data[2];
1409 
1410  default:
1411  // to do the rest, sort our indices before comparing
1412  {
1413  TableIndices<2> sorted_indices (indices);
1414  sorted_indices.sort ();
1415 
1416  for (unsigned int d=0, c=0; d<dim; ++d)
1417  for (unsigned int e=d+1; e<dim; ++e, ++c)
1418  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1419  return data[dim+c];
1420  Assert (false, ExcInternalError());
1421  }
1422  }
1423 
1424  static Number dummy_but_referenceable = Number();
1425  return dummy_but_referenceable;
1426  }
1427 
1428 
1429 
1430  template <int dim, typename Number>
1431  inline
1432  Number &
1433  symmetric_tensor_access (const TableIndices<4> &indices,
1435  {
1436  switch (dim)
1437  {
1438  case 1:
1439  return data[0][0];
1440 
1441  case 2:
1442  // each entry of the tensor can be
1443  // thought of as an entry in a
1444  // matrix that maps the rolled-out
1445  // rank-2 tensors into rolled-out
1446  // rank-2 tensors. this is the
1447  // format in which we store rank-4
1448  // tensors. determine which
1449  // position the present entry is
1450  // stored in
1451  {
1452  unsigned int base_index[2] ;
1453  if ((indices[0] == 0) && (indices[1] == 0))
1454  base_index[0] = 0;
1455  else if ((indices[0] == 1) && (indices[1] == 1))
1456  base_index[0] = 1;
1457  else
1458  base_index[0] = 2;
1459 
1460  if ((indices[2] == 0) && (indices[3] == 0))
1461  base_index[1] = 0;
1462  else if ((indices[2] == 1) && (indices[3] == 1))
1463  base_index[1] = 1;
1464  else
1465  base_index[1] = 2;
1466 
1467  return data[base_index[0]][base_index[1]];
1468  }
1469 
1470  case 3:
1471  // each entry of the tensor can be
1472  // thought of as an entry in a
1473  // matrix that maps the rolled-out
1474  // rank-2 tensors into rolled-out
1475  // rank-2 tensors. this is the
1476  // format in which we store rank-4
1477  // tensors. determine which
1478  // position the present entry is
1479  // stored in
1480  {
1481  unsigned int base_index[2] ;
1482  if ((indices[0] == 0) && (indices[1] == 0))
1483  base_index[0] = 0;
1484  else if ((indices[0] == 1) && (indices[1] == 1))
1485  base_index[0] = 1;
1486  else if ((indices[0] == 2) && (indices[1] == 2))
1487  base_index[0] = 2;
1488  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1489  ((indices[0] == 1) && (indices[1] == 0)))
1490  base_index[0] = 3;
1491  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1492  ((indices[0] == 2) && (indices[1] == 0)))
1493  base_index[0] = 4;
1494  else
1495  {
1496  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1497  ((indices[0] == 2) && (indices[1] == 1)),
1498  ExcInternalError());
1499  base_index[0] = 5;
1500  }
1501 
1502  if ((indices[2] == 0) && (indices[3] == 0))
1503  base_index[1] = 0;
1504  else if ((indices[2] == 1) && (indices[3] == 1))
1505  base_index[1] = 1;
1506  else if ((indices[2] == 2) && (indices[3] == 2))
1507  base_index[1] = 2;
1508  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1509  ((indices[2] == 1) && (indices[3] == 0)))
1510  base_index[1] = 3;
1511  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1512  ((indices[2] == 2) && (indices[3] == 0)))
1513  base_index[1] = 4;
1514  else
1515  {
1516  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1517  ((indices[2] == 2) && (indices[3] == 1)),
1518  ExcInternalError());
1519  base_index[1] = 5;
1520  }
1521 
1522  return data[base_index[0]][base_index[1]];
1523  }
1524 
1525  default:
1526  Assert (false, ExcNotImplemented());
1527  }
1528 
1529  static Number dummy;
1530  return dummy;
1531  }
1532 
1533 
1534  template <int dim, typename Number>
1535  inline
1536  Number
1537  symmetric_tensor_access (const TableIndices<4> &indices,
1539  {
1540  switch (dim)
1541  {
1542  case 1:
1543  return data[0][0];
1544 
1545  case 2:
1546  // each entry of the tensor can be
1547  // thought of as an entry in a
1548  // matrix that maps the rolled-out
1549  // rank-2 tensors into rolled-out
1550  // rank-2 tensors. this is the
1551  // format in which we store rank-4
1552  // tensors. determine which
1553  // position the present entry is
1554  // stored in
1555  {
1556  unsigned int base_index[2] ;
1557  if ((indices[0] == 0) && (indices[1] == 0))
1558  base_index[0] = 0;
1559  else if ((indices[0] == 1) && (indices[1] == 1))
1560  base_index[0] = 1;
1561  else
1562  base_index[0] = 2;
1563 
1564  if ((indices[2] == 0) && (indices[3] == 0))
1565  base_index[1] = 0;
1566  else if ((indices[2] == 1) && (indices[3] == 1))
1567  base_index[1] = 1;
1568  else
1569  base_index[1] = 2;
1570 
1571  return data[base_index[0]][base_index[1]];
1572  }
1573 
1574  case 3:
1575  // each entry of the tensor can be
1576  // thought of as an entry in a
1577  // matrix that maps the rolled-out
1578  // rank-2 tensors into rolled-out
1579  // rank-2 tensors. this is the
1580  // format in which we store rank-4
1581  // tensors. determine which
1582  // position the present entry is
1583  // stored in
1584  {
1585  unsigned int base_index[2] ;
1586  if ((indices[0] == 0) && (indices[1] == 0))
1587  base_index[0] = 0;
1588  else if ((indices[0] == 1) && (indices[1] == 1))
1589  base_index[0] = 1;
1590  else if ((indices[0] == 2) && (indices[1] == 2))
1591  base_index[0] = 2;
1592  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1593  ((indices[0] == 1) && (indices[1] == 0)))
1594  base_index[0] = 3;
1595  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1596  ((indices[0] == 2) && (indices[1] == 0)))
1597  base_index[0] = 4;
1598  else
1599  {
1600  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1601  ((indices[0] == 2) && (indices[1] == 1)),
1602  ExcInternalError());
1603  base_index[0] = 5;
1604  }
1605 
1606  if ((indices[2] == 0) && (indices[3] == 0))
1607  base_index[1] = 0;
1608  else if ((indices[2] == 1) && (indices[3] == 1))
1609  base_index[1] = 1;
1610  else if ((indices[2] == 2) && (indices[3] == 2))
1611  base_index[1] = 2;
1612  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1613  ((indices[2] == 1) && (indices[3] == 0)))
1614  base_index[1] = 3;
1615  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1616  ((indices[2] == 2) && (indices[3] == 0)))
1617  base_index[1] = 4;
1618  else
1619  {
1620  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1621  ((indices[2] == 2) && (indices[3] == 1)),
1622  ExcInternalError());
1623  base_index[1] = 5;
1624  }
1625 
1626  return data[base_index[0]][base_index[1]];
1627  }
1628 
1629  default:
1630  Assert (false, ExcNotImplemented());
1631  }
1632 
1633  static Number dummy;
1634  return dummy;
1635  }
1636 
1637 } // end of namespace internal
1638 
1639 
1640 
1641 template <int rank, int dim, typename Number>
1642 inline
1643 Number &
1645 {
1646  for (unsigned int r=0; r<rank; ++r)
1647  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1648  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1649 }
1650 
1651 
1652 
1653 template <int rank, int dim, typename Number>
1654 inline
1655 Number
1657 (const TableIndices<rank> &indices) const
1658 {
1659  for (unsigned int r=0; r<rank; ++r)
1660  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1661  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1662 }
1663 
1664 
1665 
1666 template <int rank, int dim, typename Number>
1667 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1668 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1669 {
1670  return
1671  internal::SymmetricTensorAccessors::
1672  Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
1673 }
1674 
1675 
1676 
1677 template <int rank, int dim, typename Number>
1678 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
1679 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
1680 {
1681  return
1682  internal::SymmetricTensorAccessors::
1683  Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
1684 }
1685 
1686 
1687 
1688 template <int rank, int dim, typename Number>
1689 inline
1690 Number
1692 {
1693  return operator()(indices);
1694 }
1695 
1696 
1697 
1698 template <int rank, int dim, typename Number>
1699 inline
1700 Number &
1702 {
1703  return operator()(indices);
1704 }
1705 
1706 
1707 
1708 
1709 namespace internal
1710 {
1711  namespace SymmetricTensor
1712  {
1713  template <int dim, typename Number>
1714  unsigned int
1715  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
1716  const unsigned int index)
1717  {
1718  return index;
1719  }
1720 
1721 
1722  template <int dim, typename Number>
1724  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
1725  const unsigned int index)
1726  {
1727  return
1730  }
1731 
1732  }
1733 }
1734 
1735 
1736 
1737 template <int rank, int dim, typename Number>
1738 inline
1739 Number
1740 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
1741 {
1742  AssertIndexRange (index, n_independent_components);
1743  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1744 }
1745 
1746 
1747 
1748 template <int rank, int dim, typename Number>
1749 inline
1750 Number &
1751 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
1752 {
1753  AssertIndexRange (index, n_independent_components);
1754  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1755 }
1756 
1757 
1758 
1759 namespace internal
1760 {
1761  template <int dim, typename Number>
1762  inline
1765  {
1766  switch (dim)
1767  {
1768  case 1:
1769  return numbers::NumberTraits<Number>::abs(data[0]);
1770 
1771  case 2:
1772  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1775 
1776  case 3:
1777  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1783 
1784  default:
1785  {
1786  typename numbers::NumberTraits<Number>::real_type return_value
1788 
1789  for (unsigned int d=0; d<dim; ++d)
1790  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
1791  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
1792  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
1793 
1794  return std::sqrt(return_value);
1795  }
1796  }
1797  }
1798 
1799 
1800 
1801  template <int dim, typename Number>
1802  inline
1805  {
1806  switch (dim)
1807  {
1808  case 1:
1809  return numbers::NumberTraits<Number>::abs (data[0][0]);
1810 
1811  default:
1812  {
1813  typename numbers::NumberTraits<Number>::real_type return_value
1815 
1816  const unsigned int n_independent_components = data.dimension;
1817 
1818  for (unsigned int i=0; i<dim; ++i)
1819  for (unsigned int j=0; j<dim; ++j)
1820  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
1821  for (unsigned int i=0; i<dim; ++i)
1822  for (unsigned int j=dim; j<n_independent_components; ++j)
1823  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1824  for (unsigned int i=dim; i<n_independent_components; ++i)
1825  for (unsigned int j=0; j<dim; ++j)
1826  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1827  for (unsigned int i=dim; i<n_independent_components; ++i)
1828  for (unsigned int j=dim; j<n_independent_components; ++j)
1829  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1830 
1831  return std::sqrt(return_value);
1832  }
1833  }
1834  }
1835 
1836 } // end of namespace internal
1837 
1838 
1839 
1840 template <int rank, int dim, typename Number>
1841 inline
1844 {
1845  return internal::compute_norm<dim,Number> (data);
1846 }
1847 
1848 
1849 
1850 namespace internal
1851 {
1852  namespace SymmetricTensor
1853  {
1854  namespace
1855  {
1856  // a function to do the unrolling from a set of indices to a
1857  // scalar index into the array in which we store the elements of
1858  // a symmetric tensor
1859  //
1860  // this function is for rank-2 tensors
1861  template <int dim>
1862  inline
1863  unsigned int
1864  component_to_unrolled_index
1865  (const TableIndices<2> &indices)
1866  {
1867  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
1868  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
1869 
1870  switch (dim)
1871  {
1872  case 1:
1873  {
1874  return 0;
1875  }
1876 
1877  case 2:
1878  {
1879  static const unsigned int table[2][2] = {{0, 2},
1880  {2, 1}
1881  };
1882  return table[indices[0]][indices[1]];
1883  }
1884 
1885  case 3:
1886  {
1887  static const unsigned int table[3][3] = {{0, 3, 4},
1888  {3, 1, 5},
1889  {4, 5, 2}
1890  };
1891  return table[indices[0]][indices[1]];
1892  }
1893 
1894  case 4:
1895  {
1896  static const unsigned int table[4][4] = {{0, 4, 5, 6},
1897  {4, 1, 7, 8},
1898  {5, 7, 2, 9},
1899  {6, 8, 9, 3}
1900  };
1901  return table[indices[0]][indices[1]];
1902  }
1903 
1904  default:
1905  // for the remainder, manually figure out the numbering
1906  {
1907  if (indices[0] == indices[1])
1908  return indices[0];
1909 
1910  TableIndices<2> sorted_indices (indices);
1911  sorted_indices.sort ();
1912 
1913  for (unsigned int d=0, c=0; d<dim; ++d)
1914  for (unsigned int e=d+1; e<dim; ++e, ++c)
1915  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1916  return dim+c;
1917 
1918  // should never get here:
1919  Assert(false, ExcInternalError());
1920  return 0;
1921  }
1922  }
1923  }
1924 
1925  // a function to do the unrolling from a set of indices to a
1926  // scalar index into the array in which we store the elements of
1927  // a symmetric tensor
1928  //
1929  // this function is for tensors of ranks not already handled
1930  // above
1931  template <int dim, int rank>
1932  inline
1933  unsigned int
1934  component_to_unrolled_index
1935  (const TableIndices<rank> &indices)
1936  {
1937  (void)indices;
1938  Assert (false, ExcNotImplemented());
1940  }
1941  }
1942  }
1943 }
1944 
1945 
1946 template <int rank, int dim, typename Number>
1947 inline
1948 unsigned int
1950 (const TableIndices<rank> &indices)
1951 {
1952  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
1953 }
1954 
1955 
1956 
1957 namespace internal
1958 {
1959  namespace SymmetricTensor
1960  {
1961  namespace
1962  {
1963  // a function to do the inverse of the unrolling from a set of
1964  // indices to a scalar index into the array in which we store
1965  // the elements of a symmetric tensor. in other words, it goes
1966  // from the scalar index into the array to a set of indices of
1967  // the tensor
1968  //
1969  // this function is for rank-2 tensors
1970  template <int dim>
1971  inline
1973  unrolled_to_component_indices
1974  (const unsigned int i,
1975  const int2type<2> &)
1976  {
1979  switch (dim)
1980  {
1981  case 1:
1982  {
1983  return TableIndices<2>(0,0);
1984  }
1985 
1986  case 2:
1987  {
1988  const TableIndices<2> table[3] =
1989  {
1990  TableIndices<2> (0,0),
1991  TableIndices<2> (1,1),
1992  TableIndices<2> (0,1)
1993  };
1994  return table[i];
1995  }
1996 
1997  case 3:
1998  {
1999  const TableIndices<2> table[6] =
2000  {
2001  TableIndices<2> (0,0),
2002  TableIndices<2> (1,1),
2003  TableIndices<2> (2,2),
2004  TableIndices<2> (0,1),
2005  TableIndices<2> (0,2),
2006  TableIndices<2> (1,2)
2007  };
2008  return table[i];
2009  }
2010 
2011  default:
2012  if (i<dim)
2013  return TableIndices<2> (i,i);
2014 
2015  for (unsigned int d=0, c=0; d<dim; ++d)
2016  for (unsigned int e=d+1; e<dim; ++e, ++c)
2017  if (c==i)
2018  return TableIndices<2>(d,e);
2019 
2020  // should never get here:
2021  Assert(false, ExcInternalError());
2022  return TableIndices<2>(0, 0);
2023  }
2024  }
2025 
2026  // a function to do the inverse of the unrolling from a set of
2027  // indices to a scalar index into the array in which we store
2028  // the elements of a symmetric tensor. in other words, it goes
2029  // from the scalar index into the array to a set of indices of
2030  // the tensor
2031  //
2032  // this function is for tensors of a rank not already handled
2033  // above
2034  template <int dim, int rank>
2035  inline
2037  unrolled_to_component_indices
2038  (const unsigned int i,
2039  const int2type<rank> &)
2040  {
2041  (void)i;
2044  Assert (false, ExcNotImplemented());
2045  return TableIndices<rank>();
2046  }
2047 
2048  }
2049  }
2050 }
2051 
2052 template <int rank, int dim, typename Number>
2053 inline
2056 (const unsigned int i)
2057 {
2058  return
2059  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2061 }
2062 
2063 
2064 
2065 template <int rank, int dim, typename Number>
2066 template <class Archive>
2067 inline
2068 void
2069 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2070 {
2071  ar &data;
2072 }
2073 
2074 
2075 #endif // DOXYGEN
2076 
2077 /* ----------------- Non-member functions operating on tensors. ------------ */
2078 
2079 
2086 template <int rank, int dim, typename Number, typename OtherNumber>
2087 inline
2090  const Tensor<rank, dim, OtherNumber> &right)
2091 {
2092  return Tensor<rank, dim, Number>(left) + right;
2093 }
2094 
2095 
2102 template <int rank, int dim, typename Number, typename OtherNumber>
2103 inline
2107 {
2108  return left + Tensor<rank, dim, OtherNumber>(right);
2109 }
2110 
2111 
2118 template <int rank, int dim, typename Number, typename OtherNumber>
2119 inline
2122  const Tensor<rank, dim, OtherNumber> &right)
2123 {
2124  return Tensor<rank, dim, Number>(left) - right;
2125 }
2126 
2127 
2134 template <int rank, int dim, typename Number, typename OtherNumber>
2135 inline
2139 {
2140  return left - Tensor<rank, dim, OtherNumber>(right);
2141 }
2142 
2143 
2144 
2158 template <int dim, typename Number>
2159 inline
2161 {
2162  switch (dim)
2163  {
2164  case 1:
2165  return t.data[0];
2166  case 2:
2167  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2168  case 3:
2169  {
2170  // in analogy to general tensors, but
2171  // there's something to be simplified for
2172  // the present case
2173  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2174  return ( tmp + tmp
2175  +t.data[0]*t.data[1]*t.data[2]
2176  -t.data[0]*t.data[5]*t.data[5]
2177  -t.data[1]*t.data[4]*t.data[4]
2178  -t.data[2]*t.data[3]*t.data[3]);
2179  }
2180  default:
2181  Assert (false, ExcNotImplemented());
2182  return 0;
2183  }
2184 }
2185 
2186 
2187 
2197 template <int dim, typename Number>
2198 inline
2200 {
2201  return determinant (t);
2202 }
2203 
2204 
2205 
2213 template <int dim, typename Number>
2215 {
2216  Number t = d.data[0];
2217  for (unsigned int i=1; i<dim; ++i)
2218  t += d.data[i];
2219  return t;
2220 }
2221 
2222 
2232 template <int dim, typename Number>
2233 inline
2235 {
2236  return trace (t);
2237 }
2238 
2239 
2252 template <typename Number>
2253 inline
2255 {
2256  return 0;
2257 }
2258 
2259 
2260 
2281 template <typename Number>
2282 inline
2284 {
2285  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2286 }
2287 
2288 
2289 
2299 template <typename Number>
2300 inline
2302 {
2303  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2304  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2305 }
2306 
2307 
2308 
2309 
2319 template <int rank, int dim, typename Number>
2320 inline
2323 {
2324  return t;
2325 }
2326 
2327 
2328 
2338 template <int dim, typename Number>
2339 inline
2342 {
2344 
2345  // subtract scaled trace from the diagonal
2346  const Number tr = trace(t) / dim;
2347  for (unsigned int i=0; i<dim; ++i)
2348  tmp.data[i] -= tr;
2349 
2350  return tmp;
2351 }
2352 
2353 
2354 
2362 template <int dim, typename Number>
2363 inline
2365 unit_symmetric_tensor ()
2366 {
2367  // create a default constructed matrix filled with
2368  // zeros, then set the diagonal elements to one
2370  switch (dim)
2371  {
2372  case 1:
2373  tmp.data[0] = 1;
2374  break;
2375  case 2:
2376  tmp.data[0] = tmp.data[1] = 1;
2377  break;
2378  case 3:
2379  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
2380  break;
2381  default:
2382  for (unsigned int d=0; d<dim; ++d)
2383  tmp.data[d] = 1;
2384  }
2385  return tmp;
2386 }
2387 
2388 
2389 
2398 template <int dim>
2399 inline
2401 unit_symmetric_tensor ()
2402 {
2403  return unit_symmetric_tensor<dim,double>();
2404 }
2405 
2406 
2407 
2422 template <int dim, typename Number>
2423 inline
2425 deviator_tensor ()
2426 {
2428 
2429  // fill the elements treating the diagonal
2430  for (unsigned int i=0; i<dim; ++i)
2431  for (unsigned int j=0; j<dim; ++j)
2432  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
2433 
2434  // then fill the ones that copy over the
2435  // non-diagonal elements. note that during
2436  // the double-contraction, we handle the
2437  // off-diagonal elements twice, so simply
2438  // copying requires a weight of 1/2
2439  for (unsigned int i=dim;
2440  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2441  ++i)
2442  tmp.data[i][i] = 0.5;
2443 
2444  return tmp;
2445 }
2446 
2447 
2448 
2463 template <int dim>
2464 inline
2466 deviator_tensor ()
2467 {
2468  return deviator_tensor<dim,double>();
2469 }
2470 
2471 
2472 
2495 template <int dim, typename Number>
2496 inline
2498 identity_tensor ()
2499 {
2501 
2502  // fill the elements treating the diagonal
2503  for (unsigned int i=0; i<dim; ++i)
2504  tmp.data[i][i] = 1;
2505 
2506  // then fill the ones that copy over the
2507  // non-diagonal elements. note that during
2508  // the double-contraction, we handle the
2509  // off-diagonal elements twice, so simply
2510  // copying requires a weight of 1/2
2511  for (unsigned int i=dim;
2512  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2513  ++i)
2514  tmp.data[i][i] = 0.5;
2515 
2516  return tmp;
2517 }
2518 
2519 
2520 
2542 template <int dim>
2543 inline
2545 identity_tensor ()
2546 {
2547  return identity_tensor<dim,double>();
2548 }
2549 
2550 
2551 
2562 template <int dim, typename Number>
2563 inline
2566 {
2567  // if desired, take over the
2568  // inversion of a 4x4 tensor
2569  // from the FullMatrix
2570  AssertThrow (false, ExcNotImplemented());
2571 
2573 }
2574 
2575 
2576 
2577 #ifndef DOXYGEN
2578 
2579 template <typename Number>
2580 inline
2582 invert (const SymmetricTensor<2,1,Number> &t)
2583 {
2585 
2586  tmp[0][0] = 1.0/t[0][0];
2587 
2588  return tmp;
2589 }
2590 
2591 
2592 
2593 template <typename Number>
2594 inline
2596 invert (const SymmetricTensor<2,2,Number> &t)
2597 {
2599 
2600  // Sympy result: ([
2601  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
2602  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
2603  const TableIndices<2> idx_00 (0,0);
2604  const TableIndices<2> idx_01 (0,1);
2605  const TableIndices<2> idx_11 (1,1);
2606  const Number inv_det_t
2607  = 1.0/(t[idx_00]*t[idx_11]
2608  - t[idx_01]*t[idx_01]);
2609  tmp[idx_00] = t[idx_11];
2610  tmp[idx_01] = -t[idx_01];
2611  tmp[idx_11] = t[idx_00];
2612  tmp *= inv_det_t;
2613 
2614  return tmp;
2615 }
2616 
2617 
2618 
2619 template <typename Number>
2620 inline
2622 invert (const SymmetricTensor<2,3,Number> &t)
2623 {
2625 
2626  // Sympy result: ([
2627  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2628  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2629  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
2630  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2631  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2632  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
2633  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2634  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
2635  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
2636  const TableIndices<2> idx_00 (0,0);
2637  const TableIndices<2> idx_01 (0,1);
2638  const TableIndices<2> idx_02 (0,2);
2639  const TableIndices<2> idx_11 (1,1);
2640  const TableIndices<2> idx_12 (1,2);
2641  const TableIndices<2> idx_22 (2,2);
2642  const Number inv_det_t
2643  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
2644  - t[idx_00]*t[idx_12]*t[idx_12]
2645  - t[idx_01]*t[idx_01]*t[idx_22]
2646  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
2647  - t[idx_02]*t[idx_02]*t[idx_11]);
2648  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
2649  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
2650  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
2651  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
2652  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
2653  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
2654  tmp *= inv_det_t;
2655 
2656  return tmp;
2657 }
2658 
2659 #endif /* DOXYGEN */
2660 
2661 
2662 
2676 template <int dim, typename Number>
2677 inline
2680 {
2682  switch (dim)
2683  {
2684  case 1:
2685  tmp.data[0][0] = 1./t.data[0][0];
2686  break;
2687  case 2:
2688 
2689  // inverting this tensor is a little more
2690  // complicated than necessary, since we
2691  // store the data of 't' as a 3x3 matrix
2692  // t.data, but the product between a rank-4
2693  // and a rank-2 tensor is really not the
2694  // product between this matrix and the
2695  // 3-vector of a rhs, but rather
2696  //
2697  // B.vec = t.data * mult * A.vec
2698  //
2699  // where mult is a 3x3 matrix with
2700  // entries [[1,0,0],[0,1,0],[0,0,2]] to
2701  // capture the fact that we need to add up
2702  // both the c_ij12*a_12 and the c_ij21*a_21
2703  // terms
2704  //
2705  // in addition, in this scheme, the
2706  // identity tensor has the matrix
2707  // representation mult^-1.
2708  //
2709  // the inverse of 't' therefore has the
2710  // matrix representation
2711  //
2712  // inv.data = mult^-1 * t.data^-1 * mult^-1
2713  //
2714  // in order to compute it, let's first
2715  // compute the inverse of t.data and put it
2716  // into tmp.data; at the end of the
2717  // function we then scale the last row and
2718  // column of the inverse by 1/2,
2719  // corresponding to the left and right
2720  // multiplication with mult^-1
2721  {
2722  const Number t4 = t.data[0][0]*t.data[1][1],
2723  t6 = t.data[0][0]*t.data[1][2],
2724  t8 = t.data[0][1]*t.data[1][0],
2725  t00 = t.data[0][2]*t.data[1][0],
2726  t01 = t.data[0][1]*t.data[2][0],
2727  t04 = t.data[0][2]*t.data[2][0],
2728  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
2729  t8*t.data[2][2]+t00*t.data[2][1]+
2730  t01*t.data[1][2]-t04*t.data[1][1]);
2731  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
2732  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
2733  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
2734  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
2735  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
2736  tmp.data[1][2] = -(t6-t00)*t07;
2737  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
2738  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
2739  tmp.data[2][2] = (t4-t8)*t07;
2740 
2741  // scale last row and column as mentioned
2742  // above
2743  tmp.data[2][0] /= 2;
2744  tmp.data[2][1] /= 2;
2745  tmp.data[0][2] /= 2;
2746  tmp.data[1][2] /= 2;
2747  tmp.data[2][2] /= 4;
2748  }
2749  break;
2750  default:
2751  Assert (false, ExcNotImplemented());
2752  }
2753  return tmp;
2754 }
2755 
2756 
2757 
2771 template <>
2773 invert (const SymmetricTensor<4,3,double> &t);
2774 // this function is implemented in the .cc file for double data types
2775 
2776 
2777 
2792 template <int dim, typename Number>
2793 inline
2797 {
2799 
2800  // fill only the elements really needed
2801  for (unsigned int i=0; i<dim; ++i)
2802  for (unsigned int j=i; j<dim; ++j)
2803  for (unsigned int k=0; k<dim; ++k)
2804  for (unsigned int l=k; l<dim; ++l)
2805  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
2806 
2807  return tmp;
2808 }
2809 
2810 
2811 
2820 template <int dim,typename Number>
2821 inline
2824 {
2825  Number array[(dim*dim+dim)/2];
2826  for (unsigned int d=0; d<dim; ++d)
2827  array[d] = t[d][d];
2828  for (unsigned int d=0, c=0; d<dim; ++d)
2829  for (unsigned int e=d+1; e<dim; ++e, ++c)
2830  array[dim+c] = (t[d][e]+t[e][d])*0.5;
2831  return SymmetricTensor<2,dim,Number>(array);
2832 }
2833 
2834 
2835 
2843 template <int rank, int dim, typename Number>
2844 inline
2847  const Number factor)
2848 {
2850  tt *= factor;
2851  return tt;
2852 }
2853 
2854 
2855 
2863 template <int rank, int dim, typename Number>
2864 inline
2866 operator * (const Number factor,
2868 {
2869  // simply forward to the other operator
2870  return t*factor;
2871 }
2872 
2873 
2874 
2900 template <int rank, int dim, typename Number, typename OtherNumber>
2901 inline
2904  const OtherNumber factor)
2905 {
2906  // form the product. we have to convert the two factors into the final
2907  // type via explicit casts because, for awkward reasons, the C++
2908  // standard committee saw it fit to not define an
2909  // operator*(float,std::complex<double>)
2910  // (as well as with switched arguments and double<->float).
2911  typedef typename ProductType<Number,OtherNumber>::type product_type;
2913  // we used to shorten the following by 'tt *= product_type(factor);'
2914  // which requires that a converting constructor
2915  // 'product_type::product_type(const OtherNumber) is defined.
2916  // however, a user-defined constructor is not allowed for aggregates,
2917  // e.g. VectorizedArray. therefore, we work around this issue using a
2918  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
2919  // which we assume to be defined.
2920  product_type new_factor;
2921  new_factor = factor;
2922  tt *= new_factor;
2923  return tt;
2924 }
2925 
2926 
2927 
2936 template <int rank, int dim, typename Number, typename OtherNumber>
2937 inline
2939 operator * (const Number factor,
2941 {
2942  // simply forward to the other operator with switched arguments
2943  return (t*factor);
2944 }
2945 
2946 
2947 
2953 template <int rank, int dim, typename Number>
2954 inline
2957  const Number factor)
2958 {
2960  tt /= factor;
2961  return tt;
2962 }
2963 
2964 
2965 
2972 template <int rank, int dim>
2973 inline
2975 operator * (const SymmetricTensor<rank,dim> &t,
2976  const double factor)
2977 {
2979  tt *= factor;
2980  return tt;
2981 }
2982 
2983 
2984 
2991 template <int rank, int dim>
2992 inline
2994 operator * (const double factor,
2995  const SymmetricTensor<rank,dim> &t)
2996 {
2998  tt *= factor;
2999  return tt;
3000 }
3001 
3002 
3003 
3009 template <int rank, int dim>
3010 inline
3012 operator / (const SymmetricTensor<rank,dim> &t,
3013  const double factor)
3014 {
3016  tt /= factor;
3017  return tt;
3018 }
3019 
3029 template <int dim, typename Number>
3030 inline
3031 Number
3034 {
3035  return (t1*t2);
3036 }
3037 
3038 
3048 template <int dim, typename Number>
3049 inline
3050 Number
3052  const Tensor<2,dim,Number> &t2)
3053 {
3054  Number s = 0;
3055  for (unsigned int i=0; i<dim; ++i)
3056  for (unsigned int j=0; j<dim; ++j)
3057  s += t1[i][j] * t2[i][j];
3058  return s;
3059 }
3060 
3061 
3071 template <int dim, typename Number>
3072 inline
3073 Number
3076 {
3077  return scalar_product(t2, t1);
3078 }
3079 
3080 
3096 template <typename Number>
3097 inline
3098 void
3100  const SymmetricTensor<4,1,Number> &t,
3101  const SymmetricTensor<2,1,Number> &s)
3102 {
3103  tmp[0][0] = t[0][0][0][0] * s[0][0];
3104 }
3105 
3106 
3107 
3123 template <typename Number>
3124 inline
3125 void
3127  const SymmetricTensor<2,1,Number> &s,
3128  const SymmetricTensor<4,1,Number> &t)
3129 {
3130  tmp[0][0] = t[0][0][0][0] * s[0][0];
3131 }
3132 
3133 
3134 
3149 template <typename Number>
3150 inline
3151 void
3153  const SymmetricTensor<4,2,Number> &t,
3154  const SymmetricTensor<2,2,Number> &s)
3155 {
3156  const unsigned int dim = 2;
3157 
3158  for (unsigned int i=0; i<dim; ++i)
3159  for (unsigned int j=i; j<dim; ++j)
3160  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3161  t[i][j][1][1] * s[1][1] +
3162  2 * t[i][j][0][1] * s[0][1];
3163 }
3164 
3165 
3166 
3182 template <typename Number>
3183 inline
3184 void
3186  const SymmetricTensor<2,2,Number> &s,
3187  const SymmetricTensor<4,2,Number> &t)
3188 {
3189  const unsigned int dim = 2;
3190 
3191  for (unsigned int i=0; i<dim; ++i)
3192  for (unsigned int j=i; j<dim; ++j)
3193  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3194  s[1][1] * t[1][1][i][j] +
3195  2 * s[0][1] * t[0][1][i][j];
3196 }
3197 
3198 
3199 
3215 template <typename Number>
3216 inline
3217 void
3219  const SymmetricTensor<4,3,Number> &t,
3220  const SymmetricTensor<2,3,Number> &s)
3221 {
3222  const unsigned int dim = 3;
3223 
3224  for (unsigned int i=0; i<dim; ++i)
3225  for (unsigned int j=i; j<dim; ++j)
3226  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3227  t[i][j][1][1] * s[1][1] +
3228  t[i][j][2][2] * s[2][2] +
3229  2 * t[i][j][0][1] * s[0][1] +
3230  2 * t[i][j][0][2] * s[0][2] +
3231  2 * t[i][j][1][2] * s[1][2];
3232 }
3233 
3234 
3235 
3251 template <typename Number>
3252 inline
3253 void
3255  const SymmetricTensor<2,3,Number> &s,
3256  const SymmetricTensor<4,3,Number> &t)
3257 {
3258  const unsigned int dim = 3;
3259 
3260  for (unsigned int i=0; i<dim; ++i)
3261  for (unsigned int j=i; j<dim; ++j)
3262  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3263  s[1][1] * t[1][1][i][j] +
3264  s[2][2] * t[2][2][i][j] +
3265  2 * s[0][1] * t[0][1][i][j] +
3266  2 * s[0][2] * t[0][2][i][j] +
3267  2 * s[1][2] * t[1][2][i][j];
3268 }
3269 
3270 
3271 
3279 template <int dim, typename Number>
3281 operator * (const SymmetricTensor<2,dim,Number> &src1,
3282  const Tensor<1,dim,Number> &src2)
3283 {
3284  Tensor<1,dim,Number> dest;
3285  for (unsigned int i=0; i<dim; ++i)
3286  for (unsigned int j=0; j<dim; ++j)
3287  dest[i] += src1[i][j] * src2[j];
3288  return dest;
3289 }
3290 
3291 
3299 template <int dim, typename Number>
3301 operator * (const Tensor<1,dim,Number> &src1,
3302  const SymmetricTensor<2,dim,Number> &src2)
3303 {
3304  // this is easy for symmetric tensors:
3305  return src2 * src1;
3306 }
3307 
3308 
3318 template <int dim, typename Number>
3319 inline
3320 std::ostream &operator << (std::ostream &out,
3322 {
3323  //make out lives a bit simpler by outputing
3324  //the tensor through the operator for the
3325  //general Tensor class
3327 
3328  for (unsigned int i=0; i<dim; ++i)
3329  for (unsigned int j=0; j<dim; ++j)
3330  tt[i][j] = t[i][j];
3331 
3332  return out << tt;
3333 }
3334 
3335 
3336 
3346 template <int dim, typename Number>
3347 inline
3348 std::ostream &operator << (std::ostream &out,
3350 {
3351  //make out lives a bit simpler by outputing
3352  //the tensor through the operator for the
3353  //general Tensor class
3355 
3356  for (unsigned int i=0; i<dim; ++i)
3357  for (unsigned int j=0; j<dim; ++j)
3358  for (unsigned int k=0; k<dim; ++k)
3359  for (unsigned int l=0; l<dim; ++l)
3360  tt[i][j][k][l] = t[i][j][k][l];
3361 
3362  return out << tt;
3363 }
3364 
3365 
3366 DEAL_II_NAMESPACE_CLOSE
3367 
3368 #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()
static const unsigned int invalid_unsigned_int
Definition: types.h:170
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)
void double_contract(SymmetricTensor< 2, 2, Number > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, Number > &s)
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)
void double_contract(SymmetricTensor< 2, 3, Number > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, Number > &t)
void double_contract(SymmetricTensor< 2, 1, Number > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, Number > &s)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1170
static TableIndices< rank > unrolled_to_component_indices(const unsigned int i)
base_tensor_type data
void double_contract(SymmetricTensor< 2, 1, Number > &tmp, const SymmetricTensor< 2, 1, Number > &s, const SymmetricTensor< 4, 1, Number > &t)
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
#define AssertThrow(cond, exc)
Definition: exceptions.h:369
static real_type abs(const number &x)
Definition: numbers.h:354
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void serialize(Archive &ar, const unsigned int version)
static const unsigned int dimension
static unsigned int component_to_unrolled_index(const TableIndices< rank > &indices)
void double_contract(SymmetricTensor< 2, 2, Number > &tmp, const SymmetricTensor< 2, 2, Number > &s, const SymmetricTensor< 4, 2, Number > &t)
static::ExceptionBase & ExcMessage(std::string arg1)
static std::size_t memory_consumption()
SymmetricTensor operator+(const SymmetricTensor &s) const
Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor & operator-=(const SymmetricTensor &)
friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:313
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)
internal::SymmetricTensorAccessors::double_contraction_result< rank, 2, dim, Number >::type operator*(const SymmetricTensor< 2, dim, Number > &s) const
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)
Number scalar_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)
void double_contract(SymmetricTensor< 2, 3, Number > &tmp, const SymmetricTensor< 4, 3, Number > &t, const SymmetricTensor< 2, 3, Number > &s)
Number scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, Number > &t2)
SymmetricTensor operator-() const
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: mpi.h:41
SymmetricTensor & operator/=(const Number factor)
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< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
SymmetricTensor & operator+=(const SymmetricTensor &)
SymmetricTensor & operator*=(const Number factor)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
friend SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
StreamType & operator<<(StreamType &s, UpdateFlags u)
Number scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
SymmetricTensor & operator=(const Number d)
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static::ExceptionBase & ExcInternalError()