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