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