Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
symmetric_tensor.h
1 // ---------------------------------------------------------------------
2 // @f$Id: symmetric_tensor.h 30315 2013-08-15 02:52:42Z bangerth @f$
3 //
4 // Copyright (C) 2005 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__symmetric_tensor_h
18 #define __deal2__symmetric_tensor_h
19 
20 
21 #include <deal.II/base/tensor.h>
22 #include <deal.II/base/table_indices.h>
23 
25 
26 template <int rank, int dim, typename Number=double> class SymmetricTensor;
27 
28 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
30 template <int dim, typename Number> SymmetricTensor<4,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 dim2, typename Number> Number
38 
39 template <int dim, typename Number> SymmetricTensor<2,dim,Number>
41 template <int dim, typename Number> Number
43 
44 
45 
46 namespace internal
47 {
53  namespace SymmetricTensorAccessors
54  {
66  inline
67  TableIndices<2> merge (const TableIndices<2> &previous_indices,
68  const unsigned int new_index,
69  const unsigned int position)
70  {
71  Assert (position < 2, ExcIndexRange (position, 0, 2));
72 
73  if (position == 0)
74  return TableIndices<2>(new_index);
75  else
76  return TableIndices<2>(previous_indices[0], new_index);
77  }
78 
79 
80 
92  inline
93  TableIndices<4> merge (const TableIndices<4> &previous_indices,
94  const unsigned int new_index,
95  const unsigned int position)
96  {
97  Assert (position < 4, ExcIndexRange (position, 0, 4));
98 
99  switch (position)
100  {
101  case 0:
102  return TableIndices<4>(new_index);
103  case 1:
104  return TableIndices<4>(previous_indices[0],
105  new_index);
106  case 2:
107  return TableIndices<4>(previous_indices[0],
108  previous_indices[1],
109  new_index);
110  case 3:
111  return TableIndices<4>(previous_indices[0],
112  previous_indices[1],
113  previous_indices[2],
114  new_index);
115  }
116  Assert (false, ExcInternalError());
117  return TableIndices<4>();
118  }
119 
120 
135  template <int rank1, int rank2, int dim, typename Number>
137  {
138  typedef ::SymmetricTensor<rank1+rank2-4,dim,Number> type;
139  };
140 
141 
156  template <int dim, typename Number>
157  struct double_contraction_result<2,2,dim,Number>
158  {
159  typedef Number type;
160  };
161 
162 
163 
183  template <int rank, int dim, typename Number>
184  struct StorageType;
185 
190  template <int dim, typename Number>
191  struct StorageType<2,dim,Number>
192  {
198  static const unsigned int
199  n_independent_components = (dim*dim + dim)/2;
200 
206  };
207 
208 
209 
214  template <int dim, typename Number>
215  struct StorageType<4,dim,Number>
216  {
224  static const unsigned int
225  n_rank2_components = (dim*dim + dim)/2;
226 
231  static const unsigned int
232  n_independent_components = (n_rank2_components *
234 
245  };
246 
247 
248 
255  template <int rank, int dim, bool constness, typename Number>
257 
266  template <int rank, int dim, typename Number>
267  struct AccessorTypes<rank,dim,true,Number>
268  {
269  typedef const ::SymmetricTensor<rank,dim,Number> tensor_type;
270 
271  typedef Number reference;
272  };
273 
283  template <int rank, int dim, typename Number>
284  struct AccessorTypes<rank,dim,false,Number>
285  {
286  typedef ::SymmetricTensor<rank,dim,Number> tensor_type;
287 
288  typedef Number &reference;
289  };
290 
291 
327  template <int rank, int dim, bool constness, int P, typename Number>
328  class Accessor
329  {
330  public:
336  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
337 
338  private:
379  Accessor (tensor_type &tensor,
380  const TableIndices<rank> &previous_indices);
381 
387  Accessor ();
388 
394  Accessor (const Accessor &a);
395 
396  public:
397 
401  Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
402 
403  private:
408  tensor_type &tensor;
409  const TableIndices<rank> previous_indices;
410 
411  // declare some other classes
412  // as friends. make sure to
413  // work around bugs in some
414  // compilers
415  template <int,int,typename> friend class SymmetricTensor;
416  template <int,int,bool,int,typename>
417  friend class Accessor;
418 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
419  friend class ::SymmetricTensor<rank,dim,Number>;
420  friend class Accessor<rank,dim,constness,P+1,Number>;
421 # endif
422  };
423 
424 
425 
436  template <int rank, int dim, bool constness, typename Number>
437  class Accessor<rank,dim,constness,1,Number>
438  {
439  public:
445  typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
446 
447  private:
493  Accessor (tensor_type &tensor,
494  const TableIndices<rank> &previous_indices);
495 
501  Accessor ();
502 
508  Accessor (const Accessor &a);
509 
510  public:
511 
515  reference operator [] (const unsigned int);
516 
517  private:
522  tensor_type &tensor;
523  const TableIndices<rank> previous_indices;
524 
525  // declare some other classes
526  // as friends. make sure to
527  // work around bugs in some
528  // compilers
529  template <int,int,typename> friend class SymmetricTensor;
530  template <int,int,bool,int,typename>
532 # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
533  friend class ::SymmetricTensor<rank,dim,Number>;
534  friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
535 # endif
536  };
537  }
538 }
539 
540 
541 
608 template <int rank, int dim, typename Number>
609 class SymmetricTensor
610 {
611 public:
629  static const unsigned int dimension = dim;
630 
639  static const unsigned int n_independent_components
642 
647  SymmetricTensor ();
648 
668 
696  SymmetricTensor (const Number (&array) [n_independent_components]);
697 
702 
715  SymmetricTensor &operator = (const Number d);
716 
723  operator Tensor<rank,dim,Number> () const;
724 
728  bool operator == (const SymmetricTensor &) const;
729 
733  bool operator != (const SymmetricTensor &) const;
734 
739 
744 
750  SymmetricTensor &operator *= (const Number factor);
751 
756  SymmetricTensor &operator /= (const Number factor);
757 
765 
773 
779 
833 
842 
847  Number &operator() (const TableIndices<rank> &indices);
848 
862  Number operator() (const TableIndices<rank> &indices) const;
863 
869  internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
870  operator [] (const unsigned int row) const;
871 
877  internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
878  operator [] (const unsigned int row);
879 
885  Number
886  operator [] (const TableIndices<rank> &indices) const;
887 
893  Number &
894  operator [] (const TableIndices<rank> &indices);
895 
905  Number
906  access_raw_entry (const unsigned int unrolled_index) const;
907 
917  Number &
918  access_raw_entry (const unsigned int unrolled_index);
919 
935  Number norm () const;
936 
949  static
950  unsigned int
952 
963  static
965  unrolled_to_component_indices (const unsigned int i);
966 
985  void clear ();
986 
993  static std::size_t memory_consumption ();
994 
999  template <class Archive>
1000  void serialize(Archive &ar, const unsigned int version);
1001 
1002 private:
1007  typedef
1010 
1015  typedef typename base_tensor_descriptor::base_tensor_type base_tensor_type;
1016 
1021  base_tensor_type data;
1022 
1026  template <int, int, typename> friend class SymmetricTensor;
1027 
1031  template <int dim2, typename Number2>
1032  friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
1033 
1034  template <int dim2, typename Number2>
1035  friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
1036 
1037  template <int dim2, typename Number2>
1039  deviator (const SymmetricTensor<2,dim2,Number2> &t);
1040 
1041  template <int dim2, typename Number2>
1043 
1044  template <int dim2, typename Number2>
1046 
1047  template <int dim2, typename Number2>
1049 
1050  template <int dim2, typename Number2>
1052 };
1053 
1054 
1055 
1056 // ------------------------- inline functions ------------------------
1057 
1058 #ifndef DOXYGEN
1059 
1060 namespace internal
1061 {
1062  namespace SymmetricTensorAccessors
1063  {
1064  template <int rank, int dim, bool constness, int P, typename Number>
1066  Accessor (tensor_type &tensor,
1067  const TableIndices<rank> &previous_indices)
1068  :
1069  tensor (tensor),
1070  previous_indices (previous_indices)
1071  {}
1072 
1073 
1074 
1075  template <int rank, int dim, bool constness, int P, typename Number>
1076  Accessor<rank,dim,constness,P-1,Number>
1077  Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
1078  {
1079  return Accessor<rank,dim,constness,P-1,Number> (tensor,
1080  merge (previous_indices, i, rank-P));
1081  }
1082 
1083 
1084 
1085  template <int rank, int dim, bool constness, typename Number>
1086  Accessor<rank,dim,constness,1,Number>::
1087  Accessor (tensor_type &tensor,
1088  const TableIndices<rank> &previous_indices)
1089  :
1090  tensor (tensor),
1091  previous_indices (previous_indices)
1092  {}
1093 
1094 
1095 
1096  template <int rank, int dim, bool constness, typename Number>
1097  typename Accessor<rank,dim,constness,1,Number>::reference
1098  Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
1099  {
1100  return tensor(merge (previous_indices, i, rank-1));
1101  }
1102 
1103 
1104  }
1105 }
1106 
1107 
1108 
1109 template <int rank, int dim, typename Number>
1110 inline
1112 {}
1113 
1114 
1115 
1116 template <int rank, int dim, typename Number>
1117 inline
1119 {
1120  Assert (rank == 2, ExcNotImplemented());
1121  switch (dim)
1122  {
1123  case 2:
1124  Assert (t[0][1] == t[1][0], ExcInternalError());
1125 
1126  data[0] = t[0][0];
1127  data[1] = t[1][1];
1128  data[2] = t[0][1];
1129 
1130  break;
1131  case 3:
1132  Assert (t[0][1] == t[1][0], ExcInternalError());
1133  Assert (t[0][2] == t[2][0], ExcInternalError());
1134  Assert (t[1][2] == t[2][1], ExcInternalError());
1135 
1136  data[0] = t[0][0];
1137  data[1] = t[1][1];
1138  data[2] = t[2][2];
1139  data[3] = t[0][1];
1140  data[4] = t[0][2];
1141  data[5] = t[1][2];
1142 
1143  break;
1144  default:
1145  Assert (false, ExcNotImplemented());
1146  }
1147 }
1148 
1149 
1150 
1151 template <int rank, int dim, typename Number>
1152 inline
1153 SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [n_independent_components])
1154  :
1155  data (array)
1156 {}
1157 
1158 
1159 
1160 template <int rank, int dim, typename Number>
1161 inline
1164 {
1165  data = t.data;
1166  return *this;
1167 }
1168 
1169 
1170 
1171 template <int rank, int dim, typename Number>
1172 inline
1175 {
1176  Assert (d==0, ExcMessage ("Only assignment with zero is allowed"));
1177  (void) d;
1178 
1179  data = 0;
1180 
1181  return *this;
1182 }
1183 
1184 
1185 
1186 // helper function to convert symmetric tensor
1187 // to generic tensor
1188 namespace internal
1189 {
1190  template <typename Number>
1191  inline
1193  conversion (const Tensor<1,1,Number> &data)
1194  {
1195  const Number t[1][1] = {{data[0]}};
1196  return Tensor<2,1,Number>(t);
1197  }
1198 
1199  template <typename Number>
1200  inline
1202  conversion (const Tensor<1,3,Number> &data)
1203  {
1204  const Number t[2][2] = {{data[0], data[2]},
1205  {data[2], data[1]}
1206  };
1207  return Tensor<2,2,Number>(t);
1208  }
1209 
1210  template <typename Number>
1211  inline
1213  conversion (const Tensor<1,6,Number> &data)
1214  {
1215  const Number t[3][3] = {{data[0], data[3], data[4]},
1216  {data[3], data[1], data[5]},
1217  {data[4], data[5], data[2]}
1218  };
1219  return Tensor<2,3,Number>(t);
1220  }
1221 }
1222 
1223 
1224 
1225 template <int rank, int dim, typename Number>
1226 inline
1228 operator Tensor<rank,dim,Number> () const
1229 {
1230  Assert (rank == 2, ExcNotImplemented());
1231  return internal::conversion(data);
1232 }
1233 
1234 
1235 
1236 template <int rank, int dim, typename Number>
1237 inline
1238 bool
1240 (const SymmetricTensor<rank,dim,Number> &t) const
1241 {
1242  return data == t.data;
1243 }
1244 
1245 
1246 
1247 template <int rank, int dim, typename Number>
1248 inline
1249 bool
1250 SymmetricTensor<rank,dim,Number>::operator !=
1251 (const SymmetricTensor<rank,dim,Number> &t) const
1252 {
1253  return data != t.data;
1254 }
1255 
1256 
1257 
1258 template <int rank, int dim, typename Number>
1259 inline
1261 SymmetricTensor<rank,dim,Number>::operator +=
1263 {
1264  data += t.data;
1265  return *this;
1266 }
1267 
1268 
1269 
1270 template <int rank, int dim, typename Number>
1271 inline
1273 SymmetricTensor<rank,dim,Number>::operator -=
1275 {
1276  data -= t.data;
1277  return *this;
1278 }
1279 
1280 
1281 
1282 template <int rank, int dim, typename Number>
1283 inline
1286 {
1287  data *= d;
1288  return *this;
1289 }
1290 
1291 
1292 
1293 template <int rank, int dim, typename Number>
1294 inline
1297 {
1298  data /= d;
1299  return *this;
1300 }
1301 
1302 
1303 
1304 template <int rank, int dim, typename Number>
1305 inline
1308 {
1309  SymmetricTensor tmp = *this;
1310  tmp.data += t.data;
1311  return tmp;
1312 }
1313 
1314 
1315 
1316 template <int rank, int dim, typename Number>
1317 inline
1320 {
1321  SymmetricTensor tmp = *this;
1322  tmp.data -= t.data;
1323  return tmp;
1324 }
1325 
1326 
1327 
1328 template <int rank, int dim, typename Number>
1329 inline
1332 {
1333  SymmetricTensor tmp = *this;
1334  tmp.data = -tmp.data;
1335  return tmp;
1336 }
1337 
1338 
1339 
1340 template <int rank, int dim, typename Number>
1341 inline
1342 void
1344 {
1345  data.clear ();
1346 }
1347 
1348 
1349 
1350 template <int rank, int dim, typename Number>
1351 inline
1352 std::size_t
1354 {
1355  return
1357 }
1358 
1359 
1360 
1361 namespace internal
1362 {
1363 
1364  template <int dim, typename Number>
1365  inline
1366  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
1367  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1369  {
1370  switch (dim)
1371  {
1372  case 1:
1373  return data[0] * sdata[0];
1374  case 2:
1375  return (data[0] * sdata[0] +
1376  data[1] * sdata[1] +
1377  2*data[2] * sdata[2]);
1378  case 3:
1379  return (data[0] * sdata[0] +
1380  data[1] * sdata[1] +
1381  data[2] * sdata[2] +
1382  2*data[3] * sdata[3] +
1383  2*data[4] * sdata[4] +
1384  2*data[5] * sdata[5]);
1385  default:
1386  Number sum = 0;
1387  for (unsigned int d=0; d<dim; ++d)
1388  sum += data[d] * sdata[d];
1389  for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
1390  sum += Number(2.) * data[d] * sdata[d];
1391  return sum;
1392  }
1393  }
1394 
1395 
1396 
1397  template <int dim, typename Number>
1398  inline
1399  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
1400  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1402  {
1404  switch (dim)
1405  {
1406  case 1:
1407  tmp[0] = data[0][0] * sdata[0];
1408  break;
1409  case 2:
1410  for (unsigned int i=0; i<3; ++i)
1411  tmp[i] = (data[i][0] * sdata[0] +
1412  data[i][1] * sdata[1] +
1413  2 * data[i][2] * sdata[2]);
1414  break;
1415  case 3:
1416  for (unsigned int i=0; i<6; ++i)
1417  tmp[i] = (data[i][0] * sdata[0] +
1418  data[i][1] * sdata[1] +
1419  data[i][2] * sdata[2] +
1420  2 * data[i][3] * sdata[3] +
1421  2 * data[i][4] * sdata[4] +
1422  2 * data[i][5] * sdata[5]);
1423  break;
1424  default:
1425  Assert (false, ExcNotImplemented());
1426  }
1427  return SymmetricTensor<2,dim,Number>(tmp);
1428  }
1429 
1430 
1431 
1432  template <int dim, typename Number>
1433  inline
1435  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
1437  {
1439  switch (dim)
1440  {
1441  case 1:
1442  tmp[0] = data[0] * sdata[0][0];
1443  break;
1444  case 2:
1445  for (unsigned int i=0; i<3; ++i)
1446  tmp[i] = (data[0] * sdata[0][i] +
1447  data[1] * sdata[1][i] +
1448  2 * data[2] * sdata[2][i]);
1449  break;
1450  case 3:
1451  for (unsigned int i=0; i<6; ++i)
1452  tmp[i] = (data[0] * sdata[0][i] +
1453  data[1] * sdata[1][i] +
1454  data[2] * sdata[2][i] +
1455  2 * data[3] * sdata[3][i] +
1456  2 * data[4] * sdata[4][i] +
1457  2 * data[5] * sdata[5][i]);
1458  break;
1459  default:
1460  Assert (false, ExcNotImplemented());
1461  }
1462  return tmp;
1463  }
1464 
1465 
1466 
1467  template <int dim, typename Number>
1468  inline
1470  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
1472  {
1474  switch (dim)
1475  {
1476  case 1:
1477  tmp[0][0] = data[0][0] * sdata[0][0];
1478  break;
1479  case 2:
1480  for (unsigned int i=0; i<3; ++i)
1481  for (unsigned int j=0; j<3; ++j)
1482  tmp[i][j] = (data[i][0] * sdata[0][j] +
1483  data[i][1] * sdata[1][j] +
1484  2*data[i][2] * sdata[2][j]);
1485  break;
1486  case 3:
1487  for (unsigned int i=0; i<6; ++i)
1488  for (unsigned int j=0; j<6; ++j)
1489  tmp[i][j] = (data[i][0] * sdata[0][j] +
1490  data[i][1] * sdata[1][j] +
1491  data[i][2] * sdata[2][j] +
1492  2*data[i][3] * sdata[3][j] +
1493  2*data[i][4] * sdata[4][j] +
1494  2*data[i][5] * sdata[5][j]);
1495  break;
1496  default:
1497  Assert (false, ExcNotImplemented());
1498  }
1499  return tmp;
1500  }
1501 
1502 } // end of namespace internal
1503 
1504 
1505 
1506 template <int rank, int dim, typename Number>
1507 inline
1510 {
1511  // need to have two different function calls
1512  // because a scalar and rank-2 tensor are not
1513  // the same data type (see internal function
1514  // above)
1515  return internal::perform_double_contraction<dim,Number> (data, s.data);
1516 }
1517 
1518 
1519 
1520 template <int rank, int dim, typename Number>
1521 inline
1524 {
1527  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
1528  return tmp;
1529 }
1530 
1531 
1532 
1533 // internal namespace to switch between the
1534 // access of different tensors. There used to
1535 // be explicit instantiations before for
1536 // different ranks and dimensions, but since
1537 // we now allow for templates on the data
1538 // type, and since we cannot partially
1539 // specialize the implementation, this got
1540 // into a separate namespace
1541 namespace internal
1542 {
1543  template <int dim, typename Number>
1544  inline
1545  Number &
1546  symmetric_tensor_access (const TableIndices<2> &indices,
1548  {
1549  switch (dim)
1550  {
1551  case 1:
1552  return data[0];
1553 
1554  case 2:
1555  // first treat the main diagonal
1556  // elements, which are stored
1557  // consecutively at the beginning
1558  if (indices[0] == indices[1])
1559  return data[indices[0]];
1560 
1561  // the rest is messier and requires a few
1562  // switches. at least for the 2x2 case it
1563  // is reasonably simple
1564  Assert (((indices[0]==1) && (indices[1]==0)) ||
1565  ((indices[0]==0) && (indices[1]==1)),
1566  ExcInternalError());
1567  return data[2];
1568 
1569  case 3:
1570  // first treat the main diagonal
1571  // elements, which are stored
1572  // consecutively at the beginning
1573  if (indices[0] == indices[1])
1574  return data[indices[0]];
1575 
1576  // the rest is messier and requires a few
1577  // switches, but simpler if we just sort
1578  // our indices
1579  {
1580  TableIndices<2> sorted_indices (indices);
1581  sorted_indices.sort ();
1582 
1583  if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
1584  return data[3];
1585  else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
1586  return data[4];
1587  else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
1588  return data[5];
1589  else
1590  Assert (false, ExcInternalError());
1591  }
1592  }
1593 
1594  static Number dummy_but_referenceable = Number();
1595  return dummy_but_referenceable;
1596  }
1597 
1598 
1599 
1600  template <int dim, typename Number>
1601  inline
1602  Number
1603  symmetric_tensor_access (const TableIndices<2> &indices,
1605  {
1606  switch (dim)
1607  {
1608  case 1:
1609  return data[0];
1610 
1611  case 2:
1612  // first treat the main diagonal
1613  // elements, which are stored
1614  // consecutively at the beginning
1615  if (indices[0] == indices[1])
1616  return data[indices[0]];
1617 
1618  // the rest is messier and requires a few
1619  // switches. at least for the 2x2 case it
1620  // is reasonably simple
1621  Assert (((indices[0]==1) && (indices[1]==0)) ||
1622  ((indices[0]==0) && (indices[1]==1)),
1623  ExcInternalError());
1624  return data[2];
1625 
1626  case 3:
1627  // first treat the main diagonal
1628  // elements, which are stored
1629  // consecutively at the beginning
1630  if (indices[0] == indices[1])
1631  return data[indices[0]];
1632 
1633  // the rest is messier and requires a few
1634  // switches, but simpler if we just sort
1635  // our indices
1636  {
1637  TableIndices<2> sorted_indices (indices);
1638  sorted_indices.sort ();
1639 
1640  if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
1641  return data[3];
1642  else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
1643  return data[4];
1644  else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
1645  return data[5];
1646  else
1647  Assert (false, ExcInternalError());
1648  }
1649  }
1650 
1651  static Number dummy_but_referenceable = 0;
1652  return dummy_but_referenceable;
1653  }
1654 
1655 
1656 
1657  template <int dim, typename Number>
1658  inline
1659  Number &
1660  symmetric_tensor_access (const TableIndices<4> &indices,
1662  {
1663  switch (dim)
1664  {
1665  case 1:
1666  return data[0][0];
1667 
1668  case 2:
1669  // each entry of the tensor can be
1670  // thought of as an entry in a
1671  // matrix that maps the rolled-out
1672  // rank-2 tensors into rolled-out
1673  // rank-2 tensors. this is the
1674  // format in which we store rank-4
1675  // tensors. determine which
1676  // position the present entry is
1677  // stored in
1678  {
1679  unsigned int base_index[2] ;
1680  if ((indices[0] == 0) && (indices[1] == 0))
1681  base_index[0] = 0;
1682  else if ((indices[0] == 1) && (indices[1] == 1))
1683  base_index[0] = 1;
1684  else
1685  base_index[0] = 2;
1686 
1687  if ((indices[2] == 0) && (indices[3] == 0))
1688  base_index[1] = 0;
1689  else if ((indices[2] == 1) && (indices[3] == 1))
1690  base_index[1] = 1;
1691  else
1692  base_index[1] = 2;
1693 
1694  return data[base_index[0]][base_index[1]];
1695  }
1696 
1697  case 3:
1698  // each entry of the tensor can be
1699  // thought of as an entry in a
1700  // matrix that maps the rolled-out
1701  // rank-2 tensors into rolled-out
1702  // rank-2 tensors. this is the
1703  // format in which we store rank-4
1704  // tensors. determine which
1705  // position the present entry is
1706  // stored in
1707  {
1708  unsigned int base_index[2] ;
1709  if ((indices[0] == 0) && (indices[1] == 0))
1710  base_index[0] = 0;
1711  else if ((indices[0] == 1) && (indices[1] == 1))
1712  base_index[0] = 1;
1713  else if ((indices[0] == 2) && (indices[1] == 2))
1714  base_index[0] = 2;
1715  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1716  ((indices[0] == 1) && (indices[1] == 0)))
1717  base_index[0] = 3;
1718  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1719  ((indices[0] == 2) && (indices[1] == 0)))
1720  base_index[0] = 4;
1721  else
1722  {
1723  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1724  ((indices[0] == 2) && (indices[1] == 1)),
1725  ExcInternalError());
1726  base_index[0] = 5;
1727  }
1728 
1729  if ((indices[2] == 0) && (indices[3] == 0))
1730  base_index[1] = 0;
1731  else if ((indices[2] == 1) && (indices[3] == 1))
1732  base_index[1] = 1;
1733  else if ((indices[2] == 2) && (indices[3] == 2))
1734  base_index[1] = 2;
1735  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1736  ((indices[2] == 1) && (indices[3] == 0)))
1737  base_index[1] = 3;
1738  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1739  ((indices[2] == 2) && (indices[3] == 0)))
1740  base_index[1] = 4;
1741  else
1742  {
1743  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1744  ((indices[2] == 2) && (indices[3] == 1)),
1745  ExcInternalError());
1746  base_index[1] = 5;
1747  }
1748 
1749  return data[base_index[0]][base_index[1]];
1750  }
1751 
1752  default:
1753  Assert (false, ExcNotImplemented());
1754  }
1755 
1756  static Number dummy;
1757  return dummy;
1758  }
1759 
1760 
1761  template <int dim, typename Number>
1762  inline
1763  Number
1764  symmetric_tensor_access (const TableIndices<4> &indices,
1766  {
1767  switch (dim)
1768  {
1769  case 1:
1770  return data[0][0];
1771 
1772  case 2:
1773  // each entry of the tensor can be
1774  // thought of as an entry in a
1775  // matrix that maps the rolled-out
1776  // rank-2 tensors into rolled-out
1777  // rank-2 tensors. this is the
1778  // format in which we store rank-4
1779  // tensors. determine which
1780  // position the present entry is
1781  // stored in
1782  {
1783  unsigned int base_index[2] ;
1784  if ((indices[0] == 0) && (indices[1] == 0))
1785  base_index[0] = 0;
1786  else if ((indices[0] == 1) && (indices[1] == 1))
1787  base_index[0] = 1;
1788  else
1789  base_index[0] = 2;
1790 
1791  if ((indices[2] == 0) && (indices[3] == 0))
1792  base_index[1] = 0;
1793  else if ((indices[2] == 1) && (indices[3] == 1))
1794  base_index[1] = 1;
1795  else
1796  base_index[1] = 2;
1797 
1798  return data[base_index[0]][base_index[1]];
1799  }
1800 
1801  case 3:
1802  // each entry of the tensor can be
1803  // thought of as an entry in a
1804  // matrix that maps the rolled-out
1805  // rank-2 tensors into rolled-out
1806  // rank-2 tensors. this is the
1807  // format in which we store rank-4
1808  // tensors. determine which
1809  // position the present entry is
1810  // stored in
1811  {
1812  unsigned int base_index[2] ;
1813  if ((indices[0] == 0) && (indices[1] == 0))
1814  base_index[0] = 0;
1815  else if ((indices[0] == 1) && (indices[1] == 1))
1816  base_index[0] = 1;
1817  else if ((indices[0] == 2) && (indices[1] == 2))
1818  base_index[0] = 2;
1819  else if (((indices[0] == 0) && (indices[1] == 1)) ||
1820  ((indices[0] == 1) && (indices[1] == 0)))
1821  base_index[0] = 3;
1822  else if (((indices[0] == 0) && (indices[1] == 2)) ||
1823  ((indices[0] == 2) && (indices[1] == 0)))
1824  base_index[0] = 4;
1825  else
1826  {
1827  Assert (((indices[0] == 1) && (indices[1] == 2)) ||
1828  ((indices[0] == 2) && (indices[1] == 1)),
1829  ExcInternalError());
1830  base_index[0] = 5;
1831  }
1832 
1833  if ((indices[2] == 0) && (indices[3] == 0))
1834  base_index[1] = 0;
1835  else if ((indices[2] == 1) && (indices[3] == 1))
1836  base_index[1] = 1;
1837  else if ((indices[2] == 2) && (indices[3] == 2))
1838  base_index[1] = 2;
1839  else if (((indices[2] == 0) && (indices[3] == 1)) ||
1840  ((indices[2] == 1) && (indices[3] == 0)))
1841  base_index[1] = 3;
1842  else if (((indices[2] == 0) && (indices[3] == 2)) ||
1843  ((indices[2] == 2) && (indices[3] == 0)))
1844  base_index[1] = 4;
1845  else
1846  {
1847  Assert (((indices[2] == 1) && (indices[3] == 2)) ||
1848  ((indices[2] == 2) && (indices[3] == 1)),
1849  ExcInternalError());
1850  base_index[1] = 5;
1851  }
1852 
1853  return data[base_index[0]][base_index[1]];
1854  }
1855 
1856  default:
1857  Assert (false, ExcNotImplemented());
1858  }
1859 
1860  static Number dummy;
1861  return dummy;
1862  }
1863 
1864 } // end of namespace internal
1865 
1866 
1867 
1868 template <int rank, int dim, typename Number>
1869 inline
1870 Number &
1872 {
1873  for (unsigned int r=0; r<rank; ++r)
1874  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1875  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1876 }
1877 
1878 
1879 
1880 template <int rank, int dim, typename Number>
1881 inline
1882 Number
1884 (const TableIndices<rank> &indices) const
1885 {
1886  for (unsigned int r=0; r<rank; ++r)
1887  Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
1888  return internal::symmetric_tensor_access<dim,Number> (indices, data);
1889 }
1890 
1891 
1892 
1893 template <int rank, int dim, typename Number>
1894 internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
1895 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
1896 {
1897  return
1899  Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
1900 }
1901 
1902 
1903 
1904 template <int rank, int dim, typename Number>
1905 internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
1906 SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
1907 {
1908  return
1910  Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
1911 }
1912 
1913 
1914 
1915 template <int rank, int dim, typename Number>
1916 inline
1917 Number
1919 {
1920  return data[component_to_unrolled_index(indices)];
1921 }
1922 
1923 
1924 
1925 template <int rank, int dim, typename Number>
1926 inline
1927 Number &
1929 {
1930  return data[component_to_unrolled_index(indices)];
1931 }
1932 
1933 
1934 
1935 template <int rank, int dim, typename Number>
1936 inline
1937 Number
1938 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index) const
1939 {
1940  AssertIndexRange (index, data.dimension);
1941  return data[index];
1942 }
1943 
1944 
1945 
1946 template <int rank, int dim, typename Number>
1947 inline
1948 Number &
1949 SymmetricTensor<rank,dim,Number>::access_raw_entry (const unsigned int index)
1950 {
1951  AssertIndexRange (index, data.dimension);
1952  return data[index];
1953 }
1954 
1955 
1956 
1957 namespace internal
1958 {
1959  template <int dim, typename Number>
1960  inline
1961  Number
1963  {
1964  Number return_value;
1965  switch (dim)
1966  {
1967  case 1:
1968  return_value = std::fabs(data[0]);
1969  break;
1970  case 2:
1971  return_value = std::sqrt(data[0]*data[0] + data[1]*data[1] +
1972  2*data[2]*data[2]);
1973  break;
1974  case 3:
1975  return_value = std::sqrt(data[0]*data[0] + data[1]*data[1] +
1976  data[2]*data[2] + 2*data[3]*data[3] +
1977  2*data[4]*data[4] + 2*data[5]*data[5]);
1978  break;
1979  default:
1980  return_value = 0;
1981  for (unsigned int d=0; d<dim; ++d)
1982  return_value += data[d] * data[d];
1983  for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
1984  return_value += 2 * data[d] * data[d];
1985  return_value = std::sqrt(return_value);
1986  }
1987  return return_value;
1988  }
1989 
1990 
1991 
1992  template <int dim, typename Number>
1993  inline
1994  Number
1996  {
1997  Number return_value;
1998  const unsigned int n_independent_components = data.dimension;
1999 
2000  switch (dim)
2001  {
2002  case 1:
2003  return_value = std::fabs (data[0][0]);
2004  break;
2005  default:
2006  return_value = 0;
2007  for (unsigned int i=0; i<dim; ++i)
2008  for (unsigned int j=0; j<dim; ++j)
2009  return_value += data[i][j] * data[i][j];
2010  for (unsigned int i=0; i<dim; ++i)
2011  for (unsigned int j=dim; j<n_independent_components; ++j)
2012  return_value += 2 * data[i][j] * data[i][j];
2013  for (unsigned int i=dim; i<n_independent_components; ++i)
2014  for (unsigned int j=0; j<dim; ++j)
2015  return_value += 2 * data[i][j] * data[i][j];
2016  for (unsigned int i=dim; i<n_independent_components; ++i)
2017  for (unsigned int j=dim; j<n_independent_components; ++j)
2018  return_value += 4 * data[i][j] * data[i][j];
2019  return_value = std::sqrt(return_value);
2020  }
2021 
2022  return return_value;
2023  }
2024 
2025 } // end of namespace internal
2026 
2027 
2028 
2029 template <int rank, int dim, typename Number>
2030 inline
2031 Number
2033 {
2034  return internal::compute_norm<dim,Number> (data);
2035 }
2036 
2037 
2038 
2039 template <int rank, int dim, typename Number>
2040 inline
2041 unsigned int
2043 (const TableIndices<rank> &indices)
2044 {
2045  Assert (rank == 2, ExcNotImplemented());
2046  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
2047  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
2048 
2049  switch (dim)
2050  {
2051  case 1:
2052  return 0;
2053  case 2:
2054  {
2055  static const unsigned int table[2][2] = {{0, 2},
2056  {2, 1}
2057  };
2058  return table[indices[0]][indices[1]];
2059  }
2060  case 3:
2061  {
2062  static const unsigned int table[3][3] = {{0, 3, 4},
2063  {3, 1, 5},
2064  {4, 5, 2}
2065  };
2066  return table[indices[0]][indices[1]];
2067  }
2068  case 4:
2069  {
2070  static const unsigned int table[4][4] = {{0, 4, 5, 6},
2071  {4, 1, 7, 8},
2072  {5, 7, 2, 9},
2073  {6, 8, 9, 3}
2074  };
2075  return table[indices[0]][indices[1]];
2076  }
2077  default:
2078  Assert (false, ExcNotImplemented());
2079  return 0;
2080  }
2081 }
2082 
2083 
2084 
2085 template <int rank, int dim, typename Number>
2086 inline
2089 (const unsigned int i)
2090 {
2091  Assert (rank == 2, ExcNotImplemented());
2092  Assert (i < n_independent_components, ExcIndexRange(i, 0, n_independent_components));
2093  switch (dim)
2094  {
2095  case 1:
2096  return TableIndices<2>(0,0);
2097  case 2:
2098  {
2099  static const TableIndices<2> table[3] =
2100  {
2101  TableIndices<2> (0,0),
2102  TableIndices<2> (1,1),
2103  TableIndices<2> (0,1)
2104  };
2105  return table[i];
2106  }
2107  case 3:
2108  {
2109  static const TableIndices<2> table[6] =
2110  {
2111  TableIndices<2> (0,0),
2112  TableIndices<2> (1,1),
2113  TableIndices<2> (2,2),
2114  TableIndices<2> (0,1),
2115  TableIndices<2> (0,2),
2116  TableIndices<2> (1,2)
2117  };
2118  return table[i];
2119  }
2120  default:
2121  Assert (false, ExcNotImplemented());
2122  return TableIndices<2>(0,0);
2123  }
2124 }
2125 
2126 
2127 
2128 template <int rank, int dim, typename Number>
2129 template <class Archive>
2130 inline
2131 void
2132 SymmetricTensor<rank,dim,Number>::serialize(Archive &ar, const unsigned int)
2133 {
2134  ar &data;
2135 }
2136 
2137 
2138 #endif // DOXYGEN
2139 
2140 /* ----------------- Non-member functions operating on tensors. ------------ */
2141 
2155 template <int dim, typename Number>
2156 inline
2158 {
2159  switch (dim)
2160  {
2161  case 1:
2162  return t.data[0];
2163  case 2:
2164  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
2165  case 3:
2166  // in analogy to general tensors, but
2167  // there's something to be simplified for
2168  // the present case
2169  return ( t.data[0]*t.data[1]*t.data[2]
2170  -t.data[0]*t.data[5]*t.data[5]
2171  -t.data[1]*t.data[4]*t.data[4]
2172  -t.data[2]*t.data[3]*t.data[3]
2173  +2*t.data[3]*t.data[4]*t.data[5] );
2174  default:
2175  Assert (false, ExcNotImplemented());
2176  return 0;
2177  }
2178 }
2179 
2180 
2181 
2191 template <int dim, typename Number>
2192 inline
2194 {
2195  return determinant (t);
2196 }
2197 
2198 
2199 
2208 template <int dim, typename Number>
2210 {
2211  Number t=0;
2212  for (unsigned int i=0; i<dim; ++i)
2213  t += d.data[i];
2214  return t;
2215 }
2216 
2217 
2227 template <int dim, typename Number>
2228 inline
2230 {
2231  return trace (t);
2232 }
2233 
2234 
2242 template <typename Number>
2243 inline
2245 {
2246  return 0;
2247 }
2248 
2249 
2250 
2258 template <typename Number>
2259 inline
2261 {
2262  return t[0][0]*t[1][1] - t[0][1]*t[0][1];
2263 }
2264 
2265 
2266 
2274 template <typename Number>
2275 inline
2277 {
2278  return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
2279  - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
2280 }
2281 
2282 
2283 
2284 
2294 template <int rank, int dim, typename Number>
2295 inline
2298 {
2299  return t;
2300 }
2301 
2302 
2303 
2313 template <int dim, typename Number>
2314 inline
2317 {
2319 
2320  // subtract scaled trace from the diagonal
2321  const Number tr = trace(t) / dim;
2322  for (unsigned int i=0; i<dim; ++i)
2323  tmp.data[i] -= tr;
2324 
2325  return tmp;
2326 }
2327 
2328 
2329 
2336 template <int dim, typename Number>
2337 inline
2339 unit_symmetric_tensor ()
2340 {
2342  switch (dim)
2343  {
2344  case 1:
2345  tmp.data[0] = 1;
2346  break;
2347  case 2:
2348  tmp.data[0] = tmp.data[1] = 1;
2349  break;
2350  case 3:
2351  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
2352  break;
2353  default:
2354  for (unsigned int d=0; d<dim; ++d)
2355  tmp.data[d] = 1;
2356  }
2357  return tmp;
2358 }
2359 
2360 
2361 
2368 template <int dim>
2369 inline
2371 unit_symmetric_tensor ()
2372 {
2373  return unit_symmetric_tensor<dim,double>();
2374 }
2375 
2376 
2377 
2392 template <int dim, typename Number>
2393 inline
2395 deviator_tensor ()
2396 {
2398 
2399  // fill the elements treating the diagonal
2400  for (unsigned int i=0; i<dim; ++i)
2401  for (unsigned int j=0; j<dim; ++j)
2402  tmp.data[i][j] = (i==j ? 1 : 0) - 1./dim;
2403 
2404  // then fill the ones that copy over the
2405  // non-diagonal elements. note that during
2406  // the double-contraction, we handle the
2407  // off-diagonal elements twice, so simply
2408  // copying requires a weight of 1/2
2409  for (unsigned int i=dim;
2410  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2411  ++i)
2412  tmp.data[i][i] = 0.5;
2413 
2414  return tmp;
2415 }
2416 
2417 
2418 
2433 template <int dim>
2434 inline
2436 deviator_tensor ()
2437 {
2438  return deviator_tensor<dim,double>();
2439 }
2440 
2441 
2442 
2465 template <int dim, typename Number>
2466 inline
2468 identity_tensor ()
2469 {
2471 
2472  // fill the elements treating the diagonal
2473  for (unsigned int i=0; i<dim; ++i)
2474  tmp.data[i][i] = 1;
2475 
2476  // then fill the ones that copy over the
2477  // non-diagonal elements. note that during
2478  // the double-contraction, we handle the
2479  // off-diagonal elements twice, so simply
2480  // copying requires a weight of 1/2
2481  for (unsigned int i=dim;
2482  i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
2483  ++i)
2484  tmp.data[i][i] = 0.5;
2485 
2486  return tmp;
2487 }
2488 
2489 
2490 
2512 template <int dim>
2513 inline
2515 identity_tensor ()
2516 {
2517  return identity_tensor<dim,double>();
2518 }
2519 
2520 
2521 
2535 template <int dim, typename Number>
2536 inline
2539 {
2541  switch (dim)
2542  {
2543  case 1:
2544  tmp.data[0][0] = 1./t.data[0][0];
2545  break;
2546  case 2:
2547 
2548  // inverting this tensor is a little more
2549  // complicated than necessary, since we
2550  // store the data of 't' as a 3x3 matrix
2551  // t.data, but the product between a rank-4
2552  // and a rank-2 tensor is really not the
2553  // product between this matrix and the
2554  // 3-vector of a rhs, but rather
2555  //
2556  // B.vec = t.data * mult * A.vec
2557  //
2558  // where mult is a 3x3 matrix with
2559  // entries [[1,0,0],[0,1,0],[0,0,2]] to
2560  // capture the fact that we need to add up
2561  // both the c_ij12*a_12 and the c_ij21*a_21
2562  // terms
2563  //
2564  // in addition, in this scheme, the
2565  // identity tensor has the matrix
2566  // representation mult^-1.
2567  //
2568  // the inverse of 't' therefore has the
2569  // matrix representation
2570  //
2571  // inv.data = mult^-1 * t.data^-1 * mult^-1
2572  //
2573  // in order to compute it, let's first
2574  // compute the inverse of t.data and put it
2575  // into tmp.data; at the end of the
2576  // function we then scale the last row and
2577  // column of the inverse by 1/2,
2578  // corresponding to the left and right
2579  // multiplication with mult^-1
2580  {
2581  const Number t4 = t.data[0][0]*t.data[1][1],
2582  t6 = t.data[0][0]*t.data[1][2],
2583  t8 = t.data[0][1]*t.data[1][0],
2584  t00 = t.data[0][2]*t.data[1][0],
2585  t01 = t.data[0][1]*t.data[2][0],
2586  t04 = t.data[0][2]*t.data[2][0],
2587  t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
2588  t8*t.data[2][2]+t00*t.data[2][1]+
2589  t01*t.data[1][2]-t04*t.data[1][1]);
2590  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
2591  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
2592  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
2593  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
2594  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
2595  tmp.data[1][2] = -(t6-t00)*t07;
2596  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
2597  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
2598  tmp.data[2][2] = (t4-t8)*t07;
2599 
2600  // scale last row and column as mentioned
2601  // above
2602  tmp.data[2][0] /= 2;
2603  tmp.data[2][1] /= 2;
2604  tmp.data[0][2] /= 2;
2605  tmp.data[1][2] /= 2;
2606  tmp.data[2][2] /= 4;
2607  }
2608  break;
2609  default:
2610  Assert (false, ExcNotImplemented());
2611  }
2612  return tmp;
2613 }
2614 
2615 
2616 
2630 template <>
2633 // this function is implemented in the .cc file for double data types
2634 
2635 
2636 
2651 template <int dim, typename Number>
2652 inline
2656 {
2658 
2659  // fill only the elements really needed
2660  for (unsigned int i=0; i<dim; ++i)
2661  for (unsigned int j=i; j<dim; ++j)
2662  for (unsigned int k=0; k<dim; ++k)
2663  for (unsigned int l=k; l<dim; ++l)
2664  tmp[i][j][k][l] = t1[i][j] * t2[k][l];
2665 
2666  return tmp;
2667 }
2668 
2669 
2670 
2679 template <typename Number>
2680 inline
2683 {
2684  const Number array[1]
2685  = { t[0][0] };
2686  return SymmetricTensor<2,1,Number>(array);
2687 }
2688 
2689 
2690 
2699 template <typename Number>
2700 inline
2703 {
2704  const Number array[3]
2705  = { t[0][0], t[1][1], (t[0][1] + t[1][0])/2 };
2706  return SymmetricTensor<2,2,Number>(array);
2707 }
2708 
2709 
2710 
2719 template <typename Number>
2720 inline
2723 {
2724  const Number array[6]
2725  = { t[0][0], t[1][1], t[2][2],
2726  (t[0][1] + t[1][0])/2,
2727  (t[0][2] + t[2][0])/2,
2728  (t[1][2] + t[2][1])/2
2729  };
2730  return SymmetricTensor<2,3,Number>(array);
2731 }
2732 
2733 
2734 
2741 template <int rank, int dim, typename Number>
2742 inline
2745  const Number factor)
2746 {
2748  tt *= factor;
2749  return tt;
2750 }
2751 
2752 
2753 
2760 template <int rank, int dim, typename Number>
2761 inline
2763 operator * (const Number factor,
2765 {
2767  tt *= factor;
2768  return tt;
2769 }
2770 
2771 
2772 
2778 template <int rank, int dim, typename Number>
2779 inline
2782  const Number factor)
2783 {
2785  tt /= factor;
2786  return tt;
2787 }
2788 
2789 
2790 
2797 template <int rank, int dim>
2798 inline
2801  const double factor)
2802 {
2804  tt *= factor;
2805  return tt;
2806 }
2807 
2808 
2809 
2816 template <int rank, int dim>
2817 inline
2819 operator * (const double factor,
2820  const SymmetricTensor<rank,dim> &t)
2821 {
2823  tt *= factor;
2824  return tt;
2825 }
2826 
2827 
2828 
2834 template <int rank, int dim>
2835 inline
2838  const double factor)
2839 {
2841  tt /= factor;
2842  return tt;
2843 }
2844 
2854 template <int dim, typename Number>
2855 inline
2856 Number
2859 {
2860  return (t1*t2);
2861 }
2862 
2863 
2874 template <int dim, typename Number>
2875 inline
2876 Number
2878  const Tensor<2,dim,Number> &t2)
2879 {
2880  Number s = 0;
2881  for (unsigned int i=0; i<dim; ++i)
2882  for (unsigned int j=0; j<dim; ++j)
2883  s += t1[i][j] * t2[i][j];
2884  return s;
2885 }
2886 
2887 
2898 template <int dim, typename Number>
2899 inline
2900 Number
2903 {
2904  return scalar_product(t2, t1);
2905 }
2906 
2907 
2923 template <typename Number>
2924 inline
2925 void
2927  const SymmetricTensor<4,1,Number> &t,
2928  const SymmetricTensor<2,1,Number> &s)
2929 {
2930  tmp[0][0] = t[0][0][0][0] * s[0][0];
2931 }
2932 
2933 
2934 
2950 template <typename Number>
2951 inline
2952 void
2954  const SymmetricTensor<2,1,Number> &s,
2955  const SymmetricTensor<4,1,Number> &t)
2956 {
2957  tmp[0][0] = t[0][0][0][0] * s[0][0];
2958 }
2959 
2960 
2961 
2976 template <typename Number>
2977 inline
2978 void
2980  const SymmetricTensor<4,2,Number> &t,
2981  const SymmetricTensor<2,2,Number> &s)
2982 {
2983  const unsigned int dim = 2;
2984 
2985  for (unsigned int i=0; i<dim; ++i)
2986  for (unsigned int j=i; j<dim; ++j)
2987  tmp[i][j] = t[i][j][0][0] * s[0][0] +
2988  t[i][j][1][1] * s[1][1] +
2989  2 * t[i][j][0][1] * s[0][1];
2990 }
2991 
2992 
2993 
3009 template <typename Number>
3010 inline
3011 void
3013  const SymmetricTensor<2,2,Number> &s,
3014  const SymmetricTensor<4,2,Number> &t)
3015 {
3016  const unsigned int dim = 2;
3017 
3018  for (unsigned int i=0; i<dim; ++i)
3019  for (unsigned int j=i; j<dim; ++j)
3020  tmp[i][j] = s[0][0] * t[0][0][i][j] * +
3021  s[1][1] * t[1][1][i][j] +
3022  2 * s[0][1] * t[0][1][i][j];
3023 }
3024 
3025 
3026 
3042 template <typename Number>
3043 inline
3044 void
3046  const SymmetricTensor<4,3,Number> &t,
3047  const SymmetricTensor<2,3,Number> &s)
3048 {
3049  const unsigned int dim = 3;
3050 
3051  for (unsigned int i=0; i<dim; ++i)
3052  for (unsigned int j=i; j<dim; ++j)
3053  tmp[i][j] = t[i][j][0][0] * s[0][0] +
3054  t[i][j][1][1] * s[1][1] +
3055  t[i][j][2][2] * s[2][2] +
3056  2 * t[i][j][0][1] * s[0][1] +
3057  2 * t[i][j][0][2] * s[0][2] +
3058  2 * t[i][j][1][2] * s[1][2];
3059 }
3060 
3061 
3062 
3078 template <typename Number>
3079 inline
3080 void
3082  const SymmetricTensor<2,3,Number> &s,
3083  const SymmetricTensor<4,3,Number> &t)
3084 {
3085  const unsigned int dim = 3;
3086 
3087  for (unsigned int i=0; i<dim; ++i)
3088  for (unsigned int j=i; j<dim; ++j)
3089  tmp[i][j] = s[0][0] * t[0][0][i][j] +
3090  s[1][1] * t[1][1][i][j] +
3091  s[2][2] * t[2][2][i][j] +
3092  2 * s[0][1] * t[0][1][i][j] +
3093  2 * s[0][2] * t[0][2][i][j] +
3094  2 * s[1][2] * t[1][2][i][j];
3095 }
3096 
3097 
3098 
3115 template <int dim, typename Number>
3118  const Tensor<1,dim,Number> &src2)
3119 {
3120  Tensor<1,dim,Number> dest;
3121  for (unsigned int i=0; i<dim; ++i)
3122  for (unsigned int j=0; j<dim; ++j)
3123  dest[i] += src1[i][j] * src2[j];
3124  return dest;
3125 }
3126 
3127 
3137 template <int dim, typename Number>
3138 inline
3139 std::ostream &operator << (std::ostream &out,
3141 {
3142  //make out lives a bit simpler by outputing
3143  //the tensor through the operator for the
3144  //general Tensor class
3146 
3147  for (unsigned int i=0; i<dim; ++i)
3148  for (unsigned int j=0; j<dim; ++j)
3149  tt[i][j] = t[i][j];
3150 
3151  return out << tt;
3152 }
3153 
3154 
3155 
3165 template <int dim, typename Number>
3166 inline
3167 std::ostream &operator << (std::ostream &out,
3169 {
3170  //make out lives a bit simpler by outputing
3171  //the tensor through the operator for the
3172  //general Tensor class
3174 
3175  for (unsigned int i=0; i<dim; ++i)
3176  for (unsigned int j=0; j<dim; ++j)
3177  for (unsigned int k=0; k<dim; ++k)
3178  for (unsigned int l=0; l<dim; ++l)
3179  tt[i][j][k][l] = t[i][j][k][l];
3180 
3181  return out << tt;
3182 }
3183 
3184 
3185 DEAL_II_NAMESPACE_CLOSE
3186 
3187 #endif
VectorizedArray< Number > operator/(const VectorizedArray< Number > &u, const VectorizedArray< Number > &v)
Number scalar_product(const Tensor< 2, dim, Number > &t1, const Tensor< 2, dim, Number > &t2)
Definition: tensor.h:1610
Number second_invariant(const SymmetricTensor< 2, 1, Number > &)
Number trace(const SymmetricTensor< 2, dim, Number > &d)
VectorizedArray< Number > operator*(const VectorizedArray< Number > &u, const VectorizedArray< Number > &v)
SymmetricTensor< rank, dim, Number > transpose(const SymmetricTensor< rank, dim, Number > &t)
Number second_invariant(const SymmetricTensor< 2, 3, Number > &t)
SymmetricTensor operator+(const SymmetricTensor &s) const
void serialize(Archive &ar, const unsigned int version)
AccessorTypes< rank, dim, constness, Number >::reference reference
Number second_invariant(const SymmetricTensor< 2, 2, Number > &t)
::ExceptionBase & ExcMessage(std::string arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
internal::SymmetricTensorAccessors::StorageType< rank, dim, Number > base_tensor_descriptor
friend SymmetricTensor< 4, dim2, Number2 > deviator_tensor()
friend SymmetricTensor< 2, dim2, Number2 > unit_symmetric_tensor()
SymmetricTensor & operator*=(const Number factor)
SymmetricTensor & operator/=(const Number factor)
Number trace(const Tensor< 2, dim, Number > &d)
Definition: tensor.h:1759
TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
Number scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
SymmetricTensor< 2, dim, Number > deviator(const SymmetricTensor< 2, dim, Number > &t)
static std::size_t memory_consumption()
Number determinant(const Tensor< rank, 1, Number > &t)
Definition: tensor.h:1631
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Number norm() const
Number scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const Tensor< 2, dim, Number > &t2)
void double_contract(SymmetricTensor< 2, 3, Number > &tmp, const SymmetricTensor< 2, 3, Number > &s, const SymmetricTensor< 4, 3, Number > &t)
SymmetricTensor< 2, 3, Number > symmetrize(const Tensor< 2, 3, Number > &t)
base_tensor_type data
SymmetricTensor & operator+=(const SymmetricTensor &)
void double_contract(SymmetricTensor< 2, 2, Number > &tmp, const SymmetricTensor< 4, 2, Number > &t, const SymmetricTensor< 2, 2, Number > &s)
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1781
void double_contract(SymmetricTensor< 2, 1, Number > &tmp, const SymmetricTensor< 4, 1, Number > &t, const SymmetricTensor< 2, 1, Number > &s)
SymmetricTensor< 2, 1, Number > symmetrize(const Tensor< 2, 1, Number > &t)
T sum(const T &t, const MPI_Comm &mpi_communicator)
Definition: mpi.h:447
static TableIndices< rank > unrolled_to_component_indices(const unsigned int i)
#define Assert(cond, exc)
Definition: exceptions.h:299
void double_contract(SymmetricTensor< 2, 2, Number > &tmp, const SymmetricTensor< 2, 2, Number > &s, const SymmetricTensor< 4, 2, Number > &t)
static unsigned int component_to_unrolled_index(const TableIndices< rank > &indices)
bool operator==(const SymmetricTensor &) const
Definition: tensor.h:26
friend SymmetricTensor< 4, dim2, Number2 > identity_tensor()
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Accessor< rank, dim, constness, P-1, Number > operator[](const unsigned int i)
internal::SymmetricTensorAccessors::double_contraction_result< rank, 2, dim, Number >::type operator*(const SymmetricTensor< 2, dim, Number > &s) const
void double_contract(SymmetricTensor< 2, 3, Number > &tmp, const SymmetricTensor< 4, 3, Number > &t, const SymmetricTensor< 2, 3, Number > &s)
OS & operator<<(OS &o, const Event &e)
Definition: event.h:304
static const unsigned int n_independent_components
Number first_invariant(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
SymmetricTensor & operator-=(const SymmetricTensor &)
internal::SymmetricTensorAccessors::Accessor< rank, dim, true, rank-1, Number > operator[](const unsigned int row) const
Number scalar_product(const Tensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
SymmetricTensor< 2, 2, Number > symmetrize(const Tensor< 2, 2, Number > &t)
bool operator!=(const SymmetricTensor &) const
Tensor< 1, n_independent_components, Number > base_tensor_type
::ExceptionBase & ExcNotImplemented()
void double_contract(SymmetricTensor< 2, 1, Number > &tmp, const SymmetricTensor< 2, 1, Number > &s, const SymmetricTensor< 4, 1, Number > &t)
AccessorTypes< rank, dim, constness, Number >::reference reference
SymmetricTensor operator-() const
double third_invariant(const SymmetricTensor< 2, dim, Number > &t)
SymmetricTensor & operator=(const SymmetricTensor &)
::ExceptionBase & ExcInternalError()
base_tensor_descriptor::base_tensor_type base_tensor_type
Number access_raw_entry(const unsigned int unrolled_index) const
friend Number2 trace(const SymmetricTensor< 2, dim2, Number2 > &d)
Number & operator()(const TableIndices< rank > &indices)
static const unsigned int dimension