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