deal.II version GIT relicensing-1927-g3de9220933 2024-10-03 08:40:00+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
95 } // namespace Rotations
96
113 namespace Contravariant
114 {
133 template <int dim, typename Number>
136 const Tensor<2, dim, Number> &F);
137
152 template <int dim, typename Number>
155 const Tensor<2, dim, Number> &F);
156
172 template <int dim, typename Number>
175 const Tensor<2, dim, Number> &F);
176
191 template <int dim, typename Number>
194 const Tensor<2, dim, Number> &F);
195
211 template <int dim, typename Number>
214 const Tensor<2, dim, Number> &F);
215
236 template <int dim, typename Number>
239 const Tensor<2, dim, Number> &F);
240
255 template <int dim, typename Number>
258 const Tensor<2, dim, Number> &F);
259
274 template <int dim, typename Number>
277 const Tensor<2, dim, Number> &F);
278
293 template <int dim, typename Number>
296 const Tensor<2, dim, Number> &F);
297
312 template <int dim, typename Number>
315 const Tensor<2, dim, Number> &F);
316
318 } // namespace Contravariant
319
338 namespace Covariant
339 {
358 template <int dim, typename Number>
361 const Tensor<2, dim, Number> &F);
362
377 template <int dim, typename Number>
380 const Tensor<2, dim, Number> &F);
381
397 template <int dim, typename Number>
400 const Tensor<2, dim, Number> &F);
401
416 template <int dim, typename Number>
419 const Tensor<2, dim, Number> &F);
420
436 template <int dim, typename Number>
439 const Tensor<2, dim, Number> &F);
440
461 template <int dim, typename Number>
464 const Tensor<2, dim, Number> &F);
465
480 template <int dim, typename Number>
483 const Tensor<2, dim, Number> &F);
484
499 template <int dim, typename Number>
502 const Tensor<2, dim, Number> &F);
503
518 template <int dim, typename Number>
521 const Tensor<2, dim, Number> &F);
522
537 template <int dim, typename Number>
540 const Tensor<2, dim, Number> &F);
541
543 } // namespace Covariant
544
550 namespace Piola
551 {
572 template <int dim, typename Number>
575 const Tensor<2, dim, Number> &F);
576
592 template <int dim, typename Number>
595 const Tensor<2, dim, Number> &F);
596
613 template <int dim, typename Number>
616 const Tensor<2, dim, Number> &F);
617
634 template <int dim, typename Number>
637 const Tensor<2, dim, Number> &F);
638
656 template <int dim, typename Number>
659 const Tensor<2, dim, Number> &F);
660
683 template <int dim, typename Number>
686 const Tensor<2, dim, Number> &F);
687
703 template <int dim, typename Number>
706 const Tensor<2, dim, Number> &F);
707
723 template <int dim, typename Number>
726 const Tensor<2, dim, Number> &F);
727
744 template <int dim, typename Number>
747 const Tensor<2, dim, Number> &F);
748
765 template <int dim, typename Number>
768 const Tensor<2, dim, Number> &F);
769
771 } // namespace Piola
772
800 template <int dim, typename Number>
803 const Tensor<2, dim, Number> &F);
804
822 template <int dim, typename Number>
825 const Tensor<2, dim, Number> &B);
826
838 template <int dim, typename Number>
841 const Tensor<2, dim, Number> &B);
842
854 template <int dim, typename Number>
857 const Tensor<2, dim, Number> &B);
858
869 template <int dim, typename Number>
872 const Tensor<2, dim, Number> &B);
873
885 template <int dim, typename Number>
888 const Tensor<2, dim, Number> &B);
889
892 } // namespace Transformations
893} // namespace Physics
894
895
896
897#ifndef DOXYGEN
898
899
900
901template <typename Number>
904{
905 // Make things work with AD types
906 using std::cos;
907 using std::sin;
908
909 const Number rotation[2][2] = {{cos(angle), -sin(angle)},
910 {sin(angle), cos(angle)}};
911 return Tensor<2, 2>(rotation);
912}
913
914
915
916template <typename Number>
919 const Tensor<1, 3, Number> &axis,
920 const Number &angle)
921{
922 // Make things work with AD types
923 using std::abs;
924 using std::cos;
925 using std::sin;
926
927 Assert(abs(axis.norm() - 1.0) < 1e-9,
928 ExcMessage("The supplied axial vector is not a unit vector."));
929 const Number c = cos(angle);
930 const Number s = sin(angle);
931 const Number t = 1. - c;
932 const Number rotation[3][3] = {{t * axis[0] * axis[0] + c,
933 t * axis[0] * axis[1] - s * axis[2],
934 t * axis[0] * axis[2] + s * axis[1]},
935 {t * axis[0] * axis[1] + s * axis[2],
936 t * axis[1] * axis[1] + c,
937 t * axis[1] * axis[2] - s * axis[0]},
938 {t * axis[0] * axis[2] - s * axis[1],
939 t * axis[1] * axis[2] + s * axis[0],
940 t * axis[2] * axis[2] + c}};
941 return Tensor<2, 3, Number>(rotation);
942}
943
944
945
946template <int dim, typename Number>
949 const Tensor<1, dim, Number> &V,
950 const Tensor<2, dim, Number> &F)
951{
953}
954
955
956
957template <int dim, typename Number>
960 const Tensor<2, dim, Number> &T,
961 const Tensor<2, dim, Number> &F)
962{
964}
965
966
967
968template <int dim, typename Number>
972 const Tensor<2, dim, Number> &F)
973{
975}
976
977
978
979template <int dim, typename Number>
982 const Tensor<4, dim, Number> &H,
983 const Tensor<2, dim, Number> &F)
984{
986}
987
988
989
990template <int dim, typename Number>
994 const Tensor<2, dim, Number> &F)
995{
997}
998
999
1000
1001template <int dim, typename Number>
1004 const Tensor<1, dim, Number> &v,
1005 const Tensor<2, dim, Number> &F)
1006{
1008}
1009
1010
1011
1012template <int dim, typename Number>
1015 const Tensor<2, dim, Number> &t,
1016 const Tensor<2, dim, Number> &F)
1017{
1019}
1020
1021
1022
1023template <int dim, typename Number>
1027 const Tensor<2, dim, Number> &F)
1028{
1030}
1031
1032
1033
1034template <int dim, typename Number>
1037 const Tensor<4, dim, Number> &h,
1038 const Tensor<2, dim, Number> &F)
1039{
1041}
1042
1043
1044
1045template <int dim, typename Number>
1049 const Tensor<2, dim, Number> &F)
1050{
1052}
1053
1054
1055
1056template <int dim, typename Number>
1059 const Tensor<1, dim, Number> &V,
1060 const Tensor<2, dim, Number> &F)
1061{
1063 transpose(invert(F)));
1064}
1065
1066
1067
1068template <int dim, typename Number>
1071 const Tensor<2, dim, Number> &T,
1072 const Tensor<2, dim, Number> &F)
1073{
1075 transpose(invert(F)));
1076}
1077
1078
1079
1080template <int dim, typename Number>
1084 const Tensor<2, dim, Number> &F)
1085{
1087 transpose(invert(F)));
1088}
1089
1090
1091
1092template <int dim, typename Number>
1095 const Tensor<4, dim, Number> &H,
1096 const Tensor<2, dim, Number> &F)
1097{
1099 transpose(invert(F)));
1100}
1101
1102
1103
1104template <int dim, typename Number>
1108 const Tensor<2, dim, Number> &F)
1109{
1111 transpose(invert(F)));
1112}
1113
1114
1115
1116template <int dim, typename Number>
1119 const Tensor<2, dim, Number> &F)
1120{
1122}
1123
1124
1125
1126template <int dim, typename Number>
1129 const Tensor<2, dim, Number> &F)
1130{
1132}
1133
1134
1135
1136template <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{
1173 return Number(1.0 / determinant(F)) * Contravariant::push_forward(V, F);
1174}
1175
1176
1177
1178template <int dim, typename Number>
1181 const Tensor<2, dim, Number> &F)
1182{
1183 return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
1184}
1185
1186
1187
1188template <int dim, typename Number>
1192 const Tensor<2, dim, Number> &F)
1193{
1194 return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, 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(H, 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(H, F);
1216}
1217
1218
1219
1220template <int dim, typename Number>
1223 const Tensor<2, dim, Number> &F)
1224{
1225 return Number(determinant(F)) * Contravariant::pull_back(v, F);
1226}
1227
1228
1229
1230template <int dim, typename Number>
1233 const Tensor<2, dim, Number> &F)
1234{
1235 return Number(determinant(F)) * Contravariant::pull_back(t, F);
1236}
1237
1238
1239
1240template <int dim, typename Number>
1244 const Tensor<2, dim, Number> &F)
1245{
1246 return Number(determinant(F)) * Contravariant::pull_back(t, 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(h, 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(h, F);
1268}
1269
1270
1271
1272template <int dim, typename Number>
1275 const Tensor<2, dim, Number> &F)
1276{
1277 return cofactor(F) * N;
1278}
1279
1280
1281template <int dim, typename Number>
1284 const Tensor<2, dim, Number> &B)
1285{
1286 return contract<1, 0>(B, V);
1287}
1288
1289
1290
1291template <int dim, typename Number>
1294 const Tensor<2, dim, Number> &B)
1295{
1296 return contract<1, 0>(B, contract<1, 1>(T, B));
1297}
1298
1299
1300
1301template <int dim, typename Number>
1305 const Tensor<2, dim, Number> &B)
1306{
1308 for (unsigned int i = 0; i < dim; ++i)
1309 for (unsigned int J = 0; J < dim; ++J)
1310 // Loop over I but complex.h defines a macro I, so use I_ instead
1311 for (unsigned int I_ = 0; I_ < dim; ++I_)
1312 tmp_1[i][J] += B[i][I_] * T[I_][J];
1313
1315 for (unsigned int i = 0; i < dim; ++i)
1316 for (unsigned int j = i; j < dim; ++j)
1317 for (unsigned int J = 0; J < dim; ++J)
1318 out[i][j] += B[j][J] * tmp_1[i][J];
1319
1320 return out;
1321}
1322
1323
1324
1325template <int dim, typename Number>
1328 const Tensor<2, dim, Number> &B)
1329{
1330 // This contraction order and indexing might look a bit dubious, so a
1331 // quick explanation as to what's going on is probably in order:
1332 //
1333 // When the contract() function operates on the inner indices, the
1334 // result has the inner index and outer index transposed, i.e.
1335 // contract<2,1>(H,F) implies
1336 // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
1337 // rather than T_{IJkL} (the desired result).
1338 // So, in effect, contraction of the 3rd (inner) index with F as the
1339 // second argument results in its transposition with respect to its
1340 // adjacent neighbor. This is due to the position of the argument F,
1341 // leading to the free index being on the right hand side of the result.
1342 // However, given that we can do two transformations from the LHS of H
1343 // and two from the right we can undo the otherwise erroneous
1344 // swapping of the outer indices upon application of the second
1345 // sets of contractions.
1346 //
1347 // Note: Its significantly quicker (in 3d) to push forward
1348 // each index individually
1349 return contract<1, 1>(
1350 B, contract<1, 1>(B, contract<2, 1>(contract<2, 1>(H, B), B)));
1351}
1352
1353
1354
1355template <int dim, typename Number>
1359 const Tensor<2, dim, Number> &B)
1360{
1361 // The first and last transformation operations respectively
1362 // break and recover the symmetry properties of the tensors.
1363 // We also want to perform a minimal number of operations here
1364 // and avoid some complications related to the transposition of
1365 // tensor indices when contracting inner indices using the contract()
1366 // function. (For an explanation of the contraction operations,
1367 // please see the note in the equivalent function for standard
1368 // Tensors.) So what we'll do here is manually perform the first
1369 // and last contractions that break/recover the tensor symmetries
1370 // on the inner indices, and use the contract() function only on
1371 // the outer indices.
1372 //
1373 // Note: Its significantly quicker (in 3d) to push forward
1374 // each index individually
1375
1376 // Push forward (inner) index 1
1378 // Loop over I but complex.h defines a macro I, so use I_ instead
1379 for (unsigned int I_ = 0; I_ < dim; ++I_)
1380 for (unsigned int j = 0; j < dim; ++j)
1381 for (unsigned int K = 0; K < dim; ++K)
1382 for (unsigned int L = 0; L < dim; ++L)
1383 for (unsigned int J = 0; J < dim; ++J)
1384 tmp[I_][j][K][L] += B[j][J] * H[I_][J][K][L];
1385
1386 // Push forward (outer) indices 0 and 3
1387 tmp = contract<1, 0>(B, contract<3, 1>(tmp, B));
1388
1389 // Push forward (inner) index 2
1391 for (unsigned int i = 0; i < dim; ++i)
1392 for (unsigned int j = i; j < dim; ++j)
1393 for (unsigned int k = 0; k < dim; ++k)
1394 for (unsigned int l = k; l < dim; ++l)
1395 for (unsigned int K = 0; K < dim; ++K)
1396 out[i][j][k][l] += B[k][K] * tmp[i][j][K][l];
1397
1398 return out;
1399}
1400
1401#endif // DOXYGEN
1402
1404
1405#endif
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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 > &)