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 : : #include <cblas.h>
19 : :
20 : : #include "Types.hpp"
21 : : #include "Exception.hpp"
22 : :
23 : : // ignore old-style-casts required for lapack/blas calls
24 : : #if defined(__clang__)
25 : : #pragma clang diagnostic ignored "-Wold-style-cast"
26 : : #endif
27 : :
28 : : // Lapacke forward declarations
29 : : extern "C" {
30 : :
31 : : using lapack_int = long;
32 : :
33 : : #define LAPACK_ROW_MAJOR 101
34 : : #define LAPACK_COL_MAJOR 102
35 : :
36 : : extern lapack_int LAPACKE_dgetrf( int, lapack_int, lapack_int, double*,
37 : : lapack_int, lapack_int* );
38 : : extern lapack_int LAPACKE_dgetri( int, lapack_int, double*, lapack_int,
39 : : const lapack_int* );
40 : :
41 : : }
42 : :
43 : : namespace tk {
44 : :
45 : : //! Flip sign of vector components
46 : : //! \param[in] v Vector whose components to multiply by -1.0
47 : : inline void
48 : : flip( std::array< real, 3 >& v )
49 : : {
50 : : v[0] = -v[0];
51 : : v[1] = -v[1];
52 : : v[2] = -v[2];
53 : : }
54 : :
55 : : //! Compute the cross-product of two vectors
56 : : //! \param[in] v1x x coordinate of the 1st vector
57 : : //! \param[in] v1y y coordinate of the 1st vector
58 : : //! \param[in] v1z z coordinate of the 1st vector
59 : : //! \param[in] v2x x coordinate of the 2nd vector
60 : : //! \param[in] v2y y coordinate of the 2nd vector
61 : : //! \param[in] v2z z coordinate of the 2nd vector
62 : : //! \param[out] rx x coordinate of the product vector
63 : : //! \param[out] ry y coordinate of the product vector
64 : : //! \param[out] rz z coordinate of the product vector
65 : : #pragma omp declare simd
66 : : inline void
67 : 1157859162 : cross( real v1x, real v1y, real v1z,
68 : : real v2x, real v2y, real v2z,
69 : : real& rx, real& ry, real& rz )
70 : : {
71 : 1157859162 : rx = v1y*v2z - v2y*v1z;
72 : 1157859162 : ry = v1z*v2x - v2z*v1x;
73 : 1157859162 : rz = v1x*v2y - v2x*v1y;
74 : 1157859162 : }
75 : :
76 : : //! Compute the cross-product of two vectors
77 : : //! \param[in] v1 1st vector
78 : : //! \param[in] v2 2nd vector
79 : : //! \return Cross-product
80 : : inline std::array< real, 3 >
81 : 830956914 : cross( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 )
82 : : {
83 : : real rx, ry, rz;
84 : 830956914 : cross( v1[0], v1[1], v1[2], v2[0], v2[1], v2[2], rx, ry, rz );
85 : 830956914 : return { std::move(rx), std::move(ry), std::move(rz) };
86 : : }
87 : :
88 : : //! Compute the cross-product of two vectors divided by a scalar
89 : : //! \param[in] v1x x coordinate of the 1st vector
90 : : //! \param[in] v1y y coordinate of the 1st vector
91 : : //! \param[in] v1z z coordinate of the 1st vector
92 : : //! \param[in] v2x x coordinate of the 2nd vector
93 : : //! \param[in] v2y y coordinate of the 2nd vector
94 : : //! \param[in] v2z z coordinate of the 2nd vector
95 : : //! \param[in] j The scalar to divide the product with
96 : : //! \param[out] rx x coordinate of the product vector
97 : : //! \param[out] ry y coordinate of the product vector
98 : : //! \param[out] rz z coordinate of the product vector
99 : : #pragma omp declare simd uniform(j)
100 : : inline void
101 : 221638818 : crossdiv( real v1x, real v1y, real v1z,
102 : : real v2x, real v2y, real v2z,
103 : : real j,
104 : : real& rx, real& ry, real& rz )
105 : : {
106 : 221638818 : cross( v1x, v1y, v1z, v2x, v2y, v2z, rx, ry, rz );
107 : 221638818 : rx /= j;
108 : 221638818 : ry /= j;
109 : 221638818 : rz /= j;
110 : 221638818 : }
111 : :
112 : : //! Compute the cross-product of two vectors divided by a scalar
113 : : //! \param[in] v1 1st vector
114 : : //! \param[in] v2 2nd vector
115 : : //! \param[in] j Scalar to divide each component by
116 : : //! \return Cross-product divided by scalar
117 : : inline std::array< real, 3 >
118 : 247848235 : crossdiv( const std::array< real, 3 >& v1,
119 : : const std::array< real, 3 >& v2,
120 : : real j )
121 : : {
122 : 247848235 : return {{ (v1[1]*v2[2] - v2[1]*v1[2]) / j,
123 : 247848235 : (v1[2]*v2[0] - v2[2]*v1[0]) / j,
124 : 495696470 : (v1[0]*v2[1] - v2[0]*v1[1]) / j }};
125 : : }
126 : :
127 : : //! Compute the dot-product of two vectors
128 : : //! \param[in] v1 1st vector
129 : : //! \param[in] v2 2nd vector
130 : : //! \return Dot-product
131 : : inline real
132 : 1879581379 : dot( const std::array< real, 3 >& v1, const std::array< real, 3 >& v2 )
133 : : {
134 : 1879581379 : return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
135 : : }
136 : :
137 : : //! Compute the dot-product of a matrix and a vector
138 : : //! \param[in] m Matrix
139 : : //! \param[in] v vector
140 : : //! \return Dot-product
141 : : inline std::array< real, 3 >
142 : 686000 : matvec(
143 : : const std::array< std::array< real, 3 >, 3 >& m,
144 : : const std::array< real, 3 >& v )
145 : : {
146 : 686000 : std::array< real, 3 > mv{0, 0, 0};
147 [ + + ]: 2744000 : for (std::size_t i=0; i<3; ++i) {
148 [ + + ]: 8232000 : for (std::size_t j=0; j<3; ++j)
149 : 6174000 : mv[i] += m[i][j]*v[j];
150 : : }
151 : :
152 : 686000 : return mv;
153 : : }
154 : :
155 : : //! Compute length of a vector
156 : : //! \param[in] x X coordinate of vector
157 : : //! \param[in] y Y coordinate of vector
158 : : //! \param[in] z Z coordinate of vector
159 : : //! \return length
160 : : #pragma omp declare simd
161 : : inline real
162 : 12673735 : length( real x, real y, real z )
163 : : {
164 : 12673735 : return std::sqrt( x*x + y*y + z*z );
165 : : }
166 : :
167 : : //! Compute length of a vector
168 : : //! \param[in] v vector
169 : : //! \return length
170 : : inline real
171 : 17881104 : length( const std::array< real, 3 >& v )
172 : : {
173 : 17881104 : return std::sqrt( dot(v,v) );
174 : : }
175 : :
176 : : //! Scale vector to unit length
177 : : //! \param[in,out] v Vector to normalize
178 : : inline void
179 : 2 : unit( std::array< real, 3 >& v ) noexcept(ndebug)
180 : : {
181 : 2 : auto len = length( v );
182 [ - + ][ - - ]: 2 : Assert( len > std::numeric_limits< tk::real >::epsilon(), "div by zero" );
[ - - ][ - - ]
183 : 2 : v[0] /= len;
184 : 2 : v[1] /= len;
185 : 2 : v[2] /= len;
186 : 2 : }
187 : :
188 : : //! Compute the triple-product of three vectors
189 : : //! \param[in] v1x x coordinate of the 1st vector
190 : : //! \param[in] v1y y coordinate of the 1st vector
191 : : //! \param[in] v1z z coordinate of the 1st vector
192 : : //! \param[in] v2x x coordinate of the 2nd vector
193 : : //! \param[in] v2y y coordinate of the 2nd vector
194 : : //! \param[in] v2z z coordinate of the 2nd vector
195 : : //! \param[in] v3x x coordinate of the 3rd vector
196 : : //! \param[in] v3y y coordinate of the 3rd vector
197 : : //! \param[in] v3z z coordinate of the 3rd vector
198 : : //! \return Scalar value of the triple product
199 : : #pragma omp declare simd
200 : : inline tk::real
201 : 84839598 : triple( real v1x, real v1y, real v1z,
202 : : real v2x, real v2y, real v2z,
203 : : real v3x, real v3y, real v3z )
204 : : {
205 : : real cx, cy, cz;
206 : 84839598 : cross( v2x, v2y, v2z, v3x, v3y, v3z, cx, cy, cz );
207 : 84839598 : return v1x*cx + v1y*cy + v1z*cz;
208 : : }
209 : :
210 : : //! Compute the triple-product of three vectors
211 : : //! \param[in] v1 1st vector
212 : : //! \param[in] v2 2nd vector
213 : : //! \param[in] v3 3rd vector
214 : : //! \return Triple-product
215 : : inline real
216 : 830956913 : triple( const std::array< real, 3 >& v1,
217 : : const std::array< real, 3 >& v2,
218 : : const std::array< real, 3 >& v3 )
219 : : {
220 : 830956913 : return dot( v1, cross(v2,v3) );
221 : : }
222 : :
223 : : //! Rotate vector about X axis
224 : : //! \param[in] v Vector to rotate
225 : : //! \param[in] angle Angle to use to rotate with
226 : : //! \return Rotated vector
227 : : inline std::array< real, 3 >
228 : 153116 : rotateX( const std::array< real, 3 >& v, real angle )
229 : : {
230 : : using std::cos; using std::sin;
231 : :
232 : : std::array< std::array< real, 3 >, 3 >
233 : 153116 : R{{ {{ 1.0, 0.0, 0.0 }},
234 : 153116 : {{ 0.0, cos(angle), -sin(angle) }},
235 : 153116 : {{ 0.0, sin(angle), cos(angle) }} }};
236 : :
237 : 153116 : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
238 : : }
239 : :
240 : : //! Rotate vector about Y axis
241 : : //! \param[in] v Vector to rotate
242 : : //! \param[in] angle Angle to use to rotate with
243 : : //! \return Rotated vector
244 : : inline std::array< real, 3 >
245 : 153116 : rotateY( const std::array< real, 3 >& v, real angle )
246 : : {
247 : : using std::cos; using std::sin;
248 : :
249 : : std::array< std::array< real, 3 >, 3 >
250 : 153116 : R{{ {{ cos(angle), 0.0, sin(angle) }},
251 : : {{ 0.0, 1.0, 0.0 }},
252 : 153116 : {{ -sin(angle), 0.0, cos(angle) }} }};
253 : :
254 : 153116 : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
255 : : }
256 : :
257 : : //! Rotate vector about Z axis
258 : : //! \param[in] v Vector to rotate
259 : : //! \param[in] angle Angle to use to rotate with
260 : : //! \return Rotated vector
261 : : inline std::array< real, 3 >
262 : 153116 : rotateZ( const std::array< real, 3 >& v, real angle )
263 : : {
264 : : using std::cos; using std::sin;
265 : :
266 : : std::array< std::array< real, 3 >, 3 >
267 : 153116 : R{{ {{ cos(angle), -sin(angle), 0.0 }},
268 : 153116 : {{ sin(angle), cos(angle), 0.0 }},
269 : 153116 : {{ 0.0, 0.0, 1.0 }} }};
270 : :
271 : 153116 : return {{ dot(R[0],v), dot(R[1],v), dot(R[2],v) }};
272 : : }
273 : :
274 : : //! \brief Compute the determinant of the Jacobian of a coordinate
275 : : //! transformation over a tetrahedron
276 : : //! \param[in] v1 (x,y,z) coordinates of 1st vertex of the tetrahedron
277 : : //! \param[in] v2 (x,y,z) coordinates of 2nd vertex of the tetrahedron
278 : : //! \param[in] v3 (x,y,z) coordinates of 3rd vertex of the tetrahedron
279 : : //! \param[in] v4 (x,y,z) coordinates of 4th vertex of the tetrahedron
280 : : //! \return Determinant of the Jacobian of transformation of physical
281 : : //! tetrahedron to reference (xi, eta, zeta) space
282 : : inline real
283 : 735272582 : Jacobian( const std::array< real, 3 >& v1,
284 : : const std::array< real, 3 >& v2,
285 : : const std::array< real, 3 >& v3,
286 : : const std::array< real, 3 >& v4 )
287 : : {
288 : 735272582 : std::array< real, 3 > ba{{ v2[0]-v1[0], v2[1]-v1[1], v2[2]-v1[2] }},
289 : 735272582 : ca{{ v3[0]-v1[0], v3[1]-v1[1], v3[2]-v1[2] }},
290 : 735272582 : da{{ v4[0]-v1[0], v4[1]-v1[1], v4[2]-v1[2] }};
291 : 735272582 : return triple( ba, ca, da );
292 : : }
293 : :
294 : : //! \brief Compute the inverse of the Jacobian of a coordinate transformation
295 : : //! over a tetrahedron
296 : : //! \param[in] v1 (x,y,z) coordinates of 1st vertex of the tetrahedron
297 : : //! \param[in] v2 (x,y,z) coordinates of 2nd vertex of the tetrahedron
298 : : //! \param[in] v3 (x,y,z) coordinates of 3rd vertex of the tetrahedron
299 : : //! \param[in] v4 (x,y,z) coordinates of 4th vertex of the tetrahedron
300 : : //! \return Inverse of the Jacobian of transformation of physical
301 : : //! tetrahedron to reference (xi, eta, zeta) space
302 : : inline std::array< std::array< real, 3 >, 3 >
303 : 11786039 : inverseJacobian( const std::array< real, 3 >& v1,
304 : : const std::array< real, 3 >& v2,
305 : : const std::array< real, 3 >& v3,
306 : : const std::array< real, 3 >& v4 )
307 : : {
308 : : std::array< std::array< real, 3 >, 3 > jacInv;
309 : :
310 : 11786039 : auto detJ = Jacobian( v1, v2, v3, v4 );
311 : :
312 : 23572078 : jacInv[0][0] = ( (v3[1]-v1[1])*(v4[2]-v1[2])
313 : 11786039 : - (v4[1]-v1[1])*(v3[2]-v1[2])) / detJ;
314 : 23572078 : jacInv[1][0] = -( (v2[1]-v1[1])*(v4[2]-v1[2])
315 : 11786039 : - (v4[1]-v1[1])*(v2[2]-v1[2])) / detJ;
316 : 23572078 : jacInv[2][0] = ( (v2[1]-v1[1])*(v3[2]-v1[2])
317 : 11786039 : - (v3[1]-v1[1])*(v2[2]-v1[2])) / detJ;
318 : :
319 : 23572078 : jacInv[0][1] = -( (v3[0]-v1[0])*(v4[2]-v1[2])
320 : 11786039 : - (v4[0]-v1[0])*(v3[2]-v1[2])) / detJ;
321 : 23572078 : jacInv[1][1] = ( (v2[0]-v1[0])*(v4[2]-v1[2])
322 : 11786039 : - (v4[0]-v1[0])*(v2[2]-v1[2])) / detJ;
323 : 23572078 : jacInv[2][1] = -( (v2[0]-v1[0])*(v3[2]-v1[2])
324 : 11786039 : - (v3[0]-v1[0])*(v2[2]-v1[2])) / detJ;
325 : :
326 : 23572078 : jacInv[0][2] = ( (v3[0]-v1[0])*(v4[1]-v1[1])
327 : 11786039 : - (v4[0]-v1[0])*(v3[1]-v1[1])) / detJ;
328 : 23572078 : jacInv[1][2] = -( (v2[0]-v1[0])*(v4[1]-v1[1])
329 : 11786039 : - (v4[0]-v1[0])*(v2[1]-v1[1])) / detJ;
330 : 23572078 : jacInv[2][2] = ( (v2[0]-v1[0])*(v3[1]-v1[1])
331 : 11786039 : - (v3[0]-v1[0])*(v2[1]-v1[1])) / detJ;
332 : :
333 : 11786039 : return jacInv;
334 : : }
335 : :
336 : : //! Compute the determinant of 3x3 matrix
337 : : //! \param[in] a 3x3 matrix
338 : : //! \return Determinant of the 3x3 matrix
339 : : inline tk::real
340 : 92185440 : determinant( const std::array< std::array< tk::real, 3 >, 3 >& a )
341 : : {
342 : 92185440 : return ( a[0][0] * (a[1][1]*a[2][2]-a[1][2]*a[2][1])
343 : 92185440 : - a[0][1] * (a[1][0]*a[2][2]-a[1][2]*a[2][0])
344 : 92185440 : + a[0][2] * (a[1][0]*a[2][1]-a[1][1]*a[2][0]) );
345 : : }
346 : :
347 : : //! Compute the inverse of 3x3 matrix
348 : : //! \param[in] m 3x3 matrix
349 : : //! \return Inverse of the 3x3 matrix
350 : : inline std::array< std::array< tk::real, 3 >, 3 >
351 : : inverse( const std::array< std::array< tk::real, 3 >, 3 >& m )
352 : : {
353 : : tk::real det = m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) -
354 : : m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) +
355 : : m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
356 : :
357 : : tk::real invdet = 1.0 / det;
358 : :
359 : : std::array< std::array< tk::real, 3 >, 3 > minv;
360 : : minv[0][0] = (m[1][1] * m[2][2] - m[2][1] * m[1][2]) * invdet;
361 : : minv[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) * invdet;
362 : : minv[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * invdet;
363 : : minv[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) * invdet;
364 : : minv[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * invdet;
365 : : minv[1][2] = (m[1][0] * m[0][2] - m[0][0] * m[1][2]) * invdet;
366 : : minv[2][0] = (m[1][0] * m[2][1] - m[2][0] * m[1][1]) * invdet;
367 : : minv[2][1] = (m[2][0] * m[0][1] - m[0][0] * m[2][1]) * invdet;
368 : : minv[2][2] = (m[0][0] * m[1][1] - m[1][0] * m[0][1]) * invdet;
369 : :
370 : : return minv;
371 : : }
372 : :
373 : : //! Solve a 3x3 system of equations using Cramer's rule
374 : : //! \param[in] a 3x3 lhs matrix
375 : : //! \param[in] b 3x1 rhs matrix
376 : : //! \return Array of solutions of the 3x3 system
377 : : inline std::array < tk::real, 3 >
378 : 23046360 : cramer( const std::array< std::array< tk::real, 3 >, 3>& a,
379 : : const std::array< tk::real, 3 >& b )
380 : : {
381 : 23046360 : auto de = determinant( a );
382 : :
383 : 23046360 : auto nu(0.0);
384 : : std::array < real, 3 > x;
385 : :
386 : 46092720 : nu = determinant( {{{{b[0], a[0][1], a[0][2]}},
387 : 23046360 : {{b[1], a[1][1], a[1][2]}},
388 : 23046360 : {{b[2], a[2][1], a[2][2]}}}} );
389 : 23046360 : x[0] = nu/de;
390 : :
391 : 46092720 : nu = determinant( {{{{a[0][0], b[0], a[0][2]}},
392 : 23046360 : {{a[1][0], b[1], a[1][2]}},
393 : 23046360 : {{a[2][0], b[2], a[2][2]}}}} );
394 : 23046360 : x[1] = nu/de;
395 : :
396 : 46092720 : nu = determinant( {{{{a[0][0], a[0][1], b[0]}},
397 : 23046360 : {{a[1][0], a[1][1], b[1]}},
398 : 23046360 : {{a[2][0], a[2][1], b[2]}}}} );
399 : 23046360 : x[2] = nu/de;
400 : :
401 : 23046360 : return x;
402 : : }
403 : :
404 : : //! Move a point to a reference space given coordinates of origin of that space
405 : : //! \param[in] origin Origin of reference frame to which point is to be moved
406 : : //! \param[in,out] point Point that needs to be moved to reference frame
407 : : inline void
408 : 89211 : movePoint( const std::array< tk::real, 3 >& origin,
409 : : std::array< tk::real, 3 >& point )
410 : : {
411 [ + + ]: 356844 : for (std::size_t i=0; i<3; ++i)
412 : 267633 : point[i] -= origin[i];
413 : 89211 : }
414 : :
415 : : //! Rotate a point in 3D space by specifying rotation angles in degrees
416 : : //! \param[in] angles Angles in 3D space by which point is to be rotated
417 : : //! \param[in,out] point Point that needs to be rotated
418 : : inline void
419 : 86855 : rotatePoint( const std::array< tk::real, 3 >& angles,
420 : : std::array< tk::real, 3 >& point )
421 : : {
422 : : // Convert angles to radian
423 : 86855 : tk::real pi = 4.0*std::atan(1.0);
424 : 86855 : auto a = angles[0] * pi/180.0;
425 : 86855 : auto b = angles[1] * pi/180.0;
426 : 86855 : auto c = angles[2] * pi/180.0;
427 : :
428 : : // Rotation matrix
429 : : std::array< std::array< tk::real, 3 >, 3 > rotMat;
430 : : {
431 : : using namespace std;
432 : 86855 : rotMat[0][0] = cos(b)*cos(c);
433 : 86855 : rotMat[0][1] = - cos(b)*sin(c);
434 : 86855 : rotMat[0][2] = sin(b);
435 : :
436 : 86855 : rotMat[1][0] = sin(a)*sin(b)*cos(c) + cos(a)*sin(c);
437 : 86855 : rotMat[1][1] = - sin(a)*sin(b)*sin(c) + cos(a)*cos(c);
438 : 86855 : rotMat[1][2] = - sin(a)*cos(b);
439 : :
440 : 86855 : rotMat[2][0] = - cos(a)*sin(b)*cos(c) + sin(a)*sin(c);
441 : 86855 : rotMat[2][1] = cos(a)*sin(b)*sin(c) + sin(a)*cos(c);
442 : 86855 : rotMat[2][2] = cos(a)*cos(b);
443 : : }
444 : :
445 : : // Apply rotation
446 : 86855 : std::array< tk::real, 3 > x{{0.0, 0.0, 0.0}};
447 [ + + ]: 347420 : for (std::size_t i=0; i<3; ++i) {
448 [ + + ]: 1042260 : for (std::size_t j=0; j<3; ++j) {
449 : 781695 : x[i] += rotMat[i][j]*point[j];
450 : : }
451 : : }
452 : 86855 : point = x;
453 : 86855 : }
454 : :
455 : : //! \brief Get the Right Cauchy-Green strain tensor from the inverse deformation
456 : : //! gradient tensor.
457 : : //! \param[in] g Inverse deformation gradient tensor
458 : : //! \return Right Cauchy-Green tensor
459 : : inline std::array< std::array< real, 3 >, 3 >
460 : 0 : getRightCauchyGreen(const std::array< std::array< real, 3 >, 3 >& g)
461 : : {
462 : : // allocate matrices
463 : : double G[9], C[9];
464 : :
465 : : // initialize c-matrices
466 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
467 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
468 : 0 : G[i*3+j] = g[i][j];
469 : : }
470 : :
471 : : // get g.g^T
472 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
473 : : 3, 3, 3, 1.0, G, 3, G, 3, 0.0, C, 3);
474 : :
475 : : // get inv(g.g^T)
476 : : lapack_int ipiv[3];
477 : :
478 : : #ifndef NDEBUG
479 : : lapack_int ierr =
480 : : #endif
481 [ - - ]: 0 : LAPACKE_dgetrf(LAPACK_ROW_MAJOR, 3, 3, C, 3, ipiv);
482 [ - - ][ - - ]: 0 : Assert(ierr==0, "Lapack error in LU factorization of g.g^T");
[ - - ][ - - ]
483 : :
484 : : #ifndef NDEBUG
485 : : lapack_int jerr =
486 : : #endif
487 [ - - ]: 0 : LAPACKE_dgetri(LAPACK_ROW_MAJOR, 3, C, 3, ipiv);
488 [ - - ][ - - ]: 0 : Assert(jerr==0, "Lapack error in inverting g.g^T");
[ - - ][ - - ]
489 : :
490 : : // Output C as 2D array
491 : 0 : return {{ {C[0], C[1], C[2]},
492 : 0 : {C[3], C[4], C[5]},
493 : 0 : {C[6], C[7], C[8]} }};
494 : : }
495 : :
496 : : //! \brief Get the Left Cauchy-Green strain tensor from the inverse deformation
497 : : //! gradient tensor.
498 : : //! \param[in] g Inverse deformation gradient tensor
499 : : //! \return Left Cauchy-Green tensor
500 : : inline std::array< std::array< real, 3 >, 3 >
501 : 0 : getLeftCauchyGreen(const std::array< std::array< real, 3 >, 3 >& g)
502 : : {
503 : : // allocate matrices
504 : : double G[9], b[9];
505 : :
506 : : // initialize c-matrices
507 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
508 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
509 : 0 : G[i*3+j] = g[i][j];
510 : : }
511 : :
512 : : // get g^T.g
513 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
514 : : 3, 3, 3, 1.0, G, 3, G, 3, 0.0, b, 3);
515 : :
516 : : // get inv(g^T.g)
517 : : lapack_int ipiv[3];
518 : :
519 : : #ifndef NDEBUG
520 : : lapack_int ierr =
521 : : #endif
522 [ - - ]: 0 : LAPACKE_dgetrf(LAPACK_ROW_MAJOR, 3, 3, b, 3, ipiv);
523 [ - - ][ - - ]: 0 : Assert(ierr==0, "Lapack error in LU factorization of g^T.g");
[ - - ][ - - ]
524 : :
525 : : #ifndef NDEBUG
526 : : lapack_int jerr =
527 : : #endif
528 [ - - ]: 0 : LAPACKE_dgetri(LAPACK_ROW_MAJOR, 3, b, 3, ipiv);
529 [ - - ][ - - ]: 0 : Assert(jerr==0, "Lapack error in inverting g^T.g");
[ - - ][ - - ]
530 : :
531 : : // Output b as 2D array
532 : 0 : return {{ {b[0], b[1], b[2]},
533 : 0 : {b[3], b[4], b[5]},
534 : 0 : {b[6], b[7], b[8]} }};
535 : : }
536 : :
537 : : //! \brief Get the volume-preserving part of the right Cauchy-Green strain
538 : : //! tensor from the inverse deformation gradient tensor.
539 : : //! \param[in] g Inverse deformation gradient tensor
540 : : //! \return Volume-preserving part of the right Cauchy-Green strain tensor
541 : : inline std::array< std::array< tk::real, 3 >, 3 >
542 : 0 : getIsochorRightCauchyGreen(const std::array< std::array< real, 3 >, 3 >& g)
543 : : {
544 : 0 : auto Ct = tk::getRightCauchyGreen(g);
545 : 0 : auto detC = std::pow(tk::determinant(Ct), 1.0/3.0);
546 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
547 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
548 : 0 : Ct[i][j] /= detC;
549 : : }
550 : :
551 : 0 : return Ct;
552 : : }
553 : :
554 : : //! \brief Get the deviatoric Hencky strain tensor from the inverse deformation
555 : : //! gradient tensor.
556 : : //! \param[in] g Inverse deformation gradient tensor
557 : : //! \return Deviatoric Hencky strain tensor
558 : : inline std::array< std::array< real, 3 >, 3 >
559 : 0 : getDevHencky(const std::array< std::array< real, 3 >, 3 >& g)
560 : : {
561 : : // Get right volm-preserving part of Cauchy-Green strain tensor
562 [ - - ]: 0 : auto C = getIsochorRightCauchyGreen(g);
563 : :
564 : 0 : std::array< std::array< real, 3 >, 3 > devH{{{0,0,0}, {0,0,0}, {0,0,0}}};
565 : :
566 : : // Compute approximation of Hencky strain tensor from section 3.4 of
567 : : // Barton, P. T. (2019). An interface-capturing Godunov method for the
568 : : // simulation of compressible solid-fluid problems. Journal of Computational
569 : : // Physics, 390, 25-50.
570 : :
571 : : // get g.g^T (i.e. inv(C))
572 : : double G[9], CInv[9];
573 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
574 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
575 : 0 : G[i*3+j] = g[i][j];
576 : : }
577 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
578 : : 3, 3, 3, 1.0, G, 3, G, 3, 0.0, CInv, 3);
579 : :
580 : : // volume-preserving part
581 : 0 : auto detCInv = std::pow(tk::determinant(
582 : 0 : {{ { CInv[0], CInv[1], CInv[2] },
583 : 0 : { CInv[3], CInv[4], CInv[5] },
584 : 0 : { CInv[6], CInv[7], CInv[8] } }} ), 1.0/3.0);
585 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
586 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
587 : 0 : CInv[3*i+j] /= detCInv;
588 : :
589 : : // Compute (C-CInv)/4
590 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
591 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
592 : 0 : devH[i][j] = 0.25*(C[i][j]-CInv[3*i+j]);
593 : :
594 : : // Compute deviatoric part
595 : 0 : tk::real trH = devH[0][0] + devH[1][1] + devH[2][2];
596 : :
597 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
598 : 0 : devH[i][i] -= trH/3.0;
599 : :
600 : : // Output devH
601 : 0 : return devH;
602 : : }
603 : :
604 : : //! \brief Get transpose of a 3 by 3 matrix
605 : : //! \param[in] mat matrix to be transposed
606 : : //! \return transposed matrix
607 : : inline std::array< std::array< tk::real, 3 >, 3 >
608 : 0 : transpose3by3(std::array< std::array< tk::real, 3 >, 3 > mat)
609 : : {
610 : : std::array< std::array< tk::real, 3 >, 3 > transMat;
611 [ - - ]: 0 : for (size_t i=0; i<3; ++i)
612 [ - - ]: 0 : for (size_t j=0; j<3; ++j)
613 : 0 : transMat[i][j] = mat[j][i];
614 : 0 : return transMat;
615 : : }
616 : :
617 : : //! \brief Get rotation matrix in 2D array form given a normal
618 : : //! direction. The remaining directions are given by section 5.3.1 of
619 : : //! Miller, G. H., and P. Colella. "A conservative three-dimensional
620 : : //! Eulerian method for coupled solid–fluid shock capturing."
621 : : //! Journal of Computational Physics 183.1 (2002): 26-82.
622 : : //! \param[in] r Coordinates of the first basis vector r = (rx,ry,rz).
623 : : //! \return rotation matrix
624 : : inline std::array< std::array< tk::real, 3 >, 3 >
625 : 1392147 : getRotationMatrix(const std::array< tk::real, 3 >& r)
626 : : {
627 : : std::array< std::array< tk::real, 3 >, 3 > rotMat;
628 : 1392147 : tk::real rx = r[0];
629 : 1392147 : tk::real ry = r[1];
630 : 1392147 : tk::real rz = r[2];
631 [ + + ]: 1392147 : if (std::abs(ry+rz) <= std::abs(ry-rz))
632 : : {
633 : 819426 : tk::real norm = std::sqrt(2*(1-rx*ry-rx*rz-ry*rz));
634 : 819426 : rotMat[0][0] = rx;
635 : 819426 : rotMat[0][1] = ry;
636 : 819426 : rotMat[0][2] = rz;
637 : 819426 : rotMat[1][0] = (ry-rz)/norm;
638 : 819426 : rotMat[1][1] = (rz-rx)/norm;
639 : 819426 : rotMat[1][2] = (rx-ry)/norm;
640 : 819426 : rotMat[2][0] = (rx*(ry+rz)-ry*ry-rz*rz)/norm;
641 : 819426 : rotMat[2][1] = (ry*(rx+rz)-rx*rx-rz*rz)/norm;
642 : 819426 : rotMat[2][2] = (rz*(rx+ry)-rx*rx-ry*ry)/norm;
643 : : }
644 : : else
645 : : {
646 : 572721 : tk::real norm = std::sqrt(2*(1+rz*(ry-rx)+rx*ry));
647 : 572721 : rotMat[0][0] = rx;
648 : 572721 : rotMat[0][1] = ry;
649 : 572721 : rotMat[0][2] = rz;
650 : 572721 : rotMat[1][0] = (ry+rz)/norm;
651 : 572721 : rotMat[1][1] = (rz-rx)/norm;
652 : 572721 : rotMat[1][2] = (-rx-ry)/norm;
653 : 572721 : rotMat[2][0] = (rx*(rz-ry)-ry*ry-rz*rz)/norm;
654 : 572721 : rotMat[2][1] = (ry*(rx+rz)+rx*rx+rz*rz)/norm;
655 : 572721 : rotMat[2][2] = (rz*(rx-ry)-rx*rx-ry*ry)/norm;
656 : : }
657 : 1392147 : return rotMat;
658 : : }
659 : :
660 : : //! \brief Rotate a second order tensor (e.g. a Strain/Stress matrix) from
661 : : //! the (x,y,z) to a new (r,s,t) coordinate system.
662 : : //! The directions are given by section 5.3.1 of
663 : : //! Miller, G. H., and P. Colella. "A conservative three-dimensional
664 : : //! Eulerian method for coupled solid–fluid shock capturing."
665 : : //! Journal of Computational Physics 183.1 (2002): 26-82.
666 : : //! \param[in] mat matrix to be rotated.
667 : : //! \param[in] r Coordinates of the first basis vector r = (rx,ry,rz).
668 : : //! \return rotated tensor
669 : : inline std::array< std::array< tk::real, 3 >, 3 >
670 : 1392147 : rotateTensor(const std::array< std::array< tk::real, 3 >, 3 >& mat,
671 : : const std::array< tk::real, 3 >& r )
672 : : {
673 : : // define rotation matrix
674 : 1392147 : std::array< std::array< tk::real, 3 >, 3 > rotMatrix = getRotationMatrix(r);
675 : : double rotMat[9];
676 : :
677 : : // define matrices
678 : : double matAuxIn[9], matAuxOut[9];
679 [ + + ]: 5568588 : for (std::size_t i=0; i<3; ++i)
680 [ + + ]: 16705764 : for (std::size_t j=0; j<3; ++j)
681 : : {
682 : 12529323 : matAuxIn[i*3+j] = mat[i][j];
683 : 12529323 : rotMat[i*3+j] = rotMatrix[i][j];
684 : : }
685 : :
686 : : // compute matAuxIn*rotMat and store it into matAuxOut
687 [ + - ]: 1392147 : cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
688 : : 3, 3, 3, 1.0, matAuxIn, 3, rotMat, 3, 0.0, matAuxOut, 3);
689 : :
690 : : // matAuxOut -> matAuxIn
691 [ + + ]: 13921470 : for (std::size_t i=0; i<9; i++)
692 : : {
693 : 12529323 : matAuxIn[i] = matAuxOut[i];
694 : 12529323 : matAuxOut[i] = 0.0;
695 : : }
696 : :
697 : : // compute rotMat^T*matAuxIn and store it into matAuxOut
698 [ + - ]: 1392147 : cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
699 : : 3, 3, 3, 1.0, rotMat, 3, matAuxIn, 3, 0.0, matAuxOut, 3);
700 : :
701 : : // return matAuxOut as a 2D array
702 : 1392147 : return {{ {matAuxOut[0], matAuxOut[1], matAuxOut[2]},
703 : 1392147 : {matAuxOut[3], matAuxOut[4], matAuxOut[5]},
704 : 1392147 : {matAuxOut[6], matAuxOut[7], matAuxOut[8]} }};
705 : : }
706 : :
707 : : //! \brief Rotate a second order tensor (e.g. a Strain/Stress matrix) from
708 : : //! the (x,y,z) to a new (r,s,t) coordinate system.
709 : : //! The directions are given by section 5.3.1 of
710 : : //! Miller, G. H., and P. Colella. "A conservative three-dimensional
711 : : //! Eulerian method for coupled solid–fluid shock capturing."
712 : : //! Journal of Computational Physics 183.1 (2002): 26-82.
713 : : //! \param[in] mat matrix to be rotated.
714 : : //! \param[in] r Coordinates of the first basis vector r = (rx,ry,rz).
715 : : //! \return rotated tensor
716 : : inline std::array< std::array< tk::real, 3 >, 3 >
717 : 0 : unrotateTensor(const std::array< std::array< tk::real, 3 >, 3 >& mat,
718 : : const std::array< tk::real, 3 >& r )
719 : : {
720 : : // define rotation matrix
721 : 0 : std::array< std::array< tk::real, 3 >, 3 > rotMatrix = getRotationMatrix(r);
722 : : double rotMat[9];
723 : :
724 : : // define matrices
725 : : double matAuxIn[9], matAuxOut[9];
726 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
727 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
728 : : {
729 : 0 : matAuxIn[i*3+j] = mat[i][j];
730 : 0 : rotMat[i*3+j] = rotMatrix[i][j];
731 : : }
732 : : // compute matAuxIn*rotMat and store it into matAuxOut
733 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
734 : : 3, 3, 3, 1.0, matAuxIn, 3, rotMat, 3, 0.0, matAuxOut, 3);
735 : : // matAuxOut -> matAuxIn
736 [ - - ]: 0 : for (std::size_t i=0; i<9; i++)
737 : : {
738 : 0 : matAuxIn[i] = matAuxOut[i];
739 : 0 : matAuxOut[i] = 0.0;
740 : : }
741 : : // compute rotMat^T*matAuxIn and store it into matAuxOut
742 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
743 : : 3, 3, 3, 1.0, rotMat, 3, matAuxIn, 3, 0.0, matAuxOut, 3);
744 : : // return matAuxOut as a 2D array
745 : 0 : return {{ {matAuxOut[0], matAuxOut[1], matAuxOut[2]},
746 : 0 : {matAuxOut[3], matAuxOut[4], matAuxOut[5]},
747 : 0 : {matAuxOut[6], matAuxOut[7], matAuxOut[8]} }};
748 : : }
749 : :
750 : : //! \brief Rotate a vector (e.g. a velocity) from
751 : : //! the (x,y,z) to a new (r,s,t) coordinate system.
752 : : //! The directions are given by section 5.3.1 of
753 : : //! Miller, G. H., and P. Colella. "A conservative three-dimensional
754 : : //! Eulerian method for coupled solid–fluid shock capturing."
755 : : //! Journal of Computational Physics 183.1 (2002): 26-82.
756 : : //! \param[in] v Vector to be rotated.
757 : : //! \param[in] r Coordinates of the first basis vector r = (rx,ry,rz).
758 : : //! \return rotated vector
759 : : inline std::array< tk::real, 3 >
760 : 0 : rotateVector( const std::array< tk::real, 3 >& v,
761 : : const std::array< tk::real, 3 >& r )
762 : : {
763 : : // define rotation matrix
764 : 0 : std::array< std::array< tk::real, 3 >, 3 > rotMat = getRotationMatrix(r);
765 : :
766 : : // return rotMat*v
767 : 0 : return matvec(rotMat,v);
768 : : }
769 : :
770 : : //! \brief Unrotate a vector (e.g. a velocity) from
771 : : //! the (x,y,z) to a new (r,s,t) coordinate system.
772 : : //! The directions are given by section 5.3.1 of
773 : : //! Miller, G. H., and P. Colella. "A conservative three-dimensional
774 : : //! Eulerian method for coupled solid–fluid shock capturing."
775 : : //! Journal of Computational Physics 183.1 (2002): 26-82.
776 : : //! \param[in] v Vector to be rotated.
777 : : //! \param[in] r Coordinates of the first basis vector r = (rx,ry,rz).
778 : : //! \return rotated vector
779 : : inline std::array< tk::real, 3 >
780 : 0 : unrotateVector( const std::array< tk::real, 3 >& v,
781 : : const std::array< tk::real, 3 >& r )
782 : : {
783 : : // define rotation matrix
784 : 0 : std::array< std::array< tk::real, 3 >, 3 > rotMat = getRotationMatrix(r);
785 : :
786 : : // return rotMat^T*v
787 : 0 : return matvec(transpose3by3(rotMat),v);
788 : : }
789 : :
790 : : //! \brief Reflect a second order tensor (e.g. a Strain/Stress matrix)
791 : : //! \param[in] mat matrix to be rotated.
792 : : //! \param[in] reflectMat Reflection matrix
793 : : //! \return reflected tensor
794 : : inline std::array< std::array< tk::real, 3 >, 3 >
795 : 0 : reflectTensor(const std::array< std::array< tk::real, 3 >, 3 >& mat,
796 : : const std::array< std::array< tk::real, 3 >, 3 >& reflectMat)
797 : : {
798 : : // define reflection matrix
799 : : double refMat[9];
800 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
801 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
802 : 0 : refMat[i*3+j] = reflectMat[i][j];
803 : :
804 : : // define matAux (I need matrices as row major 1D arrays)
805 : : double matAuxIn[9], matAuxOut[9];
806 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
807 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
808 : 0 : matAuxIn[i*3+j] = mat[i][j];
809 : :
810 : : // compute matAuxIn*refMat and store it into matAuxOut
811 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
812 : : 3, 3, 3, 1.0, matAuxIn, 3, refMat, 3, 0.0, matAuxOut, 3);
813 : :
814 : : // matAuxOut -> matAuxIn
815 [ - - ]: 0 : for (std::size_t i=0; i<9; i++)
816 : : {
817 : 0 : matAuxIn[i] = matAuxOut[i];
818 : 0 : matAuxOut[i] = 0.0;
819 : : }
820 : :
821 : : // compute refMat^T*matAuxIn and store it into matAuxOut
822 [ - - ]: 0 : cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
823 : : 3, 3, 3, 1.0, refMat, 3, matAuxIn, 3, 0.0, matAuxOut, 3);
824 : :
825 : : // return matAuxOut as a 2D array
826 : 0 : return {{ {matAuxOut[0], matAuxOut[1], matAuxOut[2]},
827 : 0 : {matAuxOut[3], matAuxOut[4], matAuxOut[5]},
828 : 0 : {matAuxOut[6], matAuxOut[7], matAuxOut[8]} }};
829 : : }
830 : :
831 : : } // tk::
832 : :
833 : : #endif // Vector_h
|
