MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
grad.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// Internal header, included only by .cpp files.
13// Template function implementations.
14
15#ifndef MFEM_QUADINTERP_GRAD
16#define MFEM_QUADINTERP_GRAD
17
18#include "../quadinterpolator.hpp"
19#include "../../general/forall.hpp"
20#include "../../linalg/dtensor.hpp"
21#include "../../fem/kernels.hpp"
22#include "../../linalg/kernels.hpp"
23
24namespace mfem
25{
26
27namespace internal
28{
29
30namespace quadrature_interpolator
31{
32
33template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS>
34static void Derivatives1D(const int NE,
35 const real_t *b_,
36 const real_t *g_,
37 const real_t *j_,
38 const real_t *x_,
39 real_t *y_,
40 const int sdim,
41 const int vdim,
42 const int d1d,
43 const int q1d)
44{
45 MFEM_CONTRACT_VAR(b_);
46 const int SDIM = GRAD_PHYS ? sdim : 1;
47 const auto g = Reshape(g_, q1d, d1d);
48 const auto j = Reshape(j_, q1d, SDIM, NE);
49 const auto x = Reshape(x_, d1d, vdim, NE);
50 auto y = Q_LAYOUT == QVectorLayout::byNODES ?
51 Reshape(y_, q1d, vdim, SDIM, NE):
52 Reshape(y_, vdim, SDIM, q1d, NE);
53
54 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
55 {
56 for (int c = 0; c < vdim; c++)
57 {
58 for (int q = 0; q < q1d; q++)
59 {
60 real_t du[3] = {0.0, 0.0, 0.0};
61 for (int d = 0; d < d1d; d++)
62 {
63 du[0] += g(q, d) * x(d, c, e);
64 }
65 if (GRAD_PHYS)
66 {
67 if (SDIM == 1) { du[0] /= j(q, 0, e); }
68 else if (SDIM == 2)
69 {
70 const real_t Jloc[2] = {j(q,0,e), j(q,1,e)};
71 real_t Jinv[3];
72 kernels::CalcLeftInverse<2,1>(Jloc, Jinv);
73 const real_t U = Jinv[0]*du[0];
74 const real_t V = Jinv[1]*du[0];
75 du[0] = U;
76 du[1] = V;
77 }
78 else // SDIM == 3
79 {
80 const real_t Jloc[3] = {j(q,0,e), j(q,1,e), j(q,2,e)};
81 real_t Jinv[3];
82 kernels::CalcLeftInverse<3,1>(Jloc, Jinv);
83 const real_t U = Jinv[0]*du[0];
84 const real_t V = Jinv[1]*du[0];
85 const real_t W = Jinv[2]*du[0];
86 du[0] = U;
87 du[1] = V;
88 du[2] = W;
89 }
90 }
91 for (int d = 0; d < SDIM; ++d)
92 {
93 if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c, d, q, e) = du[d]; }
94 if (Q_LAYOUT == QVectorLayout::byNODES) { y(q, c, d, e) = du[d]; }
95 }
96 }
97 }
98 });
99}
100
101// Template compute kernel for derivatives in 2D: tensor product version.
102template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
103 int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
104 int T_NBZ = 1>
105static void Derivatives2D(const int NE,
106 const real_t *b_,
107 const real_t *g_,
108 const real_t *j_,
109 const real_t *x_,
110 real_t *y_,
111 const int sdim = 2,
112 const int vdim = 0,
113 const int d1d = 0,
114 const int q1d = 0)
115{
116 const int D1D = T_D1D ? T_D1D : d1d;
117 const int Q1D = T_Q1D ? T_Q1D : q1d;
118 const int VDIM = T_VDIM ? T_VDIM : vdim;
119 const int SDIM = GRAD_PHYS ? sdim : 2;
120 static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
121
122 const auto b = Reshape(b_, Q1D, D1D);
123 const auto g = Reshape(g_, Q1D, D1D);
124 const auto j = Reshape(j_, Q1D, Q1D, SDIM, 2, NE);
125 const auto x = Reshape(x_, D1D, D1D, VDIM, NE);
126 auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
127 Reshape(y_, Q1D, Q1D, VDIM, SDIM, NE):
128 Reshape(y_, VDIM, SDIM, Q1D, Q1D, NE);
129
130 mfem::forall_2D_batch(NE, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE (int e)
131 {
132 const int D1D = T_D1D ? T_D1D : d1d;
133 const int Q1D = T_Q1D ? T_Q1D : q1d;
134 const int VDIM = T_VDIM ? T_VDIM : vdim;
135 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
136 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
137
138 const int tidz = MFEM_THREAD_ID(z);
139 MFEM_SHARED real_t BG[2][MQ1*MD1];
140 kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,BG);
141 DeviceMatrix B(BG[0], D1D, Q1D);
142 DeviceMatrix G(BG[1], D1D, Q1D);
143
144 MFEM_SHARED real_t XY[NBZ][MD1*MD1];
145 DeviceTensor<2> X((real_t*)(XY+tidz), D1D, D1D);
146
147 MFEM_SHARED real_t s_DQ[2][NBZ][MD1*MQ1];
148 DeviceTensor<2> DQ0(s_DQ[0][tidz], D1D, Q1D);
149 DeviceTensor<2> DQ1(s_DQ[1][tidz], D1D, Q1D);
150
151 for (int c = 0; c < VDIM; ++c)
152 {
153 kernels::internal::LoadX<MD1,NBZ>(e,D1D,c,x,XY);
154 MFEM_FOREACH_THREAD(dy,y,D1D)
155 {
156 MFEM_FOREACH_THREAD(qx,x,Q1D)
157 {
158 real_t u = 0.0;
159 real_t v = 0.0;
160 for (int dx = 0; dx < D1D; ++dx)
161 {
162 const real_t input = X(dx,dy);
163 u += input * B(dx,qx);
164 v += input * G(dx,qx);
165 }
166 DQ0(dy,qx) = u;
167 DQ1(dy,qx) = v;
168 }
169 }
170 MFEM_SYNC_THREAD;
171 MFEM_FOREACH_THREAD(qy,y,Q1D)
172 {
173 MFEM_FOREACH_THREAD(qx,x,Q1D)
174 {
175 real_t du[3] = {0.0, 0.0, 0.0};
176 for (int dy = 0; dy < D1D; ++dy)
177 {
178 du[0] += DQ1(dy,qx) * B(dy,qy);
179 du[1] += DQ0(dy,qx) * G(dy,qy);
180 }
181 if (GRAD_PHYS)
182 {
183 if (SDIM == 2)
184 {
185 real_t Jloc[4], Jinv[4];
186 Jloc[0] = j(qx,qy,0,0,e);
187 Jloc[1] = j(qx,qy,1,0,e);
188 Jloc[2] = j(qx,qy,0,1,e);
189 Jloc[3] = j(qx,qy,1,1,e);
190 kernels::CalcInverse<2>(Jloc, Jinv);
191 const real_t U = Jinv[0]*du[0] + Jinv[1]*du[1];
192 const real_t V = Jinv[2]*du[0] + Jinv[3]*du[1];
193 du[0] = U;
194 du[1] = V;
195 }
196 else
197 {
198 real_t Jloc[6], Jinv[6];
199 Jloc[0] = j(qx,qy,0,0,e);
200 Jloc[1] = j(qx,qy,1,0,e);
201 Jloc[2] = j(qx,qy,2,0,e);
202 Jloc[3] = j(qx,qy,0,1,e);
203 Jloc[4] = j(qx,qy,1,1,e);
204 Jloc[5] = j(qx,qy,2,1,e);
205 kernels::CalcLeftInverse<3,2>(Jloc, Jinv);
206 const real_t U = Jinv[0]*du[0] + Jinv[1]*du[1];
207 const real_t V = Jinv[2]*du[0] + Jinv[3]*du[1];
208 const real_t W = Jinv[4]*du[0] + Jinv[5]*du[1];
209 du[0] = U;
210 du[1] = V;
211 du[2] = W;
212 }
213 }
214 for (int d = 0; d < SDIM; ++d)
215 {
216 if (Q_LAYOUT == QVectorLayout::byVDIM)
217 {
218 y(c,d,qx,qy,e) = du[d];
219 }
220 else // Q_LAYOUT == QVectorLayout::byNODES
221 {
222 y(qx,qy,c,d,e) = du[d];
223 }
224 }
225 }
226 }
227 MFEM_SYNC_THREAD;
228 }
229 });
230}
231
232// Template compute kernel for derivatives in 3D: tensor product version.
233template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
234 int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0>
235static void Derivatives3D(const int NE,
236 const real_t *b_,
237 const real_t *g_,
238 const real_t *j_,
239 const real_t *x_,
240 real_t *y_,
241 const int sdim = 3,
242 const int vdim = 0,
243 const int d1d = 0,
244 const int q1d = 0)
245{
246 const int D1D = T_D1D ? T_D1D : d1d;
247 const int Q1D = T_Q1D ? T_Q1D : q1d;
248 const int VDIM = T_VDIM ? T_VDIM : vdim;
249
250 const auto b = Reshape(b_, Q1D, D1D);
251 const auto g = Reshape(g_, Q1D, D1D);
252 const auto j = Reshape(j_, Q1D, Q1D, Q1D, 3, 3, NE);
253 const auto x = Reshape(x_, D1D, D1D, D1D, VDIM, NE);
254 auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
255 Reshape(y_, Q1D, Q1D, Q1D, VDIM, 3, NE):
256 Reshape(y_, VDIM, 3, Q1D, Q1D, Q1D, NE);
257
258 mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
259 {
260 const int D1D = T_D1D ? T_D1D : d1d;
261 const int Q1D = T_Q1D ? T_Q1D : q1d;
262 const int VDIM = T_VDIM ? T_VDIM : vdim;
263 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_INTERP_1D;
264 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_INTERP_1D;
265
266 MFEM_SHARED real_t BG[2][MQ1*MD1];
267 kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,BG);
268 DeviceMatrix B(BG[0], D1D, Q1D);
269 DeviceMatrix G(BG[1], D1D, Q1D);
270
271 MFEM_SHARED real_t sm0[3][MQ1*MQ1*MQ1];
272 MFEM_SHARED real_t sm1[3][MQ1*MQ1*MQ1];
273 DeviceTensor<3> X(sm0[2], D1D, D1D, D1D);
274 DeviceTensor<3> DDQ0(sm0[0], D1D, D1D, Q1D);
275 DeviceTensor<3> DDQ1(sm0[1], D1D, D1D, Q1D);
276 DeviceTensor<3> DQQ0(sm1[0], D1D, Q1D, Q1D);
277 DeviceTensor<3> DQQ1(sm1[1], D1D, Q1D, Q1D);
278 DeviceTensor<3> DQQ2(sm1[2], D1D, Q1D, Q1D);
279
280 for (int c = 0; c < VDIM; ++c)
281 {
282 kernels::internal::LoadX(e,D1D,c,x,X);
283 MFEM_FOREACH_THREAD(dz,z,D1D)
284 {
285 MFEM_FOREACH_THREAD(dy,y,D1D)
286 {
287 MFEM_FOREACH_THREAD(qx,x,Q1D)
288 {
289 real_t u = 0.0;
290 real_t v = 0.0;
291 for (int dx = 0; dx < D1D; ++dx)
292 {
293 const real_t input = X(dx,dy,dz);
294 u += input * B(dx,qx);
295 v += input * G(dx,qx);
296 }
297 DDQ0(dz,dy,qx) = u;
298 DDQ1(dz,dy,qx) = v;
299 }
300 }
301 }
302 MFEM_SYNC_THREAD;
303 MFEM_FOREACH_THREAD(dz,z,D1D)
304 {
305 MFEM_FOREACH_THREAD(qy,y,Q1D)
306 {
307 MFEM_FOREACH_THREAD(qx,x,Q1D)
308 {
309 real_t u = 0.0;
310 real_t v = 0.0;
311 real_t w = 0.0;
312 for (int dy = 0; dy < D1D; ++dy)
313 {
314 u += DDQ1(dz,dy,qx) * B(dy,qy);
315 v += DDQ0(dz,dy,qx) * G(dy,qy);
316 w += DDQ0(dz,dy,qx) * B(dy,qy);
317 }
318 DQQ0(dz,qy,qx) = u;
319 DQQ1(dz,qy,qx) = v;
320 DQQ2(dz,qy,qx) = w;
321 }
322 }
323 }
324 MFEM_SYNC_THREAD;
325 MFEM_FOREACH_THREAD(qz,z,Q1D)
326 {
327 MFEM_FOREACH_THREAD(qy,y,Q1D)
328 {
329 MFEM_FOREACH_THREAD(qx,x,Q1D)
330 {
331 real_t u = 0.0;
332 real_t v = 0.0;
333 real_t w = 0.0;
334 for (int dz = 0; dz < D1D; ++dz)
335 {
336 u += DQQ0(dz,qy,qx) * B(dz,qz);
337 v += DQQ1(dz,qy,qx) * B(dz,qz);
338 w += DQQ2(dz,qy,qx) * G(dz,qz);
339 }
340 if (GRAD_PHYS)
341 {
342 real_t Jloc[9], Jinv[9];
343 for (int col = 0; col < 3; col++)
344 {
345 for (int row = 0; row < 3; row++)
346 {
347 Jloc[row+3*col] = j(qx,qy,qz,row,col,e);
348 }
349 }
350 kernels::CalcInverse<3>(Jloc, Jinv);
351 const real_t U = Jinv[0]*u + Jinv[1]*v + Jinv[2]*w;
352 const real_t V = Jinv[3]*u + Jinv[4]*v + Jinv[5]*w;
353 const real_t W = Jinv[6]*u + Jinv[7]*v + Jinv[8]*w;
354 u = U; v = V; w = W;
355 }
356 if (Q_LAYOUT == QVectorLayout::byVDIM)
357 {
358 y(c,0,qx,qy,qz,e) = u;
359 y(c,1,qx,qy,qz,e) = v;
360 y(c,2,qx,qy,qz,e) = w;
361 }
362 if (Q_LAYOUT == QVectorLayout::byNODES)
363 {
364 y(qx,qy,qz,c,0,e) = u;
365 y(qx,qy,qz,c,1,e) = v;
366 y(qx,qy,qz,c,2,e) = w;
367 }
368 }
369 }
370 }
371 MFEM_SYNC_THREAD;
372 }
373 });
374}
375
376template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS>
377static void CollocatedDerivatives1D(const int NE,
378 const real_t *g_,
379 const real_t *j_,
380 const real_t *x_,
381 real_t *y_,
382 const int sdim,
383 const int vdim,
384 const int d1d)
385{
386 Derivatives1D<Q_LAYOUT, GRAD_PHYS>(
387 NE, nullptr, g_, j_, x_, y_, sdim, vdim, d1d, d1d);
388}
389
390// Template compute kernel for derivatives in 2D: tensor product version.
391template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
392 int T_VDIM = 0, int T_D1D = 0,
393 int T_NBZ = 1>
394static void CollocatedDerivatives2D(const int NE,
395 const real_t *g_,
396 const real_t *j_,
397 const real_t *x_,
398 real_t *y_,
399 const int sdim = 2,
400 const int vdim = 0,
401 const int d1d = 0)
402{
403 const int D1D = T_D1D ? T_D1D : d1d;
404 const int VDIM = T_VDIM ? T_VDIM : vdim;
405 const int SDIM = GRAD_PHYS ? sdim : 2;
406 static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
407
408 const auto g = Reshape(g_, D1D, D1D);
409 const auto j = Reshape(j_, D1D, D1D, SDIM, 2, NE);
410 const auto x = Reshape(x_, D1D, D1D, VDIM, NE);
411 auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
412 Reshape(y_, D1D, D1D, VDIM, SDIM, NE):
413 Reshape(y_, VDIM, SDIM, D1D, D1D, NE);
414
415 mfem::forall_2D_batch(NE, D1D, D1D, NBZ, [=] MFEM_HOST_DEVICE (int e)
416 {
417 const int D1D = T_D1D ? T_D1D : d1d;
418 const int VDIM = T_VDIM ? T_VDIM : vdim;
419 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
420
421 const int tidz = MFEM_THREAD_ID(z);
422
423 MFEM_SHARED real_t XY[NBZ][MD1*MD1];
424 DeviceTensor<2> X((real_t*)(XY+tidz), D1D, D1D);
425
426 for (int c = 0; c < VDIM; ++c)
427 {
428 kernels::internal::LoadX<MD1,NBZ>(e,D1D,c,x,XY);
429 MFEM_FOREACH_THREAD(dy,y,D1D)
430 {
431 MFEM_FOREACH_THREAD(dx,x,D1D)
432 {
433 real_t u = 0.0;
434 real_t v = 0.0;
435 real_t w = 0.0;
436 for (int dxy = 0; dxy < D1D; ++dxy)
437 {
438 u += X(dxy, dy) * g(dx,dxy);
439 v += X(dx, dxy) * g(dy,dxy);
440 }
441
442 if (GRAD_PHYS)
443 {
444 if (SDIM == 2)
445 {
446 real_t Jloc[4], Jinv[4];
447 Jloc[0] = j(dx,dy,0,0,e);
448 Jloc[1] = j(dx,dy,1,0,e);
449 Jloc[2] = j(dx,dy,0,1,e);
450 Jloc[3] = j(dx,dy,1,1,e);
451 kernels::CalcInverse<2>(Jloc, Jinv);
452 const real_t U = Jinv[0]*u + Jinv[1]*v;
453 const real_t V = Jinv[2]*u + Jinv[3]*v;
454 u = U;
455 v = V;
456 }
457 else
458 {
459 real_t Jloc[6], Jinv[6];
460 Jloc[0] = j(dx,dy,0,0,e);
461 Jloc[1] = j(dx,dy,1,0,e);
462 Jloc[2] = j(dx,dy,2,0,e);
463 Jloc[3] = j(dx,dy,0,1,e);
464 Jloc[4] = j(dx,dy,1,1,e);
465 Jloc[5] = j(dx,dy,2,1,e);
466 kernels::CalcLeftInverse<3,2>(Jloc, Jinv);
467 const real_t U = Jinv[0]*u + Jinv[1]*v;
468 const real_t V = Jinv[2]*u + Jinv[3]*v;
469 const real_t W = Jinv[4]*u + Jinv[5]*v;
470 u = U;
471 v = V;
472 w = W;
473 }
474 }
475
476 if (Q_LAYOUT == QVectorLayout::byVDIM)
477 {
478 y(c,0,dx,dy,e) = u;
479 y(c,1,dx,dy,e) = v;
480 if (SDIM == 3) { y(c,2,dx,dy,e) = w; }
481 }
482 if (Q_LAYOUT == QVectorLayout::byNODES)
483 {
484 y(dx,dy,c,0,e) = u;
485 y(dx,dy,c,1,e) = v;
486 if (SDIM == 3) { y(dx,dy,c,2,e) = w; }
487 }
488 }
489 }
490 MFEM_SYNC_THREAD;
491 }
492 });
493}
494
495// Template compute kernel for derivatives in 3D: tensor product version.
496template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
497 int T_VDIM = 0, int T_D1D = 0>
498static void CollocatedDerivatives3D(const int NE,
499 const real_t *g_,
500 const real_t *j_,
501 const real_t *x_,
502 real_t *y_,
503 const int sdim = 3,
504 const int vdim = 0,
505 const int d1d = 0)
506{
507 MFEM_VERIFY(sdim == 3, "");
508
509 const int D1D = T_D1D ? T_D1D : d1d;
510 const int VDIM = T_VDIM ? T_VDIM : vdim;
511
512 const auto g = Reshape(g_, D1D, D1D);
513 const auto j = Reshape(j_, D1D, D1D, D1D, 3, 3, NE);
514 const auto x = Reshape(x_, D1D, D1D, D1D, VDIM, NE);
515 auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
516 Reshape(y_, D1D, D1D, D1D, VDIM, 3, NE):
517 Reshape(y_, VDIM, 3, D1D, D1D, D1D, NE);
518
519 mfem::forall_3D(NE, D1D, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
520 {
521 const int D1D = T_D1D ? T_D1D : d1d;
522 const int VDIM = T_VDIM ? T_VDIM : vdim;
523 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_INTERP_1D;
524
525 MFEM_SHARED real_t uvw[MD1*MD1*MD1];
526 DeviceTensor<3> X(uvw, D1D, D1D, D1D);
527
528 for (int c = 0; c < VDIM; ++c)
529 {
530 kernels::internal::LoadX(e,D1D,c,x,X);
531 MFEM_FOREACH_THREAD(dz,z,D1D)
532 {
533 MFEM_FOREACH_THREAD(dy,y,D1D)
534 {
535 MFEM_FOREACH_THREAD(dx,x,D1D)
536 {
537 real_t u = 0.0;
538 real_t v = 0.0;
539 real_t w = 0.0;
540 for (int dxyz = 0; dxyz < D1D; ++dxyz)
541 {
542 u += X(dxyz, dy, dz) * g(dx,dxyz);
543 v += X(dx, dxyz, dz) * g(dy,dxyz);
544 w += X(dx, dy, dxyz) * g(dz,dxyz);
545 }
546
547 if (GRAD_PHYS)
548 {
549 real_t Jloc[9], Jinv[9];
550 for (int col = 0; col < 3; col++)
551 {
552 for (int row = 0; row < 3; row++)
553 {
554 Jloc[row+3*col] = j(dx,dy,dz,row,col,e);
555 }
556 }
557 kernels::CalcInverse<3>(Jloc, Jinv);
558 const real_t U = Jinv[0]*u + Jinv[1]*v + Jinv[2]*w;
559 const real_t V = Jinv[3]*u + Jinv[4]*v + Jinv[5]*w;
560 const real_t W = Jinv[6]*u + Jinv[7]*v + Jinv[8]*w;
561 u = U; v = V; w = W;
562 }
563 if (Q_LAYOUT == QVectorLayout::byVDIM)
564 {
565 y(c,0,dx,dy,dz,e) = u;
566 y(c,1,dx,dy,dz,e) = v;
567 y(c,2,dx,dy,dz,e) = w;
568 }
569 if (Q_LAYOUT == QVectorLayout::byNODES)
570 {
571 y(dx,dy,dz,c,0,e) = u;
572 y(dx,dy,dz,c,1,e) = v;
573 y(dx,dy,dz,c,2,e) = w;
574 }
575
576 }
577 }
578 }
579 MFEM_SYNC_THREAD;
580 }
581 });
582}
583
584} // namespace quadrature_interpolator
585
586} // namespace internal
587
588/// @cond Suppress_Doxygen_warnings
589
590template<int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
591 int Q1D, int NBZ>
592QuadratureInterpolator::GradKernelType
593QuadratureInterpolator::GradKernels::Kernel()
594{
595 if (DIM == 1) { return internal::quadrature_interpolator::Derivatives1D<Q_LAYOUT, GRAD_PHYS>; }
596 else if (DIM == 2) { return internal::quadrature_interpolator::Derivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D, NBZ>; }
597 else if (DIM == 3) { return internal::quadrature_interpolator::Derivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D>; }
598 else { MFEM_ABORT(""); }
599}
600
601template<int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int VDIM, int D1D,
602 int NBZ>
603QuadratureInterpolator::CollocatedGradKernelType
604QuadratureInterpolator::CollocatedGradKernels::Kernel()
605{
606 if (DIM == 1) { return internal::quadrature_interpolator::CollocatedDerivatives1D<Q_LAYOUT, GRAD_PHYS>; }
607 else if (DIM == 2) { return internal::quadrature_interpolator::CollocatedDerivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, NBZ>; }
608 else if (DIM == 3) { return internal::quadrature_interpolator::CollocatedDerivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D>; }
609 else { MFEM_ABORT(""); }
610}
611
612/// @endcond
613
614} // namespace mfem
615
616#endif
mfem::QuadratureInterpolator::CollocatedGradKernelType
void(*)(const int, const real_t *, const real_t *, const real_t *, real_t *, const int, const int, const int) CollocatedGradKernelType
Definition quadinterpolator.hpp:160
mfem::QuadratureInterpolator::GradKernelType
void(*)(const int, const real_t *, const real_t *, const real_t *, const real_t *, real_t *, const int, const int, const int, const int) GradKernelType
Definition quadinterpolator.hpp:157
b
real_t b
Definition lissajous.cpp:42
SDIM
constexpr int SDIM
Definition minimal-surface.cpp:73
DIM
constexpr int DIM
Definition minimal-surface.cpp:72
mfem::kernels::CalcInverse
MFEM_HOST_DEVICE void CalcInverse(const T *data, T *inv_data)
Return the inverse of a matrix with given size and data into the matrix with data inv_data.
Definition kernels.hpp:246
mfem::kernels::CalcLeftInverse< 2, 1 >
MFEM_HOST_DEVICE void CalcLeftInverse< 2, 1 >(const real_t *d, real_t *left_inv)
Definition kernels.hpp:1140
mfem::kernels::CalcLeftInverse< 3, 2 >
MFEM_HOST_DEVICE void CalcLeftInverse< 3, 2 >(const real_t *d, real_t *left_inv)
Definition kernels.hpp:1157
mfem::kernels::CalcLeftInverse< 3, 1 >
MFEM_HOST_DEVICE void CalcLeftInverse< 3, 1 >(const real_t *d, real_t *left_inv)
Definition kernels.hpp:1148
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::DeviceMatrix
DeviceTensor< 2, real_t > DeviceMatrix
Definition dtensor.hpp:143
mfem::Reshape
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:131
mfem::forall_2D_batch
void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
Definition forall.hpp:768
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:774
mfem::QVectorLayout
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition fespace.hpp:53
mfem::real_t
float real_t
Definition config.hpp:43
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:753
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0