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 : #include "MooseError.h"
17 : #include "MooseUtils.h"
18 :
19 : #include "libmesh/tensor_tools.h"
20 :
21 : namespace Moose
22 : {
23 : namespace Kokkos
24 : {
25 :
26 : struct Real33
27 : {
28 : Real a[3][3];
29 :
30 : #ifdef MOOSE_KOKKOS_SCOPE
31 11805453 : KOKKOS_INLINE_FUNCTION Real33() { *this = 0; }
32 : KOKKOS_INLINE_FUNCTION Real33(const Real & scalar) { *this = scalar; }
33 400573540 : KOKKOS_INLINE_FUNCTION Real33(const Real33 & tensor) { *this = tensor; }
34 26590540 : KOKKOS_INLINE_FUNCTION Real & operator()(unsigned int i, unsigned int j) { return a[i][j]; }
35 2954174256 : KOKKOS_INLINE_FUNCTION Real operator()(unsigned int i, unsigned int j) const { return a[i][j]; }
36 421133300 : KOKKOS_INLINE_FUNCTION Real33 & operator=(const Real33 & tensor)
37 : {
38 1684533200 : for (unsigned int i = 0; i < 3; ++i)
39 5053599600 : for (unsigned int j = 0; j < 3; ++j)
40 3790199700 : a[i][j] = tensor.a[i][j];
41 :
42 421133300 : return *this;
43 : }
44 11805453 : KOKKOS_INLINE_FUNCTION Real33 & operator=(const Real scalar)
45 : {
46 47221812 : for (unsigned int i = 0; i < 3; ++i)
47 141665436 : for (unsigned int j = 0; j < 3; ++j)
48 106249077 : a[i][j] = scalar;
49 :
50 11805453 : return *this;
51 : }
52 1926282 : KOKKOS_INLINE_FUNCTION void operator+=(const Real33 tensor)
53 : {
54 7705128 : for (unsigned int i = 0; i < 3; ++i)
55 23115384 : for (unsigned int j = 0; j < 3; ++j)
56 17336538 : a[i][j] += tensor.a[i][j];
57 1926282 : }
58 0 : KOKKOS_INLINE_FUNCTION void identity(const unsigned int dim = 3)
59 : {
60 0 : *this = 0;
61 :
62 0 : for (unsigned int i = 0; i < dim; ++i)
63 0 : a[i][i] = 1;
64 0 : }
65 716240 : KOKKOS_INLINE_FUNCTION Real determinant(const unsigned int dim = 3)
66 : {
67 716240 : Real det = 0;
68 :
69 716240 : if (dim == 0)
70 526 : det = 1;
71 715714 : else if (dim == 1)
72 228040 : det = a[0][0];
73 487674 : else if (dim == 2)
74 466170 : det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
75 21504 : else if (dim == 3)
76 21504 : det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
77 21504 : a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
78 21504 : a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
79 :
80 716240 : return det;
81 : }
82 358120 : KOKKOS_INLINE_FUNCTION Real33 inverse(const unsigned int dim = 3)
83 : {
84 358120 : Real inv_det = 1.0 / determinant(dim);
85 358120 : Real33 inv_mat;
86 :
87 358120 : if (dim == 1)
88 : {
89 6394 : inv_mat(0, 0) = inv_det;
90 : }
91 351726 : else if (dim == 2)
92 : {
93 340974 : inv_mat(0, 0) = a[1][1] * inv_det;
94 340974 : inv_mat(0, 1) = -a[0][1] * inv_det;
95 340974 : inv_mat(1, 0) = -a[1][0] * inv_det;
96 340974 : inv_mat(1, 1) = a[0][0] * inv_det;
97 : }
98 10752 : else if (dim == 3)
99 : {
100 10752 : inv_mat(0, 0) = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) * inv_det;
101 10752 : inv_mat(0, 1) = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) * inv_det;
102 10752 : inv_mat(0, 2) = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * inv_det;
103 10752 : inv_mat(1, 0) = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) * inv_det;
104 10752 : inv_mat(1, 1) = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) * inv_det;
105 10752 : inv_mat(1, 2) = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * inv_det;
106 10752 : inv_mat(2, 0) = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) * inv_det;
107 10752 : inv_mat(2, 1) = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) * inv_det;
108 10752 : inv_mat(2, 2) = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * inv_det;
109 : }
110 :
111 358120 : return inv_mat;
112 : }
113 216304 : KOKKOS_INLINE_FUNCTION Real33 transpose()
114 : {
115 216304 : Real33 tr_mat;
116 :
117 865216 : for (unsigned int i = 0; i < 3; ++i)
118 2595648 : for (unsigned int j = 0; j < 3; ++j)
119 1946736 : tr_mat(i, j) = a[j][i];
120 :
121 216304 : return tr_mat;
122 : }
123 : #endif
124 : };
125 :
126 : struct Real3
127 : {
128 : Real v[3];
129 :
130 : #ifdef MOOSE_KOKKOS_SCOPE
131 9110471 : KOKKOS_INLINE_FUNCTION Real3()
132 9110471 : {
133 9110471 : v[0] = 0;
134 9110471 : v[1] = 0;
135 9110471 : v[2] = 0;
136 9110471 : }
137 53211839 : KOKKOS_INLINE_FUNCTION Real3(const Real & scalar)
138 53211839 : {
139 53211839 : v[0] = scalar;
140 53211839 : v[1] = scalar;
141 53211839 : v[2] = scalar;
142 53211839 : }
143 428217896 : KOKKOS_INLINE_FUNCTION Real3(const Real3 & vector)
144 428217896 : {
145 428217896 : v[0] = vector.v[0];
146 428217896 : v[1] = vector.v[1];
147 428217896 : v[2] = vector.v[2];
148 428217896 : }
149 622009528 : KOKKOS_INLINE_FUNCTION Real3(const Real & x, const Real & y, const Real & z)
150 622009528 : {
151 622009528 : v[0] = x;
152 622009528 : v[1] = y;
153 622009528 : v[2] = z;
154 622009528 : }
155 12 : Real3(const libMesh::TypeVector<Real> & vector)
156 12 : {
157 12 : v[0] = vector(0);
158 12 : v[1] = vector(1);
159 12 : v[2] = vector(2);
160 12 : }
161 :
162 166064 : KOKKOS_INLINE_FUNCTION Real & operator()(unsigned int i) { return v[i]; }
163 3080 : KOKKOS_INLINE_FUNCTION Real operator()(unsigned int i) const { return v[i]; }
164 :
165 74185316 : KOKKOS_INLINE_FUNCTION Real3 & operator=(const Real3 & vector)
166 : {
167 74185316 : v[0] = vector.v[0];
168 74185316 : v[1] = vector.v[1];
169 74185316 : v[2] = vector.v[2];
170 :
171 74185316 : return *this;
172 : }
173 : KOKKOS_INLINE_FUNCTION Real3 & operator=(const Real scalar)
174 : {
175 : v[0] = scalar;
176 : v[1] = scalar;
177 : v[2] = scalar;
178 :
179 : return *this;
180 : }
181 223208 : Real3 & operator=(const libMesh::TypeVector<Real> & vector)
182 : {
183 223208 : v[0] = vector(0);
184 223208 : v[1] = vector(1);
185 223208 : v[2] = vector(2);
186 :
187 223208 : return *this;
188 : }
189 : KOKKOS_INLINE_FUNCTION void operator+=(const Real scalar)
190 : {
191 : v[0] += scalar;
192 : v[1] += scalar;
193 : v[2] += scalar;
194 : }
195 248860376 : KOKKOS_INLINE_FUNCTION void operator+=(const Real3 vector)
196 : {
197 248860376 : v[0] += vector.v[0];
198 248860376 : v[1] += vector.v[1];
199 248860376 : v[2] += vector.v[2];
200 248860376 : }
201 : KOKKOS_INLINE_FUNCTION void operator-=(const Real scalar)
202 : {
203 : v[0] -= scalar;
204 : v[1] -= scalar;
205 : v[2] -= scalar;
206 : }
207 : KOKKOS_INLINE_FUNCTION void operator-=(const Real3 vector)
208 : {
209 : v[0] -= vector.v[0];
210 : v[1] -= vector.v[1];
211 : v[2] -= vector.v[2];
212 : }
213 : KOKKOS_INLINE_FUNCTION void operator*=(const Real scalar)
214 : {
215 : v[0] *= scalar;
216 : v[1] *= scalar;
217 : v[2] *= scalar;
218 : }
219 0 : KOKKOS_INLINE_FUNCTION Real norm() { return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); }
220 : KOKKOS_INLINE_FUNCTION Real dot_product(const Real3 vector)
221 : {
222 : return v[0] * vector.v[0] + v[1] * vector.v[1] + v[2] * vector.v[2];
223 : }
224 0 : KOKKOS_INLINE_FUNCTION Real3 cross_product(const Real3 vector)
225 : {
226 0 : Real3 cross;
227 :
228 0 : cross.v[0] = v[1] * vector.v[2] - v[2] * vector.v[1];
229 0 : cross.v[1] = v[2] * vector.v[0] - v[0] * vector.v[2];
230 0 : cross.v[2] = v[0] * vector.v[1] - v[1] * vector.v[0];
231 :
232 0 : return cross;
233 : }
234 1926282 : KOKKOS_INLINE_FUNCTION Real33 cartesian_product(const Real3 vector)
235 : {
236 1926282 : Real33 tensor;
237 :
238 7705128 : for (unsigned int i = 0; i < 3; ++i)
239 23115384 : for (unsigned int j = 0; j < 3; ++j)
240 17336538 : tensor(i, j) = v[i] * vector.v[j];
241 :
242 1926282 : return tensor;
243 : }
244 : #endif
245 : };
246 :
247 : #ifdef MOOSE_KOKKOS_SCOPE
248 : KOKKOS_INLINE_FUNCTION Real3
249 295065768 : operator*(const Real left, const Real3 right)
250 : {
251 295065768 : return {left * right.v[0], left * right.v[1], left * right.v[2]};
252 : }
253 : KOKKOS_INLINE_FUNCTION Real3
254 : operator*(const Real3 left, const Real right)
255 : {
256 : return {left.v[0] * right, left.v[1] * right, left.v[2] * right};
257 : }
258 : KOKKOS_INLINE_FUNCTION Real
259 57182000 : operator*(const Real3 left, const Real3 right)
260 : {
261 57182000 : return left.v[0] * right.v[0] + left.v[1] * right.v[1] + left.v[2] * right.v[2];
262 : }
263 : KOKKOS_INLINE_FUNCTION Real3
264 326943760 : operator*(const Real33 left, const Real3 right)
265 : {
266 326943760 : return {left(0, 0) * right.v[0] + left(0, 1) * right.v[1] + left(0, 2) * right.v[2],
267 326943760 : left(1, 0) * right.v[0] + left(1, 1) * right.v[1] + left(1, 2) * right.v[2],
268 326943760 : left(2, 0) * right.v[0] + left(2, 1) * right.v[1] + left(2, 2) * right.v[2]};
269 : }
270 : KOKKOS_INLINE_FUNCTION Real33
271 216304 : operator*(const Real33 left, const Real33 right)
272 : {
273 216304 : Real33 mul;
274 :
275 865216 : for (unsigned int i = 0; i < 3; ++i)
276 2595648 : for (unsigned int j = 0; j < 3; ++j)
277 7786944 : for (unsigned int k = 0; k < 3; ++k)
278 5840208 : mul(i, j) += left(i, k) * right(k, j);
279 :
280 216304 : return mul;
281 : }
282 : KOKKOS_INLINE_FUNCTION Real3
283 : operator+(const Real left, const Real3 right)
284 : {
285 : return {left + right.v[0], left + right.v[1], left + right.v[2]};
286 : }
287 : KOKKOS_INLINE_FUNCTION Real3
288 : operator+(const Real3 left, const Real right)
289 : {
290 : return {left.v[0] + right, left.v[1] + right, left.v[2] + right};
291 : }
292 : KOKKOS_INLINE_FUNCTION Real3
293 : operator+(const Real3 left, const Real3 right)
294 : {
295 : return {left.v[0] + right.v[0], left.v[1] + right.v[1], left.v[2] + right.v[2]};
296 : }
297 : KOKKOS_INLINE_FUNCTION Real3
298 : operator-(const Real left, const Real3 right)
299 : {
300 : return {left - right.v[0], left - right.v[1], left - right.v[2]};
301 : }
302 : KOKKOS_INLINE_FUNCTION Real3
303 : operator-(const Real3 left, const Real right)
304 : {
305 : return {left.v[0] - right, left.v[1] - right, left.v[2] - right};
306 : }
307 : KOKKOS_INLINE_FUNCTION Real3
308 0 : operator-(const Real3 left, const Real3 right)
309 : {
310 0 : return {left.v[0] - right.v[0], left.v[1] - right.v[1], left.v[2] - right.v[2]};
311 : }
312 : #endif
313 :
314 : template <typename T1, typename T2>
315 : struct Pair
316 : {
317 : T1 first;
318 : T2 second;
319 :
320 : template <typename T3, typename T4>
321 0 : auto & operator=(const std::pair<T3, T4> pair)
322 : {
323 0 : first = pair.first;
324 0 : second = pair.second;
325 :
326 0 : return *this;
327 : }
328 : };
329 :
330 : template <typename T1, typename T2>
331 : bool
332 40526 : operator<(const Pair<T1, T2> & left, const Pair<T1, T2> & right)
333 : {
334 40526 : return std::make_pair(left.first, left.second) < std::make_pair(right.first, right.second);
335 : }
336 :
337 : } // namespace Kokkos
338 : } // namespace Moose
|
