1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #ifndef MFEM_FEM_KERNELS_HPP
13 #define MFEM_FEM_KERNELS_HPP
14
15 #include "../config/config.hpp"
16 #include "../linalg/dtensor.hpp"
17
18 namespace mfem
19 {
20
21 namespace kernels
22 {
23
24 // Experimental helper functions for MFEM_FORALL FEM kernels
25 // For the 2D functions, NBZ should be tied to '1' for now
26 namespace internal
27 {
28
29 /// Load B1d matrice into shared memory
30 template<int MD1, int MQ1>
LoadB(const int D1D,const int Q1D,const ConstDeviceMatrix & b,double (& sB)[MQ1 * MD1])31 MFEM_HOST_DEVICE inline void LoadB(const int D1D, const int Q1D,
32 const ConstDeviceMatrix &b,
33 double (&sB)[MQ1*MD1])
34 {
35 const int tidz = MFEM_THREAD_ID(z);
36 DeviceMatrix B(sB, D1D, Q1D);
37
38 if (tidz == 0)
39 {
40 MFEM_FOREACH_THREAD(d,y,D1D)
41 {
42 MFEM_FOREACH_THREAD(q,x,Q1D)
43 {
44 B(d,q) = b(q,d);
45 }
46 }
47 }
48 MFEM_SYNC_THREAD;
49 }
50
51 /// Load Bt1d matrices into shared memory
52 template<int MD1, int MQ1>
LoadBt(const int D1D,const int Q1D,const ConstDeviceMatrix & b,double (& sB)[MQ1 * MD1])53 MFEM_HOST_DEVICE inline void LoadBt(const int D1D, const int Q1D,
54 const ConstDeviceMatrix &b,
55 double (&sB)[MQ1*MD1])
56 {
57 const int tidz = MFEM_THREAD_ID(z);
58 DeviceMatrix Bt(sB, Q1D, D1D);
59
60 if (tidz == 0)
61 {
62 MFEM_FOREACH_THREAD(d,y,D1D)
63 {
64 MFEM_FOREACH_THREAD(q,x,Q1D)
65 {
66 Bt(q,d) = b(q,d);
67 }
68 }
69 }
70 MFEM_SYNC_THREAD;
71 }
72
73 /// Load B1d & G1d matrices into shared memory
74 template<int MD1, int MQ1>
LoadBG(const int D1D,const int Q1D,const ConstDeviceMatrix & b,const ConstDeviceMatrix & g,double (& sBG)[2][MQ1 * MD1])75 MFEM_HOST_DEVICE inline void LoadBG(const int D1D, const int Q1D,
76 const ConstDeviceMatrix &b,
77 const ConstDeviceMatrix &g,
78 double (&sBG)[2][MQ1*MD1])
79 {
80 const int tidz = MFEM_THREAD_ID(z);
81 DeviceMatrix B(sBG[0], D1D, Q1D);
82 DeviceMatrix G(sBG[1], D1D, Q1D);
83
84 if (tidz == 0)
85 {
86 MFEM_FOREACH_THREAD(d,y,D1D)
87 {
88 MFEM_FOREACH_THREAD(q,x,Q1D)
89 {
90 B(d,q) = b(q,d);
91 G(d,q) = g(q,d);
92 }
93 }
94 }
95 MFEM_SYNC_THREAD;
96 }
97
98 /// Load Bt1d & Gt1d matrices into shared memory
99 template<int MD1, int MQ1>
LoadBGt(const int D1D,const int Q1D,const ConstDeviceMatrix & b,const ConstDeviceMatrix & g,double (& sBG)[2][MQ1 * MD1])100 MFEM_HOST_DEVICE inline void LoadBGt(const int D1D, const int Q1D,
101 const ConstDeviceMatrix &b,
102 const ConstDeviceMatrix &g,
103 double (&sBG)[2][MQ1*MD1])
104 {
105 const int tidz = MFEM_THREAD_ID(z);
106 DeviceMatrix Bt(sBG[0], Q1D, D1D);
107 DeviceMatrix Gt(sBG[1], Q1D, D1D);
108
109 if (tidz == 0)
110 {
111 MFEM_FOREACH_THREAD(d,y,D1D)
112 {
113 MFEM_FOREACH_THREAD(q,x,Q1D)
114 {
115 Bt(q,d) = b(q,d);
116 Gt(q,d) = g(q,d);
117 }
118 }
119 }
120 MFEM_SYNC_THREAD;
121 }
122
123 /// Load 2D input scalar into shared memory
124 template<int MD1, int NBZ>
LoadX(const int e,const int D1D,const DeviceTensor<3,const double> & x,double (& sX)[NBZ][MD1 * MD1])125 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
126 const DeviceTensor<3, const double> &x,
127 double (&sX)[NBZ][MD1*MD1])
128 {
129 const int tidz = MFEM_THREAD_ID(z);
130 DeviceMatrix X(sX[tidz], D1D, D1D);
131
132 MFEM_FOREACH_THREAD(dy,y,D1D)
133 {
134 MFEM_FOREACH_THREAD(dx,x,D1D)
135 {
136 X(dx,dy) = x(dx,dy,e);
137 }
138 }
139 MFEM_SYNC_THREAD;
140 }
141
142 /// Load 2D input scalar into shared memory, with comp
LoadX(const int e,const int D1D,const int c,const DeviceTensor<4,const double> & x,DeviceMatrix & DD)143 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
144 const DeviceTensor<4, const double> &x,
145 DeviceMatrix &DD)
146 {
147 MFEM_FOREACH_THREAD(dy,y,D1D)
148 {
149 MFEM_FOREACH_THREAD(dx,x,D1D)
150 {
151 DD(dx,dy) = x(dx,dy,c,e);
152 }
153 }
154 MFEM_SYNC_THREAD;
155 }
156
157 template<int MD1, int NBZ>
LoadX(const int e,const int D1D,const int c,const DeviceTensor<4,const double> & x,double (& sm)[NBZ][MD1 * MD1])158 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
159 const DeviceTensor<4, const double> &x,
160 double (&sm)[NBZ][MD1*MD1])
161 {
162 const int tidz = MFEM_THREAD_ID(z);
163 DeviceMatrix DD(sm[tidz], D1D, D1D);
164 LoadX(e,D1D,c,x,DD);
165 }
166
167 /// 2D Scalar Evaluation, 1/2
EvalX(const int D1D,const int Q1D,ConstDeviceMatrix & B,DeviceMatrix & DD,DeviceMatrix & DQ)168 MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
169 ConstDeviceMatrix &B,
170 DeviceMatrix &DD,
171 DeviceMatrix &DQ)
172 {
173 MFEM_FOREACH_THREAD(dy,y,D1D)
174 {
175 MFEM_FOREACH_THREAD(qx,x,Q1D)
176 {
177 double u = 0.0;
178 for (int dx = 0; dx < D1D; ++dx)
179 {
180 u += B(dx,qx) * DD(dx,dy);
181 }
182 DQ(dy,qx) = u;
183 }
184 }
185 MFEM_SYNC_THREAD;
186 }
187
188 template<int MD1, int MQ1, int NBZ>
EvalX(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],double (& sDD)[NBZ][MD1 * MD1],double (& sDQ)[NBZ][MD1 * MQ1])189 MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
190 const double (&sB)[MQ1*MD1],
191 double (&sDD)[NBZ][MD1*MD1],
192 double (&sDQ)[NBZ][MD1*MQ1])
193 {
194 const int tidz = MFEM_THREAD_ID(z);
195 ConstDeviceMatrix B(sB, D1D, Q1D);
196 DeviceMatrix DD(sDD[tidz], D1D, D1D);
197 DeviceMatrix DQ(sDQ[tidz], D1D, Q1D);
198 EvalX(D1D,Q1D,B,DD,DQ);
199 }
200
201 /// 2D Scalar Evaluation, 2/2
EvalY(const int D1D,const int Q1D,ConstDeviceMatrix & B,DeviceMatrix & DQ,DeviceMatrix & QQ)202 MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
203 ConstDeviceMatrix &B,
204 DeviceMatrix &DQ,
205 DeviceMatrix &QQ)
206 {
207 MFEM_FOREACH_THREAD(qy,y,Q1D)
208 {
209 MFEM_FOREACH_THREAD(qx,x,Q1D)
210 {
211 double u = 0.0;
212 for (int dy = 0; dy < D1D; ++dy)
213 {
214 u += DQ(dy,qx) * B(dy,qy);
215 }
216 QQ(qx,qy) = u;
217 }
218 }
219 MFEM_SYNC_THREAD;
220 }
221
222 template<int MD1, int MQ1, int NBZ>
EvalY(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],double (& sDQ)[NBZ][MD1 * MQ1],double (& sQQ)[NBZ][MQ1 * MQ1])223 MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
224 const double (&sB)[MQ1*MD1],
225 double (&sDQ)[NBZ][MD1*MQ1],
226 double (&sQQ)[NBZ][MQ1*MQ1])
227 {
228 const int tidz = MFEM_THREAD_ID(z);
229 ConstDeviceMatrix B(sB, D1D, Q1D);
230 DeviceMatrix DQ(sDQ[tidz], D1D, Q1D);
231 DeviceMatrix QQ(sQQ[tidz], Q1D, Q1D);
232 EvalY(D1D,Q1D,B,DQ,QQ);
233 }
234
235 /// Pull 2D Scalar Evaluation
PullEval(const int qx,const int qy,DeviceMatrix & QQ,double & P)236 MFEM_HOST_DEVICE inline void PullEval(const int qx, const int qy,
237 DeviceMatrix &QQ,
238 double &P)
239 {
240 P = QQ(qx,qy);
241 }
242
243 template<int MQ1, int NBZ>
PullEval(const int Q1D,const int qx,const int qy,double (& sQQ)[NBZ][MQ1 * MQ1],double & P)244 MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
245 const int qx, const int qy,
246 double (&sQQ)[NBZ][MQ1*MQ1],
247 double &P)
248 {
249 const int tidz = MFEM_THREAD_ID(z);
250 DeviceMatrix QQ(sQQ[tidz], Q1D, Q1D);
251 PullEval(qx,qy,QQ,P);
252 }
253
254 /// Load 2D input vector into shared memory
255 template<int MD1, int NBZ>
LoadX(const int e,const int D1D,const DeviceTensor<4,const double> & X,double (& sX)[2][NBZ][MD1 * MD1])256 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
257 const DeviceTensor<4, const double> &X,
258 double (&sX)[2][NBZ][MD1*MD1])
259 {
260 const int tidz = MFEM_THREAD_ID(z);
261 DeviceMatrix X0(sX[0][tidz], D1D, D1D);
262 DeviceMatrix X1(sX[1][tidz], D1D, D1D);
263
264 MFEM_FOREACH_THREAD(dy,y,D1D)
265 {
266 MFEM_FOREACH_THREAD(dx,x,D1D)
267 {
268 X0(dx,dy) = X(dx,dy,0,e);
269 X1(dx,dy) = X(dx,dy,1,e);
270 }
271 }
272 MFEM_SYNC_THREAD;
273 }
274
275 /// 2D Evaluation, 1/2 (only B)
276 template<int MD1, int MQ1, int NBZ>
EvalX(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sX)[2][NBZ][MD1 * MD1],double (& sDQ)[2][NBZ][MD1 * MQ1])277 MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
278 const double (&sB)[MQ1*MD1],
279 const double (&sX)[2][NBZ][MD1*MD1],
280 double (&sDQ)[2][NBZ][MD1*MQ1])
281 {
282 const int tidz = MFEM_THREAD_ID(z);
283 ConstDeviceMatrix B(sB, D1D, Q1D);
284 ConstDeviceMatrix X0(sX[0][tidz], D1D, D1D);
285 ConstDeviceMatrix X1(sX[1][tidz], D1D, D1D);
286 DeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
287 DeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
288
289 MFEM_FOREACH_THREAD(dy,y,D1D)
290 {
291 MFEM_FOREACH_THREAD(qx,x,Q1D)
292 {
293 double u[2] = {0.0, 0.0};
294 for (int dx = 0; dx < D1D; ++dx)
295 {
296 const double xx = X0(dx,dy);
297 const double xy = X1(dx,dy);
298 u[0] += B(dx,qx) * xx;
299 u[1] += B(dx,qx) * xy;
300 }
301 DQ0(qx,dy) = u[0];
302 DQ1(qx,dy) = u[1];
303 }
304 }
305 MFEM_SYNC_THREAD;
306 }
307
308 /// 2D Evaluation, 2/2 (only B)
309 template<int MD1, int MQ1, int NBZ>
EvalY(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDQ)[2][NBZ][MD1 * MQ1],double (& sQQ)[2][NBZ][MQ1 * MQ1])310 MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
311 const double (&sB)[MQ1*MD1],
312 const double (&sDQ)[2][NBZ][MD1*MQ1],
313 double (&sQQ)[2][NBZ][MQ1*MQ1])
314 {
315 const int tidz = MFEM_THREAD_ID(z);
316 ConstDeviceMatrix B(sB, D1D, Q1D);
317 ConstDeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
318 ConstDeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
319 DeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
320 DeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
321
322 MFEM_FOREACH_THREAD(qy,y,Q1D)
323 {
324 MFEM_FOREACH_THREAD(qx,x,Q1D)
325 {
326 double u[2] = {0.0, 0.0};
327 for (int dy = 0; dy < D1D; ++dy)
328 {
329 u[0] += DQ0(qx,dy) * B(dy,qy);
330 u[1] += DQ1(qx,dy) * B(dy,qy);
331 }
332 QQ0(qx,qy) = u[0];
333 QQ1(qx,qy) = u[1];
334 }
335 }
336 MFEM_SYNC_THREAD;
337 }
338
339 /// Pull 2D Evaluation
340 template<int MQ1, int NBZ>
PullEval(const int Q1D,const int qx,const int qy,const double (& sQQ)[2][NBZ][MQ1 * MQ1],double (& P)[2])341 MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
342 const int qx, const int qy,
343 const double (&sQQ)[2][NBZ][MQ1*MQ1],
344 double (&P)[2])
345 {
346 const int tidz = MFEM_THREAD_ID(z);
347 ConstDeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
348 ConstDeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
349
350 P[0] = QQ0(qx,qy);
351 P[1] = QQ1(qx,qy);
352 }
353
354 /// Push 2D Evaluation
355 template<int MQ1, int NBZ>
PushEval(const int Q1D,const int qx,const int qy,const double * P,double (& sQQ)[2][NBZ][MQ1 * MQ1])356 MFEM_HOST_DEVICE inline void PushEval(const int Q1D,
357 const int qx, const int qy,
358 const double *P,
359 double (&sQQ)[2][NBZ][MQ1*MQ1])
360 {
361 const int tidz = MFEM_THREAD_ID(z);
362 DeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
363 DeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
364
365 QQ0(qx,qy) = P[0];
366 QQ1(qx,qy) = P[1];
367 }
368
369 /// 2D Transposed evaluation, 1/2
370 template<int MD1, int MQ1, int NBZ>
EvalXt(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sQQ)[2][NBZ][MQ1 * MQ1],double (& sDQ)[2][NBZ][MD1 * MQ1])371 MFEM_HOST_DEVICE inline void EvalXt(const int D1D, const int Q1D,
372 const double (&sB)[MQ1*MD1],
373 const double (&sQQ)[2][NBZ][MQ1*MQ1],
374 double (&sDQ)[2][NBZ][MD1*MQ1])
375 {
376 const int tidz = MFEM_THREAD_ID(z);
377 ConstDeviceMatrix Bt(sB, Q1D, D1D);
378 ConstDeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
379 ConstDeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
380 DeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
381 DeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
382
383 MFEM_FOREACH_THREAD(qy,y,Q1D)
384 {
385 MFEM_FOREACH_THREAD(dx,x,D1D)
386 {
387 double u[2] = {0.0, 0.0};
388 for (int qx = 0; qx < Q1D; ++qx)
389 {
390 u[0] += QQ0(qx,qy) * Bt(qx,dx);
391 u[1] += QQ1(qx,qy) * Bt(qx,dx);
392 }
393 DQ0(qy,dx) = u[0];
394 DQ1(qy,dx) = u[1];
395 }
396 }
397 MFEM_SYNC_THREAD;
398 }
399
400 /// 2D Transposed evaluation, 2/2
401 template<int MD1, int MQ1, int NBZ>
EvalYt(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDQ)[2][NBZ][MD1 * MQ1],const DeviceTensor<4> & Y,const int e)402 MFEM_HOST_DEVICE inline void EvalYt(const int D1D, const int Q1D,
403 const double (&sB)[MQ1*MD1],
404 const double (&sDQ)[2][NBZ][MD1*MQ1],
405 const DeviceTensor<4> &Y, // output
406 const int e)
407 {
408 const int tidz = MFEM_THREAD_ID(z);
409 ConstDeviceMatrix Bt(sB, Q1D, D1D);
410 ConstDeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
411 ConstDeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
412
413 MFEM_FOREACH_THREAD(dy,y,D1D)
414 {
415 MFEM_FOREACH_THREAD(dx,x,D1D)
416 {
417 double u[2] = {0.0, 0.0};
418 for (int qy = 0; qy < Q1D; ++qy)
419 {
420 u[0] += Bt(qy,dy) * DQ0(qy,dx);
421 u[1] += Bt(qy,dy) * DQ1(qy,dx);
422 }
423 Y(dx,dy,0,e) += u[0];
424 Y(dx,dy,1,e) += u[1];
425 }
426 }
427 MFEM_SYNC_THREAD;
428 }
429
430 /// 2D Gradient, 1/2
431 template<int MD1, int MQ1, int NBZ>
GradX(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& sX)[2][NBZ][MD1 * MD1],double (& sDQ)[4][NBZ][MD1 * MQ1])432 MFEM_HOST_DEVICE inline void GradX(const int D1D, const int Q1D,
433 const double (&sBG)[2][MQ1*MD1],
434 const double (&sX)[2][NBZ][MD1*MD1],
435 double (&sDQ)[4][NBZ][MD1*MQ1])
436 {
437 const int tidz = MFEM_THREAD_ID(z);
438 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
439 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
440 ConstDeviceMatrix X0(sX[0][tidz], D1D, D1D);
441 ConstDeviceMatrix X1(sX[1][tidz], D1D, D1D);
442 DeviceMatrix X0B(sDQ[0][tidz], Q1D, D1D);
443 DeviceMatrix X0G(sDQ[1][tidz], Q1D, D1D);
444 DeviceMatrix X1B(sDQ[2][tidz], Q1D, D1D);
445 DeviceMatrix X1G(sDQ[3][tidz], Q1D, D1D);
446
447 MFEM_FOREACH_THREAD(dy,y,D1D)
448 {
449 MFEM_FOREACH_THREAD(qx,x,Q1D)
450 {
451 double u[2] = {0.0, 0.0};
452 double v[2] = {0.0, 0.0};
453 for (int dx = 0; dx < D1D; ++dx)
454 {
455 const double Bx = B(dx,qx);
456 const double Gx = G(dx,qx);
457 const double x0 = X0(dx,dy);
458 const double x1 = X1(dx,dy);
459 u[0] += Bx * x0;
460 v[0] += Gx * x0;
461 u[1] += Bx * x1;
462 v[1] += Gx * x1;
463 }
464 X0B(qx,dy) = u[0];
465 X0G(qx,dy) = v[0];
466 X1B(qx,dy) = u[1];
467 X1G(qx,dy) = v[1];
468 }
469 }
470 MFEM_SYNC_THREAD;
471 }
472
473 /// 2D Gradient, 2/2
474 template<int MD1, int MQ1, int NBZ>
GradY(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& sDQ)[4][NBZ][MD1 * MQ1],double (& sQQ)[4][NBZ][MQ1 * MQ1])475 MFEM_HOST_DEVICE inline void GradY(const int D1D, const int Q1D,
476 const double (&sBG)[2][MQ1*MD1],
477 const double (&sDQ)[4][NBZ][MD1*MQ1],
478 double (&sQQ)[4][NBZ][MQ1*MQ1])
479 {
480 const int tidz = MFEM_THREAD_ID(z);
481 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
482 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
483 ConstDeviceMatrix X0B(sDQ[0][tidz], Q1D, D1D);
484 ConstDeviceMatrix X0G(sDQ[1][tidz], Q1D, D1D);
485 ConstDeviceMatrix X1B(sDQ[2][tidz], Q1D, D1D);
486 ConstDeviceMatrix X1G(sDQ[3][tidz], Q1D, D1D);
487 DeviceMatrix X0GB(sQQ[0][tidz], Q1D, Q1D);
488 DeviceMatrix X0BG(sQQ[1][tidz], Q1D, Q1D);
489 DeviceMatrix X1GB(sQQ[2][tidz], Q1D, Q1D);
490 DeviceMatrix X1BG(sQQ[3][tidz], Q1D, Q1D);
491
492 MFEM_FOREACH_THREAD(qy,y,Q1D)
493 {
494 MFEM_FOREACH_THREAD(qx,x,Q1D)
495 {
496 double u[2] = {0.0, 0.0};
497 double v[2] = {0.0, 0.0};
498 for (int dy = 0; dy < D1D; ++dy)
499 {
500 const double By = B(dy,qy);
501 const double Gy = G(dy,qy);
502 u[0] += X0G(qx,dy) * By;
503 v[0] += X0B(qx,dy) * Gy;
504 u[1] += X1G(qx,dy) * By;
505 v[1] += X1B(qx,dy) * Gy;
506 }
507 X0GB(qx,qy) = u[0];
508 X0BG(qx,qy) = v[0];
509 X1GB(qx,qy) = u[1];
510 X1BG(qx,qy) = v[1];
511 }
512 }
513 MFEM_SYNC_THREAD;
514 }
515
516 /// Pull 2D Gradient
517 template<int MQ1, int NBZ>
PullGrad(const int Q1D,const int qx,const int qy,const double (& sQQ)[4][NBZ][MQ1 * MQ1],double * Jpr)518 MFEM_HOST_DEVICE inline void PullGrad(const int Q1D,
519 const int qx, const int qy,
520 const double (&sQQ)[4][NBZ][MQ1*MQ1],
521 double *Jpr)
522 {
523 const int tidz = MFEM_THREAD_ID(z);
524 ConstDeviceMatrix X0GB(sQQ[0][tidz], Q1D, Q1D);
525 ConstDeviceMatrix X0BG(sQQ[1][tidz], Q1D, Q1D);
526 ConstDeviceMatrix X1GB(sQQ[2][tidz], Q1D, Q1D);
527 ConstDeviceMatrix X1BG(sQQ[3][tidz], Q1D, Q1D);
528
529 Jpr[0] = X0GB(qx,qy);
530 Jpr[1] = X1GB(qx,qy);
531 Jpr[2] = X0BG(qx,qy);
532 Jpr[3] = X1BG(qx,qy);
533 }
534
535 /// Push 2D Gradient
536 template<int MQ1, int NBZ>
PushGrad(const int Q1D,const int qx,const int qy,const double * A,double (& sQQ)[4][NBZ][MQ1 * MQ1])537 MFEM_HOST_DEVICE inline void PushGrad(const int Q1D,
538 const int qx, const int qy,
539 const double *A,
540 double (&sQQ)[4][NBZ][MQ1*MQ1])
541 {
542 const int tidz = MFEM_THREAD_ID(z);
543 DeviceMatrix X0GB(sQQ[0][tidz], Q1D, Q1D);
544 DeviceMatrix X0BG(sQQ[1][tidz], Q1D, Q1D);
545 DeviceMatrix X1GB(sQQ[2][tidz], Q1D, Q1D);
546 DeviceMatrix X1BG(sQQ[3][tidz], Q1D, Q1D);
547
548 X0GB(qx,qy) = A[0];
549 X1GB(qx,qy) = A[2];
550 X0BG(qx,qy) = A[1];
551 X1BG(qx,qy) = A[3];
552 }
553
554 /// 2D Transposed gradient, 1/2
555 template<int MD1, int MQ1, int NBZ>
GradYt(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& GQ)[4][NBZ][MQ1 * MQ1],double (& GD)[4][NBZ][MD1 * MQ1])556 MFEM_HOST_DEVICE inline void GradYt(const int D1D, const int Q1D,
557 const double (&sBG)[2][MQ1*MD1],
558 const double (&GQ)[4][NBZ][MQ1*MQ1],
559 double (&GD)[4][NBZ][MD1*MQ1])
560 {
561 const int tidz = MFEM_THREAD_ID(z);
562 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
563 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
564 ConstDeviceMatrix QQx0(GQ[0][tidz], Q1D, Q1D);
565 ConstDeviceMatrix QQx1(GQ[1][tidz], Q1D, Q1D);
566 ConstDeviceMatrix QQy0(GQ[2][tidz], Q1D, Q1D);
567 ConstDeviceMatrix QQy1(GQ[3][tidz], Q1D, Q1D);
568 DeviceMatrix DQxB(GD[0][tidz], Q1D, D1D);
569 DeviceMatrix DQxG(GD[1][tidz], Q1D, D1D);
570 DeviceMatrix DQyB(GD[2][tidz], Q1D, D1D);
571 DeviceMatrix DQyG(GD[3][tidz], Q1D, D1D);
572
573 MFEM_FOREACH_THREAD(qy,y,Q1D)
574 {
575 MFEM_FOREACH_THREAD(dx,x,D1D)
576 {
577 double u[2] = {0.0, 0.0};
578 double v[2] = {0.0, 0.0};
579 for (int qx = 0; qx < Q1D; ++qx)
580 {
581 u[0] += Gt(qx,dx) * QQx0(qx,qy);
582 u[1] += Gt(qx,dx) * QQy0(qx,qy);
583 v[0] += Bt(qx,dx) * QQx1(qx,qy);
584 v[1] += Bt(qx,dx) * QQy1(qx,qy);
585 }
586 DQxB(qy,dx) = u[0];
587 DQyB(qy,dx) = u[1];
588 DQxG(qy,dx) = v[0];
589 DQyG(qy,dx) = v[1];
590 }
591 }
592 MFEM_SYNC_THREAD;
593 }
594
595 /// 2D Transposed gradient, 2/2
596 template<int MD1, int MQ1, int NBZ>
GradXt(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& GD)[4][NBZ][MD1 * MQ1],const DeviceTensor<4> & Y,const int e)597 MFEM_HOST_DEVICE inline void GradXt(const int D1D, const int Q1D,
598 const double (&sBG)[2][MQ1*MD1],
599 const double (&GD)[4][NBZ][MD1*MQ1],
600 const DeviceTensor<4> &Y, // output
601 const int e)
602 {
603 const int tidz = MFEM_THREAD_ID(z);
604 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
605 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
606 ConstDeviceMatrix DQxB(GD[0][tidz], Q1D, D1D);
607 ConstDeviceMatrix DQxG(GD[1][tidz], Q1D, D1D);
608 ConstDeviceMatrix DQyB(GD[2][tidz], Q1D, D1D);
609 ConstDeviceMatrix DQyG(GD[3][tidz], Q1D, D1D);
610
611 MFEM_FOREACH_THREAD(dy,y,D1D)
612 {
613 MFEM_FOREACH_THREAD(dx,x,D1D)
614 {
615 double u[2] = {0.0, 0.0};
616 double v[2] = {0.0, 0.0};
617 for (int qy = 0; qy < Q1D; ++qy)
618 {
619 u[0] += DQxB(qy,dx) * Bt(qy,dy);
620 u[1] += DQyB(qy,dx) * Bt(qy,dy);
621 v[0] += DQxG(qy,dx) * Gt(qy,dy);
622 v[1] += DQyG(qy,dx) * Gt(qy,dy);
623 }
624 Y(dx,dy,0,e) += u[0] + v[0];
625 Y(dx,dy,1,e) += u[1] + v[1];
626 }
627 }
628 MFEM_SYNC_THREAD;
629 }
630
631 /// Load 3D scalar input vector into shared memory
LoadX(const int e,const int D1D,const DeviceTensor<4,const double> & x,DeviceCube & X)632 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
633 const DeviceTensor<4, const double> &x,
634 DeviceCube &X)
635 {
636 MFEM_FOREACH_THREAD(dz,z,D1D)
637 {
638 MFEM_FOREACH_THREAD(dy,y,D1D)
639 {
640 MFEM_FOREACH_THREAD(dx,x,D1D)
641 {
642 X(dx,dy,dz) = x(dx,dy,dz,e);
643 }
644 }
645 }
646 MFEM_SYNC_THREAD;
647 }
648
649 template<int MD1>
LoadX(const int e,const int D1D,const DeviceTensor<4,const double> & x,double (& sm)[MD1 * MD1 * MD1])650 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
651 const DeviceTensor<4, const double> &x,
652 double (&sm)[MD1*MD1*MD1])
653 {
654 DeviceCube X(sm, D1D,D1D,D1D);
655 LoadX(e,D1D,x,X);
656 }
657
658 /// Load 3D scalar input vector into shared memory, with comp & DeviceTensor
LoadX(const int e,const int D1D,const int c,const DeviceTensor<5,const double> & x,DeviceTensor<3> & X)659 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
660 const DeviceTensor<5, const double> &x,
661 DeviceTensor<3> &X)
662 {
663 MFEM_FOREACH_THREAD(dz,z,D1D)
664 {
665 MFEM_FOREACH_THREAD(dy,y,D1D)
666 {
667 MFEM_FOREACH_THREAD(dx,x,D1D)
668 {
669 X(dx,dy,dz) = x(dx,dy,dz,c,e);
670 }
671 }
672 }
673 MFEM_SYNC_THREAD;
674 }
675
676 /// Load 3D scalar input vector into shared memory, with comp & pointer
677 template<int MD1>
LoadX(const int e,const int D1D,const int c,const DeviceTensor<5,const double> & x,double (& sm)[MD1 * MD1 * MD1])678 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
679 const DeviceTensor<5, const double> &x,
680 double (&sm)[MD1*MD1*MD1])
681 {
682 DeviceCube X(sm, D1D, D1D, D1D);
683 return LoadX<MD1>(e,D1D,c,x,X);
684 }
685
686 /// 3D Scalar Evaluation, 1/3
EvalX(const int D1D,const int Q1D,ConstDeviceMatrix & B,const DeviceCube & DDD,DeviceCube & DDQ)687 MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
688 ConstDeviceMatrix &B,
689 const DeviceCube &DDD,
690 DeviceCube &DDQ)
691 {
692 MFEM_FOREACH_THREAD(dz,z,D1D)
693 {
694 MFEM_FOREACH_THREAD(dy,y,D1D)
695 {
696 MFEM_FOREACH_THREAD(qx,x,Q1D)
697 {
698 double u = 0.0;
699 for (int dx = 0; dx < D1D; ++dx)
700 {
701 const double Bx = B(dx,qx);
702 u += Bx * DDD(dx,dy,dz);
703 }
704 DDQ(dz,dy,qx) = u;
705 }
706 }
707 }
708 MFEM_SYNC_THREAD;
709 }
710
711 template<int MD1, int MQ1>
EvalX(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDDD)[MD1 * MD1 * MD1],double (& sDDQ)[MD1 * MD1 * MQ1])712 MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
713 const double (&sB)[MQ1*MD1],
714 const double (&sDDD)[MD1*MD1*MD1],
715 double (&sDDQ)[MD1*MD1*MQ1])
716 {
717 ConstDeviceMatrix B(sB, D1D, Q1D);
718 const DeviceCube DDD(sDDD, D1D, D1D, D1D);
719 DeviceCube DDQ(sDDQ, Q1D, D1D, D1D);
720 EvalX(D1D,Q1D,B,DDD,DDQ);
721 }
722
723 /// 3D Scalar Evaluation, 2/3
EvalY(const int D1D,const int Q1D,ConstDeviceMatrix & B,const DeviceCube & DDQ,DeviceCube & DQQ)724 MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
725 ConstDeviceMatrix &B,
726 const DeviceCube &DDQ,
727 DeviceCube &DQQ)
728 {
729 MFEM_FOREACH_THREAD(dz,z,D1D)
730 {
731 MFEM_FOREACH_THREAD(qy,y,Q1D)
732 {
733 MFEM_FOREACH_THREAD(qx,x,Q1D)
734 {
735 double u = 0.0;
736 for (int dy = 0; dy < D1D; ++dy)
737 {
738 const double By = B(dy,qy);
739 u += DDQ(dz,dy,qx) * By;
740 }
741 DQQ(dz,qy,qx) = u;
742 }
743 }
744 }
745 MFEM_SYNC_THREAD;
746 }
747
748 template<int MD1, int MQ1>
EvalY(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDDQ)[MD1 * MD1 * MQ1],double (& sDQQ)[MD1 * MQ1 * MQ1])749 MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
750 const double (&sB)[MQ1*MD1],
751 const double (&sDDQ)[MD1*MD1*MQ1],
752 double (&sDQQ)[MD1*MQ1*MQ1])
753 {
754 ConstDeviceMatrix B(sB, D1D, Q1D);
755 const DeviceCube DDQ(sDDQ, Q1D, D1D, D1D);
756 DeviceCube DQQ(sDQQ, Q1D, Q1D, D1D);
757 EvalY(D1D,Q1D,B,DDQ,DQQ);
758 }
759
760 /// 3D Scalar Evaluation, 3/3
EvalZ(const int D1D,const int Q1D,ConstDeviceMatrix & B,const DeviceCube & DQQ,DeviceCube & QQQ)761 MFEM_HOST_DEVICE inline void EvalZ(const int D1D, const int Q1D,
762 ConstDeviceMatrix &B,
763 const DeviceCube &DQQ,
764 DeviceCube &QQQ)
765 {
766 MFEM_FOREACH_THREAD(qz,z,Q1D)
767 {
768 MFEM_FOREACH_THREAD(qy,y,Q1D)
769 {
770 MFEM_FOREACH_THREAD(qx,x,Q1D)
771 {
772 double u = 0.0;
773 for (int dz = 0; dz < D1D; ++dz)
774 {
775 const double Bz = B(dz,qz);
776 u += DQQ(dz,qy,qx) * Bz;
777 }
778 QQQ(qz,qy,qx) = u;
779 }
780 }
781 }
782 MFEM_SYNC_THREAD;
783 }
784
785 template<int MD1, int MQ1>
EvalZ(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDQQ)[MD1 * MQ1 * MQ1],double (& sQQQ)[MQ1 * MQ1 * MQ1])786 MFEM_HOST_DEVICE inline void EvalZ(const int D1D, const int Q1D,
787 const double (&sB)[MQ1*MD1],
788 const double (&sDQQ)[MD1*MQ1*MQ1],
789 double (&sQQQ)[MQ1*MQ1*MQ1])
790 {
791 ConstDeviceMatrix B(sB, D1D, Q1D);
792 const DeviceCube DQQ(sDQQ, Q1D, Q1D, D1D);
793 DeviceCube QQQ(sQQQ, Q1D, Q1D, Q1D);
794 EvalZ(D1D,Q1D,B,DQQ,QQQ);
795 }
796
797 /// Pull 3D Scalar Evaluation
PullEval(const int x,const int y,const int z,const DeviceCube & QQQ,double & X)798 MFEM_HOST_DEVICE inline void PullEval(const int x, const int y, const int z,
799 const DeviceCube &QQQ,
800 double &X)
801 {
802 X = QQQ(z,y,x);
803 }
804
805 template<int MQ1>
PullEval(const int Q1D,const int x,const int y,const int z,const double (& sQQQ)[MQ1 * MQ1 * MQ1],double & X)806 MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
807 const int x, const int y, const int z,
808 const double (&sQQQ)[MQ1*MQ1*MQ1],
809 double &X)
810 {
811 const DeviceCube QQQ(sQQQ, Q1D, Q1D, Q1D);
812 PullEval(x,y,z,QQQ,X);
813 }
814
815 /// Load 3D input vector into shared memory
816 template<int MD1>
LoadX(const int e,const int D1D,const DeviceTensor<5,const double> & X,double (* sm)[MD1 * MD1 * MD1])817 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
818 const DeviceTensor<5, const double> &X,
819 double (*sm)[MD1*MD1*MD1])
820 {
821 DeviceCube Xx(sm[0], D1D, D1D, D1D);
822 DeviceCube Xy(sm[1], D1D, D1D, D1D);
823 DeviceCube Xz(sm[2], D1D, D1D, D1D);
824
825 MFEM_FOREACH_THREAD(dz,z,D1D)
826 {
827 MFEM_FOREACH_THREAD(dy,y,D1D)
828 {
829 MFEM_FOREACH_THREAD(dx,x,D1D)
830 {
831 Xx(dx,dy,dz) = X(dx,dy,dz,0,e);
832 Xy(dx,dy,dz) = X(dx,dy,dz,1,e);
833 Xz(dx,dy,dz) = X(dx,dy,dz,2,e);
834 }
835 }
836 }
837 MFEM_SYNC_THREAD;
838 }
839
840 /// 3D Vector Evaluation, 1/3 (only B)
841 template<int MD1, int MQ1>
EvalX(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDDD)[3][MD1 * MD1 * MD1],double (& sDDQ)[3][MD1 * MD1 * MQ1])842 MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
843 const double (&sB)[MQ1*MD1],
844 const double (&sDDD)[3][MD1*MD1*MD1],
845 double (&sDDQ)[3][MD1*MD1*MQ1])
846 {
847 ConstDeviceMatrix B(sB, D1D, Q1D);
848 ConstDeviceCube Xx(sDDD[0], D1D, D1D, D1D);
849 ConstDeviceCube Xy(sDDD[1], D1D, D1D, D1D);
850 ConstDeviceCube Xz(sDDD[2], D1D, D1D, D1D);
851 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
852 DeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
853 DeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
854
855 MFEM_FOREACH_THREAD(dz,z,D1D)
856 {
857 MFEM_FOREACH_THREAD(dy,y,D1D)
858 {
859 MFEM_FOREACH_THREAD(qx,x,Q1D)
860 {
861 double u[3] = {0.0, 0.0, 0.0};
862 for (int dx = 0; dx < D1D; ++dx)
863 {
864 const double Bx = B(dx,qx);
865 u[0] += Bx * Xx(dx,dy,dz);
866 u[1] += Bx * Xy(dx,dy,dz);
867 u[2] += Bx * Xz(dx,dy,dz);
868 }
869 XxB(qx,dy,dz) = u[0];
870 XyB(qx,dy,dz) = u[1];
871 XzB(qx,dy,dz) = u[2];
872 }
873 }
874 }
875 MFEM_SYNC_THREAD;
876 }
877
878 /// 3D Vector Evaluation, 2/3 (only B)
879 template<int MD1, int MQ1>
EvalY(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDDQ)[3][MD1 * MD1 * MQ1],double (& sDQQ)[3][MD1 * MQ1 * MQ1])880 MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
881 const double (&sB)[MQ1*MD1],
882 const double (&sDDQ)[3][MD1*MD1*MQ1],
883 double (&sDQQ)[3][MD1*MQ1*MQ1])
884 {
885 ConstDeviceMatrix B(sB, D1D, Q1D);
886 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
887 ConstDeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
888 ConstDeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
889 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
890 DeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
891 DeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
892
893 MFEM_FOREACH_THREAD(dz,z,D1D)
894 {
895 MFEM_FOREACH_THREAD(qy,y,Q1D)
896 {
897 MFEM_FOREACH_THREAD(qx,x,Q1D)
898 {
899 double u[3] = {0.0, 0.0, 0.0};
900 for (int dy = 0; dy < D1D; ++dy)
901 {
902 const double By = B(dy,qy);
903 u[0] += XxB(qx,dy,dz) * By;
904 u[1] += XyB(qx,dy,dz) * By;
905 u[2] += XzB(qx,dy,dz) * By;
906 }
907 XxBB(qx,qy,dz) = u[0];
908 XyBB(qx,qy,dz) = u[1];
909 XzBB(qx,qy,dz) = u[2];
910 }
911 }
912 }
913 MFEM_SYNC_THREAD;
914 }
915
916 /// 3D Vector Evaluation, 3/3 (only B)
917 template<int MD1, int MQ1>
EvalZ(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDQQ)[3][MD1 * MQ1 * MQ1],double (& sQQQ)[3][MQ1 * MQ1 * MQ1])918 MFEM_HOST_DEVICE inline void EvalZ(const int D1D, const int Q1D,
919 const double (&sB)[MQ1*MD1],
920 const double (&sDQQ)[3][MD1*MQ1*MQ1],
921 double (&sQQQ)[3][MQ1*MQ1*MQ1])
922 {
923 ConstDeviceMatrix B(sB, D1D, Q1D);
924 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
925 ConstDeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
926 ConstDeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
927 DeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
928 DeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
929 DeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
930
931 MFEM_FOREACH_THREAD(qz,z,Q1D)
932 {
933 MFEM_FOREACH_THREAD(qy,y,Q1D)
934 {
935 MFEM_FOREACH_THREAD(qx,x,Q1D)
936 {
937 double u[3] = {0.0, 0.0, 0.0};
938 for (int dz = 0; dz < D1D; ++dz)
939 {
940 const double Bz = B(dz,qz);
941 u[0] += XxBB(qx,qy,dz) * Bz;
942 u[1] += XyBB(qx,qy,dz) * Bz;
943 u[2] += XzBB(qx,qy,dz) * Bz;
944 }
945 XxBBB(qx,qy,qz) = u[0];
946 XyBBB(qx,qy,qz) = u[1];
947 XzBBB(qx,qy,qz) = u[2];
948 }
949 }
950 }
951 MFEM_SYNC_THREAD;
952 }
953
954 /// Pull 3D Vector Evaluation
955 template<int MQ1>
PullEval(const int Q1D,const int x,const int y,const int z,const double (& sQQQ)[3][MQ1 * MQ1 * MQ1],double (& X)[3])956 MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
957 const int x, const int y, const int z,
958 const double (&sQQQ)[3][MQ1*MQ1*MQ1],
959 double (&X)[3])
960 {
961 ConstDeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
962 ConstDeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
963 ConstDeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
964
965 X[0] = XxBBB(x,y,z);
966 X[1] = XyBBB(x,y,z);
967 X[2] = XzBBB(x,y,z);
968 }
969
970 /// Push 3D Vector Evaluation
971 template<int MQ1>
PushEval(const int Q1D,const int x,const int y,const int z,const double (& A)[3],double (& sQQQ)[3][MQ1 * MQ1 * MQ1])972 MFEM_HOST_DEVICE inline void PushEval(const int Q1D,
973 const int x, const int y, const int z,
974 const double (&A)[3],
975 double (&sQQQ)[3][MQ1*MQ1*MQ1])
976 {
977 DeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
978 DeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
979 DeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
980
981 XxBBB(x,y,z) = A[0];
982 XyBBB(x,y,z) = A[1];
983 XzBBB(x,y,z) = A[2];
984 }
985
986 /// 3D Transposed Vector Evaluation, 1/3
987 template<int MD1, int MQ1>
EvalXt(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sQQQ)[3][MQ1 * MQ1 * MQ1],double (& sDQQ)[3][MD1 * MQ1 * MQ1])988 MFEM_HOST_DEVICE inline void EvalXt(const int D1D, const int Q1D,
989 const double (&sB)[MQ1*MD1],
990 const double (&sQQQ)[3][MQ1*MQ1*MQ1],
991 double (&sDQQ)[3][MD1*MQ1*MQ1])
992 {
993 ConstDeviceMatrix Bt(sB, Q1D, D1D);
994 ConstDeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
995 ConstDeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
996 ConstDeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
997 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
998 DeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
999 DeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
1000
1001 MFEM_FOREACH_THREAD(qz,z,Q1D)
1002 {
1003 MFEM_FOREACH_THREAD(qy,y,Q1D)
1004 {
1005 MFEM_FOREACH_THREAD(dx,x,D1D)
1006 {
1007 double u[3] = {0.0, 0.0, 0.0};
1008 for (int qx = 0; qx < Q1D; ++qx)
1009 {
1010 const double Btx = Bt(qx,dx);
1011 u[0] += XxBBB(qx,qy,qz) * Btx;
1012 u[1] += XyBBB(qx,qy,qz) * Btx;
1013 u[2] += XzBBB(qx,qy,qz) * Btx;
1014 }
1015 XxBB(qz,qy,dx) = u[0];
1016 XyBB(qz,qy,dx) = u[1];
1017 XzBB(qz,qy,dx) = u[2];
1018 }
1019 }
1020 }
1021 MFEM_SYNC_THREAD;
1022 }
1023
1024 /// 3D Transposed Vector Evaluation, 2/3
1025 template<int MD1, int MQ1>
EvalYt(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDQQ)[3][MD1 * MQ1 * MQ1],double (& sDDQ)[3][MD1 * MD1 * MQ1])1026 MFEM_HOST_DEVICE inline void EvalYt(const int D1D, const int Q1D,
1027 const double (&sB)[MQ1*MD1],
1028 const double (&sDQQ)[3][MD1*MQ1*MQ1],
1029 double (&sDDQ)[3][MD1*MD1*MQ1])
1030 {
1031 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1032 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1033 ConstDeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
1034 ConstDeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
1035 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1036 DeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
1037 DeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
1038
1039 MFEM_FOREACH_THREAD(qz,z,Q1D)
1040 {
1041 MFEM_FOREACH_THREAD(dy,y,D1D)
1042 {
1043 MFEM_FOREACH_THREAD(dx,x,D1D)
1044 {
1045 double u[3] = {0.0, 0.0, 0.0};
1046 for (int qy = 0; qy < Q1D; ++qy)
1047 {
1048 const double Bty = Bt(qy,dy);
1049 u[0] += XxBB(qz,qy,dx) * Bty;
1050 u[1] += XyBB(qz,qy,dx) * Bty;
1051 u[2] += XzBB(qz,qy,dx) * Bty;
1052
1053 }
1054 XxB(qz,dy,dx) = u[0];
1055 XyB(qz,dy,dx) = u[1];
1056 XzB(qz,dy,dx)= u[2];
1057 }
1058 }
1059 }
1060 MFEM_SYNC_THREAD;
1061 }
1062
1063 /// 3D Transposed Vector Evaluation, 3/3
1064 template<int MD1, int MQ1>
EvalZt(const int D1D,const int Q1D,const double (& sB)[MQ1 * MD1],const double (& sDDQ)[3][MD1 * MD1 * MQ1],const DeviceTensor<5> & Y,const int e)1065 MFEM_HOST_DEVICE inline void EvalZt(const int D1D, const int Q1D,
1066 const double (&sB)[MQ1*MD1],
1067 const double (&sDDQ)[3][MD1*MD1*MQ1],
1068 const DeviceTensor<5> &Y, // output
1069 const int e)
1070 {
1071 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1072 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1073 ConstDeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
1074 ConstDeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
1075
1076 MFEM_FOREACH_THREAD(dz,z,D1D)
1077 {
1078 MFEM_FOREACH_THREAD(dy,y,D1D)
1079 {
1080 MFEM_FOREACH_THREAD(dx,x,D1D)
1081 {
1082 double u[3] = {0.0, 0.0, 0.0};
1083 for (int qz = 0; qz < Q1D; ++qz)
1084 {
1085 const double Btz = Bt(qz,dz);
1086 u[0] += XxB(qz,dy,dx) * Btz;
1087 u[1] += XyB(qz,dy,dx) * Btz;
1088 u[2] += XzB(qz,dy,dx) * Btz;
1089 }
1090 Y(dx,dy,dz,0,e) += u[0];
1091 Y(dx,dy,dz,1,e) += u[1];
1092 Y(dx,dy,dz,2,e) += u[2];
1093 }
1094 }
1095 }
1096 }
1097
1098 /// 3D Gradient, 1/3
1099 template<int MD1, int MQ1>
GradX(const int D1D,const int Q1D,const double (* sBG)[MQ1 * MD1],const double (* sDDD)[MD1 * MD1 * MD1],double (* sDDQ)[MD1 * MD1 * MQ1])1100 MFEM_HOST_DEVICE inline void GradX(const int D1D, const int Q1D,
1101 const double (*sBG)[MQ1*MD1],
1102 const double (*sDDD)[MD1*MD1*MD1],
1103 double (*sDDQ)[MD1*MD1*MQ1])
1104 {
1105 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1106 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1107 ConstDeviceCube Xx(sDDD[0], D1D, D1D, D1D);
1108 ConstDeviceCube Xy(sDDD[1], D1D, D1D, D1D);
1109 ConstDeviceCube Xz(sDDD[2], D1D, D1D, D1D);
1110 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1111 DeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
1112 DeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
1113 DeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
1114 DeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
1115 DeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
1116
1117 MFEM_FOREACH_THREAD(dz,z,D1D)
1118 {
1119 MFEM_FOREACH_THREAD(dy,y,D1D)
1120 {
1121 MFEM_FOREACH_THREAD(qx,x,Q1D)
1122 {
1123 double u[3] = {0.0, 0.0, 0.0};
1124 double v[3] = {0.0, 0.0, 0.0};
1125 for (int dx = 0; dx < D1D; ++dx)
1126 {
1127 const double xx = Xx(dx,dy,dz);
1128 const double xy = Xy(dx,dy,dz);
1129 const double xz = Xz(dx,dy,dz);
1130 const double Bx = B(dx,qx);
1131 const double Gx = G(dx,qx);
1132 u[0] += Bx * xx;
1133 u[1] += Bx * xy;
1134 u[2] += Bx * xz;
1135
1136 v[0] += Gx * xx;
1137 v[1] += Gx * xy;
1138 v[2] += Gx * xz;
1139 }
1140 XxB(qx,dy,dz) = u[0];
1141 XyB(qx,dy,dz) = u[1];
1142 XzB(qx,dy,dz) = u[2];
1143
1144 XxG(qx,dy,dz) = v[0];
1145 XyG(qx,dy,dz) = v[1];
1146 XzG(qx,dy,dz) = v[2];
1147 }
1148 }
1149 }
1150 MFEM_SYNC_THREAD;
1151 }
1152
1153 /// 3D Gradient, 2/3
1154 template<int MD1, int MQ1>
GradY(const int D1D,const int Q1D,const double (* sBG)[MQ1 * MD1],const double (* sDDQ)[MD1 * MD1 * MQ1],double (* sDQQ)[MD1 * MQ1 * MQ1])1155 MFEM_HOST_DEVICE inline void GradY(const int D1D, const int Q1D,
1156 const double (*sBG)[MQ1*MD1],
1157 const double (*sDDQ)[MD1*MD1*MQ1],
1158 double (*sDQQ)[MD1*MQ1*MQ1])
1159 {
1160 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1161 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1162 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1163 ConstDeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
1164 ConstDeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
1165 ConstDeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
1166 ConstDeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
1167 ConstDeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
1168 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1169 DeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
1170 DeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
1171 DeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
1172 DeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
1173 DeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
1174 DeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
1175 DeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
1176 DeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
1177
1178 MFEM_FOREACH_THREAD(dz,z,D1D)
1179 {
1180 MFEM_FOREACH_THREAD(qy,y,Q1D)
1181 {
1182 MFEM_FOREACH_THREAD(qx,x,Q1D)
1183 {
1184 double u[3] = {0.0, 0.0, 0.0};
1185 double v[3] = {0.0, 0.0, 0.0};
1186 double w[3] = {0.0, 0.0, 0.0};
1187 for (int dy = 0; dy < D1D; ++dy)
1188 {
1189 const double By = B(dy,qy);
1190 const double Gy = G(dy,qy);
1191
1192 u[0] += XxB(qx,dy,dz) * By;
1193 u[1] += XyB(qx,dy,dz) * By;
1194 u[2] += XzB(qx,dy,dz) * By;
1195
1196 v[0] += XxG(qx,dy,dz) * By;
1197 v[1] += XyG(qx,dy,dz) * By;
1198 v[2] += XzG(qx,dy,dz) * By;
1199
1200 w[0] += XxB(qx,dy,dz) * Gy;
1201 w[1] += XyB(qx,dy,dz) * Gy;
1202 w[2] += XzB(qx,dy,dz) * Gy;
1203 }
1204 XxBB(qx,qy,dz) = u[0];
1205 XyBB(qx,qy,dz) = u[1];
1206 XzBB(qx,qy,dz) = u[2];
1207
1208 XxBG(qx,qy,dz) = v[0];
1209 XyBG(qx,qy,dz) = v[1];
1210 XzBG(qx,qy,dz) = v[2];
1211
1212 XxGB(qx,qy,dz) = w[0];
1213 XyGB(qx,qy,dz) = w[1];
1214 XzGB(qx,qy,dz) = w[2];
1215 }
1216 }
1217 }
1218 MFEM_SYNC_THREAD;
1219 }
1220
1221 /// 3D Gradient, 3/3
1222 template<int MD1, int MQ1>
GradZ(const int D1D,const int Q1D,const double (* sBG)[MQ1 * MD1],const double (* sDQQ)[MD1 * MQ1 * MQ1],double (* sQQQ)[MQ1 * MQ1 * MQ1])1223 MFEM_HOST_DEVICE inline void GradZ(const int D1D, const int Q1D,
1224 const double (*sBG)[MQ1*MD1],
1225 const double (*sDQQ)[MD1*MQ1*MQ1],
1226 double (*sQQQ)[MQ1*MQ1*MQ1])
1227 {
1228 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1229 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1230 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1231 ConstDeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
1232 ConstDeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
1233 ConstDeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
1234 ConstDeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
1235 ConstDeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
1236 ConstDeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
1237 ConstDeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
1238 ConstDeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
1239 DeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
1240 DeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
1241 DeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
1242 DeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
1243 DeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
1244 DeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
1245 DeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
1246 DeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
1247 DeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
1248
1249 MFEM_FOREACH_THREAD(qz,z,Q1D)
1250 {
1251 MFEM_FOREACH_THREAD(qy,y,Q1D)
1252 {
1253 MFEM_FOREACH_THREAD(qx,x,Q1D)
1254 {
1255 double u[3] = {0.0, 0.0, 0.0};
1256 double v[3] = {0.0, 0.0, 0.0};
1257 double w[3] = {0.0, 0.0, 0.0};
1258 for (int dz = 0; dz < D1D; ++dz)
1259 {
1260 const double Bz = B(dz,qz);
1261 const double Gz = G(dz,qz);
1262
1263 u[0] += XxBG(qx,qy,dz) * Bz;
1264 u[1] += XyBG(qx,qy,dz) * Bz;
1265 u[2] += XzBG(qx,qy,dz) * Bz;
1266
1267 v[0] += XxGB(qx,qy,dz) * Bz;
1268 v[1] += XyGB(qx,qy,dz) * Bz;
1269 v[2] += XzGB(qx,qy,dz) * Bz;
1270
1271 w[0] += XxBB(qx,qy,dz) * Gz;
1272 w[1] += XyBB(qx,qy,dz) * Gz;
1273 w[2] += XzBB(qx,qy,dz) * Gz;
1274 }
1275 XxBBG(qx,qy,qz) = u[0];
1276 XyBBG(qx,qy,qz) = u[1];
1277 XzBBG(qx,qy,qz) = u[2];
1278
1279 XxBGB(qx,qy,qz) = v[0];
1280 XyBGB(qx,qy,qz) = v[1];
1281 XzBGB(qx,qy,qz) = v[2];
1282
1283 XxGBB(qx,qy,qz)= w[0];
1284 XyGBB(qx,qy,qz) = w[1];
1285 XzGBB(qx,qy,qz) = w[2];
1286 }
1287 }
1288 }
1289 MFEM_SYNC_THREAD;
1290 }
1291
1292 /// Pull 3D Gradient
1293 template<int MQ1>
PullGrad(const int Q1D,const int x,const int y,const int z,const double (* sQQQ)[MQ1 * MQ1 * MQ1],double * Jpr)1294 MFEM_HOST_DEVICE inline void PullGrad(const int Q1D,
1295 const int x, const int y, const int z,
1296 const double (*sQQQ)[MQ1*MQ1*MQ1],
1297 double *Jpr)
1298 {
1299 ConstDeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
1300 ConstDeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
1301 ConstDeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
1302 ConstDeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
1303 ConstDeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
1304 ConstDeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
1305 ConstDeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
1306 ConstDeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
1307 ConstDeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
1308
1309 Jpr[0] = XxBBG(x,y,z);
1310 Jpr[3] = XxBGB(x,y,z);
1311 Jpr[6] = XxGBB(x,y,z);
1312 Jpr[1] = XyBBG(x,y,z);
1313 Jpr[4] = XyBGB(x,y,z);
1314 Jpr[7] = XyGBB(x,y,z);
1315 Jpr[2] = XzBBG(x,y,z);
1316 Jpr[5] = XzBGB(x,y,z);
1317 Jpr[8] = XzGBB(x,y,z);
1318 }
1319
1320 /// Push 3D Gradient
1321 template<int MQ1>
PushGrad(const int Q1D,const int x,const int y,const int z,const double * A,double (& sQQQ)[9][MQ1 * MQ1 * MQ1])1322 MFEM_HOST_DEVICE inline void PushGrad(const int Q1D,
1323 const int x, const int y, const int z,
1324 const double *A,
1325 double (&sQQQ)[9][MQ1*MQ1*MQ1])
1326 {
1327 DeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
1328 DeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
1329 DeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
1330 DeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
1331 DeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
1332 DeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
1333 DeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
1334 DeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
1335 DeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
1336
1337 XxBBG(x,y,z) = A[0];
1338 XxBGB(x,y,z) = A[1];
1339 XxGBB(x,y,z) = A[2];
1340 XyBBG(x,y,z) = A[3];
1341 XyBGB(x,y,z) = A[4];
1342 XyGBB(x,y,z) = A[5];
1343 XzBBG(x,y,z) = A[6];
1344 XzBGB(x,y,z) = A[7];
1345 XzGBB(x,y,z) = A[8];
1346 }
1347
1348 /// 3D Transposed Gradient, 1/3
1349 template<int MD1, int MQ1>
GradZt(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& sQQQ)[9][MQ1 * MQ1 * MQ1],double (& sDQQ)[9][MD1 * MQ1 * MQ1])1350 MFEM_HOST_DEVICE inline void GradZt(const int D1D, const int Q1D,
1351 const double (&sBG)[2][MQ1*MD1],
1352 const double (&sQQQ)[9][MQ1*MQ1*MQ1],
1353 double (&sDQQ)[9][MD1*MQ1*MQ1])
1354 {
1355
1356 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
1357 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
1358 ConstDeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
1359 ConstDeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
1360 ConstDeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
1361 ConstDeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
1362 ConstDeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
1363 ConstDeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
1364 ConstDeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
1365 ConstDeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
1366 ConstDeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
1367 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1368 DeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
1369 DeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
1370 DeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
1371 DeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
1372 DeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
1373 DeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
1374 DeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
1375 DeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
1376
1377 MFEM_FOREACH_THREAD(qz,z,Q1D)
1378 {
1379 MFEM_FOREACH_THREAD(qy,y,Q1D)
1380 {
1381 MFEM_FOREACH_THREAD(dx,x,D1D)
1382 {
1383 double u[3] = {0.0, 0.0, 0.0};
1384 double v[3] = {0.0, 0.0, 0.0};
1385 double w[3] = {0.0, 0.0, 0.0};
1386 for (int qx = 0; qx < Q1D; ++qx)
1387 {
1388 const double Btx = Bt(qx,dx);
1389 const double Gtx = Gt(qx,dx);
1390
1391 u[0] += XxBBG(qx,qy,qz) * Gtx;
1392 v[0] += XxBGB(qx,qy,qz) * Btx;
1393 w[0] += XxGBB(qx,qy,qz) * Btx;
1394
1395 u[1] += XyBBG(qx,qy,qz) * Gtx;
1396 v[1] += XyBGB(qx,qy,qz) * Btx;
1397 w[1] += XyGBB(qx,qy,qz) * Btx;
1398
1399 u[2] += XzBBG(qx,qy,qz) * Gtx;
1400 v[2] += XzBGB(qx,qy,qz) * Btx;
1401 w[2] += XzGBB(qx,qy,qz) * Btx;
1402 }
1403 XxBB(qz,qy,dx) = u[0];
1404 XxBG(qz,qy,dx) = v[0];
1405 XxGB(qz,qy,dx) = w[0];
1406
1407 XyBB(qz,qy,dx) = u[1];
1408 XyBG(qz,qy,dx) = v[1];
1409 XyGB(qz,qy,dx) = w[1];
1410
1411 XzBB(qz,qy,dx) = u[2];
1412 XzBG(qz,qy,dx) = v[2];
1413 XzGB(qz,qy,dx) = w[2];
1414 }
1415 }
1416 }
1417 MFEM_SYNC_THREAD;
1418 }
1419
1420 /// 3D Transposed Gradient, 2/3
1421 template<int MD1, int MQ1>
GradYt(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& sDQQ)[9][MD1 * MQ1 * MQ1],double (& sDDQ)[9][MD1 * MD1 * MQ1])1422 MFEM_HOST_DEVICE inline void GradYt(const int D1D, const int Q1D,
1423 const double (&sBG)[2][MQ1*MD1],
1424 const double (&sDQQ)[9][MD1*MQ1*MQ1],
1425 double (&sDDQ)[9][MD1*MD1*MQ1])
1426 {
1427 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
1428 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
1429 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1430 ConstDeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
1431 ConstDeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
1432 ConstDeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
1433 ConstDeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
1434 ConstDeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
1435 ConstDeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
1436 ConstDeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
1437 ConstDeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
1438 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1439 DeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
1440 DeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
1441 DeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
1442 DeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
1443 DeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
1444 DeviceCube XxC(sDDQ[6], Q1D, D1D, D1D);
1445 DeviceCube XyC(sDDQ[7], Q1D, D1D, D1D);
1446 DeviceCube XzC(sDDQ[8], Q1D, D1D, D1D);
1447
1448 MFEM_FOREACH_THREAD(qz,z,Q1D)
1449 {
1450 MFEM_FOREACH_THREAD(dy,y,D1D)
1451 {
1452 MFEM_FOREACH_THREAD(dx,x,D1D)
1453 {
1454 double u[3] = {0.0, 0.0, 0.0};
1455 double v[3] = {0.0, 0.0, 0.0};
1456 double w[3] = {0.0, 0.0, 0.0};
1457 for (int qy = 0; qy < Q1D; ++qy)
1458 {
1459 const double Bty = Bt(qy,dy);
1460 const double Gty = Gt(qy,dy);
1461
1462 u[0] += XxBB(qz,qy,dx) * Bty;
1463 v[0] += XxBG(qz,qy,dx) * Gty;
1464 w[0] += XxGB(qz,qy,dx) * Bty;
1465
1466 u[1] += XyBB(qz,qy,dx) * Bty;
1467 v[1] += XyBG(qz,qy,dx) * Gty;
1468 w[1] += XyGB(qz,qy,dx) * Bty;
1469
1470 u[2] += XzBB(qz,qy,dx) * Bty;
1471 v[2] += XzBG(qz,qy,dx) * Gty;
1472 w[2] += XzGB(qz,qy,dx) * Bty;
1473
1474 }
1475 XxB(qz,dy,dx) = u[0];
1476 XxC(qz,dy,dx) = v[0];
1477 XxG(qz,dy,dx) = w[0];
1478
1479 XyB(qz,dy,dx) = u[1];
1480 XyC(qz,dy,dx) = v[1];
1481 XyG(qz,dy,dx) = w[1];
1482
1483 XzB(qz,dy,dx) = u[2];
1484 XzC(qz,dy,dx) = v[2];
1485 XzG(qz,dy,dx) = w[2];
1486 }
1487 }
1488 }
1489 MFEM_SYNC_THREAD;
1490 }
1491
1492 /// 3D Transposed Gradient, 3/3
1493 template<int MD1, int MQ1>
GradXt(const int D1D,const int Q1D,const double (& sBG)[2][MQ1 * MD1],const double (& sDDQ)[9][MD1 * MD1 * MQ1],const DeviceTensor<5> & Y,const int e)1494 MFEM_HOST_DEVICE inline void GradXt(const int D1D, const int Q1D,
1495 const double (&sBG)[2][MQ1*MD1],
1496 const double (&sDDQ)[9][MD1*MD1*MQ1],
1497 const DeviceTensor<5> &Y, // output
1498 const int e)
1499 {
1500 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
1501 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
1502 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1503 ConstDeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
1504 ConstDeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
1505 ConstDeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
1506 ConstDeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
1507 ConstDeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
1508 ConstDeviceCube XxC(sDDQ[6], Q1D, D1D, D1D);
1509 ConstDeviceCube XyC(sDDQ[7], Q1D, D1D, D1D);
1510 ConstDeviceCube XzC(sDDQ[8], Q1D, D1D, D1D);
1511
1512 MFEM_FOREACH_THREAD(dz,z,D1D)
1513 {
1514 MFEM_FOREACH_THREAD(dy,y,D1D)
1515 {
1516 MFEM_FOREACH_THREAD(dx,x,D1D)
1517 {
1518 double u[3] = {0.0, 0.0, 0.0};
1519 double v[3] = {0.0, 0.0, 0.0};
1520 double w[3] = {0.0, 0.0, 0.0};
1521 for (int qz = 0; qz < Q1D; ++qz)
1522 {
1523 const double Btz = Bt(qz,dz);
1524 const double Gtz = Gt(qz,dz);
1525
1526 u[0] += XxB(qz,dy,dx) * Btz;
1527 v[0] += XxC(qz,dy,dx) * Btz;
1528 w[0] += XxG(qz,dy,dx) * Gtz;
1529
1530 u[1] += XyB(qz,dy,dx) * Btz;
1531 v[1] += XyC(qz,dy,dx)* Btz;
1532 w[1] += XyG(qz,dy,dx) * Gtz;
1533
1534 u[2] += XzB(qz,dy,dx) * Btz;
1535 v[2] += XzC(qz,dy,dx) * Btz;
1536 w[2] += XzG(qz,dy,dx) * Gtz;
1537 }
1538 Y(dx,dy,dz,0,e) += u[0] + v[0] + w[0];
1539 Y(dx,dy,dz,1,e) += u[1] + v[1] + w[1];
1540 Y(dx,dy,dz,2,e) += u[2] + v[2] + w[2];
1541 }
1542 }
1543 }
1544 }
1545
1546 } // namespace kernels::internal
1547
1548 } // namespace kernels
1549
1550 } // namespace mfem
1551
1552 #endif // MFEM_FEM_KERNELS_HPP
1553