MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_hdiv_kernels.cpp
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#include "bilininteg_hdiv_kernels.hpp"
13
14namespace mfem
15{
16
17namespace internal
18{
19
20void PAHdivMassSetup2D(const int Q1D,
21 const int coeffDim,
22 const int NE,
23 const Array<real_t> &w,
24 const Vector &j,
25 Vector &coeff_,
26 Vector &op)
27{
28 const bool symmetric = (coeffDim != 4);
29 const int NQ = Q1D*Q1D;
30 auto W = w.Read();
31
32 auto J = Reshape(j.Read(), NQ, 2, 2, NE);
33 auto C = Reshape(coeff_.Read(), coeffDim, NQ, NE);
34 auto y = Reshape(op.Write(), NQ, symmetric ? 3 : 4, NE);
35
36 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
37 {
38 for (int q = 0; q < NQ; ++q)
39 {
40 const real_t J11 = J(q,0,0,e);
41 const real_t J21 = J(q,1,0,e);
42 const real_t J12 = J(q,0,1,e);
43 const real_t J22 = J(q,1,1,e);
44 const real_t c_detJ = W[q] / ((J11*J22)-(J21*J12));
45
46 // (1/detJ) J^T C J
47 if (coeffDim == 3 || coeffDim == 4) // Matrix coefficient
48 {
49 const real_t C11 = C(0,q,e);
50 const real_t C12 = C(1,q,e);
51 const real_t C21 = symmetric ? C12 : C(2,q,e);
52 const real_t C22 = symmetric ? C(2,q,e) : C(3,q,e);
53 const real_t R11 = C11*J11 + C12*J21;
54 const real_t R21 = C21*J11 + C22*J21;
55 const real_t R12 = C11*J12 + C12*J22;
56 const real_t R22 = C21*J12 + C22*J22;
57
58 y(q,0,e) = c_detJ * (J11*R11 + J21*R21); // 1,1
59 y(q,1,e) = c_detJ * (J11*R12 + J21*R22); // 1,2
60
61 if (symmetric)
62 {
63 y(q,2,e) = c_detJ * (J12*R12 + J22*R22); // 2,2
64 }
65 else
66 {
67 y(q,2,e) = c_detJ * (J12*R11 + J22*R21); // 2,1
68 y(q,3,e) = c_detJ * (J12*R12 + J22*R22); // 2,2
69 }
70 }
71 else // Vector or scalar coefficient
72 {
73 const real_t C1 = C(0,q,e);
74 const real_t C2 = (coeffDim == 2 ? C(1,q,e) : C1);
75 y(q,0,e) = c_detJ * (J11*C1*J11 + J21*C2*J21); // 1,1
76 y(q,1,e) = c_detJ * (J11*C1*J12 + J21*C2*J22); // 1,2
77 y(q,2,e) = c_detJ * (J12*C1*J12 + J22*C2*J22); // 2,2
78 }
79 }
80 });
81}
82
83void PAHdivMassSetup3D(const int Q1D,
84 const int coeffDim,
85 const int NE,
86 const Array<real_t> &w,
87 const Vector &j,
88 Vector &coeff_,
89 Vector &op)
90{
91 const bool symmetric = (coeffDim != 9);
92 const int NQ = Q1D*Q1D*Q1D;
93 auto W = w.Read();
94 auto J = Reshape(j.Read(), NQ, 3, 3, NE);
95 auto C = Reshape(coeff_.Read(), coeffDim, NQ, NE);
96 auto y = Reshape(op.Write(), NQ, symmetric ? 6 : 9, NE);
97
98 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
99 {
100 for (int q = 0; q < NQ; ++q)
101 {
102 const real_t J11 = J(q,0,0,e);
103 const real_t J21 = J(q,1,0,e);
104 const real_t J31 = J(q,2,0,e);
105 const real_t J12 = J(q,0,1,e);
106 const real_t J22 = J(q,1,1,e);
107 const real_t J32 = J(q,2,1,e);
108 const real_t J13 = J(q,0,2,e);
109 const real_t J23 = J(q,1,2,e);
110 const real_t J33 = J(q,2,2,e);
111 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
112 J21 * (J12 * J33 - J32 * J13) +
113 J31 * (J12 * J23 - J22 * J13);
114 const real_t c_detJ = W[q] / detJ;
115
116 // (1/detJ) J^T C J
117 if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
118 {
119 real_t M[3][3];
120 M[0][0] = C(0, q, e);
121 M[0][1] = C(1, q, e);
122 M[0][2] = C(2, q, e);
123 M[1][0] = (!symmetric) ? C(3, q, e) : M[0][1];
124 M[1][1] = (!symmetric) ? C(4, q, e) : C(3, q, e);
125 M[1][2] = (!symmetric) ? C(5, q, e) : C(4, q, e);
126 M[2][0] = (!symmetric) ? C(6, q, e) : M[0][2];
127 M[2][1] = (!symmetric) ? C(7, q, e) : M[1][2];
128 M[2][2] = (!symmetric) ? C(8, q, e) : C(5, q, e);
129
130 int idx = 0;
131 for (int i=0; i<3; ++i)
132 for (int j = (symmetric ? i : 0); j<3; ++j)
133 {
134 y(q,idx,e) = 0.0;
135 for (int k=0; k<3; ++k)
136 {
137 real_t MJ_kj = 0.0;
138 for (int l=0; l<3; ++l)
139 {
140 MJ_kj += M[k][l] * J(q,l,j,e);
141 }
142
143 y(q,idx,e) += J(q,k,i,e) * MJ_kj;
144 }
145
146 y(q,idx,e) *= c_detJ;
147 idx++;
148 }
149 }
150 else // Vector or scalar coefficient version
151 {
152 int idx = 0;
153 for (int i=0; i<3; ++i)
154 for (int j=i; j<3; ++j)
155 {
156 y(q,idx,e) = 0.0;
157 for (int k=0; k<3; ++k)
158 {
159 y(q,idx,e) += J(q,k,i,e) * C(coeffDim == 3 ? k : 0, q, e) * J(q,k,j,e);
160 }
161
162 y(q,idx,e) *= c_detJ;
163 idx++;
164 }
165 }
166 }
167 });
168}
169
170void PAHdivMassAssembleDiagonal2D(const int D1D,
171 const int Q1D,
172 const int NE,
173 const bool symmetric,
174 const Array<real_t> &Bo_,
175 const Array<real_t> &Bc_,
176 const Vector &op_,
177 Vector &diag_)
178{
179 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
180 auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
181 auto op = Reshape(op_.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
182 auto diag = Reshape(diag_.ReadWrite(), 2*(D1D-1)*D1D, NE);
183
184 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
185 {
186 constexpr static int VDIM = 2;
187 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
188
189 int osc = 0;
190
191 for (int c = 0; c < VDIM; ++c) // loop over x, y components
192 {
193 const int D1Dx = (c == 1) ? D1D - 1 : D1D;
194 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
195
196 for (int dy = 0; dy < D1Dy; ++dy)
197 {
198 real_t mass[MAX_Q1D];
199 for (int qx = 0; qx < Q1D; ++qx)
200 {
201 mass[qx] = 0.0;
202 for (int qy = 0; qy < Q1D; ++qy)
203 {
204 const real_t wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
205 mass[qx] += wy*wy*((c == 0) ? op(qx,qy,0,e) : op(qx,qy,symmetric ? 2 : 3,e));
206 }
207 }
208
209 for (int dx = 0; dx < D1Dx; ++dx)
210 {
211 real_t val = 0.0;
212 for (int qx = 0; qx < Q1D; ++qx)
213 {
214 const real_t wx = (c == 0) ? Bc(qx,dx) : Bo(qx,dx);
215 val += mass[qx] * wx * wx;
216 }
217 diag(dx + (dy * D1Dx) + osc, e) += val;
218 }
219 }
220
221 osc += D1Dx * D1Dy;
222 } // loop (c) over components
223 }); // end of element loop
224}
225
226void PAHdivMassAssembleDiagonal3D(const int D1D,
227 const int Q1D,
228 const int NE,
229 const bool symmetric,
230 const Array<real_t> &Bo_,
231 const Array<real_t> &Bc_,
232 const Vector &op_,
233 Vector &diag_)
234{
235 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
236 "Error: D1D > HDIV_MAX_D1D");
237 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
238 "Error: Q1D > HDIV_MAX_Q1D");
239 constexpr static int VDIM = 3;
240
241 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
242 auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
243 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
244 auto diag = Reshape(diag_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
245
246 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
247 {
248 int osc = 0;
249
250 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
251 {
252 const int D1Dz = (c == 2) ? D1D : D1D - 1;
253 const int D1Dy = (c == 1) ? D1D : D1D - 1;
254 const int D1Dx = (c == 0) ? D1D : D1D - 1;
255
256 const int opc = (c == 0) ? 0 : ((c == 1) ? (symmetric ? 3 : 4) :
257 (symmetric ? 5 : 8));
258
259 real_t mass[DofQuadLimits::HDIV_MAX_Q1D];
260
261 for (int dz = 0; dz < D1Dz; ++dz)
262 {
263 for (int dy = 0; dy < D1Dy; ++dy)
264 {
265 for (int qx = 0; qx < Q1D; ++qx)
266 {
267 mass[qx] = 0.0;
268 for (int qy = 0; qy < Q1D; ++qy)
269 {
270 const real_t wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
271 for (int qz = 0; qz < Q1D; ++qz)
272 {
273 const real_t wz = (c == 2) ? Bc(qz,dz) : Bo(qz,dz);
274 mass[qx] += wy * wy * wz * wz * op(qx,qy,qz,opc,e);
275 }
276 }
277 }
278
279 for (int dx = 0; dx < D1Dx; ++dx)
280 {
281 real_t val = 0.0;
282 for (int qx = 0; qx < Q1D; ++qx)
283 {
284 const real_t wx = (c == 0) ? Bc(qx,dx) : Bo(qx,dx);
285 val += mass[qx] * wx * wx;
286 }
287 diag(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += val;
288 }
289 }
290 }
291
292 osc += D1Dx * D1Dy * D1Dz;
293 } // loop c
294 }); // end of element loop
295}
296
297void PAHdivMassApply(const int dim,
298 const int D1D,
299 const int Q1D,
300 const int NE,
301 const bool symmetric,
302 const Array<real_t> &Bo,
303 const Array<real_t> &Bc,
304 const Array<real_t> &Bot,
305 const Array<real_t> &Bct,
306 const Vector &op,
307 const Vector &x,
308 Vector &y)
309{
310 const int id = (D1D << 4) | Q1D;
311
312 if (dim == 2)
313 {
314 switch (id)
315 {
316 case 0x22: return SmemPAHdivMassApply2D<2,2>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
317 case 0x33: return SmemPAHdivMassApply2D<3,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
318 case 0x44: return SmemPAHdivMassApply2D<4,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
319 case 0x55: return SmemPAHdivMassApply2D<5,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
320 default: // fallback
321 return PAHdivMassApply2D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
322 }
323 }
324 else if (dim == 3)
325 {
326 switch (id)
327 {
328 case 0x23: return SmemPAHdivMassApply3D<2,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
329 case 0x34: return SmemPAHdivMassApply3D<3,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
330 case 0x45: return SmemPAHdivMassApply3D<4,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
331 case 0x56: return SmemPAHdivMassApply3D<5,6>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
332 case 0x67: return SmemPAHdivMassApply3D<6,7>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
333 case 0x78: return SmemPAHdivMassApply3D<7,8>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
334 default: // fallback
335 return PAHdivMassApply3D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
336 }
337 }
338}
339
340void PAHdivMassApply2D(const int D1D,
341 const int Q1D,
342 const int NE,
343 const bool symmetric,
344 const Array<real_t> &Bo_,
345 const Array<real_t> &Bc_,
346 const Array<real_t> &Bot_,
347 const Array<real_t> &Bct_,
348 const Vector &op_,
349 const Vector &x_,
350 Vector &y_)
351{
352 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
353 auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
354 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
355 auto Bct = Reshape(Bct_.Read(), D1D, Q1D);
356 auto op = Reshape(op_.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
357 auto x = Reshape(x_.Read(), 2*(D1D-1)*D1D, NE);
358 auto y = Reshape(y_.ReadWrite(), 2*(D1D-1)*D1D, NE);
359
360 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
361 {
362 constexpr static int VDIM = 2;
363 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
364 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
365
366 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
367
368 for (int qy = 0; qy < Q1D; ++qy)
369 {
370 for (int qx = 0; qx < Q1D; ++qx)
371 {
372 for (int c = 0; c < VDIM; ++c)
373 {
374 mass[qy][qx][c] = 0.0;
375 }
376 }
377 }
378
379 int osc = 0;
380
381 for (int c = 0; c < VDIM; ++c) // loop over x, y components
382 {
383 const int D1Dx = (c == 1) ? D1D - 1 : D1D;
384 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
385
386 for (int dy = 0; dy < D1Dy; ++dy)
387 {
388 real_t massX[MAX_Q1D];
389 for (int qx = 0; qx < Q1D; ++qx)
390 {
391 massX[qx] = 0.0;
392 }
393
394 for (int dx = 0; dx < D1Dx; ++dx)
395 {
396 const real_t t = x(dx + (dy * D1Dx) + osc, e);
397 for (int qx = 0; qx < Q1D; ++qx)
398 {
399 massX[qx] += t * ((c == 0) ? Bc(qx,dx) : Bo(qx,dx));
400 }
401 }
402
403 for (int qy = 0; qy < Q1D; ++qy)
404 {
405 const real_t wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
406 for (int qx = 0; qx < Q1D; ++qx)
407 {
408 mass[qy][qx][c] += massX[qx] * wy;
409 }
410 }
411 }
412
413 osc += D1Dx * D1Dy;
414 } // loop (c) over components
415
416 // Apply D operator.
417 for (int qy = 0; qy < Q1D; ++qy)
418 {
419 for (int qx = 0; qx < Q1D; ++qx)
420 {
421 const real_t O11 = op(qx,qy,0,e);
422 const real_t O12 = op(qx,qy,1,e);
423 const real_t O21 = symmetric ? O12 : op(qx,qy,2,e);
424 const real_t O22 = symmetric ? op(qx,qy,2,e) : op(qx,qy,3,e);
425 const real_t massX = mass[qy][qx][0];
426 const real_t massY = mass[qy][qx][1];
427 mass[qy][qx][0] = (O11*massX)+(O12*massY);
428 mass[qy][qx][1] = (O21*massX)+(O22*massY);
429 }
430 }
431
432 for (int qy = 0; qy < Q1D; ++qy)
433 {
434 osc = 0;
435
436 for (int c = 0; c < VDIM; ++c) // loop over x, y components
437 {
438 const int D1Dx = (c == 1) ? D1D - 1 : D1D;
439 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
440
441 real_t massX[MAX_D1D];
442 for (int dx = 0; dx < D1Dx; ++dx)
443 {
444 massX[dx] = 0;
445 }
446 for (int qx = 0; qx < Q1D; ++qx)
447 {
448 for (int dx = 0; dx < D1Dx; ++dx)
449 {
450 massX[dx] += mass[qy][qx][c] * ((c == 0) ? Bct(dx,qx) :
451 Bot(dx,qx));
452 }
453 }
454
455 for (int dy = 0; dy < D1Dy; ++dy)
456 {
457 const real_t wy = (c == 1) ? Bct(dy,qy) : Bot(dy,qy);
458
459 for (int dx = 0; dx < D1Dx; ++dx)
460 {
461 y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
462 }
463 }
464
465 osc += D1Dx * D1Dy;
466 } // loop c
467 } // loop qy
468 }); // end of element loop
469}
470
471void PAHdivMassApply3D(const int D1D,
472 const int Q1D,
473 const int NE,
474 const bool symmetric,
475 const Array<real_t> &Bo_,
476 const Array<real_t> &Bc_,
477 const Array<real_t> &Bot_,
478 const Array<real_t> &Bct_,
479 const Vector &op_,
480 const Vector &x_,
481 Vector &y_)
482{
483 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
484 "Error: D1D > HDIV_MAX_D1D");
485 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
486 "Error: Q1D > HDIV_MAX_Q1D");
487 constexpr static int VDIM = 3;
488
489 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
490 auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
491 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
492 auto Bct = Reshape(Bct_.Read(), D1D, Q1D);
493 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
494 auto x = Reshape(x_.Read(), 3*(D1D-1)*(D1D-1)*D1D, NE);
495 auto y = Reshape(y_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
496
497 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
498 {
499 real_t mass[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D][VDIM];
500
501 for (int qz = 0; qz < Q1D; ++qz)
502 {
503 for (int qy = 0; qy < Q1D; ++qy)
504 {
505 for (int qx = 0; qx < Q1D; ++qx)
506 {
507 for (int c = 0; c < VDIM; ++c)
508 {
509 mass[qz][qy][qx][c] = 0.0;
510 }
511 }
512 }
513 }
514
515 int osc = 0;
516
517 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
518 {
519 const int D1Dz = (c == 2) ? D1D : D1D - 1;
520 const int D1Dy = (c == 1) ? D1D : D1D - 1;
521 const int D1Dx = (c == 0) ? D1D : D1D - 1;
522
523 for (int dz = 0; dz < D1Dz; ++dz)
524 {
525 real_t massXY[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
526 for (int qy = 0; qy < Q1D; ++qy)
527 {
528 for (int qx = 0; qx < Q1D; ++qx)
529 {
530 massXY[qy][qx] = 0.0;
531 }
532 }
533
534 for (int dy = 0; dy < D1Dy; ++dy)
535 {
536 real_t massX[DofQuadLimits::HDIV_MAX_Q1D];
537 for (int qx = 0; qx < Q1D; ++qx)
538 {
539 massX[qx] = 0.0;
540 }
541
542 for (int dx = 0; dx < D1Dx; ++dx)
543 {
544 const real_t t = x(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
545 for (int qx = 0; qx < Q1D; ++qx)
546 {
547 massX[qx] += t * ((c == 0) ? Bc(qx,dx) : Bo(qx,dx));
548 }
549 }
550
551 for (int qy = 0; qy < Q1D; ++qy)
552 {
553 const real_t wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
554 for (int qx = 0; qx < Q1D; ++qx)
555 {
556 const real_t wx = massX[qx];
557 massXY[qy][qx] += wx * wy;
558 }
559 }
560 }
561
562 for (int qz = 0; qz < Q1D; ++qz)
563 {
564 const real_t wz = (c == 2) ? Bc(qz,dz) : Bo(qz,dz);
565 for (int qy = 0; qy < Q1D; ++qy)
566 {
567 for (int qx = 0; qx < Q1D; ++qx)
568 {
569 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
570 }
571 }
572 }
573 }
574
575 osc += D1Dx * D1Dy * D1Dz;
576 } // loop (c) over components
577
578 // Apply D operator.
579 for (int qz = 0; qz < Q1D; ++qz)
580 {
581 for (int qy = 0; qy < Q1D; ++qy)
582 {
583 for (int qx = 0; qx < Q1D; ++qx)
584 {
585 const real_t O11 = op(qx,qy,qz,0,e);
586 const real_t O12 = op(qx,qy,qz,1,e);
587 const real_t O13 = op(qx,qy,qz,2,e);
588 const real_t O21 = symmetric ? O12 : op(qx,qy,qz,3,e);
589 const real_t O22 = symmetric ? op(qx,qy,qz,3,e) : op(qx,qy,qz,4,e);
590 const real_t O23 = symmetric ? op(qx,qy,qz,4,e) : op(qx,qy,qz,5,e);
591 const real_t O31 = symmetric ? O13 : op(qx,qy,qz,6,e);
592 const real_t O32 = symmetric ? O23 : op(qx,qy,qz,7,e);
593 const real_t O33 = symmetric ? op(qx,qy,qz,5,e) : op(qx,qy,qz,8,e);
594
595 const real_t massX = mass[qz][qy][qx][0];
596 const real_t massY = mass[qz][qy][qx][1];
597 const real_t massZ = mass[qz][qy][qx][2];
598 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
599 mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
600 mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
601 }
602 }
603 }
604
605 for (int qz = 0; qz < Q1D; ++qz)
606 {
607 real_t massXY[DofQuadLimits::HDIV_MAX_D1D][DofQuadLimits::HDIV_MAX_D1D];
608
609 osc = 0;
610
611 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
612 {
613 const int D1Dz = (c == 2) ? D1D : D1D - 1;
614 const int D1Dy = (c == 1) ? D1D : D1D - 1;
615 const int D1Dx = (c == 0) ? D1D : D1D - 1;
616
617 for (int dy = 0; dy < D1Dy; ++dy)
618 {
619 for (int dx = 0; dx < D1Dx; ++dx)
620 {
621 massXY[dy][dx] = 0;
622 }
623 }
624 for (int qy = 0; qy < Q1D; ++qy)
625 {
626 real_t massX[DofQuadLimits::HDIV_MAX_D1D];
627 for (int dx = 0; dx < D1Dx; ++dx)
628 {
629 massX[dx] = 0;
630 }
631 for (int qx = 0; qx < Q1D; ++qx)
632 {
633 for (int dx = 0; dx < D1Dx; ++dx)
634 {
635 massX[dx] += mass[qz][qy][qx][c] *
636 ((c == 0) ? Bct(dx,qx) : Bot(dx,qx));
637 }
638 }
639 for (int dy = 0; dy < D1Dy; ++dy)
640 {
641 const real_t wy = (c == 1) ? Bct(dy,qy) : Bot(dy,qy);
642 for (int dx = 0; dx < D1Dx; ++dx)
643 {
644 massXY[dy][dx] += massX[dx] * wy;
645 }
646 }
647 }
648
649 for (int dz = 0; dz < D1Dz; ++dz)
650 {
651 const real_t wz = (c == 2) ? Bct(dz,qz) : Bot(dz,qz);
652 for (int dy = 0; dy < D1Dy; ++dy)
653 {
654 for (int dx = 0; dx < D1Dx; ++dx)
655 {
656 y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) +=
657 massXY[dy][dx] * wz;
658 }
659 }
660 }
661
662 osc += D1Dx * D1Dy * D1Dz;
663 } // loop c
664 } // loop qz
665 }); // end of element loop
666}
667
668// NOTE: this is identical to PACurlCurlSetup2D
669void PADivDivSetup2D(const int Q1D,
670 const int NE,
671 const Array<real_t> &w,
672 const Vector &j,
673 Vector &coeff_,
674 Vector &op)
675{
676 const int NQ = Q1D*Q1D;
677 auto W = w.Read();
678 auto J = Reshape(j.Read(), NQ, 2, 2, NE);
679 auto coeff = Reshape(coeff_.Read(), NQ, NE);
680 auto y = Reshape(op.Write(), NQ, NE);
681 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
682 {
683 for (int q = 0; q < NQ; ++q)
684 {
685 const real_t J11 = J(q,0,0,e);
686 const real_t J21 = J(q,1,0,e);
687 const real_t J12 = J(q,0,1,e);
688 const real_t J22 = J(q,1,1,e);
689 const real_t detJ = (J11*J22)-(J21*J12);
690 y(q,e) = W[q] * coeff(q,e) / detJ;
691 }
692 });
693}
694
695void PADivDivSetup3D(const int Q1D,
696 const int NE,
697 const Array<real_t> &w,
698 const Vector &j,
699 Vector &coeff_,
700 Vector &op)
701{
702 const int NQ = Q1D*Q1D*Q1D;
703 auto W = w.Read();
704 auto J = Reshape(j.Read(), NQ, 3, 3, NE);
705 auto coeff = Reshape(coeff_.Read(), NQ, NE);
706 auto y = Reshape(op.Write(), NQ, NE);
707
708 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
709 {
710 for (int q = 0; q < NQ; ++q)
711 {
712 const real_t J11 = J(q,0,0,e);
713 const real_t J21 = J(q,1,0,e);
714 const real_t J31 = J(q,2,0,e);
715 const real_t J12 = J(q,0,1,e);
716 const real_t J22 = J(q,1,1,e);
717 const real_t J32 = J(q,2,1,e);
718 const real_t J13 = J(q,0,2,e);
719 const real_t J23 = J(q,1,2,e);
720 const real_t J33 = J(q,2,2,e);
721 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
722 J21 * (J12 * J33 - J32 * J13) +
723 J31 * (J12 * J23 - J22 * J13);
724 y(q,e) = W[q] * coeff(q, e) / detJ;
725 }
726 });
727}
728
729void PADivDivAssembleDiagonal2D(const int D1D,
730 const int Q1D,
731 const int NE,
732 const Array<real_t> &Bo_,
733 const Array<real_t> &Gc_,
734 const Vector &op_,
735 Vector &diag_)
736{
737 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
738 auto Gc = Reshape(Gc_.Read(), Q1D, D1D);
739 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
740 auto diag = Reshape(diag_.ReadWrite(), 2*(D1D-1)*D1D, NE);
741
742 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
743 {
744 constexpr static int VDIM = 2;
745 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
746
747 int osc = 0;
748
749 for (int c = 0; c < VDIM; ++c) // loop over x, y components
750 {
751 const int D1Dx = (c == 1) ? D1D - 1 : D1D;
752 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
753
754 real_t div[MAX_Q1D];
755
756 for (int dy = 0; dy < D1Dy; ++dy)
757 {
758 for (int qx = 0; qx < Q1D; ++qx)
759 {
760 div[qx] = 0.0;
761 for (int qy = 0; qy < Q1D; ++qy)
762 {
763 const real_t wy = (c == 0) ? Bo(qy,dy) : Gc(qy,dy);
764 div[qx] += wy * wy * op(qx,qy,e);
765 }
766 }
767
768 for (int dx = 0; dx < D1Dx; ++dx)
769 {
770 real_t val = 0.0;
771 for (int qx = 0; qx < Q1D; ++qx)
772 {
773 const real_t wx = (c == 0) ? Gc(qx,dx) : Bo(qx,dx);
774 val += div[qx] * wx * wx;
775 }
776 diag(dx + (dy * D1Dx) + osc, e) += val;
777 }
778 }
779
780 osc += D1Dx * D1Dy;
781 } // loop c
782 });
783}
784
785void PADivDivAssembleDiagonal3D(const int D1D,
786 const int Q1D,
787 const int NE,
788 const Array<real_t> &Bo_,
789 const Array<real_t> &Gc_,
790 const Vector &op_,
791 Vector &diag_)
792{
793 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
794 "Error: D1D > HDIV_MAX_D1D");
795 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
796 "Error: Q1D > HDIV_MAX_Q1D");
797 constexpr static int VDIM = 3;
798
799 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
800 auto Gc = Reshape(Gc_.Read(), Q1D, D1D);
801 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
802 auto diag = Reshape(diag_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
803
804 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
805 {
806 int osc = 0;
807
808 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
809 {
810 const int D1Dz = (c == 2) ? D1D : D1D - 1;
811 const int D1Dy = (c == 1) ? D1D : D1D - 1;
812 const int D1Dx = (c == 0) ? D1D : D1D - 1;
813
814 for (int dz = 0; dz < D1Dz; ++dz)
815 {
816 for (int dy = 0; dy < D1Dy; ++dy)
817 {
818 real_t a[DofQuadLimits::HDIV_MAX_Q1D];
819
820 for (int qx = 0; qx < Q1D; ++qx)
821 {
822 a[qx] = 0.0;
823 for (int qy = 0; qy < Q1D; ++qy)
824 {
825 const real_t wy = (c == 1) ? Gc(qy,dy) : Bo(qy,dy);
826
827 for (int qz = 0; qz < Q1D; ++qz)
828 {
829 const real_t wz = (c == 2) ? Gc(qz,dz) : Bo(qz,dz);
830 a[qx] += wy * wy * wz * wz * op(qx,qy,qz,e);
831 }
832 }
833 }
834
835 for (int dx = 0; dx < D1Dx; ++dx)
836 {
837 real_t val = 0.0;
838 for (int qx = 0; qx < Q1D; ++qx)
839 {
840 const real_t wx = (c == 0) ? Gc(qx,dx) : Bo(qx,dx);
841 val += a[qx] * wx * wx;
842 }
843 diag(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += val;
844 }
845 }
846 }
847
848 osc += D1Dx * D1Dy * D1Dz;
849 } // loop c
850 }); // end of element loop
851}
852
853void PADivDivApply2D(const int D1D,
854 const int Q1D,
855 const int NE,
856 const Array<real_t> &Bo_,
857 const Array<real_t> &Gc_,
858 const Array<real_t> &Bot_,
859 const Array<real_t> &Gct_,
860 const Vector &op_,
861 const Vector &x_,
862 Vector &y_)
863{
864 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
865 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
866 auto Gc = Reshape(Gc_.Read(), Q1D, D1D);
867 auto Gct = Reshape(Gct_.Read(), D1D, Q1D);
868 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
869 auto x = Reshape(x_.Read(), 2*(D1D-1)*D1D, NE);
870 auto y = Reshape(y_.ReadWrite(), 2*(D1D-1)*D1D, NE);
871
872 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
873 {
874 constexpr static int VDIM = 2;
875 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
876 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
877
878 real_t div[MAX_Q1D][MAX_Q1D];
879
880 // div[qy][qx] will be computed as du_x/dx + du_y/dy
881
882 for (int qy = 0; qy < Q1D; ++qy)
883 {
884 for (int qx = 0; qx < Q1D; ++qx)
885 {
886 div[qy][qx] = 0;
887 }
888 }
889
890 int osc = 0;
891
892 for (int c = 0; c < VDIM; ++c) // loop over x, y components
893 {
894 const int D1Dx = (c == 1) ? D1D - 1 : D1D;
895 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
896
897 for (int dy = 0; dy < D1Dy; ++dy)
898 {
899 real_t gradX[MAX_Q1D];
900 for (int qx = 0; qx < Q1D; ++qx)
901 {
902 gradX[qx] = 0;
903 }
904
905 for (int dx = 0; dx < D1Dx; ++dx)
906 {
907 const real_t t = x(dx + (dy * D1Dx) + osc, e);
908 for (int qx = 0; qx < Q1D; ++qx)
909 {
910 gradX[qx] += t * ((c == 0) ? Gc(qx,dx) : Bo(qx,dx));
911 }
912 }
913
914 for (int qy = 0; qy < Q1D; ++qy)
915 {
916 const real_t wy = (c == 0) ? Bo(qy,dy) : Gc(qy,dy);
917 for (int qx = 0; qx < Q1D; ++qx)
918 {
919 div[qy][qx] += gradX[qx] * wy;
920 }
921 }
922 }
923
924 osc += D1Dx * D1Dy;
925 } // loop (c) over components
926
927 // Apply D operator.
928 for (int qy = 0; qy < Q1D; ++qy)
929 {
930 for (int qx = 0; qx < Q1D; ++qx)
931 {
932 div[qy][qx] *= op(qx,qy,e);
933 }
934 }
935
936 for (int qy = 0; qy < Q1D; ++qy)
937 {
938 osc = 0;
939
940 for (int c = 0; c < VDIM; ++c) // loop over x, y components
941 {
942 const int D1Dx = (c == 1) ? D1D - 1 : D1D;
943 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
944
945 real_t gradX[MAX_D1D];
946 for (int dx = 0; dx < D1Dx; ++dx)
947 {
948 gradX[dx] = 0;
949 }
950 for (int qx = 0; qx < Q1D; ++qx)
951 {
952 for (int dx = 0; dx < D1Dx; ++dx)
953 {
954 gradX[dx] += div[qy][qx] * (c == 0 ? Gct(dx,qx) : Bot(dx,qx));
955 }
956 }
957 for (int dy = 0; dy < D1Dy; ++dy)
958 {
959 const real_t wy = (c == 0) ? Bot(dy,qy) : Gct(dy,qy);
960 for (int dx = 0; dx < D1Dx; ++dx)
961 {
962 y(dx + (dy * D1Dx) + osc, e) += gradX[dx] * wy;
963 }
964 }
965
966 osc += D1Dx * D1Dy;
967 } // loop c
968 } // loop qy
969 }); // end of element loop
970}
971
972void PADivDivApply3D(const int D1D,
973 const int Q1D,
974 const int NE,
975 const Array<real_t> &Bo_,
976 const Array<real_t> &Gc_,
977 const Array<real_t> &Bot_,
978 const Array<real_t> &Gct_,
979 const Vector &op_,
980 const Vector &x_,
981 Vector &y_)
982{
983 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
984 "Error: D1D > HDIV_MAX_D1D");
985 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
986 "Error: Q1D > HDIV_MAX_Q1D");
987 constexpr static int VDIM = 3;
988
989 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
990 auto Gc = Reshape(Gc_.Read(), Q1D, D1D);
991 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
992 auto Gct = Reshape(Gct_.Read(), D1D, Q1D);
993 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
994 auto x = Reshape(x_.Read(), 3*(D1D-1)*(D1D-1)*D1D, NE);
995 auto y = Reshape(y_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
996
997 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
998 {
999 real_t div[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1000
1001 for (int qz = 0; qz < Q1D; ++qz)
1002 {
1003 for (int qy = 0; qy < Q1D; ++qy)
1004 {
1005 for (int qx = 0; qx < Q1D; ++qx)
1006 {
1007 div[qz][qy][qx] = 0.0;
1008 }
1009 }
1010 }
1011
1012 int osc = 0;
1013
1014 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1015 {
1016 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1017 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1018 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1019
1020 for (int dz = 0; dz < D1Dz; ++dz)
1021 {
1022 real_t aXY[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1023 for (int qy = 0; qy < Q1D; ++qy)
1024 {
1025 for (int qx = 0; qx < Q1D; ++qx)
1026 {
1027 aXY[qy][qx] = 0.0;
1028 }
1029 }
1030
1031 for (int dy = 0; dy < D1Dy; ++dy)
1032 {
1033 real_t aX[DofQuadLimits::HDIV_MAX_Q1D];
1034 for (int qx = 0; qx < Q1D; ++qx)
1035 {
1036 aX[qx] = 0.0;
1037 }
1038
1039 for (int dx = 0; dx < D1Dx; ++dx)
1040 {
1041 const real_t t = x(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1042 for (int qx = 0; qx < Q1D; ++qx)
1043 {
1044 aX[qx] += t * ((c == 0) ? Gc(qx,dx) : Bo(qx,dx));
1045 }
1046 }
1047
1048 for (int qy = 0; qy < Q1D; ++qy)
1049 {
1050 const real_t wy = (c == 1) ? Gc(qy,dy) : Bo(qy,dy);
1051 for (int qx = 0; qx < Q1D; ++qx)
1052 {
1053 const real_t wx = aX[qx];
1054 aXY[qy][qx] += wx * wy;
1055 }
1056 }
1057 }
1058
1059 for (int qz = 0; qz < Q1D; ++qz)
1060 {
1061 const real_t wz = (c == 2) ? Gc(qz,dz) : Bo(qz,dz);
1062 for (int qy = 0; qy < Q1D; ++qy)
1063 {
1064 for (int qx = 0; qx < Q1D; ++qx)
1065 {
1066 div[qz][qy][qx] += aXY[qy][qx] * wz;
1067 }
1068 }
1069 }
1070 }
1071
1072 osc += D1Dx * D1Dy * D1Dz;
1073 } // loop (c) over components
1074
1075 // Apply D operator.
1076 for (int qz = 0; qz < Q1D; ++qz)
1077 {
1078 for (int qy = 0; qy < Q1D; ++qy)
1079 {
1080 for (int qx = 0; qx < Q1D; ++qx)
1081 {
1082 div[qz][qy][qx] *= op(qx,qy,qz,e);
1083 }
1084 }
1085 }
1086
1087 for (int qz = 0; qz < Q1D; ++qz)
1088 {
1089 real_t aXY[DofQuadLimits::HDIV_MAX_D1D][DofQuadLimits::HDIV_MAX_D1D];
1090
1091 osc = 0;
1092
1093 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1094 {
1095 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1096 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1097 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1098
1099 for (int dy = 0; dy < D1Dy; ++dy)
1100 {
1101 for (int dx = 0; dx < D1Dx; ++dx)
1102 {
1103 aXY[dy][dx] = 0;
1104 }
1105 }
1106 for (int qy = 0; qy < Q1D; ++qy)
1107 {
1108 real_t aX[DofQuadLimits::HDIV_MAX_D1D];
1109 for (int dx = 0; dx < D1Dx; ++dx)
1110 {
1111 aX[dx] = 0;
1112 }
1113 for (int qx = 0; qx < Q1D; ++qx)
1114 {
1115 for (int dx = 0; dx < D1Dx; ++dx)
1116 {
1117 aX[dx] += div[qz][qy][qx] *
1118 (c == 0 ? Gct(dx,qx) : Bot(dx,qx));
1119 }
1120 }
1121 for (int dy = 0; dy < D1Dy; ++dy)
1122 {
1123 const real_t wy = (c == 1) ? Gct(dy,qy) : Bot(dy,qy);
1124 for (int dx = 0; dx < D1Dx; ++dx)
1125 {
1126 aXY[dy][dx] += aX[dx] * wy;
1127 }
1128 }
1129 }
1130
1131 for (int dz = 0; dz < D1Dz; ++dz)
1132 {
1133 const real_t wz = (c == 2) ? Gct(dz,qz) : Bot(dz,qz);
1134 for (int dy = 0; dy < D1Dy; ++dy)
1135 {
1136 for (int dx = 0; dx < D1Dx; ++dx)
1137 {
1138 y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) +=
1139 aXY[dy][dx] * wz;
1140 }
1141 }
1142 }
1143
1144 osc += D1Dx * D1Dy * D1Dz;
1145 } // loop c
1146 } // loop qz
1147 }); // end of element loop
1148}
1149
1150void PAHdivL2Setup2D(const int Q1D,
1151 const int NE,
1152 const Array<real_t> &w,
1153 Vector &coeff_,
1154 Vector &op)
1155{
1156 const int NQ = Q1D*Q1D;
1157 auto W = w.Read();
1158 auto coeff = Reshape(coeff_.Read(), NQ, NE);
1159 auto y = Reshape(op.Write(), NQ, NE);
1160 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1161 {
1162 for (int q = 0; q < NQ; ++q)
1163 {
1164 y(q,e) = W[q] * coeff(q,e);
1165 }
1166 });
1167}
1168
1169void PAHdivL2Setup3D(const int Q1D,
1170 const int NE,
1171 const Array<real_t> &w,
1172 Vector &coeff_,
1173 Vector &op)
1174{
1175 const int NQ = Q1D*Q1D*Q1D;
1176 auto W = w.Read();
1177 auto coeff = Reshape(coeff_.Read(), NQ, NE);
1178 auto y = Reshape(op.Write(), NQ, NE);
1179
1180 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1181 {
1182 for (int q = 0; q < NQ; ++q)
1183 {
1184 y(q,e) = W[q] * coeff(q, e);
1185 }
1186 });
1187}
1188
1189void PAHdivL2AssembleDiagonal_ADAt_2D(const int D1D,
1190 const int Q1D,
1191 const int L2D1D,
1192 const int NE,
1193 const Array<real_t> &L2Bo_,
1194 const Array<real_t> &Gct_,
1195 const Array<real_t> &Bot_,
1196 const Vector &op_,
1197 const Vector &D_,
1198 Vector &diag_)
1199{
1200 constexpr static int VDIM = 2;
1201
1202 auto L2Bo = Reshape(L2Bo_.Read(), Q1D, L2D1D);
1203 auto Gct = Reshape(Gct_.Read(), D1D, Q1D);
1204 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
1205 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
1206 auto D = Reshape(D_.Read(), 2*(D1D-1)*D1D, NE);
1207 auto diag = Reshape(diag_.ReadWrite(), L2D1D, L2D1D, NE);
1208
1209 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1210 {
1211 for (int ry = 0; ry < L2D1D; ++ry)
1212 {
1213 for (int rx = 0; rx < L2D1D; ++rx)
1214 {
1215 // Compute row (rx,ry), assuming all contributions are from
1216 // a single element.
1217
1218 real_t row[2*DofQuadLimits::HDIV_MAX_D1D*(DofQuadLimits::HDIV_MAX_D1D-1)];
1219 real_t div[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1220
1221 for (int i=0; i<2*D1D*(D1D - 1); ++i)
1222 {
1223 row[i] = 0;
1224 }
1225
1226 for (int qy = 0; qy < Q1D; ++qy)
1227 {
1228 for (int qx = 0; qx < Q1D; ++qx)
1229 {
1230 div[qy][qx] = op(qx,qy,e) * L2Bo(qx,rx) * L2Bo(qy,ry);
1231 }
1232 }
1233
1234 for (int qy = 0; qy < Q1D; ++qy)
1235 {
1236 int osc = 0;
1237 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1238 {
1239 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1240 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1241
1242 real_t aX[DofQuadLimits::HDIV_MAX_D1D];
1243 for (int dx = 0; dx < D1Dx; ++dx)
1244 {
1245 aX[dx] = 0;
1246 }
1247 for (int qx = 0; qx < Q1D; ++qx)
1248 {
1249 for (int dx = 0; dx < D1Dx; ++dx)
1250 {
1251 aX[dx] += div[qy][qx] * ((c == 0) ? Gct(dx,qx) :
1252 Bot(dx,qx));
1253 }
1254 }
1255
1256 for (int dy = 0; dy < D1Dy; ++dy)
1257 {
1258 const real_t wy = (c == 1) ? Gct(dy,qy) : Bot(dy,qy);
1259
1260 for (int dx = 0; dx < D1Dx; ++dx)
1261 {
1262 row[dx + (dy * D1Dx) + osc] += aX[dx] * wy;
1263 }
1264 }
1265
1266 osc += D1Dx * D1Dy;
1267 } // loop c
1268 } // loop qy
1269
1270 real_t val = 0.0;
1271 for (int i=0; i<2*D1D*(D1D - 1); ++i)
1272 {
1273 val += row[i] * row[i] * D(i,e);
1274 }
1275 diag(rx,ry,e) += val;
1276 } // loop rx
1277 } // loop ry
1278 }); // end of element loop
1279}
1280
1281void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
1282 const int Q1D,
1283 const int L2D1D,
1284 const int NE,
1285 const Array<real_t> &L2Bo_,
1286 const Array<real_t> &Gct_,
1287 const Array<real_t> &Bot_,
1288 const Vector &op_,
1289 const Vector &D_,
1290 Vector &diag_)
1291{
1292 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
1293 "Error: D1D > HDIV_MAX_D1D");
1294 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
1295 "Error: Q1D > HDIV_MAX_Q1D");
1296 constexpr static int VDIM = 3;
1297
1298 auto L2Bo = Reshape(L2Bo_.Read(), Q1D, L2D1D);
1299 auto Gct = Reshape(Gct_.Read(), D1D, Q1D);
1300 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
1301 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
1302 auto D = Reshape(D_.Read(), 3*(D1D-1)*(D1D-1)*D1D, NE);
1303 auto diag = Reshape(diag_.ReadWrite(), L2D1D, L2D1D, L2D1D, NE);
1304
1305 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1306 {
1307 for (int rz = 0; rz < L2D1D; ++rz)
1308 {
1309 for (int ry = 0; ry < L2D1D; ++ry)
1310 {
1311 for (int rx = 0; rx < L2D1D; ++rx)
1312 {
1313 // Compute row (rx,ry,rz), assuming all contributions are from
1314 // a single element.
1315
1316 real_t row[3*DofQuadLimits::HDIV_MAX_D1D*(DofQuadLimits::HDIV_MAX_D1D-1)*
1317 (DofQuadLimits::HDIV_MAX_D1D-1)];
1318 real_t div[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1319
1320 for (int i=0; i<3*D1D*(D1D - 1)*(D1D - 1); ++i)
1321 {
1322 row[i] = 0;
1323 }
1324
1325 for (int qz = 0; qz < Q1D; ++qz)
1326 {
1327 for (int qy = 0; qy < Q1D; ++qy)
1328 {
1329 for (int qx = 0; qx < Q1D; ++qx)
1330 {
1331 div[qz][qy][qx] = op(qx,qy,qz,e) * L2Bo(qx,rx) *
1332 L2Bo(qy,ry) * L2Bo(qz,rz);
1333 }
1334 }
1335 }
1336
1337 for (int qz = 0; qz < Q1D; ++qz)
1338 {
1339 real_t aXY[DofQuadLimits::HDIV_MAX_D1D][DofQuadLimits::HDIV_MAX_D1D];
1340
1341 int osc = 0;
1342 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1343 {
1344 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1345 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1346 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1347
1348 for (int dy = 0; dy < D1Dy; ++dy)
1349 {
1350 for (int dx = 0; dx < D1Dx; ++dx)
1351 {
1352 aXY[dy][dx] = 0;
1353 }
1354 }
1355 for (int qy = 0; qy < Q1D; ++qy)
1356 {
1357 real_t aX[DofQuadLimits::HDIV_MAX_D1D];
1358 for (int dx = 0; dx < D1Dx; ++dx)
1359 {
1360 aX[dx] = 0;
1361 }
1362 for (int qx = 0; qx < Q1D; ++qx)
1363 {
1364 for (int dx = 0; dx < D1Dx; ++dx)
1365 {
1366 aX[dx] += div[qz][qy][qx] * ((c == 0) ? Gct(dx,qx)
1367 : Bot(dx,qx));
1368 }
1369 }
1370 for (int dy = 0; dy < D1Dy; ++dy)
1371 {
1372 const real_t wy = (c == 1) ? Gct(dy,qy) : Bot(dy,qy);
1373 for (int dx = 0; dx < D1Dx; ++dx)
1374 {
1375 aXY[dy][dx] += aX[dx] * wy;
1376 }
1377 }
1378 }
1379
1380 for (int dz = 0; dz < D1Dz; ++dz)
1381 {
1382 const real_t wz = (c == 2) ? Gct(dz,qz) : Bot(dz,qz);
1383 for (int dy = 0; dy < D1Dy; ++dy)
1384 {
1385 for (int dx = 0; dx < D1Dx; ++dx)
1386 {
1387 row[dx + ((dy + (dz * D1Dy)) * D1Dx) + osc] +=
1388 aXY[dy][dx] * wz;
1389 }
1390 }
1391 }
1392
1393 osc += D1Dx * D1Dy * D1Dz;
1394 } // loop c
1395 } // loop qz
1396
1397 real_t val = 0.0;
1398 for (int i=0; i<3*D1D*(D1D - 1)*(D1D - 1); ++i)
1399 {
1400 val += row[i] * row[i] * D(i,e);
1401 }
1402 diag(rx,ry,rz,e) += val;
1403 } // loop rx
1404 } // loop ry
1405 } // loop rz
1406 }); // end of element loop
1407}
1408
1409// Apply to x corresponding to DOFs in H(div) (trial), whose divergence is
1410// integrated against L_2 test functions corresponding to y.
1411void PAHdivL2Apply2D(const int D1D,
1412 const int Q1D,
1413 const int L2D1D,
1414 const int NE,
1415 const Array<real_t> &Bo_,
1416 const Array<real_t> &Gc_,
1417 const Array<real_t> &L2Bot_,
1418 const Vector &op_,
1419 const Vector &x_,
1420 Vector &y_)
1421{
1422 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
1423 auto Gc = Reshape(Gc_.Read(), Q1D, D1D);
1424 auto L2Bot = Reshape(L2Bot_.Read(), L2D1D, Q1D);
1425 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
1426 auto x = Reshape(x_.Read(), 2*(D1D-1)*D1D, NE);
1427 auto y = Reshape(y_.ReadWrite(), L2D1D, L2D1D, NE);
1428
1429 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1430 {
1431 constexpr static int VDIM = 2;
1432 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
1433 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
1434
1435 real_t div[MAX_Q1D][MAX_Q1D];
1436
1437 for (int qy = 0; qy < Q1D; ++qy)
1438 {
1439 for (int qx = 0; qx < Q1D; ++qx)
1440 {
1441 div[qy][qx] = 0.0;
1442 }
1443 }
1444
1445 int osc = 0;
1446
1447 for (int c = 0; c < VDIM; ++c) // loop over x, y components
1448 {
1449 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1450 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1451
1452 for (int dy = 0; dy < D1Dy; ++dy)
1453 {
1454 real_t aX[MAX_Q1D];
1455 for (int qx = 0; qx < Q1D; ++qx)
1456 {
1457 aX[qx] = 0.0;
1458 }
1459
1460 for (int dx = 0; dx < D1Dx; ++dx)
1461 {
1462 const real_t t = x(dx + (dy * D1Dx) + osc, e);
1463 for (int qx = 0; qx < Q1D; ++qx)
1464 {
1465 aX[qx] += t * ((c == 0) ? Gc(qx,dx) : Bo(qx,dx));
1466 }
1467 }
1468
1469 for (int qy = 0; qy < Q1D; ++qy)
1470 {
1471 const real_t wy = (c == 1) ? Gc(qy,dy) : Bo(qy,dy);
1472 for (int qx = 0; qx < Q1D; ++qx)
1473 {
1474 div[qy][qx] += aX[qx] * wy;
1475 }
1476 }
1477 }
1478
1479 osc += D1Dx * D1Dy;
1480 } // loop (c) over components
1481
1482 // Apply D operator.
1483 for (int qy = 0; qy < Q1D; ++qy)
1484 {
1485 for (int qx = 0; qx < Q1D; ++qx)
1486 {
1487 div[qy][qx] *= op(qx,qy,e);
1488 }
1489 }
1490
1491 for (int qy = 0; qy < Q1D; ++qy)
1492 {
1493 real_t aX[MAX_D1D];
1494 for (int dx = 0; dx < L2D1D; ++dx)
1495 {
1496 aX[dx] = 0;
1497 }
1498 for (int qx = 0; qx < Q1D; ++qx)
1499 {
1500 for (int dx = 0; dx < L2D1D; ++dx)
1501 {
1502 aX[dx] += div[qy][qx] * L2Bot(dx,qx);
1503 }
1504 }
1505 for (int dy = 0; dy < L2D1D; ++dy)
1506 {
1507 const real_t wy = L2Bot(dy,qy);
1508 for (int dx = 0; dx < L2D1D; ++dx)
1509 {
1510 y(dx,dy,e) += aX[dx] * wy;
1511 }
1512 }
1513 }
1514 }); // end of element loop
1515}
1516
1517void PAHdivL2ApplyTranspose2D(const int D1D,
1518 const int Q1D,
1519 const int L2D1D,
1520 const int NE,
1521 const Array<real_t> &L2Bo_,
1522 const Array<real_t> &Gct_,
1523 const Array<real_t> &Bot_,
1524 const Vector &op_,
1525 const Vector &x_,
1526 Vector &y_)
1527{
1528 auto L2Bo = Reshape(L2Bo_.Read(), Q1D, L2D1D);
1529 auto Gct = Reshape(Gct_.Read(), D1D, Q1D);
1530 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
1531 auto op = Reshape(op_.Read(), Q1D, Q1D, NE);
1532 auto x = Reshape(x_.Read(), L2D1D, L2D1D, NE);
1533 auto y = Reshape(y_.ReadWrite(), 2*(D1D-1)*D1D, NE);
1534
1535 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1536 {
1537 constexpr static int VDIM = 2;
1538 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
1539 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
1540
1541 real_t div[MAX_Q1D][MAX_Q1D];
1542
1543 for (int qy = 0; qy < Q1D; ++qy)
1544 {
1545 for (int qx = 0; qx < Q1D; ++qx)
1546 {
1547 div[qy][qx] = 0.0;
1548 }
1549 }
1550
1551 for (int dy = 0; dy < L2D1D; ++dy)
1552 {
1553 real_t aX[MAX_Q1D];
1554 for (int qx = 0; qx < Q1D; ++qx)
1555 {
1556 aX[qx] = 0.0;
1557 }
1558
1559 for (int dx = 0; dx < L2D1D; ++dx)
1560 {
1561 const real_t t = x(dx,dy,e);
1562 for (int qx = 0; qx < Q1D; ++qx)
1563 {
1564 aX[qx] += t * L2Bo(qx,dx);
1565 }
1566 }
1567
1568 for (int qy = 0; qy < Q1D; ++qy)
1569 {
1570 const real_t wy = L2Bo(qy,dy);
1571 for (int qx = 0; qx < Q1D; ++qx)
1572 {
1573 div[qy][qx] += aX[qx] * wy;
1574 }
1575 }
1576 }
1577
1578 // Apply D operator.
1579 for (int qy = 0; qy < Q1D; ++qy)
1580 {
1581 for (int qx = 0; qx < Q1D; ++qx)
1582 {
1583 div[qy][qx] *= op(qx,qy,e);
1584 }
1585 }
1586
1587 for (int qy = 0; qy < Q1D; ++qy)
1588 {
1589 real_t aX[MAX_D1D];
1590
1591 int osc = 0;
1592 for (int c = 0; c < VDIM; ++c) // loop over x, y components
1593 {
1594 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1595 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1596
1597 for (int dx = 0; dx < D1Dx; ++dx)
1598 {
1599 aX[dx] = 0;
1600 }
1601 for (int qx = 0; qx < Q1D; ++qx)
1602 {
1603 for (int dx = 0; dx < D1Dx; ++dx)
1604 {
1605 aX[dx] += div[qy][qx] * ((c == 0) ? Gct(dx,qx) : Bot(dx,qx));
1606 }
1607 }
1608 for (int dy = 0; dy < D1Dy; ++dy)
1609 {
1610 const real_t wy = (c == 0) ? Bot(dy,qy) : Gct(dy,qy);
1611 for (int dx = 0; dx < D1Dx; ++dx)
1612 {
1613 y(dx + (dy * D1Dx) + osc, e) += aX[dx] * wy;
1614 }
1615 }
1616
1617 osc += D1Dx * D1Dy;
1618 } // loop c
1619 } // loop qy
1620 }); // end of element loop
1621}
1622
1623// Apply to x corresponding to DOFs in H(div) (trial), whose divergence is
1624// integrated against L_2 test functions corresponding to y.
1625void PAHdivL2Apply3D(const int D1D,
1626 const int Q1D,
1627 const int L2D1D,
1628 const int NE,
1629 const Array<real_t> &Bo_,
1630 const Array<real_t> &Gc_,
1631 const Array<real_t> &L2Bot_,
1632 const Vector &op_,
1633 const Vector &x_,
1634 Vector &y_)
1635{
1636 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
1637 "Error: D1D > HDIV_MAX_D1D");
1638 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
1639 "Error: Q1D > HDIV_MAX_Q1D");
1640 constexpr static int VDIM = 3;
1641
1642 auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
1643 auto Gc = Reshape(Gc_.Read(), Q1D, D1D);
1644 auto L2Bot = Reshape(L2Bot_.Read(), L2D1D, Q1D);
1645 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
1646 auto x = Reshape(x_.Read(), 3*(D1D-1)*(D1D-1)*D1D, NE);
1647 auto y = Reshape(y_.ReadWrite(), L2D1D, L2D1D, L2D1D, NE);
1648
1649 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1650 {
1651 real_t div[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1652
1653 for (int qz = 0; qz < Q1D; ++qz)
1654 {
1655 for (int qy = 0; qy < Q1D; ++qy)
1656 {
1657 for (int qx = 0; qx < Q1D; ++qx)
1658 {
1659 div[qz][qy][qx] = 0.0;
1660 }
1661 }
1662 }
1663
1664 int osc = 0;
1665
1666 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1667 {
1668 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1669 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1670 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1671
1672 for (int dz = 0; dz < D1Dz; ++dz)
1673 {
1674 real_t aXY[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1675 for (int qy = 0; qy < Q1D; ++qy)
1676 {
1677 for (int qx = 0; qx < Q1D; ++qx)
1678 {
1679 aXY[qy][qx] = 0.0;
1680 }
1681 }
1682
1683 for (int dy = 0; dy < D1Dy; ++dy)
1684 {
1685 real_t aX[DofQuadLimits::HDIV_MAX_Q1D];
1686 for (int qx = 0; qx < Q1D; ++qx)
1687 {
1688 aX[qx] = 0.0;
1689 }
1690
1691 for (int dx = 0; dx < D1Dx; ++dx)
1692 {
1693 const real_t t = x(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1694 for (int qx = 0; qx < Q1D; ++qx)
1695 {
1696 aX[qx] += t * ((c == 0) ? Gc(qx,dx) : Bo(qx,dx));
1697 }
1698 }
1699
1700 for (int qy = 0; qy < Q1D; ++qy)
1701 {
1702 const real_t wy = (c == 1) ? Gc(qy,dy) : Bo(qy,dy);
1703 for (int qx = 0; qx < Q1D; ++qx)
1704 {
1705 aXY[qy][qx] += aX[qx] * wy;
1706 }
1707 }
1708 }
1709
1710 for (int qz = 0; qz < Q1D; ++qz)
1711 {
1712 const real_t wz = (c == 2) ? Gc(qz,dz) : Bo(qz,dz);
1713 for (int qy = 0; qy < Q1D; ++qy)
1714 {
1715 for (int qx = 0; qx < Q1D; ++qx)
1716 {
1717 div[qz][qy][qx] += aXY[qy][qx] * wz;
1718 }
1719 }
1720 }
1721 }
1722
1723 osc += D1Dx * D1Dy * D1Dz;
1724 } // loop (c) over components
1725
1726 // Apply D operator.
1727 for (int qz = 0; qz < Q1D; ++qz)
1728 {
1729 for (int qy = 0; qy < Q1D; ++qy)
1730 {
1731 for (int qx = 0; qx < Q1D; ++qx)
1732 {
1733 div[qz][qy][qx] *= op(qx,qy,qz,e);
1734 }
1735 }
1736 }
1737
1738 for (int qz = 0; qz < Q1D; ++qz)
1739 {
1740 real_t aXY[DofQuadLimits::HDIV_MAX_D1D][DofQuadLimits::HDIV_MAX_D1D];
1741
1742 for (int dy = 0; dy < L2D1D; ++dy)
1743 {
1744 for (int dx = 0; dx < L2D1D; ++dx)
1745 {
1746 aXY[dy][dx] = 0;
1747 }
1748 }
1749 for (int qy = 0; qy < Q1D; ++qy)
1750 {
1751 real_t aX[DofQuadLimits::HDIV_MAX_D1D];
1752 for (int dx = 0; dx < L2D1D; ++dx)
1753 {
1754 aX[dx] = 0;
1755 }
1756 for (int qx = 0; qx < Q1D; ++qx)
1757 {
1758 for (int dx = 0; dx < L2D1D; ++dx)
1759 {
1760 aX[dx] += div[qz][qy][qx] * L2Bot(dx,qx);
1761 }
1762 }
1763 for (int dy = 0; dy < L2D1D; ++dy)
1764 {
1765 const real_t wy = L2Bot(dy,qy);
1766 for (int dx = 0; dx < L2D1D; ++dx)
1767 {
1768 aXY[dy][dx] += aX[dx] * wy;
1769 }
1770 }
1771 }
1772
1773 for (int dz = 0; dz < L2D1D; ++dz)
1774 {
1775 const real_t wz = L2Bot(dz,qz);
1776 for (int dy = 0; dy < L2D1D; ++dy)
1777 {
1778 for (int dx = 0; dx < L2D1D; ++dx)
1779 {
1780 y(dx,dy,dz,e) += aXY[dy][dx] * wz;
1781 }
1782 }
1783 }
1784 } // loop qz
1785 }); // end of element loop
1786}
1787
1788void PAHdivL2ApplyTranspose3D(const int D1D,
1789 const int Q1D,
1790 const int L2D1D,
1791 const int NE,
1792 const Array<real_t> &L2Bo_,
1793 const Array<real_t> &Gct_,
1794 const Array<real_t> &Bot_,
1795 const Vector &op_,
1796 const Vector &x_,
1797 Vector &y_)
1798{
1799 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
1800 "Error: D1D > HDIV_MAX_D1D");
1801 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
1802 "Error: Q1D > HDIV_MAX_Q1D");
1803 constexpr static int VDIM = 3;
1804
1805 auto L2Bo = Reshape(L2Bo_.Read(), Q1D, L2D1D);
1806 auto Gct = Reshape(Gct_.Read(), D1D, Q1D);
1807 auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
1808 auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, NE);
1809 auto x = Reshape(x_.Read(), L2D1D, L2D1D, L2D1D, NE);
1810 auto y = Reshape(y_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
1811
1812 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1813 {
1814 real_t div[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1815
1816 for (int qz = 0; qz < Q1D; ++qz)
1817 {
1818 for (int qy = 0; qy < Q1D; ++qy)
1819 {
1820 for (int qx = 0; qx < Q1D; ++qx)
1821 {
1822 div[qz][qy][qx] = 0.0;
1823 }
1824 }
1825 }
1826
1827 for (int dz = 0; dz < L2D1D; ++dz)
1828 {
1829 real_t aXY[DofQuadLimits::HDIV_MAX_Q1D][DofQuadLimits::HDIV_MAX_Q1D];
1830 for (int qy = 0; qy < Q1D; ++qy)
1831 {
1832 for (int qx = 0; qx < Q1D; ++qx)
1833 {
1834 aXY[qy][qx] = 0.0;
1835 }
1836 }
1837
1838 for (int dy = 0; dy < L2D1D; ++dy)
1839 {
1840 real_t aX[DofQuadLimits::HDIV_MAX_Q1D];
1841 for (int qx = 0; qx < Q1D; ++qx)
1842 {
1843 aX[qx] = 0.0;
1844 }
1845
1846 for (int dx = 0; dx < L2D1D; ++dx)
1847 {
1848 const real_t t = x(dx,dy,dz,e);
1849 for (int qx = 0; qx < Q1D; ++qx)
1850 {
1851 aX[qx] += t * L2Bo(qx,dx);
1852 }
1853 }
1854
1855 for (int qy = 0; qy < Q1D; ++qy)
1856 {
1857 const real_t wy = L2Bo(qy,dy);
1858 for (int qx = 0; qx < Q1D; ++qx)
1859 {
1860 aXY[qy][qx] += aX[qx] * wy;
1861 }
1862 }
1863 }
1864
1865 for (int qz = 0; qz < Q1D; ++qz)
1866 {
1867 const real_t wz = L2Bo(qz,dz);
1868 for (int qy = 0; qy < Q1D; ++qy)
1869 {
1870 for (int qx = 0; qx < Q1D; ++qx)
1871 {
1872 div[qz][qy][qx] += aXY[qy][qx] * wz;
1873 }
1874 }
1875 }
1876 }
1877
1878 // Apply D operator.
1879 for (int qz = 0; qz < Q1D; ++qz)
1880 {
1881 for (int qy = 0; qy < Q1D; ++qy)
1882 {
1883 for (int qx = 0; qx < Q1D; ++qx)
1884 {
1885 div[qz][qy][qx] *= op(qx,qy,qz,e);
1886 }
1887 }
1888 }
1889
1890 for (int qz = 0; qz < Q1D; ++qz)
1891 {
1892 real_t aXY[DofQuadLimits::HDIV_MAX_D1D][DofQuadLimits::HDIV_MAX_D1D];
1893
1894 int osc = 0;
1895 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1896 {
1897 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1898 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1899 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1900
1901 for (int dy = 0; dy < D1Dy; ++dy)
1902 {
1903 for (int dx = 0; dx < D1Dx; ++dx)
1904 {
1905 aXY[dy][dx] = 0;
1906 }
1907 }
1908 for (int qy = 0; qy < Q1D; ++qy)
1909 {
1910 real_t aX[DofQuadLimits::HDIV_MAX_D1D];
1911 for (int dx = 0; dx < D1Dx; ++dx)
1912 {
1913 aX[dx] = 0;
1914 }
1915 for (int qx = 0; qx < Q1D; ++qx)
1916 {
1917 for (int dx = 0; dx < D1Dx; ++dx)
1918 {
1919 aX[dx] += div[qz][qy][qx] * ((c == 0) ? Gct(dx,qx) :
1920 Bot(dx,qx));
1921 }
1922 }
1923 for (int dy = 0; dy < D1Dy; ++dy)
1924 {
1925 const real_t wy = (c == 1) ? Gct(dy,qy) : Bot(dy,qy);
1926 for (int dx = 0; dx < D1Dx; ++dx)
1927 {
1928 aXY[dy][dx] += aX[dx] * wy;
1929 }
1930 }
1931 }
1932
1933 for (int dz = 0; dz < D1Dz; ++dz)
1934 {
1935 const real_t wz = (c == 2) ? Gct(dz,qz) : Bot(dz,qz);
1936 for (int dy = 0; dy < D1Dy; ++dy)
1937 {
1938 for (int dx = 0; dx < D1Dx; ++dx)
1939 {
1940 y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) +=
1941 aXY[dy][dx] * wz;
1942 }
1943 }
1944 }
1945
1946 osc += D1Dx * D1Dy * D1Dz;
1947 } // loop c
1948 } // loop qz
1949 }); // end of element loop
1950}
1951
1952} // namespace internal
1953
1954} // namespace mfem
bilininteg_hdiv_kernels.hpp
dim
int dim
Definition ex24.cpp:53
a
real_t a
Definition lissajous.cpp:41
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::real_t
float real_t
Definition config.hpp:43
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:753
mfem::DeviceDofQuadLimits::Get
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:128
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0