Reference documentation for deal.II version Git 4a360c6 2017-05-26 05:30:59 +0200
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 () = default;
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  Assert (rank == 2, ExcNotImplemented());
930  switch (dim)
931  {
932  case 2:
933  Assert (t[0][1] == t[1][0], ExcInternalError());
934 
935  data[0] = t[0][0];
936  data[1] = t[1][1];
937  data[2] = t[0][1];
938 
939  break;
940  case 3:
941  Assert (t[0][1] == t[1][0], ExcInternalError());
942  Assert (t[0][2] == t[2][0], ExcInternalError());
943  Assert (t[1][2] == t[2][1], ExcInternalError());
944 
945  data[0] = t[0][0];
946  data[1] = t[1][1];
947  data[2] = t[2][2];
948  data[3] = t[0][1];
949  data[4] = t[0][2];
950  data[5] = t[1][2];
951 
952  break;
953  default:
954  for (unsigned int d=0; d<dim; ++d)
955  for (unsigned int e=0; e<d; ++e)
956  Assert(t[d][e] == t[e][d], ExcInternalError());
957 
958  for (unsigned int d=0; d<dim; ++d)
959  data[d] = t[d][d];
960 
961  for (unsigned int d=0, c=0; d<dim; ++d)
962  for (unsigned int e=d+1; e<dim; ++e, ++c)
963  data[dim+c] = t[d][e];
964  }
965 }
966 
967 
968 
969 template <int rank, int dim, typename Number>
970 template <typename OtherNumber>
971 inline
974 {
975  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
976  data[i] = initializer.data[i];
977 }
978 
979 
980 
981 
982 template <int rank, int dim, typename Number>
983 inline
984 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
985  :
986  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
987 {
988  // ensure that the reinterpret_cast above actually works
989  Assert (sizeof(typename base_tensor_type::array_type)
990  == sizeof(array),
991  ExcInternalError());
992 }
993 
994 
995 
996 template <int rank, int dim, typename Number>
997 inline
1000 {
1001  Assert (d==Number(), ExcMessage ("Only assignment with zero is allowed"));
1002  (void) d;
1003 
1004  data = 0;
1005 
1006  return *this;
1007 }
1008 
1009 
1010 namespace internal
1011 {
1012  namespace SymmetricTensor
1013  {
1014  template <int dim, typename Number>
1016  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1017  {
1019 
1020  // diagonal entries are stored first
1021  for (unsigned int d=0; d<dim; ++d)
1022  t[d][d] = s.access_raw_entry(d);
1023 
1024  // off-diagonal entries come next, row by row
1025  for (unsigned int d=0, c=0; d<dim; ++d)
1026  for (unsigned int e=d+1; e<dim; ++e, ++c)
1027  {
1028  t[d][e] = s.access_raw_entry(dim+c);
1029  t[e][d] = s.access_raw_entry(dim+c);
1030  }
1031  return t;
1032  }
1033 
1034 
1035  template <int dim, typename Number>
1037  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1038  {
1039  // utilize the symmetry properties of SymmetricTensor<4,dim>
1040  // discussed in the class documentation to avoid accessing all
1041  // independent elements of the input tensor more than once
1043 
1044  for (unsigned int i=0; i<dim; ++i)
1045  for (unsigned int j=i; j<dim; ++j)
1046  for (unsigned int k=0; k<dim; ++k)
1047  for (unsigned int l=k; l<dim; ++l)
1048  t[TableIndices<4>(i,j,k,l)]
1049  = t[TableIndices<4>(i,j,l,k)]
1050  = t[TableIndices<4>(j,i,k,l)]
1051  = t[TableIndices<4>(j,i,l,k)]
1052  = st[TableIndices<4>(i,j,k,l)];
1053 
1054  return t;
1055  }
1056  }
1057 }
1058 
1059 
1060 
1061 template <int rank, int dim, typename Number>
1062 inline
1064 operator Tensor<rank,dim,Number> () const
1065 {
1066  return internal::SymmetricTensor::convert_to_tensor (*this);
1067 }
1068 
1069 
1070 
1071 template <int rank, int dim, typename Number>
1072 inline
1073 bool
1075 (const SymmetricTensor<rank,dim,Number> &t) const
1076 {
1077  return data == t.data;
1078 }
1079 
1080 
1081 
1082 template <int rank, int dim, typename Number>
1083 inline
1084 bool
1085 SymmetricTensor<rank,dim,Number>::operator !=
1086 (const SymmetricTensor<rank,dim,Number> &t) const
1087 {
1088  return data != t.data;
1089 }
1090 
1091 
1092 
1093 template <int rank, int dim, typename Number>
1094 inline
1096 SymmetricTensor<rank,dim,Number>::operator +=
1098 {
1099  data += t.data;
1100  return *this;
1101 }
1102 
1103 
1104 
1105 template <int rank, int dim, typename Number>
1106 inline
1108 SymmetricTensor<rank,dim,Number>::operator -=
1110 {
1111  data -= t.data;
1112  return *this;
1113 }
1114 
1115 
1116 
1117 template <int rank, int dim, typename Number>
1118 inline
1121 {
1122  data *= d;
1123  return *this;
1124 }
1125 
1126 
1127 
1128 template <int rank, int dim, typename Number>
1129 inline
1132 {
1133  data /= d;
1134  return *this;
1135 }
1136 
1137 
1138 
1139 template <int rank, int dim, typename Number>
1140 inline
1143 {
1144  SymmetricTensor tmp = *this;
1145  tmp.data += t.data;
1146  return tmp;
1147 }
1148 
1149 
1150 
1151 template <int rank, int dim, typename Number>
1152 inline
1155 {
1156  SymmetricTensor tmp = *this;
1157  tmp.data -= t.data;
1158  return tmp;
1159 }
1160 
1161 
1162 
1163 template <int rank, int dim, typename Number>
1164 inline
1167 {
1168  SymmetricTensor tmp = *this;
1169  tmp.data = -tmp.data;
1170  return tmp;
1171 }
1172 
1173 
1174 
1175 template <int rank, int dim, typename Number>
1176 inline
1177 void
1179 {
1180  data.clear ();
1181 }
1182 
1183 
1184 
1185 template <int rank, int dim, typename Number>
1186 inline
1187 std::size_t
1189 {
1190  // all memory consists of statically allocated memory of the current
1191  // object, no pointers
1192  return sizeof(SymmetricTensor<rank,dim,Number>);
1193 }
1194 
1195 
1196 
1197 namespace internal
1198 {
1199 
1200  template <int dim, typename Number>
1201  inline
1202  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
1203  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1205  {
1206  switch (dim)
1207  {
1208  case 1:
1209  return data[0] * sdata[0];
1210  default:
1211  // Start with the non-diagonal part to avoid some multiplications by
1212  // 2.
1213  Number sum = data[dim] * sdata[dim];
1214  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1215  sum += data[d] * sdata[d];
1216  sum += sum; // sum = sum * 2.;
1217 
1218  // Now add the contributions from the diagonal
1219  for (unsigned int d=0; d<dim; ++d)
1220  sum += data[d] * sdata[d];
1221  return sum;
1222  }
1223  }
1224 
1225 
1226 
1227  template <int dim, typename Number>
1228  inline
1229  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
1230  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1232  {
1233  const unsigned int data_dim =
1235  Number tmp [data_dim];
1236  for (unsigned int i=0; i<data_dim; ++i)
1237  tmp[i] = perform_double_contraction<dim,Number>(data[i], sdata);
1238  return ::SymmetricTensor<2,dim,Number>(tmp);
1239  }
1240 
1241 
1242 
1243  template <int dim, typename Number>
1244  inline
1246  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1248  {
1250  for (unsigned int i=0; i<tmp.dimension; ++i)
1251  {
1252  // Start with the non-diagonal part
1253  Number sum = data[dim] * sdata[dim][i];
1254  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1255  sum += data[d] * sdata[d][i];
1256  sum += sum; // sum = sum * 2.;
1257 
1258  // Now add the contributions from the diagonal
1259  for (unsigned int d=0; d<dim; ++d)
1260  sum += data[d] * sdata[d][i];
1261  tmp[i] = sum;
1262  }
1263  return tmp;
1264  }
1265 
1266 
1267 
1268  template <int dim, typename Number>
1269  inline
1271  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1273  {
1274  const unsigned int data_dim =
1277  for (unsigned int i=0; i<data_dim; ++i)
1278  for (unsigned int j=0; j<data_dim; ++j)
1279  {
1280  // Start with the non-diagonal part
1281  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1282  tmp[i][j] += data[i][d] * sdata[d][j];
1283  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1284 
1285  // Now add the contributions from the diagonal
1286  for (unsigned int d=0; d<dim; ++d)
1287  tmp[i][j] += data[i][d] * sdata[d][j];
1288  }
1289  return tmp;
1290  }
1291 
1292 } // end of namespace internal
1293 
1294 
1295 
1296 template <int rank, int dim, typename Number>
1297 inline
1300 {
1301  // need to have two different function calls
1302  // because a scalar and rank-2 tensor are not
1303  // the same data type (see internal function
1304  // above)
1305  return internal::perform_double_contraction<dim,Number> (data, s.data);
1306 }
1307 
1308 
1309 
1310 template <int rank, int dim, typename Number>
1311 inline
1314 {
1317  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
1318  return tmp;
1319 }
1320 
1321 
1322 
1323 // internal namespace to switch between the
1324 // access of different tensors. There used to
1325 // be explicit instantiations before for
1326 // different ranks and dimensions, but since
1327 // we now allow for templates on the data
1328 // type, and since we cannot partially
1329 // specialize the implementation, this got
1330 // into a separate namespace
1331 namespace internal
1332 {
1333  template <int dim, typename Number>
1334  inline
1335  Number &
1336  symmetric_tensor_access (const TableIndices<2> &indices,
1338  {
1339  // 1d is very simple and done first
1340  if (dim == 1)
1341  return data[0];
1342 
1343  // first treat the main diagonal elements, which are stored consecutively
1344  // at the beginning
1345  if (indices[0] == indices[1])
1346  return data[indices[0]];
1347 
1348  // the rest is messier and requires a few switches.
1349  switch (dim)
1350  {
1351  case 2:
1352  // at least for the 2x2 case it is reasonably simple
1353  Assert (((indices[0]==1) && (indices[1]==0)) ||
1354  ((indices[0]==0) && (indices[1]==1)),
1355  ExcInternalError());
1356  return data[2];
1357 
1358  default:
1359  // to do the rest, sort our indices before comparing
1360  {
1361  TableIndices<2> sorted_indices (indices);
1362  sorted_indices.sort ();
1363 
1364  for (unsigned int d=0, c=0; d<dim; ++d)
1365  for (unsigned int e=d+1; e<dim; ++e, ++c)
1366  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1367  return data[dim+c];
1368  Assert (false, ExcInternalError());
1369  }
1370  }
1371 
1372  static Number dummy_but_referenceable = Number();
1373  return dummy_but_referenceable;
1374  }
1375 
1376 
1377 
1378  template <int dim, typename Number>
1379  inline
1380  Number
1381  symmetric_tensor_access (const TableIndices<2> &indices,
1383  {
1384  // 1d is very simple and done first
1385  if (dim == 1)
1386  return data[0];
1387 
1388  // first treat the main diagonal elements, which are stored consecutively
1389  // at the beginning
1390  if (indices[0] == indices[1])
1391  return data[indices[0]];
1392 
1393  // the rest is messier and requires a few switches.
1394  switch (dim)
1395  {
1396  case 2:
1397  // at least for the 2x2 case it is reasonably simple
1398  Assert (((indices[0]==1) && (indices[1]==0)) ||
1399  ((indices[0]==0) && (indices[1]==1)),
1400  ExcInternalError());
1401  return data[2];
1402 
1403  default:
1404  // to do the rest, sort our indices before comparing
1405  {
1406  TableIndices<2> sorted_indices (indices);
1407  sorted_indices.sort ();
1408 
1409  for (unsigned int d=0, c=0; d<dim; ++d)
1410  for (unsigned int e=d+1; e<dim; ++e, ++c)
1411  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1412  return data[dim+c];
1413  Assert (false, ExcInternalError());
1414  }
1415  }
1416 
1417  static Number dummy_but_referenceable = Number();
1418  return dummy_but_referenceable;
1419  }
1420 
1421 
1422 
1423  template <int dim, typename Number>
1424  inline
1425  Number &
1426  symmetric_tensor_access (const TableIndices<4> &indices,
1428  {
1429  switch (dim)
1430  {
1431  case 1:
1432  return data[0][0];
1433 
1434  case 2:
1435  // each entry of the tensor can be
1436  // thought of as an entry in a
1437  // matrix that maps the rolled-out
1438  // rank-2 tensors into rolled-out
1439  // rank-2 tensors. this is the
1440  // format in which we store rank-4
1441  // tensors. determine which
1442  // position the present entry is
1443  // stored in
1444  {
1445  unsigned int base_index[2] ;
1446  if ((indices[0] == 0) && (indices[1] == 0))
1447  base_index[0] = 0;
1448  else if ((indices[0] == 1) && (indices[1] == 1))
1449  base_index[0] = 1;
1450  else
1451  base_index[0] = 2;
1452 
1453  if ((indices[2] == 0) && (indices[3] == 0))
1454  base_index[1] = 0;
1455  else if ((indices[2] == 1) && (indices[3] == 1))
1456  base_index[1] = 1;
1457  else
1458  base_index[1] = 2;
1459 
1460  return data[base_index[0]][base_index[1]];
1461  }
1462 
1463  case 3:
1464  // each entry of the tensor can be
1465  // thought of as an entry in a
1466  // matrix that maps the rolled-out
1467  // rank-2 tensors into rolled-out
1468  // rank-2 tensors. this is the
1469  // format in which we store rank-4
1470  // tensors. determine which
1471  // position the present entry is
1472  // stored in
1473  {
1474  unsigned int base_index[2] ;
1475  if ((indices[0] == 0) && (indices[1] == 0))
1476  base_index[0] = 0;
1477  else if ((indices[0] == 1) && (indices[1] == 1))
1478  base_index[0] = 1;
1479  else if ((indices[0] == 2) && (indices[1] == 2))
1480  base_index[0] = 2;
1481  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1482  ((indices[0] == 1) && (indices[1] == 0)))
1483  base_index[0] = 3;
1484  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1485  ((indices[0] == 2) && (indices[1] == 0)))
1486  base_index[0] = 4;
1487  else
1488  {
1489  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1490  ((indices[0] == 2) && (indices[1] == 1)),
1491  ExcInternalError());
1492  base_index[0] = 5;
1493  }
1494 
1495  if ((indices[2] == 0) && (indices[3] == 0))
1496  base_index[1] = 0;
1497  else if ((indices[2] == 1) && (indices[3] == 1))
1498  base_index[1] = 1;
1499  else if ((indices[2] == 2) && (indices[3] == 2))
1500  base_index[1] = 2;
1501  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1502  ((indices[2] == 1) && (indices[3] == 0)))
1503  base_index[1] = 3;
1504  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1505  ((indices[2] == 2) && (indices[3] == 0)))
1506  base_index[1] = 4;
1507  else
1508  {
1509  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1510  ((indices[2] == 2) && (indices[3] == 1)),
1511  ExcInternalError());
1512  base_index[1] = 5;
1513  }
1514 
1515  return data[base_index[0]][base_index[1]];
1516  }
1517 
1518  default:
1519  Assert (false, ExcNotImplemented());
1520  }
1521 
1522  static Number dummy;
1523  return dummy;
1524  }
1525 
1526 
1527  template <int dim, typename Number>
1528  inline
1529  Number
1530  symmetric_tensor_access (const TableIndices<4> &indices,
1532  {
1533  switch (dim)
1534  {
1535  case 1:
1536  return data[0][0];
1537 
1538  case 2:
1539  // each entry of the tensor can be
1540  // thought of as an entry in a
1541  // matrix that maps the rolled-out
1542  // rank-2 tensors into rolled-out
1543  // rank-2 tensors. this is the
1544  // format in which we store rank-4
1545  // tensors. determine which
1546  // position the present entry is
1547  // stored in
1548  {
1549  unsigned int base_index[2] ;
1550  if ((indices[0] == 0) && (indices[1] == 0))
1551  base_index[0] = 0;
1552  else if ((indices[0] == 1) && (indices[1] == 1))
1553  base_index[0] = 1;
1554  else
1555  base_index[0] = 2;
1556 
1557  if ((indices[2] == 0) && (indices[3] == 0))
1558  base_index[1] = 0;
1559  else if ((indices[2] == 1) && (indices[3] == 1))
1560  base_index[1] = 1;
1561  else
1562  base_index[1] = 2;
1563 
1564  return data[base_index[0]][base_index[1]];
1565  }
1566 
1567  case 3:
1568  // each entry of the tensor can be
1569  // thought of as an entry in a
1570  // matrix that maps the rolled-out
1571  // rank-2 tensors into rolled-out
1572  // rank-2 tensors. this is the
1573  // format in which we store rank-4
1574  // tensors. determine which
1575  // position the present entry is
1576  // stored in
1577  {
1578  unsigned int base_index[2] ;
1579  if ((indices[0] == 0) && (indices[1] == 0))
1580  base_index[0] = 0;
1581  else if ((indices[0] == 1) && (indices[1] == 1))
1582  base_index[0] = 1;
1583  else if ((indices[0] == 2) && (indices[1] == 2))
1584  base_index[0] = 2;
1585  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1586  ((indices[0] == 1) && (indices[1] == 0)))
1587  base_index[0] = 3;
1588  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1589  ((indices[0] == 2) && (indices[1] == 0)))
1590  base_index[0] = 4;
1591  else
1592  {
1593  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1594  ((indices[0] == 2) && (indices[1] == 1)),
1595  ExcInternalError());
1596  base_index[0] = 5;
1597  }
1598 
1599  if ((indices[2] == 0) && (indices[3] == 0))
1600  base_index[1] = 0;
1601  else if ((indices[2] == 1) && (indices[3] == 1))
1602  base_index[1] = 1;
1603  else if ((indices[2] == 2) && (indices[3] == 2))
1604  base_index[1] = 2;
1605  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1606  ((indices[2] == 1) && (indices[3] == 0)))
1607  base_index[1] = 3;
1608  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1609  ((indices[2] == 2) && (indices[3] == 0)))
1610  base_index[1] = 4;
1611  else
1612  {
1613  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1614  ((indices[2] == 2) && (indices[3] == 1)),
1615  ExcInternalError());
1616  base_index[1] = 5;
1617  }
1618 
1619  return data[base_index[0]][base_index[1]];
1620  }
1621 
1622  default:
1623  Assert (false, ExcNotImplemented());
1624  }
1625 
1626  static Number dummy;
1627  return dummy;
1628  }
1629 
1630 } // end of namespace internal
1631 
1632 
1633 
1634 template <int rank, int dim, typename Number>
1635 inline
1636 Number &
1638 {
1639  for (unsigned int r=0; r<rank; ++r)
1640  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1641  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1642 }
1643 
1644 
1645 
1646 template <int rank, int dim, typename Number>
1647 inline
1648 Number
1650 (const TableIndices<rank> &indices) const
1651 {
1652  for (unsigned int r=0; r<rank; ++r)
1653  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1654  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1655 }
1656 
1657 
1658 
1659 template <int rank, int dim, typename Number>
1660 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1661 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1662 {
1663  return
1664  internal::SymmetricTensorAccessors::
1665  Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
1666 }
1667 
1668 
1669 
1670 template <int rank, int dim, typename Number>
1671 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
1672 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
1673 {
1674  return
1675  internal::SymmetricTensorAccessors::
1676  Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
1677 }
1678 
1679 
1680 
1681 template <int rank, int dim, typename Number>
1682 inline
1683 Number
1685 {
1686  return operator()(indices);
1687 }
1688 
1689 
1690 
1691 template <int rank, int dim, typename Number>
1692 inline
1693 Number &
1695 {
1696  return operator()(indices);
1697 }
1698 
1699 
1700 
1701 
1702 namespace internal
1703 {
1704  namespace SymmetricTensor
1705  {
1706  template <int dim, typename Number>
1707  unsigned int
1708  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
1709  const unsigned int index)
1710  {
1711  return index;
1712  }
1713 
1714 
1715  template <int dim, typename Number>
1717  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
1718  const unsigned int index)
1719  {
1720  return
1723  }
1724 
1725  }
1726 }
1727 
1728 
1729 
1730 template <int rank, int dim, typename Number>
1731 inline
1732 Number
1733 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
1734 {
1735  AssertIndexRange (index, n_independent_components);
1736  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1737 }
1738 
1739 
1740 
1741 template <int rank, int dim, typename Number>
1742 inline
1743 Number &
1744 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
1745 {
1746  AssertIndexRange (index, n_independent_components);
1747  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1748 }
1749 
1750 
1751 
1752 namespace internal
1753 {
1754  template <int dim, typename Number>
1755  inline
1758  {
1759  switch (dim)
1760  {
1761  case 1:
1762  return numbers::NumberTraits<Number>::abs(data[0]);
1763 
1764  case 2:
1765  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1768 
1769  case 3:
1770  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1776 
1777  default:
1778  {
1779  typename numbers::NumberTraits<Number>::real_type return_value
1781 
1782  for (unsigned int d=0; d<dim; ++d)
1783  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
1784  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
1785  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
1786 
1787  return std::sqrt(return_value);
1788  }
1789  }
1790  }
1791 
1792 
1793 
1794  template <int dim, typename Number>
1795  inline
1798  {
1799  switch (dim)
1800  {
1801  case 1:
1802  return numbers::NumberTraits<Number>::abs (data[0][0]);
1803 
1804  default:
1805  {
1806  typename numbers::NumberTraits<Number>::real_type return_value
1808 
1809  const unsigned int n_independent_components = data.dimension;
1810 
1811  for (unsigned int i=0; i<dim; ++i)
1812  for (unsigned int j=0; j<dim; ++j)
1813  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
1814  for (unsigned int i=0; i<dim; ++i)
1815  for (unsigned int j=dim; j<n_independent_components; ++j)
1816  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1817  for (unsigned int i=dim; i<n_independent_components; ++i)
1818  for (unsigned int j=0; j<dim; ++j)
1819  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1820  for (unsigned int i=dim; i<n_independent_components; ++i)
1821  for (unsigned int j=dim; j<n_independent_components; ++j)
1822  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1823 
1824  return std::sqrt(return_value);
1825  }
1826  }
1827  }
1828 
1829 } // end of namespace internal
1830 
1831 
1832 
1833 template <int rank, int dim, typename Number>
1834 inline
1837 {
1838  return internal::compute_norm<dim,Number> (data);
1839 }
1840 
1841 
1842 
1843 namespace internal
1844 {
1845  namespace SymmetricTensor
1846  {
1847  namespace
1848  {
1849  // a function to do the unrolling from a set of indices to a
1850  // scalar index into the array in which we store the elements of
1851  // a symmetric tensor
1852  //
1853  // this function is for rank-2 tensors
1854  template <int dim>
1855  inline
1856  unsigned int
1857  component_to_unrolled_index
1858  (const TableIndices<2> &indices)
1859  {
1860  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
1861  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
1862 
1863  switch (dim)
1864  {
1865  case 1:
1866  {
1867  return 0;
1868  }
1869 
1870  case 2:
1871  {
1872  static const unsigned int table[2][2] = {{0, 2},
1873  {2, 1}
1874  };
1875  return table[indices[0]][indices[1]];
1876  }
1877 
1878  case 3:
1879  {
1880  static const unsigned int table[3][3] = {{0, 3, 4},
1881  {3, 1, 5},
1882  {4, 5, 2}
1883  };
1884  return table[indices[0]][indices[1]];
1885  }
1886 
1887  case 4:
1888  {
1889  static const unsigned int table[4][4] = {{0, 4, 5, 6},
1890  {4, 1, 7, 8},
1891  {5, 7, 2, 9},
1892  {6, 8, 9, 3}
1893  };
1894  return table[indices[0]][indices[1]];
1895  }
1896 
1897  default:
1898  // for the remainder, manually figure out the numbering
1899  {
1900  if (indices[0] == indices[1])
1901  return indices[0];
1902 
1903  TableIndices<2> sorted_indices (indices);
1904  sorted_indices.sort ();
1905 
1906  for (unsigned int d=0, c=0; d<dim; ++d)
1907  for (unsigned int e=d+1; e<dim; ++e, ++c)
1908  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1909  return dim+c;
1910 
1911  // should never get here:
1912  Assert(false, ExcInternalError());
1913  return 0;
1914  }
1915  }
1916  }
1917 
1918  // a function to do the unrolling from a set of indices to a
1919  // scalar index into the array in which we store the elements of
1920  // a symmetric tensor
1921  //
1922  // this function is for tensors of ranks not already handled
1923  // above
1924  template <int dim, int rank>
1925  inline
1926  unsigned int
1927  component_to_unrolled_index
1928  (const TableIndices<rank> &indices)
1929  {
1930  (void)indices;
1931  Assert (false, ExcNotImplemented());
1933  }
1934  }
1935  }
1936 }
1937 
1938 
1939 template <int rank, int dim, typename Number>
1940 inline
1941 unsigned int
1943 (const TableIndices<rank> &indices)
1944 {
1945  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
1946 }
1947 
1948 
1949 
1950 namespace internal
1951 {
1952  namespace SymmetricTensor
1953  {
1954  namespace
1955  {
1956  // a function to do the inverse of the unrolling from a set of
1957  // indices to a scalar index into the array in which we store
1958  // the elements of a symmetric tensor. in other words, it goes
1959  // from the scalar index into the array to a set of indices of
1960  // the tensor
1961  //
1962  // this function is for rank-2 tensors
1963  template <int dim>
1964  inline
1966  unrolled_to_component_indices
1967  (const unsigned int i,
1968  const int2type<2> &)
1969  {
1972  switch (dim)
1973  {
1974  case 1:
1975  {
1976  return TableIndices<2>(0,0);
1977  }
1978 
1979  case 2:
1980  {
1981  const TableIndices<2> table[3] =
1982  {
1983  TableIndices<2> (0,0),
1984  TableIndices<2> (1,1),
1985  TableIndices<2> (0,1)
1986  };
1987  return table[i];
1988  }
1989 
1990  case 3:
1991  {
1992  const TableIndices<2> table[6] =
1993  {
1994  TableIndices<2> (0,0),
1995  TableIndices<2> (1,1),
1996  TableIndices<2> (2,2),
1997  TableIndices<2> (0,1),
1998  TableIndices<2> (0,2),
1999  TableIndices<2> (1,2)
2000  };
2001  return table[i];
2002  }
2003 
2004  default:
2005  if (i<dim)
2006  return TableIndices<2> (i,i);
2007 
2008  for (unsigned int d=0, c=0; d<dim; ++d)
2009  for (unsigned int e=d+1; e<dim; ++e, ++c)
2010  if (c==i)
2011  return TableIndices<2>(d,e);
2012 
2013  // should never get here:
2014  Assert(false, ExcInternalError());
2015  return TableIndices<2>(0, 0);
2016  }
2017  }
2018 
2019  // a function to do the inverse of the unrolling from a set of
2020  // indices to a scalar index into the array in which we store
2021  // the elements of a symmetric tensor. in other words, it goes
2022  // from the scalar index into the array to a set of indices of
2023  // the tensor
2024  //
2025  // this function is for tensors of a rank not already handled
2026  // above
2027  template <int dim, int rank>
2028  inline
2030  unrolled_to_component_indices
2031  (const unsigned int i,
2032  const int2type<rank> &)
2033  {
2034  (void)i;
2037  Assert (false, ExcNotImplemented());
2038  return TableIndices<rank>();
2039  }
2040 
2041  }
2042  }
2043 }
2044 
2045 template <int rank, int dim, typename Number>
2046 inline
2049 (const unsigned int i)
2050 {
2051  return
2052  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2054 }
2055 
2056 
2057 
2058 template <int rank, int dim, typename Number>
2059 template <class Archive>
2060 inline
2061 void
2062 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2063 {
2064  ar &data;
2065 }
2066 
2067 
2068 #endif // DOXYGEN
2069 
2070 /* ----------------- Non-member functions operating on tensors. ------------ */
2071 
2072 
2079 template <int rank, int dim, typename Number, typename OtherNumber>
2080 inline
2083  const Tensor<rank, dim, OtherNumber> &right)
2084 {
2085  return Tensor<rank, dim, Number>(left) + right;
2086 }
2087 
2088 
2095 template <int rank, int dim, typename Number, typename OtherNumber>
2096 inline
2100 {
2101  return left + Tensor<rank, dim, OtherNumber>(right);
2102 }
2103 
2104 
2111 template <int rank, int dim, typename Number, typename OtherNumber>
2112 inline
2115  const Tensor<rank, dim, OtherNumber> &right)
2116 {
2117  return Tensor<rank, dim, Number>(left) - right;
2118 }
2119 
2120 
2127 template <int rank, int dim, typename Number, typename OtherNumber>
2128 inline
2132 {
2133  return left - Tensor<rank, dim, OtherNumber>(right);
2134 }
2135 
2136 
2137 
2151 template <int dim, typename Number>
2152 inline
2154 {
2155  switch (dim)
2156  {
2157  case 1:
2158  return t.data[0];
2159  case 2:
2160  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2161  case 3:
2162  {
2163  // in analogy to general tensors, but
2164  // there's something to be simplified for
2165  // the present case
2166  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2167  return ( tmp + tmp
2168  +t.data[0]*t.data[1]*t.data[2]
2169  -t.data[0]*t.data[5]*t.data[5]
2170  -t.data[1]*t.data[4]*t.data[4]
2171  -t.data[2]*t.data[3]*t.data[3]);
2172  }
2173  default:
2174  Assert (false, ExcNotImplemented());
2175  return 0;
2176  }
2177 }
2178 
2179 
2180 
2190 template <int dim, typename Number>
2191 inline
2193 {
2194  return determinant (t);
2195 }
2196 
2197 
2198 
2206 template <int dim, typename Number>
2208 {
2209  Number t = d.data[0];
2210  for (unsigned int i=1; i<dim; ++i)
2211  t += d.data[i];
2212  return t;
2213 }
2214 
2215 
2225 template <int dim, typename Number>
2226 inline
2228 {
2229  return trace (t);
2230 }
2231 
2232 
2245 template <typename Number>
2246 inline
2248 {
2249  return 0;
2250 }
2251 
2252 
2253 
2274 template <typename Number>
2275 inline
2277 {
2278  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2279 }
2280 
2281 
2282 
2292 template <typename Number>
2293 inline
2295 {
2296  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2297  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2298 }
2299 
2300 
2301 
2302 
2312 template <int rank, int dim, typename Number>
2313 inline
2316 {
2317  return t;
2318 }
2319 
2320 
2321 
2331 template <int dim, typename Number>
2332 inline
2335 {
2337 
2338  // subtract scaled trace from the diagonal
2339  const Number tr = trace(t) / dim;
2340  for (unsigned int i=0; i<dim; ++i)
2341  tmp.data[i] -= tr;
2342 
2343  return tmp;
2344 }
2345 
2346 
2347 
2355 template <int dim, typename Number>
2356 inline
2358 unit_symmetric_tensor ()
2359 {
2360  // create a default constructed matrix filled with
2361  // zeros, then set the diagonal elements to one
2363  switch (dim)
2364  {
2365  case 1:
2366  tmp.data[0] = 1;
2367  break;
2368  case 2:
2369  tmp.data[0] = tmp.data[1] = 1;
2370  break;
2371  case 3:
2372  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
2373  break;
2374  default:
2375  for (unsigned int d=0; d<dim; ++d)
2376  tmp.data[d] = 1;
2377  }
2378  return tmp;
2379 }
2380 
2381 
2382 
2391 template <int dim>
2392 inline
2394 unit_symmetric_tensor ()
2395 {
2396  return unit_symmetric_tensor<dim,double>();
2397 }
2398 
2399 
2400 
2415 template <int dim, typename Number>
2416 inline
2418 deviator_tensor ()
2419 {
2421 
2422  // fill the elements treating the diagonal
2423  for (unsigned int i=0; i<dim; ++i)
2424  for (unsigned int j=0; j<dim; ++j)
2425  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
2426 
2427  // then fill the ones that copy over the
2428  // non-diagonal elements. note that during
2429  // the double-contraction, we handle the
2430  // off-diagonal elements twice, so simply
2431  // copying requires a weight of 1/2
2432  for (unsigned int i=dim;
2433  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2434  ++i)
2435  tmp.data[i][i] = 0.5;
2436 
2437  return tmp;
2438 }
2439 
2440 
2441 
2456 template <int dim>
2457 inline
2459 deviator_tensor ()
2460 {
2461  return deviator_tensor<dim,double>();
2462 }
2463 
2464 
2465 
2488 template <int dim, typename Number>
2489 inline
2491 identity_tensor ()
2492 {
2494 
2495  // fill the elements treating the diagonal
2496  for (unsigned int i=0; i<dim; ++i)
2497  tmp.data[i][i] = 1;
2498 
2499  // then fill the ones that copy over the
2500  // non-diagonal elements. note that during
2501  // the double-contraction, we handle the
2502  // off-diagonal elements twice, so simply
2503  // copying requires a weight of 1/2
2504  for (unsigned int i=dim;
2505  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2506  ++i)
2507  tmp.data[i][i] = 0.5;
2508 
2509  return tmp;
2510 }
2511 
2512 
2513 
2535 template <int dim>
2536 inline
2538 identity_tensor ()
2539 {
2540  return identity_tensor<dim,double>();
2541 }
2542 
2543 
2544 
2555 template <int dim, typename Number>
2556 inline
2559 {
2560  // if desired, take over the
2561  // inversion of a 4x4 tensor
2562  // from the FullMatrix
2563  AssertThrow (false, ExcNotImplemented());
2564 
2566 }
2567 
2568 
2569 
2570 #ifndef DOXYGEN
2571 
2572 template <typename Number>
2573 inline
2575 invert (const SymmetricTensor<2,1,Number> &t)
2576 {
2578 
2579  tmp[0][0] = 1.0/t[0][0];
2580 
2581  return tmp;
2582 }
2583 
2584 
2585 
2586 template <typename Number>
2587 inline
2589 invert (const SymmetricTensor<2,2,Number> &t)
2590 {
2592 
2593  // Sympy result: ([
2594  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
2595  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
2596  const TableIndices<2> idx_00 (0,0);
2597  const TableIndices<2> idx_01 (0,1);
2598  const TableIndices<2> idx_11 (1,1);
2599  const Number inv_det_t
2600  = 1.0/(t[idx_00]*t[idx_11]
2601  - t[idx_01]*t[idx_01]);
2602  tmp[idx_00] = t[idx_11];
2603  tmp[idx_01] = -t[idx_01];
2604  tmp[idx_11] = t[idx_00];
2605  tmp *= inv_det_t;
2606 
2607  return tmp;
2608 }
2609 
2610 
2611 
2612 template <typename Number>
2613 inline
2615 invert (const SymmetricTensor<2,3,Number> &t)
2616 {
2618 
2619  // Sympy result: ([
2620  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2621  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2622  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
2623  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2624  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2625  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
2626  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2627  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
2628  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
2629  const TableIndices<2> idx_00 (0,0);
2630  const TableIndices<2> idx_01 (0,1);
2631  const TableIndices<2> idx_02 (0,2);
2632  const TableIndices<2> idx_11 (1,1);
2633  const TableIndices<2> idx_12 (1,2);
2634  const TableIndices<2> idx_22 (2,2);
2635  const Number inv_det_t
2636  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
2637  - t[idx_00]*t[idx_12]*t[idx_12]
2638  - t[idx_01]*t[idx_01]*t[idx_22]
2639  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
2640  - t[idx_02]*t[idx_02]*t[idx_11]);
2641  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
2642  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
2643  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
2644  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
2645  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
2646  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
2647  tmp *= inv_det_t;
2648 
2649  return tmp;
2650 }
2651 
2652 #endif /* DOXYGEN */
2653 
2654 
2655 
2669 template <int dim, typename Number>
2670 inline
2673 {
2675  switch (dim)
2676  {
2677  case 1:
2678  tmp.data[0][0] = 1./t.data[0][0];
2679  break;
2680  case 2:
2681 
2682  // inverting this tensor is a little more
2683  // complicated than necessary, since we
2684  // store the data of 't' as a 3x3 matrix
2685  // t.data, but the product between a rank-4
2686  // and a rank-2 tensor is really not the
2687  // product between this matrix and the
2688  // 3-vector of a rhs, but rather
2689  //
2690  // B.vec = t.data * mult * A.vec
2691  //
2692  // where mult is a 3x3 matrix with
2693  // entries [[1,0,0],[0,1,0],[0,0,2]] to
2694  // capture the fact that we need to add up
2695  // both the c_ij12*a_12 and the c_ij21*a_21
2696  // terms
2697  //
2698  // in addition, in this scheme, the
2699  // identity tensor has the matrix
2700  // representation mult^-1.
2701  //
2702  // the inverse of 't' therefore has the
2703  // matrix representation
2704  //
2705  // inv.data = mult^-1 * t.data^-1 * mult^-1
2706  //
2707  // in order to compute it, let's first
2708  // compute the inverse of t.data and put it
2709  // into tmp.data; at the end of the
2710  // function we then scale the last row and
2711  // column of the inverse by 1/2,
2712  // corresponding to the left and right
2713  // multiplication with mult^-1
2714  {
2715  const Number t4 = t.data[0][0]*t.data[1][1],
2716  t6 = t.data[0][0]*t.data[1][2],
2717  t8 = t.data[0][1]*t.data[1][0],
2718  t00 = t.data[0][2]*t.data[1][0],
2719  t01 = t.data[0][1]*t.data[2][0],
2720  t04 = t.data[0][2]*t.data[2][0],
2721  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
2722  t8*t.data[2][2]+t00*t.data[2][1]+
2723  t01*t.data[1][2]-t04*t.data[1][1]);
2724  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
2725  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
2726  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
2727  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
2728  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
2729  tmp.data[1][2] = -(t6-t00)*t07;
2730  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
2731  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
2732  tmp.data[2][2] = (t4-t8)*t07;
2733 
2734  // scale last row and column as mentioned
2735  // above
2736  tmp.data[2][0] /= 2;
2737  tmp.data[2][1] /= 2;
2738  tmp.data[0][2] /= 2;
2739  tmp.data[1][2] /= 2;
2740  tmp.data[2][2] /= 4;
2741  }
2742  break;
2743  default:
2744  Assert (false, ExcNotImplemented());
2745  }
2746  return tmp;
2747 }
2748 
2749 
2750 
2764 template <>
2766 invert (const SymmetricTensor<4,3,double> &t);
2767 // this function is implemented in the .cc file for double data types
2768 
2769 
2770 
2785 template <int dim, typename Number>
2786 inline
2790 {
2792 
2793  // fill only the elements really needed
2794  for (unsigned int i=0; i<dim; ++i)
2795  for (unsigned int j=i; j<dim; ++j)
2796  for (unsigned int k=0; k<dim; ++k)
2797  for (unsigned int l=k; l<dim; ++l)
2798  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
2799 
2800  return tmp;
2801 }
2802 
2803 
2804 
2813 template <int dim,typename Number>
2814 inline
2817 {
2818  Number array[(dim*dim+dim)/2];
2819  for (unsigned int d=0; d<dim; ++d)
2820  array[d] = t[d][d];
2821  for (unsigned int d=0, c=0; d<dim; ++d)
2822  for (unsigned int e=d+1; e<dim; ++e, ++c)
2823  array[dim+c] = (t[d][e]+t[e][d])*0.5;
2824  return SymmetricTensor<2,dim,Number>(array);
2825 }
2826 
2827 
2828 
2836 template <int rank, int dim, typename Number>
2837 inline
2840  const Number factor)
2841 {
2843  tt *= factor;
2844  return tt;
2845 }
2846 
2847 
2848 
2856 template <int rank, int dim, typename Number>
2857 inline
2859 operator * (const Number factor,
2861 {
2862  // simply forward to the other operator
2863  return t*factor;
2864 }
2865 
2866 
2867 
2893 template <int rank, int dim, typename Number, typename OtherNumber>
2894 inline
2897  const OtherNumber factor)
2898 {
2899  // form the product. we have to convert the two factors into the final
2900  // type via explicit casts because, for awkward reasons, the C++
2901  // standard committee saw it fit to not define an
2902  // operator*(float,std::complex<double>)
2903  // (as well as with switched arguments and double<->float).
2904  typedef typename ProductType<Number,OtherNumber>::type product_type;
2906  // we used to shorten the following by 'tt *= product_type(factor);'
2907  // which requires that a converting constructor
2908  // 'product_type::product_type(const OtherNumber) is defined.
2909  // however, a user-defined constructor is not allowed for aggregates,
2910  // e.g. VectorizedArray. therefore, we work around this issue using a
2911  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
2912  // which we assume to be defined.
2913  product_type new_factor;
2914  new_factor = factor;
2915  tt *= new_factor;
2916  return tt;
2917 }
2918 
2919 
2920 
2929 template <int rank, int dim, typename Number, typename OtherNumber>
2930 inline
2932 operator * (const Number factor,
2934 {
2935  // simply forward to the other operator with switched arguments
2936  return (t*factor);
2937 }
2938 
2939 
2940 
2946 template <int rank, int dim, typename Number>
2947 inline
2950  const Number factor)
2951 {
2953  tt /= factor;
2954  return tt;
2955 }
2956 
2957 
2958 
2965 template <int rank, int dim>
2966 inline
2968 operator * (const SymmetricTensor<rank,dim> &t,
2969  const double factor)
2970 {
2972  tt *= factor;
2973  return tt;
2974 }
2975 
2976 
2977 
2984 template <int rank, int dim>
2985 inline
2987 operator * (const double factor,
2988  const SymmetricTensor<rank,dim> &t)
2989 {
2991  tt *= factor;
2992  return tt;
2993 }
2994 
2995 
2996 
3002 template <int rank, int dim>
3003 inline
3005 operator / (const SymmetricTensor<rank,dim> &t,
3006  const double factor)
3007 {
3009  tt /= factor;
3010  return tt;
3011 }
3012 
3022 template <int dim, typename Number>
3023 inline
3024 Number
3027 {
3028  return (t1*t2);
3029 }
3030 
3031 
3041 template <int dim, typename Number>
3042 inline
3043 Number
3045  const Tensor<2,dim,Number> &t2)
3046 {
3047  Number s = 0;
3048  for (unsigned int i=0; i<dim; ++i)
3049  for (unsigned int j=0; j<dim; ++j)
3050  s += t1[i][j] * t2[i][j];
3051  return s;
3052 }
3053 
3054 
3064 template <int dim, typename Number>
3065 inline
3066 Number
3069 {
3070  return scalar_product(t2, t1);
3071 }
3072 
3073 
3089 template <typename Number>
3090 inline
3091 void
3093  const SymmetricTensor<4,1,Number> &t,
3094  const SymmetricTensor<2,1,Number> &s)
3095 {
3096  tmp[0][0] = t[0][0][0][0] * s[0][0];
3097 }
3098 
3099 
3100 
3116 template <typename Number>
3117 inline
3118 void
3120  const SymmetricTensor<2,1,Number> &s,
3121  const SymmetricTensor<4,1,Number> &t)
3122 {
3123  tmp[0][0] = t[0][0][0][0] * s[0][0];
3124 }
3125 
3126 
3127 
3142 template <typename Number>
3143 inline
3144 void
3146  const SymmetricTensor<4,2,Number> &t,
3147  const SymmetricTensor<2,2,Number> &s)
3148 {
3149  const unsigned int dim = 2;
3150 
3151  for (unsigned int i=0; i<dim; ++i)
3152  for (unsigned int j=i; j<dim; ++j)
3153  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3154  t[i][j][1][1] * s[1][1] +
3155  2 * t[i][j][0][1] * s[0][1];
3156 }
3157 
3158 
3159 
3175 template <typename Number>
3176 inline
3177 void
3179  const SymmetricTensor<2,2,Number> &s,
3180  const SymmetricTensor<4,2,Number> &t)
3181 {
3182  const unsigned int dim = 2;
3183 
3184  for (unsigned int i=0; i<dim; ++i)
3185  for (unsigned int j=i; j<dim; ++j)
3186  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3187  s[1][1] * t[1][1][i][j] +
3188  2 * s[0][1] * t[0][1][i][j];
3189 }
3190 
3191 
3192 
3208 template <typename Number>
3209 inline
3210 void
3212  const SymmetricTensor<4,3,Number> &t,
3213  const SymmetricTensor<2,3,Number> &s)
3214 {
3215  const unsigned int dim = 3;
3216 
3217  for (unsigned int i=0; i<dim; ++i)
3218  for (unsigned int j=i; j<dim; ++j)
3219  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3220  t[i][j][1][1] * s[1][1] +
3221  t[i][j][2][2] * s[2][2] +
3222  2 * t[i][j][0][1] * s[0][1] +
3223  2 * t[i][j][0][2] * s[0][2] +
3224  2 * t[i][j][1][2] * s[1][2];
3225 }
3226 
3227 
3228 
3244 template <typename Number>
3245 inline
3246 void
3248  const SymmetricTensor<2,3,Number> &s,
3249  const SymmetricTensor<4,3,Number> &t)
3250 {
3251  const unsigned int dim = 3;
3252 
3253  for (unsigned int i=0; i<dim; ++i)
3254  for (unsigned int j=i; j<dim; ++j)
3255  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3256  s[1][1] * t[1][1][i][j] +
3257  s[2][2] * t[2][2][i][j] +
3258  2 * s[0][1] * t[0][1][i][j] +
3259  2 * s[0][2] * t[0][2][i][j] +
3260  2 * s[1][2] * t[1][2][i][j];
3261 }
3262 
3263 
3264 
3272 template <int dim, typename Number>
3274 operator * (const SymmetricTensor<2,dim,Number> &src1,
3275  const Tensor<1,dim,Number> &src2)
3276 {
3277  Tensor<1,dim,Number> dest;
3278  for (unsigned int i=0; i<dim; ++i)
3279  for (unsigned int j=0; j<dim; ++j)
3280  dest[i] += src1[i][j] * src2[j];
3281  return dest;
3282 }
3283 
3284 
3292 template <int dim, typename Number>
3294 operator * (const Tensor<1,dim,Number> &src1,
3295  const SymmetricTensor<2,dim,Number> &src2)
3296 {
3297  // this is easy for symmetric tensors:
3298  return src2 * src1;
3299 }
3300 
3301 
3311 template <int dim, typename Number>
3312 inline
3313 std::ostream &operator << (std::ostream &out,
3315 {
3316  //make out lives a bit simpler by outputing
3317  //the tensor through the operator for the
3318  //general Tensor class
3320 
3321  for (unsigned int i=0; i<dim; ++i)
3322  for (unsigned int j=0; j<dim; ++j)
3323  tt[i][j] = t[i][j];
3324 
3325  return out << tt;
3326 }
3327 
3328 
3329 
3339 template <int dim, typename Number>
3340 inline
3341 std::ostream &operator << (std::ostream &out,
3343 {
3344  //make out lives a bit simpler by outputing
3345  //the tensor through the operator for the
3346  //general Tensor class
3348 
3349  for (unsigned int i=0; i<dim; ++i)
3350  for (unsigned int j=0; j<dim; ++j)
3351  for (unsigned int k=0; k<dim; ++k)
3352  for (unsigned int l=0; l<dim; ++l)
3353  tt[i][j][k][l] = t[i][j][k][l];
3354 
3355  return out << tt;
3356 }
3357 
3358 
3359 DEAL_II_NAMESPACE_CLOSE
3360 
3361 #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:365
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)
SymmetricTensor()=default
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()