Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Base/Vector.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Vector algebra
10 : : \details Vector algebra.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef Vector_h
14 : : #define Vector_h
15 : :
16 : : #include <array>
17 : : #include <cmath>
18 : : #include <vector>
19 : :
20 : : #include "Types.hpp"
21 : : #include "Exception.hpp"
22 : :
23 : : namespace tk {
24 : :
25 : : //! Compute the cross-product of two vectors
26 : : //! \param[in] v1x x coordinate of the 1st vector
27 : : //! \param[in] v1y y coordinate of the 1st vector
28 : : //! \param[in] v1z z coordinate of the 1st vector
29 : : //! \param[in] v2x x coordinate of the 2nd vector
30 : : //! \param[in] v2y y coordinate of the 2nd vector
31 : : //! \param[in] v2z z coordinate of the 2nd vector
32 : : //! \param[out] rx x coordinate of the product vector
33 : : //! \param[out] ry y coordinate of the product vector
34 : : //! \param[out] rz z coordinate of the product vector
35 : : inline void
36 : : cross( real v1x, real v1y, real v1z,
37 : : real v2x, real v2y, real v2z,
38 : : real& rx, real& ry, real& rz )
39 : : {
40 : 27756330 : rx = v1y*v2z - v2y*v1z;
41 : 27756330 : ry = v1z*v2x - v2z*v1x;
42 [ + + ]: 27662340 : rz = v1x*v2y - v2x*v1y;
43 : : }
44 : :
45 : : //! Compute the cross-product of two vectors
46 : : //! \param[in] v1 1st vector
47 : : //! \param[in] v2 2nd vector
48 : : //! \return Cross-product
49 : : inline std::array< real, 3 >
50 : : cross( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 ) {
51 : : real rx, ry, rz;
52 : : cross( v1[0], v1[1], v1[2], v2[0], v2[1], v2[2], rx, ry, rz );
53 : 1 : return { std::move(rx), std::move(ry), std::move(rz) };
54 : : }
55 : :
56 : : //! Compute the cross-product of two vectors divided by a scalar
57 : : //! \param[in] v1x x coordinate of the 1st vector
58 : : //! \param[in] v1y y coordinate of the 1st vector
59 : : //! \param[in] v1z z coordinate of the 1st vector
60 : : //! \param[in] v2x x coordinate of the 2nd vector
61 : : //! \param[in] v2y y coordinate of the 2nd vector
62 : : //! \param[in] v2z z coordinate of the 2nd vector
63 : : //! \param[in] j The scalar to divide the product with
64 : : //! \param[out] rx x coordinate of the product vector
65 : : //! \param[out] ry y coordinate of the product vector
66 : : //! \param[out] rz z coordinate of the product vector
67 : : inline void
68 : : crossdiv( real v1x, real v1y, real v1z,
69 : : real v2x, real v2y, real v2z,
70 : : real j,
71 : : real& rx, real& ry, real& rz )
72 : : {
73 : : cross( v1x, v1y, v1z, v2x, v2y, v2z, rx, ry, rz );
74 : 4711868 : rx /= j;
75 : 4711868 : ry /= j;
76 [ - + ]: 4711868 : rz /= j;
77 : : }
78 : :
79 : : //! Compute the cross-product of two vectors divided by a scalar
80 : : //! \param[in] v1 1st vector
81 : : //! \param[in] v2 2nd vector
82 : : //! \param[in] j Scalar to divide each component by
83 : : //! \return Cross-product divided by scalar
84 : : inline std::array< real, 3 >
85 : : crossdiv( const std::array< real, 3 >& v1,
86 : : const std::array< real, 3 >& v2,
87 : : real j )
88 : : {
89 : 6059068 : return {{ (v1[1]*v2[2] - v2[1]*v1[2]) / j,
90 : 6059068 : (v1[2]*v2[0] - v2[2]*v1[0]) / j,
91 : 6059069 : (v1[0]*v2[1] - v2[0]*v1[1]) / j }};
92 : : }
93 : :
94 : : //! Compute the dot-product of two vectors
95 : : //! \param[in] v1x x coordinate of 1st vector
96 : : //! \param[in] v1y y coordinate of 1st vector
97 : : //! \param[in] v1z z coordinate of 1st vector
98 : : //! \param[in] v2x x coordinate of 2nd vector
99 : : //! \param[in] v2y y coordinate of 2nd vector
100 : : //! \param[in] v2z z coordinate of 2ndt vector
101 : : //! \return Dot-product
102 : : inline real
103 : : dot( real v1x, real v1y, real v1z, real v2x, real v2y, real v2z ) {
104 : : return v1x*v2x + v1y*v2y + v1z*v2z;
105 : : }
106 : :
107 : : //! Compute the dot-product of two vectors
108 : : //! \param[in] v1 1st vector
109 : : //! \param[in] v2 2nd vector
110 : : //! \return Dot-product
111 : : inline real
112 : : dot( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 ) {
113 [ + + ][ - + ]: 13770336 : return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
114 : : }
115 : :
116 : : //! Compute length of a vector
117 : : //! \param[in] x X coordinate of vector
118 : : //! \param[in] y Y coordinate of vector
119 : : //! \param[in] z Z coordinate of vector
120 : : //! \return length
121 : : inline real
122 : : length( real x, real y, real z ) {
123 [ + + ][ + + ]: 75135980 : return std::sqrt( x*x + y*y + z*z );
124 : : }
125 : :
126 : : //! Compute length of a vector
127 : : //! \param[in] v vector
128 : : //! \return length
129 : : inline real
130 : : length( const std::array< real, 3 >& v ) {
131 : 239447 : return std::sqrt( dot(v,v) );
132 : : }
133 : :
134 : : //! Scale vector to unit length
135 : : //! \param[in,out] v Vector to normalize
136 : : inline void
137 : 2 : unit( std::array< real, 3 >& v ) noexcept(ndebug) {
138 : : auto len = length( v );
139 : : // cppcheck-suppress throwInNoexceptFunction
140 : : Assert( len > std::numeric_limits< tk::real >::epsilon(), "div by zero" );
141 : 2 : v[0] /= len;
142 : 2 : v[1] /= len;
143 : 2 : v[2] /= len;
144 : 2 : }
145 : :
146 : : //! Compute the triple-product of three vectors
147 : : //! \param[in] v1x x coordinate of the 1st vector
148 : : //! \param[in] v1y y coordinate of the 1st vector
149 : : //! \param[in] v1z z coordinate of the 1st vector
150 : : //! \param[in] v2x x coordinate of the 2nd vector
151 : : //! \param[in] v2y y coordinate of the 2nd vector
152 : : //! \param[in] v2z z coordinate of the 2nd vector
153 : : //! \param[in] v3x x coordinate of the 3rd vector
154 : : //! \param[in] v3y y coordinate of the 3rd vector
155 : : //! \param[in] v3z z coordinate of the 3rd vector
156 : : //! \return Scalar value of the triple product
157 : : inline tk::real
158 : : triple( real v1x, real v1y, real v1z,
159 : : real v2x, real v2y, real v2z,
160 : : real v3x, real v3y, real v3z )
161 : : {
162 : : real cx, cy, cz;
163 : : cross( v2x, v2y, v2z, v3x, v3y, v3z, cx, cy, cz );
164 : : return v1x*cx + v1y*cy + v1z*cz;
165 : : }
166 : :
167 : : //! Compute the triple-product of three vectors
168 : : //! \param[in] v1 1st vector
169 : : //! \param[in] v2 2nd vector
170 : : //! \param[in] v3 3rd vector
171 : : //! \return Triple-product
172 : : inline real
173 : : triple( const std::array< real, 3 >& v1,
174 : : const std::array< real, 3 >& v2,
175 : : const std::array< real, 3 >& v3 )
176 : : {
177 : : return dot( v1, cross(v2,v3) );
178 : : }
179 : :
180 : : //! Rotate vector about X axis
181 : : //! \param[in] v Vector to rotate
182 : : //! \param[in] angle Angle to use to rotate with
183 : : //! \return Rotated vector
184 : : inline std::array< real, 3 >
185 : : rotateX( const std::array< real, 3 >& v, real angle )
186 : : {
187 : : using std::cos; using std::sin;
188 : :
189 : : std::array< std::array< real, 3 >, 3 >
190 : : R{{ {{ 1.0, 0.0, 0.0 }},
191 : : {{ 0.0, cos(angle), -sin(angle) }},
192 : : {{ 0.0, sin(angle), cos(angle) }} }};
193 : :
194 : : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
195 : : }
196 : :
197 : : //! Rotate vector about Y axis
198 : : //! \param[in] v Vector to rotate
199 : : //! \param[in] angle Angle to use to rotate with
200 : : //! \return Rotated vector
201 : : inline std::array< real, 3 >
202 : : rotateY( const std::array< real, 3 >& v, real angle )
203 : : {
204 : : using std::cos; using std::sin;
205 : :
206 : : std::array< std::array< real, 3 >, 3 >
207 : : R{{ {{ cos(angle), 0.0, sin(angle) }},
208 : : {{ 0.0, 1.0, 0.0 }},
209 : : {{ -sin(angle), 0.0, cos(angle) }} }};
210 : :
211 : : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
212 : : }
213 : :
214 : : //! Rotate vector about Z axis
215 : : //! \param[in] v Vector to rotate
216 : : //! \param[in] angle Angle to use to rotate with
217 : : //! \return Rotated vector
218 : : inline std::array< real, 3 >
219 : : rotateZ( const std::array< real, 3 >& v, real angle )
220 : : {
221 : : using std::cos; using std::sin;
222 : :
223 : : std::array< std::array< real, 3 >, 3 >
224 : : R{{ {{ cos(angle), -sin(angle), 0.0 }},
225 : : {{ sin(angle), cos(angle), 0.0 }},
226 : : {{ 0.0, 0.0, 1.0 }} }};
227 : :
228 : : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
229 : : }
230 : :
231 : : //! Compute the unit normal vector of a triangle
232 : : //! \param[in] x1 x coordinate of the 1st vertex of the triangle
233 : : //! \param[in] x2 x coordinate of the 2nd vertex of the triangle
234 : : //! \param[in] x3 x coordinate of the 3rd vertex of the triangle
235 : : //! \param[in] y1 y coordinate of the 1st vertex of the triangle
236 : : //! \param[in] y2 y coordinate of the 2nd vertex of the triangle
237 : : //! \param[in] y3 y coordinate of the 3rd vertex of the triangle
238 : : //! \param[in] z1 z coordinate of the 1st vertex of the triangle
239 : : //! \param[in] z2 z coordinate of the 2nd vertex of the triangle
240 : : //! \param[in] z3 z coordinate of the 3rd vertex of the triangle
241 : : //! \param[out] nx x coordinate of the unit normal
242 : : //! \param[out] ny y coordinate of the unit normal
243 : : //! \param[out] nz z coordinate of the unit normal
244 : : //! \return Triangle area
245 : : inline real
246 : 0 : normal( real x1, real x2, real x3,
247 : : real y1, real y2, real y3,
248 : : real z1, real z2, real z3,
249 : : real& nx, real& ny, real& nz )
250 : : {
251 : 0 : real ax = x2 - x1;
252 : 0 : real ay = y2 - y1;
253 : 0 : real az = z2 - z1;
254 : :
255 : 0 : real bx = x3 - x1;
256 : 0 : real by = y3 - y1;
257 : 0 : real bz = z3 - z1;
258 : :
259 : 0 : real n1 = ay*bz - az*by;
260 : 0 : real n2 = -(ax*bz - az*bx);
261 : 0 : real n3 = ax*by - ay*bx;
262 : :
263 : 0 : auto farea = std::sqrt( n1*n1 + n2*n2 + n3*n3 );
264 : :
265 : 0 : nx = n1/farea;
266 : 0 : ny = n2/farea;
267 : 0 : nz = n3/farea;
268 : :
269 : 0 : return farea;
270 : : }
271 : :
272 : : } // tk::
273 : :
274 : : #endif // Vector_h
|
