@@ -22,52 +22,61 @@ struct CalculateRotationControlConstraint {
22
22
23
23
KOKKOS_FUNCTION
24
24
void operator ()(const int i_constraint) const {
25
- const auto i_node1 = base_node_index (i_constraint);
26
- const auto i_node2 = target_node_index (i_constraint);
25
+ const auto i_base = base_node_index (i_constraint);
26
+ const auto i_target = target_node_index (i_constraint);
27
27
28
28
// Initial difference between nodes
29
29
const auto X0_data = Kokkos::Array<double , 3 >{
30
30
X0_ (i_constraint, 0 ), X0_ (i_constraint, 1 ), X0_ (i_constraint, 2 )
31
31
};
32
32
const auto X0 = View_3::const_type{X0_data.data ()};
33
33
34
- // Base node displacement
35
- const auto u1_data =
36
- Kokkos::Array<double , 3 >{node_u (i_node1 , 0 ), node_u (i_node1 , 1 ), node_u (i_node1 , 2 )};
37
- const auto R1_data = Kokkos::Array<double , 4 >{
38
- node_u (i_node1 , 3 ), node_u (i_node1 , 4 ), node_u (i_node1 , 5 ), node_u (i_node1 , 6 )
34
+ // Base node displacement (translation and rotation)
35
+ const auto ub_data =
36
+ Kokkos::Array<double , 3 >{node_u (i_base , 0 ), node_u (i_base , 1 ), node_u (i_base , 2 )};
37
+ const auto Rb_data = Kokkos::Array<double , 4 >{
38
+ node_u (i_base , 3 ), node_u (i_base , 4 ), node_u (i_base , 5 ), node_u (i_base , 6 )
39
39
};
40
- const auto u1 = View_3::const_type{u1_data.data ()};
41
- const auto R1 = Kokkos::View<double [4 ]>::const_type{R1_data.data ()};
42
-
43
- // Target node displacement
44
- const auto R2_data = Kokkos::Array<double , 4 >{
45
- node_u (i_node2, 3 ), node_u (i_node2, 4 ), node_u (i_node2, 5 ), node_u (i_node2, 6 )
40
+ const auto ub = View_3::const_type{ub_data.data ()};
41
+ const auto Rb = Kokkos::View<double [4 ]>::const_type{Rb_data.data ()};
42
+
43
+ // Target node displacement (translation and rotation)
44
+ const auto ut_data =
45
+ Kokkos::Array<double , 3 >{node_u (i_target, 0 ), node_u (i_target, 1 ), node_u (i_target, 2 )};
46
+ const auto ut = View_3::const_type{ut_data.data ()};
47
+ const auto Rt_data = Kokkos::Array<double , 4 >{
48
+ node_u (i_target, 3 ), node_u (i_target, 4 ), node_u (i_target, 5 ), node_u (i_target, 6 )
46
49
};
47
- const auto R2 = Kokkos::View<double [4 ]>::const_type{R2_data.data ()};
48
- const auto u2_data =
49
- Kokkos::Array<double , 3 >{node_u (i_node2, 0 ), node_u (i_node2, 1 ), node_u (i_node2, 2 )};
50
- const auto u2 = View_3::const_type{u2_data.data ()};
50
+ const auto Rt = Kokkos::View<double [4 ]>::const_type{Rt_data.data ()};
51
+
52
+ auto AX_data = Kokkos::Array<double , 3 >{};
53
+ const auto AX = Kokkos::View<double [3 ]>{AX_data.data ()};
54
+
55
+ // Control rotation vector
56
+ auto RV_data = Kokkos::Array<double , 3 >{};
57
+ const auto RV = Kokkos::View<double [3 ]>{RV_data.data ()};
51
58
52
59
// Rotation control
53
- auto RC_data = Kokkos::Array<double , 4 >{};
54
- const auto RC = Kokkos::View<double [4 ]>{RC_data.data ()};
55
- auto RCt_data = Kokkos::Array<double , 4 >{};
56
- const auto RCt = Kokkos::View<double [4 ]>{RCt_data.data ()};
57
- auto RV_data = Kokkos::Array<double , 4 >{};
58
- const auto RV = Kokkos::View<double [4 ]>{RV_data.data ()};
60
+ auto Rc_data = Kokkos::Array<double , 4 >{};
61
+ const auto Rc = Kokkos::View<double [4 ]>{Rc_data.data ()};
62
+ auto RcT_data = Kokkos::Array<double , 4 >{};
63
+ const auto RcT = Kokkos::View<double [4 ]>{RcT_data.data ()};
59
64
60
- auto R1t_data = Kokkos::Array<double , 4 >{};
61
- const auto R1t = Kokkos::View<double [4 ]>{R1t_data.data ()};
65
+ // Base rotation transpose
66
+ auto RbT_data = Kokkos::Array<double , 4 >{};
67
+ const auto RbT = Kokkos::View<double [4 ]>{RbT_data.data ()};
62
68
63
- auto R1_X0_data = Kokkos::Array<double , 4 >{};
64
- const auto R1_X0 = Kokkos::View<double [4 ]>{R1_X0_data.data ()};
69
+ // Base rotation * X0
70
+ auto Rb_X0_data = Kokkos::Array<double , 4 >{};
71
+ const auto Rb_X0 = Kokkos::View<double [4 ]>{Rb_X0_data.data ()};
65
72
66
- auto R2_RCt_data = Kokkos::Array<double , 4 >{};
67
- const auto R2_RCt = Kokkos::View<double [4 ]>{R2_RCt_data.data ()};
73
+ // Target rotation * Control rotation transpose
74
+ auto Rt_RcT_data = Kokkos::Array<double , 4 >{};
75
+ const auto Rt_RcT = Kokkos::View<double [4 ]>{Rt_RcT_data.data ()};
68
76
69
- auto R2_RCt_R1t_data = Kokkos::Array<double , 4 >{};
70
- const auto R2_RCt_R1t = Kokkos::View<double [4 ]>{R2_RCt_R1t_data.data ()};
77
+ // Target rotation * Control rotation transpose * Base rotation transpose
78
+ auto Rt_RcT_RbT_data = Kokkos::Array<double , 4 >{};
79
+ const auto Rt_RcT_RbT = Kokkos::View<double [4 ]>{Rt_RcT_RbT_data.data ()};
71
80
72
81
auto A_data = Kokkos::Array<double , 9 >{};
73
82
const auto A = View_3x3{A_data.data ()};
@@ -79,30 +88,38 @@ struct CalculateRotationControlConstraint {
79
88
const auto V3 = View_3{V3_data.data ()};
80
89
81
90
// ----------------------------------------------------------------------
82
- // Residual Vector
91
+ // Position residual
83
92
// ----------------------------------------------------------------------
84
93
85
- // Phi(0:3) = u2 + X0 - u1 - R1*X0
86
- QuaternionInverse (R1, R1t);
87
- RotateVectorByQuaternion (R1, X0, R1_X0);
94
+ // Phi(0:3) = ut + X0 - ub - Rb*X0
95
+ RotateVectorByQuaternion (Rb, X0, Rb_X0);
88
96
for (int i = 0 ; i < 3 ; ++i) {
89
- residual_terms (i_constraint, i) = u2 (i) + X0 (i) - u1 (i) - R1_X0 (i);
97
+ residual_terms (i_constraint, i) = ut (i) + X0 (i) - ub (i) - Rb_X0 (i);
90
98
}
91
99
92
- // Angular residual
93
- // If this is a rotation control constraint, calculate RC from control and axis
100
+ // ----------------------------------------------------------------------
101
+ // Rotation residual
102
+ // ----------------------------------------------------------------------
103
+
104
+ auto rotation_command = constraint_inputs (i_constraint, 0 );
105
+
106
+ // Copy rotation axis for this constraint
94
107
for (auto i = 0U ; i < 3U ; ++i) {
95
- RV (i) = axes (i_constraint, 0 , i) * constraint_inputs (i_constraint, 0 );
108
+ AX (i) = axes (i_constraint, 0 , i);
109
+ RV (i) = AX (i) * rotation_command;
96
110
}
97
- RotationVectorToQuaternion (RV, RC);
98
- QuaternionInverse (RC, RCt);
99
111
100
- // Phi(3:6) = axial(R2*inv(RC)*inv(R1))
101
- QuaternionCompose (R2, RCt, R2_RCt);
102
- QuaternionCompose (R2_RCt, R1t, R2_RCt_R1t);
103
- QuaternionToRotationMatrix (R2_RCt_R1t, C);
112
+ // Convert scaled axis to quaternion and calculate inverse
113
+ RotationVectorToQuaternion (RV, Rc);
114
+ QuaternionInverse (Rc, RcT);
115
+
116
+ // Phi(3:6) = axial(Rt*inv(Rc)*inv(Rb))
117
+ QuaternionInverse (Rb, RbT);
118
+ QuaternionCompose (Rt, RcT, Rt_RcT);
119
+ QuaternionCompose (Rt_RcT, RbT, Rt_RcT_RbT);
120
+ QuaternionToRotationMatrix (Rt_RcT_RbT, C);
104
121
AxialVectorOfMatrix (C, V3);
105
- for (int i = 0 ; i < 3 ; ++i) {
122
+ for (auto i = 0U ; i < 3U ; ++i) {
106
123
residual_terms (i_constraint, i + 3 ) = V3 (i);
107
124
}
108
125
@@ -119,7 +136,7 @@ struct CalculateRotationControlConstraint {
119
136
target_gradient_terms (i_constraint, i, i) = 1 .;
120
137
}
121
138
122
- // B(3:6,3:6) = AX(R1*RC *inv(R2 )) = transpose(AX(R2 *inv(RC )*inv(R1 )))
139
+ // B(3:6,3:6) = AX(Rb*Rc *inv(Rt )) = transpose(AX(Rt *inv(Rc )*inv(Rb )))
123
140
AX_Matrix (C, A);
124
141
for (int i = 0 ; i < 3 ; ++i) {
125
142
for (int j = 0 ; j < 3 ; ++j) {
@@ -136,15 +153,15 @@ struct CalculateRotationControlConstraint {
136
153
base_gradient_terms (i_constraint, i, i) = -1 .;
137
154
}
138
155
139
- // B(0:3,3:6) = tilde(R1 *X0)
140
- VecTilde (R1_X0 , A);
156
+ // B(0:3,3:6) = tilde(Rb *X0)
157
+ VecTilde (Rb_X0 , A);
141
158
for (int i = 0 ; i < 3 ; ++i) {
142
159
for (int j = 0 ; j < 3 ; ++j) {
143
160
base_gradient_terms (i_constraint, i, j + 3 ) = A (i, j);
144
161
}
145
162
}
146
163
147
- // B(3:6,3:6) = -AX(R2 *inv(RC )*inv(R1 ))
164
+ // B(3:6,3:6) = -AX(Rt *inv(Rc )*inv(Rb ))
148
165
AX_Matrix (C, A);
149
166
for (int i = 0 ; i < 3 ; ++i) {
150
167
for (int j = 0 ; j < 3 ; ++j) {
0 commit comments