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 "KokkosArray.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 5208290 : KOKKOS_INLINE_FUNCTION Real33() { *this = 0; }
32 : KOKKOS_INLINE_FUNCTION Real33(const Real & scalar) { *this = scalar; }
33 162077182 : KOKKOS_INLINE_FUNCTION Real33(const Real33 & tensor) { *this = tensor; }
34 12482020 : KOKKOS_INLINE_FUNCTION Real & operator()(unsigned int i, unsigned int j) { return a[i][j]; }
35 1194323364 : KOKKOS_INLINE_FUNCTION Real operator()(unsigned int i, unsigned int j) const { return a[i][j]; }
36 173854028 : KOKKOS_INLINE_FUNCTION Real33 & operator=(const Real33 & tensor)
37 : {
38 695416112 : for (unsigned int i = 0; i < 3; ++i)
39 2086248336 : for (unsigned int j = 0; j < 3; ++j)
40 1564686252 : a[i][j] = tensor.a[i][j];
41 :
42 173854028 : return *this;
43 : }
44 5208290 : KOKKOS_INLINE_FUNCTION Real33 & operator=(const Real scalar)
45 : {
46 20833160 : for (unsigned int i = 0; i < 3; ++i)
47 62499480 : for (unsigned int j = 0; j < 3; ++j)
48 46874610 : a[i][j] = scalar;
49 :
50 5208290 : return *this;
51 : }
52 943570 : KOKKOS_INLINE_FUNCTION void operator+=(const Real33 tensor)
53 : {
54 3774280 : for (unsigned int i = 0; i < 3; ++i)
55 11322840 : for (unsigned int j = 0; j < 3; ++j)
56 8492130 : a[i][j] += tensor.a[i][j];
57 943570 : }
58 367484 : KOKKOS_INLINE_FUNCTION Real determinant(const unsigned int dim = 3)
59 : {
60 367484 : Real det = 0;
61 :
62 367484 : if (dim == 0)
63 526 : det = 1;
64 366958 : else if (dim == 1)
65 101190 : det = a[0][0];
66 265768 : else if (dim == 2)
67 244264 : det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
68 21504 : else if (dim == 3)
69 21504 : det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
70 21504 : a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
71 21504 : a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
72 :
73 367484 : return det;
74 : }
75 183742 : KOKKOS_INLINE_FUNCTION Real33 inverse(const unsigned int dim = 3)
76 : {
77 183742 : Real inv_det = 1.0 / determinant(dim);
78 183742 : Real33 inv_mat;
79 :
80 183742 : if (dim == 1)
81 : {
82 6394 : inv_mat(0, 0) = inv_det;
83 : }
84 177348 : else if (dim == 2)
85 : {
86 166596 : inv_mat(0, 0) = a[1][1] * inv_det;
87 166596 : inv_mat(0, 1) = -a[0][1] * inv_det;
88 166596 : inv_mat(1, 0) = -a[1][0] * inv_det;
89 166596 : inv_mat(1, 1) = a[0][0] * inv_det;
90 : }
91 10752 : else if (dim == 3)
92 : {
93 10752 : inv_mat(0, 0) = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) * inv_det;
94 10752 : inv_mat(0, 1) = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) * inv_det;
95 10752 : inv_mat(0, 2) = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * inv_det;
96 10752 : inv_mat(1, 0) = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) * inv_det;
97 10752 : inv_mat(1, 1) = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) * inv_det;
98 10752 : inv_mat(1, 2) = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * inv_det;
99 10752 : inv_mat(2, 0) = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) * inv_det;
100 10752 : inv_mat(2, 1) = (a[0][1] * a[2][0] - a[0][0] * a[2][1]) * inv_det;
101 10752 : inv_mat(2, 2) = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * inv_det;
102 : }
103 :
104 183742 : return inv_mat;
105 : }
106 89454 : KOKKOS_INLINE_FUNCTION Real33 transpose()
107 : {
108 89454 : Real33 tr_mat;
109 :
110 357816 : for (unsigned int i = 0; i < 3; ++i)
111 1073448 : for (unsigned int j = 0; j < 3; ++j)
112 805086 : tr_mat(i, j) = a[j][i];
113 :
114 89454 : return tr_mat;
115 : }
116 : #endif
117 : };
118 :
119 : struct Real3
120 : {
121 : Real v[3];
122 :
123 : #ifdef MOOSE_KOKKOS_SCOPE
124 3993184 : KOKKOS_INLINE_FUNCTION Real3()
125 3993184 : {
126 3993184 : v[0] = 0;
127 3993184 : v[1] = 0;
128 3993184 : v[2] = 0;
129 3993184 : }
130 18228752 : KOKKOS_INLINE_FUNCTION Real3(const Real & scalar)
131 18228752 : {
132 18228752 : v[0] = scalar;
133 18228752 : v[1] = scalar;
134 18228752 : v[2] = scalar;
135 18228752 : }
136 169600728 : KOKKOS_INLINE_FUNCTION Real3(const Real3 & vector)
137 169600728 : {
138 169600728 : v[0] = vector.v[0];
139 169600728 : v[1] = vector.v[1];
140 169600728 : v[2] = vector.v[2];
141 169600728 : }
142 238374340 : KOKKOS_INLINE_FUNCTION Real3(const Real & x, const Real & y, const Real & z)
143 238374340 : {
144 238374340 : v[0] = x;
145 238374340 : v[1] = y;
146 238374340 : v[2] = z;
147 238374340 : }
148 12 : Real3(const libMesh::TypeVector<Real> & vector)
149 12 : {
150 12 : v[0] = vector(0);
151 12 : v[1] = vector(1);
152 12 : v[2] = vector(2);
153 12 : }
154 :
155 122672 : KOKKOS_INLINE_FUNCTION Real & operator()(unsigned int i) { return v[i]; }
156 3080 : KOKKOS_INLINE_FUNCTION Real operator()(unsigned int i) const { return v[i]; }
157 :
158 30005982 : KOKKOS_INLINE_FUNCTION Real3 & operator=(const Real3 & vector)
159 : {
160 30005982 : v[0] = vector.v[0];
161 30005982 : v[1] = vector.v[1];
162 30005982 : v[2] = vector.v[2];
163 :
164 30005982 : return *this;
165 : }
166 : KOKKOS_INLINE_FUNCTION Real3 & operator=(const Real scalar)
167 : {
168 : v[0] = scalar;
169 : v[1] = scalar;
170 : v[2] = scalar;
171 :
172 : return *this;
173 : }
174 167876 : Real3 & operator=(const libMesh::TypeVector<Real> & vector)
175 : {
176 167876 : v[0] = vector(0);
177 167876 : v[1] = vector(1);
178 167876 : v[2] = vector(2);
179 :
180 167876 : return *this;
181 : }
182 : KOKKOS_INLINE_FUNCTION void operator+=(const Real scalar)
183 : {
184 : v[0] += scalar;
185 : v[1] += scalar;
186 : v[2] += scalar;
187 : }
188 78541076 : KOKKOS_INLINE_FUNCTION void operator+=(const Real3 vector)
189 : {
190 78541076 : v[0] += vector.v[0];
191 78541076 : v[1] += vector.v[1];
192 78541076 : v[2] += vector.v[2];
193 78541076 : }
194 : KOKKOS_INLINE_FUNCTION void operator-=(const Real scalar)
195 : {
196 : v[0] -= scalar;
197 : v[1] -= scalar;
198 : v[2] -= scalar;
199 : }
200 : KOKKOS_INLINE_FUNCTION void operator-=(const Real3 vector)
201 : {
202 : v[0] -= vector.v[0];
203 : v[1] -= vector.v[1];
204 : v[2] -= vector.v[2];
205 : }
206 : KOKKOS_INLINE_FUNCTION void operator*=(const Real scalar)
207 : {
208 : v[0] *= scalar;
209 : v[1] *= scalar;
210 : v[2] *= scalar;
211 : }
212 0 : KOKKOS_INLINE_FUNCTION Real norm() { return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); }
213 : KOKKOS_INLINE_FUNCTION Real dot_product(const Real3 vector)
214 : {
215 : return v[0] * vector.v[0] + v[1] * vector.v[1] + v[2] * vector.v[2];
216 : }
217 0 : KOKKOS_INLINE_FUNCTION Real3 cross_product(const Real3 vector)
218 : {
219 0 : Real3 cross;
220 :
221 0 : cross.v[0] = v[1] * vector.v[2] - v[2] * vector.v[1];
222 0 : cross.v[1] = v[2] * vector.v[0] - v[0] * vector.v[2];
223 0 : cross.v[2] = v[0] * vector.v[1] - v[1] * vector.v[0];
224 :
225 0 : return cross;
226 : }
227 943570 : KOKKOS_INLINE_FUNCTION Real33 cartesian_product(const Real3 vector)
228 : {
229 943570 : Real33 tensor;
230 :
231 3774280 : for (unsigned int i = 0; i < 3; ++i)
232 11322840 : for (unsigned int j = 0; j < 3; ++j)
233 8492130 : tensor(i, j) = v[i] * vector.v[j];
234 :
235 943570 : return tensor;
236 : }
237 : #endif
238 : };
239 :
240 : #ifdef MOOSE_KOKKOS_SCOPE
241 : KOKKOS_INLINE_FUNCTION Real3
242 106208468 : operator*(const Real left, const Real3 right)
243 : {
244 106208468 : return {left * right.v[0], left * right.v[1], left * right.v[2]};
245 : }
246 : KOKKOS_INLINE_FUNCTION Real3
247 : operator*(const Real3 left, const Real right)
248 : {
249 : return {left.v[0] * right, left.v[1] * right, left.v[2] * right};
250 : }
251 : KOKKOS_INLINE_FUNCTION Real
252 38644000 : operator*(const Real3 left, const Real3 right)
253 : {
254 38644000 : return left.v[0] * right.v[0] + left.v[1] * right.v[1] + left.v[2] * right.v[2];
255 : }
256 : KOKKOS_INLINE_FUNCTION Real3
257 132165872 : operator*(const Real33 left, const Real3 right)
258 : {
259 132165872 : return {left(0, 0) * right.v[0] + left(0, 1) * right.v[1] + left(0, 2) * right.v[2],
260 132165872 : left(1, 0) * right.v[0] + left(1, 1) * right.v[1] + left(1, 2) * right.v[2],
261 132165872 : left(2, 0) * right.v[0] + left(2, 1) * right.v[1] + left(2, 2) * right.v[2]};
262 : }
263 : KOKKOS_INLINE_FUNCTION Real33
264 89454 : operator*(const Real33 left, const Real33 right)
265 : {
266 89454 : Real33 mul;
267 :
268 357816 : for (unsigned int i = 0; i < 3; ++i)
269 1073448 : for (unsigned int j = 0; j < 3; ++j)
270 3220344 : for (unsigned int k = 0; k < 3; ++k)
271 2415258 : mul(i, j) += left(i, k) * right(k, j);
272 :
273 89454 : return mul;
274 : }
275 : KOKKOS_INLINE_FUNCTION Real3
276 : operator+(const Real left, const Real3 right)
277 : {
278 : return {left + right.v[0], left + right.v[1], left + right.v[2]};
279 : }
280 : KOKKOS_INLINE_FUNCTION Real3
281 : operator+(const Real3 left, const Real right)
282 : {
283 : return {left.v[0] + right, left.v[1] + right, left.v[2] + right};
284 : }
285 : KOKKOS_INLINE_FUNCTION Real3
286 : operator+(const Real3 left, const Real3 right)
287 : {
288 : return {left.v[0] + right.v[0], left.v[1] + right.v[1], left.v[2] + right.v[2]};
289 : }
290 : KOKKOS_INLINE_FUNCTION Real3
291 : operator-(const Real left, const Real3 right)
292 : {
293 : return {left - right.v[0], left - right.v[1], left - right.v[2]};
294 : }
295 : KOKKOS_INLINE_FUNCTION Real3
296 : operator-(const Real3 left, const Real right)
297 : {
298 : return {left.v[0] - right, left.v[1] - right, left.v[2] - right};
299 : }
300 : KOKKOS_INLINE_FUNCTION Real3
301 0 : operator-(const Real3 left, const Real3 right)
302 : {
303 0 : return {left.v[0] - right.v[0], left.v[1] - right.v[1], left.v[2] - right.v[2]};
304 : }
305 : #endif
306 :
307 : template <typename T1, typename T2>
308 : struct Pair
309 : {
310 : T1 first;
311 : T2 second;
312 :
313 : template <typename T3, typename T4>
314 0 : auto & operator=(const std::pair<T3, T4> pair)
315 : {
316 0 : first = pair.first;
317 0 : second = pair.second;
318 :
319 0 : return *this;
320 : }
321 : };
322 :
323 : template <typename T1, typename T2>
324 : bool
325 8788 : operator<(const Pair<T1, T2> & left, const Pair<T1, T2> & right)
326 : {
327 8788 : return std::make_pair(left.first, left.second) < std::make_pair(right.first, right.second);
328 : }
329 :
330 : } // namespace Kokkos
331 : } // namespace Moose
|
