MFEM v4.9.0
Finite element discretization library
Loading...
Searching...
No Matches
kernels.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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#include "../linalg/tensor.hpp"
18
19namespace mfem
20{
21
22namespace kernels
23{
24
25// Experimental helper functions for mfem::forall FEM kernels
26// For the 2D functions, NBZ should be tied to '1' for now
27namespace internal
28{
29
30// Types for tensors mapped to registers
31// - N is the number of threads in each of the x and y dimensions
32// - N should not be greater than 32, to have a maximum of 1024 threads
33// On GPU, the last two dimensions are set to 0 to match a 2D tile of threads
34#if ((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
35 (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
36template <int N = 0>
37using s_regs2d_t = mfem::future::tensor<real_t, 0, 0>;
38
39template <int VDIM, int N>
40using v_regs2d_t = mfem::future::tensor<real_t, VDIM, 0, 0>;
41
42template <int VDIM, int DIM, int N = 0>
43using vd_regs2d_t = mfem::future::tensor<real_t, VDIM, DIM, 0, 0>;
44
45template <int N>
46using s_regs3d_t = mfem::future::tensor<real_t, N, 0, 0>;
47
48template <int VDIM, int N>
49using v_regs3d_t = mfem::future::tensor<real_t, VDIM, N, 0, 0>;
50
51template <int VDIM, int DIM, int N>
52using vd_regs3d_t = mfem::future::tensor<real_t, VDIM, DIM, N, 0, 0>;
53
54// on GPU, SetMaxOf is a no-op, for minimal register usage
55constexpr int SetMaxOf(int n) { return n; }
56#else
57template <int N>
58using s_regs2d_t = mfem::future::tensor<real_t, N, N>;
59
60template <int VDIM, int N>
61using v_regs2d_t = mfem::future::tensor<real_t, VDIM, N, N>;
62
63template <int VDIM, int DIM, int N>
64using vd_regs2d_t = mfem::future::tensor<real_t, VDIM, DIM, N, N>;
65
66template <int N>
67using s_regs3d_t = mfem::future::tensor<real_t, N, N, N>;
68
69template <int VDIM, int N>
70using v_regs3d_t = mfem::future::tensor<real_t, VDIM, N, N, N>;
71
72template <int VDIM, int DIM, int N>
73using vd_regs3d_t = mfem::future::tensor<real_t, VDIM, DIM, N, N, N>;
74
75// on CPU, get next multiple of 4, allowing better alignments
76template <int N>
77constexpr int NextMultipleOf(int n)
78{
79 static_assert(N > 0 && (N & (N - 1)) == 0, "N must be a power of 2");
80 return (n + (N - 1)) & ~(N - 1);
81}
82constexpr int SetMaxOf(int n) { return NextMultipleOf<4>(n); }
83#endif // CUDA/HIP && DEVICE_COMPILE
84
85/// Load 2D matrix into shared memory
86template <int MQ1>
87inline MFEM_HOST_DEVICE void LoadMatrix(const int d1d, const int q1d,
88 const real_t *M, real_t (*N)[MQ1])
89{
90 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
91 {
92 MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
93 {
94 N[dy][qx] = M[dy * q1d + qx];
95 }
96 }
97 MFEM_SYNC_THREAD;
98}
99
100/// Load 2D input VDIM*DIM vector into given register tensor, specific component
101template <int VDIM, int DIM, int MQ1 = 0>
102inline MFEM_HOST_DEVICE void LoadDofs2d(const int e, const int d1d, const int c,
103 const DeviceTensor<4, const real_t> &X,
104 vd_regs2d_t<VDIM, DIM, MQ1> &Y)
105{
106 for (int d = 0; d < DIM; d++)
107 {
108 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
109 {
110 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
111 {
112 Y[c][d][dy][dx] = X(dx, dy, c, e);
113 }
114 }
115 }
116 MFEM_SYNC_THREAD;
117}
118
119/// Load 2D input VDIM*DIM vector into given register tensor
120template <int VDIM, int DIM, int MQ1 = 0>
121inline MFEM_HOST_DEVICE void LoadDofs2d(const int e, const int d1d,
122 const DeviceTensor<4, const real_t> &X,
123 vd_regs2d_t<VDIM, DIM, MQ1> &Y)
124{
125 for (int c = 0; c < VDIM; ++c) { LoadDofs2d(e, d1d, c, X, Y); }
126}
127
128/// Load 2D input VDIM vector into given register tensor
129template <int VDIM, int MQ1 = 0>
130inline MFEM_HOST_DEVICE void LoadDofs2d(const int e, const int d1d,
131 const DeviceTensor<4, const real_t> &X,
132 v_regs2d_t<VDIM, MQ1> &Y)
133{
134 for (int c = 0; c < VDIM; ++c)
135 {
136 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
137 {
138 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
139 {
140 Y[c][dy][dx] = X(dx, dy, c, e);
141 }
142 }
143 }
144 MFEM_SYNC_THREAD;
145}
146
147/// Load 2D input scalar into given register tensor
148template <int MQ1 = 0>
149inline MFEM_HOST_DEVICE void LoadDofs2d(const int e, const int d1d,
150 const DeviceTensor<3, const real_t> &X,
151 s_regs2d_t<MQ1> &Y)
152{
153 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
154 {
155 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
156 {
157 Y[dy][dx] = X(dx, dy, e);
158 }
159 }
160 MFEM_SYNC_THREAD;
161}
162
163/// Write 2D vector into given device tensor, with read (i) write (j) indices
164template <int VDIM, int DIM, int MQ1 = 0>
165inline MFEM_HOST_DEVICE void WriteDofs2d(const int e, const int d1d,
166 const int i, const int j,
167 vd_regs2d_t<VDIM, DIM, MQ1> &X,
168 const DeviceTensor<4, real_t> &Y)
169{
170 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
171 {
172 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
173 {
174 real_t y = 0.0;
175 for (int d = 0; d < DIM; d++) { y += X(i, d, dy, dx); }
176 Y(dx, dy, j, e) += y;
177 }
178 }
179 MFEM_SYNC_THREAD;
180}
181
182/// Write 2D VDIM*DIM vector into given device tensor
183template <int VDIM, int DIM, int MQ1 = 0>
184inline MFEM_HOST_DEVICE void WriteDofs2d(const int e, const int d1d,
185 vd_regs2d_t<VDIM, DIM, MQ1> &X,
186 const DeviceTensor<4, real_t> &Y)
187{
188 for (int c = 0; c < VDIM; ++c) { WriteDofs2d(e, d1d, c, c, X, Y); }
189}
190
191/// Write 2D VDIM vector into given device tensor
192template <int VDIM, int MQ1 = 0>
193inline MFEM_HOST_DEVICE void WriteDofs2d(const int e, const int d1d,
194 v_regs2d_t<VDIM, MQ1> &X,
195 const DeviceTensor<4, real_t> &Y)
196{
197 for (int c = 0; c < VDIM; ++c)
198 {
199 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
200 {
201 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
202 {
203 Y(dx, dy, c, e) += X(c, dy, dx);
204 }
205 }
206 }
207 MFEM_SYNC_THREAD;
208}
209
210/// Load 3D input VDIM*DIM vector into given register tensor, specific component
211template <int VDIM, int DIM, int MQ1>
212inline MFEM_HOST_DEVICE void LoadDofs3d(const int e, const int d1d, const int c,
213 const DeviceTensor<5, const real_t> &X,
214 vd_regs3d_t<VDIM, DIM, MQ1> &Y)
215{
216 for (int d = 0; d < DIM; d++)
217 {
218 for (int dz = 0; dz < d1d; ++dz)
219 {
220 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
221 {
222 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
223 {
224 Y[c][d][dz][dy][dx] = X(dx, dy, dz, c, e);
225 }
226 }
227 }
228 }
229 MFEM_SYNC_THREAD;
230}
231
232/// Load 3D input VDIM*DIM vector into given register tensor
233template <int VDIM, int DIM, int MQ1>
234inline MFEM_HOST_DEVICE void LoadDofs3d(const int e, const int d1d,
235 const DeviceTensor<5, const real_t> &X,
236 vd_regs3d_t<VDIM, DIM, MQ1> &Y)
237{
238 for (int c = 0; c < VDIM; ++c) { LoadDofs3d(e, d1d, c, X, Y); }
239}
240
241/// Load 3D input VDIM vector into given register tensor
242template <int VDIM, int MQ1>
243inline MFEM_HOST_DEVICE void LoadDofs3d(const int e, const int d1d,
244 const DeviceTensor<5, const real_t> &X,
245 v_regs3d_t<VDIM, MQ1> &Y)
246{
247 for (int c = 0; c < VDIM; ++c)
248 {
249 for (int dz = 0; dz < d1d; ++dz)
250 {
251 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
252 {
253 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
254 {
255 Y[c][dz][dy][dx] = X(dx,dy,dz,c,e);
256 }
257 }
258 }
259 }
260 MFEM_SYNC_THREAD;
261}
262
263/// Load 3D input scalar into given register tensor
264template <int MQ1>
265inline MFEM_HOST_DEVICE void LoadDofs3d(const int e, const int d1d,
266 const DeviceTensor<4, const real_t> &X,
267 s_regs3d_t<MQ1> &Y)
268{
269 for (int dz = 0; dz < d1d; ++dz)
270 {
271 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
272 {
273 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
274 {
275 Y[dz][dy][dx] = X(dx,dy,dz,e);
276 }
277 }
278 }
279 MFEM_SYNC_THREAD;
280}
281
282/// Write 3D scalar into given device tensor, with read (i) write (j) indices
283template <int VDIM, int DIM, int MQ1>
284inline MFEM_HOST_DEVICE void WriteDofs3d(const int e, const int d1d,
285 const int i, const int j,
286 vd_regs3d_t<VDIM, DIM, MQ1> &X,
287 const DeviceTensor<5, real_t> &Y)
288{
289 for (int dz = 0; dz < d1d; ++dz)
290 {
291 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
292 {
293 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
294 {
295 real_t value = 0.0;
296 for (int d = 0; d < DIM; d++) { value += X(i, d, dz, dy, dx); }
297 Y(dx, dy, dz, j, e) += value;
298 }
299 }
300 }
301 MFEM_SYNC_THREAD;
302}
303
304/// Write 3D VDIM*DIM vector into given device tensor
305template <int VDIM, int DIM, int MQ1>
306inline MFEM_HOST_DEVICE void WriteDofs3d(const int e, const int d1d,
307 vd_regs3d_t<VDIM, DIM, MQ1> &X,
308 const DeviceTensor<5, real_t> &Y)
309{
310 for (int c = 0; c < VDIM; ++c) { WriteDofs3d(e, d1d, c, c, X, Y); }
311}
312
313/// Write 3D VDIM vector into given device tensor
314template <int VDIM, int MQ1>
315inline MFEM_HOST_DEVICE void WriteDofs3d(const int e, const int d1d,
316 v_regs3d_t<VDIM, MQ1> &X,
317 const DeviceTensor<5, real_t> &Y)
318{
319 for (int c = 0; c < VDIM; ++c)
320 {
321 for (int dz = 0; dz < d1d; ++dz)
322 {
323 MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
324 {
325 MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
326 {
327 Y(dx, dy, dz, c, e) += X(c, dz, dy, dx);
328 }
329 }
330 }
331 }
332 MFEM_SYNC_THREAD;
333}
334
335/// 2D scalar contraction, X direction
336template <bool Transpose, int MQ1>
337inline MFEM_HOST_DEVICE void ContractX2d(const int d1d, const int q1d,
338 real_t (&smem)[MQ1][MQ1],
339 const real_t (*B)[MQ1],
340 const s_regs2d_t<MQ1> &X,
341 s_regs2d_t<MQ1> &Y)
342{
343 MFEM_FOREACH_THREAD_DIRECT(y, y, d1d)
344 {
345 MFEM_FOREACH_THREAD_DIRECT(x, x, (Transpose ? q1d : d1d))
346 {
347 smem[y][x] = X[y][x];
348 }
349 }
350 MFEM_SYNC_THREAD;
351 MFEM_FOREACH_THREAD_DIRECT(y, y, d1d)
352 {
353 MFEM_FOREACH_THREAD_DIRECT(x, x, (Transpose ? d1d : q1d))
354 {
355 real_t u = 0.0;
356 for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
357 {
358 u += (Transpose ? B[x][k] : B[k][x]) * smem[y][k];
359 }
360 Y[y][x] = u;
361 }
362 }
363 MFEM_SYNC_THREAD;
364}
365
366/// 2D scalar contraction, Y direction
367template <bool Transpose, int MQ1>
368inline MFEM_HOST_DEVICE void ContractY2d(const int d1d, const int q1d,
369 real_t (&smem)[MQ1][MQ1],
370 const real_t (*B)[MQ1],
371 const s_regs2d_t<MQ1> &X,
372 s_regs2d_t<MQ1> &Y)
373{
374 MFEM_FOREACH_THREAD_DIRECT(y, y, (Transpose ? q1d : d1d))
375 {
376 MFEM_FOREACH_THREAD_DIRECT(x, x, q1d) { smem[y][x] = X[y][x]; }
377 }
378 MFEM_SYNC_THREAD;
379 MFEM_FOREACH_THREAD_DIRECT(y, y, (Transpose ? d1d : q1d))
380 {
381 MFEM_FOREACH_THREAD_DIRECT(x, x, q1d)
382 {
383 real_t u = 0.0;
384 for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
385 {
386 u += (Transpose ? B[y][k] : B[k][y]) * smem[k][x];
387 }
388 Y[y][x] = u;
389 }
390 }
391 MFEM_SYNC_THREAD;
392}
393
394/// 2D scalar copy
395template <int MQ1 = 0>
396inline MFEM_HOST_DEVICE void Copy2d(const int q1d,
397 s_regs2d_t<MQ1> &X,
398 s_regs2d_t<MQ1> &Y)
399{
400 MFEM_FOREACH_THREAD_DIRECT(y, y, q1d)
401 {
402 MFEM_FOREACH_THREAD_DIRECT(x, x, q1d) { Y[y][x] = X[y][x]; }
403 }
404 MFEM_SYNC_THREAD;
405}
406
407/// 2D scalar contraction: X & Y directions, with additional copy
408template <bool Transpose, int MQ1>
409inline MFEM_HOST_DEVICE void Contract2d(const int d1d, const int q1d,
410 real_t (&smem)[MQ1][MQ1],
411 const real_t (*Bx)[MQ1],
412 const real_t (*By)[MQ1],
413 s_regs2d_t<MQ1> &X,
414 s_regs2d_t<MQ1> &Y)
415{
416 if (!Transpose)
417 {
418 ContractX2d<false>(d1d, q1d, smem, Bx, X, Y);
419 ContractY2d<false>(d1d, q1d, smem, By, Y, X);
420 Copy2d(q1d, X, Y);
421 }
422 else
423 {
424 Copy2d(q1d, X, Y);
425 ContractY2d<true>(d1d, q1d, smem, By, Y, X);
426 ContractX2d<true>(d1d, q1d, smem, Bx, X, Y);
427 }
428}
429
430/// 2D scalar evaluation
431template <int MQ1, bool Transpose = false>
432inline MFEM_HOST_DEVICE void Eval2d(const int d1d, const int q1d,
433 real_t (&smem)[MQ1][MQ1],
434 const real_t (*B)[MQ1],
435 s_regs2d_t<MQ1> &X,
436 s_regs2d_t<MQ1> &Y)
437{
438 Contract2d<Transpose, MQ1>(d1d, q1d, smem, B, B, X, Y);
439}
440
441/// 2D vector evaluation
442template <int VDIM, int MQ1, bool Transpose = false>
443inline MFEM_HOST_DEVICE void Eval2d(const int d1d, const int q1d,
444 real_t (&smem)[MQ1][MQ1],
445 const real_t (*B)[MQ1],
446 v_regs2d_t<VDIM, MQ1> &X,
447 v_regs2d_t<VDIM, MQ1> &Y)
448{
449 for (int c = 0; c < VDIM; c++)
450 {
451 Eval2d<MQ1, Transpose>(d1d, q1d, smem, B, X[c], Y[c]);
452 }
453}
454
455/// 2D vector transposed evaluation
456template <int VDIM, int MQ1>
457inline MFEM_HOST_DEVICE void EvalTranspose2d(const int d1d, const int q1d,
458 real_t (&smem)[MQ1][MQ1],
459 const real_t (*B)[MQ1],
460 v_regs2d_t<VDIM, MQ1> &X,
461 v_regs2d_t<VDIM, MQ1> &Y)
462{
463 Eval2d<VDIM, MQ1, true>(d1d, q1d, smem, B, X, Y);
464}
465
466/// 2D vector gradient, with component
467template <int VDIM, int DIM, int MQ1, bool Transpose = false>
468inline MFEM_HOST_DEVICE void Grad2d(const int d1d, const int q1d,
469 real_t (&smem)[MQ1][MQ1],
470 const real_t (*B)[MQ1],
471 const real_t (*G)[MQ1],
472 vd_regs2d_t<VDIM, DIM, MQ1> &X,
473 vd_regs2d_t<VDIM, DIM, MQ1> &Y,
474 const int c)
475{
476
477 for (int d = 0; d < DIM; d++)
478 {
479 const real_t (*Bx)[MQ1] = (d == 0) ? G : B;
480 const real_t (*By)[MQ1] = (d == 1) ? G : B;
481 Contract2d<Transpose>(d1d, q1d, smem, Bx, By, X[c][d], Y[c][d]);
482 }
483
484}
485
486/// 2D vector gradient
487template <int VDIM, int DIM, int MQ1, bool Transpose = false>
488inline MFEM_HOST_DEVICE void Grad2d(const int d1d, const int q1d,
489 real_t (&smem)[MQ1][MQ1],
490 const real_t (*B)[MQ1],
491 const real_t (*G)[MQ1],
492 vd_regs2d_t<VDIM, DIM, MQ1> &X,
493 vd_regs2d_t<VDIM, DIM, MQ1> &Y)
494{
495 for (int c = 0; c < VDIM; ++c)
496 {
497 Grad2d<VDIM, DIM, MQ1, Transpose>(d1d, q1d, smem, B, G, X, Y, c);
498 }
499}
500
501/// 2D vector transposed gradient
502template <int VDIM, int DIM, int MQ1>
503inline MFEM_HOST_DEVICE void GradTranspose2d(const int d1d, const int q1d,
504 real_t (&smem)[MQ1][MQ1],
505 const real_t (*B)[MQ1],
506 const real_t (*G)[MQ1],
507 vd_regs2d_t<VDIM, DIM, MQ1> &X,
508 vd_regs2d_t<VDIM, DIM, MQ1> &Y)
509{
510 constexpr bool Transpose = true;
511 Grad2d<VDIM, DIM, MQ1, Transpose>(d1d, q1d, smem, B, G, X, Y);
512}
513
514/// 2D scalar contraction, with component
515template <int VDIM, int DIM, int MQ1>
516inline MFEM_HOST_DEVICE void GradTranspose2d(const int d1d, const int q1d,
517 real_t (&smem)[MQ1][MQ1],
518 const real_t (*B)[MQ1],
519 const real_t (*G)[MQ1],
520 vd_regs2d_t<VDIM, DIM, MQ1> &X,
521 vd_regs2d_t<VDIM, DIM, MQ1> &Y,
522 const int c)
523{
524 constexpr bool Transpose = true;
525 Grad2d<VDIM, DIM, MQ1, Transpose>(d1d, q1d, smem, B, G, X, Y, c);
526}
527
528/// 3D scalar contraction, X direction
529template <bool Transpose, int MQ1>
530inline MFEM_HOST_DEVICE void ContractX3d(const int d1d, const int q1d,
531 real_t (&smem)[MQ1][MQ1],
532 const real_t (*B)[MQ1],
533 const s_regs3d_t<MQ1> &X,
534 s_regs3d_t<MQ1> &Y)
535{
536 for (int z = 0; z < d1d; ++z)
537 {
538 MFEM_FOREACH_THREAD_DIRECT(y, y, d1d)
539 {
540 MFEM_FOREACH_THREAD_DIRECT(x, x, (Transpose ? q1d : d1d))
541 {
542 smem[y][x] = X[z][y][x];
543 }
544 }
545 MFEM_SYNC_THREAD;
546 MFEM_FOREACH_THREAD_DIRECT(y, y, d1d)
547 {
548 MFEM_FOREACH_THREAD_DIRECT(x, x, (Transpose ? d1d : q1d))
549 {
550 real_t u = 0.0;
551 for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
552 {
553 u += (Transpose ? B[x][k] : B[k][x]) * smem[y][k];
554 }
555 Y[z][y][x] = u;
556 }
557 }
558 MFEM_SYNC_THREAD;
559 }
560}
561
562/// 3D scalar contraction, Y direction
563template <bool Transpose, int MQ1>
564inline MFEM_HOST_DEVICE void ContractY3d(const int d1d, const int q1d,
565 real_t (&smem)[MQ1][MQ1],
566 const real_t (*B)[MQ1],
567 const s_regs3d_t<MQ1> &X,
568 s_regs3d_t<MQ1> &Y)
569{
570 for (int z = 0; z < d1d; ++z)
571 {
572 MFEM_FOREACH_THREAD_DIRECT(y, y, (Transpose ? q1d : d1d))
573 {
574 MFEM_FOREACH_THREAD_DIRECT(x, x, q1d) { smem[y][x] = X[z][y][x]; }
575 }
576 MFEM_SYNC_THREAD;
577 MFEM_FOREACH_THREAD_DIRECT(y, y, (Transpose ? d1d : q1d))
578 {
579 MFEM_FOREACH_THREAD_DIRECT(x, x, q1d)
580 {
581 real_t u = 0.0;
582 for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
583 {
584 u += (Transpose ? B[y][k] : B[k][y]) * smem[k][x];
585 }
586 Y[z][y][x] = u;
587 }
588 }
589 MFEM_SYNC_THREAD;
590 }
591}
592
593/// 3D scalar contraction, Z direction
594template <bool Transpose, int MQ1>
595inline MFEM_HOST_DEVICE void ContractZ3d(const int d1d, const int q1d,
596 const real_t (*B)[MQ1],
597 const s_regs3d_t<MQ1> &X,
598 s_regs3d_t<MQ1> &Y)
599{
600 for (int z = 0; z < (Transpose ? d1d : q1d); ++z)
601 {
602 MFEM_FOREACH_THREAD_DIRECT(y, y, q1d)
603 {
604 MFEM_FOREACH_THREAD_DIRECT(x, x, q1d)
605 {
606 real_t u = 0.0;
607 for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
608 {
609 u += (Transpose ? B[z][k] : B[k][z]) * X[k][y][x];
610 }
611 Y[z][y][x] = u;
612 }
613 }
614 }
615}
616
617/// 3D scalar contraction: X, Y & Z directions
618template <bool Transpose, int MQ1>
619inline MFEM_HOST_DEVICE void Contract3d(const int d1d, const int q1d,
620 real_t (&smem)[MQ1][MQ1],
621 const real_t (*Bx)[MQ1],
622 const real_t (*By)[MQ1],
623 const real_t (*Bz)[MQ1],
624 s_regs3d_t<MQ1> &X,
625 s_regs3d_t<MQ1> &Y)
626{
627 if (!Transpose)
628 {
629 ContractX3d<false>(d1d, q1d, smem, Bx, X, Y);
630 ContractY3d<false>(d1d, q1d, smem, By, Y, X);
631 ContractZ3d<false>(d1d, q1d, Bz, X, Y);
632 }
633 else
634 {
635 ContractZ3d<true>(d1d, q1d, Bz, X, Y);
636 ContractY3d<true>(d1d, q1d, smem, By, Y, X);
637 ContractX3d<true>(d1d, q1d, smem, Bx, X, Y);
638 }
639}
640
641/// 3D scalar evaluation
642template <int MQ1, bool Transpose = false>
643inline MFEM_HOST_DEVICE void Eval3d(const int d1d, const int q1d,
644 real_t (&smem)[MQ1][MQ1],
645 const real_t (*B)[MQ1],
646 s_regs3d_t<MQ1> &X,
647 s_regs3d_t<MQ1> &Y)
648{
649 Contract3d<Transpose>(d1d, q1d, smem, B, B, B, X, Y);
650}
651
652/// 3D vector evaluation
653template <int VDIM, int MQ1, bool Transpose = false>
654inline MFEM_HOST_DEVICE void Eval3d(const int d1d, const int q1d,
655 real_t (&smem)[MQ1][MQ1],
656 const real_t (*B)[MQ1],
657 v_regs3d_t<VDIM, MQ1> &X,
658 v_regs3d_t<VDIM, MQ1> &Y)
659{
660 for (int c = 0; c < VDIM; c++)
661 {
662 Eval3d<MQ1, Transpose>(d1d, q1d, smem, B, X[c], Y[c]);
663 }
664}
665
666/// 3D vector transposed evaluation
667template <int VDIM, int MQ1>
668inline MFEM_HOST_DEVICE void EvalTranspose3d(const int d1d, const int q1d,
669 real_t (&smem)[MQ1][MQ1],
670 const real_t (*B)[MQ1],
671 v_regs3d_t<VDIM, MQ1> &X,
672 v_regs3d_t<VDIM, MQ1> &Y)
673{
674 Eval3d<VDIM, MQ1, true>(d1d, q1d, smem, B, X, Y);
675}
676
677/// 3D vector gradient, with component
678template <int VDIM, int DIM, int MQ1, bool Transpose = false>
679inline MFEM_HOST_DEVICE void Grad3d(const int d1d, const int q1d,
680 real_t (&smem)[MQ1][MQ1],
681 const real_t (*B)[MQ1],
682 const real_t (*G)[MQ1],
683 vd_regs3d_t<VDIM, DIM, MQ1> &X,
684 vd_regs3d_t<VDIM, DIM, MQ1> &Y,
685 const int c)
686{
687 for (int d = 0; d < DIM; d++)
688 {
689 const real_t (*Bx)[MQ1] = (d == 0) ? G : B;
690 const real_t (*By)[MQ1] = (d == 1) ? G : B;
691 const real_t (*Bz)[MQ1] = (d == 2) ? G : B;
692 Contract3d<Transpose>(d1d, q1d, smem, Bx, By, Bz, X[c][d], Y[c][d]);
693 }
694}
695
696/// 3D vector gradient
697template <int VDIM, int DIM, int MQ1, bool Transpose = false>
698inline MFEM_HOST_DEVICE void Grad3d(const int d1d, const int q1d,
699 real_t (&smem)[MQ1][MQ1],
700 const real_t (*B)[MQ1],
701 const real_t (*G)[MQ1],
702 vd_regs3d_t<VDIM, DIM, MQ1> &X,
703 vd_regs3d_t<VDIM, DIM, MQ1> &Y)
704{
705 for (int c = 0; c < VDIM; c++)
706 {
707 Grad3d<VDIM, DIM, MQ1, Transpose>(d1d, q1d, smem, B, G, X, Y, c);
708 }
709}
710
711/// 3D vector transposed gradient
712template <int VDIM, int DIM, int MQ1>
713inline MFEM_HOST_DEVICE void GradTranspose3d(const int d1d, const int q1d,
714 real_t (&smem)[MQ1][MQ1],
715 const real_t (*B)[MQ1],
716 const real_t (*G)[MQ1],
717 vd_regs3d_t<VDIM, DIM, MQ1> &X,
718 vd_regs3d_t<VDIM, DIM, MQ1> &Y)
719{
720 Grad3d<VDIM, DIM, MQ1, true>(d1d, q1d, smem, B, G, X, Y);
721}
722
723/// 3D vector transposed gradient, with component
724template <int VDIM, int DIM, int MQ1>
725inline MFEM_HOST_DEVICE void GradTranspose3d(const int d1d, const int q1d,
726 real_t (&smem)[MQ1][MQ1],
727 const real_t (*B)[MQ1],
728 const real_t (*G)[MQ1],
729 vd_regs3d_t<VDIM, DIM, MQ1> &X,
730 vd_regs3d_t<VDIM, DIM, MQ1> &Y,
731 const int c)
732{
733 Grad3d<VDIM, DIM, MQ1, true>(d1d, q1d, smem, B, G, X, Y, c);
734}
735
736/// Load B1d matrix into shared memory
737template<int MD1, int MQ1>
738MFEM_HOST_DEVICE inline void LoadB(const int D1D, const int Q1D,
739 const ConstDeviceMatrix &b,
740 real_t (&sB)[MQ1*MD1])
741{
742 const int tidz = MFEM_THREAD_ID(z);
743 DeviceMatrix B(sB, D1D, Q1D);
744
745 if (tidz == 0)
746 {
747 MFEM_FOREACH_THREAD(d,y,D1D)
748 {
749 MFEM_FOREACH_THREAD(q,x,Q1D)
750 {
751 B(d,q) = b(q,d);
752 }
753 }
754 }
755 MFEM_SYNC_THREAD;
756}
757
758/// Load Bt1d matrix into shared memory
759template<int MD1, int MQ1>
760MFEM_HOST_DEVICE inline void LoadBt(const int D1D, const int Q1D,
761 const ConstDeviceMatrix &b,
762 real_t (&sB)[MQ1*MD1])
763{
764 const int tidz = MFEM_THREAD_ID(z);
765 DeviceMatrix Bt(sB, Q1D, D1D);
766
767 if (tidz == 0)
768 {
769 MFEM_FOREACH_THREAD(d,y,D1D)
770 {
771 MFEM_FOREACH_THREAD(q,x,Q1D)
772 {
773 Bt(q,d) = b(q,d);
774 }
775 }
776 }
777 MFEM_SYNC_THREAD;
778}
779
780/// Load B1d & G1d matrices into shared memory
781template<int MD1, int MQ1>
782MFEM_HOST_DEVICE inline void LoadBG(const int D1D, const int Q1D,
783 const ConstDeviceMatrix &b,
784 const ConstDeviceMatrix &g,
785 real_t (&sBG)[2][MQ1*MD1])
786{
787 const int tidz = MFEM_THREAD_ID(z);
788 DeviceMatrix B(sBG[0], D1D, Q1D);
789 DeviceMatrix G(sBG[1], D1D, Q1D);
790
791 if (tidz == 0)
792 {
793 MFEM_FOREACH_THREAD(d,y,D1D)
794 {
795 MFEM_FOREACH_THREAD(q,x,Q1D)
796 {
797 B(d,q) = b(q,d);
798 G(d,q) = g(q,d);
799 }
800 }
801 }
802 MFEM_SYNC_THREAD;
803}
804
805/// Load Bt1d & Gt1d matrices into shared memory
806template<int MD1, int MQ1>
807MFEM_HOST_DEVICE inline void LoadBGt(const int D1D, const int Q1D,
808 const ConstDeviceMatrix &b,
809 const ConstDeviceMatrix &g,
810 real_t (&sBG)[2][MQ1*MD1])
811{
812 const int tidz = MFEM_THREAD_ID(z);
813 DeviceMatrix Bt(sBG[0], Q1D, D1D);
814 DeviceMatrix Gt(sBG[1], Q1D, D1D);
815
816 if (tidz == 0)
817 {
818 MFEM_FOREACH_THREAD(d,y,D1D)
819 {
820 MFEM_FOREACH_THREAD(q,x,Q1D)
821 {
822 Bt(q,d) = b(q,d);
823 Gt(q,d) = g(q,d);
824 }
825 }
826 }
827 MFEM_SYNC_THREAD;
828}
829
830/// Load 2D input scalar into given DeviceMatrix
831MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
832 const DeviceTensor<3, const real_t> &x,
833 DeviceMatrix &DD)
834{
835 MFEM_FOREACH_THREAD(dy,y,D1D)
836 {
837 MFEM_FOREACH_THREAD(dx,x,D1D)
838 {
839 DD(dx,dy) = x(dx,dy,e);
840 }
841 }
842 MFEM_SYNC_THREAD;
843}
844
845
846/// Load 2D input scalar into shared memory
847template<int MD1, int NBZ>
848MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
849 const DeviceTensor<3, const real_t> &x,
850 real_t (&sX)[NBZ][MD1*MD1])
851{
852 const int tidz = MFEM_THREAD_ID(z);
853 DeviceMatrix X(sX[tidz], D1D, D1D);
854 LoadX(e, D1D, x, X);
855}
856
857/// Load 2D input scalar into shared memory, with comp
858MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
859 const DeviceTensor<4, const real_t> &x,
860 DeviceMatrix &DD)
861{
862 MFEM_FOREACH_THREAD(dy,y,D1D)
863 {
864 MFEM_FOREACH_THREAD(dx,x,D1D)
865 {
866 DD(dx,dy) = x(dx,dy,c,e);
867 }
868 }
869 MFEM_SYNC_THREAD;
870}
871
872template<int MD1, int NBZ>
873MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
874 const DeviceTensor<4, const real_t> &x,
875 real_t (&sm)[NBZ][MD1*MD1])
876{
877 const int tidz = MFEM_THREAD_ID(z);
878 DeviceMatrix DD(sm[tidz], D1D, D1D);
879 LoadX(e,D1D,c,x,DD);
880}
881
882/// 2D Scalar Evaluation, 1/2
883MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
884 ConstDeviceMatrix &B,
885 DeviceMatrix &DD,
886 DeviceMatrix &DQ)
887{
888 MFEM_FOREACH_THREAD(dy,y,D1D)
889 {
890 MFEM_FOREACH_THREAD(qx,x,Q1D)
891 {
892 real_t u = 0.0;
893 for (int dx = 0; dx < D1D; ++dx)
894 {
895 u += B(dx,qx) * DD(dx,dy);
896 }
897 DQ(dy,qx) = u;
898 }
899 }
900 MFEM_SYNC_THREAD;
901}
902
903template<int MD1, int MQ1, int NBZ>
904MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
905 const real_t (&sB)[MQ1*MD1],
906 real_t (&sDD)[NBZ][MD1*MD1],
907 real_t (&sDQ)[NBZ][MD1*MQ1])
908{
909 const int tidz = MFEM_THREAD_ID(z);
910 ConstDeviceMatrix B(sB, D1D, Q1D);
911 DeviceMatrix DD(sDD[tidz], D1D, D1D);
912 DeviceMatrix DQ(sDQ[tidz], D1D, Q1D);
913 EvalX(D1D,Q1D,B,DD,DQ);
914}
915
916/// 2D Scalar Evaluation, 2/2
917MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
918 ConstDeviceMatrix &B,
919 DeviceMatrix &DQ,
920 DeviceMatrix &QQ)
921{
922 MFEM_FOREACH_THREAD(qy,y,Q1D)
923 {
924 MFEM_FOREACH_THREAD(qx,x,Q1D)
925 {
926 real_t u = 0.0;
927 for (int dy = 0; dy < D1D; ++dy)
928 {
929 u += DQ(dy,qx) * B(dy,qy);
930 }
931 QQ(qx,qy) = u;
932 }
933 }
934 MFEM_SYNC_THREAD;
935}
936
937template<int MD1, int MQ1, int NBZ>
938MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
939 const real_t (&sB)[MQ1*MD1],
940 real_t (&sDQ)[NBZ][MD1*MQ1],
941 real_t (&sQQ)[NBZ][MQ1*MQ1])
942{
943 const int tidz = MFEM_THREAD_ID(z);
944 ConstDeviceMatrix B(sB, D1D, Q1D);
945 DeviceMatrix DQ(sDQ[tidz], D1D, Q1D);
946 DeviceMatrix QQ(sQQ[tidz], Q1D, Q1D);
947 EvalY(D1D,Q1D,B,DQ,QQ);
948}
949
950/// Pull 2D Scalar Evaluation
951MFEM_HOST_DEVICE inline void PullEval(const int qx, const int qy,
952 DeviceMatrix &QQ,
953 real_t &P)
954{
955 P = QQ(qx,qy);
956}
957
958template<int MQ1, int NBZ>
959MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
960 const int qx, const int qy,
961 real_t (&sQQ)[NBZ][MQ1*MQ1],
962 real_t &P)
963{
964 const int tidz = MFEM_THREAD_ID(z);
965 DeviceMatrix QQ(sQQ[tidz], Q1D, Q1D);
966 PullEval(qx,qy,QQ,P);
967}
968
969/// Load 2D input vector into shared memory
970template<int MD1, int NBZ>
971MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
972 const DeviceTensor<4, const real_t> &X,
973 real_t (&sX)[2][NBZ][MD1*MD1])
974{
975 const int tidz = MFEM_THREAD_ID(z);
976 DeviceMatrix X0(sX[0][tidz], D1D, D1D);
977 DeviceMatrix X1(sX[1][tidz], D1D, D1D);
978
979 MFEM_FOREACH_THREAD(dy,y,D1D)
980 {
981 MFEM_FOREACH_THREAD(dx,x,D1D)
982 {
983 X0(dx,dy) = X(dx,dy,0,e);
984 X1(dx,dy) = X(dx,dy,1,e);
985 }
986 }
987 MFEM_SYNC_THREAD;
988}
989
990/// 2D Evaluation, 1/2 (only B)
991template<int MD1, int MQ1, int NBZ>
992MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
993 const real_t (&sB)[MQ1*MD1],
994 const real_t (&sX)[2][NBZ][MD1*MD1],
995 real_t (&sDQ)[2][NBZ][MD1*MQ1])
996{
997 const int tidz = MFEM_THREAD_ID(z);
998 ConstDeviceMatrix B(sB, D1D, Q1D);
999 ConstDeviceMatrix X0(sX[0][tidz], D1D, D1D);
1000 ConstDeviceMatrix X1(sX[1][tidz], D1D, D1D);
1001 DeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
1002 DeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
1003
1004 MFEM_FOREACH_THREAD(dy,y,D1D)
1005 {
1006 MFEM_FOREACH_THREAD(qx,x,Q1D)
1007 {
1008 real_t u[2] = {0.0, 0.0};
1009 for (int dx = 0; dx < D1D; ++dx)
1010 {
1011 const real_t xx = X0(dx,dy);
1012 const real_t xy = X1(dx,dy);
1013 u[0] += B(dx,qx) * xx;
1014 u[1] += B(dx,qx) * xy;
1015 }
1016 DQ0(qx,dy) = u[0];
1017 DQ1(qx,dy) = u[1];
1018 }
1019 }
1020 MFEM_SYNC_THREAD;
1021}
1022
1023/// 2D Evaluation, 2/2 (only B)
1024template<int MD1, int MQ1, int NBZ>
1025MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
1026 const real_t (&sB)[MQ1*MD1],
1027 const real_t (&sDQ)[2][NBZ][MD1*MQ1],
1028 real_t (&sQQ)[2][NBZ][MQ1*MQ1])
1029{
1030 const int tidz = MFEM_THREAD_ID(z);
1031 ConstDeviceMatrix B(sB, D1D, Q1D);
1032 ConstDeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
1033 ConstDeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
1034 DeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
1035 DeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
1036
1037 MFEM_FOREACH_THREAD(qy,y,Q1D)
1038 {
1039 MFEM_FOREACH_THREAD(qx,x,Q1D)
1040 {
1041 real_t u[2] = {0.0, 0.0};
1042 for (int dy = 0; dy < D1D; ++dy)
1043 {
1044 u[0] += DQ0(qx,dy) * B(dy,qy);
1045 u[1] += DQ1(qx,dy) * B(dy,qy);
1046 }
1047 QQ0(qx,qy) = u[0];
1048 QQ1(qx,qy) = u[1];
1049 }
1050 }
1051 MFEM_SYNC_THREAD;
1052}
1053
1054/// Pull 2D Evaluation
1055template<int MQ1, int NBZ>
1056MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
1057 const int qx, const int qy,
1058 const real_t (&sQQ)[2][NBZ][MQ1*MQ1],
1059 real_t (&P)[2])
1060{
1061 const int tidz = MFEM_THREAD_ID(z);
1062 ConstDeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
1063 ConstDeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
1064
1065 P[0] = QQ0(qx,qy);
1066 P[1] = QQ1(qx,qy);
1067}
1068
1069/// Push 2D Evaluation
1070template<int MQ1, int NBZ>
1071MFEM_HOST_DEVICE inline void PushEval(const int Q1D,
1072 const int qx, const int qy,
1073 const real_t *P,
1074 real_t (&sQQ)[2][NBZ][MQ1*MQ1])
1075{
1076 const int tidz = MFEM_THREAD_ID(z);
1077 DeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
1078 DeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
1079
1080 QQ0(qx,qy) = P[0];
1081 QQ1(qx,qy) = P[1];
1082}
1083
1084/// 2D Transposed evaluation, 1/2
1085template<int MD1, int MQ1, int NBZ>
1086MFEM_HOST_DEVICE inline void EvalXt(const int D1D, const int Q1D,
1087 const real_t (&sB)[MQ1*MD1],
1088 const real_t (&sQQ)[2][NBZ][MQ1*MQ1],
1089 real_t (&sDQ)[2][NBZ][MD1*MQ1])
1090{
1091 const int tidz = MFEM_THREAD_ID(z);
1092 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1093 ConstDeviceMatrix QQ0(sQQ[0][tidz], Q1D, Q1D);
1094 ConstDeviceMatrix QQ1(sQQ[1][tidz], Q1D, Q1D);
1095 DeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
1096 DeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
1097
1098 MFEM_FOREACH_THREAD(qy,y,Q1D)
1099 {
1100 MFEM_FOREACH_THREAD(dx,x,D1D)
1101 {
1102 real_t u[2] = {0.0, 0.0};
1103 for (int qx = 0; qx < Q1D; ++qx)
1104 {
1105 u[0] += QQ0(qx,qy) * Bt(qx,dx);
1106 u[1] += QQ1(qx,qy) * Bt(qx,dx);
1107 }
1108 DQ0(qy,dx) = u[0];
1109 DQ1(qy,dx) = u[1];
1110 }
1111 }
1112 MFEM_SYNC_THREAD;
1113}
1114
1115/// 2D Transposed evaluation, 2/2
1116template<int MD1, int MQ1, int NBZ>
1117MFEM_HOST_DEVICE inline void EvalYt(const int D1D, const int Q1D,
1118 const real_t (&sB)[MQ1*MD1],
1119 const real_t (&sDQ)[2][NBZ][MD1*MQ1],
1120 const DeviceTensor<4> &Y, // output
1121 const int e)
1122{
1123 const int tidz = MFEM_THREAD_ID(z);
1124 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1125 ConstDeviceMatrix DQ0(sDQ[0][tidz], Q1D, D1D);
1126 ConstDeviceMatrix DQ1(sDQ[1][tidz], Q1D, D1D);
1127
1128 MFEM_FOREACH_THREAD(dy,y,D1D)
1129 {
1130 MFEM_FOREACH_THREAD(dx,x,D1D)
1131 {
1132 real_t u[2] = {0.0, 0.0};
1133 for (int qy = 0; qy < Q1D; ++qy)
1134 {
1135 u[0] += Bt(qy,dy) * DQ0(qy,dx);
1136 u[1] += Bt(qy,dy) * DQ1(qy,dx);
1137 }
1138 Y(dx,dy,0,e) += u[0];
1139 Y(dx,dy,1,e) += u[1];
1140 }
1141 }
1142 MFEM_SYNC_THREAD;
1143}
1144
1145/// 2D Gradient, 1/2
1146template<int MD1, int MQ1, int NBZ>
1147MFEM_HOST_DEVICE inline void GradX(const int D1D, const int Q1D,
1148 const real_t (&sBG)[2][MQ1*MD1],
1149 const real_t (&sX)[2][NBZ][MD1*MD1],
1150 real_t (&sDQ)[4][NBZ][MD1*MQ1])
1151{
1152 const int tidz = MFEM_THREAD_ID(z);
1153 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1154 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1155 ConstDeviceMatrix X0(sX[0][tidz], D1D, D1D);
1156 ConstDeviceMatrix X1(sX[1][tidz], D1D, D1D);
1157 DeviceMatrix X0B(sDQ[0][tidz], Q1D, D1D);
1158 DeviceMatrix X0G(sDQ[1][tidz], Q1D, D1D);
1159 DeviceMatrix X1B(sDQ[2][tidz], Q1D, D1D);
1160 DeviceMatrix X1G(sDQ[3][tidz], Q1D, D1D);
1161
1162 MFEM_FOREACH_THREAD(dy,y,D1D)
1163 {
1164 MFEM_FOREACH_THREAD(qx,x,Q1D)
1165 {
1166 real_t u[2] = {0.0, 0.0};
1167 real_t v[2] = {0.0, 0.0};
1168 for (int dx = 0; dx < D1D; ++dx)
1169 {
1170 const real_t Bx = B(dx,qx);
1171 const real_t Gx = G(dx,qx);
1172 const real_t x0 = X0(dx,dy);
1173 const real_t x1 = X1(dx,dy);
1174 u[0] += Bx * x0;
1175 v[0] += Gx * x0;
1176 u[1] += Bx * x1;
1177 v[1] += Gx * x1;
1178 }
1179 X0B(qx,dy) = u[0];
1180 X0G(qx,dy) = v[0];
1181 X1B(qx,dy) = u[1];
1182 X1G(qx,dy) = v[1];
1183 }
1184 }
1185 MFEM_SYNC_THREAD;
1186}
1187
1188/// 2D Gradient, 2/2
1189template<int MD1, int MQ1, int NBZ>
1190MFEM_HOST_DEVICE inline void GradY(const int D1D, const int Q1D,
1191 const real_t (&sBG)[2][MQ1*MD1],
1192 const real_t (&sDQ)[4][NBZ][MD1*MQ1],
1193 real_t (&sQQ)[4][NBZ][MQ1*MQ1])
1194{
1195 const int tidz = MFEM_THREAD_ID(z);
1196 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1197 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1198 ConstDeviceMatrix X0B(sDQ[0][tidz], Q1D, D1D);
1199 ConstDeviceMatrix X0G(sDQ[1][tidz], Q1D, D1D);
1200 ConstDeviceMatrix X1B(sDQ[2][tidz], Q1D, D1D);
1201 ConstDeviceMatrix X1G(sDQ[3][tidz], Q1D, D1D);
1202 DeviceMatrix X0GB(sQQ[0][tidz], Q1D, Q1D);
1203 DeviceMatrix X0BG(sQQ[1][tidz], Q1D, Q1D);
1204 DeviceMatrix X1GB(sQQ[2][tidz], Q1D, Q1D);
1205 DeviceMatrix X1BG(sQQ[3][tidz], Q1D, Q1D);
1206
1207 MFEM_FOREACH_THREAD(qy,y,Q1D)
1208 {
1209 MFEM_FOREACH_THREAD(qx,x,Q1D)
1210 {
1211 real_t u[2] = {0.0, 0.0};
1212 real_t v[2] = {0.0, 0.0};
1213 for (int dy = 0; dy < D1D; ++dy)
1214 {
1215 const real_t By = B(dy,qy);
1216 const real_t Gy = G(dy,qy);
1217 u[0] += X0G(qx,dy) * By;
1218 v[0] += X0B(qx,dy) * Gy;
1219 u[1] += X1G(qx,dy) * By;
1220 v[1] += X1B(qx,dy) * Gy;
1221 }
1222 X0GB(qx,qy) = u[0];
1223 X0BG(qx,qy) = v[0];
1224 X1GB(qx,qy) = u[1];
1225 X1BG(qx,qy) = v[1];
1226 }
1227 }
1228 MFEM_SYNC_THREAD;
1229}
1230
1231/// Pull 2D Gradient
1232template<int MQ1, int NBZ>
1233MFEM_HOST_DEVICE inline void PullGrad(const int Q1D,
1234 const int qx, const int qy,
1235 const real_t (&sQQ)[4][NBZ][MQ1*MQ1],
1236 real_t *Jpr)
1237{
1238 const int tidz = MFEM_THREAD_ID(z);
1239 ConstDeviceMatrix X0GB(sQQ[0][tidz], Q1D, Q1D);
1240 ConstDeviceMatrix X0BG(sQQ[1][tidz], Q1D, Q1D);
1241 ConstDeviceMatrix X1GB(sQQ[2][tidz], Q1D, Q1D);
1242 ConstDeviceMatrix X1BG(sQQ[3][tidz], Q1D, Q1D);
1243
1244 Jpr[0] = X0GB(qx,qy);
1245 Jpr[1] = X1GB(qx,qy);
1246 Jpr[2] = X0BG(qx,qy);
1247 Jpr[3] = X1BG(qx,qy);
1248}
1249
1250/// Push 2D Gradient
1251template<int MQ1, int NBZ>
1252MFEM_HOST_DEVICE inline void PushGrad(const int Q1D,
1253 const int qx, const int qy,
1254 const real_t *A,
1255 real_t (&sQQ)[4][NBZ][MQ1*MQ1])
1256{
1257 const int tidz = MFEM_THREAD_ID(z);
1258 DeviceMatrix X0GB(sQQ[0][tidz], Q1D, Q1D);
1259 DeviceMatrix X0BG(sQQ[1][tidz], Q1D, Q1D);
1260 DeviceMatrix X1GB(sQQ[2][tidz], Q1D, Q1D);
1261 DeviceMatrix X1BG(sQQ[3][tidz], Q1D, Q1D);
1262
1263 X0GB(qx,qy) = A[0];
1264 X1GB(qx,qy) = A[2];
1265 X0BG(qx,qy) = A[1];
1266 X1BG(qx,qy) = A[3];
1267}
1268
1269/// 2D Transposed gradient, 1/2
1270template<int MD1, int MQ1, int NBZ>
1271MFEM_HOST_DEVICE inline void GradYt(const int D1D, const int Q1D,
1272 const real_t (&sBG)[2][MQ1*MD1],
1273 const real_t (&GQ)[4][NBZ][MQ1*MQ1],
1274 real_t (&GD)[4][NBZ][MD1*MQ1])
1275{
1276 const int tidz = MFEM_THREAD_ID(z);
1277 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
1278 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
1279 ConstDeviceMatrix QQx0(GQ[0][tidz], Q1D, Q1D);
1280 ConstDeviceMatrix QQx1(GQ[1][tidz], Q1D, Q1D);
1281 ConstDeviceMatrix QQy0(GQ[2][tidz], Q1D, Q1D);
1282 ConstDeviceMatrix QQy1(GQ[3][tidz], Q1D, Q1D);
1283 DeviceMatrix DQxB(GD[0][tidz], Q1D, D1D);
1284 DeviceMatrix DQxG(GD[1][tidz], Q1D, D1D);
1285 DeviceMatrix DQyB(GD[2][tidz], Q1D, D1D);
1286 DeviceMatrix DQyG(GD[3][tidz], Q1D, D1D);
1287
1288 MFEM_FOREACH_THREAD(qy,y,Q1D)
1289 {
1290 MFEM_FOREACH_THREAD(dx,x,D1D)
1291 {
1292 real_t u[2] = {0.0, 0.0};
1293 real_t v[2] = {0.0, 0.0};
1294 for (int qx = 0; qx < Q1D; ++qx)
1295 {
1296 u[0] += Gt(qx,dx) * QQx0(qx,qy);
1297 u[1] += Gt(qx,dx) * QQy0(qx,qy);
1298 v[0] += Bt(qx,dx) * QQx1(qx,qy);
1299 v[1] += Bt(qx,dx) * QQy1(qx,qy);
1300 }
1301 DQxB(qy,dx) = u[0];
1302 DQyB(qy,dx) = u[1];
1303 DQxG(qy,dx) = v[0];
1304 DQyG(qy,dx) = v[1];
1305 }
1306 }
1307 MFEM_SYNC_THREAD;
1308}
1309
1310/// 2D Transposed gradient, 2/2
1311template<int MD1, int MQ1, int NBZ>
1312MFEM_HOST_DEVICE inline void GradXt(const int D1D, const int Q1D,
1313 const real_t (&sBG)[2][MQ1*MD1],
1314 const real_t (&GD)[4][NBZ][MD1*MQ1],
1315 const DeviceTensor<4> &Y, // output
1316 const int e)
1317{
1318 const int tidz = MFEM_THREAD_ID(z);
1319 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
1320 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
1321 ConstDeviceMatrix DQxB(GD[0][tidz], Q1D, D1D);
1322 ConstDeviceMatrix DQxG(GD[1][tidz], Q1D, D1D);
1323 ConstDeviceMatrix DQyB(GD[2][tidz], Q1D, D1D);
1324 ConstDeviceMatrix DQyG(GD[3][tidz], Q1D, D1D);
1325
1326 MFEM_FOREACH_THREAD(dy,y,D1D)
1327 {
1328 MFEM_FOREACH_THREAD(dx,x,D1D)
1329 {
1330 real_t u[2] = {0.0, 0.0};
1331 real_t v[2] = {0.0, 0.0};
1332 for (int qy = 0; qy < Q1D; ++qy)
1333 {
1334 u[0] += DQxB(qy,dx) * Bt(qy,dy);
1335 u[1] += DQyB(qy,dx) * Bt(qy,dy);
1336 v[0] += DQxG(qy,dx) * Gt(qy,dy);
1337 v[1] += DQyG(qy,dx) * Gt(qy,dy);
1338 }
1339 Y(dx,dy,0,e) += u[0] + v[0];
1340 Y(dx,dy,1,e) += u[1] + v[1];
1341 }
1342 }
1343 MFEM_SYNC_THREAD;
1344}
1345
1346/// Load 3D scalar input vector into shared memory
1347MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
1348 const DeviceTensor<4, const real_t> &x,
1349 DeviceCube &X)
1350{
1351 MFEM_FOREACH_THREAD(dz,z,D1D)
1352 {
1353 MFEM_FOREACH_THREAD(dy,y,D1D)
1354 {
1355 MFEM_FOREACH_THREAD(dx,x,D1D)
1356 {
1357 X(dx,dy,dz) = x(dx,dy,dz,e);
1358 }
1359 }
1360 }
1361 MFEM_SYNC_THREAD;
1362}
1363
1364template<int MD1>
1365MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
1366 const DeviceTensor<4, const real_t> &x,
1367 real_t (&sm)[MD1*MD1*MD1])
1368{
1369 DeviceCube X(sm, D1D,D1D,D1D);
1370 LoadX(e,D1D,x,X);
1371}
1372
1373/// Load 3D scalar input vector into shared memory, with comp & DeviceTensor
1374MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
1375 const DeviceTensor<5, const real_t> &x,
1376 DeviceTensor<3> &X)
1377{
1378 MFEM_FOREACH_THREAD(dz,z,D1D)
1379 {
1380 MFEM_FOREACH_THREAD(dy,y,D1D)
1381 {
1382 MFEM_FOREACH_THREAD(dx,x,D1D)
1383 {
1384 X(dx,dy,dz) = x(dx,dy,dz,c,e);
1385 }
1386 }
1387 }
1388 MFEM_SYNC_THREAD;
1389}
1390
1391/// Load 3D scalar input vector into shared memory, with comp & pointer
1392template<int MD1>
1393MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
1394 const DeviceTensor<5, const real_t> &x,
1395 real_t (&sm)[MD1*MD1*MD1])
1396{
1397 DeviceCube X(sm, D1D, D1D, D1D);
1398 return LoadX<MD1>(e,D1D,c,x,X);
1399}
1400
1401/// 3D Scalar Evaluation, 1/3
1402MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
1403 ConstDeviceMatrix &B,
1404 const DeviceCube &DDD,
1405 DeviceCube &DDQ)
1406{
1407 MFEM_FOREACH_THREAD(dz,z,D1D)
1408 {
1409 MFEM_FOREACH_THREAD(dy,y,D1D)
1410 {
1411 MFEM_FOREACH_THREAD(qx,x,Q1D)
1412 {
1413 real_t u = 0.0;
1414 for (int dx = 0; dx < D1D; ++dx)
1415 {
1416 const real_t Bx = B(dx,qx);
1417 u += Bx * DDD(dx,dy,dz);
1418 }
1419 DDQ(dz,dy,qx) = u;
1420 }
1421 }
1422 }
1423 MFEM_SYNC_THREAD;
1424}
1425
1426template<int MD1, int MQ1>
1427MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
1428 const real_t (&sB)[MQ1*MD1],
1429 const real_t (&sDDD)[MD1*MD1*MD1],
1430 real_t (&sDDQ)[MD1*MD1*MQ1])
1431{
1432 ConstDeviceMatrix B(sB, D1D, Q1D);
1433 const DeviceCube DDD(sDDD, D1D, D1D, D1D);
1434 DeviceCube DDQ(sDDQ, Q1D, D1D, D1D);
1435 EvalX(D1D,Q1D,B,DDD,DDQ);
1436}
1437
1438/// 3D Scalar Evaluation, 2/3
1439MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
1440 ConstDeviceMatrix &B,
1441 const DeviceCube &DDQ,
1442 DeviceCube &DQQ)
1443{
1444 MFEM_FOREACH_THREAD(dz,z,D1D)
1445 {
1446 MFEM_FOREACH_THREAD(qy,y,Q1D)
1447 {
1448 MFEM_FOREACH_THREAD(qx,x,Q1D)
1449 {
1450 real_t u = 0.0;
1451 for (int dy = 0; dy < D1D; ++dy)
1452 {
1453 const real_t By = B(dy,qy);
1454 u += DDQ(dz,dy,qx) * By;
1455 }
1456 DQQ(dz,qy,qx) = u;
1457 }
1458 }
1459 }
1460 MFEM_SYNC_THREAD;
1461}
1462
1463template<int MD1, int MQ1>
1464MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
1465 const real_t (&sB)[MQ1*MD1],
1466 const real_t (&sDDQ)[MD1*MD1*MQ1],
1467 real_t (&sDQQ)[MD1*MQ1*MQ1])
1468{
1469 ConstDeviceMatrix B(sB, D1D, Q1D);
1470 const DeviceCube DDQ(sDDQ, Q1D, D1D, D1D);
1471 DeviceCube DQQ(sDQQ, Q1D, Q1D, D1D);
1472 EvalY(D1D,Q1D,B,DDQ,DQQ);
1473}
1474
1475/// 3D Scalar Evaluation, 3/3
1476MFEM_HOST_DEVICE inline void EvalZ(const int D1D, const int Q1D,
1477 ConstDeviceMatrix &B,
1478 const DeviceCube &DQQ,
1479 DeviceCube &QQQ)
1480{
1481 MFEM_FOREACH_THREAD(qz,z,Q1D)
1482 {
1483 MFEM_FOREACH_THREAD(qy,y,Q1D)
1484 {
1485 MFEM_FOREACH_THREAD(qx,x,Q1D)
1486 {
1487 real_t u = 0.0;
1488 for (int dz = 0; dz < D1D; ++dz)
1489 {
1490 const real_t Bz = B(dz,qz);
1491 u += DQQ(dz,qy,qx) * Bz;
1492 }
1493 QQQ(qz,qy,qx) = u;
1494 }
1495 }
1496 }
1497 MFEM_SYNC_THREAD;
1498}
1499
1500template<int MD1, int MQ1>
1501MFEM_HOST_DEVICE inline void EvalZ(const int D1D, const int Q1D,
1502 const real_t (&sB)[MQ1*MD1],
1503 const real_t (&sDQQ)[MD1*MQ1*MQ1],
1504 real_t (&sQQQ)[MQ1*MQ1*MQ1])
1505{
1506 ConstDeviceMatrix B(sB, D1D, Q1D);
1507 const DeviceCube DQQ(sDQQ, Q1D, Q1D, D1D);
1508 DeviceCube QQQ(sQQQ, Q1D, Q1D, Q1D);
1509 EvalZ(D1D,Q1D,B,DQQ,QQQ);
1510}
1511
1512/// Pull 3D Scalar Evaluation
1513MFEM_HOST_DEVICE inline void PullEval(const int x, const int y, const int z,
1514 const DeviceCube &QQQ,
1515 real_t &X)
1516{
1517 X = QQQ(z,y,x);
1518}
1519
1520template<int MQ1>
1521MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
1522 const int x, const int y, const int z,
1523 const real_t (&sQQQ)[MQ1*MQ1*MQ1],
1524 real_t &X)
1525{
1526 const DeviceCube QQQ(sQQQ, Q1D, Q1D, Q1D);
1527 PullEval(x,y,z,QQQ,X);
1528}
1529
1530/// Load 3D input vector into shared memory
1531template<int MD1>
1532MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
1533 const DeviceTensor<5, const real_t> &X,
1534 real_t (*sm)[MD1*MD1*MD1])
1535{
1536 DeviceCube Xx(sm[0], D1D, D1D, D1D);
1537 DeviceCube Xy(sm[1], D1D, D1D, D1D);
1538 DeviceCube Xz(sm[2], D1D, D1D, D1D);
1539
1540 MFEM_FOREACH_THREAD(dz,z,D1D)
1541 {
1542 MFEM_FOREACH_THREAD(dy,y,D1D)
1543 {
1544 MFEM_FOREACH_THREAD(dx,x,D1D)
1545 {
1546 Xx(dx,dy,dz) = X(dx,dy,dz,0,e);
1547 Xy(dx,dy,dz) = X(dx,dy,dz,1,e);
1548 Xz(dx,dy,dz) = X(dx,dy,dz,2,e);
1549 }
1550 }
1551 }
1552 MFEM_SYNC_THREAD;
1553}
1554
1555/// 3D Vector Evaluation, 1/3 (only B)
1556template<int MD1, int MQ1>
1557MFEM_HOST_DEVICE inline void EvalX(const int D1D, const int Q1D,
1558 const real_t (&sB)[MQ1*MD1],
1559 const real_t (&sDDD)[3][MD1*MD1*MD1],
1560 real_t (&sDDQ)[3][MD1*MD1*MQ1])
1561{
1562 ConstDeviceMatrix B(sB, D1D, Q1D);
1563 ConstDeviceCube Xx(sDDD[0], D1D, D1D, D1D);
1564 ConstDeviceCube Xy(sDDD[1], D1D, D1D, D1D);
1565 ConstDeviceCube Xz(sDDD[2], D1D, D1D, D1D);
1566 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1567 DeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
1568 DeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
1569
1570 MFEM_FOREACH_THREAD(dz,z,D1D)
1571 {
1572 MFEM_FOREACH_THREAD(dy,y,D1D)
1573 {
1574 MFEM_FOREACH_THREAD(qx,x,Q1D)
1575 {
1576 real_t u[3] = {0.0, 0.0, 0.0};
1577 for (int dx = 0; dx < D1D; ++dx)
1578 {
1579 const real_t Bx = B(dx,qx);
1580 u[0] += Bx * Xx(dx,dy,dz);
1581 u[1] += Bx * Xy(dx,dy,dz);
1582 u[2] += Bx * Xz(dx,dy,dz);
1583 }
1584 XxB(qx,dy,dz) = u[0];
1585 XyB(qx,dy,dz) = u[1];
1586 XzB(qx,dy,dz) = u[2];
1587 }
1588 }
1589 }
1590 MFEM_SYNC_THREAD;
1591}
1592
1593/// 3D Vector Evaluation, 2/3 (only B)
1594template<int MD1, int MQ1>
1595MFEM_HOST_DEVICE inline void EvalY(const int D1D, const int Q1D,
1596 const real_t (&sB)[MQ1*MD1],
1597 const real_t (&sDDQ)[3][MD1*MD1*MQ1],
1598 real_t (&sDQQ)[3][MD1*MQ1*MQ1])
1599{
1600 ConstDeviceMatrix B(sB, D1D, Q1D);
1601 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1602 ConstDeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
1603 ConstDeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
1604 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1605 DeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
1606 DeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
1607
1608 MFEM_FOREACH_THREAD(dz,z,D1D)
1609 {
1610 MFEM_FOREACH_THREAD(qy,y,Q1D)
1611 {
1612 MFEM_FOREACH_THREAD(qx,x,Q1D)
1613 {
1614 real_t u[3] = {0.0, 0.0, 0.0};
1615 for (int dy = 0; dy < D1D; ++dy)
1616 {
1617 const real_t By = B(dy,qy);
1618 u[0] += XxB(qx,dy,dz) * By;
1619 u[1] += XyB(qx,dy,dz) * By;
1620 u[2] += XzB(qx,dy,dz) * By;
1621 }
1622 XxBB(qx,qy,dz) = u[0];
1623 XyBB(qx,qy,dz) = u[1];
1624 XzBB(qx,qy,dz) = u[2];
1625 }
1626 }
1627 }
1628 MFEM_SYNC_THREAD;
1629}
1630
1631/// 3D Vector Evaluation, 3/3 (only B)
1632template<int MD1, int MQ1>
1633MFEM_HOST_DEVICE inline void EvalZ(const int D1D, const int Q1D,
1634 const real_t (&sB)[MQ1*MD1],
1635 const real_t (&sDQQ)[3][MD1*MQ1*MQ1],
1636 real_t (&sQQQ)[3][MQ1*MQ1*MQ1])
1637{
1638 ConstDeviceMatrix B(sB, D1D, Q1D);
1639 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1640 ConstDeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
1641 ConstDeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
1642 DeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
1643 DeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
1644 DeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
1645
1646 MFEM_FOREACH_THREAD(qz,z,Q1D)
1647 {
1648 MFEM_FOREACH_THREAD(qy,y,Q1D)
1649 {
1650 MFEM_FOREACH_THREAD(qx,x,Q1D)
1651 {
1652 real_t u[3] = {0.0, 0.0, 0.0};
1653 for (int dz = 0; dz < D1D; ++dz)
1654 {
1655 const real_t Bz = B(dz,qz);
1656 u[0] += XxBB(qx,qy,dz) * Bz;
1657 u[1] += XyBB(qx,qy,dz) * Bz;
1658 u[2] += XzBB(qx,qy,dz) * Bz;
1659 }
1660 XxBBB(qx,qy,qz) = u[0];
1661 XyBBB(qx,qy,qz) = u[1];
1662 XzBBB(qx,qy,qz) = u[2];
1663 }
1664 }
1665 }
1666 MFEM_SYNC_THREAD;
1667}
1668
1669/// Pull 3D Vector Evaluation
1670template<int MQ1>
1671MFEM_HOST_DEVICE inline void PullEval(const int Q1D,
1672 const int x, const int y, const int z,
1673 const real_t (&sQQQ)[3][MQ1*MQ1*MQ1],
1674 real_t (&X)[3])
1675{
1676 ConstDeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
1677 ConstDeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
1678 ConstDeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
1679
1680 X[0] = XxBBB(x,y,z);
1681 X[1] = XyBBB(x,y,z);
1682 X[2] = XzBBB(x,y,z);
1683}
1684
1685/// Push 3D Vector Evaluation
1686template<int MQ1>
1687MFEM_HOST_DEVICE inline void PushEval(const int Q1D,
1688 const int x, const int y, const int z,
1689 const real_t (&A)[3],
1690 real_t (&sQQQ)[3][MQ1*MQ1*MQ1])
1691{
1692 DeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
1693 DeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
1694 DeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
1695
1696 XxBBB(x,y,z) = A[0];
1697 XyBBB(x,y,z) = A[1];
1698 XzBBB(x,y,z) = A[2];
1699}
1700
1701/// 3D Transposed Vector Evaluation, 1/3
1702template<int MD1, int MQ1>
1703MFEM_HOST_DEVICE inline void EvalXt(const int D1D, const int Q1D,
1704 const real_t (&sB)[MQ1*MD1],
1705 const real_t (&sQQQ)[3][MQ1*MQ1*MQ1],
1706 real_t (&sDQQ)[3][MD1*MQ1*MQ1])
1707{
1708 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1709 ConstDeviceCube XxBBB(sQQQ[0], Q1D, Q1D, Q1D);
1710 ConstDeviceCube XyBBB(sQQQ[1], Q1D, Q1D, Q1D);
1711 ConstDeviceCube XzBBB(sQQQ[2], Q1D, Q1D, Q1D);
1712 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1713 DeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
1714 DeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
1715
1716 MFEM_FOREACH_THREAD(qz,z,Q1D)
1717 {
1718 MFEM_FOREACH_THREAD(qy,y,Q1D)
1719 {
1720 MFEM_FOREACH_THREAD(dx,x,D1D)
1721 {
1722 real_t u[3] = {0.0, 0.0, 0.0};
1723 for (int qx = 0; qx < Q1D; ++qx)
1724 {
1725 const real_t Btx = Bt(qx,dx);
1726 u[0] += XxBBB(qx,qy,qz) * Btx;
1727 u[1] += XyBBB(qx,qy,qz) * Btx;
1728 u[2] += XzBBB(qx,qy,qz) * Btx;
1729 }
1730 XxBB(qz,qy,dx) = u[0];
1731 XyBB(qz,qy,dx) = u[1];
1732 XzBB(qz,qy,dx) = u[2];
1733 }
1734 }
1735 }
1736 MFEM_SYNC_THREAD;
1737}
1738
1739/// 3D Transposed Vector Evaluation, 2/3
1740template<int MD1, int MQ1>
1741MFEM_HOST_DEVICE inline void EvalYt(const int D1D, const int Q1D,
1742 const real_t (&sB)[MQ1*MD1],
1743 const real_t (&sDQQ)[3][MD1*MQ1*MQ1],
1744 real_t (&sDDQ)[3][MD1*MD1*MQ1])
1745{
1746 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1747 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1748 ConstDeviceCube XyBB(sDQQ[1], Q1D, Q1D, D1D);
1749 ConstDeviceCube XzBB(sDQQ[2], Q1D, Q1D, D1D);
1750 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1751 DeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
1752 DeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
1753
1754 MFEM_FOREACH_THREAD(qz,z,Q1D)
1755 {
1756 MFEM_FOREACH_THREAD(dy,y,D1D)
1757 {
1758 MFEM_FOREACH_THREAD(dx,x,D1D)
1759 {
1760 real_t u[3] = {0.0, 0.0, 0.0};
1761 for (int qy = 0; qy < Q1D; ++qy)
1762 {
1763 const real_t Bty = Bt(qy,dy);
1764 u[0] += XxBB(qz,qy,dx) * Bty;
1765 u[1] += XyBB(qz,qy,dx) * Bty;
1766 u[2] += XzBB(qz,qy,dx) * Bty;
1767
1768 }
1769 XxB(qz,dy,dx) = u[0];
1770 XyB(qz,dy,dx) = u[1];
1771 XzB(qz,dy,dx)= u[2];
1772 }
1773 }
1774 }
1775 MFEM_SYNC_THREAD;
1776}
1777
1778/// 3D Transposed Vector Evaluation, 3/3
1779template<int MD1, int MQ1>
1780MFEM_HOST_DEVICE inline void EvalZt(const int D1D, const int Q1D,
1781 const real_t (&sB)[MQ1*MD1],
1782 const real_t (&sDDQ)[3][MD1*MD1*MQ1],
1783 const DeviceTensor<5> &Y, // output
1784 const int e)
1785{
1786 ConstDeviceMatrix Bt(sB, Q1D, D1D);
1787 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1788 ConstDeviceCube XyB(sDDQ[1], Q1D, D1D, D1D);
1789 ConstDeviceCube XzB(sDDQ[2], Q1D, D1D, D1D);
1790
1791 MFEM_FOREACH_THREAD(dz,z,D1D)
1792 {
1793 MFEM_FOREACH_THREAD(dy,y,D1D)
1794 {
1795 MFEM_FOREACH_THREAD(dx,x,D1D)
1796 {
1797 real_t u[3] = {0.0, 0.0, 0.0};
1798 for (int qz = 0; qz < Q1D; ++qz)
1799 {
1800 const real_t Btz = Bt(qz,dz);
1801 u[0] += XxB(qz,dy,dx) * Btz;
1802 u[1] += XyB(qz,dy,dx) * Btz;
1803 u[2] += XzB(qz,dy,dx) * Btz;
1804 }
1805 Y(dx,dy,dz,0,e) += u[0];
1806 Y(dx,dy,dz,1,e) += u[1];
1807 Y(dx,dy,dz,2,e) += u[2];
1808 }
1809 }
1810 }
1811}
1812
1813/// 3D Gradient, 1/3
1814template<int MD1, int MQ1>
1815MFEM_HOST_DEVICE inline void GradX(const int D1D, const int Q1D,
1816 const real_t (*sBG)[MQ1*MD1],
1817 const real_t (*sDDD)[MD1*MD1*MD1],
1818 real_t (*sDDQ)[MD1*MD1*MQ1])
1819{
1820 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1821 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1822 ConstDeviceCube Xx(sDDD[0], D1D, D1D, D1D);
1823 ConstDeviceCube Xy(sDDD[1], D1D, D1D, D1D);
1824 ConstDeviceCube Xz(sDDD[2], D1D, D1D, D1D);
1825 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1826 DeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
1827 DeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
1828 DeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
1829 DeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
1830 DeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
1831
1832 MFEM_FOREACH_THREAD(dz,z,D1D)
1833 {
1834 MFEM_FOREACH_THREAD(dy,y,D1D)
1835 {
1836 MFEM_FOREACH_THREAD(qx,x,Q1D)
1837 {
1838 real_t u[3] = {0.0, 0.0, 0.0};
1839 real_t v[3] = {0.0, 0.0, 0.0};
1840 for (int dx = 0; dx < D1D; ++dx)
1841 {
1842 const real_t xx = Xx(dx,dy,dz);
1843 const real_t xy = Xy(dx,dy,dz);
1844 const real_t xz = Xz(dx,dy,dz);
1845 const real_t Bx = B(dx,qx);
1846 const real_t Gx = G(dx,qx);
1847 u[0] += Bx * xx;
1848 u[1] += Bx * xy;
1849 u[2] += Bx * xz;
1850
1851 v[0] += Gx * xx;
1852 v[1] += Gx * xy;
1853 v[2] += Gx * xz;
1854 }
1855 XxB(qx,dy,dz) = u[0];
1856 XyB(qx,dy,dz) = u[1];
1857 XzB(qx,dy,dz) = u[2];
1858
1859 XxG(qx,dy,dz) = v[0];
1860 XyG(qx,dy,dz) = v[1];
1861 XzG(qx,dy,dz) = v[2];
1862 }
1863 }
1864 }
1865 MFEM_SYNC_THREAD;
1866}
1867
1868/// 3D Gradient, 2/3
1869template<int MD1, int MQ1>
1870MFEM_HOST_DEVICE inline void GradY(const int D1D, const int Q1D,
1871 const real_t (*sBG)[MQ1*MD1],
1872 const real_t (*sDDQ)[MD1*MD1*MQ1],
1873 real_t (*sDQQ)[MD1*MQ1*MQ1])
1874{
1875 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1876 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1877 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
1878 ConstDeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
1879 ConstDeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
1880 ConstDeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
1881 ConstDeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
1882 ConstDeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
1883 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1884 DeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
1885 DeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
1886 DeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
1887 DeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
1888 DeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
1889 DeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
1890 DeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
1891 DeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
1892
1893 MFEM_FOREACH_THREAD(dz,z,D1D)
1894 {
1895 MFEM_FOREACH_THREAD(qy,y,Q1D)
1896 {
1897 MFEM_FOREACH_THREAD(qx,x,Q1D)
1898 {
1899 real_t u[3] = {0.0, 0.0, 0.0};
1900 real_t v[3] = {0.0, 0.0, 0.0};
1901 real_t w[3] = {0.0, 0.0, 0.0};
1902 for (int dy = 0; dy < D1D; ++dy)
1903 {
1904 const real_t By = B(dy,qy);
1905 const real_t Gy = G(dy,qy);
1906
1907 u[0] += XxB(qx,dy,dz) * By;
1908 u[1] += XyB(qx,dy,dz) * By;
1909 u[2] += XzB(qx,dy,dz) * By;
1910
1911 v[0] += XxG(qx,dy,dz) * By;
1912 v[1] += XyG(qx,dy,dz) * By;
1913 v[2] += XzG(qx,dy,dz) * By;
1914
1915 w[0] += XxB(qx,dy,dz) * Gy;
1916 w[1] += XyB(qx,dy,dz) * Gy;
1917 w[2] += XzB(qx,dy,dz) * Gy;
1918 }
1919 XxBB(qx,qy,dz) = u[0];
1920 XyBB(qx,qy,dz) = u[1];
1921 XzBB(qx,qy,dz) = u[2];
1922
1923 XxBG(qx,qy,dz) = v[0];
1924 XyBG(qx,qy,dz) = v[1];
1925 XzBG(qx,qy,dz) = v[2];
1926
1927 XxGB(qx,qy,dz) = w[0];
1928 XyGB(qx,qy,dz) = w[1];
1929 XzGB(qx,qy,dz) = w[2];
1930 }
1931 }
1932 }
1933 MFEM_SYNC_THREAD;
1934}
1935
1936/// 3D Gradient, 3/3
1937template<int MD1, int MQ1>
1938MFEM_HOST_DEVICE inline void GradZ(const int D1D, const int Q1D,
1939 const real_t (*sBG)[MQ1*MD1],
1940 const real_t (*sDQQ)[MD1*MQ1*MQ1],
1941 real_t (*sQQQ)[MQ1*MQ1*MQ1])
1942{
1943 ConstDeviceMatrix B(sBG[0], D1D, Q1D);
1944 ConstDeviceMatrix G(sBG[1], D1D, Q1D);
1945 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
1946 ConstDeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
1947 ConstDeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
1948 ConstDeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
1949 ConstDeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
1950 ConstDeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
1951 ConstDeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
1952 ConstDeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
1953 ConstDeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
1954 DeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
1955 DeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
1956 DeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
1957 DeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
1958 DeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
1959 DeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
1960 DeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
1961 DeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
1962 DeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
1963
1964 MFEM_FOREACH_THREAD(qz,z,Q1D)
1965 {
1966 MFEM_FOREACH_THREAD(qy,y,Q1D)
1967 {
1968 MFEM_FOREACH_THREAD(qx,x,Q1D)
1969 {
1970 real_t u[3] = {0.0, 0.0, 0.0};
1971 real_t v[3] = {0.0, 0.0, 0.0};
1972 real_t w[3] = {0.0, 0.0, 0.0};
1973 for (int dz = 0; dz < D1D; ++dz)
1974 {
1975 const real_t Bz = B(dz,qz);
1976 const real_t Gz = G(dz,qz);
1977
1978 u[0] += XxBG(qx,qy,dz) * Bz;
1979 u[1] += XyBG(qx,qy,dz) * Bz;
1980 u[2] += XzBG(qx,qy,dz) * Bz;
1981
1982 v[0] += XxGB(qx,qy,dz) * Bz;
1983 v[1] += XyGB(qx,qy,dz) * Bz;
1984 v[2] += XzGB(qx,qy,dz) * Bz;
1985
1986 w[0] += XxBB(qx,qy,dz) * Gz;
1987 w[1] += XyBB(qx,qy,dz) * Gz;
1988 w[2] += XzBB(qx,qy,dz) * Gz;
1989 }
1990 XxBBG(qx,qy,qz) = u[0];
1991 XyBBG(qx,qy,qz) = u[1];
1992 XzBBG(qx,qy,qz) = u[2];
1993
1994 XxBGB(qx,qy,qz) = v[0];
1995 XyBGB(qx,qy,qz) = v[1];
1996 XzBGB(qx,qy,qz) = v[2];
1997
1998 XxGBB(qx,qy,qz)= w[0];
1999 XyGBB(qx,qy,qz) = w[1];
2000 XzGBB(qx,qy,qz) = w[2];
2001 }
2002 }
2003 }
2004 MFEM_SYNC_THREAD;
2005}
2006
2007/// Pull 3D Gradient
2008template<int MQ1>
2009MFEM_HOST_DEVICE inline void PullGrad(const int Q1D,
2010 const int x, const int y, const int z,
2011 const real_t (*sQQQ)[MQ1*MQ1*MQ1],
2012 real_t *Jpr)
2013{
2014 ConstDeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
2015 ConstDeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
2016 ConstDeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
2017 ConstDeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
2018 ConstDeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
2019 ConstDeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
2020 ConstDeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
2021 ConstDeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
2022 ConstDeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
2023
2024 Jpr[0] = XxBBG(x,y,z);
2025 Jpr[3] = XxBGB(x,y,z);
2026 Jpr[6] = XxGBB(x,y,z);
2027 Jpr[1] = XyBBG(x,y,z);
2028 Jpr[4] = XyBGB(x,y,z);
2029 Jpr[7] = XyGBB(x,y,z);
2030 Jpr[2] = XzBBG(x,y,z);
2031 Jpr[5] = XzBGB(x,y,z);
2032 Jpr[8] = XzGBB(x,y,z);
2033}
2034
2035/// Push 3D Gradient
2036template<int MQ1>
2037MFEM_HOST_DEVICE inline void PushGrad(const int Q1D,
2038 const int x, const int y, const int z,
2039 const real_t *A,
2040 real_t (&sQQQ)[9][MQ1*MQ1*MQ1])
2041{
2042 DeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
2043 DeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
2044 DeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
2045 DeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
2046 DeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
2047 DeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
2048 DeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
2049 DeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
2050 DeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
2051
2052 XxBBG(x,y,z) = A[0];
2053 XxBGB(x,y,z) = A[1];
2054 XxGBB(x,y,z) = A[2];
2055 XyBBG(x,y,z) = A[3];
2056 XyBGB(x,y,z) = A[4];
2057 XyGBB(x,y,z) = A[5];
2058 XzBBG(x,y,z) = A[6];
2059 XzBGB(x,y,z) = A[7];
2060 XzGBB(x,y,z) = A[8];
2061}
2062
2063/// 3D Transposed Gradient, 1/3
2064template<int MD1, int MQ1>
2065MFEM_HOST_DEVICE inline void GradZt(const int D1D, const int Q1D,
2066 const real_t (&sBG)[2][MQ1*MD1],
2067 const real_t (&sQQQ)[9][MQ1*MQ1*MQ1],
2068 real_t (&sDQQ)[9][MD1*MQ1*MQ1])
2069{
2070
2071 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
2072 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
2073 ConstDeviceCube XxBBG(sQQQ[0], Q1D, Q1D, Q1D);
2074 ConstDeviceCube XxBGB(sQQQ[1], Q1D, Q1D, Q1D);
2075 ConstDeviceCube XxGBB(sQQQ[2], Q1D, Q1D, Q1D);
2076 ConstDeviceCube XyBBG(sQQQ[3], Q1D, Q1D, Q1D);
2077 ConstDeviceCube XyBGB(sQQQ[4], Q1D, Q1D, Q1D);
2078 ConstDeviceCube XyGBB(sQQQ[5], Q1D, Q1D, Q1D);
2079 ConstDeviceCube XzBBG(sQQQ[6], Q1D, Q1D, Q1D);
2080 ConstDeviceCube XzBGB(sQQQ[7], Q1D, Q1D, Q1D);
2081 ConstDeviceCube XzGBB(sQQQ[8], Q1D, Q1D, Q1D);
2082 DeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
2083 DeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
2084 DeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
2085 DeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
2086 DeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
2087 DeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
2088 DeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
2089 DeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
2090 DeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
2091
2092 MFEM_FOREACH_THREAD(qz,z,Q1D)
2093 {
2094 MFEM_FOREACH_THREAD(qy,y,Q1D)
2095 {
2096 MFEM_FOREACH_THREAD(dx,x,D1D)
2097 {
2098 real_t u[3] = {0.0, 0.0, 0.0};
2099 real_t v[3] = {0.0, 0.0, 0.0};
2100 real_t w[3] = {0.0, 0.0, 0.0};
2101 for (int qx = 0; qx < Q1D; ++qx)
2102 {
2103 const real_t Btx = Bt(qx,dx);
2104 const real_t Gtx = Gt(qx,dx);
2105
2106 u[0] += XxBBG(qx,qy,qz) * Gtx;
2107 v[0] += XxBGB(qx,qy,qz) * Btx;
2108 w[0] += XxGBB(qx,qy,qz) * Btx;
2109
2110 u[1] += XyBBG(qx,qy,qz) * Gtx;
2111 v[1] += XyBGB(qx,qy,qz) * Btx;
2112 w[1] += XyGBB(qx,qy,qz) * Btx;
2113
2114 u[2] += XzBBG(qx,qy,qz) * Gtx;
2115 v[2] += XzBGB(qx,qy,qz) * Btx;
2116 w[2] += XzGBB(qx,qy,qz) * Btx;
2117 }
2118 XxBB(qz,qy,dx) = u[0];
2119 XxBG(qz,qy,dx) = v[0];
2120 XxGB(qz,qy,dx) = w[0];
2121
2122 XyBB(qz,qy,dx) = u[1];
2123 XyBG(qz,qy,dx) = v[1];
2124 XyGB(qz,qy,dx) = w[1];
2125
2126 XzBB(qz,qy,dx) = u[2];
2127 XzBG(qz,qy,dx) = v[2];
2128 XzGB(qz,qy,dx) = w[2];
2129 }
2130 }
2131 }
2132 MFEM_SYNC_THREAD;
2133}
2134
2135/// 3D Transposed Gradient, 2/3
2136template<int MD1, int MQ1>
2137MFEM_HOST_DEVICE inline void GradYt(const int D1D, const int Q1D,
2138 const real_t (&sBG)[2][MQ1*MD1],
2139 const real_t (&sDQQ)[9][MD1*MQ1*MQ1],
2140 real_t (&sDDQ)[9][MD1*MD1*MQ1])
2141{
2142 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
2143 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
2144 ConstDeviceCube XxBB(sDQQ[0], Q1D, Q1D, D1D);
2145 ConstDeviceCube XxBG(sDQQ[1], Q1D, Q1D, D1D);
2146 ConstDeviceCube XxGB(sDQQ[2], Q1D, Q1D, D1D);
2147 ConstDeviceCube XyBB(sDQQ[3], Q1D, Q1D, D1D);
2148 ConstDeviceCube XyBG(sDQQ[4], Q1D, Q1D, D1D);
2149 ConstDeviceCube XyGB(sDQQ[5], Q1D, Q1D, D1D);
2150 ConstDeviceCube XzBB(sDQQ[6], Q1D, Q1D, D1D);
2151 ConstDeviceCube XzBG(sDQQ[7], Q1D, Q1D, D1D);
2152 ConstDeviceCube XzGB(sDQQ[8], Q1D, Q1D, D1D);
2153 DeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
2154 DeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
2155 DeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
2156 DeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
2157 DeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
2158 DeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
2159 DeviceCube XxC(sDDQ[6], Q1D, D1D, D1D);
2160 DeviceCube XyC(sDDQ[7], Q1D, D1D, D1D);
2161 DeviceCube XzC(sDDQ[8], Q1D, D1D, D1D);
2162
2163 MFEM_FOREACH_THREAD(qz,z,Q1D)
2164 {
2165 MFEM_FOREACH_THREAD(dy,y,D1D)
2166 {
2167 MFEM_FOREACH_THREAD(dx,x,D1D)
2168 {
2169 real_t u[3] = {0.0, 0.0, 0.0};
2170 real_t v[3] = {0.0, 0.0, 0.0};
2171 real_t w[3] = {0.0, 0.0, 0.0};
2172 for (int qy = 0; qy < Q1D; ++qy)
2173 {
2174 const real_t Bty = Bt(qy,dy);
2175 const real_t Gty = Gt(qy,dy);
2176
2177 u[0] += XxBB(qz,qy,dx) * Bty;
2178 v[0] += XxBG(qz,qy,dx) * Gty;
2179 w[0] += XxGB(qz,qy,dx) * Bty;
2180
2181 u[1] += XyBB(qz,qy,dx) * Bty;
2182 v[1] += XyBG(qz,qy,dx) * Gty;
2183 w[1] += XyGB(qz,qy,dx) * Bty;
2184
2185 u[2] += XzBB(qz,qy,dx) * Bty;
2186 v[2] += XzBG(qz,qy,dx) * Gty;
2187 w[2] += XzGB(qz,qy,dx) * Bty;
2188
2189 }
2190 XxB(qz,dy,dx) = u[0];
2191 XxC(qz,dy,dx) = v[0];
2192 XxG(qz,dy,dx) = w[0];
2193
2194 XyB(qz,dy,dx) = u[1];
2195 XyC(qz,dy,dx) = v[1];
2196 XyG(qz,dy,dx) = w[1];
2197
2198 XzB(qz,dy,dx) = u[2];
2199 XzC(qz,dy,dx) = v[2];
2200 XzG(qz,dy,dx) = w[2];
2201 }
2202 }
2203 }
2204 MFEM_SYNC_THREAD;
2205}
2206
2207/// 3D Transposed Gradient, 3/3
2208template<int MD1, int MQ1>
2209MFEM_HOST_DEVICE inline void GradXt(const int D1D, const int Q1D,
2210 const real_t (&sBG)[2][MQ1*MD1],
2211 const real_t (&sDDQ)[9][MD1*MD1*MQ1],
2212 const DeviceTensor<5> &Y, // output
2213 const int e)
2214{
2215 ConstDeviceMatrix Bt(sBG[0], Q1D, D1D);
2216 ConstDeviceMatrix Gt(sBG[1], Q1D, D1D);
2217 ConstDeviceCube XxB(sDDQ[0], Q1D, D1D, D1D);
2218 ConstDeviceCube XxG(sDDQ[1], Q1D, D1D, D1D);
2219 ConstDeviceCube XyB(sDDQ[2], Q1D, D1D, D1D);
2220 ConstDeviceCube XyG(sDDQ[3], Q1D, D1D, D1D);
2221 ConstDeviceCube XzB(sDDQ[4], Q1D, D1D, D1D);
2222 ConstDeviceCube XzG(sDDQ[5], Q1D, D1D, D1D);
2223 ConstDeviceCube XxC(sDDQ[6], Q1D, D1D, D1D);
2224 ConstDeviceCube XyC(sDDQ[7], Q1D, D1D, D1D);
2225 ConstDeviceCube XzC(sDDQ[8], Q1D, D1D, D1D);
2226
2227 MFEM_FOREACH_THREAD(dz,z,D1D)
2228 {
2229 MFEM_FOREACH_THREAD(dy,y,D1D)
2230 {
2231 MFEM_FOREACH_THREAD(dx,x,D1D)
2232 {
2233 real_t u[3] = {0.0, 0.0, 0.0};
2234 real_t v[3] = {0.0, 0.0, 0.0};
2235 real_t w[3] = {0.0, 0.0, 0.0};
2236 for (int qz = 0; qz < Q1D; ++qz)
2237 {
2238 const real_t Btz = Bt(qz,dz);
2239 const real_t Gtz = Gt(qz,dz);
2240
2241 u[0] += XxB(qz,dy,dx) * Btz;
2242 v[0] += XxC(qz,dy,dx) * Btz;
2243 w[0] += XxG(qz,dy,dx) * Gtz;
2244
2245 u[1] += XyB(qz,dy,dx) * Btz;
2246 v[1] += XyC(qz,dy,dx)* Btz;
2247 w[1] += XyG(qz,dy,dx) * Gtz;
2248
2249 u[2] += XzB(qz,dy,dx) * Btz;
2250 v[2] += XzC(qz,dy,dx) * Btz;
2251 w[2] += XzG(qz,dy,dx) * Gtz;
2252 }
2253 Y(dx,dy,dz,0,e) += u[0] + v[0] + w[0];
2254 Y(dx,dy,dz,1,e) += u[1] + v[1] + w[1];
2255 Y(dx,dy,dz,2,e) += u[2] + v[2] + w[2];
2256 }
2257 }
2258 }
2259}
2260
2261} // namespace internal
2262
2263} // namespace kernels
2264
2265} // namespace mfem
2266
2267#endif // MFEM_FEM_KERNELS_HPP
b
real_t b
Definition lissajous.cpp:42
DIM
constexpr int DIM
Definition minimal-surface.cpp:71
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:17
mfem::DeviceCube
DeviceTensor< 3, real_t > DeviceCube
Definition dtensor.hpp:153
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::ConstDeviceCube
DeviceTensor< 3, const real_t > ConstDeviceCube
Definition dtensor.hpp:154
mfem::Transpose
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:443
mfem::real_t
float real_t
Definition config.hpp:46
mfem::ConstDeviceMatrix
DeviceTensor< 2, const real_t > ConstDeviceMatrix
Definition dtensor.hpp:151
mfem::DeviceMatrix
DeviceTensor< 2, real_t > DeviceMatrix
Definition dtensor.hpp:150
tensor.hpp
Implementation of the tensor class.
Generated on Thu Dec 11 2025 11:52:55 for MFEM by doxygen 1.11.0