Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
transformations.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2022 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/config.h>
20
22#include <deal.II/base/tensor.h>
23
25
26
27namespace Physics
28{
29 namespace Transformations
30 {
35 namespace Rotations
36 {
57 template <typename Number>
59 rotation_matrix_2d(const Number &angle);
60
61
90 template <typename Number>
92 rotation_matrix_3d(const Tensor<1, 3, Number> &axis, const Number &angle);
93
99 template <typename Number>
101 rotation_matrix_3d(const Point<3, Number> &axis, const Number &angle);
102
105 } // namespace Rotations
106
123 namespace Contravariant
124 {
143 template <int dim, typename Number>
146 const Tensor<2, dim, Number> &F);
147
162 template <int dim, typename Number>
165 const Tensor<2, dim, Number> &F);
166
182 template <int dim, typename Number>
185 const Tensor<2, dim, Number> & F);
186
201 template <int dim, typename Number>
204 const Tensor<2, dim, Number> &F);
205
221 template <int dim, typename Number>
224 const Tensor<2, dim, Number> & F);
225
246 template <int dim, typename Number>
249 const Tensor<2, dim, Number> &F);
250
265 template <int dim, typename Number>
268 const Tensor<2, dim, Number> &F);
269
284 template <int dim, typename Number>
287 const Tensor<2, dim, Number> & F);
288
303 template <int dim, typename Number>
306 const Tensor<2, dim, Number> &F);
307
322 template <int dim, typename Number>
325 const Tensor<2, dim, Number> & F);
326
328 } // namespace Contravariant
329
348 namespace Covariant
349 {
368 template <int dim, typename Number>
371 const Tensor<2, dim, Number> &F);
372
387 template <int dim, typename Number>
390 const Tensor<2, dim, Number> &F);
391
407 template <int dim, typename Number>
410 const Tensor<2, dim, Number> & F);
411
426 template <int dim, typename Number>
429 const Tensor<2, dim, Number> &F);
430
446 template <int dim, typename Number>
449 const Tensor<2, dim, Number> & F);
450
471 template <int dim, typename Number>
474 const Tensor<2, dim, Number> &F);
475
490 template <int dim, typename Number>
493 const Tensor<2, dim, Number> &F);
494
509 template <int dim, typename Number>
512 const Tensor<2, dim, Number> & F);
513
528 template <int dim, typename Number>
531 const Tensor<2, dim, Number> &F);
532
547 template <int dim, typename Number>
550 const Tensor<2, dim, Number> & F);
551
553 } // namespace Covariant
554
560 namespace Piola
561 {
582 template <int dim, typename Number>
585 const Tensor<2, dim, Number> &F);
586
602 template <int dim, typename Number>
605 const Tensor<2, dim, Number> &F);
606
623 template <int dim, typename Number>
626 const Tensor<2, dim, Number> & F);
627
644 template <int dim, typename Number>
647 const Tensor<2, dim, Number> &F);
648
666 template <int dim, typename Number>
669 const Tensor<2, dim, Number> & F);
670
693 template <int dim, typename Number>
696 const Tensor<2, dim, Number> &F);
697
713 template <int dim, typename Number>
716 const Tensor<2, dim, Number> &F);
717
733 template <int dim, typename Number>
736 const Tensor<2, dim, Number> & F);
737
754 template <int dim, typename Number>
757 const Tensor<2, dim, Number> &F);
758
775 template <int dim, typename Number>
778 const Tensor<2, dim, Number> & F);
779
781 } // namespace Piola
782
810 template <int dim, typename Number>
813 const Tensor<2, dim, Number> &F);
814
832 template <int dim, typename Number>
835 const Tensor<2, dim, Number> &B);
836
848 template <int dim, typename Number>
851 const Tensor<2, dim, Number> &B);
852
864 template <int dim, typename Number>
867 const Tensor<2, dim, Number> & B);
868
879 template <int dim, typename Number>
882 const Tensor<2, dim, Number> &B);
883
895 template <int dim, typename Number>
898 const Tensor<2, dim, Number> & B);
899
902 } // namespace Transformations
903} // namespace Physics
904
905
906
907#ifndef DOXYGEN
908
909
910
911template <typename Number>
914{
915 // Make things work with AD types
916 using std::cos;
917 using std::sin;
918
919 const Number rotation[2][2] = {{cos(angle), -sin(angle)},
920 {sin(angle), cos(angle)}};
921 return Tensor<2, 2>(rotation);
922}
923
924
925
926template <typename Number>
929 const Tensor<1, 3, Number> &axis,
930 const Number & angle)
931{
932 // Make things work with AD types
933 using std::abs;
934 using std::cos;
935 using std::sin;
936
937 Assert(abs(axis.norm() - 1.0) < 1e-9,
938 ExcMessage("The supplied axial vector is not a unit vector."));
939 const Number c = cos(angle);
940 const Number s = sin(angle);
941 const Number t = 1. - c;
942 const Number rotation[3][3] = {{t * axis[0] * axis[0] + c,
943 t * axis[0] * axis[1] - s * axis[2],
944 t * axis[0] * axis[2] + s * axis[1]},
945 {t * axis[0] * axis[1] + s * axis[2],
946 t * axis[1] * axis[1] + c,
947 t * axis[1] * axis[2] - s * axis[0]},
948 {t * axis[0] * axis[2] - s * axis[1],
949 t * axis[1] * axis[2] + s * axis[0],
950 t * axis[2] * axis[2] + c}};
951 return Tensor<2, 3, Number>(rotation);
952}
953
954
955
956template <typename Number>
959 const Point<3, Number> &axis,
960 const Number & angle)
961{
962 return rotation_matrix_3d(static_cast<Tensor<1, 3, Number>>(axis), angle);
963}
964
965
966
967template <int dim, typename Number>
970 const Tensor<1, dim, Number> &V,
971 const Tensor<2, dim, Number> &F)
972{
974}
975
976
977
978template <int dim, typename Number>
981 const Tensor<2, dim, Number> &T,
982 const Tensor<2, dim, Number> &F)
983{
985}
986
987
988
989template <int dim, typename Number>
993 const Tensor<2, dim, Number> & F)
994{
996}
997
998
999
1000template <int dim, typename Number>
1003 const Tensor<4, dim, Number> &H,
1004 const Tensor<2, dim, Number> &F)
1005{
1007}
1008
1009
1010
1011template <int dim, typename Number>
1015 const Tensor<2, dim, Number> & F)
1016{
1018}
1019
1020
1021
1022template <int dim, typename Number>
1025 const Tensor<1, dim, Number> &v,
1026 const Tensor<2, dim, Number> &F)
1027{
1029}
1030
1031
1032
1033template <int dim, typename Number>
1036 const Tensor<2, dim, Number> &t,
1037 const Tensor<2, dim, Number> &F)
1038{
1040}
1041
1042
1043
1044template <int dim, typename Number>
1048 const Tensor<2, dim, Number> & F)
1049{
1051}
1052
1053
1054
1055template <int dim, typename Number>
1058 const Tensor<4, dim, Number> &h,
1059 const Tensor<2, dim, Number> &F)
1060{
1062}
1063
1064
1065
1066template <int dim, typename Number>
1070 const Tensor<2, dim, Number> & F)
1071{
1073}
1074
1075
1076
1077template <int dim, typename Number>
1080 const Tensor<1, dim, Number> &V,
1081 const Tensor<2, dim, Number> &F)
1082{
1084 transpose(invert(F)));
1085}
1086
1087
1088
1089template <int dim, typename Number>
1092 const Tensor<2, dim, Number> &T,
1093 const Tensor<2, dim, Number> &F)
1094{
1096 transpose(invert(F)));
1097}
1098
1099
1100
1101template <int dim, typename Number>
1105 const Tensor<2, dim, Number> & F)
1106{
1108 transpose(invert(F)));
1109}
1110
1111
1112
1113template <int dim, typename Number>
1116 const Tensor<4, dim, Number> &H,
1117 const Tensor<2, dim, Number> &F)
1118{
1120 transpose(invert(F)));
1121}
1122
1123
1124
1125template <int dim, typename Number>
1129 const Tensor<2, dim, Number> & F)
1130{
1132 transpose(invert(F)));
1133}
1134
1135
1136
1137template <int dim, typename Number>
1140 const Tensor<2, dim, Number> &F)
1141{
1143}
1144
1145
1146
1147template <int dim, typename Number>
1150 const Tensor<2, dim, Number> &F)
1151{
1153}
1154
1155
1156
1157template <int dim, typename Number>
1161 const Tensor<2, dim, Number> & F)
1162{
1164}
1165
1166
1167
1168template <int dim, typename Number>
1171 const Tensor<2, dim, Number> &F)
1172{
1174}
1175
1176
1177
1178template <int dim, typename Number>
1182 const Tensor<2, dim, Number> & F)
1183{
1185}
1186
1187
1188
1189template <int dim, typename Number>
1192 const Tensor<2, dim, Number> &F)
1193{
1194 return Number(1.0 / determinant(F)) * Contravariant::push_forward(V, F);
1195}
1196
1197
1198
1199template <int dim, typename Number>
1202 const Tensor<2, dim, Number> &F)
1203{
1204 return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1205}
1206
1207
1208
1209template <int dim, typename Number>
1213 const Tensor<2, dim, Number> & F)
1214{
1215 return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1216}
1217
1218
1219
1220template <int dim, typename Number>
1223 const Tensor<2, dim, Number> &F)
1224{
1225 return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
1226}
1227
1228
1229
1230template <int dim, typename Number>
1234 const Tensor<2, dim, Number> & F)
1235{
1236 return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
1237}
1238
1239
1240
1241template <int dim, typename Number>
1244 const Tensor<2, dim, Number> &F)
1245{
1246 return Number(determinant(F)) * Contravariant::pull_back(v, F);
1247}
1248
1249
1250
1251template <int dim, typename Number>
1254 const Tensor<2, dim, Number> &F)
1255{
1256 return Number(determinant(F)) * Contravariant::pull_back(t, F);
1257}
1258
1259
1260
1261template <int dim, typename Number>
1265 const Tensor<2, dim, Number> & F)
1266{
1267 return Number(determinant(F)) * Contravariant::pull_back(t, F);
1268}
1269
1270
1271
1272template <int dim, typename Number>
1275 const Tensor<2, dim, Number> &F)
1276{
1277 return Number(determinant(F)) * Contravariant::pull_back(h, F);
1278}
1279
1280
1281
1282template <int dim, typename Number>
1286 const Tensor<2, dim, Number> & F)
1287{
1288 return Number(determinant(F)) * Contravariant::pull_back(h, F);
1289}
1290
1291
1292
1293template <int dim, typename Number>
1296 const Tensor<2, dim, Number> &F)
1297{
1298 return cofactor(F) * N;
1299}
1300
1301
1302template <int dim, typename Number>
1305 const Tensor<2, dim, Number> &B)
1306{
1307 return contract<1, 0>(B, V);
1308}
1309
1310
1311
1312template <int dim, typename Number>
1315 const Tensor<2, dim, Number> &B)
1316{
1317 return contract<1, 0>(B, contract<1, 1>(T, B));
1318}
1319
1320
1321
1322template <int dim, typename Number>
1326 const Tensor<2, dim, Number> & B)
1327{
1329 for (unsigned int i = 0; i < dim; ++i)
1330 for (unsigned int J = 0; J < dim; ++J)
1331 // Loop over I but complex.h defines a macro I, so use I_ instead
1332 for (unsigned int I_ = 0; I_ < dim; ++I_)
1333 tmp_1[i][J] += B[i][I_] * T[I_][J];
1334
1336 for (unsigned int i = 0; i < dim; ++i)
1337 for (unsigned int j = i; j < dim; ++j)
1338 for (unsigned int J = 0; J < dim; ++J)
1339 out[i][j] += B[j][J] * tmp_1[i][J];
1340
1341 return out;
1342}
1343
1344
1345
1346template <int dim, typename Number>
1349 const Tensor<2, dim, Number> &B)
1350{
1351 // This contraction order and indexing might look a bit dubious, so a
1352 // quick explanation as to what's going on is probably in order:
1353 //
1354 // When the contract() function operates on the inner indices, the
1355 // result has the inner index and outer index transposed, i.e.
1356 // contract<2,1>(H,F) implies
1357 // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
1358 // rather than T_{IJkL} (the desired result).
1359 // So, in effect, contraction of the 3rd (inner) index with F as the
1360 // second argument results in its transposition with respect to its
1361 // adjacent neighbor. This is due to the position of the argument F,
1362 // leading to the free index being on the right hand side of the result.
1363 // However, given that we can do two transformations from the LHS of H
1364 // and two from the right we can undo the otherwise erroneous
1365 // swapping of the outer indices upon application of the second
1366 // sets of contractions.
1367 //
1368 // Note: Its significantly quicker (in 3d) to push forward
1369 // each index individually
1370 return contract<1, 1>(
1371 B, contract<1, 1>(B, contract<2, 1>(contract<2, 1>(H, B), B)));
1372}
1373
1374
1375
1376template <int dim, typename Number>
1380 const Tensor<2, dim, Number> & B)
1381{
1382 // The first and last transformation operations respectively
1383 // break and recover the symmetry properties of the tensors.
1384 // We also want to perform a minimal number of operations here
1385 // and avoid some complications related to the transposition of
1386 // tensor indices when contracting inner indices using the contract()
1387 // function. (For an explanation of the contraction operations,
1388 // please see the note in the equivalent function for standard
1389 // Tensors.) So what we'll do here is manually perform the first
1390 // and last contractions that break/recover the tensor symmetries
1391 // on the inner indices, and use the contract() function only on
1392 // the outer indices.
1393 //
1394 // Note: Its significantly quicker (in 3d) to push forward
1395 // each index individually
1396
1397 // Push forward (inner) index 1
1399 // Loop over I but complex.h defines a macro I, so use I_ instead
1400 for (unsigned int I_ = 0; I_ < dim; ++I_)
1401 for (unsigned int j = 0; j < dim; ++j)
1402 for (unsigned int K = 0; K < dim; ++K)
1403 for (unsigned int L = 0; L < dim; ++L)
1404 for (unsigned int J = 0; J < dim; ++J)
1405 tmp[I_][j][K][L] += B[j][J] * H[I_][J][K][L];
1406
1407 // Push forward (outer) indices 0 and 3
1408 tmp = contract<1, 0>(B, contract<3, 1>(tmp, B));
1409
1410 // Push forward (inner) index 2
1412 for (unsigned int i = 0; i < dim; ++i)
1413 for (unsigned int j = i; j < dim; ++j)
1414 for (unsigned int k = 0; k < dim; ++k)
1415 for (unsigned int l = k; l < dim; ++l)
1416 for (unsigned int K = 0; K < dim; ++K)
1417 out[i][j][k][l] += B[k][K] * tmp[i][j][K][l];
1418
1419 return out;
1420}
1421
1422#endif // DOXYGEN
1423
1425
1426#endif
Definition point.h:112
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_DEPRECATED
Definition config.h:172
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
Tensor< 1, dim, Number > push_forward(const Tensor< 1, dim, Number > &V, const Tensor< 2, dim, Number > &F)
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)
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)
Tensor< 2, 3, Number > rotation_matrix_3d(const Tensor< 1, 3, Number > &axis, const Number &angle)
Tensor< 2, 2, Number > rotation_matrix_2d(const Number &angle)
Tensor< 1, dim, Number > nansons_formula(const Tensor< 1, dim, Number > &N, const Tensor< 2, dim, Number > &F)
Tensor< 1, dim, Number > basis_transformation(const Tensor< 1, dim, Number > &V, const Tensor< 2, dim, Number > &B)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)