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