Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_transformations_h
16#define dealii_transformations_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/tensor.h>
22
24
25
26namespace Physics
27{
28 namespace Transformations
29 {
34 namespace Rotations
35 {
56 template <typename Number>
58 rotation_matrix_2d(const Number &angle);
59
60
89 template <typename Number>
91 rotation_matrix_3d(const Tensor<1, 3, Number> &axis, const Number &angle);
92
98 template <typename Number>
100 rotation_matrix_3d(const Point<3, Number> &axis, const Number &angle);
101
104 } // namespace Rotations
105
122 namespace Contravariant
123 {
142 template <int dim, typename Number>
145 const Tensor<2, dim, Number> &F);
146
161 template <int dim, typename Number>
164 const Tensor<2, dim, Number> &F);
165
181 template <int dim, typename Number>
184 const Tensor<2, dim, Number> &F);
185
200 template <int dim, typename Number>
203 const Tensor<2, dim, Number> &F);
204
220 template <int dim, typename Number>
223 const Tensor<2, dim, Number> &F);
224
245 template <int dim, typename Number>
248 const Tensor<2, dim, Number> &F);
249
264 template <int dim, typename Number>
267 const Tensor<2, dim, Number> &F);
268
283 template <int dim, typename Number>
286 const Tensor<2, dim, Number> &F);
287
302 template <int dim, typename Number>
305 const Tensor<2, dim, Number> &F);
306
321 template <int dim, typename Number>
324 const Tensor<2, dim, Number> &F);
325
327 } // namespace Contravariant
328
347 namespace Covariant
348 {
367 template <int dim, typename Number>
370 const Tensor<2, dim, Number> &F);
371
386 template <int dim, typename Number>
389 const Tensor<2, dim, Number> &F);
390
406 template <int dim, typename Number>
409 const Tensor<2, dim, Number> &F);
410
425 template <int dim, typename Number>
428 const Tensor<2, dim, Number> &F);
429
445 template <int dim, typename Number>
448 const Tensor<2, dim, Number> &F);
449
470 template <int dim, typename Number>
473 const Tensor<2, dim, Number> &F);
474
489 template <int dim, typename Number>
492 const Tensor<2, dim, Number> &F);
493
508 template <int dim, typename Number>
511 const Tensor<2, dim, Number> &F);
512
527 template <int dim, typename Number>
530 const Tensor<2, dim, Number> &F);
531
546 template <int dim, typename Number>
549 const Tensor<2, dim, Number> &F);
550
552 } // namespace Covariant
553
559 namespace Piola
560 {
581 template <int dim, typename Number>
584 const Tensor<2, dim, Number> &F);
585
601 template <int dim, typename Number>
604 const Tensor<2, dim, Number> &F);
605
622 template <int dim, typename Number>
625 const Tensor<2, dim, Number> &F);
626
643 template <int dim, typename Number>
646 const Tensor<2, dim, Number> &F);
647
665 template <int dim, typename Number>
668 const Tensor<2, dim, Number> &F);
669
692 template <int dim, typename Number>
695 const Tensor<2, dim, Number> &F);
696
712 template <int dim, typename Number>
715 const Tensor<2, dim, Number> &F);
716
732 template <int dim, typename Number>
735 const Tensor<2, dim, Number> &F);
736
753 template <int dim, typename Number>
756 const Tensor<2, dim, Number> &F);
757
774 template <int dim, typename Number>
777 const Tensor<2, dim, Number> &F);
778
780 } // namespace Piola
781
809 template <int dim, typename Number>
812 const Tensor<2, dim, Number> &F);
813
831 template <int dim, typename Number>
834 const Tensor<2, dim, Number> &B);
835
847 template <int dim, typename Number>
850 const Tensor<2, dim, Number> &B);
851
863 template <int dim, typename Number>
866 const Tensor<2, dim, Number> &B);
867
878 template <int dim, typename Number>
881 const Tensor<2, dim, Number> &B);
882
894 template <int dim, typename Number>
897 const Tensor<2, dim, Number> &B);
898
901 } // namespace Transformations
902} // namespace Physics
903
904
905
906#ifndef DOXYGEN
907
908
909
910template <typename Number>
913{
914 // Make things work with AD types
915 using std::cos;
916 using std::sin;
917
918 const Number rotation[2][2] = {{cos(angle), -sin(angle)},
919 {sin(angle), cos(angle)}};
920 return Tensor<2, 2>(rotation);
921}
922
923
924
925template <typename Number>
928 const Tensor<1, 3, Number> &axis,
929 const Number &angle)
930{
931 // Make things work with AD types
932 using std::abs;
933 using std::cos;
934 using std::sin;
935
936 Assert(abs(axis.norm() - 1.0) < 1e-9,
937 ExcMessage("The supplied axial vector is not a unit vector."));
938 const Number c = cos(angle);
939 const Number s = sin(angle);
940 const Number t = 1. - c;
941 const Number rotation[3][3] = {{t * axis[0] * axis[0] + c,
942 t * axis[0] * axis[1] - s * axis[2],
943 t * axis[0] * axis[2] + s * axis[1]},
944 {t * axis[0] * axis[1] + s * axis[2],
945 t * axis[1] * axis[1] + c,
946 t * axis[1] * axis[2] - s * axis[0]},
947 {t * axis[0] * axis[2] - s * axis[1],
948 t * axis[1] * axis[2] + s * axis[0],
949 t * axis[2] * axis[2] + c}};
950 return Tensor<2, 3, Number>(rotation);
951}
952
953
954
955template <typename Number>
958 const Point<3, Number> &axis,
959 const Number &angle)
960{
961 return rotation_matrix_3d(static_cast<Tensor<1, 3, Number>>(axis), angle);
962}
963
964
965
966template <int dim, typename Number>
969 const Tensor<1, dim, Number> &V,
970 const Tensor<2, dim, Number> &F)
971{
973}
974
975
976
977template <int dim, typename Number>
980 const Tensor<2, dim, Number> &T,
981 const Tensor<2, dim, Number> &F)
982{
984}
985
986
987
988template <int dim, typename Number>
992 const Tensor<2, dim, Number> &F)
993{
995}
996
997
998
999template <int dim, typename Number>
1002 const Tensor<4, dim, Number> &H,
1003 const Tensor<2, dim, Number> &F)
1004{
1006}
1007
1008
1009
1010template <int dim, typename Number>
1014 const Tensor<2, dim, Number> &F)
1015{
1017}
1018
1019
1020
1021template <int dim, typename Number>
1024 const Tensor<1, dim, Number> &v,
1025 const Tensor<2, dim, Number> &F)
1026{
1028}
1029
1030
1031
1032template <int dim, typename Number>
1035 const Tensor<2, dim, Number> &t,
1036 const Tensor<2, dim, Number> &F)
1037{
1039}
1040
1041
1042
1043template <int dim, typename Number>
1047 const Tensor<2, dim, Number> &F)
1048{
1050}
1051
1052
1053
1054template <int dim, typename Number>
1057 const Tensor<4, dim, Number> &h,
1058 const Tensor<2, dim, Number> &F)
1059{
1061}
1062
1063
1064
1065template <int dim, typename Number>
1069 const Tensor<2, dim, Number> &F)
1070{
1072}
1073
1074
1075
1076template <int dim, typename Number>
1079 const Tensor<1, dim, Number> &V,
1080 const Tensor<2, dim, Number> &F)
1081{
1083 transpose(invert(F)));
1084}
1085
1086
1087
1088template <int dim, typename Number>
1091 const Tensor<2, dim, Number> &T,
1092 const Tensor<2, dim, Number> &F)
1093{
1095 transpose(invert(F)));
1096}
1097
1098
1099
1100template <int dim, typename Number>
1104 const Tensor<2, dim, Number> &F)
1105{
1107 transpose(invert(F)));
1108}
1109
1110
1111
1112template <int dim, typename Number>
1115 const Tensor<4, dim, Number> &H,
1116 const Tensor<2, dim, Number> &F)
1117{
1119 transpose(invert(F)));
1120}
1121
1122
1123
1124template <int dim, typename Number>
1128 const Tensor<2, dim, Number> &F)
1129{
1131 transpose(invert(F)));
1132}
1133
1134
1135
1136template <int dim, typename Number>
1139 const Tensor<2, dim, Number> &F)
1140{
1142}
1143
1144
1145
1146template <int dim, typename Number>
1149 const Tensor<2, dim, Number> &F)
1150{
1152}
1153
1154
1155
1156template <int dim, typename Number>
1160 const Tensor<2, dim, Number> &F)
1161{
1163}
1164
1165
1166
1167template <int dim, typename Number>
1170 const Tensor<2, dim, Number> &F)
1171{
1173}
1174
1175
1176
1177template <int dim, typename Number>
1181 const Tensor<2, dim, Number> &F)
1182{
1184}
1185
1186
1187
1188template <int dim, typename Number>
1191 const Tensor<2, dim, Number> &F)
1192{
1193 return Number(1.0 / determinant(F)) * Contravariant::push_forward(V, F);
1194}
1195
1196
1197
1198template <int dim, typename Number>
1201 const Tensor<2, dim, Number> &F)
1202{
1203 return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1204}
1205
1206
1207
1208template <int dim, typename Number>
1212 const Tensor<2, dim, Number> &F)
1213{
1214 return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1215}
1216
1217
1218
1219template <int dim, typename Number>
1222 const Tensor<2, dim, Number> &F)
1223{
1224 return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
1225}
1226
1227
1228
1229template <int dim, typename Number>
1233 const Tensor<2, dim, Number> &F)
1234{
1235 return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
1236}
1237
1238
1239
1240template <int dim, typename Number>
1243 const Tensor<2, dim, Number> &F)
1244{
1245 return Number(determinant(F)) * Contravariant::pull_back(v, F);
1246}
1247
1248
1249
1250template <int dim, typename Number>
1253 const Tensor<2, dim, Number> &F)
1254{
1255 return Number(determinant(F)) * Contravariant::pull_back(t, F);
1256}
1257
1258
1259
1260template <int dim, typename Number>
1264 const Tensor<2, dim, Number> &F)
1265{
1266 return Number(determinant(F)) * Contravariant::pull_back(t, F);
1267}
1268
1269
1270
1271template <int dim, typename Number>
1274 const Tensor<2, dim, Number> &F)
1275{
1276 return Number(determinant(F)) * Contravariant::pull_back(h, F);
1277}
1278
1279
1280
1281template <int dim, typename Number>
1285 const Tensor<2, dim, Number> &F)
1286{
1287 return Number(determinant(F)) * Contravariant::pull_back(h, F);
1288}
1289
1290
1291
1292template <int dim, typename Number>
1295 const Tensor<2, dim, Number> &F)
1296{
1297 return cofactor(F) * N;
1298}
1299
1300
1301template <int dim, typename Number>
1304 const Tensor<2, dim, Number> &B)
1305{
1306 return contract<1, 0>(B, V);
1307}
1308
1309
1310
1311template <int dim, typename Number>
1314 const Tensor<2, dim, Number> &B)
1315{
1316 return contract<1, 0>(B, contract<1, 1>(T, B));
1317}
1318
1319
1320
1321template <int dim, typename Number>
1325 const Tensor<2, dim, Number> &B)
1326{
1328 for (unsigned int i = 0; i < dim; ++i)
1329 for (unsigned int J = 0; J < dim; ++J)
1330 // Loop over I but complex.h defines a macro I, so use I_ instead
1331 for (unsigned int I_ = 0; I_ < dim; ++I_)
1332 tmp_1[i][J] += B[i][I_] * T[I_][J];
1333
1335 for (unsigned int i = 0; i < dim; ++i)
1336 for (unsigned int j = i; j < dim; ++j)
1337 for (unsigned int J = 0; J < dim; ++J)
1338 out[i][j] += B[j][J] * tmp_1[i][J];
1339
1340 return out;
1341}
1342
1343
1344
1345template <int dim, typename Number>
1348 const Tensor<2, dim, Number> &B)
1349{
1350 // This contraction order and indexing might look a bit dubious, so a
1351 // quick explanation as to what's going on is probably in order:
1352 //
1353 // When the contract() function operates on the inner indices, the
1354 // result has the inner index and outer index transposed, i.e.
1355 // contract<2,1>(H,F) implies
1356 // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
1357 // rather than T_{IJkL} (the desired result).
1358 // So, in effect, contraction of the 3rd (inner) index with F as the
1359 // second argument results in its transposition with respect to its
1360 // adjacent neighbor. This is due to the position of the argument F,
1361 // leading to the free index being on the right hand side of the result.
1362 // However, given that we can do two transformations from the LHS of H
1363 // and two from the right we can undo the otherwise erroneous
1364 // swapping of the outer indices upon application of the second
1365 // sets of contractions.
1366 //
1367 // Note: Its significantly quicker (in 3d) to push forward
1368 // each index individually
1369 return contract<1, 1>(
1370 B, contract<1, 1>(B, contract<2, 1>(contract<2, 1>(H, B), B)));
1371}
1372
1373
1374
1375template <int dim, typename Number>
1379 const Tensor<2, dim, Number> &B)
1380{
1381 // The first and last transformation operations respectively
1382 // break and recover the symmetry properties of the tensors.
1383 // We also want to perform a minimal number of operations here
1384 // and avoid some complications related to the transposition of
1385 // tensor indices when contracting inner indices using the contract()
1386 // function. (For an explanation of the contraction operations,
1387 // please see the note in the equivalent function for standard
1388 // Tensors.) So what we'll do here is manually perform the first
1389 // and last contractions that break/recover the tensor symmetries
1390 // on the inner indices, and use the contract() function only on
1391 // the outer indices.
1392 //
1393 // Note: Its significantly quicker (in 3d) to push forward
1394 // each index individually
1395
1396 // Push forward (inner) index 1
1398 // Loop over I but complex.h defines a macro I, so use I_ instead
1399 for (unsigned int I_ = 0; I_ < dim; ++I_)
1400 for (unsigned int j = 0; j < dim; ++j)
1401 for (unsigned int K = 0; K < dim; ++K)
1402 for (unsigned int L = 0; L < dim; ++L)
1403 for (unsigned int J = 0; J < dim; ++J)
1404 tmp[I_][j][K][L] += B[j][J] * H[I_][J][K][L];
1405
1406 // Push forward (outer) indices 0 and 3
1407 tmp = contract<1, 0>(B, contract<3, 1>(tmp, B));
1408
1409 // Push forward (inner) index 2
1411 for (unsigned int i = 0; i < dim; ++i)
1412 for (unsigned int j = i; j < dim; ++j)
1413 for (unsigned int k = 0; k < dim; ++k)
1414 for (unsigned int l = k; l < dim; ++l)
1415 for (unsigned int K = 0; K < dim; ++K)
1416 out[i][j][k][l] += B[k][K] * tmp[i][j][K][l];
1417
1418 return out;
1419}
1420
1421#endif // DOXYGEN
1422
1424
1425#endif
Definition point.h:111
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_DEPRECATED
Definition config.h:207
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
const InputIterator OutputIterator out
Definition parallel.h:167
::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 > &)