Reference documentation for deal.II version Git ab1cc5b 2017-03-24 06:23:49 -0600
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 #include <deal.II/base/vectorization.h>
25 
26 DEAL_II_NAMESPACE_OPEN
27 
28 template <int rank, int dim, typename Number=double> class SymmetricTensor;
29 
30 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
32 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
34 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
36 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
38 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
40 template <int dim2, typename Number> Number
42 
43 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
45 template <int dim, typename Number> Number
47 
48 
49 
50 namespace internal
51 {
56  namespace SymmetricTensorAccessors
57  {
64  inline
65  TableIndices<2> merge (const TableIndices<2> &previous_indices,
66  const unsigned int new_index,
67  const unsigned int position)
68  {
69  Assert (position < 2, ExcIndexRange (position, 0, 2));
70 
71  if (position == 0)
72  return TableIndices<2>(new_index);
73  else
74  return TableIndices<2>(previous_indices[0], new_index);
75  }
76 
77 
78 
85  inline
86  TableIndices<4> merge (const TableIndices<4> &previous_indices,
87  const unsigned int new_index,
88  const unsigned int position)
89  {
90  Assert (position < 4, ExcIndexRange (position, 0, 4));
91 
92  switch (position)
93  {
94  case 0:
95  return TableIndices<4>(new_index);
96  case 1:
97  return TableIndices<4>(previous_indices[0],
98  new_index);
99  case 2:
100  return TableIndices<4>(previous_indices[0],
101  previous_indices[1],
102  new_index);
103  case 3:
104  return TableIndices<4>(previous_indices[0],
105  previous_indices[1],
106  previous_indices[2],
107  new_index);
108  }
109  Assert (false, ExcInternalError());
110  return TableIndices<4>();
111  }
112 
113 
122  template <int rank1, int rank2, int dim, typename Number>
124  {
125  typedef ::SymmetricTensor<rank1+rank2-4,dim,Number> type;
126  };
127 
128 
137  template <int dim, typename Number>
138  struct double_contraction_result<2,2,dim,Number>
139  {
140  typedef Number type;
141  };
142 
143 
144 
157  template <int rank, int dim, typename Number>
158  struct StorageType;
159 
163  template <int dim, typename Number>
164  struct StorageType<2,dim,Number>
165  {
170  static const unsigned int
171  n_independent_components = (dim*dim + dim)/2;
172 
177  };
178 
179 
180 
184  template <int dim, typename Number>
185  struct StorageType<4,dim,Number>
186  {
192  static const unsigned int
193  n_rank2_components = (dim*dim + dim)/2;
194 
198  static const unsigned int
199  n_independent_components = (n_rank2_components *
201 
209  };
210 
211 
212 
217  template <int rank, int dim, bool constness, typename Number>
219 
226  template <int rank, int dim, typename Number>
227  struct AccessorTypes<rank,dim,true,Number>
228  {
229  typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
230 
231  typedef Number reference;
232  };
233 
240  template <int rank, int dim, typename Number>
241  struct AccessorTypes<rank,dim,false,Number>
242  {
243  typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
244 
245  typedef Number &reference;
246  };
247 
248 
283  template <int rank, int dim, bool constness, int P, typename Number>
284  class Accessor
285  {
286  public:
290  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
291  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
292 
293  private:
312  Accessor (tensor_type &tensor,
313  const TableIndices<rank> &previous_indices);
314 
318  Accessor (const Accessor &a);
319 
320  public:
321 
325  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
326 
330  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;
331 
332  private:
336  tensor_type &tensor;
337  const TableIndices<rank> previous_indices;
338 
339  // declare some other classes
340  // as friends. make sure to
341  // work around bugs in some
342  // compilers
343  template <int,int,typename> friend class ::SymmetricTensor;
344  template <int,int,bool,int,typename>
345  friend class Accessor;
346 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
347  friend class ::SymmetricTensor<rank,dim,Number>;
348  friend class Accessor<rank,dim,constness,P+1,Number>;
349 # endif
350  };
351 
352 
353 
363  template <int rank, int dim, bool constness, typename Number>
364  class Accessor<rank,dim,constness,1,Number>
365  {
366  public:
370  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
371  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
372 
373  private:
395  Accessor (tensor_type &tensor,
396  const TableIndices<rank> &previous_indices);
397 
401  Accessor ();
402 
406  Accessor (const Accessor &a);
407 
408  public:
409 
413  reference operator [] (const unsigned int);
414 
418  reference operator [] (const unsigned int) const;
419 
420  private:
424  tensor_type &tensor;
425  const TableIndices<rank> previous_indices;
426 
427  // declare some other classes
428  // as friends. make sure to
429  // work around bugs in some
430  // compilers
431  template <int,int,typename> friend class ::SymmetricTensor;
432  template <int,int,bool,int,typename>
433  friend class SymmetricTensorAccessors::Accessor;
434 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
435  friend class ::SymmetricTensor<rank,dim,Number>;
436  friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
437 # endif
438  };
439  }
440 }
441 
442 
443 
507 template <int rank, int dim, typename Number>
508 class SymmetricTensor
509 {
510 public:
519  static const unsigned int dimension = dim;
520 
526  static const unsigned int n_independent_components
529 
533  SymmetricTensor ();
534 
546 
562  SymmetricTensor (const Number (&array) [n_independent_components]);
563 
569  template <typename OtherNumber>
570  explicit
572 
579  SymmetricTensor &operator = (const Number d);
580 
585  operator Tensor<rank,dim,Number> () const;
586 
590  bool operator == (const SymmetricTensor &) const;
591 
595  bool operator != (const SymmetricTensor &) const;
596 
601 
606 
611  SymmetricTensor &operator *= (const Number factor);
612 
616  SymmetricTensor &operator /= (const Number factor);
617 
623 
629 
634 
661 
668 
672  Number &operator() (const TableIndices<rank> &indices);
673 
677  Number operator() (const TableIndices<rank> &indices) const;
678 
683  internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
684  operator [] (const unsigned int row) const;
685 
690  internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
691  operator [] (const unsigned int row);
692 
698  Number
699  operator [] (const TableIndices<rank> &indices) const;
700 
706  Number &
707  operator [] (const TableIndices<rank> &indices);
708 
714  Number
715  access_raw_entry (const unsigned int unrolled_index) const;
716 
722  Number &
723  access_raw_entry (const unsigned int unrolled_index);
724 
735  norm () const;
736 
744  static
745  unsigned int
747 
753  static
755  unrolled_to_component_indices (const unsigned int i);
756 
769  void clear ();
770 
775  static std::size_t memory_consumption ();
776 
781  template <class Archive>
782  void serialize(Archive &ar, const unsigned int version);
783 
784 private:
788  typedef
791 
795  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
796 
800  base_tensor_type data;
801 
805  template <int, int, typename> friend class SymmetricTensor;
806 
810  template <int dim2, typename Number2>
811  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
812 
813  template <int dim2, typename Number2>
814  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
815 
816  template <int dim2, typename Number2>
818  deviator (const SymmetricTensor<2,dim2,Number2> &t);
819 
820  template <int dim2, typename Number2>
822 
823  template <int dim2, typename Number2>
825 
826  template <int dim2, typename Number2>
828 
829  template <int dim2, typename Number2>
831 
832  template <int dim2, typename Number2>
834 };
835 
836 
837 
838 // ------------------------- inline functions ------------------------
839 
840 #ifndef DOXYGEN
841 
842 namespace internal
843 {
844  namespace SymmetricTensorAccessors
845  {
846  template <int rank, int dim, bool constness, int P, typename Number>
847  Accessor<rank,dim,constness,P,Number>::
848  Accessor (tensor_type &tensor,
849  const TableIndices<rank> &previous_indices)
850  :
851  tensor (tensor),
852  previous_indices (previous_indices)
853  {}
854 
855 
856  template <int rank, int dim, bool constness, int P, typename Number>
857  Accessor<rank,dim,constness,P,Number>::
858  Accessor (const Accessor &a)
859  :
860  tensor (a.tensor),
861  previous_indices (a.previous_indices)
862  {}
863 
864 
865 
866  template <int rank, int dim, bool constness, int P, typename Number>
867  Accessor<rank,dim,constness,P-1,Number>
868  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
869  {
870  return Accessor<rank,dim,constness,P-1,Number> (tensor,
871  merge (previous_indices, i, rank-P));
872  }
873 
874 
875 
876  template <int rank, int dim, bool constness, int P, typename Number>
877  Accessor<rank,dim,constness,P-1,Number>
878  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i) const
879  {
880  return Accessor<rank,dim,constness,P-1,Number> (tensor,
881  merge (previous_indices, i, rank-P));
882  }
883 
884 
885 
886  template <int rank, int dim, bool constness, typename Number>
887  Accessor<rank,dim,constness,1,Number>::
888  Accessor (tensor_type &tensor,
889  const TableIndices<rank> &previous_indices)
890  :
891  tensor (tensor),
892  previous_indices (previous_indices)
893  {}
894 
895 
896 
897  template <int rank, int dim, bool constness, typename Number>
898  Accessor<rank,dim,constness,1,Number>::
899  Accessor (const Accessor &a)
900  :
901  tensor (a.tensor),
902  previous_indices (a.previous_indices)
903  {}
904 
905 
906 
907  template <int rank, int dim, bool constness, typename Number>
908  typename Accessor<rank,dim,constness,1,Number>::reference
909  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
910  {
911  return tensor(merge (previous_indices, i, rank-1));
912  }
913 
914 
915  template <int rank, int dim, bool constness, typename Number>
916  typename Accessor<rank,dim,constness,1,Number>::reference
917  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i) const
918  {
919  return tensor(merge (previous_indices, i, rank-1));
920  }
921  }
922 }
923 
924 
925 
926 template <int rank, int dim, typename Number>
927 inline
929 {}
930 
931 
932 
933 template <int rank, int dim, typename Number>
934 inline
936 {
937  Assert (rank == 2, ExcNotImplemented());
938  switch (dim)
939  {
940  case 2:
941  Assert (t[0][1] == t[1][0], ExcInternalError());
942 
943  data[0] = t[0][0];
944  data[1] = t[1][1];
945  data[2] = t[0][1];
946 
947  break;
948  case 3:
949  Assert (t[0][1] == t[1][0], ExcInternalError());
950  Assert (t[0][2] == t[2][0], ExcInternalError());
951  Assert (t[1][2] == t[2][1], ExcInternalError());
952 
953  data[0] = t[0][0];
954  data[1] = t[1][1];
955  data[2] = t[2][2];
956  data[3] = t[0][1];
957  data[4] = t[0][2];
958  data[5] = t[1][2];
959 
960  break;
961  default:
962  for (unsigned int d=0; d<dim; ++d)
963  for (unsigned int e=0; e<d; ++e)
964  Assert(t[d][e] == t[e][d], ExcInternalError());
965 
966  for (unsigned int d=0; d<dim; ++d)
967  data[d] = t[d][d];
968 
969  for (unsigned int d=0, c=0; d<dim; ++d)
970  for (unsigned int e=d+1; e<dim; ++e, ++c)
971  data[dim+c] = t[d][e];
972  }
973 }
974 
975 
976 
977 template <int rank, int dim, typename Number>
978 template <typename OtherNumber>
979 inline
982 {
983  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
984  data[i] = initializer.data[i];
985 }
986 
987 
988 
989 
990 template <int rank, int dim, typename Number>
991 inline
992 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
993  :
994  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
995 {
996  // ensure that the reinterpret_cast above actually works
997  Assert (sizeof(typename base_tensor_type::array_type)
998  == sizeof(array),
999  ExcInternalError());
1000 }
1001 
1002 
1003 
1004 template <int rank, int dim, typename Number>
1005 inline
1008 {
1009  Assert (d==Number(), ExcMessage ("Only assignment with zero is allowed"));
1010  (void) d;
1011 
1012  data = 0;
1013 
1014  return *this;
1015 }
1016 
1017 
1018 namespace internal
1019 {
1020  namespace SymmetricTensor
1021  {
1022  template <int dim, typename Number>
1024  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1025  {
1026  Number t[dim][dim];
1027 
1028  // diagonal entries are stored first
1029  for (unsigned int d=0; d<dim; ++d)
1030  t[d][d] = s.access_raw_entry(d);
1031 
1032  // off-diagonal entries come next, row by row
1033  for (unsigned int d=0, c=0; d<dim; ++d)
1034  for (unsigned int e=d+1; e<dim; ++e, ++c)
1035  {
1036  t[d][e] = s.access_raw_entry(dim+c);
1037  t[e][d] = s.access_raw_entry(dim+c);
1038  }
1039  return ::Tensor<2,dim,Number>(t);
1040  }
1041 
1042 
1043  template <int dim, typename Number>
1045  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1046  {
1047  // utilize the symmetry properties of SymmetricTensor<4,dim>
1048  // discussed in the class documentation to avoid accessing all
1049  // independent elements of the input tensor more than once
1051 
1052  for (unsigned int i=0; i<dim; ++i)
1053  for (unsigned int j=i; j<dim; ++j)
1054  for (unsigned int k=0; k<dim; ++k)
1055  for (unsigned int l=k; l<dim; ++l)
1056  t[TableIndices<4>(i,j,k,l)]
1057  = t[TableIndices<4>(i,j,l,k)]
1058  = t[TableIndices<4>(j,i,k,l)]
1059  = t[TableIndices<4>(j,i,l,k)]
1060  = st[TableIndices<4>(i,j,k,l)];
1061 
1062  return t;
1063  }
1064  }
1065 }
1066 
1067 
1068 
1069 template <int rank, int dim, typename Number>
1070 inline
1072 operator Tensor<rank,dim,Number> () const
1073 {
1074  return internal::SymmetricTensor::convert_to_tensor (*this);
1075 }
1076 
1077 
1078 
1079 template <int rank, int dim, typename Number>
1080 inline
1081 bool
1083 (const SymmetricTensor<rank,dim,Number> &t) const
1084 {
1085  return data == t.data;
1086 }
1087 
1088 
1089 
1090 template <int rank, int dim, typename Number>
1091 inline
1092 bool
1093 SymmetricTensor<rank,dim,Number>::operator !=
1094 (const SymmetricTensor<rank,dim,Number> &t) const
1095 {
1096  return data != t.data;
1097 }
1098 
1099 
1100 
1101 template <int rank, int dim, typename Number>
1102 inline
1104 SymmetricTensor<rank,dim,Number>::operator +=
1106 {
1107  data += t.data;
1108  return *this;
1109 }
1110 
1111 
1112 
1113 template <int rank, int dim, typename Number>
1114 inline
1116 SymmetricTensor<rank,dim,Number>::operator -=
1118 {
1119  data -= t.data;
1120  return *this;
1121 }
1122 
1123 
1124 
1125 template <int rank, int dim, typename Number>
1126 inline
1129 {
1130  data *= d;
1131  return *this;
1132 }
1133 
1134 
1135 
1136 template <int rank, int dim, typename Number>
1137 inline
1140 {
1141  data /= d;
1142  return *this;
1143 }
1144 
1145 
1146 
1147 template <int rank, int dim, typename Number>
1148 inline
1151 {
1152  SymmetricTensor tmp = *this;
1153  tmp.data += t.data;
1154  return tmp;
1155 }
1156 
1157 
1158 
1159 template <int rank, int dim, typename Number>
1160 inline
1163 {
1164  SymmetricTensor tmp = *this;
1165  tmp.data -= t.data;
1166  return tmp;
1167 }
1168 
1169 
1170 
1171 template <int rank, int dim, typename Number>
1172 inline
1175 {
1176  SymmetricTensor tmp = *this;
1177  tmp.data = -tmp.data;
1178  return tmp;
1179 }
1180 
1181 
1182 
1183 template <int rank, int dim, typename Number>
1184 inline
1185 void
1187 {
1188  data.clear ();
1189 }
1190 
1191 
1192 
1193 template <int rank, int dim, typename Number>
1194 inline
1195 std::size_t
1197 {
1198  // all memory consists of statically allocated memory of the current
1199  // object, no pointers
1200  return sizeof(SymmetricTensor<rank,dim,Number>);
1201 }
1202 
1203 
1204 
1205 namespace internal
1206 {
1207 
1208  template <int dim, typename Number>
1209  inline
1210  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
1211  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1213  {
1214  switch (dim)
1215  {
1216  case 1:
1217  return data[0] * sdata[0];
1218  default:
1219  // Start with the non-diagonal part to avoid some multiplications by
1220  // 2.
1221  Number sum = data[dim] * sdata[dim];
1222  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1223  sum += data[d] * sdata[d];
1224  sum += sum; // sum = sum * 2.;
1225 
1226  // Now add the contributions from the diagonal
1227  for (unsigned int d=0; d<dim; ++d)
1228  sum += data[d] * sdata[d];
1229  return sum;
1230  }
1231  }
1232 
1233 
1234 
1235  template <int dim, typename Number>
1236  inline
1237  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
1238  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1240  {
1241  const unsigned int data_dim =
1243  Number tmp [data_dim];
1244  for (unsigned int i=0; i<data_dim; ++i)
1245  tmp[i] = perform_double_contraction<dim,Number>(data[i], sdata);
1246  return ::SymmetricTensor<2,dim,Number>(tmp);
1247  }
1248 
1249 
1250 
1251  template <int dim, typename Number>
1252  inline
1254  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1256  {
1258  for (unsigned int i=0; i<tmp.dimension; ++i)
1259  {
1260  // Start with the non-diagonal part
1261  Number sum = data[dim] * sdata[dim][i];
1262  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1263  sum += data[d] * sdata[d][i];
1264  sum += sum; // sum = sum * 2.;
1265 
1266  // Now add the contributions from the diagonal
1267  for (unsigned int d=0; d<dim; ++d)
1268  sum += data[d] * sdata[d][i];
1269  tmp[i] = sum;
1270  }
1271  return tmp;
1272  }
1273 
1274 
1275 
1276  template <int dim, typename Number>
1277  inline
1279  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1281  {
1282  const unsigned int data_dim =
1285  for (unsigned int i=0; i<data_dim; ++i)
1286  for (unsigned int j=0; j<data_dim; ++j)
1287  {
1288  // Start with the non-diagonal part
1289  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1290  tmp[i][j] += data[i][d] * sdata[d][j];
1291  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1292 
1293  // Now add the contributions from the diagonal
1294  for (unsigned int d=0; d<dim; ++d)
1295  tmp[i][j] += data[i][d] * sdata[d][j];
1296  }
1297  return tmp;
1298  }
1299 
1300 } // end of namespace internal
1301 
1302 
1303 
1304 template <int rank, int dim, typename Number>
1305 inline
1308 {
1309  // need to have two different function calls
1310  // because a scalar and rank-2 tensor are not
1311  // the same data type (see internal function
1312  // above)
1313  return internal::perform_double_contraction<dim,Number> (data, s.data);
1314 }
1315 
1316 
1317 
1318 template <int rank, int dim, typename Number>
1319 inline
1322 {
1325  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
1326  return tmp;
1327 }
1328 
1329 
1330 
1331 // internal namespace to switch between the
1332 // access of different tensors. There used to
1333 // be explicit instantiations before for
1334 // different ranks and dimensions, but since
1335 // we now allow for templates on the data
1336 // type, and since we cannot partially
1337 // specialize the implementation, this got
1338 // into a separate namespace
1339 namespace internal
1340 {
1341  template <int dim, typename Number>
1342  inline
1343  Number &
1344  symmetric_tensor_access (const TableIndices<2> &indices,
1346  {
1347  // 1d is very simple and done first
1348  if (dim == 1)
1349  return data[0];
1350 
1351  // first treat the main diagonal elements, which are stored consecutively
1352  // at the beginning
1353  if (indices[0] == indices[1])
1354  return data[indices[0]];
1355 
1356  // the rest is messier and requires a few switches.
1357  switch (dim)
1358  {
1359  case 2:
1360  // at least for the 2x2 case it is reasonably simple
1361  Assert (((indices[0]==1) && (indices[1]==0)) ||
1362  ((indices[0]==0) && (indices[1]==1)),
1363  ExcInternalError());
1364  return data[2];
1365 
1366  default:
1367  // to do the rest, sort our indices before comparing
1368  {
1369  TableIndices<2> sorted_indices (indices);
1370  sorted_indices.sort ();
1371 
1372  for (unsigned int d=0, c=0; d<dim; ++d)
1373  for (unsigned int e=d+1; e<dim; ++e, ++c)
1374  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1375  return data[dim+c];
1376  Assert (false, ExcInternalError());
1377  }
1378  }
1379 
1380  static Number dummy_but_referenceable = Number();
1381  return dummy_but_referenceable;
1382  }
1383 
1384 
1385 
1386  template <int dim, typename Number>
1387  inline
1388  Number
1389  symmetric_tensor_access (const TableIndices<2> &indices,
1391  {
1392  // 1d is very simple and done first
1393  if (dim == 1)
1394  return data[0];
1395 
1396  // first treat the main diagonal elements, which are stored consecutively
1397  // at the beginning
1398  if (indices[0] == indices[1])
1399  return data[indices[0]];
1400 
1401  // the rest is messier and requires a few switches.
1402  switch (dim)
1403  {
1404  case 2:
1405  // at least for the 2x2 case it is reasonably simple
1406  Assert (((indices[0]==1) && (indices[1]==0)) ||
1407  ((indices[0]==0) && (indices[1]==1)),
1408  ExcInternalError());
1409  return data[2];
1410 
1411  default:
1412  // to do the rest, sort our indices before comparing
1413  {
1414  TableIndices<2> sorted_indices (indices);
1415  sorted_indices.sort ();
1416 
1417  for (unsigned int d=0, c=0; d<dim; ++d)
1418  for (unsigned int e=d+1; e<dim; ++e, ++c)
1419  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1420  return data[dim+c];
1421  Assert (false, ExcInternalError());
1422  }
1423  }
1424 
1425  static Number dummy_but_referenceable = Number();
1426  return dummy_but_referenceable;
1427  }
1428 
1429 
1430 
1431  template <int dim, typename Number>
1432  inline
1433  Number &
1434  symmetric_tensor_access (const TableIndices<4> &indices,
1436  {
1437  switch (dim)
1438  {
1439  case 1:
1440  return data[0][0];
1441 
1442  case 2:
1443  // each entry of the tensor can be
1444  // thought of as an entry in a
1445  // matrix that maps the rolled-out
1446  // rank-2 tensors into rolled-out
1447  // rank-2 tensors. this is the
1448  // format in which we store rank-4
1449  // tensors. determine which
1450  // position the present entry is
1451  // stored in
1452  {
1453  unsigned int base_index[2] ;
1454  if ((indices[0] == 0) && (indices[1] == 0))
1455  base_index[0] = 0;
1456  else if ((indices[0] == 1) && (indices[1] == 1))
1457  base_index[0] = 1;
1458  else
1459  base_index[0] = 2;
1460 
1461  if ((indices[2] == 0) && (indices[3] == 0))
1462  base_index[1] = 0;
1463  else if ((indices[2] == 1) && (indices[3] == 1))
1464  base_index[1] = 1;
1465  else
1466  base_index[1] = 2;
1467 
1468  return data[base_index[0]][base_index[1]];
1469  }
1470 
1471  case 3:
1472  // each entry of the tensor can be
1473  // thought of as an entry in a
1474  // matrix that maps the rolled-out
1475  // rank-2 tensors into rolled-out
1476  // rank-2 tensors. this is the
1477  // format in which we store rank-4
1478  // tensors. determine which
1479  // position the present entry is
1480  // stored in
1481  {
1482  unsigned int base_index[2] ;
1483  if ((indices[0] == 0) && (indices[1] == 0))
1484  base_index[0] = 0;
1485  else if ((indices[0] == 1) && (indices[1] == 1))
1486  base_index[0] = 1;
1487  else if ((indices[0] == 2) && (indices[1] == 2))
1488  base_index[0] = 2;
1489  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1490  ((indices[0] == 1) && (indices[1] == 0)))
1491  base_index[0] = 3;
1492  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1493  ((indices[0] == 2) && (indices[1] == 0)))
1494  base_index[0] = 4;
1495  else
1496  {
1497  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1498  ((indices[0] == 2) && (indices[1] == 1)),
1499  ExcInternalError());
1500  base_index[0] = 5;
1501  }
1502 
1503  if ((indices[2] == 0) && (indices[3] == 0))
1504  base_index[1] = 0;
1505  else if ((indices[2] == 1) && (indices[3] == 1))
1506  base_index[1] = 1;
1507  else if ((indices[2] == 2) && (indices[3] == 2))
1508  base_index[1] = 2;
1509  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1510  ((indices[2] == 1) && (indices[3] == 0)))
1511  base_index[1] = 3;
1512  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1513  ((indices[2] == 2) && (indices[3] == 0)))
1514  base_index[1] = 4;
1515  else
1516  {
1517  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1518  ((indices[2] == 2) && (indices[3] == 1)),
1519  ExcInternalError());
1520  base_index[1] = 5;
1521  }
1522 
1523  return data[base_index[0]][base_index[1]];
1524  }
1525 
1526  default:
1527  Assert (false, ExcNotImplemented());
1528  }
1529 
1530  static Number dummy;
1531  return dummy;
1532  }
1533 
1534 
1535  template <int dim, typename Number>
1536  inline
1537  Number
1538  symmetric_tensor_access (const TableIndices<4> &indices,
1540  {
1541  switch (dim)
1542  {
1543  case 1:
1544  return data[0][0];
1545 
1546  case 2:
1547  // each entry of the tensor can be
1548  // thought of as an entry in a
1549  // matrix that maps the rolled-out
1550  // rank-2 tensors into rolled-out
1551  // rank-2 tensors. this is the
1552  // format in which we store rank-4
1553  // tensors. determine which
1554  // position the present entry is
1555  // stored in
1556  {
1557  unsigned int base_index[2] ;
1558  if ((indices[0] == 0) && (indices[1] == 0))
1559  base_index[0] = 0;
1560  else if ((indices[0] == 1) && (indices[1] == 1))
1561  base_index[0] = 1;
1562  else
1563  base_index[0] = 2;
1564 
1565  if ((indices[2] == 0) && (indices[3] == 0))
1566  base_index[1] = 0;
1567  else if ((indices[2] == 1) && (indices[3] == 1))
1568  base_index[1] = 1;
1569  else
1570  base_index[1] = 2;
1571 
1572  return data[base_index[0]][base_index[1]];
1573  }
1574 
1575  case 3:
1576  // each entry of the tensor can be
1577  // thought of as an entry in a
1578  // matrix that maps the rolled-out
1579  // rank-2 tensors into rolled-out
1580  // rank-2 tensors. this is the
1581  // format in which we store rank-4
1582  // tensors. determine which
1583  // position the present entry is
1584  // stored in
1585  {
1586  unsigned int base_index[2] ;
1587  if ((indices[0] == 0) && (indices[1] == 0))
1588  base_index[0] = 0;
1589  else if ((indices[0] == 1) && (indices[1] == 1))
1590  base_index[0] = 1;
1591  else if ((indices[0] == 2) && (indices[1] == 2))
1592  base_index[0] = 2;
1593  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1594  ((indices[0] == 1) && (indices[1] == 0)))
1595  base_index[0] = 3;
1596  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1597  ((indices[0] == 2) && (indices[1] == 0)))
1598  base_index[0] = 4;
1599  else
1600  {
1601  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1602  ((indices[0] == 2) && (indices[1] == 1)),
1603  ExcInternalError());
1604  base_index[0] = 5;
1605  }
1606 
1607  if ((indices[2] == 0) && (indices[3] == 0))
1608  base_index[1] = 0;
1609  else if ((indices[2] == 1) && (indices[3] == 1))
1610  base_index[1] = 1;
1611  else if ((indices[2] == 2) && (indices[3] == 2))
1612  base_index[1] = 2;
1613  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1614  ((indices[2] == 1) && (indices[3] == 0)))
1615  base_index[1] = 3;
1616  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1617  ((indices[2] == 2) && (indices[3] == 0)))
1618  base_index[1] = 4;
1619  else
1620  {
1621  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1622  ((indices[2] == 2) && (indices[3] == 1)),
1623  ExcInternalError());
1624  base_index[1] = 5;
1625  }
1626 
1627  return data[base_index[0]][base_index[1]];
1628  }
1629 
1630  default:
1631  Assert (false, ExcNotImplemented());
1632  }
1633 
1634  static Number dummy;
1635  return dummy;
1636  }
1637 
1638 } // end of namespace internal
1639 
1640 
1641 
1642 template <int rank, int dim, typename Number>
1643 inline
1644 Number &
1646 {
1647  for (unsigned int r=0; r<rank; ++r)
1648  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1649  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1650 }
1651 
1652 
1653 
1654 template <int rank, int dim, typename Number>
1655 inline
1656 Number
1658 (const TableIndices<rank> &indices) const
1659 {
1660  for (unsigned int r=0; r<rank; ++r)
1661  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1662  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1663 }
1664 
1665 
1666 
1667 template <int rank, int dim, typename Number>
1668 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1669 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1670 {
1671  return
1672  internal::SymmetricTensorAccessors::
1673  Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
1674 }
1675 
1676 
1677 
1678 template <int rank, int dim, typename Number>
1679 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
1680 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
1681 {
1682  return
1683  internal::SymmetricTensorAccessors::
1684  Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
1685 }
1686 
1687 
1688 
1689 template <int rank, int dim, typename Number>
1690 inline
1691 Number
1693 {
1694  return operator()(indices);
1695 }
1696 
1697 
1698 
1699 template <int rank, int dim, typename Number>
1700 inline
1701 Number &
1703 {
1704  return operator()(indices);
1705 }
1706 
1707 
1708 
1709 
1710 namespace internal
1711 {
1712  namespace SymmetricTensor
1713  {
1714  template <int dim, typename Number>
1715  unsigned int
1716  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
1717  const unsigned int index)
1718  {
1719  return index;
1720  }
1721 
1722 
1723  template <int dim, typename Number>
1725  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
1726  const unsigned int index)
1727  {
1728  return
1731  }
1732 
1733  }
1734 }
1735 
1736 
1737 
1738 template <int rank, int dim, typename Number>
1739 inline
1740 Number
1741 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
1742 {
1743  AssertIndexRange (index, n_independent_components);
1744  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1745 }
1746 
1747 
1748 
1749 template <int rank, int dim, typename Number>
1750 inline
1751 Number &
1752 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
1753 {
1754  AssertIndexRange (index, n_independent_components);
1755  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1756 }
1757 
1758 
1759 
1760 namespace internal
1761 {
1762  template <int dim, typename Number>
1763  inline
1766  {
1767  switch (dim)
1768  {
1769  case 1:
1770  return numbers::NumberTraits<Number>::abs(data[0]);
1771 
1772  case 2:
1773  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1776 
1777  case 3:
1778  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1784 
1785  default:
1786  {
1787  typename numbers::NumberTraits<Number>::real_type return_value
1789 
1790  for (unsigned int d=0; d<dim; ++d)
1791  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
1792  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
1793  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
1794 
1795  return std::sqrt(return_value);
1796  }
1797  }
1798  }
1799 
1800 
1801 
1802  template <int dim, typename Number>
1803  inline
1806  {
1807  switch (dim)
1808  {
1809  case 1:
1810  return numbers::NumberTraits<Number>::abs (data[0][0]);
1811 
1812  default:
1813  {
1814  typename numbers::NumberTraits<Number>::real_type return_value
1816 
1817  const unsigned int n_independent_components = data.dimension;
1818 
1819  for (unsigned int i=0; i<dim; ++i)
1820  for (unsigned int j=0; j<dim; ++j)
1821  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
1822  for (unsigned int i=0; i<dim; ++i)
1823  for (unsigned int j=dim; j<n_independent_components; ++j)
1824  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1825  for (unsigned int i=dim; i<n_independent_components; ++i)
1826  for (unsigned int j=0; j<dim; ++j)
1827  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1828  for (unsigned int i=dim; i<n_independent_components; ++i)
1829  for (unsigned int j=dim; j<n_independent_components; ++j)
1830  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1831 
1832  return std::sqrt(return_value);
1833  }
1834  }
1835  }
1836 
1837 } // end of namespace internal
1838 
1839 
1840 
1841 template <int rank, int dim, typename Number>
1842 inline
1845 {
1846  return internal::compute_norm<dim,Number> (data);
1847 }
1848 
1849 
1850 
1851 namespace internal
1852 {
1853  namespace SymmetricTensor
1854  {
1855  namespace
1856  {
1857  // a function to do the unrolling from a set of indices to a
1858  // scalar index into the array in which we store the elements of
1859  // a symmetric tensor
1860  //
1861  // this function is for rank-2 tensors
1862  template <int dim>
1863  inline
1864  unsigned int
1865  component_to_unrolled_index
1866  (const TableIndices<2> &indices)
1867  {
1868  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
1869  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
1870 
1871  switch (dim)
1872  {
1873  case 1:
1874  {
1875  return 0;
1876  }
1877 
1878  case 2:
1879  {
1880  static const unsigned int table[2][2] = {{0, 2},
1881  {2, 1}
1882  };
1883  return table[indices[0]][indices[1]];
1884  }
1885 
1886  case 3:
1887  {
1888  static const unsigned int table[3][3] = {{0, 3, 4},
1889  {3, 1, 5},
1890  {4, 5, 2}
1891  };
1892  return table[indices[0]][indices[1]];
1893  }
1894 
1895  case 4:
1896  {
1897  static const unsigned int table[4][4] = {{0, 4, 5, 6},
1898  {4, 1, 7, 8},
1899  {5, 7, 2, 9},
1900  {6, 8, 9, 3}
1901  };
1902  return table[indices[0]][indices[1]];
1903  }
1904 
1905  default:
1906  // for the remainder, manually figure out the numbering
1907  {
1908  if (indices[0] == indices[1])
1909  return indices[0];
1910 
1911  TableIndices<2> sorted_indices (indices);
1912  sorted_indices.sort ();
1913 
1914  for (unsigned int d=0, c=0; d<dim; ++d)
1915  for (unsigned int e=d+1; e<dim; ++e, ++c)
1916  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1917  return dim+c;
1918 
1919  // should never get here:
1920  Assert(false, ExcInternalError());
1921  return 0;
1922  }
1923  }
1924  }
1925 
1926  // a function to do the unrolling from a set of indices to a
1927  // scalar index into the array in which we store the elements of
1928  // a symmetric tensor
1929  //
1930  // this function is for tensors of ranks not already handled
1931  // above
1932  template <int dim, int rank>
1933  inline
1934  unsigned int
1935  component_to_unrolled_index
1936  (const TableIndices<rank> &indices)
1937  {
1938  (void)indices;
1939  Assert (false, ExcNotImplemented());
1941  }
1942  }
1943  }
1944 }
1945 
1946 
1947 template <int rank, int dim, typename Number>
1948 inline
1949 unsigned int
1951 (const TableIndices<rank> &indices)
1952 {
1953  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
1954 }
1955 
1956 
1957 
1958 namespace internal
1959 {
1960  namespace SymmetricTensor
1961  {
1962  namespace
1963  {
1964  // a function to do the inverse of the unrolling from a set of
1965  // indices to a scalar index into the array in which we store
1966  // the elements of a symmetric tensor. in other words, it goes
1967  // from the scalar index into the array to a set of indices of
1968  // the tensor
1969  //
1970  // this function is for rank-2 tensors
1971  template <int dim>
1972  inline
1974  unrolled_to_component_indices
1975  (const unsigned int i,
1976  const int2type<2> &)
1977  {
1980  switch (dim)
1981  {
1982  case 1:
1983  {
1984  return TableIndices<2>(0,0);
1985  }
1986 
1987  case 2:
1988  {
1989  const TableIndices<2> table[3] =
1990  {
1991  TableIndices<2> (0,0),
1992  TableIndices<2> (1,1),
1993  TableIndices<2> (0,1)
1994  };
1995  return table[i];
1996  }
1997 
1998  case 3:
1999  {
2000  const TableIndices<2> table[6] =
2001  {
2002  TableIndices<2> (0,0),
2003  TableIndices<2> (1,1),
2004  TableIndices<2> (2,2),
2005  TableIndices<2> (0,1),
2006  TableIndices<2> (0,2),
2007  TableIndices<2> (1,2)
2008  };
2009  return table[i];
2010  }
2011 
2012  default:
2013  if (i<dim)
2014  return TableIndices<2> (i,i);
2015 
2016  for (unsigned int d=0, c=0; d<dim; ++d)
2017  for (unsigned int e=d+1; e<dim; ++e, ++c)
2018  if (c==i)
2019  return TableIndices<2>(d,e);
2020 
2021  // should never get here:
2022  Assert(false, ExcInternalError());
2023  return TableIndices<2>(0, 0);
2024  }
2025  }
2026 
2027  // a function to do the inverse of the unrolling from a set of
2028  // indices to a scalar index into the array in which we store
2029  // the elements of a symmetric tensor. in other words, it goes
2030  // from the scalar index into the array to a set of indices of
2031  // the tensor
2032  //
2033  // this function is for tensors of a rank not already handled
2034  // above
2035  template <int dim, int rank>
2036  inline
2038  unrolled_to_component_indices
2039  (const unsigned int i,
2040  const int2type<rank> &)
2041  {
2042  (void)i;
2045  Assert (false, ExcNotImplemented());
2046  return TableIndices<rank>();
2047  }
2048 
2049  }
2050  }
2051 }
2052 
2053 template <int rank, int dim, typename Number>
2054 inline
2057 (const unsigned int i)
2058 {
2059  return
2060  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2062 }
2063 
2064 
2065 
2066 template <int rank, int dim, typename Number>
2067 template <class Archive>
2068 inline
2069 void
2070 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2071 {
2072  ar &data;
2073 }
2074 
2075 
2076 #endif // DOXYGEN
2077 
2078 /* ----------------- Non-member functions operating on tensors. ------------ */
2079 
2080 
2087 template <int rank, int dim, typename Number, typename OtherNumber>
2088 inline
2091  const Tensor<rank, dim, OtherNumber> &right)
2092 {
2093  return Tensor<rank, dim, Number>(left) + right;
2094 }
2095 
2096 
2103 template <int rank, int dim, typename Number, typename OtherNumber>
2104 inline
2108 {
2109  return left + Tensor<rank, dim, OtherNumber>(right);
2110 }
2111 
2112 
2119 template <int rank, int dim, typename Number, typename OtherNumber>
2120 inline
2123  const Tensor<rank, dim, OtherNumber> &right)
2124 {
2125  return Tensor<rank, dim, Number>(left) - right;
2126 }
2127 
2128 
2135 template <int rank, int dim, typename Number, typename OtherNumber>
2136 inline
2140 {
2141  return left - Tensor<rank, dim, OtherNumber>(right);
2142 }
2143 
2144 
2145 
2159 template <int dim, typename Number>
2160 inline
2162 {
2163  switch (dim)
2164  {
2165  case 1:
2166  return t.data[0];
2167  case 2:
2168  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2169  case 3:
2170  {
2171  // in analogy to general tensors, but
2172  // there's something to be simplified for
2173  // the present case
2174  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2175  return ( tmp + tmp
2176  +t.data[0]*t.data[1]*t.data[2]
2177  -t.data[0]*t.data[5]*t.data[5]
2178  -t.data[1]*t.data[4]*t.data[4]
2179  -t.data[2]*t.data[3]*t.data[3]);
2180  }
2181  default:
2182  Assert (false, ExcNotImplemented());
2183  return 0;
2184  }
2185 }
2186 
2187 
2188 
2198 template <int dim, typename Number>
2199 inline
2201 {
2202  return determinant (t);
2203 }
2204 
2205 
2206 
2214 template <int dim, typename Number>
2216 {
2217  Number t = d.data[0];
2218  for (unsigned int i=1; i<dim; ++i)
2219  t += d.data[i];
2220  return t;
2221 }
2222 
2223 
2233 template <int dim, typename Number>
2234 inline
2236 {
2237  return trace (t);
2238 }
2239 
2240 
2253 template <typename Number>
2254 inline
2256 {
2257  return 0;
2258 }
2259 
2260 
2261 
2282 template <typename Number>
2283 inline
2285 {
2286  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2287 }
2288 
2289 
2290 
2300 template <typename Number>
2301 inline
2303 {
2304  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2305  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2306 }
2307 
2308 
2309 
2310 
2320 template <int rank, int dim, typename Number>
2321 inline
2324 {
2325  return t;
2326 }
2327 
2328 
2329 
2339 template <int dim, typename Number>
2340 inline
2343 {
2345 
2346  // subtract scaled trace from the diagonal
2347  const Number tr = trace(t) / dim;
2348  for (unsigned int i=0; i<dim; ++i)
2349  tmp.data[i] -= tr;
2350 
2351  return tmp;
2352 }
2353 
2354 
2355 
2363 template <int dim, typename Number>
2364 inline
2366 unit_symmetric_tensor ()
2367 {
2368  // create a default constructed matrix filled with
2369  // zeros, then set the diagonal elements to one
2371  switch (dim)
2372  {
2373  case 1:
2374  tmp.data[0] = 1;
2375  break;
2376  case 2:
2377  tmp.data[0] = tmp.data[1] = 1;
2378  break;
2379  case 3:
2380  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
2381  break;
2382  default:
2383  for (unsigned int d=0; d<dim; ++d)
2384  tmp.data[d] = 1;
2385  }
2386  return tmp;
2387 }
2388 
2389 
2390 
2399 template <int dim>
2400 inline
2402 unit_symmetric_tensor ()
2403 {
2404  return unit_symmetric_tensor<dim,double>();
2405 }
2406 
2407 
2408 
2423 template <int dim, typename Number>
2424 inline
2426 deviator_tensor ()
2427 {
2429 
2430  // fill the elements treating the diagonal
2431  for (unsigned int i=0; i<dim; ++i)
2432  for (unsigned int j=0; j<dim; ++j)
2433  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
2434 
2435  // then fill the ones that copy over the
2436  // non-diagonal elements. note that during
2437  // the double-contraction, we handle the
2438  // off-diagonal elements twice, so simply
2439  // copying requires a weight of 1/2
2440  for (unsigned int i=dim;
2441  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2442  ++i)
2443  tmp.data[i][i] = 0.5;
2444 
2445  return tmp;
2446 }
2447 
2448 
2449 
2464 template <int dim>
2465 inline
2467 deviator_tensor ()
2468 {
2469  return deviator_tensor<dim,double>();
2470 }
2471 
2472 
2473 
2496 template <int dim, typename Number>
2497 inline
2499 identity_tensor ()
2500 {
2502 
2503  // fill the elements treating the diagonal
2504  for (unsigned int i=0; i<dim; ++i)
2505  tmp.data[i][i] = 1;
2506 
2507  // then fill the ones that copy over the
2508  // non-diagonal elements. note that during
2509  // the double-contraction, we handle the
2510  // off-diagonal elements twice, so simply
2511  // copying requires a weight of 1/2
2512  for (unsigned int i=dim;
2513  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2514  ++i)
2515  tmp.data[i][i] = 0.5;
2516 
2517  return tmp;
2518 }
2519 
2520 
2521 
2543 template <int dim>
2544 inline
2546 identity_tensor ()
2547 {
2548  return identity_tensor<dim,double>();
2549 }
2550 
2551 
2552 
2563 template <int dim, typename Number>
2564 inline
2567 {
2568  // if desired, take over the
2569  // inversion of a 4x4 tensor
2570  // from the FullMatrix
2571  AssertThrow (false, ExcNotImplemented());
2572 
2574 }
2575 
2576 
2577 
2578 #ifndef DOXYGEN
2579 
2580 template <typename Number>
2581 inline
2583 invert (const SymmetricTensor<2,1,Number> &t)
2584 {
2586 
2587  tmp[0][0] = 1.0/t[0][0];
2588 
2589  return tmp;
2590 }
2591 
2592 
2593 
2594 template <typename Number>
2595 inline
2597 invert (const SymmetricTensor<2,2,Number> &t)
2598 {
2600 
2601  // Sympy result: ([
2602  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
2603  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
2604  const TableIndices<2> idx_00 (0,0);
2605  const TableIndices<2> idx_01 (0,1);
2606  const TableIndices<2> idx_11 (1,1);
2607  const Number inv_det_t
2608  = 1.0/(t[idx_00]*t[idx_11]
2609  - t[idx_01]*t[idx_01]);
2610  tmp[idx_00] = t[idx_11];
2611  tmp[idx_01] = -t[idx_01];
2612  tmp[idx_11] = t[idx_00];
2613  tmp *= inv_det_t;
2614 
2615  return tmp;
2616 }
2617 
2618 
2619 
2620 template <typename Number>
2621 inline
2623 invert (const SymmetricTensor<2,3,Number> &t)
2624 {
2626 
2627  // Sympy result: ([
2628  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2629  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2630  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
2631  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2632  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2633  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
2634  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2635  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
2636  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
2637  const TableIndices<2> idx_00 (0,0);
2638  const TableIndices<2> idx_01 (0,1);
2639  const TableIndices<2> idx_02 (0,2);
2640  const TableIndices<2> idx_11 (1,1);
2641  const TableIndices<2> idx_12 (1,2);
2642  const TableIndices<2> idx_22 (2,2);
2643  const Number inv_det_t
2644  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
2645  - t[idx_00]*t[idx_12]*t[idx_12]
2646  - t[idx_01]*t[idx_01]*t[idx_22]
2647  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
2648  - t[idx_02]*t[idx_02]*t[idx_11]);
2649  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
2650  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
2651  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
2652  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
2653  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
2654  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
2655  tmp *= inv_det_t;
2656 
2657  return tmp;
2658 }
2659 
2660 #endif /* DOXYGEN */
2661 
2662 
2663 
2677 template <int dim, typename Number>
2678 inline
2681 {
2683  switch (dim)
2684  {
2685  case 1:
2686  tmp.data[0][0] = 1./t.data[0][0];
2687  break;
2688  case 2:
2689 
2690  // inverting this tensor is a little more
2691  // complicated than necessary, since we
2692  // store the data of 't' as a 3x3 matrix
2693  // t.data, but the product between a rank-4
2694  // and a rank-2 tensor is really not the
2695  // product between this matrix and the
2696  // 3-vector of a rhs, but rather
2697  //
2698  // B.vec = t.data * mult * A.vec
2699  //
2700  // where mult is a 3x3 matrix with
2701  // entries [[1,0,0],[0,1,0],[0,0,2]] to
2702  // capture the fact that we need to add up
2703  // both the c_ij12*a_12 and the c_ij21*a_21
2704  // terms
2705  //
2706  // in addition, in this scheme, the
2707  // identity tensor has the matrix
2708  // representation mult^-1.
2709  //
2710  // the inverse of 't' therefore has the
2711  // matrix representation
2712  //
2713  // inv.data = mult^-1 * t.data^-1 * mult^-1
2714  //
2715  // in order to compute it, let's first
2716  // compute the inverse of t.data and put it
2717  // into tmp.data; at the end of the
2718  // function we then scale the last row and
2719  // column of the inverse by 1/2,
2720  // corresponding to the left and right
2721  // multiplication with mult^-1
2722  {
2723  const Number t4 = t.data[0][0]*t.data[1][1],
2724  t6 = t.data[0][0]*t.data[1][2],
2725  t8 = t.data[0][1]*t.data[1][0],
2726  t00 = t.data[0][2]*t.data[1][0],
2727  t01 = t.data[0][1]*t.data[2][0],
2728  t04 = t.data[0][2]*t.data[2][0],
2729  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
2730  t8*t.data[2][2]+t00*t.data[2][1]+
2731  t01*t.data[1][2]-t04*t.data[1][1]);
2732  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
2733  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
2734  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
2735  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
2736  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
2737  tmp.data[1][2] = -(t6-t00)*t07;
2738  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
2739  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
2740  tmp.data[2][2] = (t4-t8)*t07;
2741 
2742  // scale last row and column as mentioned
2743  // above
2744  tmp.data[2][0] /= 2;
2745  tmp.data[2][1] /= 2;
2746  tmp.data[0][2] /= 2;
2747  tmp.data[1][2] /= 2;
2748  tmp.data[2][2] /= 4;
2749  }
2750  break;
2751  default:
2752  Assert (false, ExcNotImplemented());
2753  }
2754  return tmp;
2755 }
2756 
2757 
2758 
2772 template <>
2774 invert (const SymmetricTensor<4,3,double> &t);
2775 // this function is implemented in the .cc file for double data types
2776 
2777 
2778 
2793 template <int dim, typename Number>
2794 inline
2798 {
2800 
2801  // fill only the elements really needed
2802  for (unsigned int i=0; i<dim; ++i)
2803  for (unsigned int j=i; j<dim; ++j)
2804  for (unsigned int k=0; k<dim; ++k)
2805  for (unsigned int l=k; l<dim; ++l)
2806  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
2807 
2808  return tmp;
2809 }
2810 
2811 
2812 
2821 template <int dim,typename Number>
2822 inline
2825 {
2826  Number array[(dim*dim+dim)/2];
2827  for (unsigned int d=0; d<dim; ++d)
2828  array[d] = t[d][d];
2829  for (unsigned int d=0, c=0; d<dim; ++d)
2830  for (unsigned int e=d+1; e<dim; ++e, ++c)
2831  array[dim+c] = (t[d][e]+t[e][d])*0.5;
2832  return SymmetricTensor<2,dim,Number>(array);
2833 }
2834 
2835 
2836 
2844 template <int rank, int dim, typename Number>
2845 inline
2848  const Number factor)
2849 {
2851  tt *= factor;
2852  return tt;
2853 }
2854 
2855 
2856 
2864 template <int rank, int dim, typename Number>
2865 inline
2867 operator * (const Number factor,
2869 {
2870  // simply forward to the other operator
2871  return t*factor;
2872 }
2873 
2874 
2875 #ifndef DEAL_II_WITH_CXX11
2876 
2877 template <typename T, typename U, int rank, int dim>
2878 struct ProductType<T,SymmetricTensor<rank,dim,U> >
2879 {
2881 };
2882 
2883 template <typename T, typename U, int rank, int dim>
2884 struct ProductType<SymmetricTensor<rank,dim,T>,U>
2885 {
2887 };
2888 
2889 #endif
2890 
2891 
2892 
2918 template <int rank, int dim, typename Number, typename OtherNumber>
2919 inline
2922  const OtherNumber factor)
2923 {
2924  // form the product. we have to convert the two factors into the final
2925  // type via explicit casts because, for awkward reasons, the C++
2926  // standard committee saw it fit to not define an
2927  // operator*(float,std::complex<double>)
2928  // (as well as with switched arguments and double<->float).
2929  typedef typename ProductType<Number,OtherNumber>::type product_type;
2931  // we used to shorten the following by 'tt *= product_type(factor);'
2932  // which requires that a converting constructor
2933  // 'product_type::product_type(const OtherNumber) is defined.
2934  // however, a user-defined constructor is not allowed for aggregates,
2935  // e.g. VectorizedArray. therefore, we work around this issue using a
2936  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
2937  // which we assume to be defined.
2938  product_type new_factor;
2939  new_factor = factor;
2940  tt *= new_factor;
2941  return tt;
2942 }
2943 
2944 
2945 
2954 template <int rank, int dim, typename Number, typename OtherNumber>
2955 inline
2957 operator * (const Number factor,
2959 {
2960  // simply forward to the other operator with switched arguments
2961  return (t*factor);
2962 }
2963 
2964 
2965 
2971 template <int rank, int dim, typename Number>
2972 inline
2975  const Number factor)
2976 {
2978  tt /= factor;
2979  return tt;
2980 }
2981 
2982 
2983 
2990 template <int rank, int dim>
2991 inline
2993 operator * (const SymmetricTensor<rank,dim> &t,
2994  const double factor)
2995 {
2997  tt *= factor;
2998  return tt;
2999 }
3000 
3001 
3002 
3009 template <int rank, int dim>
3010 inline
3012 operator * (const double factor,
3013  const SymmetricTensor<rank,dim> &t)
3014 {
3016  tt *= factor;
3017  return tt;
3018 }
3019 
3020 
3021 
3027 template <int rank, int dim>
3028 inline
3030 operator / (const SymmetricTensor<rank,dim> &t,
3031  const double factor)
3032 {
3034  tt /= factor;
3035  return tt;
3036 }
3037 
3047 template <int dim, typename Number>
3048 inline
3049 Number
3052 {
3053  return (t1*t2);
3054 }
3055 
3056 
3066 template <int dim, typename Number>
3067 inline
3068 Number
3070  const Tensor<2,dim,Number> &t2)
3071 {
3072  Number s = 0;
3073  for (unsigned int i=0; i<dim; ++i)
3074  for (unsigned int j=0; j<dim; ++j)
3075  s += t1[i][j] * t2[i][j];
3076  return s;
3077 }
3078 
3079 
3089 template <int dim, typename Number>
3090 inline
3091 Number
3094 {
3095  return scalar_product(t2, t1);
3096 }
3097 
3098 
3114 template <typename Number>
3115 inline
3116 void
3118  const SymmetricTensor<4,1,Number> &t,
3119  const SymmetricTensor<2,1,Number> &s)
3120 {
3121  tmp[0][0] = t[0][0][0][0] * s[0][0];
3122 }
3123 
3124 
3125 
3141 template <typename Number>
3142 inline
3143 void
3145  const SymmetricTensor<2,1,Number> &s,
3146  const SymmetricTensor<4,1,Number> &t)
3147 {
3148  tmp[0][0] = t[0][0][0][0] * s[0][0];
3149 }
3150 
3151 
3152 
3167 template <typename Number>
3168 inline
3169 void
3171  const SymmetricTensor<4,2,Number> &t,
3172  const SymmetricTensor<2,2,Number> &s)
3173 {
3174  const unsigned int dim = 2;
3175 
3176  for (unsigned int i=0; i<dim; ++i)
3177  for (unsigned int j=i; j<dim; ++j)
3178  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3179  t[i][j][1][1] * s[1][1] +
3180  2 * t[i][j][0][1] * s[0][1];
3181 }
3182 
3183 
3184 
3200 template <typename Number>
3201 inline
3202 void
3204  const SymmetricTensor<2,2,Number> &s,
3205  const SymmetricTensor<4,2,Number> &t)
3206 {
3207  const unsigned int dim = 2;
3208 
3209  for (unsigned int i=0; i<dim; ++i)
3210  for (unsigned int j=i; j<dim; ++j)
3211  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3212  s[1][1] * t[1][1][i][j] +
3213  2 * s[0][1] * t[0][1][i][j];
3214 }
3215 
3216 
3217 
3233 template <typename Number>
3234 inline
3235 void
3237  const SymmetricTensor<4,3,Number> &t,
3238  const SymmetricTensor<2,3,Number> &s)
3239 {
3240  const unsigned int dim = 3;
3241 
3242  for (unsigned int i=0; i<dim; ++i)
3243  for (unsigned int j=i; j<dim; ++j)
3244  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3245  t[i][j][1][1] * s[1][1] +
3246  t[i][j][2][2] * s[2][2] +
3247  2 * t[i][j][0][1] * s[0][1] +
3248  2 * t[i][j][0][2] * s[0][2] +
3249  2 * t[i][j][1][2] * s[1][2];
3250 }
3251 
3252 
3253 
3269 template <typename Number>
3270 inline
3271 void
3273  const SymmetricTensor<2,3,Number> &s,
3274  const SymmetricTensor<4,3,Number> &t)
3275 {
3276  const unsigned int dim = 3;
3277 
3278  for (unsigned int i=0; i<dim; ++i)
3279  for (unsigned int j=i; j<dim; ++j)
3280  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3281  s[1][1] * t[1][1][i][j] +
3282  s[2][2] * t[2][2][i][j] +
3283  2 * s[0][1] * t[0][1][i][j] +
3284  2 * s[0][2] * t[0][2][i][j] +
3285  2 * s[1][2] * t[1][2][i][j];
3286 }
3287 
3288 
3289 
3297 template <int dim, typename Number>
3299 operator * (const SymmetricTensor<2,dim,Number> &src1,
3300  const Tensor<1,dim,Number> &src2)
3301 {
3302  Tensor<1,dim,Number> dest;
3303  for (unsigned int i=0; i<dim; ++i)
3304  for (unsigned int j=0; j<dim; ++j)
3305  dest[i] += src1[i][j] * src2[j];
3306  return dest;
3307 }
3308 
3309 
3317 template <int dim, typename Number>
3319 operator * (const Tensor<1,dim,Number> &src1,
3320  const SymmetricTensor<2,dim,Number> &src2)
3321 {
3322  // this is easy for symmetric tensors:
3323  return src2 * src1;
3324 }
3325 
3326 
3336 template <int dim, typename Number>
3337 inline
3338 std::ostream &operator << (std::ostream &out,
3340 {
3341  //make out lives a bit simpler by outputing
3342  //the tensor through the operator for the
3343  //general Tensor class
3345 
3346  for (unsigned int i=0; i<dim; ++i)
3347  for (unsigned int j=0; j<dim; ++j)
3348  tt[i][j] = t[i][j];
3349 
3350  return out << tt;
3351 }
3352 
3353 
3354 
3364 template <int dim, typename Number>
3365 inline
3366 std::ostream &operator << (std::ostream &out,
3368 {
3369  //make out lives a bit simpler by outputing
3370  //the tensor through the operator for the
3371  //general Tensor class
3373 
3374  for (unsigned int i=0; i<dim; ++i)
3375  for (unsigned int j=0; j<dim; ++j)
3376  for (unsigned int k=0; k<dim; ++k)
3377  for (unsigned int l=0; l<dim; ++l)
3378  tt[i][j][k][l] = t[i][j][k][l];
3379 
3380  return out << tt;
3381 }
3382 
3383 
3384 DEAL_II_NAMESPACE_CLOSE
3385 
3386 #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:342
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()