Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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 "KokkosThread.h"
13 : #include "KokkosScalar.h"
14 : #include "KokkosJaggedArray.h"
15 :
16 : #ifdef MOOSE_KOKKOS_SCOPE
17 : #include "KokkosADReal.h"
18 : #endif
19 :
20 : #include "MooseError.h"
21 : #include "MooseUtils.h"
22 :
23 : #include "libmesh/tensor_tools.h"
24 :
25 : namespace Moose::Kokkos
26 : {
27 :
28 : template <typename T>
29 : struct Vector3;
30 :
31 : using Real3 = Vector3<Real>;
32 : using ADReal3 = Vector3<ADReal>;
33 :
34 : struct Real33;
35 :
36 : template <typename T>
37 : struct Vector3
38 : {
39 : T v[3];
40 :
41 : #ifdef MOOSE_KOKKOS_SCOPE
42 : Vector3(const libMesh::TypeVector<T> & vector);
43 83968047 : KOKKOS_INLINE_FUNCTION Vector3() { *this = T{}; }
44 67866365 : KOKKOS_INLINE_FUNCTION Vector3(const T & scalar) { *this = scalar; }
45 555113461 : KOKKOS_INLINE_FUNCTION Vector3(const Vector3<T> & vector) { *this = vector; }
46 : KOKKOS_INLINE_FUNCTION Vector3(const T & x, const T & y, const T & z);
47 :
48 : KOKKOS_INLINE_FUNCTION Vector3<T> operator-() const;
49 136508 : KOKKOS_INLINE_FUNCTION T & operator()(unsigned int i) { return v[i]; }
50 47539154 : KOKKOS_INLINE_FUNCTION const T & operator()(unsigned int i) const { return v[i]; }
51 :
52 : Vector3<T> & operator=(const libMesh::TypeVector<T> & vector);
53 :
54 : template <typename U>
55 : KOKKOS_INLINE_FUNCTION Vector3<T> & operator=(const Vector3<U> & vector);
56 : KOKKOS_INLINE_FUNCTION Vector3<T> & operator=(const Vector3<T> & vector);
57 : KOKKOS_INLINE_FUNCTION Vector3<T> & operator=(const T & scalar);
58 :
59 : template <typename U>
60 : KOKKOS_INLINE_FUNCTION void operator+=(const Vector3<U> & vector);
61 : KOKKOS_INLINE_FUNCTION void operator+=(const T & scalar);
62 : template <typename U>
63 : KOKKOS_INLINE_FUNCTION void operator-=(const Vector3<U> & vector);
64 : KOKKOS_INLINE_FUNCTION void operator-=(const T & scalar);
65 : KOKKOS_INLINE_FUNCTION void operator*=(const T & scalar);
66 :
67 : KOKKOS_INLINE_FUNCTION Real norm() const;
68 : KOKKOS_INLINE_FUNCTION Real dot_product(const Real3 vector) const;
69 : KOKKOS_INLINE_FUNCTION Real3 cross_product(const Real3 vector) const;
70 : KOKKOS_INLINE_FUNCTION Real33 cartesian_product(const Real3 vector) const;
71 : #endif
72 : };
73 :
74 : struct Real33
75 : {
76 : Real a[3][3];
77 :
78 : #ifdef MOOSE_KOKKOS_SCOPE
79 17729931 : KOKKOS_INLINE_FUNCTION Real33() { *this = 0; }
80 : KOKKOS_INLINE_FUNCTION Real33(const Real scalar) { *this = scalar; }
81 512146284 : KOKKOS_INLINE_FUNCTION Real33(const Real33 & tensor) { *this = tensor; }
82 :
83 56694032 : KOKKOS_INLINE_FUNCTION Real & operator()(unsigned int i, unsigned int j) { return a[i][j]; }
84 3666976524 : KOKKOS_INLINE_FUNCTION Real operator()(unsigned int i, unsigned int j) const { return a[i][j]; }
85 :
86 : KOKKOS_INLINE_FUNCTION Real33 & operator=(const Real33 & tensor);
87 : KOKKOS_INLINE_FUNCTION Real33 & operator=(const Real scalar);
88 : KOKKOS_INLINE_FUNCTION void operator+=(const Real33 tensor);
89 :
90 : KOKKOS_INLINE_FUNCTION void identity(const unsigned int dim = 3);
91 : KOKKOS_INLINE_FUNCTION Real determinant(const unsigned int dim = 3) const;
92 : KOKKOS_INLINE_FUNCTION Real33 inverse(const unsigned int dim = 3) const;
93 : KOKKOS_INLINE_FUNCTION Real33 transpose() const;
94 : KOKKOS_INLINE_FUNCTION Real3 row(const unsigned int i) const;
95 : KOKKOS_INLINE_FUNCTION Real3 col(const unsigned int j) const;
96 : #endif
97 : };
98 :
99 : #ifdef MOOSE_KOKKOS_SCOPE
100 :
101 : template <typename T>
102 119 : Vector3<T>::Vector3(const libMesh::TypeVector<T> & vector)
103 : {
104 119 : v[0] = vector(0);
105 119 : v[1] = vector(1);
106 119 : v[2] = vector(2);
107 119 : }
108 :
109 : template <typename T>
110 : KOKKOS_INLINE_FUNCTION
111 798213884 : Vector3<T>::Vector3(const T & x, const T & y, const T & z)
112 : {
113 773467196 : v[0] = x;
114 773467196 : v[1] = y;
115 773467196 : v[2] = z;
116 773467196 : }
117 :
118 : template <typename T>
119 : Vector3<T> &
120 1071317 : Vector3<T>::operator=(const libMesh::TypeVector<T> & vector)
121 : {
122 1071317 : v[0] = vector(0);
123 1071317 : v[1] = vector(1);
124 1071317 : v[2] = vector(2);
125 :
126 1071317 : return *this;
127 : }
128 :
129 : template <typename T>
130 : KOKKOS_INLINE_FUNCTION Vector3<T>
131 2104 : Vector3<T>::operator-() const
132 : {
133 2104 : Vector3<T> vector(*this);
134 2104 : vector *= -1;
135 :
136 2104 : return vector;
137 : }
138 :
139 : template <typename T>
140 : KOKKOS_INLINE_FUNCTION Vector3<T> &
141 703560975 : Vector3<T>::operator=(const Vector3<T> & vector)
142 : {
143 703560975 : v[0] = vector.v[0];
144 703560975 : v[1] = vector.v[1];
145 703560975 : v[2] = vector.v[2];
146 :
147 703560975 : return *this;
148 : }
149 :
150 : template <typename T>
151 : template <typename U>
152 : KOKKOS_INLINE_FUNCTION Vector3<T> &
153 276560 : Vector3<T>::operator=(const Vector3<U> & vector)
154 : {
155 276560 : v[0] = vector.v[0];
156 276560 : v[1] = vector.v[1];
157 276560 : v[2] = vector.v[2];
158 :
159 276560 : return *this;
160 : }
161 :
162 : template <typename T>
163 : KOKKOS_INLINE_FUNCTION Vector3<T> &
164 137618348 : Vector3<T>::operator=(const T & scalar)
165 : {
166 137618348 : v[0] = scalar;
167 137618348 : v[1] = scalar;
168 137618348 : v[2] = scalar;
169 :
170 137618348 : return *this;
171 : }
172 :
173 : template <typename T>
174 : template <typename U>
175 : KOKKOS_INLINE_FUNCTION void
176 318424992 : Vector3<T>::operator+=(const Vector3<U> & vector)
177 : {
178 318424992 : v[0] += vector.v[0];
179 318424992 : v[1] += vector.v[1];
180 318424992 : v[2] += vector.v[2];
181 318424992 : }
182 :
183 : template <typename T>
184 : KOKKOS_INLINE_FUNCTION void
185 : Vector3<T>::operator+=(const T & scalar)
186 : {
187 : v[0] += scalar;
188 : v[1] += scalar;
189 : v[2] += scalar;
190 : }
191 :
192 : template <typename T>
193 : template <typename U>
194 : KOKKOS_INLINE_FUNCTION void
195 : Vector3<T>::operator-=(const Vector3<U> & vector)
196 : {
197 : v[0] -= vector.v[0];
198 : v[1] -= vector.v[1];
199 : v[2] -= vector.v[2];
200 : }
201 :
202 : template <typename T>
203 : KOKKOS_INLINE_FUNCTION void
204 : Vector3<T>::operator-=(const T & scalar)
205 : {
206 : v[0] -= scalar;
207 : v[1] -= scalar;
208 : v[2] -= scalar;
209 : }
210 :
211 : template <typename T>
212 : KOKKOS_INLINE_FUNCTION void
213 320744 : Vector3<T>::operator*=(const T & scalar)
214 : {
215 320744 : v[0] *= scalar;
216 320744 : v[1] *= scalar;
217 320744 : v[2] *= scalar;
218 320744 : }
219 :
220 : template <typename T>
221 : KOKKOS_INLINE_FUNCTION Vector3<T>
222 : operator+(const T & left, const Vector3<T> & right)
223 : {
224 : return {left + right.v[0], left + right.v[1], left + right.v[2]};
225 : }
226 :
227 : template <typename T>
228 : KOKKOS_INLINE_FUNCTION Vector3<T>
229 : operator+(const Vector3<T> & left, const T & right)
230 : {
231 : return {left.v[0] + right, left.v[1] + right, left.v[2] + right};
232 : }
233 :
234 : template <typename T>
235 : KOKKOS_INLINE_FUNCTION Vector3<T>
236 : operator+(const Vector3<T> & left, const Vector3<T> & right)
237 : {
238 : return {left.v[0] + right.v[0], left.v[1] + right.v[1], left.v[2] + right.v[2]};
239 : }
240 :
241 : template <typename T>
242 : KOKKOS_INLINE_FUNCTION Vector3<T>
243 : operator-(const T & left, const Vector3<T> & right)
244 : {
245 : return {left - right.v[0], left - right.v[1], left - right.v[2]};
246 : }
247 :
248 : template <typename T>
249 : KOKKOS_INLINE_FUNCTION Vector3<T>
250 : operator-(const Vector3<T> & left, const T & right)
251 : {
252 : return {left.v[0] - right, left.v[1] - right, left.v[2] - right};
253 : }
254 :
255 : template <typename T>
256 : KOKKOS_INLINE_FUNCTION Vector3<T>
257 0 : operator-(const Vector3<T> & left, const Vector3<T> & right)
258 : {
259 0 : return {left.v[0] - right.v[0], left.v[1] - right.v[1], left.v[2] - right.v[2]};
260 : }
261 :
262 : template <typename T>
263 : KOKKOS_INLINE_FUNCTION Vector3<T>
264 : operator*(const T & left, const Vector3<T> & right)
265 : {
266 : return {left * right.v[0], left * right.v[1], left * right.v[2]};
267 : }
268 :
269 : template <typename T>
270 : KOKKOS_INLINE_FUNCTION Vector3<T>
271 : operator*(const Vector3<T> & left, const T & right)
272 : {
273 : return {left.v[0] * right, left.v[1] * right, left.v[2] * right};
274 : }
275 :
276 : template <typename T>
277 : KOKKOS_INLINE_FUNCTION T
278 65843096 : operator*(const Vector3<T> & left, const Vector3<T> & right)
279 : {
280 65843096 : return left.v[0] * right.v[0] + left.v[1] * right.v[1] + left.v[2] * right.v[2];
281 : }
282 :
283 : template <>
284 : KOKKOS_INLINE_FUNCTION Real
285 318640 : Vector3<Real>::norm() const
286 : {
287 318640 : return ::Kokkos::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
288 : }
289 :
290 : template <>
291 : KOKKOS_INLINE_FUNCTION Real
292 : Vector3<Real>::dot_product(const Real3 vector) const
293 : {
294 : return v[0] * vector.v[0] + v[1] * vector.v[1] + v[2] * vector.v[2];
295 : }
296 :
297 : template <>
298 : KOKKOS_INLINE_FUNCTION Real3
299 633504 : Vector3<Real>::cross_product(const Real3 vector) const
300 : {
301 633504 : Real3 cross;
302 :
303 633504 : cross.v[0] = v[1] * vector.v[2] - v[2] * vector.v[1];
304 633504 : cross.v[1] = v[2] * vector.v[0] - v[0] * vector.v[2];
305 633504 : cross.v[2] = v[0] * vector.v[1] - v[1] * vector.v[0];
306 :
307 633504 : return cross;
308 : }
309 :
310 : template <>
311 : KOKKOS_INLINE_FUNCTION Real33
312 4601108 : Vector3<Real>::cartesian_product(const Real3 vector) const
313 : {
314 4601108 : Real33 tensor;
315 :
316 18404432 : for (unsigned int i = 0; i < 3; ++i)
317 55213296 : for (unsigned int j = 0; j < 3; ++j)
318 41409972 : tensor(i, j) = v[i] * vector.v[j];
319 :
320 4601108 : return tensor;
321 : }
322 :
323 : KOKKOS_INLINE_FUNCTION Real33 &
324 549146472 : Real33::operator=(const Real33 & tensor)
325 : {
326 2196585888 : for (unsigned int i = 0; i < 3; ++i)
327 6589757664 : for (unsigned int j = 0; j < 3; ++j)
328 4942318248 : a[i][j] = tensor.a[i][j];
329 :
330 549146472 : return *this;
331 : }
332 :
333 : KOKKOS_INLINE_FUNCTION Real33 &
334 18048571 : Real33::operator=(const Real scalar)
335 : {
336 72194284 : for (unsigned int i = 0; i < 3; ++i)
337 216582852 : for (unsigned int j = 0; j < 3; ++j)
338 162437139 : a[i][j] = scalar;
339 :
340 18048571 : return *this;
341 : }
342 :
343 : KOKKOS_INLINE_FUNCTION void
344 4601108 : Real33::operator+=(const Real33 tensor)
345 : {
346 18404432 : for (unsigned int i = 0; i < 3; ++i)
347 55213296 : for (unsigned int j = 0; j < 3; ++j)
348 41409972 : a[i][j] += tensor.a[i][j];
349 4601108 : }
350 :
351 : KOKKOS_INLINE_FUNCTION void
352 0 : Real33::identity(const unsigned int dim)
353 : {
354 0 : *this = 0;
355 :
356 0 : for (unsigned int i = 0; i < dim; ++i)
357 0 : a[i][i] = 1;
358 0 : }
359 :
360 : KOKKOS_INLINE_FUNCTION Real
361 1876312 : Real33::determinant(const unsigned int dim) const
362 : {
363 1876312 : Real det = 0;
364 :
365 1876312 : if (dim == 0)
366 1024 : det = 1;
367 1875288 : else if (dim == 1)
368 329640 : det = a[0][0];
369 1545648 : else if (dim == 2)
370 1514352 : det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
371 31296 : else if (dim == 3)
372 31296 : det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
373 31296 : a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
374 31296 : a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
375 :
376 1876312 : return det;
377 : }
378 :
379 : KOKKOS_INLINE_FUNCTION Real33
380 938156 : Real33::inverse(const unsigned int dim) const
381 : {
382 938156 : Real inv_det = 1.0 / determinant(dim);
383 938156 : Real33 inv_mat;
384 :
385 938156 : if (dim == 1)
386 : {
387 7388 : inv_mat(0, 0) = inv_det;
388 : }
389 930768 : else if (dim == 2)
390 : {
391 914256 : inv_mat(0, 0) = a[1][1] * inv_det;
392 914256 : inv_mat(0, 1) = -a[0][1] * inv_det;
393 914256 : inv_mat(1, 0) = -a[1][0] * inv_det;
394 914256 : inv_mat(1, 1) = a[0][0] * inv_det;
395 : }
396 16512 : else if (dim == 3)
397 : {
398 16512 : inv_mat(0, 0) = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) * inv_det;
399 16512 : inv_mat(0, 1) = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) * inv_det;
400 16512 : inv_mat(0, 2) = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * inv_det;
401 16512 : inv_mat(1, 0) = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) * inv_det;
402 16512 : inv_mat(1, 1) = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) * inv_det;
403 16512 : inv_mat(1, 2) = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * inv_det;
404 16512 : inv_mat(2, 0) = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) * inv_det;
405 16512 : inv_mat(2, 1) = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) * inv_det;
406 16512 : inv_mat(2, 2) = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * inv_det;
407 : }
408 :
409 938156 : return inv_mat;
410 : }
411 :
412 : KOKKOS_INLINE_FUNCTION Real33
413 318640 : Real33::transpose() const
414 : {
415 318640 : Real33 tr_mat;
416 :
417 1274560 : for (unsigned int i = 0; i < 3; ++i)
418 3823680 : for (unsigned int j = 0; j < 3; ++j)
419 2867760 : tr_mat(i, j) = a[j][i];
420 :
421 318640 : return tr_mat;
422 : }
423 :
424 : KOKKOS_INLINE_FUNCTION Real3
425 319344 : Real33::row(const unsigned int i) const
426 : {
427 319344 : return Real3(a[i][0], a[i][1], a[i][2]);
428 : }
429 :
430 : KOKKOS_INLINE_FUNCTION Real3
431 : Real33::col(const unsigned int j) const
432 : {
433 : return Real3(a[0][j], a[1][j], a[2][j]);
434 : }
435 :
436 : KOKKOS_INLINE_FUNCTION Real3
437 405529996 : operator*(const Real33 left, const Real3 right)
438 : {
439 405529996 : return {left(0, 0) * right.v[0] + left(0, 1) * right.v[1] + left(0, 2) * right.v[2],
440 405529996 : left(1, 0) * right.v[0] + left(1, 1) * right.v[1] + left(1, 2) * right.v[2],
441 405529996 : left(2, 0) * right.v[0] + left(2, 1) * right.v[1] + left(2, 2) * right.v[2]};
442 : }
443 :
444 : KOKKOS_INLINE_FUNCTION Real33
445 318640 : operator*(const Real33 left, const Real33 right)
446 : {
447 318640 : Real33 mul;
448 :
449 1274560 : for (unsigned int i = 0; i < 3; ++i)
450 3823680 : for (unsigned int j = 0; j < 3; ++j)
451 11471040 : for (unsigned int k = 0; k < 3; ++k)
452 8603280 : mul(i, j) += left(i, k) * right(k, j);
453 :
454 318640 : return mul;
455 : }
456 :
457 : KOKKOS_INLINE_FUNCTION Real3
458 : operator+(const Real left, const Real3 right)
459 : {
460 : return {left + right.v[0], left + right.v[1], left + right.v[2]};
461 : }
462 :
463 : KOKKOS_INLINE_FUNCTION Real3
464 : operator+(const Real3 left, const Real right)
465 : {
466 : return {left.v[0] + right, left.v[1] + right, left.v[2] + right};
467 : }
468 :
469 : KOKKOS_INLINE_FUNCTION Real3
470 : operator-(const Real left, const Real3 right)
471 : {
472 : return {left - right.v[0], left - right.v[1], left - right.v[2]};
473 : }
474 :
475 : KOKKOS_INLINE_FUNCTION Real3
476 : operator-(const Real3 left, const Real right)
477 : {
478 : return {left.v[0] - right, left.v[1] - right, left.v[2] - right};
479 : }
480 :
481 : KOKKOS_INLINE_FUNCTION Real3
482 359368960 : operator*(const Real left, const Real3 right)
483 : {
484 359368960 : return {left * right.v[0], left * right.v[1], left * right.v[2]};
485 : }
486 :
487 : KOKKOS_INLINE_FUNCTION Real3
488 : operator*(const Real3 left, const Real right)
489 : {
490 : return {left.v[0] * right, left.v[1] * right, left.v[2] * right};
491 : }
492 :
493 : template <typename T,
494 : typename = typename std::enable_if<
495 : std::is_same<typename std::decay<T>::type, ADReal>::value>::type>
496 : KOKKOS_INLINE_FUNCTION ADReal3
497 : operator*(const Real3 left, const T & right)
498 : {
499 : return {left(0) * right, left(1) * right, left(2) * right};
500 : }
501 :
502 : template <typename T,
503 : typename = typename std::enable_if<
504 : std::is_same<typename std::decay<T>::type, ADReal>::value>::type>
505 : KOKKOS_INLINE_FUNCTION ADReal3
506 8248896 : operator*(const T & left, const Real3 right)
507 : {
508 8248896 : return {left * right(0), left * right(1), left * right(2)};
509 : }
510 :
511 : KOKKOS_INLINE_FUNCTION ADReal
512 : operator*(const Real3 left, const ADReal3 & right)
513 : {
514 : return left(0) * right(0) + left(1) * right(1) + left(2) * right(2);
515 : }
516 :
517 : KOKKOS_INLINE_FUNCTION ADReal
518 2369344 : operator*(const ADReal3 & left, const Real3 right)
519 : {
520 2369344 : return left(0) * right(0) + left(1) * right(1) + left(2) * right(2);
521 : }
522 :
523 : #endif
524 :
525 : template <typename T1, typename T2>
526 : struct Pair
527 : {
528 : T1 first;
529 : T2 second;
530 :
531 : template <typename T3, typename T4>
532 0 : auto & operator=(const std::pair<T3, T4> pair)
533 : {
534 0 : first = pair.first;
535 0 : second = pair.second;
536 :
537 0 : return *this;
538 : }
539 : };
540 :
541 : template <typename T1, typename T2>
542 : bool
543 89094 : operator<(const Pair<T1, T2> & left, const Pair<T1, T2> & right)
544 : {
545 89094 : return std::make_pair(left.first, left.second) < std::make_pair(right.first, right.second);
546 : }
547 :
548 : } // namespace Moose::Kokkos
|
