Reference documentation for deal.II version Git 7edbc71 2017-07-25 01:57:43 -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 
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)
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,
98  case 1:
99  return TableIndices<4>(previous_indices[0],
100  new_index,
103  case 2:
104  return TableIndices<4>(previous_indices[0],
105  previous_indices[1],
106  new_index,
108  case 3:
109  return TableIndices<4>(previous_indices[0],
110  previous_indices[1],
111  previous_indices[2],
112  new_index);
113  }
114  Assert (false, ExcInternalError());
115  return TableIndices<4>();
116  }
117 
118 
127  template <int rank1, int rank2, int dim, typename Number>
129  {
130  typedef ::SymmetricTensor<rank1+rank2-4,dim,Number> type;
131  };
132 
133 
142  template <int dim, typename Number>
143  struct double_contraction_result<2,2,dim,Number>
144  {
145  typedef Number type;
146  };
147 
148 
149 
162  template <int rank, int dim, typename Number>
163  struct StorageType;
164 
168  template <int dim, typename Number>
169  struct StorageType<2,dim,Number>
170  {
175  static const unsigned int
176  n_independent_components = (dim*dim + dim)/2;
177 
182  };
183 
184 
185 
189  template <int dim, typename Number>
190  struct StorageType<4,dim,Number>
191  {
197  static const unsigned int
198  n_rank2_components = (dim*dim + dim)/2;
199 
203  static const unsigned int
204  n_independent_components = (n_rank2_components *
206 
214  };
215 
216 
217 
222  template <int rank, int dim, bool constness, typename Number>
224 
231  template <int rank, int dim, typename Number>
232  struct AccessorTypes<rank,dim,true,Number>
233  {
234  typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
235 
236  typedef Number reference;
237  };
238 
245  template <int rank, int dim, typename Number>
246  struct AccessorTypes<rank,dim,false,Number>
247  {
248  typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
249 
250  typedef Number &reference;
251  };
252 
253 
288  template <int rank, int dim, bool constness, int P, typename Number>
289  class Accessor
290  {
291  public:
295  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
296  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
297 
298  private:
317  Accessor (tensor_type &tensor,
318  const TableIndices<rank> &previous_indices);
319 
323  Accessor (const Accessor &) = default;
324 
325  public:
326 
330  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
331 
335  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i) const;
336 
337  private:
341  tensor_type &tensor;
342  const TableIndices<rank> previous_indices;
343 
344  // declare some other classes
345  // as friends. make sure to
346  // work around bugs in some
347  // compilers
348  template <int,int,typename> friend class ::SymmetricTensor;
349  template <int,int,bool,int,typename>
350  friend class Accessor;
351 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
352  friend class ::SymmetricTensor<rank,dim,Number>;
353  friend class Accessor<rank,dim,constness,P+1,Number>;
354 # endif
355  };
356 
357 
358 
368  template <int rank, int dim, bool constness, typename Number>
369  class Accessor<rank,dim,constness,1,Number>
370  {
371  public:
375  typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
376  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
377 
378  private:
400  Accessor (tensor_type &tensor,
401  const TableIndices<rank> &previous_indices);
402 
406  Accessor () = delete;
407 
411  Accessor (const Accessor &) = default;
412 
413  public:
414 
418  reference operator [] (const unsigned int);
419 
423  reference operator [] (const unsigned int) const;
424 
425  private:
429  tensor_type &tensor;
430  const TableIndices<rank> previous_indices;
431 
432  // declare some other classes
433  // as friends. make sure to
434  // work around bugs in some
435  // compilers
436  template <int,int,typename> friend class ::SymmetricTensor;
437  template <int,int,bool,int,typename>
438  friend class SymmetricTensorAccessors::Accessor;
439 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
440  friend class ::SymmetricTensor<rank,dim,Number>;
441  friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
442 # endif
443  };
444  }
445 }
446 
447 
448 
512 template <int rank, int dim, typename Number>
513 class SymmetricTensor
514 {
515 public:
524  static const unsigned int dimension = dim;
525 
531  static const unsigned int n_independent_components
534 
538  SymmetricTensor () = default;
539 
551 
567  SymmetricTensor (const Number (&array) [n_independent_components]);
568 
574  template <typename OtherNumber>
575  explicit
577 
584  SymmetricTensor &operator = (const Number d);
585 
590  operator Tensor<rank,dim,Number> () const;
591 
595  bool operator == (const SymmetricTensor &) const;
596 
600  bool operator != (const SymmetricTensor &) const;
601 
606 
611 
616  SymmetricTensor &operator *= (const Number factor);
617 
621  SymmetricTensor &operator /= (const Number factor);
622 
628 
634 
639 
666 
673 
677  Number &operator() (const TableIndices<rank> &indices);
678 
682  Number operator() (const TableIndices<rank> &indices) const;
683 
688  internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
689  operator [] (const unsigned int row) const;
690 
695  internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
696  operator [] (const unsigned int row);
697 
703  Number
704  operator [] (const TableIndices<rank> &indices) const;
705 
711  Number &
712  operator [] (const TableIndices<rank> &indices);
713 
719  Number
720  access_raw_entry (const unsigned int unrolled_index) const;
721 
727  Number &
728  access_raw_entry (const unsigned int unrolled_index);
729 
740  norm () const;
741 
749  static
750  unsigned int
752 
758  static
760  unrolled_to_component_indices (const unsigned int i);
761 
774  void clear ();
775 
780  static std::size_t memory_consumption ();
781 
786  template <class Archive>
787  void serialize(Archive &ar, const unsigned int version);
788 
789 private:
793  typedef
796 
800  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
801 
805  base_tensor_type data;
806 
810  template <int, int, typename> friend class SymmetricTensor;
811 
815  template <int dim2, typename Number2>
816  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
817 
818  template <int dim2, typename Number2>
819  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
820 
821  template <int dim2, typename Number2>
823  deviator (const SymmetricTensor<2,dim2,Number2> &t);
824 
825  template <int dim2, typename Number2>
827 
828  template <int dim2, typename Number2>
830 
831  template <int dim2, typename Number2>
833 
834  template <int dim2, typename Number2>
836 
837  template <int dim2, typename Number2>
839 };
840 
841 
842 
843 // ------------------------- inline functions ------------------------
844 
845 #ifndef DOXYGEN
846 
847 namespace internal
848 {
849  namespace SymmetricTensorAccessors
850  {
851  template <int rank, int dim, bool constness, int P, typename Number>
852  Accessor<rank,dim,constness,P,Number>::
853  Accessor (tensor_type &tensor,
854  const TableIndices<rank> &previous_indices)
855  :
856  tensor (tensor),
857  previous_indices (previous_indices)
858  {}
859 
860 
861 
862  template <int rank, int dim, bool constness, int P, typename Number>
863  Accessor<rank,dim,constness,P-1,Number>
864  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
865  {
866  return Accessor<rank,dim,constness,P-1,Number> (tensor,
867  merge (previous_indices, i, rank-P));
868  }
869 
870 
871 
872  template <int rank, int dim, bool constness, int P, typename Number>
873  Accessor<rank,dim,constness,P-1,Number>
874  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i) const
875  {
876  return Accessor<rank,dim,constness,P-1,Number> (tensor,
877  merge (previous_indices, i, rank-P));
878  }
879 
880 
881 
882  template <int rank, int dim, bool constness, typename Number>
883  Accessor<rank,dim,constness,1,Number>::
884  Accessor (tensor_type &tensor,
885  const TableIndices<rank> &previous_indices)
886  :
887  tensor (tensor),
888  previous_indices (previous_indices)
889  {}
890 
891 
892 
893  template <int rank, int dim, bool constness, typename Number>
894  typename Accessor<rank,dim,constness,1,Number>::reference
895  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
896  {
897  return tensor(merge (previous_indices, i, rank-1));
898  }
899 
900 
901  template <int rank, int dim, bool constness, typename Number>
902  typename Accessor<rank,dim,constness,1,Number>::reference
903  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i) const
904  {
905  return tensor(merge (previous_indices, i, rank-1));
906  }
907  }
908 }
909 
910 
911 
912 template <int rank, int dim, typename Number>
913 inline
915 {
916  Assert (rank == 2, ExcNotImplemented());
917  switch (dim)
918  {
919  case 2:
920  Assert (t[0][1] == t[1][0], ExcInternalError());
921 
922  data[0] = t[0][0];
923  data[1] = t[1][1];
924  data[2] = t[0][1];
925 
926  break;
927  case 3:
928  Assert (t[0][1] == t[1][0], ExcInternalError());
929  Assert (t[0][2] == t[2][0], ExcInternalError());
930  Assert (t[1][2] == t[2][1], ExcInternalError());
931 
932  data[0] = t[0][0];
933  data[1] = t[1][1];
934  data[2] = t[2][2];
935  data[3] = t[0][1];
936  data[4] = t[0][2];
937  data[5] = t[1][2];
938 
939  break;
940  default:
941  for (unsigned int d=0; d<dim; ++d)
942  for (unsigned int e=0; e<d; ++e)
943  Assert(t[d][e] == t[e][d], ExcInternalError());
944 
945  for (unsigned int d=0; d<dim; ++d)
946  data[d] = t[d][d];
947 
948  for (unsigned int d=0, c=0; d<dim; ++d)
949  for (unsigned int e=d+1; e<dim; ++e, ++c)
950  data[dim+c] = t[d][e];
951  }
952 }
953 
954 
955 
956 template <int rank, int dim, typename Number>
957 template <typename OtherNumber>
958 inline
961 {
962  for (unsigned int i=0; i<base_tensor_type::dimension; ++i)
963  data[i] = initializer.data[i];
964 }
965 
966 
967 
968 
969 template <int rank, int dim, typename Number>
970 inline
971 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
972  :
973  data (*reinterpret_cast<const typename base_tensor_type::array_type *>(array))
974 {
975  // ensure that the reinterpret_cast above actually works
976  Assert (sizeof(typename base_tensor_type::array_type)
977  == sizeof(array),
978  ExcInternalError());
979 }
980 
981 
982 
983 template <int rank, int dim, typename Number>
984 inline
987 {
988  Assert (d==Number(), ExcMessage ("Only assignment with zero is allowed"));
989  (void) d;
990 
991  data = 0;
992 
993  return *this;
994 }
995 
996 
997 namespace internal
998 {
999  namespace SymmetricTensor
1000  {
1001  template <int dim, typename Number>
1003  convert_to_tensor (const ::SymmetricTensor<2,dim,Number> &s)
1004  {
1006 
1007  // diagonal entries are stored first
1008  for (unsigned int d=0; d<dim; ++d)
1009  t[d][d] = s.access_raw_entry(d);
1010 
1011  // off-diagonal entries come next, row by row
1012  for (unsigned int d=0, c=0; d<dim; ++d)
1013  for (unsigned int e=d+1; e<dim; ++e, ++c)
1014  {
1015  t[d][e] = s.access_raw_entry(dim+c);
1016  t[e][d] = s.access_raw_entry(dim+c);
1017  }
1018  return t;
1019  }
1020 
1021 
1022  template <int dim, typename Number>
1024  convert_to_tensor (const ::SymmetricTensor<4,dim,Number> &st)
1025  {
1026  // utilize the symmetry properties of SymmetricTensor<4,dim>
1027  // discussed in the class documentation to avoid accessing all
1028  // independent elements of the input tensor more than once
1030 
1031  for (unsigned int i=0; i<dim; ++i)
1032  for (unsigned int j=i; j<dim; ++j)
1033  for (unsigned int k=0; k<dim; ++k)
1034  for (unsigned int l=k; l<dim; ++l)
1035  t[TableIndices<4>(i,j,k,l)]
1036  = t[TableIndices<4>(i,j,l,k)]
1037  = t[TableIndices<4>(j,i,k,l)]
1038  = t[TableIndices<4>(j,i,l,k)]
1039  = st[TableIndices<4>(i,j,k,l)];
1040 
1041  return t;
1042  }
1043  }
1044 }
1045 
1046 
1047 
1048 template <int rank, int dim, typename Number>
1049 inline
1051 operator Tensor<rank,dim,Number> () const
1052 {
1053  return internal::SymmetricTensor::convert_to_tensor (*this);
1054 }
1055 
1056 
1057 
1058 template <int rank, int dim, typename Number>
1059 inline
1060 bool
1062 (const SymmetricTensor<rank,dim,Number> &t) const
1063 {
1064  return data == t.data;
1065 }
1066 
1067 
1068 
1069 template <int rank, int dim, typename Number>
1070 inline
1071 bool
1072 SymmetricTensor<rank,dim,Number>::operator !=
1073 (const SymmetricTensor<rank,dim,Number> &t) const
1074 {
1075  return data != t.data;
1076 }
1077 
1078 
1079 
1080 template <int rank, int dim, typename Number>
1081 inline
1083 SymmetricTensor<rank,dim,Number>::operator +=
1085 {
1086  data += t.data;
1087  return *this;
1088 }
1089 
1090 
1091 
1092 template <int rank, int dim, typename Number>
1093 inline
1095 SymmetricTensor<rank,dim,Number>::operator -=
1097 {
1098  data -= t.data;
1099  return *this;
1100 }
1101 
1102 
1103 
1104 template <int rank, int dim, typename Number>
1105 inline
1108 {
1109  data *= d;
1110  return *this;
1111 }
1112 
1113 
1114 
1115 template <int rank, int dim, typename Number>
1116 inline
1119 {
1120  data /= d;
1121  return *this;
1122 }
1123 
1124 
1125 
1126 template <int rank, int dim, typename Number>
1127 inline
1130 {
1131  SymmetricTensor tmp = *this;
1132  tmp.data += t.data;
1133  return tmp;
1134 }
1135 
1136 
1137 
1138 template <int rank, int dim, typename Number>
1139 inline
1142 {
1143  SymmetricTensor tmp = *this;
1144  tmp.data -= t.data;
1145  return tmp;
1146 }
1147 
1148 
1149 
1150 template <int rank, int dim, typename Number>
1151 inline
1154 {
1155  SymmetricTensor tmp = *this;
1156  tmp.data = -tmp.data;
1157  return tmp;
1158 }
1159 
1160 
1161 
1162 template <int rank, int dim, typename Number>
1163 inline
1164 void
1166 {
1167  data.clear ();
1168 }
1169 
1170 
1171 
1172 template <int rank, int dim, typename Number>
1173 inline
1174 std::size_t
1176 {
1177  // all memory consists of statically allocated memory of the current
1178  // object, no pointers
1179  return sizeof(SymmetricTensor<rank,dim,Number>);
1180 }
1181 
1182 
1183 
1184 namespace internal
1185 {
1186 
1187  template <int dim, typename Number>
1188  inline
1189  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
1190  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1192  {
1193  switch (dim)
1194  {
1195  case 1:
1196  return data[0] * sdata[0];
1197  default:
1198  // Start with the non-diagonal part to avoid some multiplications by
1199  // 2.
1200  Number sum = data[dim] * sdata[dim];
1201  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1202  sum += data[d] * sdata[d];
1203  sum += sum; // sum = sum * 2.;
1204 
1205  // Now add the contributions from the diagonal
1206  for (unsigned int d=0; d<dim; ++d)
1207  sum += data[d] * sdata[d];
1208  return sum;
1209  }
1210  }
1211 
1212 
1213 
1214  template <int dim, typename Number>
1215  inline
1216  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
1217  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1219  {
1220  const unsigned int data_dim =
1222  Number tmp [data_dim];
1223  for (unsigned int i=0; i<data_dim; ++i)
1224  tmp[i] = perform_double_contraction<dim,Number>(data[i], sdata);
1225  return ::SymmetricTensor<2,dim,Number>(tmp);
1226  }
1227 
1228 
1229 
1230  template <int dim, typename Number>
1231  inline
1233  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1235  {
1237  for (unsigned int i=0; i<tmp.dimension; ++i)
1238  {
1239  // Start with the non-diagonal part
1240  Number sum = data[dim] * sdata[dim][i];
1241  for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
1242  sum += data[d] * sdata[d][i];
1243  sum += sum; // sum = sum * 2.;
1244 
1245  // Now add the contributions from the diagonal
1246  for (unsigned int d=0; d<dim; ++d)
1247  sum += data[d] * sdata[d][i];
1248  tmp[i] = sum;
1249  }
1250  return tmp;
1251  }
1252 
1253 
1254 
1255  template <int dim, typename Number>
1256  inline
1258  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1260  {
1261  const unsigned int data_dim =
1264  for (unsigned int i=0; i<data_dim; ++i)
1265  for (unsigned int j=0; j<data_dim; ++j)
1266  {
1267  // Start with the non-diagonal part
1268  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1269  tmp[i][j] += data[i][d] * sdata[d][j];
1270  tmp[i][j] += tmp[i][j]; // tmp[i][j] = tmp[i][j] * 2;
1271 
1272  // Now add the contributions from the diagonal
1273  for (unsigned int d=0; d<dim; ++d)
1274  tmp[i][j] += data[i][d] * sdata[d][j];
1275  }
1276  return tmp;
1277  }
1278 
1279 } // end of namespace internal
1280 
1281 
1282 
1283 template <int rank, int dim, typename Number>
1284 inline
1287 {
1288  // need to have two different function calls
1289  // because a scalar and rank-2 tensor are not
1290  // the same data type (see internal function
1291  // above)
1292  return internal::perform_double_contraction<dim,Number> (data, s.data);
1293 }
1294 
1295 
1296 
1297 template <int rank, int dim, typename Number>
1298 inline
1301 {
1304  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
1305  return tmp;
1306 }
1307 
1308 
1309 
1310 // internal namespace to switch between the
1311 // access of different tensors. There used to
1312 // be explicit instantiations before for
1313 // different ranks and dimensions, but since
1314 // we now allow for templates on the data
1315 // type, and since we cannot partially
1316 // specialize the implementation, this got
1317 // into a separate namespace
1318 namespace internal
1319 {
1320  template <int dim, typename Number>
1321  inline
1322  Number &
1323  symmetric_tensor_access (const TableIndices<2> &indices,
1325  {
1326  // 1d is very simple and done first
1327  if (dim == 1)
1328  return data[0];
1329 
1330  // first treat the main diagonal elements, which are stored consecutively
1331  // at the beginning
1332  if (indices[0] == indices[1])
1333  return data[indices[0]];
1334 
1335  // the rest is messier and requires a few switches.
1336  switch (dim)
1337  {
1338  case 2:
1339  // at least for the 2x2 case it is reasonably simple
1340  Assert (((indices[0]==1) && (indices[1]==0)) ||
1341  ((indices[0]==0) && (indices[1]==1)),
1342  ExcInternalError());
1343  return data[2];
1344 
1345  default:
1346  // to do the rest, sort our indices before comparing
1347  {
1348  TableIndices<2> sorted_indices (indices);
1349  sorted_indices.sort ();
1350 
1351  for (unsigned int d=0, c=0; d<dim; ++d)
1352  for (unsigned int e=d+1; e<dim; ++e, ++c)
1353  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1354  return data[dim+c];
1355  Assert (false, ExcInternalError());
1356  }
1357  }
1358 
1359  static Number dummy_but_referenceable = Number();
1360  return dummy_but_referenceable;
1361  }
1362 
1363 
1364 
1365  template <int dim, typename Number>
1366  inline
1367  Number
1368  symmetric_tensor_access (const TableIndices<2> &indices,
1370  {
1371  // 1d is very simple and done first
1372  if (dim == 1)
1373  return data[0];
1374 
1375  // first treat the main diagonal elements, which are stored consecutively
1376  // at the beginning
1377  if (indices[0] == indices[1])
1378  return data[indices[0]];
1379 
1380  // the rest is messier and requires a few switches.
1381  switch (dim)
1382  {
1383  case 2:
1384  // at least for the 2x2 case it is reasonably simple
1385  Assert (((indices[0]==1) && (indices[1]==0)) ||
1386  ((indices[0]==0) && (indices[1]==1)),
1387  ExcInternalError());
1388  return data[2];
1389 
1390  default:
1391  // to do the rest, sort our indices before comparing
1392  {
1393  TableIndices<2> sorted_indices (indices);
1394  sorted_indices.sort ();
1395 
1396  for (unsigned int d=0, c=0; d<dim; ++d)
1397  for (unsigned int e=d+1; e<dim; ++e, ++c)
1398  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1399  return data[dim+c];
1400  Assert (false, ExcInternalError());
1401  }
1402  }
1403 
1404  static Number dummy_but_referenceable = Number();
1405  return dummy_but_referenceable;
1406  }
1407 
1408 
1409 
1410  template <int dim, typename Number>
1411  inline
1412  Number &
1413  symmetric_tensor_access (const TableIndices<4> &indices,
1415  {
1416  switch (dim)
1417  {
1418  case 1:
1419  return data[0][0];
1420 
1421  case 2:
1422  // each entry of the tensor can be
1423  // thought of as an entry in a
1424  // matrix that maps the rolled-out
1425  // rank-2 tensors into rolled-out
1426  // rank-2 tensors. this is the
1427  // format in which we store rank-4
1428  // tensors. determine which
1429  // position the present entry is
1430  // stored in
1431  {
1432  unsigned int base_index[2] ;
1433  if ((indices[0] == 0) && (indices[1] == 0))
1434  base_index[0] = 0;
1435  else if ((indices[0] == 1) && (indices[1] == 1))
1436  base_index[0] = 1;
1437  else
1438  base_index[0] = 2;
1439 
1440  if ((indices[2] == 0) && (indices[3] == 0))
1441  base_index[1] = 0;
1442  else if ((indices[2] == 1) && (indices[3] == 1))
1443  base_index[1] = 1;
1444  else
1445  base_index[1] = 2;
1446 
1447  return data[base_index[0]][base_index[1]];
1448  }
1449 
1450  case 3:
1451  // each entry of the tensor can be
1452  // thought of as an entry in a
1453  // matrix that maps the rolled-out
1454  // rank-2 tensors into rolled-out
1455  // rank-2 tensors. this is the
1456  // format in which we store rank-4
1457  // tensors. determine which
1458  // position the present entry is
1459  // stored in
1460  {
1461  unsigned int base_index[2] ;
1462  if ((indices[0] == 0) && (indices[1] == 0))
1463  base_index[0] = 0;
1464  else if ((indices[0] == 1) && (indices[1] == 1))
1465  base_index[0] = 1;
1466  else if ((indices[0] == 2) && (indices[1] == 2))
1467  base_index[0] = 2;
1468  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1469  ((indices[0] == 1) && (indices[1] == 0)))
1470  base_index[0] = 3;
1471  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1472  ((indices[0] == 2) && (indices[1] == 0)))
1473  base_index[0] = 4;
1474  else
1475  {
1476  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1477  ((indices[0] == 2) && (indices[1] == 1)),
1478  ExcInternalError());
1479  base_index[0] = 5;
1480  }
1481 
1482  if ((indices[2] == 0) && (indices[3] == 0))
1483  base_index[1] = 0;
1484  else if ((indices[2] == 1) && (indices[3] == 1))
1485  base_index[1] = 1;
1486  else if ((indices[2] == 2) && (indices[3] == 2))
1487  base_index[1] = 2;
1488  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1489  ((indices[2] == 1) && (indices[3] == 0)))
1490  base_index[1] = 3;
1491  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1492  ((indices[2] == 2) && (indices[3] == 0)))
1493  base_index[1] = 4;
1494  else
1495  {
1496  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1497  ((indices[2] == 2) && (indices[3] == 1)),
1498  ExcInternalError());
1499  base_index[1] = 5;
1500  }
1501 
1502  return data[base_index[0]][base_index[1]];
1503  }
1504 
1505  default:
1506  Assert (false, ExcNotImplemented());
1507  }
1508 
1509  static Number dummy;
1510  return dummy;
1511  }
1512 
1513 
1514  template <int dim, typename Number>
1515  inline
1516  Number
1517  symmetric_tensor_access (const TableIndices<4> &indices,
1519  {
1520  switch (dim)
1521  {
1522  case 1:
1523  return data[0][0];
1524 
1525  case 2:
1526  // each entry of the tensor can be
1527  // thought of as an entry in a
1528  // matrix that maps the rolled-out
1529  // rank-2 tensors into rolled-out
1530  // rank-2 tensors. this is the
1531  // format in which we store rank-4
1532  // tensors. determine which
1533  // position the present entry is
1534  // stored in
1535  {
1536  unsigned int base_index[2] ;
1537  if ((indices[0] == 0) && (indices[1] == 0))
1538  base_index[0] = 0;
1539  else if ((indices[0] == 1) && (indices[1] == 1))
1540  base_index[0] = 1;
1541  else
1542  base_index[0] = 2;
1543 
1544  if ((indices[2] == 0) && (indices[3] == 0))
1545  base_index[1] = 0;
1546  else if ((indices[2] == 1) && (indices[3] == 1))
1547  base_index[1] = 1;
1548  else
1549  base_index[1] = 2;
1550 
1551  return data[base_index[0]][base_index[1]];
1552  }
1553 
1554  case 3:
1555  // each entry of the tensor can be
1556  // thought of as an entry in a
1557  // matrix that maps the rolled-out
1558  // rank-2 tensors into rolled-out
1559  // rank-2 tensors. this is the
1560  // format in which we store rank-4
1561  // tensors. determine which
1562  // position the present entry is
1563  // stored in
1564  {
1565  unsigned int base_index[2] ;
1566  if ((indices[0] == 0) && (indices[1] == 0))
1567  base_index[0] = 0;
1568  else if ((indices[0] == 1) && (indices[1] == 1))
1569  base_index[0] = 1;
1570  else if ((indices[0] == 2) && (indices[1] == 2))
1571  base_index[0] = 2;
1572  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1573  ((indices[0] == 1) && (indices[1] == 0)))
1574  base_index[0] = 3;
1575  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1576  ((indices[0] == 2) && (indices[1] == 0)))
1577  base_index[0] = 4;
1578  else
1579  {
1580  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1581  ((indices[0] == 2) && (indices[1] == 1)),
1582  ExcInternalError());
1583  base_index[0] = 5;
1584  }
1585 
1586  if ((indices[2] == 0) && (indices[3] == 0))
1587  base_index[1] = 0;
1588  else if ((indices[2] == 1) && (indices[3] == 1))
1589  base_index[1] = 1;
1590  else if ((indices[2] == 2) && (indices[3] == 2))
1591  base_index[1] = 2;
1592  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1593  ((indices[2] == 1) && (indices[3] == 0)))
1594  base_index[1] = 3;
1595  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1596  ((indices[2] == 2) && (indices[3] == 0)))
1597  base_index[1] = 4;
1598  else
1599  {
1600  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1601  ((indices[2] == 2) && (indices[3] == 1)),
1602  ExcInternalError());
1603  base_index[1] = 5;
1604  }
1605 
1606  return data[base_index[0]][base_index[1]];
1607  }
1608 
1609  default:
1610  Assert (false, ExcNotImplemented());
1611  }
1612 
1613  static Number dummy;
1614  return dummy;
1615  }
1616 
1617 } // end of namespace internal
1618 
1619 
1620 
1621 template <int rank, int dim, typename Number>
1622 inline
1623 Number &
1625 {
1626  for (unsigned int r=0; r<rank; ++r)
1627  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1628  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1629 }
1630 
1631 
1632 
1633 template <int rank, int dim, typename Number>
1634 inline
1635 Number
1637 (const TableIndices<rank> &indices) const
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 namespace internal
1647 {
1648  namespace SymmetricTensor
1649  {
1650  template <int rank>
1652  get_partially_filled_indices (const unsigned int row,
1653  const internal::int2type<2> &)
1654  {
1655  return TableIndices<rank> (row,
1657 
1658  }
1659 
1660 
1661  template <int rank>
1663  get_partially_filled_indices (const unsigned int row,
1664  const internal::int2type<4> &)
1665  {
1666  return TableIndices<rank> (row,
1670 
1671  }
1672  }
1673 }
1674 
1675 
1676 template <int rank, int dim, typename Number>
1677 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1678 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1679 {
1680  return
1681  internal::SymmetricTensorAccessors::
1682  Accessor<rank,dim,true,rank-1,Number> (*this,
1683  internal::SymmetricTensor::get_partially_filled_indices<rank> (row,
1685 }
1686 
1687 
1688 
1689 template <int rank, int dim, typename Number>
1690 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
1691 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
1692 {
1693  return
1694  internal::SymmetricTensorAccessors::
1695  Accessor<rank,dim,false,rank-1,Number> (*this,
1696  internal::SymmetricTensor::get_partially_filled_indices<rank> (row,
1698 }
1699 
1700 
1701 
1702 template <int rank, int dim, typename Number>
1703 inline
1704 Number
1706 {
1707  return operator()(indices);
1708 }
1709 
1710 
1711 
1712 template <int rank, int dim, typename Number>
1713 inline
1714 Number &
1716 {
1717  return operator()(indices);
1718 }
1719 
1720 
1721 
1722 
1723 namespace internal
1724 {
1725  namespace SymmetricTensor
1726  {
1727  template <int dim, typename Number>
1728  unsigned int
1729  entry_to_indices (const ::SymmetricTensor<2,dim,Number> &,
1730  const unsigned int index)
1731  {
1732  return index;
1733  }
1734 
1735 
1736  template <int dim, typename Number>
1738  entry_to_indices (const ::SymmetricTensor<4,dim,Number> &,
1739  const unsigned int index)
1740  {
1741  return
1744  }
1745 
1746  }
1747 }
1748 
1749 
1750 
1751 template <int rank, int dim, typename Number>
1752 inline
1753 Number
1754 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
1755 {
1756  AssertIndexRange (index, n_independent_components);
1757  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1758 }
1759 
1760 
1761 
1762 template <int rank, int dim, typename Number>
1763 inline
1764 Number &
1765 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
1766 {
1767  AssertIndexRange (index, n_independent_components);
1768  return data[internal::SymmetricTensor::entry_to_indices(*this, index)];
1769 }
1770 
1771 
1772 
1773 namespace internal
1774 {
1775  template <int dim, typename Number>
1776  inline
1779  {
1780  switch (dim)
1781  {
1782  case 1:
1783  return numbers::NumberTraits<Number>::abs(data[0]);
1784 
1785  case 2:
1786  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1789 
1790  case 3:
1791  return std::sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
1797 
1798  default:
1799  {
1800  typename numbers::NumberTraits<Number>::real_type return_value
1802 
1803  for (unsigned int d=0; d<dim; ++d)
1804  return_value += numbers::NumberTraits<Number>::abs_square(data[d]);
1805  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
1806  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
1807 
1808  return std::sqrt(return_value);
1809  }
1810  }
1811  }
1812 
1813 
1814 
1815  template <int dim, typename Number>
1816  inline
1819  {
1820  switch (dim)
1821  {
1822  case 1:
1823  return numbers::NumberTraits<Number>::abs (data[0][0]);
1824 
1825  default:
1826  {
1827  typename numbers::NumberTraits<Number>::real_type return_value
1829 
1830  const unsigned int n_independent_components = data.dimension;
1831 
1832  for (unsigned int i=0; i<dim; ++i)
1833  for (unsigned int j=0; j<dim; ++j)
1834  return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
1835  for (unsigned int i=0; i<dim; ++i)
1836  for (unsigned int j=dim; j<n_independent_components; ++j)
1837  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1838  for (unsigned int i=dim; i<n_independent_components; ++i)
1839  for (unsigned int j=0; j<dim; ++j)
1840  return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1841  for (unsigned int i=dim; i<n_independent_components; ++i)
1842  for (unsigned int j=dim; j<n_independent_components; ++j)
1843  return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
1844 
1845  return std::sqrt(return_value);
1846  }
1847  }
1848  }
1849 
1850 } // end of namespace internal
1851 
1852 
1853 
1854 template <int rank, int dim, typename Number>
1855 inline
1858 {
1859  return internal::compute_norm<dim,Number> (data);
1860 }
1861 
1862 
1863 
1864 namespace internal
1865 {
1866  namespace SymmetricTensor
1867  {
1868  namespace
1869  {
1870  // a function to do the unrolling from a set of indices to a
1871  // scalar index into the array in which we store the elements of
1872  // a symmetric tensor
1873  //
1874  // this function is for rank-2 tensors
1875  template <int dim>
1876  inline
1877  unsigned int
1878  component_to_unrolled_index
1879  (const TableIndices<2> &indices)
1880  {
1881  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
1882  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
1883 
1884  switch (dim)
1885  {
1886  case 1:
1887  {
1888  return 0;
1889  }
1890 
1891  case 2:
1892  {
1893  static const unsigned int table[2][2] = {{0, 2},
1894  {2, 1}
1895  };
1896  return table[indices[0]][indices[1]];
1897  }
1898 
1899  case 3:
1900  {
1901  static const unsigned int table[3][3] = {{0, 3, 4},
1902  {3, 1, 5},
1903  {4, 5, 2}
1904  };
1905  return table[indices[0]][indices[1]];
1906  }
1907 
1908  case 4:
1909  {
1910  static const unsigned int table[4][4] = {{0, 4, 5, 6},
1911  {4, 1, 7, 8},
1912  {5, 7, 2, 9},
1913  {6, 8, 9, 3}
1914  };
1915  return table[indices[0]][indices[1]];
1916  }
1917 
1918  default:
1919  // for the remainder, manually figure out the numbering
1920  {
1921  if (indices[0] == indices[1])
1922  return indices[0];
1923 
1924  TableIndices<2> sorted_indices (indices);
1925  sorted_indices.sort ();
1926 
1927  for (unsigned int d=0, c=0; d<dim; ++d)
1928  for (unsigned int e=d+1; e<dim; ++e, ++c)
1929  if ((sorted_indices[0]==d) && (sorted_indices[1]==e))
1930  return dim+c;
1931 
1932  // should never get here:
1933  Assert(false, ExcInternalError());
1934  return 0;
1935  }
1936  }
1937  }
1938 
1939  // a function to do the unrolling from a set of indices to a
1940  // scalar index into the array in which we store the elements of
1941  // a symmetric tensor
1942  //
1943  // this function is for tensors of ranks not already handled
1944  // above
1945  template <int dim, int rank>
1946  inline
1947  unsigned int
1948  component_to_unrolled_index
1949  (const TableIndices<rank> &indices)
1950  {
1951  (void)indices;
1952  Assert (false, ExcNotImplemented());
1954  }
1955  }
1956  }
1957 }
1958 
1959 
1960 template <int rank, int dim, typename Number>
1961 inline
1962 unsigned int
1964 (const TableIndices<rank> &indices)
1965 {
1966  return internal::SymmetricTensor::component_to_unrolled_index<dim> (indices);
1967 }
1968 
1969 
1970 
1971 namespace internal
1972 {
1973  namespace SymmetricTensor
1974  {
1975  namespace
1976  {
1977  // a function to do the inverse of the unrolling from a set of
1978  // indices to a scalar index into the array in which we store
1979  // the elements of a symmetric tensor. in other words, it goes
1980  // from the scalar index into the array to a set of indices of
1981  // the tensor
1982  //
1983  // this function is for rank-2 tensors
1984  template <int dim>
1985  inline
1987  unrolled_to_component_indices
1988  (const unsigned int i,
1989  const int2type<2> &)
1990  {
1993  switch (dim)
1994  {
1995  case 1:
1996  {
1997  return TableIndices<2>(0,0);
1998  }
1999 
2000  case 2:
2001  {
2002  const TableIndices<2> table[3] =
2003  {
2004  TableIndices<2> (0,0),
2005  TableIndices<2> (1,1),
2006  TableIndices<2> (0,1)
2007  };
2008  return table[i];
2009  }
2010 
2011  case 3:
2012  {
2013  const TableIndices<2> table[6] =
2014  {
2015  TableIndices<2> (0,0),
2016  TableIndices<2> (1,1),
2017  TableIndices<2> (2,2),
2018  TableIndices<2> (0,1),
2019  TableIndices<2> (0,2),
2020  TableIndices<2> (1,2)
2021  };
2022  return table[i];
2023  }
2024 
2025  default:
2026  if (i<dim)
2027  return TableIndices<2> (i,i);
2028 
2029  for (unsigned int d=0, c=0; d<dim; ++d)
2030  for (unsigned int e=d+1; e<dim; ++e, ++c)
2031  if (c==i)
2032  return TableIndices<2>(d,e);
2033 
2034  // should never get here:
2035  Assert(false, ExcInternalError());
2036  return TableIndices<2>(0, 0);
2037  }
2038  }
2039 
2040  // a function to do the inverse of the unrolling from a set of
2041  // indices to a scalar index into the array in which we store
2042  // the elements of a symmetric tensor. in other words, it goes
2043  // from the scalar index into the array to a set of indices of
2044  // the tensor
2045  //
2046  // this function is for tensors of a rank not already handled
2047  // above
2048  template <int dim, int rank>
2049  inline
2051  unrolled_to_component_indices
2052  (const unsigned int i,
2053  const int2type<rank> &)
2054  {
2055  (void)i;
2058  Assert (false, ExcNotImplemented());
2059  return TableIndices<rank>();
2060  }
2061 
2062  }
2063  }
2064 }
2065 
2066 template <int rank, int dim, typename Number>
2067 inline
2070 (const unsigned int i)
2071 {
2072  return
2073  internal::SymmetricTensor::unrolled_to_component_indices<dim> (i,
2075 }
2076 
2077 
2078 
2079 template <int rank, int dim, typename Number>
2080 template <class Archive>
2081 inline
2082 void
2083 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2084 {
2085  ar &data;
2086 }
2087 
2088 
2089 #endif // DOXYGEN
2090 
2091 /* ----------------- Non-member functions operating on tensors. ------------ */
2092 
2093 
2100 template <int rank, int dim, typename Number, typename OtherNumber>
2101 inline
2104  const Tensor<rank, dim, OtherNumber> &right)
2105 {
2106  return Tensor<rank, dim, Number>(left) + right;
2107 }
2108 
2109 
2116 template <int rank, int dim, typename Number, typename OtherNumber>
2117 inline
2121 {
2122  return left + Tensor<rank, dim, OtherNumber>(right);
2123 }
2124 
2125 
2132 template <int rank, int dim, typename Number, typename OtherNumber>
2133 inline
2136  const Tensor<rank, dim, OtherNumber> &right)
2137 {
2138  return Tensor<rank, dim, Number>(left) - right;
2139 }
2140 
2141 
2148 template <int rank, int dim, typename Number, typename OtherNumber>
2149 inline
2153 {
2154  return left - Tensor<rank, dim, OtherNumber>(right);
2155 }
2156 
2157 
2158 
2172 template <int dim, typename Number>
2173 inline
2175 {
2176  switch (dim)
2177  {
2178  case 1:
2179  return t.data[0];
2180  case 2:
2181  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2182  case 3:
2183  {
2184  // in analogy to general tensors, but
2185  // there's something to be simplified for
2186  // the present case
2187  const Number tmp = t.data[3]*t.data[4]*t.data[5];
2188  return ( tmp + tmp
2189  +t.data[0]*t.data[1]*t.data[2]
2190  -t.data[0]*t.data[5]*t.data[5]
2191  -t.data[1]*t.data[4]*t.data[4]
2192  -t.data[2]*t.data[3]*t.data[3]);
2193  }
2194  default:
2195  Assert (false, ExcNotImplemented());
2196  return 0;
2197  }
2198 }
2199 
2200 
2201 
2211 template <int dim, typename Number>
2212 inline
2214 {
2215  return determinant (t);
2216 }
2217 
2218 
2219 
2227 template <int dim, typename Number>
2229 {
2230  Number t = d.data[0];
2231  for (unsigned int i=1; i<dim; ++i)
2232  t += d.data[i];
2233  return t;
2234 }
2235 
2236 
2246 template <int dim, typename Number>
2247 inline
2249 {
2250  return trace (t);
2251 }
2252 
2253 
2266 template <typename Number>
2267 inline
2269 {
2270  return 0;
2271 }
2272 
2273 
2274 
2295 template <typename Number>
2296 inline
2298 {
2299  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2300 }
2301 
2302 
2303 
2313 template <typename Number>
2314 inline
2316 {
2317  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2318  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2319 }
2320 
2321 
2322 
2323 
2333 template <int rank, int dim, typename Number>
2334 inline
2337 {
2338  return t;
2339 }
2340 
2341 
2342 
2352 template <int dim, typename Number>
2353 inline
2356 {
2358 
2359  // subtract scaled trace from the diagonal
2360  const Number tr = trace(t) / dim;
2361  for (unsigned int i=0; i<dim; ++i)
2362  tmp.data[i] -= tr;
2363 
2364  return tmp;
2365 }
2366 
2367 
2368 
2376 template <int dim, typename Number>
2377 inline
2379 unit_symmetric_tensor ()
2380 {
2381  // create a default constructed matrix filled with
2382  // zeros, then set the diagonal elements to one
2384  switch (dim)
2385  {
2386  case 1:
2387  tmp.data[0] = 1;
2388  break;
2389  case 2:
2390  tmp.data[0] = tmp.data[1] = 1;
2391  break;
2392  case 3:
2393  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
2394  break;
2395  default:
2396  for (unsigned int d=0; d<dim; ++d)
2397  tmp.data[d] = 1;
2398  }
2399  return tmp;
2400 }
2401 
2402 
2403 
2412 template <int dim>
2413 inline
2415 unit_symmetric_tensor ()
2416 {
2417  return unit_symmetric_tensor<dim,double>();
2418 }
2419 
2420 
2421 
2436 template <int dim, typename Number>
2437 inline
2439 deviator_tensor ()
2440 {
2442 
2443  // fill the elements treating the diagonal
2444  for (unsigned int i=0; i<dim; ++i)
2445  for (unsigned int j=0; j<dim; ++j)
2446  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
2447 
2448  // then fill the ones that copy over the
2449  // non-diagonal elements. note that during
2450  // the double-contraction, we handle the
2451  // off-diagonal elements twice, so simply
2452  // copying requires a weight of 1/2
2453  for (unsigned int i=dim;
2454  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2455  ++i)
2456  tmp.data[i][i] = 0.5;
2457 
2458  return tmp;
2459 }
2460 
2461 
2462 
2477 template <int dim>
2478 inline
2480 deviator_tensor ()
2481 {
2482  return deviator_tensor<dim,double>();
2483 }
2484 
2485 
2486 
2509 template <int dim, typename Number>
2510 inline
2512 identity_tensor ()
2513 {
2515 
2516  // fill the elements treating the diagonal
2517  for (unsigned int i=0; i<dim; ++i)
2518  tmp.data[i][i] = 1;
2519 
2520  // then fill the ones that copy over the
2521  // non-diagonal elements. note that during
2522  // the double-contraction, we handle the
2523  // off-diagonal elements twice, so simply
2524  // copying requires a weight of 1/2
2525  for (unsigned int i=dim;
2526  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2527  ++i)
2528  tmp.data[i][i] = 0.5;
2529 
2530  return tmp;
2531 }
2532 
2533 
2534 
2556 template <int dim>
2557 inline
2559 identity_tensor ()
2560 {
2561  return identity_tensor<dim,double>();
2562 }
2563 
2564 
2565 
2576 template <int dim, typename Number>
2577 inline
2580 {
2581  // if desired, take over the
2582  // inversion of a 4x4 tensor
2583  // from the FullMatrix
2584  AssertThrow (false, ExcNotImplemented());
2585 
2587 }
2588 
2589 
2590 
2591 #ifndef DOXYGEN
2592 
2593 template <typename Number>
2594 inline
2596 invert (const SymmetricTensor<2,1,Number> &t)
2597 {
2599 
2600  tmp[0][0] = 1.0/t[0][0];
2601 
2602  return tmp;
2603 }
2604 
2605 
2606 
2607 template <typename Number>
2608 inline
2610 invert (const SymmetricTensor<2,2,Number> &t)
2611 {
2613 
2614  // Sympy result: ([
2615  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
2616  // [-t01/(t00*t11 - t01**2), t00/(t00*t11 - t01**2)] ])
2617  const TableIndices<2> idx_00 (0,0);
2618  const TableIndices<2> idx_01 (0,1);
2619  const TableIndices<2> idx_11 (1,1);
2620  const Number inv_det_t
2621  = 1.0/(t[idx_00]*t[idx_11]
2622  - t[idx_01]*t[idx_01]);
2623  tmp[idx_00] = t[idx_11];
2624  tmp[idx_01] = -t[idx_01];
2625  tmp[idx_11] = t[idx_00];
2626  tmp *= inv_det_t;
2627 
2628  return tmp;
2629 }
2630 
2631 
2632 
2633 template <typename Number>
2634 inline
2636 invert (const SymmetricTensor<2,3,Number> &t)
2637 {
2639 
2640  // Sympy result: ([
2641  // [ (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2642  // (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2643  // (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
2644  // [ (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2645  // (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2646  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
2647  // [ (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
2648  // (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
2649  // (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)] ])
2650  const TableIndices<2> idx_00 (0,0);
2651  const TableIndices<2> idx_01 (0,1);
2652  const TableIndices<2> idx_02 (0,2);
2653  const TableIndices<2> idx_11 (1,1);
2654  const TableIndices<2> idx_12 (1,2);
2655  const TableIndices<2> idx_22 (2,2);
2656  const Number inv_det_t
2657  = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
2658  - t[idx_00]*t[idx_12]*t[idx_12]
2659  - t[idx_01]*t[idx_01]*t[idx_22]
2660  + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
2661  - t[idx_02]*t[idx_02]*t[idx_11]);
2662  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
2663  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
2664  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
2665  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
2666  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
2667  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
2668  tmp *= inv_det_t;
2669 
2670  return tmp;
2671 }
2672 
2673 #endif /* DOXYGEN */
2674 
2675 
2676 
2690 template <int dim, typename Number>
2691 inline
2694 {
2696  switch (dim)
2697  {
2698  case 1:
2699  tmp.data[0][0] = 1./t.data[0][0];
2700  break;
2701  case 2:
2702 
2703  // inverting this tensor is a little more
2704  // complicated than necessary, since we
2705  // store the data of 't' as a 3x3 matrix
2706  // t.data, but the product between a rank-4
2707  // and a rank-2 tensor is really not the
2708  // product between this matrix and the
2709  // 3-vector of a rhs, but rather
2710  //
2711  // B.vec = t.data * mult * A.vec
2712  //
2713  // where mult is a 3x3 matrix with
2714  // entries [[1,0,0],[0,1,0],[0,0,2]] to
2715  // capture the fact that we need to add up
2716  // both the c_ij12*a_12 and the c_ij21*a_21
2717  // terms
2718  //
2719  // in addition, in this scheme, the
2720  // identity tensor has the matrix
2721  // representation mult^-1.
2722  //
2723  // the inverse of 't' therefore has the
2724  // matrix representation
2725  //
2726  // inv.data = mult^-1 * t.data^-1 * mult^-1
2727  //
2728  // in order to compute it, let's first
2729  // compute the inverse of t.data and put it
2730  // into tmp.data; at the end of the
2731  // function we then scale the last row and
2732  // column of the inverse by 1/2,
2733  // corresponding to the left and right
2734  // multiplication with mult^-1
2735  {
2736  const Number t4 = t.data[0][0]*t.data[1][1],
2737  t6 = t.data[0][0]*t.data[1][2],
2738  t8 = t.data[0][1]*t.data[1][0],
2739  t00 = t.data[0][2]*t.data[1][0],
2740  t01 = t.data[0][1]*t.data[2][0],
2741  t04 = t.data[0][2]*t.data[2][0],
2742  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
2743  t8*t.data[2][2]+t00*t.data[2][1]+
2744  t01*t.data[1][2]-t04*t.data[1][1]);
2745  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
2746  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
2747  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
2748  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
2749  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
2750  tmp.data[1][2] = -(t6-t00)*t07;
2751  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
2752  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
2753  tmp.data[2][2] = (t4-t8)*t07;
2754 
2755  // scale last row and column as mentioned
2756  // above
2757  tmp.data[2][0] /= 2;
2758  tmp.data[2][1] /= 2;
2759  tmp.data[0][2] /= 2;
2760  tmp.data[1][2] /= 2;
2761  tmp.data[2][2] /= 4;
2762  }
2763  break;
2764  default:
2765  Assert (false, ExcNotImplemented());
2766  }
2767  return tmp;
2768 }
2769 
2770 
2771 
2785 template <>
2787 invert (const SymmetricTensor<4,3,double> &t);
2788 // this function is implemented in the .cc file for double data types
2789 
2790 
2791 
2806 template <int dim, typename Number>
2807 inline
2811 {
2813 
2814  // fill only the elements really needed
2815  for (unsigned int i=0; i<dim; ++i)
2816  for (unsigned int j=i; j<dim; ++j)
2817  for (unsigned int k=0; k<dim; ++k)
2818  for (unsigned int l=k; l<dim; ++l)
2819  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
2820 
2821  return tmp;
2822 }
2823 
2824 
2825 
2834 template <int dim,typename Number>
2835 inline
2838 {
2839  Number array[(dim*dim+dim)/2];
2840  for (unsigned int d=0; d<dim; ++d)
2841  array[d] = t[d][d];
2842  for (unsigned int d=0, c=0; d<dim; ++d)
2843  for (unsigned int e=d+1; e<dim; ++e, ++c)
2844  array[dim+c] = (t[d][e]+t[e][d])*0.5;
2845  return SymmetricTensor<2,dim,Number>(array);
2846 }
2847 
2848 
2849 
2857 template <int rank, int dim, typename Number>
2858 inline
2861  const Number factor)
2862 {
2864  tt *= factor;
2865  return tt;
2866 }
2867 
2868 
2869 
2877 template <int rank, int dim, typename Number>
2878 inline
2880 operator * (const Number factor,
2882 {
2883  // simply forward to the other operator
2884  return t*factor;
2885 }
2886 
2887 
2888 
2914 template <int rank, int dim, typename Number, typename OtherNumber>
2915 inline
2918  const OtherNumber factor)
2919 {
2920  // form the product. we have to convert the two factors into the final
2921  // type via explicit casts because, for awkward reasons, the C++
2922  // standard committee saw it fit to not define an
2923  // operator*(float,std::complex<double>)
2924  // (as well as with switched arguments and double<->float).
2925  typedef typename ProductType<Number,OtherNumber>::type product_type;
2927  // we used to shorten the following by 'tt *= product_type(factor);'
2928  // which requires that a converting constructor
2929  // 'product_type::product_type(const OtherNumber) is defined.
2930  // however, a user-defined constructor is not allowed for aggregates,
2931  // e.g. VectorizedArray. therefore, we work around this issue using a
2932  // copy-assignment operator 'product_type::operator=(const OtherNumber)'
2933  // which we assume to be defined.
2934  product_type new_factor;
2935  new_factor = factor;
2936  tt *= new_factor;
2937  return tt;
2938 }
2939 
2940 
2941 
2950 template <int rank, int dim, typename Number, typename OtherNumber>
2951 inline
2953 operator * (const Number factor,
2955 {
2956  // simply forward to the other operator with switched arguments
2957  return (t*factor);
2958 }
2959 
2960 
2961 
2967 template <int rank, int dim, typename Number>
2968 inline
2971  const Number factor)
2972 {
2974  tt /= factor;
2975  return tt;
2976 }
2977 
2978 
2979 
2986 template <int rank, int dim>
2987 inline
2989 operator * (const SymmetricTensor<rank,dim> &t,
2990  const double factor)
2991 {
2993  tt *= factor;
2994  return tt;
2995 }
2996 
2997 
2998 
3005 template <int rank, int dim>
3006 inline
3008 operator * (const double factor,
3009  const SymmetricTensor<rank,dim> &t)
3010 {
3012  tt *= factor;
3013  return tt;
3014 }
3015 
3016 
3017 
3023 template <int rank, int dim>
3024 inline
3026 operator / (const SymmetricTensor<rank,dim> &t,
3027  const double factor)
3028 {
3030  tt /= factor;
3031  return tt;
3032 }
3033 
3043 template <int dim, typename Number>
3044 inline
3045 Number
3048 {
3049  return (t1*t2);
3050 }
3051 
3052 
3062 template <int dim, typename Number>
3063 inline
3064 Number
3066  const Tensor<2,dim,Number> &t2)
3067 {
3068  Number s = 0;
3069  for (unsigned int i=0; i<dim; ++i)
3070  for (unsigned int j=0; j<dim; ++j)
3071  s += t1[i][j] * t2[i][j];
3072  return s;
3073 }
3074 
3075 
3085 template <int dim, typename Number>
3086 inline
3087 Number
3090 {
3091  return scalar_product(t2, t1);
3092 }
3093 
3094 
3110 template <typename Number>
3111 inline
3112 void
3114  const SymmetricTensor<4,1,Number> &t,
3115  const SymmetricTensor<2,1,Number> &s)
3116 {
3117  tmp[0][0] = t[0][0][0][0] * s[0][0];
3118 }
3119 
3120 
3121 
3137 template <typename Number>
3138 inline
3139 void
3141  const SymmetricTensor<2,1,Number> &s,
3142  const SymmetricTensor<4,1,Number> &t)
3143 {
3144  tmp[0][0] = t[0][0][0][0] * s[0][0];
3145 }
3146 
3147 
3148 
3163 template <typename Number>
3164 inline
3165 void
3167  const SymmetricTensor<4,2,Number> &t,
3168  const SymmetricTensor<2,2,Number> &s)
3169 {
3170  const unsigned int dim = 2;
3171 
3172  for (unsigned int i=0; i<dim; ++i)
3173  for (unsigned int j=i; j<dim; ++j)
3174  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3175  t[i][j][1][1] * s[1][1] +
3176  2 * t[i][j][0][1] * s[0][1];
3177 }
3178 
3179 
3180 
3196 template <typename Number>
3197 inline
3198 void
3200  const SymmetricTensor<2,2,Number> &s,
3201  const SymmetricTensor<4,2,Number> &t)
3202 {
3203  const unsigned int dim = 2;
3204 
3205  for (unsigned int i=0; i<dim; ++i)
3206  for (unsigned int j=i; j<dim; ++j)
3207  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3208  s[1][1] * t[1][1][i][j] +
3209  2 * s[0][1] * t[0][1][i][j];
3210 }
3211 
3212 
3213 
3229 template <typename Number>
3230 inline
3231 void
3233  const SymmetricTensor<4,3,Number> &t,
3234  const SymmetricTensor<2,3,Number> &s)
3235 {
3236  const unsigned int dim = 3;
3237 
3238  for (unsigned int i=0; i<dim; ++i)
3239  for (unsigned int j=i; j<dim; ++j)
3240  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3241  t[i][j][1][1] * s[1][1] +
3242  t[i][j][2][2] * s[2][2] +
3243  2 * t[i][j][0][1] * s[0][1] +
3244  2 * t[i][j][0][2] * s[0][2] +
3245  2 * t[i][j][1][2] * s[1][2];
3246 }
3247 
3248 
3249 
3265 template <typename Number>
3266 inline
3267 void
3269  const SymmetricTensor<2,3,Number> &s,
3270  const SymmetricTensor<4,3,Number> &t)
3271 {
3272  const unsigned int dim = 3;
3273 
3274  for (unsigned int i=0; i<dim; ++i)
3275  for (unsigned int j=i; j<dim; ++j)
3276  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3277  s[1][1] * t[1][1][i][j] +
3278  s[2][2] * t[2][2][i][j] +
3279  2 * s[0][1] * t[0][1][i][j] +
3280  2 * s[0][2] * t[0][2][i][j] +
3281  2 * s[1][2] * t[1][2][i][j];
3282 }
3283 
3284 
3285 
3293 template <int dim, typename Number>
3295 operator * (const SymmetricTensor<2,dim,Number> &src1,
3296  const Tensor<1,dim,Number> &src2)
3297 {
3298  Tensor<1,dim,Number> dest;
3299  for (unsigned int i=0; i<dim; ++i)
3300  for (unsigned int j=0; j<dim; ++j)
3301  dest[i] += src1[i][j] * src2[j];
3302  return dest;
3303 }
3304 
3305 
3313 template <int dim, typename Number>
3315 operator * (const Tensor<1,dim,Number> &src1,
3316  const SymmetricTensor<2,dim,Number> &src2)
3317 {
3318  // this is easy for symmetric tensors:
3319  return src2 * src1;
3320 }
3321 
3322 
3332 template <int dim, typename Number>
3333 inline
3334 std::ostream &operator << (std::ostream &out,
3336 {
3337  //make out lives a bit simpler by outputing
3338  //the tensor through the operator for the
3339  //general Tensor class
3341 
3342  for (unsigned int i=0; i<dim; ++i)
3343  for (unsigned int j=0; j<dim; ++j)
3344  tt[i][j] = t[i][j];
3345 
3346  return out << tt;
3347 }
3348 
3349 
3350 
3360 template <int dim, typename Number>
3361 inline
3362 std::ostream &operator << (std::ostream &out,
3364 {
3365  //make out lives a bit simpler by outputing
3366  //the tensor through the operator for the
3367  //general Tensor class
3369 
3370  for (unsigned int i=0; i<dim; ++i)
3371  for (unsigned int j=0; j<dim; ++j)
3372  for (unsigned int k=0; k<dim; ++k)
3373  for (unsigned int l=0; l<dim; ++l)
3374  tt[i][j][k][l] = t[i][j][k][l];
3375 
3376  return out << tt;
3377 }
3378 
3379 
3380 DEAL_II_NAMESPACE_CLOSE
3381 
3382 #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:173
bool operator!=(const SymmetricTensor &) const
Tensor< rank, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank, dim, Number > &left, const Tensor< rank, dim, OtherNumber > &right)
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:1190
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:389
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:328
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()