Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "Moose.h"
13 : #include "libmesh/vector_value.h"
14 : #include "libmesh/tensor_value.h"
15 : #include "MooseTypes.h"
16 : /**
17 : * Utility functions to return rotations matrics
18 : */
19 : namespace RotationMatrix
20 : {
21 : /// provides a rotation matrix that will rotate the vector vec to the z axis (the "2" direction)
22 : template <bool is_ad = false>
23 : GenericRealTensorValue<is_ad>
24 1110 : rotVecToZ(GenericRealVectorValue<is_ad> vec)
25 : // provides a rotation matrix that will rotate the vector vec to the z axis (the "2" direction)
26 : {
27 : using std::sqrt, std::abs;
28 : // ensure that vec is normalised
29 1110 : vec /= sqrt(vec * vec);
30 :
31 : // construct v0 and v1 to be orthonormal to vec
32 : // and form a RH basis, that is, so v1 x vec = v0
33 :
34 : // Use Gram-Schmidt method to find v1.
35 1110 : GenericRealVectorValue<is_ad> v1;
36 : // Need a prototype for v1 first, and this is done by looking at the smallest component of vec
37 1110 : GenericRealVectorValue<is_ad> w(abs(vec(0)), abs(vec(1)), abs(vec(2)));
38 1110 : if ((w(2) >= w(1) && w(1) >= w(0)) || (w(1) >= w(2) && w(2) >= w(0)))
39 : // vec(0) is the smallest component
40 1044 : v1(0) = 1;
41 66 : else if ((w(2) >= w(0) && w(0) >= w(1)) || (w(0) >= w(2) && w(2) >= w(1)))
42 : // vec(1) is the smallest component
43 46 : v1(1) = 1;
44 : else
45 : // vec(2) is the smallest component
46 20 : v1(2) = 1;
47 : // now Gram-Schmidt
48 1110 : v1 -= (v1 * vec) * vec;
49 1110 : v1 /= sqrt(v1 * v1);
50 :
51 : // now use v0 = v1 x vec
52 1110 : GenericRealVectorValue<is_ad> v0;
53 1110 : v0(0) = v1(1) * vec(2) - v1(2) * vec(1);
54 1110 : v0(1) = v1(2) * vec(0) - v1(0) * vec(2);
55 1110 : v0(2) = v1(0) * vec(1) - v1(1) * vec(0);
56 :
57 : // the desired rotation matrix is just
58 1110 : GenericRealTensorValue<is_ad> rot(
59 1110 : v0(0), v0(1), v0(2), v1(0), v1(1), v1(2), vec(0), vec(1), vec(2));
60 2220 : return rot;
61 : }
62 :
63 : /// provides a rotation matrix that will rotate the vector vec1 to vec2
64 : template <bool is_ad = false>
65 : GenericRealTensorValue<is_ad>
66 52 : rotVec1ToVec2(GenericRealVectorValue<is_ad> vec1, GenericRealVectorValue<is_ad> vec2)
67 : // provides a rotation matrix that will rotate the vector vec1 to the vector vec2
68 : {
69 52 : GenericRealTensorValue<is_ad> rot1_to_z = rotVecToZ<is_ad>(vec1);
70 52 : GenericRealTensorValue<is_ad> rot2_to_z = rotVecToZ<is_ad>(vec2);
71 104 : return rot2_to_z.transpose() * rot1_to_z;
72 52 : }
73 :
74 : /// provides a rotation matrix that will rotate the vector vec1 to the [1,0,0], assuming vec1[2]==0
75 : template <bool is_ad = false>
76 : GenericRealTensorValue<is_ad>
77 24 : rotVec2DToX(const GenericRealVectorValue<is_ad> & vec)
78 : // provides a rotation matrix that will rotate the vector `vec` to the [1,0,0], assuming vec[2]==0
79 : {
80 : using std::atan2, std::sin, std::cos;
81 24 : const GenericReal<is_ad> theta = atan2(vec(1), vec(0));
82 24 : const GenericReal<is_ad> st = sin(theta);
83 24 : const GenericReal<is_ad> ct = cos(theta);
84 24 : return GenericRealTensorValue<is_ad>(ct, st, 0., -st, ct, 0., 0., 0., 1.);
85 : }
86 :
87 : /**
88 : * Provides rotatiom matrix for rotating from vec1 to vec2 using Rodrigues' rotation forumula.
89 : * See https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Matrix_notation
90 : * @param vec1 starting vector -- must have 3 components!
91 : * @param vec2 ending vector -- must have 3 components!
92 : * @return 3x3 rotation tensor (matrix)
93 : */
94 : template <bool is_ad = false>
95 : GenericRealTensorValue<is_ad>
96 386 : rodriguesRotationMatrix(GenericRealVectorValue<is_ad> vec1, GenericRealVectorValue<is_ad> vec2)
97 : {
98 : // normalize input vectors
99 386 : GenericRealVectorValue<is_ad> u = vec1 / vec1.norm();
100 386 : GenericRealVectorValue<is_ad> v = vec2 / vec2.norm();
101 :
102 386 : if ((u - v).norm() < libMesh::TOLERANCE)
103 232 : return GenericRealTensorValue<is_ad>(1, 0, 0, 0, 1, 0, 0, 0, 1); // identity matrix
104 :
105 270 : GenericRealVectorValue<is_ad> k_vec = u.cross(v); // calculate rotation axis
106 :
107 : // test for exactly opposite vectors to avoid divide by zero
108 270 : if (k_vec.norm() < libMesh::TOLERANCE)
109 : {
110 18 : GenericRealTensorValue<is_ad> rot_matrix;
111 18 : if ((u - GenericRealVectorValue<is_ad>(1, 0, 0)).norm() < TOLERANCE)
112 18 : rot_matrix = GenericRealTensorValue<is_ad>(-1, 0, 0, 0, -1, 0, 0, 0, 1);
113 : else
114 0 : mooseError("Rotation matrix cannot be generated for opposite-facing vectors at this time!");
115 18 : return rot_matrix;
116 18 : }
117 :
118 252 : k_vec /= k_vec.norm(); // normalize
119 252 : Real cos_theta = u * v;
120 252 : Real theta = std::acos(cos_theta);
121 252 : Real sin_theta = std::sin(theta);
122 :
123 252 : GenericRealTensorValue<is_ad> K_matrix(
124 252 : 0, -k_vec(2), k_vec(1), k_vec(2), 0, -k_vec(0), -k_vec(1), k_vec(0), 0);
125 252 : GenericRealTensorValue<is_ad> I(1, 0, 0, 0, 1, 0, 0, 0, 1); // identity matrix
126 :
127 : // construct rotation matrix
128 252 : GenericRealTensorValue<is_ad> rot_matrix;
129 252 : rot_matrix = I + sin_theta * K_matrix + (1 - cos_theta) * K_matrix * K_matrix;
130 252 : return rot_matrix;
131 252 : }
132 : }
|
