Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.1.1
\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
transformations.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2019 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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_transformations_h
17 #define dealii_transformations_h
18 
19 #include <deal.II/base/point.h>
21 #include <deal.II/base/tensor.h>
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 
26 namespace Physics
27 {
28  namespace Transformations
29  {
36  namespace Rotations
37  {
42 
58  template <typename Number>
60  rotation_matrix_2d(const Number &angle);
61 
62 
91  template <typename Number>
93  rotation_matrix_3d(const Point<3, Number> &axis, const Number &angle);
94 
96 
97  } // namespace Rotations
98 
117  namespace Contravariant
118  {
123 
137  template <int dim, typename Number>
140  const Tensor<2, dim, Number> &F);
141 
156  template <int dim, typename Number>
159  const Tensor<2, dim, Number> &F);
160 
176  template <int dim, typename Number>
179  const Tensor<2, dim, Number> & F);
180 
195  template <int dim, typename Number>
198  const Tensor<2, dim, Number> &F);
199 
215  template <int dim, typename Number>
218  const Tensor<2, dim, Number> & F);
219 
221 
226 
240  template <int dim, typename Number>
243  const Tensor<2, dim, Number> &F);
244 
259  template <int dim, typename Number>
262  const Tensor<2, dim, Number> &F);
263 
278  template <int dim, typename Number>
281  const Tensor<2, dim, Number> & F);
282 
297  template <int dim, typename Number>
300  const Tensor<2, dim, Number> &F);
301 
316  template <int dim, typename Number>
319  const Tensor<2, dim, Number> & F);
320 
322  } // namespace Contravariant
323 
344  namespace Covariant
345  {
350 
364  template <int dim, typename Number>
367  const Tensor<2, dim, Number> &F);
368 
383  template <int dim, typename Number>
386  const Tensor<2, dim, Number> &F);
387 
403  template <int dim, typename Number>
406  const Tensor<2, dim, Number> & F);
407 
422  template <int dim, typename Number>
425  const Tensor<2, dim, Number> &F);
426 
442  template <int dim, typename Number>
445  const Tensor<2, dim, Number> & F);
446 
448 
453 
467  template <int dim, typename Number>
470  const Tensor<2, dim, Number> &F);
471 
486  template <int dim, typename Number>
489  const Tensor<2, dim, Number> &F);
490 
505  template <int dim, typename Number>
508  const Tensor<2, dim, Number> & F);
509 
524  template <int dim, typename Number>
527  const Tensor<2, dim, Number> &F);
528 
543  template <int dim, typename Number>
546  const Tensor<2, dim, Number> & F);
547 
549  } // namespace Covariant
550 
558  namespace Piola
559  {
564 
580  template <int dim, typename Number>
583  const Tensor<2, dim, Number> &F);
584 
600  template <int dim, typename Number>
603  const Tensor<2, dim, Number> &F);
604 
621  template <int dim, typename Number>
624  const Tensor<2, dim, Number> & F);
625 
642  template <int dim, typename Number>
645  const Tensor<2, dim, Number> &F);
646 
664  template <int dim, typename Number>
667  const Tensor<2, dim, Number> & F);
668 
670 
675 
691  template <int dim, typename Number>
694  const Tensor<2, dim, Number> &F);
695 
711  template <int dim, typename Number>
714  const Tensor<2, dim, Number> &F);
715 
731  template <int dim, typename Number>
734  const Tensor<2, dim, Number> & F);
735 
752  template <int dim, typename Number>
755  const Tensor<2, dim, Number> &F);
756 
773  template <int dim, typename Number>
776  const Tensor<2, dim, Number> & F);
777 
779  } // namespace Piola
780 
785 
808  template <int dim, typename Number>
811  const Tensor<2, dim, Number> &F);
812 
814  } // namespace Transformations
815 } // namespace Physics
816 
817 
818 
819 #ifndef DOXYGEN
820 
821 // ------------------------- inline functions ------------------------
822 
823 namespace internal
824 {
825  namespace Physics
826  {
827  template <int dim, typename Number>
829  transformation_contraction(const Tensor<1, dim, Number> &V,
830  const Tensor<2, dim, Number> &F)
831  {
832  return contract<1, 0>(F, V);
833  }
834 
835 
836 
837  template <int dim, typename Number>
839  transformation_contraction(const Tensor<2, dim, Number> &T,
840  const Tensor<2, dim, Number> &F)
841  {
842  return contract<1, 0>(F, contract<1, 1>(T, F));
843  }
844 
845 
846 
847  template <int dim, typename Number>
848  inline ::SymmetricTensor<2, dim, Number>
849  transformation_contraction(const ::SymmetricTensor<2, dim, Number> &T,
850  const Tensor<2, dim, Number> & F)
851  {
853  for (unsigned int i = 0; i < dim; ++i)
854  for (unsigned int J = 0; J < dim; ++J)
855  // Loop over I but complex.h defines a macro I, so use I_ instead
856  for (unsigned int I_ = 0; I_ < dim; ++I_)
857  tmp_1[i][J] += F[i][I_] * T[I_][J];
858 
860  for (unsigned int i = 0; i < dim; ++i)
861  for (unsigned int j = i; j < dim; ++j)
862  for (unsigned int J = 0; J < dim; ++J)
863  out[i][j] += F[j][J] * tmp_1[i][J];
864 
865  return out;
866  }
867 
868 
869 
870  template <int dim, typename Number>
872  transformation_contraction(const Tensor<4, dim, Number> &H,
873  const Tensor<2, dim, Number> &F)
874  {
875  // This contraction order and indexing might look a bit dubious, so a
876  // quick explanation as to what's going on is probably in order:
877  //
878  // When the contract() function operates on the inner indices, the
879  // result has the inner index and outer index transposed, i.e.
880  // contract<2,1>(H,F) implies
881  // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
882  // rather than T_{IJkL} (the desired result).
883  // So, in effect, contraction of the 3rd (inner) index with F as the
884  // second argument results in its transposition with respect to its
885  // adjacent neighbor. This is due to the position of the argument F,
886  // leading to the free index being on the right hand side of the result.
887  // However, given that we can do two transformations from the LHS of H
888  // and two from the right we can undo the otherwise erroneous
889  // swapping of the outer indices upon application of the second
890  // sets of contractions.
891  //
892  // Note: Its significantly quicker (in 3d) to push forward
893  // each index individually
894  return contract<1, 1>(
895  F, contract<1, 1>(F, contract<2, 1>(contract<2, 1>(H, F), F)));
896  }
897 
898 
899 
900  template <int dim, typename Number>
901  inline ::SymmetricTensor<4, dim, Number>
902  transformation_contraction(const ::SymmetricTensor<4, dim, Number> &H,
903  const Tensor<2, dim, Number> & F)
904  {
905  // The first and last transformation operations respectively
906  // break and recover the symmetry properties of the tensors.
907  // We also want to perform a minimal number of operations here
908  // and avoid some complications related to the transposition of
909  // tensor indices when contracting inner indices using the contract()
910  // function. (For an explanation of the contraction operations,
911  // please see the note in the equivalent function for standard
912  // Tensors.) So what we'll do here is manually perform the first
913  // and last contractions that break/recover the tensor symmetries
914  // on the inner indices, and use the contract() function only on
915  // the outer indices.
916  //
917  // Note: Its significantly quicker (in 3d) to push forward
918  // each index individually
919 
920  // Push forward (inner) index 1
922  // Loop over I but complex.h defines a macro I, so use I_ instead
923  for (unsigned int I_ = 0; I_ < dim; ++I_)
924  for (unsigned int j = 0; j < dim; ++j)
925  for (unsigned int K = 0; K < dim; ++K)
926  for (unsigned int L = 0; L < dim; ++L)
927  for (unsigned int J = 0; J < dim; ++J)
928  tmp[I_][j][K][L] += F[j][J] * H[I_][J][K][L];
929 
930  // Push forward (outer) indices 0 and 3
931  tmp = contract<1, 0>(F, contract<3, 1>(tmp, F));
932 
933  // Push forward (inner) index 2
935  for (unsigned int i = 0; i < dim; ++i)
936  for (unsigned int j = i; j < dim; ++j)
937  for (unsigned int k = 0; k < dim; ++k)
938  for (unsigned int l = k; l < dim; ++l)
939  for (unsigned int K = 0; K < dim; ++K)
940  out[i][j][k][l] += F[k][K] * tmp[i][j][K][l];
941 
942  return out;
943  }
944  } // namespace Physics
945 } // namespace internal
946 
947 
948 
949 template <typename Number>
952 {
953  const double rotation[2][2] = {{std::cos(angle), -std::sin(angle)},
954  {std::sin(angle), std::cos(angle)}};
955  return Tensor<2, 2>(rotation);
956 }
957 
958 
959 
960 template <typename Number>
963  const Point<3, Number> &axis,
964  const Number & angle)
965 {
966  Assert(std::abs(axis.norm() - 1.0) < 1e-9,
967  ExcMessage("The supplied axial vector is not a unit vector."));
968  const Number c = std::cos(angle);
969  const Number s = std::sin(angle);
970  const Number t = 1. - c;
971  const double rotation[3][3] = {{t * axis[0] * axis[0] + c,
972  t * axis[0] * axis[1] - s * axis[2],
973  t * axis[0] * axis[2] + s * axis[1]},
974  {t * axis[0] * axis[1] + s * axis[2],
975  t * axis[1] * axis[1] + c,
976  t * axis[1] * axis[2] - s * axis[0]},
977  {t * axis[0] * axis[2] - s * axis[1],
978  t * axis[1] * axis[2] + s * axis[0],
979  t * axis[2] * axis[2] + c}};
980  return Tensor<2, 3, Number>(rotation);
981 }
982 
983 
984 
985 template <int dim, typename Number>
988  const Tensor<1, dim, Number> &V,
989  const Tensor<2, dim, Number> &F)
990 {
991  return internal::Physics::transformation_contraction(V, F);
992 }
993 
994 
995 
996 template <int dim, typename Number>
999  const Tensor<2, dim, Number> &T,
1000  const Tensor<2, dim, Number> &F)
1001 {
1002  return internal::Physics::transformation_contraction(T, F);
1003 }
1004 
1005 
1006 
1007 template <int dim, typename Number>
1011  const Tensor<2, dim, Number> & F)
1012 {
1013  return internal::Physics::transformation_contraction(T, F);
1014 }
1015 
1016 
1017 
1018 template <int dim, typename Number>
1021  const Tensor<4, dim, Number> &H,
1022  const Tensor<2, dim, Number> &F)
1023 {
1024  return internal::Physics::transformation_contraction(H, F);
1025 }
1026 
1027 
1028 
1029 template <int dim, typename Number>
1033  const Tensor<2, dim, Number> & F)
1034 {
1035  return internal::Physics::transformation_contraction(H, F);
1036 }
1037 
1038 
1039 
1040 template <int dim, typename Number>
1043  const Tensor<1, dim, Number> &v,
1044  const Tensor<2, dim, Number> &F)
1045 {
1046  return internal::Physics::transformation_contraction(v, invert(F));
1047 }
1048 
1049 
1050 
1051 template <int dim, typename Number>
1054  const Tensor<2, dim, Number> &t,
1055  const Tensor<2, dim, Number> &F)
1056 {
1057  return internal::Physics::transformation_contraction(t, invert(F));
1058 }
1059 
1060 
1061 
1062 template <int dim, typename Number>
1066  const Tensor<2, dim, Number> & F)
1067 {
1068  return internal::Physics::transformation_contraction(t, invert(F));
1069 }
1070 
1071 
1072 
1073 template <int dim, typename Number>
1076  const Tensor<4, dim, Number> &h,
1077  const Tensor<2, dim, Number> &F)
1078 {
1079  return internal::Physics::transformation_contraction(h, invert(F));
1080 }
1081 
1082 
1083 
1084 template <int dim, typename Number>
1088  const Tensor<2, dim, Number> & F)
1089 {
1090  return internal::Physics::transformation_contraction(h, invert(F));
1091 }
1092 
1093 
1094 
1095 template <int dim, typename Number>
1098  const Tensor<1, dim, Number> &V,
1099  const Tensor<2, dim, Number> &F)
1100 {
1101  return internal::Physics::transformation_contraction(V, transpose(invert(F)));
1102 }
1103 
1104 
1105 
1106 template <int dim, typename Number>
1109  const Tensor<2, dim, Number> &T,
1110  const Tensor<2, dim, Number> &F)
1111 {
1112  return internal::Physics::transformation_contraction(T, transpose(invert(F)));
1113 }
1114 
1115 
1116 
1117 template <int dim, typename Number>
1121  const Tensor<2, dim, Number> & F)
1122 {
1123  return internal::Physics::transformation_contraction(T, transpose(invert(F)));
1124 }
1125 
1126 
1127 
1128 template <int dim, typename Number>
1131  const Tensor<4, dim, Number> &H,
1132  const Tensor<2, dim, Number> &F)
1133 {
1134  return internal::Physics::transformation_contraction(H, transpose(invert(F)));
1135 }
1136 
1137 
1138 
1139 template <int dim, typename Number>
1143  const Tensor<2, dim, Number> & F)
1144 {
1145  return internal::Physics::transformation_contraction(H, transpose(invert(F)));
1146 }
1147 
1148 
1149 
1150 template <int dim, typename Number>
1153  const Tensor<2, dim, Number> &F)
1154 {
1155  return internal::Physics::transformation_contraction(v, transpose(F));
1156 }
1157 
1158 
1159 
1160 template <int dim, typename Number>
1163  const Tensor<2, dim, Number> &F)
1164 {
1165  return internal::Physics::transformation_contraction(t, transpose(F));
1166 }
1167 
1168 
1169 
1170 template <int dim, typename Number>
1174  const Tensor<2, dim, Number> & F)
1175 {
1176  return internal::Physics::transformation_contraction(t, transpose(F));
1177 }
1178 
1179 
1180 
1181 template <int dim, typename Number>
1184  const Tensor<2, dim, Number> &F)
1185 {
1186  return internal::Physics::transformation_contraction(h, transpose(F));
1187 }
1188 
1189 
1190 
1191 template <int dim, typename Number>
1195  const Tensor<2, dim, Number> & F)
1196 {
1197  return internal::Physics::transformation_contraction(h, transpose(F));
1198 }
1199 
1200 
1201 
1202 template <int dim, typename Number>
1205  const Tensor<2, dim, Number> &F)
1206 {
1207  return Number(1.0 / determinant(F)) * Contravariant::push_forward(V, F);
1208 }
1209 
1210 
1211 
1212 template <int dim, typename Number>
1215  const Tensor<2, dim, Number> &F)
1216 {
1217  return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1218 }
1219 
1220 
1221 
1222 template <int dim, typename Number>
1226  const Tensor<2, dim, Number> & F)
1227 {
1228  return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1229 }
1230 
1231 
1232 
1233 template <int dim, typename Number>
1236  const Tensor<2, dim, Number> &F)
1237 {
1238  return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
1239 }
1240 
1241 
1242 
1243 template <int dim, typename Number>
1247  const Tensor<2, dim, Number> & F)
1248 {
1249  return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
1250 }
1251 
1252 
1253 
1254 template <int dim, typename Number>
1257  const Tensor<2, dim, Number> &F)
1258 {
1259  return Number(determinant(F)) * Contravariant::pull_back(v, F);
1260 }
1261 
1262 
1263 
1264 template <int dim, typename Number>
1267  const Tensor<2, dim, Number> &F)
1268 {
1269  return Number(determinant(F)) * Contravariant::pull_back(t, F);
1270 }
1271 
1272 
1273 
1274 template <int dim, typename Number>
1278  const Tensor<2, dim, Number> & F)
1279 {
1280  return Number(determinant(F)) * Contravariant::pull_back(t, F);
1281 }
1282 
1283 
1284 
1285 template <int dim, typename Number>
1288  const Tensor<2, dim, Number> &F)
1289 {
1290  return Number(determinant(F)) * Contravariant::pull_back(h, F);
1291 }
1292 
1293 
1294 
1295 template <int dim, typename Number>
1299  const Tensor<2, dim, Number> & F)
1300 {
1301  return Number(determinant(F)) * Contravariant::pull_back(h, F);
1302 }
1303 
1304 
1305 
1306 template <int dim, typename Number>
1309  const Tensor<2, dim, Number> &F)
1310 {
1311  return cofactor(F) * N;
1312 }
1313 
1314 #endif // DOXYGEN
1315 
1316 DEAL_II_NAMESPACE_CLOSE
1317 
1318 #endif
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
constexpr Tensor< 2, dim, Number > cofactor(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2341
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Tensor< 2, 2, Number > rotation_matrix_2d(const Number &angle)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Definition: point.h:110
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 1, dim, Number > nansons_formula(const Tensor< 1, dim, Number > &N, const Tensor< 2, dim, Number > &F)
Tensor< 1, dim, Number > push_forward(const Tensor< 1, dim, Number > &V, const Tensor< 2, dim, Number > &F)
static ::ExceptionBase & ExcMessage(std::string arg1)
Tensor< 2, 3, Number > rotation_matrix_3d(const Point< 3, Number > &axis, const Number &angle)
Tensor< 1, dim, Number > push_forward(const Tensor< 1, dim, Number > &V, const Tensor< 2, dim, Number > &F)
#define Assert(cond, exc)
Definition: exceptions.h:1407
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Tensor< 1, dim, Number > push_forward(const Tensor< 1, dim, Number > &V, const Tensor< 2, dim, Number > &F)
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
Definition: mpi.h:90
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)