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 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Vector algebra
9 : : \details Vector algebra.
10 : : */
11 : : // *****************************************************************************
12 : : #ifndef Vector_h
13 : : #define Vector_h
14 : :
15 : : #include <array>
16 : : #include <cmath>
17 : : #include <vector>
18 : :
19 : : #include "Types.hpp"
20 : : #include "Exception.hpp"
21 : :
22 : : namespace tk {
23 : :
24 : : //! Flip sign of vector components
25 : : //! \param[in] v Vector whose components to multiply by -1.0
26 : : inline void
27 : : flip( std::array< real, 3 >& v )
28 : : {
29 : : v[0] = -v[0];
30 : : v[1] = -v[1];
31 : : v[2] = -v[2];
32 : : }
33 : :
34 : : //! Compute the cross-product of two vectors
35 : : //! \param[in] v1x x coordinate of the 1st vector
36 : : //! \param[in] v1y y coordinate of the 1st vector
37 : : //! \param[in] v1z z coordinate of the 1st vector
38 : : //! \param[in] v2x x coordinate of the 2nd vector
39 : : //! \param[in] v2y y coordinate of the 2nd vector
40 : : //! \param[in] v2z z coordinate of the 2nd vector
41 : : //! \param[out] rx x coordinate of the product vector
42 : : //! \param[out] ry y coordinate of the product vector
43 : : //! \param[out] rz z coordinate of the product vector
44 : : #pragma omp declare simd
45 : : inline void
46 : : cross( real v1x, real v1y, real v1z,
47 : : real v2x, real v2y, real v2z,
48 : : real& rx, real& ry, real& rz )
49 : : {
50 : 1551775699 : rx = v1y*v2z - v2y*v1z;
51 : 1551775699 : ry = v1z*v2x - v2z*v1x;
52 : 1538914125 : rz = v1x*v2y - v2x*v1y;
53 : : }
54 : :
55 : : //! Compute the cross-product of two vectors
56 : : //! \param[in] v1 1st vector
57 : : //! \param[in] v2 2nd vector
58 : : //! \return Cross-product
59 : : inline std::array< real, 3 >
60 : : cross( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 )
61 : : {
62 : : real rx, ry, rz;
63 : : cross( v1[0], v1[1], v1[2], v2[0], v2[1], v2[2], rx, ry, rz );
64 : : return { std::move(rx), std::move(ry), std::move(rz) };
65 : : }
66 : :
67 : : //! Compute the cross-product of two vectors divided by a scalar
68 : : //! \param[in] v1x x coordinate of the 1st vector
69 : : //! \param[in] v1y y coordinate of the 1st vector
70 : : //! \param[in] v1z z coordinate of the 1st vector
71 : : //! \param[in] v2x x coordinate of the 2nd vector
72 : : //! \param[in] v2y y coordinate of the 2nd vector
73 : : //! \param[in] v2z z coordinate of the 2nd vector
74 : : //! \param[in] j The scalar to divide the product with
75 : : //! \param[out] rx x coordinate of the product vector
76 : : //! \param[out] ry y coordinate of the product vector
77 : : //! \param[out] rz z coordinate of the product vector
78 : : #pragma omp declare simd uniform(j)
79 : : inline void
80 : : crossdiv( real v1x, real v1y, real v1z,
81 : : real v2x, real v2y, real v2z,
82 : : real j,
83 : : real& rx, real& ry, real& rz )
84 : : {
85 : : cross( v1x, v1y, v1z, v2x, v2y, v2z, rx, ry, rz );
86 : 147809125 : rx /= j;
87 : 147809125 : ry /= j;
88 : 147809125 : rz /= j;
89 : : }
90 : :
91 : : //! Compute the cross-product of two vectors divided by a scalar
92 : : //! \param[in] v1 1st vector
93 : : //! \param[in] v2 2nd vector
94 : : //! \param[in] j Scalar to divide each component by
95 : : //! \return Cross-product divided by scalar
96 : : inline std::array< real, 3 >
97 : : crossdiv( const std::array< real, 3 >& v1,
98 : : const std::array< real, 3 >& v2,
99 : : real j )
100 : : {
101 : 299383700 : return {{ (v1[1]*v2[2] - v2[1]*v1[2]) / j,
102 : 299383700 : (v1[2]*v2[0] - v2[2]*v1[0]) / j,
103 : 299290428 : (v1[0]*v2[1] - v2[0]*v1[1]) / j }};
104 : : }
105 : :
106 : : //! Compute the dot-product of two vectors
107 : : //! \param[in] v1 1st vector
108 : : //! \param[in] v2 2nd vector
109 : : //! \return Dot-product
110 : : inline real
111 : : dot( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 )
112 : : {
113 [ + + ][ + + ]: 5528040762 : 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 : : #pragma omp declare simd
122 : : inline real
123 : : length( real x, real y, real z )
124 : : {
125 [ + + ][ - - ]: 23578029 : return std::sqrt( x*x + y*y + z*z );
[ + + ][ - - ]
[ + + ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
126 : : }
127 : :
128 : : //! Compute length of a vector
129 : : //! \param[in] v vector
130 : : //! \return length
131 : : inline real
132 : : length( const std::array< real, 3 >& v )
133 : : {
134 [ + + ][ + + ]: 27339491 : return std::sqrt( dot(v,v) );
[ - - ][ - - ]
[ - - ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
135 : : }
136 : :
137 : : //! Scale vector to unit length
138 : : //! \param[in,out] v Vector to normalize
139 : : inline void
140 : : unit( std::array< real, 3 >& v ) noexcept(ndebug)
141 : : {
142 : : auto len = length( v );
143 : : Assert( len > std::numeric_limits< tk::real >::epsilon(), "div by zero" );
144 : : v[0] /= len;
145 : : v[1] /= len;
146 : : v[2] /= len;
147 : : }
148 : :
149 : : //! Compute the triple-product of three vectors
150 : : //! \param[in] v1x x coordinate of the 1st vector
151 : : //! \param[in] v1y y coordinate of the 1st vector
152 : : //! \param[in] v1z z coordinate of the 1st vector
153 : : //! \param[in] v2x x coordinate of the 2nd vector
154 : : //! \param[in] v2y y coordinate of the 2nd vector
155 : : //! \param[in] v2z z coordinate of the 2nd vector
156 : : //! \param[in] v3x x coordinate of the 3rd vector
157 : : //! \param[in] v3y y coordinate of the 3rd vector
158 : : //! \param[in] v3z z coordinate of the 3rd vector
159 : : //! \return Scalar value of the triple product
160 : : #pragma omp declare simd
161 : : inline tk::real
162 : : triple( real v1x, real v1y, real v1z,
163 : : real v2x, real v2y, real v2z,
164 : : real v3x, real v3y, real v3z )
165 : : {
166 : : real cx, cy, cz;
167 : : cross( v2x, v2y, v2z, v3x, v3y, v3z, cx, cy, cz );
168 [ - + ][ - + ]: 159962398 : return v1x*cx + v1y*cy + v1z*cz;
[ - + ][ - - ]
[ - + ][ - + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
169 : : }
170 : :
171 : : //! Compute the triple-product of three vectors
172 : : //! \param[in] v1 1st vector
173 : : //! \param[in] v2 2nd vector
174 : : //! \param[in] v3 3rd vector
175 : : //! \return Triple-product
176 : : inline real
177 : : triple( const std::array< real, 3 >& v1,
178 : : const std::array< real, 3 >& v2,
179 : : const std::array< real, 3 >& v3 )
180 : : {
181 : : return dot( v1, cross(v2,v3) );
182 : : }
183 : :
184 : : //! Rotate vector about X axis
185 : : //! \param[in] v Vector to rotate
186 : : //! \param[in] angle Angle to use to rotate with
187 : : //! \return Rotated vector
188 : : inline std::array< real, 3 >
189 : 153116 : rotateX( const std::array< real, 3 >& v, real angle )
190 : : {
191 : : using std::cos; using std::sin;
192 : :
193 : : std::array< std::array< real, 3 >, 3 >
194 : : R{{ {{ 1.0, 0.0, 0.0 }},
195 : 153116 : {{ 0.0, cos(angle), -sin(angle) }},
196 : : {{ 0.0, sin(angle), cos(angle) }} }};
197 : :
198 : 153116 : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
199 : : }
200 : :
201 : : //! Rotate vector about Y axis
202 : : //! \param[in] v Vector to rotate
203 : : //! \param[in] angle Angle to use to rotate with
204 : : //! \return Rotated vector
205 : : inline std::array< real, 3 >
206 : 153116 : rotateY( const std::array< real, 3 >& v, real angle )
207 : : {
208 : : using std::cos; using std::sin;
209 : :
210 : : std::array< std::array< real, 3 >, 3 >
211 : 153116 : R{{ {{ cos(angle), 0.0, sin(angle) }},
212 : : {{ 0.0, 1.0, 0.0 }},
213 : 153116 : {{ -sin(angle), 0.0, cos(angle) }} }};
214 : :
215 : 153116 : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
216 : : }
217 : :
218 : : //! Rotate vector about Z axis
219 : : //! \param[in] v Vector to rotate
220 : : //! \param[in] angle Angle to use to rotate with
221 : : //! \return Rotated vector
222 : : inline std::array< real, 3 >
223 : 153116 : rotateZ( const std::array< real, 3 >& v, real angle )
224 : : {
225 : : using std::cos; using std::sin;
226 : :
227 : : std::array< std::array< real, 3 >, 3 >
228 : 153116 : R{{ {{ cos(angle), -sin(angle), 0.0 }},
229 : : {{ sin(angle), cos(angle), 0.0 }},
230 : : {{ 0.0, 0.0, 1.0 }} }};
231 : :
232 : 153116 : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
233 : : }
234 : :
235 : : //! \brief Compute the determinant of the Jacobian of a coordinate
236 : : //! transformation over a tetrahedron
237 : : //! \param[in] v1 (x,y,z) coordinates of 1st vertex of the tetrahedron
238 : : //! \param[in] v2 (x,y,z) coordinates of 2nd vertex of the tetrahedron
239 : : //! \param[in] v3 (x,y,z) coordinates of 3rd vertex of the tetrahedron
240 : : //! \param[in] v4 (x,y,z) coordinates of 4th vertex of the tetrahedron
241 : : //! \return Determinant of the Jacobian of transformation of physical
242 : : //! tetrahedron to reference (xi, eta, zeta) space
243 : : inline real
244 : 1027702892 : Jacobian( const std::array< real, 3 >& v1,
245 : : const std::array< real, 3 >& v2,
246 : : const std::array< real, 3 >& v3,
247 : : const std::array< real, 3 >& v4 )
248 : : {
249 : 1027702892 : std::array< real, 3 > ba{{ v2[0]-v1[0], v2[1]-v1[1], v2[2]-v1[2] }},
250 : 1027702892 : ca{{ v3[0]-v1[0], v3[1]-v1[1], v3[2]-v1[2] }},
251 : 1027702892 : da{{ v4[0]-v1[0], v4[1]-v1[1], v4[2]-v1[2] }};
252 : 1027702892 : return triple( ba, ca, da );
253 : : }
254 : :
255 : : //! \brief Compute the inverse of the Jacobian of a coordinate transformation
256 : : //! over a tetrahedron
257 : : //! \param[in] v1 (x,y,z) coordinates of 1st vertex of the tetrahedron
258 : : //! \param[in] v2 (x,y,z) coordinates of 2nd vertex of the tetrahedron
259 : : //! \param[in] v3 (x,y,z) coordinates of 3rd vertex of the tetrahedron
260 : : //! \param[in] v4 (x,y,z) coordinates of 4th vertex of the tetrahedron
261 : : //! \return Inverse of the Jacobian of transformation of physical
262 : : //! tetrahedron to reference (xi, eta, zeta) space
263 : : inline std::array< std::array< real, 3 >, 3 >
264 : 15370366 : inverseJacobian( const std::array< real, 3 >& v1,
265 : : const std::array< real, 3 >& v2,
266 : : const std::array< real, 3 >& v3,
267 : : const std::array< real, 3 >& v4 )
268 : : {
269 : : std::array< std::array< real, 3 >, 3 > jacInv;
270 : :
271 : 15370366 : auto detJ = Jacobian( v1, v2, v3, v4 );
272 : :
273 : 15370366 : jacInv[0][0] = ( (v3[1]-v1[1])*(v4[2]-v1[2])
274 : 15370366 : - (v4[1]-v1[1])*(v3[2]-v1[2])) / detJ;
275 : 15370366 : jacInv[1][0] = -( (v2[1]-v1[1])*(v4[2]-v1[2])
276 : 15370366 : - (v4[1]-v1[1])*(v2[2]-v1[2])) / detJ;
277 : 15370366 : jacInv[2][0] = ( (v2[1]-v1[1])*(v3[2]-v1[2])
278 : 15370366 : - (v3[1]-v1[1])*(v2[2]-v1[2])) / detJ;
279 : :
280 : 15370366 : jacInv[0][1] = -( (v3[0]-v1[0])*(v4[2]-v1[2])
281 : 15370366 : - (v4[0]-v1[0])*(v3[2]-v1[2])) / detJ;
282 : 15370366 : jacInv[1][1] = ( (v2[0]-v1[0])*(v4[2]-v1[2])
283 : 15370366 : - (v4[0]-v1[0])*(v2[2]-v1[2])) / detJ;
284 : 15370366 : jacInv[2][1] = -( (v2[0]-v1[0])*(v3[2]-v1[2])
285 : 15370366 : - (v3[0]-v1[0])*(v2[2]-v1[2])) / detJ;
286 : :
287 : 15370366 : jacInv[0][2] = ( (v3[0]-v1[0])*(v4[1]-v1[1])
288 : 15370366 : - (v4[0]-v1[0])*(v3[1]-v1[1])) / detJ;
289 : 15370366 : jacInv[1][2] = -( (v2[0]-v1[0])*(v4[1]-v1[1])
290 : 15370366 : - (v4[0]-v1[0])*(v2[1]-v1[1])) / detJ;
291 : 15370366 : jacInv[2][2] = ( (v2[0]-v1[0])*(v3[1]-v1[1])
292 : 15370366 : - (v3[0]-v1[0])*(v2[1]-v1[1])) / detJ;
293 : :
294 : 15370366 : return jacInv;
295 : : }
296 : :
297 : : //! Compute the determinant of 3x3 matrix
298 : : //! \param[in] a 3x3 matrix
299 : : //! \return Determinant of the 3x3 matrix
300 : : inline tk::real
301 : 85882800 : determinant( const std::array< std::array< tk::real, 3 >, 3 >& a )
302 : : {
303 : 85882800 : return ( a[0][0] * (a[1][1]*a[2][2]-a[1][2]*a[2][1])
304 : 85882800 : - a[0][1] * (a[1][0]*a[2][2]-a[1][2]*a[2][0])
305 : 85882800 : + a[0][2] * (a[1][0]*a[2][1]-a[1][1]*a[2][0]) );
306 : : }
307 : :
308 : : //! Solve a 3x3 system of equations using Cramer's rule
309 : : //! \param[in] a 3x3 lhs matrix
310 : : //! \param[in] b 3x1 rhs matrix
311 : : //! \return Array of solutions of the 3x3 system
312 : : inline std::array < tk::real, 3 >
313 : 21470700 : cramer( const std::array< std::array< tk::real, 3 >, 3>& a,
314 : : const std::array< tk::real, 3 >& b )
315 : : {
316 : 21470700 : auto de = determinant( a );
317 : :
318 : : auto nu(0.0);
319 : : std::array < real, 3 > x;
320 : :
321 : 21470700 : nu = determinant( {{{{b[0], a[0][1], a[0][2]}},
322 : : {{b[1], a[1][1], a[1][2]}},
323 : : {{b[2], a[2][1], a[2][2]}}}} );
324 : 21470700 : x[0] = nu/de;
325 : :
326 : 21470700 : nu = determinant( {{{{a[0][0], b[0], a[0][2]}},
327 : : {{a[1][0], b[1], a[1][2]}},
328 : : {{a[2][0], b[2], a[2][2]}}}} );
329 : 21470700 : x[1] = nu/de;
330 : :
331 : 21470700 : nu = determinant( {{{{a[0][0], a[0][1], b[0]}},
332 : : {{a[1][0], a[1][1], b[1]}},
333 : : {{a[2][0], a[2][1], b[2]}}}} );
334 : 21470700 : x[2] = nu/de;
335 : :
336 : 21470700 : return x;
337 : : }
338 : :
339 : : } // tk::
340 : :
341 : : #endif // Vector_h
|
