MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_hdiv_kernels.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2026, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_BILININTEG_HDIV_KERNELS_HPP
13#define MFEM_BILININTEG_HDIV_KERNELS_HPP
14
15#include "../../config/config.hpp"
16#include "../../general/array.hpp"
17#include "../../general/forall.hpp"
18#include "../../linalg/dtensor.hpp"
19#include "../../linalg/vector.hpp"
20#include "../bilininteg.hpp"
21
22// Piola transformation in H(div): w = (1 / det (dF)) dF \hat{w}
23// div w = (1 / det (dF)) \hat{div} \hat{w}
24
25namespace mfem
26{
27
28namespace internal
29{
30
31// PA H(div) Mass Assemble 2D kernel
32void PAHdivMassSetup2D(const int Q1D,
33 const int coeffDim,
34 const int NE,
35 const Array<real_t> &w,
36 const Vector &j,
37 Vector &coeff_,
38 Vector &op);
39
40// PA H(div) Mass Assemble 3D kernel
41void PAHdivMassSetup3D(const int Q1D,
42 const int coeffDim,
43 const int NE,
44 const Array<real_t> &w,
45 const Vector &j,
46 Vector &coeff_,
47 Vector &op);
48
49// PA H(div) Mass Diagonal 2D kernel
50void PAHdivMassAssembleDiagonal2D(const int D1D,
51 const int Q1D,
52 const int NE,
53 const bool symmetric,
54 const Array<real_t> &Bo_,
55 const Array<real_t> &Bc_,
56 const Vector &op_,
57 Vector &diag_);
58
59// PA H(div) Mass Diagonal 3D kernel
60void PAHdivMassAssembleDiagonal3D(const int D1D,
61 const int Q1D,
62 const int NE,
63 const bool symmetric,
64 const Array<real_t> &Bo_,
65 const Array<real_t> &Bc_,
66 const Vector &op_,
67 Vector &diag_);
68
69// PA H(div) Mass Apply 2D kernel
70void PAHdivMassApply2D(const int NE, const bool symmetric,
71 const bool scalar_coeff, const Array<real_t> &Bo_,
72 const Array<real_t> &Bc_, const Array<real_t> &Bot_,
73 const Array<real_t> &Bct_, const Vector &op_,
74 const Vector &x_, Vector &y_, const int D1D,
75 const int TestD1D, const int Q1D);
76
77// PA H(div) Mass Apply 3D kernel
78void PAHdivMassApply3D(const int NE, const bool symmetric,
79 const bool scalar_coeff, const Array<real_t> &Bo_,
80 const Array<real_t> &Bc_, const Array<real_t> &Bot_,
81 const Array<real_t> &Bct_, const Vector &op_,
82 const Vector &x_, Vector &y_, const int D1D,
83 const int TestD1D, const int Q1D);
84
85// Shared memory PA H(div) Mass Apply 2D kernel
86template <int T_D1D = 0, int T_Q1D = 0>
87inline void SmemPAHdivMassApply2D(
88 const int NE, const bool symmetric, const bool, const Array<real_t> &Bo_,
89 const Array<real_t> &Bc_, const Array<real_t> &Bot_,
90 const Array<real_t> &Bct_, const Vector &op_, const Vector &x_, Vector &y_,
91 const int d1d = 0, const int = 0, const int q1d = 0)
92{
93 MFEM_CONTRACT_VAR(Bot_);
94 MFEM_CONTRACT_VAR(Bct_);
95
96 static constexpr int VDIM = 2;
97
98 const int D1D = T_D1D ? T_D1D : d1d;
99 const int Q1D = T_Q1D ? T_Q1D : q1d;
100
101 const auto bo = Reshape(Bo_.Read(), Q1D, D1D-1);
102 const auto bc = Reshape(Bc_.Read(), Q1D, D1D);
103 const auto D = Reshape(op_.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
104 const auto x = Reshape(x_.Read(), D1D*(D1D-1), VDIM, NE);
105 auto y = y_.ReadWrite();
106
107 mfem::forall_3D(NE, Q1D, Q1D, VDIM, [=] MFEM_HOST_DEVICE (int e)
108 {
109 const int tidz = MFEM_THREAD_ID(z);
110
111 const int D1D = T_D1D ? T_D1D : d1d;
112 const int Q1D = T_Q1D ? T_Q1D : q1d;
113
114 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
115 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
116 constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
117
118 MFEM_SHARED real_t smo[MQ1*(MD1-1)];
119 DeviceMatrix Bo(smo, D1D-1, Q1D);
120
121 MFEM_SHARED real_t smc[MQ1*MD1];
122 DeviceMatrix Bc(smc, D1D, Q1D);
123
124 MFEM_SHARED real_t sm0[VDIM*MDQ*MDQ];
125 MFEM_SHARED real_t sm1[VDIM*MDQ*MDQ];
126 DeviceMatrix X(sm0, D1D*(D1D-1), VDIM);
127 DeviceCube QD(sm1, Q1D, D1D, VDIM);
128 DeviceCube QQ(sm0, Q1D, Q1D, VDIM);
129 DeviceCube DQ(sm1, D1D, Q1D, VDIM);
130
131 // Load X, Bo and Bc into shared memory
132 MFEM_FOREACH_THREAD(vd,z,VDIM)
133 {
134 MFEM_FOREACH_THREAD(dy,y,D1D)
135 {
136 MFEM_FOREACH_THREAD(qx,x,Q1D)
137 {
138 if (qx < D1D && dy < (D1D-1))
139 {
140 X(qx + dy*D1D,vd) = x(qx+dy*D1D,vd,e);
141 }
142 if (tidz == 0)
143 {
144 if (dy < (D1D-1)) { Bo(dy,qx) = bo(qx,dy); }
145 Bc(dy,qx) = bc(qx,dy);
146 }
147 }
148 }
149 }
150 MFEM_SYNC_THREAD;
151 // Apply B operator
152 MFEM_FOREACH_THREAD(vd,z,VDIM)
153 {
154 const int nx = (vd == 0) ? D1D : D1D-1;
155 const int ny = (vd == 1) ? D1D : D1D-1;
156 DeviceCube Xxy(X, nx, ny, VDIM);
157 DeviceMatrix Bx = (vd == 0) ? Bc : Bo;
158 MFEM_FOREACH_THREAD(dy,y,ny)
159 {
160 MFEM_FOREACH_THREAD(qx,x,Q1D)
161 {
162 real_t dq = 0.0;
163 for (int dx = 0; dx < nx; ++dx)
164 {
165 dq += Xxy(dx,dy,vd) * Bx(dx,qx);
166 }
167 QD(qx,dy,vd) = dq;
168 }
169 }
170 }
171 MFEM_SYNC_THREAD;
172 MFEM_FOREACH_THREAD(vd,z,VDIM)
173 {
174 const int ny = (vd == 1) ? D1D : D1D-1;
175 DeviceMatrix By = (vd == 1) ? Bc : Bo;
176 MFEM_FOREACH_THREAD(qy,y,Q1D)
177 {
178 MFEM_FOREACH_THREAD(qx,x,Q1D)
179 {
180 real_t qq = 0.0;
181 for (int dy = 0; dy < ny; ++dy)
182 {
183 qq += QD(qx,dy,vd) * By(dy,qy);
184 }
185 QQ(qx,qy,vd) = qq;
186 }
187 }
188 }
189 MFEM_SYNC_THREAD;
190 // Apply D operator
191 if (tidz == 0)
192 {
193 MFEM_FOREACH_THREAD(qy,y,Q1D)
194 {
195 MFEM_FOREACH_THREAD(qx,x,Q1D)
196 {
197 const real_t Qx = QQ(qx,qy,0);
198 const real_t Qy = QQ(qx,qy,1);
199
200 const real_t D11 = D(qx,qy,0,e);
201 const real_t D12 = D(qx,qy,1,e);
202 const real_t D21 = symmetric ? D12 : D(qx,qy,2,e);
203 const real_t D22 = symmetric ? D(qx,qy,2,e) : D(qx,qy,3,e);
204
205 QQ(qx,qy,0) = D11*Qx + D12*Qy;
206 QQ(qx,qy,1) = D21*Qx + D22*Qy;
207 }
208 }
209 }
210 MFEM_SYNC_THREAD;
211 // Apply Bt operator
212 MFEM_FOREACH_THREAD(vd,z,VDIM)
213 {
214 const int nx = (vd == 0) ? D1D : D1D-1;
215 DeviceMatrix Btx = (vd == 0) ? Bc : Bo;
216 MFEM_FOREACH_THREAD(qy,y,Q1D)
217 {
218 MFEM_FOREACH_THREAD(dx,x,nx)
219 {
220 real_t qd = 0.0;
221 for (int qx = 0; qx < Q1D; ++qx)
222 {
223 qd += QQ(qx,qy,vd) * Btx(dx,qx);
224 }
225 DQ(dx,qy,vd) = qd;
226 }
227 }
228 }
229 MFEM_SYNC_THREAD;
230 MFEM_FOREACH_THREAD(vd,z,VDIM)
231 {
232 const int nx = (vd == 0) ? D1D : D1D-1;
233 const int ny = (vd == 1) ? D1D : D1D-1;
234 DeviceMatrix Bty = (vd == 1) ? Bc : Bo;
235 DeviceTensor<4> Yxy(y, nx, ny, VDIM, NE);
236 MFEM_FOREACH_THREAD(dy,y,ny)
237 {
238 MFEM_FOREACH_THREAD(dx,x,nx)
239 {
240 real_t dd = 0.0;
241 for (int qy = 0; qy < Q1D; ++qy)
242 {
243 dd += DQ(dx,qy,vd) * Bty(dy,qy);
244 }
245 Yxy(dx,dy,vd,e) += dd;
246 }
247 }
248 }
249 MFEM_SYNC_THREAD;
250 });
251}
252
253// Shared memory PA H(div) Mass Apply 3D kernel
254template <int T_D1D = 0, int T_Q1D = 0>
255inline void
256SmemPAHdivMassApply3D(const int NE, const bool symmetric, const bool,
257 const Array<real_t> &Bo_, const Array<real_t> &Bc_,
258 const Array<real_t> &Bot_, const Array<real_t> &Bct_,
259 const Vector &op_, const Vector &x_, Vector &y_,
260 const int d1d = 0, const int = 0, const int q1d = 0)
261{
262 MFEM_CONTRACT_VAR(Bot_);
263 MFEM_CONTRACT_VAR(Bct_);
264
265 static constexpr int VDIM = 3;
266
267 const int D1D = T_D1D ? T_D1D : d1d;
268 const int Q1D = T_Q1D ? T_Q1D : q1d;
269
270 const auto bo = Reshape(Bo_.Read(), Q1D, D1D-1);
271 const auto bc = Reshape(Bc_.Read(), Q1D, D1D);
272 const auto D = Reshape(op_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
273 const auto x = Reshape(x_.Read(), D1D*(D1D-1)*(D1D-1), VDIM, NE);
274 auto y = y_.ReadWrite();
275
276 mfem::forall_3D(NE, Q1D, Q1D, VDIM, [=] MFEM_HOST_DEVICE (int e)
277 {
278 const int tidz = MFEM_THREAD_ID(z);
279
280 const int D1D = T_D1D ? T_D1D : d1d;
281 const int Q1D = T_Q1D ? T_Q1D : q1d;
282
283 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
284 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
285 constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
286
287 MFEM_SHARED real_t smo[MQ1*(MD1-1)];
288 DeviceMatrix Bo(smo, D1D-1, Q1D);
289
290 MFEM_SHARED real_t smc[MQ1*MD1];
291 DeviceMatrix Bc(smc, D1D, Q1D);
292
293 MFEM_SHARED real_t sm0[VDIM*MDQ*MDQ*MDQ];
294 MFEM_SHARED real_t sm1[VDIM*MDQ*MDQ*MDQ];
295 DeviceMatrix X(sm0, D1D*(D1D-1)*(D1D-1), VDIM);
296 DeviceTensor<4> QDD(sm1, Q1D, D1D, D1D, VDIM);
297 DeviceTensor<4> QQD(sm0, Q1D, Q1D, D1D, VDIM);
298 DeviceTensor<4> QQQ(sm1, Q1D, Q1D, Q1D, VDIM);
299 DeviceTensor<4> DQQ(sm0, D1D, Q1D, Q1D, VDIM);
300 DeviceTensor<4> DDQ(sm1, D1D, D1D, Q1D, VDIM);
301
302 // Load X into shared memory
303 MFEM_FOREACH_THREAD(vd,z,VDIM)
304 {
305 MFEM_FOREACH_THREAD(dz,y,D1D-1)
306 {
307 MFEM_FOREACH_THREAD(dy,x,D1D-1)
308 {
309 MFEM_UNROLL(MD1)
310 for (int dx = 0; dx < D1D; ++dx)
311 {
312 X(dx+(dy+dz*(D1D-1))*D1D,vd) = x(dx+(dy+dz*(D1D-1))*D1D,vd,e);
313 }
314 }
315 }
316 }
317 // Load Bo and Bc into shared memory
318 if (tidz == 0)
319 {
320 MFEM_FOREACH_THREAD(d,y,D1D-1)
321 {
322 MFEM_FOREACH_THREAD(q,x,Q1D)
323 {
324 Bo(d,q) = bo(q,d);
325 }
326 }
327 MFEM_FOREACH_THREAD(d,y,D1D)
328 {
329 MFEM_FOREACH_THREAD(q,x,Q1D)
330 {
331 Bc(d,q) = bc(q,d);
332 }
333 }
334 }
335 MFEM_SYNC_THREAD;
336 // Apply B operator
337 MFEM_FOREACH_THREAD(vd,z,VDIM)
338 {
339 const int nx = (vd == 0) ? D1D : D1D-1;
340 const int ny = (vd == 1) ? D1D : D1D-1;
341 const int nz = (vd == 2) ? D1D : D1D-1;
342 DeviceTensor<4> Xxyz(X, nx, ny, nz, VDIM);
343 DeviceMatrix Bx = (vd == 0) ? Bc : Bo;
344 MFEM_FOREACH_THREAD(dy,y,ny)
345 {
346 MFEM_FOREACH_THREAD(qx,x,Q1D)
347 {
348 real_t u[D1D];
349 MFEM_UNROLL(MD1)
350 for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
351 MFEM_UNROLL(MD1)
352 for (int dx = 0; dx < nx; ++dx)
353 {
354 MFEM_UNROLL(MD1)
355 for (int dz = 0; dz < nz; ++dz)
356 {
357 u[dz] += Xxyz(dx,dy,dz,vd) * Bx(dx,qx);
358 }
359 }
360 MFEM_UNROLL(MD1)
361 for (int dz = 0; dz < nz; ++dz) { QDD(qx,dy,dz,vd) = u[dz]; }
362 }
363 }
364 }
365 MFEM_SYNC_THREAD;
366 MFEM_FOREACH_THREAD(vd,z,VDIM)
367 {
368 const int ny = (vd == 1) ? D1D : D1D-1;
369 const int nz = (vd == 2) ? D1D : D1D-1;
370 DeviceMatrix By = (vd == 1) ? Bc : Bo;
371 MFEM_FOREACH_THREAD(qy,y,Q1D)
372 {
373 MFEM_FOREACH_THREAD(qx,x,Q1D)
374 {
375 real_t u[D1D];
376 MFEM_UNROLL(MD1)
377 for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
378 MFEM_UNROLL(MD1)
379 for (int dy = 0; dy < ny; ++dy)
380 {
381 MFEM_UNROLL(MD1)
382 for (int dz = 0; dz < nz; ++dz)
383 {
384 u[dz] += QDD(qx,dy,dz,vd) * By(dy,qy);
385 }
386 }
387 MFEM_UNROLL(MD1)
388 for (int dz = 0; dz < nz; ++dz) { QQD(qx,qy,dz,vd) = u[dz]; }
389 }
390 }
391 }
392 MFEM_SYNC_THREAD;
393 MFEM_FOREACH_THREAD(vd,z,VDIM)
394 {
395 const int nz = (vd == 2) ? D1D : D1D-1;
396 DeviceMatrix Bz = (vd == 2) ? Bc : Bo;
397 MFEM_FOREACH_THREAD(qy,y,Q1D)
398 {
399 MFEM_FOREACH_THREAD(qx,x,Q1D)
400 {
401 real_t u[Q1D];
402 MFEM_UNROLL(MQ1)
403 for (int qz = 0; qz < Q1D; ++qz) { u[qz] = 0.0; }
404 MFEM_UNROLL(MD1)
405 for (int dz = 0; dz < nz; ++dz)
406 {
407 MFEM_UNROLL(MQ1)
408 for (int qz = 0; qz < Q1D; ++qz)
409 {
410 u[qz] += QQD(qx,qy,dz,vd) * Bz(dz,qz);
411 }
412 }
413 MFEM_UNROLL(MQ1)
414 for (int qz = 0; qz < Q1D; ++qz) { QQQ(qx,qy,qz,vd) = u[qz]; }
415 }
416 }
417 }
418 MFEM_SYNC_THREAD;
419 // Apply D operator
420 if (tidz == 0)
421 {
422 MFEM_FOREACH_THREAD(qy,y,Q1D)
423 {
424 MFEM_FOREACH_THREAD(qx,x,Q1D)
425 {
426 MFEM_UNROLL(MQ1)
427 for (int qz = 0; qz < Q1D; ++qz)
428 {
429 const real_t Qx = QQQ(qx,qy,qz,0);
430 const real_t Qy = QQQ(qx,qy,qz,1);
431 const real_t Qz = QQQ(qx,qy,qz,2);
432
433 const real_t D11 = D(qx,qy,qz,0,e);
434 const real_t D12 = D(qx,qy,qz,1,e);
435 const real_t D13 = D(qx,qy,qz,2,e);
436 const real_t D21 = symmetric ? D12 : D(qx,qy,qz,3,e);
437 const real_t D22 = symmetric ? D(qx,qy,qz,3,e) : D(qx,qy,qz,4,e);
438 const real_t D23 = symmetric ? D(qx,qy,qz,4,e) : D(qx,qy,qz,5,e);
439 const real_t D31 = symmetric ? D13 : D(qx,qy,qz,6,e);
440 const real_t D32 = symmetric ? D23 : D(qx,qy,qz,7,e);
441 const real_t D33 = symmetric ? D(qx,qy,qz,5,e) : D(qx,qy,qz,8,e);
442
443 QQQ(qx,qy,qz,0) = D11*Qx + D12*Qy + D13*Qz;
444 QQQ(qx,qy,qz,1) = D21*Qx + D22*Qy + D23*Qz;
445 QQQ(qx,qy,qz,2) = D31*Qx + D32*Qy + D33*Qz;
446 }
447 }
448 }
449 }
450 MFEM_SYNC_THREAD;
451 // Apply Bt operator
452 MFEM_FOREACH_THREAD(vd,z,VDIM)
453 {
454 const int nx = (vd == 0) ? D1D : D1D-1;
455 DeviceMatrix Btx = (vd == 0) ? Bc : Bo;
456 MFEM_FOREACH_THREAD(qy,y,Q1D)
457 {
458 MFEM_FOREACH_THREAD(dx,x,nx)
459 {
460 real_t u[Q1D];
461 MFEM_UNROLL(MQ1)
462 for (int qz = 0; qz < Q1D; ++qz) { u[qz] = 0.0; }
463 MFEM_UNROLL(MQ1)
464 for (int qx = 0; qx < Q1D; ++qx)
465 {
466 MFEM_UNROLL(MQ1)
467 for (int qz = 0; qz < Q1D; ++qz)
468 {
469 u[qz] += QQQ(qx,qy,qz,vd) * Btx(dx,qx);
470 }
471 }
472 MFEM_UNROLL(MQ1)
473 for (int qz = 0; qz < Q1D; ++qz) { DQQ(dx,qy,qz,vd) = u[qz]; }
474 }
475 }
476 }
477 MFEM_SYNC_THREAD;
478 MFEM_FOREACH_THREAD(vd,z,VDIM)
479 {
480 const int nx = (vd == 0) ? D1D : D1D-1;
481 const int ny = (vd == 1) ? D1D : D1D-1;
482 DeviceMatrix Bty = (vd == 1) ? Bc : Bo;
483 MFEM_FOREACH_THREAD(dy,y,ny)
484 {
485 MFEM_FOREACH_THREAD(dx,x,nx)
486 {
487 real_t u[Q1D];
488 MFEM_UNROLL(MQ1)
489 for (int qz = 0; qz < Q1D; ++qz) { u[qz] = 0.0; }
490 MFEM_UNROLL(MQ1)
491 for (int qy = 0; qy < Q1D; ++qy)
492 {
493 MFEM_UNROLL(MQ1)
494 for (int qz = 0; qz < Q1D; ++qz)
495 {
496 u[qz] += DQQ(dx,qy,qz,vd) * Bty(dy,qy);
497 }
498 }
499 MFEM_UNROLL(MQ1)
500 for (int qz = 0; qz < Q1D; ++qz) { DDQ(dx,dy,qz,vd) = u[qz]; }
501 }
502 }
503 }
504 MFEM_SYNC_THREAD;
505 MFEM_FOREACH_THREAD(vd,z,VDIM)
506 {
507 const int nx = (vd == 0) ? D1D : D1D-1;
508 const int ny = (vd == 1) ? D1D : D1D-1;
509 const int nz = (vd == 2) ? D1D : D1D-1;
510 DeviceTensor<5> Yxyz(y, nx, ny, nz, VDIM, NE);
511 DeviceMatrix Btz = (vd == 2) ? Bc : Bo;
512 MFEM_FOREACH_THREAD(dy,y,ny)
513 {
514 MFEM_FOREACH_THREAD(dx,x,nx)
515 {
516 real_t u[D1D];
517 MFEM_UNROLL(MD1)
518 for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
519 MFEM_UNROLL(MQ1)
520 for (int qz = 0; qz < Q1D; ++qz)
521 {
522 MFEM_UNROLL(MD1)
523 for (int dz = 0; dz < nz; ++dz)
524 {
525 u[dz] += DDQ(dx,dy,qz,vd) * Btz(dz,qz);
526 }
527 }
528 MFEM_UNROLL(MD1)
529 for (int dz = 0; dz < nz; ++dz) { Yxyz(dx,dy,dz,vd,e) += u[dz]; }
530 }
531 }
532 }
533 MFEM_SYNC_THREAD;
534 });
535}
536
537// PA H(div) div-div Assemble 2D kernel
538void PADivDivSetup2D(const int Q1D,
539 const int NE,
540 const Array<real_t> &w,
541 const Vector &j,
542 Vector &coeff_,
543 Vector &op);
544
545// PA H(div) div-div Assemble 3D kernel
546void PADivDivSetup3D(const int Q1D,
547 const int NE,
548 const Array<real_t> &w,
549 const Vector &j,
550 Vector &coeff_,
551 Vector &op);
552
553// PA H(div) div-div Diagonal 2D kernel
554void PADivDivAssembleDiagonal2D(const int D1D,
555 const int Q1D,
556 const int NE,
557 const Array<real_t> &Bo_,
558 const Array<real_t> &Gc_,
559 const Vector &op_,
560 Vector &diag_);
561
562// PA H(div) div-div Diagonal 3D kernel
563void PADivDivAssembleDiagonal3D(const int D1D,
564 const int Q1D,
565 const int NE,
566 const Array<real_t> &Bo_,
567 const Array<real_t> &Gc_,
568 const Vector &op_,
569 Vector &diag_);
570
571// PA H(div) div-div Apply 2D kernel
572void PADivDivApply2D(const int D1D,
573 const int Q1D,
574 const int NE,
575 const Array<real_t> &Bo_,
576 const Array<real_t> &Gc_,
577 const Array<real_t> &Bot_,
578 const Array<real_t> &Gct_,
579 const Vector &op_,
580 const Vector &x_,
581 Vector &y_);
582
583// PA H(div) div-div Apply 3D kernel
584void PADivDivApply3D(const int D1D,
585 const int Q1D,
586 const int NE,
587 const Array<real_t> &Bo_,
588 const Array<real_t> &Gc_,
589 const Array<real_t> &Bot_,
590 const Array<real_t> &Gct_,
591 const Vector &op_,
592 const Vector &x_,
593 Vector &y_);
594
595// PA H(div)-L2 Assemble 2D kernel
596// if geom != nullptr, then coeff_ is divided by detJ
597void PAHdivL2Setup2D(const int Q1D, const int NE, const Array<real_t> &w,
598 Vector &coeff_, Vector &op, const GeometricFactors *geom);
599
600// PA H(div)-L2 Assemble 3D kernel
601// if geom != nullptr, then coeff_ is divided by detJ
602void PAHdivL2Setup3D(const int Q1D, const int NE, const Array<real_t> &w,
603 Vector &coeff_, Vector &op, const GeometricFactors *geom);
604
605// PA H(div)-L2 Diagonal 2D kernel
606void PAHdivL2AssembleDiagonal_ADAt_2D(const int D1D,
607 const int Q1D,
608 const int L2D1D,
609 const int NE,
610 const Array<real_t> &L2Bo_,
611 const Array<real_t> &Gct_,
612 const Array<real_t> &Bot_,
613 const Vector &op_,
614 const Vector &D_,
615 Vector &diag_);
616
617// PA H(div)-L2 Diagonal 3D kernel
618void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
619 const int Q1D,
620 const int L2D1D,
621 const int NE,
622 const Array<real_t> &L2Bo_,
623 const Array<real_t> &Gct_,
624 const Array<real_t> &Bot_,
625 const Vector &op_,
626 const Vector &D_,
627 Vector &diag_);
628
629// PA H(div)-L2 Apply 2D kernel
630void PAHdivL2Apply2D(const int D1D,
631 const int Q1D,
632 const int L2D1D,
633 const int NE,
634 const Array<real_t> &Bo_,
635 const Array<real_t> &Gc_,
636 const Array<real_t> &L2Bot_,
637 const Vector &op_,
638 const Vector &x_,
639 Vector &y_);
640
641// PA H(div)-L2 Apply Transpose 2D kernel
642void PAHdivL2ApplyTranspose2D(const int D1D,
643 const int Q1D,
644 const int L2D1D,
645 const int NE,
646 const Array<real_t> &L2Bo_,
647 const Array<real_t> &Gct_,
648 const Array<real_t> &Bot_,
649 const Vector &op_,
650 const Vector &x_,
651 Vector &y_);
652
653// PA H(div)-L2 Apply 3D kernel
654void PAHdivL2Apply3D(const int D1D,
655 const int Q1D,
656 const int L2D1D,
657 const int NE,
658 const Array<real_t> &Bo_,
659 const Array<real_t> &Gc_,
660 const Array<real_t> &L2Bot_,
661 const Vector &op_,
662 const Vector &x_,
663 Vector &y_);
664
665// PA H(div)-L2 Apply Transpose 3D kernel
666void PAHdivL2ApplyTranspose3D(const int D1D,
667 const int Q1D,
668 const int L2D1D,
669 const int NE,
670 const Array<real_t> &L2Bo_,
671 const Array<real_t> &Gct_,
672 const Array<real_t> &Bot_,
673 const Vector &op_,
674 const Vector &x_,
675 Vector &y_);
676
677} // namespace internal
678
679} // namespace mfem
680
681#endif
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
mfem::DeviceCube
DeviceTensor< 3, real_t > DeviceCube
Definition dtensor.hpp:153
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
mfem::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:138
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:1244
mfem::DeviceMatrix
DeviceTensor< 2, real_t > DeviceMatrix
Definition dtensor.hpp:150
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0