12 #include "../general/forall.hpp"
26 const Array<double> &w_,
43 constexpr
static int VDIM = 2;
51 auto op =
Reshape(pa_data.
Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
59 for (
int qy = 0; qy < Q1D; ++qy)
61 for (
int qx = 0; qx < Q1D; ++qx)
63 for (
int c = 0; c < VDIM; ++c)
65 mass[qy][qx][c] = 0.0;
72 for (
int c = 0; c < VDIM; ++c)
74 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
75 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
77 for (
int dy = 0; dy < D1Dy; ++dy)
80 for (
int qx = 0; qx < Q1D; ++qx)
85 for (
int dx = 0; dx < D1Dx; ++dx)
87 const double t = X(dx + (dy * D1Dx) + osc, e);
88 for (
int qx = 0; qx < Q1D; ++qx)
90 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
94 for (
int qy = 0; qy < Q1D; ++qy)
96 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
97 for (
int qx = 0; qx < Q1D; ++qx)
99 mass[qy][qx][c] += massX[qx] * wy;
108 for (
int qy = 0; qy < Q1D; ++qy)
110 for (
int qx = 0; qx < Q1D; ++qx)
112 const double O11 = op(qx,qy,0,e);
113 const double O21 = op(qx,qy,1,e);
114 const double O12 = symmetric ? O21 : op(qx,qy,2,e);
115 const double O22 = symmetric ? op(qx,qy,2,e) : op(qx,qy,3,e);
116 const double massX = mass[qy][qx][0];
117 const double massY = mass[qy][qx][1];
118 mass[qy][qx][0] = (O11*massX)+(O12*massY);
119 mass[qy][qx][1] = (O21*massX)+(O22*massY);
123 for (
int qy = 0; qy < Q1D; ++qy)
127 for (
int c = 0; c < VDIM; ++c)
129 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
130 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
133 for (
int dx = 0; dx < D1Dx; ++dx)
137 for (
int qx = 0; qx < Q1D; ++qx)
139 for (
int dx = 0; dx < D1Dx; ++dx)
141 massX[dx] += mass[qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
145 for (
int dy = 0; dy < D1Dy; ++dy)
147 const double wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
149 for (
int dx = 0; dx < D1Dx; ++dx)
151 Y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
164 const bool symmetric,
170 constexpr
static int VDIM = 2;
175 auto op =
Reshape(pa_data.
Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
182 for (
int c = 0; c < VDIM; ++c)
184 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
185 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
189 for (
int dy = 0; dy < D1Dy; ++dy)
191 for (
int qx = 0; qx < Q1D; ++qx)
194 for (
int qy = 0; qy < Q1D; ++qy)
196 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
198 mass[qx] += wy * wy * ((c == 0) ? op(qx,qy,0,e) :
199 op(qx,qy,symmetric ? 2 : 3, e));
203 for (
int dx = 0; dx < D1Dx; ++dx)
205 for (
int qx = 0; qx < Q1D; ++qx)
207 const double wx = ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
208 D(dx + (dy * D1Dx) + osc, e) += mass[qx] * wx * wx;
221 const bool symmetric,
230 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
231 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
232 constexpr static int VDIM = 3;
234 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
235 auto Bc = Reshape(bc.Read(), Q1D, D1D);
236 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
237 auto D = Reshape(diag.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
243 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
245 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
246 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
247 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
249 const int opc = (c == 0) ? 0 : ((c == 1) ? (symmetric ? 3 : 4) :
250 (symmetric ? 5 : 8));
252 double mass[MAX_Q1D];
254 for (int dz = 0; dz < D1Dz; ++dz)
256 for (int dy = 0; dy < D1Dy; ++dy)
258 for (int qx = 0; qx < Q1D; ++qx)
261 for (int qy = 0; qy < Q1D; ++qy)
263 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
265 for (int qz = 0; qz < Q1D; ++qz)
267 const double wz = (c == 2) ? Bo(qz,dz) : Bc(qz,dz);
269 mass[qx] += wy * wy * wz * wz * op(qx,qy,qz,opc,e);
274 for (int dx = 0; dx < D1Dx; ++dx)
276 for (int qx = 0; qx < Q1D; ++qx)
278 const double wx = ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
279 D(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += mass[qx] * wx * wx;
285 osc += D1Dx * D1Dy * D1Dz;
287 }); // end of element loop
290 template<int T_D1D, int T_Q1D>
291 void SmemPAHcurlMassAssembleDiagonal3D(const int D1D,
294 const bool symmetric,
295 const Array<double> &bo,
296 const Array<double> &bc,
297 const Vector &pa_data,
300 MFEM_VERIFY(D1D <= HCURL_MAX_D1D, "Error: D1D > MAX_D1D
");
301 MFEM_VERIFY(Q1D <= HCURL_MAX_Q1D, "Error: Q1D > MAX_Q1D
");
303 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
304 auto Bc = Reshape(bc.Read(), Q1D, D1D);
305 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
306 auto D = Reshape(diag.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
308 MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
310 constexpr int VDIM = 3;
311 constexpr int tD1D = T_D1D ? T_D1D : HCURL_MAX_D1D;
312 constexpr int tQ1D = T_Q1D ? T_Q1D : HCURL_MAX_Q1D;
314 MFEM_SHARED double sBo[tQ1D][tD1D];
315 MFEM_SHARED double sBc[tQ1D][tD1D];
318 MFEM_SHARED double sop[3][tQ1D][tQ1D];
320 MFEM_FOREACH_THREAD(qx,x,Q1D)
322 MFEM_FOREACH_THREAD(qy,y,Q1D)
324 MFEM_FOREACH_THREAD(qz,z,Q1D)
326 op3[0] = op(qx,qy,qz,0,e);
327 op3[1] = op(qx,qy,qz,symmetric ? 3 : 4,e);
328 op3[2] = op(qx,qy,qz,symmetric ? 5 : 8,e);
333 const int tidx = MFEM_THREAD_ID(x);
334 const int tidy = MFEM_THREAD_ID(y);
335 const int tidz = MFEM_THREAD_ID(z);
339 MFEM_FOREACH_THREAD(d,y,D1D)
341 MFEM_FOREACH_THREAD(q,x,Q1D)
354 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
356 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
357 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
358 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
362 for (int qz=0; qz < Q1D; ++qz)
366 for (int i=0; i<3; ++i)
368 sop[i][tidx][tidy] = op3[i];
374 MFEM_FOREACH_THREAD(dz,z,D1Dz)
376 const double wz = ((c == 2) ? sBo[qz][dz] : sBc[qz][dz]);
378 MFEM_FOREACH_THREAD(dy,y,D1Dy)
380 MFEM_FOREACH_THREAD(dx,x,D1Dx)
382 for (int qy = 0; qy < Q1D; ++qy)
384 const double wy = ((c == 1) ? sBo[qy][dy] : sBc[qy][dy]);
386 for (int qx = 0; qx < Q1D; ++qx)
388 const double wx = ((c == 0) ? sBo[qx][dx] : sBc[qx][dx]);
389 dxyz += sop[c][qx][qy] * wx * wx * wy * wy * wz * wz;
399 MFEM_FOREACH_THREAD(dz,z,D1Dz)
401 MFEM_FOREACH_THREAD(dy,y,D1Dy)
403 MFEM_FOREACH_THREAD(dx,x,D1Dx)
405 D(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += dxyz;
410 osc += D1Dx * D1Dy * D1Dz;
412 }); // end of element loop
415 void PAHcurlMassApply3D(const int D1D,
418 const bool symmetric,
419 const Array<double> &bo,
420 const Array<double> &bc,
421 const Array<double> &bot,
422 const Array<double> &bct,
423 const Vector &pa_data,
427 constexpr static int MAX_D1D = HCURL_MAX_D1D;
428 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
430 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
431 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
432 constexpr static int VDIM = 3;
434 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
435 auto Bc = Reshape(bc.Read(), Q1D, D1D);
436 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
437 auto Bct = Reshape(bct.Read(), D1D, Q1D);
438 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
439 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
440 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
444 double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
446 for (int qz = 0; qz < Q1D; ++qz)
448 for (int qy = 0; qy < Q1D; ++qy)
450 for (int qx = 0; qx < Q1D; ++qx)
452 for (int c = 0; c < VDIM; ++c)
454 mass[qz][qy][qx][c] = 0.0;
462 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
464 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
465 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
466 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
468 for (int dz = 0; dz < D1Dz; ++dz)
470 double massXY[MAX_Q1D][MAX_Q1D];
471 for (int qy = 0; qy < Q1D; ++qy)
473 for (int qx = 0; qx < Q1D; ++qx)
475 massXY[qy][qx] = 0.0;
479 for (int dy = 0; dy < D1Dy; ++dy)
481 double massX[MAX_Q1D];
482 for (int qx = 0; qx < Q1D; ++qx)
487 for (int dx = 0; dx < D1Dx; ++dx)
489 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
490 for (int qx = 0; qx < Q1D; ++qx)
492 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
496 for (int qy = 0; qy < Q1D; ++qy)
498 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
499 for (int qx = 0; qx < Q1D; ++qx)
501 const double wx = massX[qx];
502 massXY[qy][qx] += wx * wy;
507 for (int qz = 0; qz < Q1D; ++qz)
509 const double wz = (c == 2) ? Bo(qz,dz) : Bc(qz,dz);
510 for (int qy = 0; qy < Q1D; ++qy)
512 for (int qx = 0; qx < Q1D; ++qx)
514 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
520 osc += D1Dx * D1Dy * D1Dz;
521 } // loop (c) over components
524 for (int qz = 0; qz < Q1D; ++qz)
526 for (int qy = 0; qy < Q1D; ++qy)
528 for (int qx = 0; qx < Q1D; ++qx)
530 const double O11 = op(qx,qy,qz,0,e);
531 const double O12 = op(qx,qy,qz,1,e);
532 const double O13 = op(qx,qy,qz,2,e);
533 const double O21 = symmetric ? O12 : op(qx,qy,qz,3,e);
534 const double O22 = symmetric ? op(qx,qy,qz,3,e) : op(qx,qy,qz,4,e);
535 const double O23 = symmetric ? op(qx,qy,qz,4,e) : op(qx,qy,qz,5,e);
536 const double O31 = symmetric ? O13 : op(qx,qy,qz,6,e);
537 const double O32 = symmetric ? O23 : op(qx,qy,qz,7,e);
538 const double O33 = symmetric ? op(qx,qy,qz,5,e) : op(qx,qy,qz,8,e);
539 const double massX = mass[qz][qy][qx][0];
540 const double massY = mass[qz][qy][qx][1];
541 const double massZ = mass[qz][qy][qx][2];
542 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
543 mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
544 mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
549 for (int qz = 0; qz < Q1D; ++qz)
551 double massXY[MAX_D1D][MAX_D1D];
555 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
557 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
558 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
559 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
561 for (int dy = 0; dy < D1Dy; ++dy)
563 for (int dx = 0; dx < D1Dx; ++dx)
565 massXY[dy][dx] = 0.0;
568 for (int qy = 0; qy < Q1D; ++qy)
570 double massX[MAX_D1D];
571 for (int dx = 0; dx < D1Dx; ++dx)
575 for (int qx = 0; qx < Q1D; ++qx)
577 for (int dx = 0; dx < D1Dx; ++dx)
579 massX[dx] += mass[qz][qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
582 for (int dy = 0; dy < D1Dy; ++dy)
584 const double wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
585 for (int dx = 0; dx < D1Dx; ++dx)
587 massXY[dy][dx] += massX[dx] * wy;
592 for (int dz = 0; dz < D1Dz; ++dz)
594 const double wz = (c == 2) ? Bot(dz,qz) : Bct(dz,qz);
595 for (int dy = 0; dy < D1Dy; ++dy)
597 for (int dx = 0; dx < D1Dx; ++dx)
599 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += massXY[dy][dx] * wz;
604 osc += D1Dx * D1Dy * D1Dz;
607 }); // end of element loop
610 template<int T_D1D, int T_Q1D>
611 void SmemPAHcurlMassApply3D(const int D1D,
614 const bool symmetric,
615 const Array<double> &bo,
616 const Array<double> &bc,
617 const Array<double> &bot,
618 const Array<double> &bct,
619 const Vector &pa_data,
623 MFEM_VERIFY(D1D <= HCURL_MAX_D1D, "Error: D1D > MAX_D1D
");
624 MFEM_VERIFY(Q1D <= HCURL_MAX_Q1D, "Error: Q1D > MAX_Q1D
");
626 const int dataSize = symmetric ? 6 : 9;
628 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
629 auto Bc = Reshape(bc.Read(), Q1D, D1D);
630 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, dataSize, NE);
631 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
632 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
634 MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
636 constexpr int VDIM = 3;
637 constexpr int tD1D = T_D1D ? T_D1D : HCURL_MAX_D1D;
638 constexpr int tQ1D = T_Q1D ? T_Q1D : HCURL_MAX_Q1D;
640 MFEM_SHARED double sBo[tQ1D][tD1D];
641 MFEM_SHARED double sBc[tQ1D][tD1D];
644 MFEM_SHARED double sop[9*tQ1D*tQ1D];
645 MFEM_SHARED double mass[tQ1D][tQ1D][3];
647 MFEM_SHARED double sX[tD1D][tD1D][tD1D];
649 MFEM_FOREACH_THREAD(qx,x,Q1D)
651 MFEM_FOREACH_THREAD(qy,y,Q1D)
653 MFEM_FOREACH_THREAD(qz,z,Q1D)
655 for (int i=0; i<dataSize; ++i)
657 op9[i] = op(qx,qy,qz,i,e);
663 const int tidx = MFEM_THREAD_ID(x);
664 const int tidy = MFEM_THREAD_ID(y);
665 const int tidz = MFEM_THREAD_ID(z);
669 MFEM_FOREACH_THREAD(d,y,D1D)
671 MFEM_FOREACH_THREAD(q,x,Q1D)
683 for (int qz=0; qz < Q1D; ++qz)
686 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
688 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
689 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
690 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
692 MFEM_FOREACH_THREAD(dz,z,D1Dz)
694 MFEM_FOREACH_THREAD(dy,y,D1Dy)
696 MFEM_FOREACH_THREAD(dx,x,D1Dx)
698 sX[dz][dy][dx] = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
706 for (int i=0; i<dataSize; ++i)
708 sop[i + (dataSize*tidx) + (dataSize*Q1D*tidy)] = op9[i];
711 MFEM_FOREACH_THREAD(qy,y,Q1D)
713 MFEM_FOREACH_THREAD(qx,x,Q1D)
717 for (int dz = 0; dz < D1Dz; ++dz)
719 const double wz = (c == 2) ? sBo[qz][dz] : sBc[qz][dz];
720 for (int dy = 0; dy < D1Dy; ++dy)
722 const double wy = (c == 1) ? sBo[qy][dy] : sBc[qy][dy];
723 for (int dx = 0; dx < D1Dx; ++dx)
725 const double t = sX[dz][dy][dx];
726 const double wx = (c == 0) ? sBo[qx][dx] : sBc[qx][dx];
727 u += t * wx * wy * wz;
737 osc += D1Dx * D1Dy * D1Dz;
741 MFEM_SYNC_THREAD; // Sync mass[qy][qx][d] and sop
744 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
746 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
747 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
748 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
752 MFEM_FOREACH_THREAD(dz,z,D1Dz)
754 const double wz = (c == 2) ? sBo[qz][dz] : sBc[qz][dz];
756 MFEM_FOREACH_THREAD(dy,y,D1Dy)
758 MFEM_FOREACH_THREAD(dx,x,D1Dx)
760 for (int qy = 0; qy < Q1D; ++qy)
762 const double wy = (c == 1) ? sBo[qy][dy] : sBc[qy][dy];
763 for (int qx = 0; qx < Q1D; ++qx)
765 const int os = (dataSize*qx) + (dataSize*Q1D*qy);
766 const int id1 = os + ((c == 0) ? 0 : ((c == 1) ? (symmetric ? 1 : 3) :
767 (symmetric ? 2 : 6))); // O11, O21, O31
768 const int id2 = os + ((c == 0) ? 1 : ((c == 1) ? (symmetric ? 3 : 4) :
769 (symmetric ? 4 : 7))); // O12, O22, O32
770 const int id3 = os + ((c == 0) ? 2 : ((c == 1) ? (symmetric ? 4 : 5) :
771 (symmetric ? 5 : 8))); // O13, O23, O33
773 const double m_c = (sop[id1] * mass[qy][qx][0]) + (sop[id2] * mass[qy][qx][1]) +
774 (sop[id3] * mass[qy][qx][2]);
776 const double wx = (c == 0) ? sBo[qx][dx] : sBc[qx][dx];
777 dxyz += m_c * wx * wy * wz;
786 MFEM_FOREACH_THREAD(dz,z,D1Dz)
788 MFEM_FOREACH_THREAD(dy,y,D1Dy)
790 MFEM_FOREACH_THREAD(dx,x,D1Dx)
792 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += dxyz;
797 osc += D1Dx * D1Dy * D1Dz;
800 }); // end of element loop
803 // PA H(curl) curl-curl assemble 2D kernel
804 static void PACurlCurlSetup2D(const int Q1D,
806 const Array<double> &w,
811 const int NQ = Q1D*Q1D;
813 auto J = Reshape(j.Read(), NQ, 2, 2, NE);
814 auto C = Reshape(coeff.Read(), NQ, NE);
815 auto y = Reshape(op.Write(), NQ, NE);
818 for (int q = 0; q < NQ; ++q)
820 const double J11 = J(q,0,0,e);
821 const double J21 = J(q,1,0,e);
822 const double J12 = J(q,0,1,e);
823 const double J22 = J(q,1,1,e);
824 const double detJ = (J11*J22)-(J21*J12);
825 y(q,e) = W[q] * C(q,e) / detJ;
830 // PA H(curl) curl-curl assemble 3D kernel
831 static void PACurlCurlSetup3D(const int Q1D,
834 const Array<double> &w,
839 const int NQ = Q1D*Q1D*Q1D;
840 const bool symmetric = (coeffDim != 9);
842 auto J = Reshape(j.Read(), NQ, 3, 3, NE);
843 auto C = Reshape(coeff.Read(), coeffDim, NQ, NE);
844 auto y = Reshape(op.Write(), NQ, symmetric ? 6 : 9, NE);
848 for (int q = 0; q < NQ; ++q)
850 const double J11 = J(q,0,0,e);
851 const double J21 = J(q,1,0,e);
852 const double J31 = J(q,2,0,e);
853 const double J12 = J(q,0,1,e);
854 const double J22 = J(q,1,1,e);
855 const double J32 = J(q,2,1,e);
856 const double J13 = J(q,0,2,e);
857 const double J23 = J(q,1,2,e);
858 const double J33 = J(q,2,2,e);
859 const double detJ = J11 * (J22 * J33 - J32 * J23) -
860 /* */ J21 * (J12 * J33 - J32 * J13) +
861 /* */ J31 * (J12 * J23 - J22 * J13);
863 const double c_detJ = W[q] / detJ;
865 if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
867 // Set y to the 6 or 9 entries of J^T M J / det
868 const double M11 = C(0, q, e);
869 const double M12 = C(1, q, e);
870 const double M13 = C(2, q, e);
871 const double M21 = (!symmetric) ? C(3, q, e) : M12;
872 const double M22 = (!symmetric) ? C(4, q, e) : C(3, q, e);
873 const double M23 = (!symmetric) ? C(5, q, e) : C(4, q, e);
874 const double M31 = (!symmetric) ? C(6, q, e) : M13;
875 const double M32 = (!symmetric) ? C(7, q, e) : M23;
876 const double M33 = (!symmetric) ? C(8, q, e) : C(5, q, e);
878 // First compute R = MJ
879 const double R11 = M11*J11 + M12*J21 + M13*J31;
880 const double R12 = M11*J12 + M12*J22 + M13*J32;
881 const double R13 = M11*J13 + M12*J23 + M13*J33;
882 const double R21 = M21*J11 + M22*J21 + M23*J31;
883 const double R22 = M21*J12 + M22*J22 + M23*J32;
884 const double R23 = M21*J13 + M22*J23 + M23*J33;
885 const double R31 = M31*J11 + M32*J21 + M33*J31;
886 const double R32 = M31*J12 + M32*J22 + M33*J32;
887 const double R33 = M31*J13 + M32*J23 + M33*J33;
889 // Now set y to J^T R / det
890 y(q,0,e) = c_detJ * (J11*R11 + J21*R21 + J31*R31); // 1,1
891 const double Y12 = c_detJ * (J11*R12 + J21*R22 + J31*R32);
892 y(q,1,e) = Y12; // 1,2
893 y(q,2,e) = c_detJ * (J11*R13 + J21*R23 + J31*R33); // 1,3
895 const double Y21 = c_detJ * (J12*R11 + J22*R21 + J32*R31);
896 const double Y22 = c_detJ * (J12*R12 + J22*R22 + J32*R32);
897 const double Y23 = c_detJ * (J12*R13 + J22*R23 + J32*R33);
899 const double Y33 = c_detJ * (J13*R13 + J23*R23 + J33*R33);
901 y(q,3,e) = symmetric ? Y22 : Y21; // 2,2 or 2,1
902 y(q,4,e) = symmetric ? Y23 : Y22; // 2,3 or 2,2
903 y(q,5,e) = symmetric ? Y33 : Y23; // 3,3 or 2,3
907 y(q,6,e) = c_detJ * (J13*R11 + J23*R21 + J33*R31); // 3,1
908 y(q,7,e) = c_detJ * (J13*R12 + J23*R22 + J33*R32); // 3,2
909 y(q,8,e) = Y33; // 3,3
912 else // Vector or scalar coefficient version
914 // Set y to the 6 entries of J^T D J / det^2
915 const double D1 = C(0, q, e);
916 const double D2 = coeffDim == 3 ? C(1, q, e) : D1;
917 const double D3 = coeffDim == 3 ? C(2, q, e) : D1;
919 y(q,0,e) = c_detJ * (D1*J11*J11 + D2*J21*J21 + D3*J31*J31); // 1,1
920 y(q,1,e) = c_detJ * (D1*J11*J12 + D2*J21*J22 + D3*J31*J32); // 1,2
921 y(q,2,e) = c_detJ * (D1*J11*J13 + D2*J21*J23 + D3*J31*J33); // 1,3
922 y(q,3,e) = c_detJ * (D1*J12*J12 + D2*J22*J22 + D3*J32*J32); // 2,2
923 y(q,4,e) = c_detJ * (D1*J12*J13 + D2*J22*J23 + D3*J32*J33); // 2,3
924 y(q,5,e) = c_detJ * (D1*J13*J13 + D2*J23*J23 + D3*J33*J33); // 3,3
930 // PA H(curl)-L2 assemble 2D kernel
931 static void PACurlL2Setup2D(const int Q1D,
933 const Array<double> &w,
937 const int NQ = Q1D*Q1D;
939 auto C = Reshape(coeff.Read(), NQ, NE);
940 auto y = Reshape(op.Write(), NQ, NE);
943 for (int q = 0; q < NQ; ++q)
945 y(q,e) = W[q] * C(q,e);
950 void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
952 // Assumes tensor-product elements
953 Mesh *mesh = fes.GetMesh();
954 const FiniteElement *fel = fes.GetFE(0);
956 const VectorTensorFiniteElement *el =
957 dynamic_cast<const VectorTensorFiniteElement*>(fel);
960 const IntegrationRule *ir
961 = IntRule ? IntRule : &MassIntegrator::GetRule(*el, *el,
962 *mesh->GetElementTransformation(0));
964 const int dims = el->GetDim();
965 MFEM_VERIFY(dims == 2 || dims == 3, "");
967 nq = ir->GetNPoints();
968 dim = mesh->Dimension();
969 MFEM_VERIFY(dim == 2 || dim == 3, "");
972 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
973 mapsC = &el->GetDofToQuad(*ir, DofToQuad::TENSOR);
974 mapsO = &el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
975 dofs1D = mapsC->ndof;
976 quad1D = mapsC->nqpt;
978 MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
980 QuadratureSpace qs(*mesh, *ir);
981 CoefficientVector coeff(qs, CoefficientStorage::SYMMETRIC);
982 if (Q) { coeff.Project(*Q); }
983 else if (MQ) { coeff.ProjectTranspose(*MQ); }
984 else if (DQ) { coeff.Project(*DQ); }
985 else { coeff.SetConstant(1.0); }
987 const int coeff_dim = coeff.GetVDim();
988 symmetric = (coeff_dim != dim*dim);
989 const int sym_dims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
990 const int ndata = (dim == 2) ? 1 : (symmetric ? sym_dims : dim*dim);
991 pa_data.SetSize(ndata * nq * ne, Device::GetMemoryType());
993 if (el->GetDerivType() != mfem::FiniteElement::CURL)
995 MFEM_ABORT("Unknown kernel.
");
1000 PACurlCurlSetup3D(quad1D, coeff_dim, ne, ir->GetWeights(), geom->J, coeff,
1005 PACurlCurlSetup2D(quad1D, ne, ir->GetWeights(), geom->J, coeff, pa_data);
1009 static void PACurlCurlApply2D(const int D1D,
1012 const Array<double> &bo,
1013 const Array<double> &bot,
1014 const Array<double> &gc,
1015 const Array<double> &gct,
1016 const Vector &pa_data,
1020 constexpr static int VDIM = 2;
1021 constexpr static int MAX_D1D = HCURL_MAX_D1D;
1022 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
1024 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
1025 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
1026 auto Gc = Reshape(gc.Read(), Q1D, D1D);
1027 auto Gct = Reshape(gct.Read(), D1D, Q1D);
1028 auto op = Reshape(pa_data.Read(), Q1D, Q1D, NE);
1029 auto X = Reshape(x.Read(), 2*(D1D-1)*D1D, NE);
1030 auto Y = Reshape(y.ReadWrite(), 2*(D1D-1)*D1D, NE);
1034 double curl[MAX_Q1D][MAX_Q1D];
1036 // curl[qy][qx] will be computed as du_y/dx - du_x/dy
1038 for (int qy = 0; qy < Q1D; ++qy)
1040 for (int qx = 0; qx < Q1D; ++qx)
1048 for (int c = 0; c < VDIM; ++c) // loop over x, y components
1050 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
1051 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
1053 for (int dy = 0; dy < D1Dy; ++dy)
1055 double gradX[MAX_Q1D];
1056 for (int qx = 0; qx < Q1D; ++qx)
1061 for (int dx = 0; dx < D1Dx; ++dx)
1063 const double t = X(dx + (dy * D1Dx) + osc, e);
1064 for (int qx = 0; qx < Q1D; ++qx)
1066 gradX[qx] += t * ((c == 0) ? Bo(qx,dx) : Gc(qx,dx));
1070 for (int qy = 0; qy < Q1D; ++qy)
1072 const double wy = (c == 0) ? -Gc(qy,dy) : Bo(qy,dy);
1073 for (int qx = 0; qx < Q1D; ++qx)
1075 curl[qy][qx] += gradX[qx] * wy;
1081 } // loop (c) over components
1083 // Apply D operator.
1084 for (int qy = 0; qy < Q1D; ++qy)
1086 for (int qx = 0; qx < Q1D; ++qx)
1088 curl[qy][qx] *= op(qx,qy,e);
1092 for (int qy = 0; qy < Q1D; ++qy)
1096 for (int c = 0; c < VDIM; ++c) // loop over x, y components
1098 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
1099 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
1101 double gradX[MAX_D1D];
1102 for (int dx = 0; dx < D1Dx; ++dx)
1106 for (int qx = 0; qx < Q1D; ++qx)
1108 for (int dx = 0; dx < D1Dx; ++dx)
1110 gradX[dx] += curl[qy][qx] * ((c == 0) ? Bot(dx,qx) : Gct(dx,qx));
1113 for (int dy = 0; dy < D1Dy; ++dy)
1115 const double wy = (c == 0) ? -Gct(dy,qy) : Bot(dy,qy);
1117 for (int dx = 0; dx < D1Dx; ++dx)
1119 Y(dx + (dy * D1Dx) + osc, e) += gradX[dx] * wy;
1126 }); // end of element loop
1129 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
1130 static void PACurlCurlApply3D(const int D1D,
1132 const bool symmetric,
1134 const Array<double> &bo,
1135 const Array<double> &bc,
1136 const Array<double> &bot,
1137 const Array<double> &bct,
1138 const Array<double> &gc,
1139 const Array<double> &gct,
1140 const Vector &pa_data,
1144 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
1145 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
1146 // Using (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u} (p. 78 of Monk), we get
1147 // (\nabla\times u) \cdot (\nabla\times v) = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{\nabla}\times\hat{v}
1148 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
1149 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
1150 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
1152 constexpr static int VDIM = 3;
1154 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
1155 auto Bc = Reshape(bc.Read(), Q1D, D1D);
1156 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
1157 auto Bct = Reshape(bct.Read(), D1D, Q1D);
1158 auto Gc = Reshape(gc.Read(), Q1D, D1D);
1159 auto Gct = Reshape(gct.Read(), D1D, Q1D);
1160 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, (symmetric ? 6 : 9), NE);
1161 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
1162 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
1166 double curl[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
1167 // curl[qz][qy][qx] will be computed as the vector curl at each quadrature point.
1169 for (int qz = 0; qz < Q1D; ++qz)
1171 for (int qy = 0; qy < Q1D; ++qy)
1173 for (int qx = 0; qx < Q1D; ++qx)
1175 for (int c = 0; c < VDIM; ++c)
1177 curl[qz][qy][qx][c] = 0.0;
1183 // We treat x, y, z components separately for optimization specific to each.
1189 const int D1Dz = D1D;
1190 const int D1Dy = D1D;
1191 const int D1Dx = D1D - 1;
1193 for (int dz = 0; dz < D1Dz; ++dz)
1195 double gradXY[MAX_Q1D][MAX_Q1D][2];
1196 for (int qy = 0; qy < Q1D; ++qy)
1198 for (int qx = 0; qx < Q1D; ++qx)
1200 for (int d = 0; d < 2; ++d)
1202 gradXY[qy][qx][d] = 0.0;
1207 for (int dy = 0; dy < D1Dy; ++dy)
1209 double massX[MAX_Q1D];
1210 for (int qx = 0; qx < Q1D; ++qx)
1215 for (int dx = 0; dx < D1Dx; ++dx)
1217 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1218 for (int qx = 0; qx < Q1D; ++qx)
1220 massX[qx] += t * Bo(qx,dx);
1224 for (int qy = 0; qy < Q1D; ++qy)
1226 const double wy = Bc(qy,dy);
1227 const double wDy = Gc(qy,dy);
1228 for (int qx = 0; qx < Q1D; ++qx)
1230 const double wx = massX[qx];
1231 gradXY[qy][qx][0] += wx * wDy;
1232 gradXY[qy][qx][1] += wx * wy;
1237 for (int qz = 0; qz < Q1D; ++qz)
1239 const double wz = Bc(qz,dz);
1240 const double wDz = Gc(qz,dz);
1241 for (int qy = 0; qy < Q1D; ++qy)
1243 for (int qx = 0; qx < Q1D; ++qx)
1245 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
1246 curl[qz][qy][qx][1] += gradXY[qy][qx][1] * wDz; // (u_0)_{x_2}
1247 curl[qz][qy][qx][2] -= gradXY[qy][qx][0] * wz; // -(u_0)_{x_1}
1253 osc += D1Dx * D1Dy * D1Dz;
1258 const int D1Dz = D1D;
1259 const int D1Dy = D1D - 1;
1260 const int D1Dx = D1D;
1262 for (int dz = 0; dz < D1Dz; ++dz)
1264 double gradXY[MAX_Q1D][MAX_Q1D][2];
1265 for (int qy = 0; qy < Q1D; ++qy)
1267 for (int qx = 0; qx < Q1D; ++qx)
1269 for (int d = 0; d < 2; ++d)
1271 gradXY[qy][qx][d] = 0.0;
1276 for (int dx = 0; dx < D1Dx; ++dx)
1278 double massY[MAX_Q1D];
1279 for (int qy = 0; qy < Q1D; ++qy)
1284 for (int dy = 0; dy < D1Dy; ++dy)
1286 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1287 for (int qy = 0; qy < Q1D; ++qy)
1289 massY[qy] += t * Bo(qy,dy);
1293 for (int qx = 0; qx < Q1D; ++qx)
1295 const double wx = Bc(qx,dx);
1296 const double wDx = Gc(qx,dx);
1297 for (int qy = 0; qy < Q1D; ++qy)
1299 const double wy = massY[qy];
1300 gradXY[qy][qx][0] += wDx * wy;
1301 gradXY[qy][qx][1] += wx * wy;
1306 for (int qz = 0; qz < Q1D; ++qz)
1308 const double wz = Bc(qz,dz);
1309 const double wDz = Gc(qz,dz);
1310 for (int qy = 0; qy < Q1D; ++qy)
1312 for (int qx = 0; qx < Q1D; ++qx)
1314 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
1315 curl[qz][qy][qx][0] -= gradXY[qy][qx][1] * wDz; // -(u_1)_{x_2}
1316 curl[qz][qy][qx][2] += gradXY[qy][qx][0] * wz; // (u_1)_{x_0}
1322 osc += D1Dx * D1Dy * D1Dz;
1327 const int D1Dz = D1D - 1;
1328 const int D1Dy = D1D;
1329 const int D1Dx = D1D;
1331 for (int dx = 0; dx < D1Dx; ++dx)
1333 double gradYZ[MAX_Q1D][MAX_Q1D][2];
1334 for (int qz = 0; qz < Q1D; ++qz)
1336 for (int qy = 0; qy < Q1D; ++qy)
1338 for (int d = 0; d < 2; ++d)
1340 gradYZ[qz][qy][d] = 0.0;
1345 for (int dy = 0; dy < D1Dy; ++dy)
1347 double massZ[MAX_Q1D];
1348 for (int qz = 0; qz < Q1D; ++qz)
1353 for (int dz = 0; dz < D1Dz; ++dz)
1355 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1356 for (int qz = 0; qz < Q1D; ++qz)
1358 massZ[qz] += t * Bo(qz,dz);
1362 for (int qy = 0; qy < Q1D; ++qy)
1364 const double wy = Bc(qy,dy);
1365 const double wDy = Gc(qy,dy);
1366 for (int qz = 0; qz < Q1D; ++qz)
1368 const double wz = massZ[qz];
1369 gradYZ[qz][qy][0] += wz * wy;
1370 gradYZ[qz][qy][1] += wz * wDy;
1375 for (int qx = 0; qx < Q1D; ++qx)
1377 const double wx = Bc(qx,dx);
1378 const double wDx = Gc(qx,dx);
1380 for (int qy = 0; qy < Q1D; ++qy)
1382 for (int qz = 0; qz < Q1D; ++qz)
1384 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
1385 curl[qz][qy][qx][0] += gradYZ[qz][qy][1] * wx; // (u_2)_{x_1}
1386 curl[qz][qy][qx][1] -= gradYZ[qz][qy][0] * wDx; // -(u_2)_{x_0}
1393 // Apply D operator.
1394 for (int qz = 0; qz < Q1D; ++qz)
1396 for (int qy = 0; qy < Q1D; ++qy)
1398 for (int qx = 0; qx < Q1D; ++qx)
1400 const double O11 = op(qx,qy,qz,0,e);
1401 const double O12 = op(qx,qy,qz,1,e);
1402 const double O13 = op(qx,qy,qz,2,e);
1403 const double O21 = symmetric ? O12 : op(qx,qy,qz,3,e);
1404 const double O22 = symmetric ? op(qx,qy,qz,3,e) : op(qx,qy,qz,4,e);
1405 const double O23 = symmetric ? op(qx,qy,qz,4,e) : op(qx,qy,qz,5,e);
1406 const double O31 = symmetric ? O13 : op(qx,qy,qz,6,e);
1407 const double O32 = symmetric ? O23 : op(qx,qy,qz,7,e);
1408 const double O33 = symmetric ? op(qx,qy,qz,5,e) : op(qx,qy,qz,8,e);
1410 const double c1 = (O11 * curl[qz][qy][qx][0]) + (O12 * curl[qz][qy][qx][1]) +
1411 (O13 * curl[qz][qy][qx][2]);
1412 const double c2 = (O21 * curl[qz][qy][qx][0]) + (O22 * curl[qz][qy][qx][1]) +
1413 (O23 * curl[qz][qy][qx][2]);
1414 const double c3 = (O31 * curl[qz][qy][qx][0]) + (O32 * curl[qz][qy][qx][1]) +
1415 (O33 * curl[qz][qy][qx][2]);
1417 curl[qz][qy][qx][0] = c1;
1418 curl[qz][qy][qx][1] = c2;
1419 curl[qz][qy][qx][2] = c3;
1427 const int D1Dz = D1D;
1428 const int D1Dy = D1D;
1429 const int D1Dx = D1D - 1;
1431 for (int qz = 0; qz < Q1D; ++qz)
1433 double gradXY12[MAX_D1D][MAX_D1D];
1434 double gradXY21[MAX_D1D][MAX_D1D];
1436 for (int dy = 0; dy < D1Dy; ++dy)
1438 for (int dx = 0; dx < D1Dx; ++dx)
1440 gradXY12[dy][dx] = 0.0;
1441 gradXY21[dy][dx] = 0.0;
1444 for (int qy = 0; qy < Q1D; ++qy)
1446 double massX[MAX_D1D][2];
1447 for (int dx = 0; dx < D1Dx; ++dx)
1449 for (int n = 0; n < 2; ++n)
1454 for (int qx = 0; qx < Q1D; ++qx)
1456 for (int dx = 0; dx < D1Dx; ++dx)
1458 const double wx = Bot(dx,qx);
1460 massX[dx][0] += wx * curl[qz][qy][qx][1];
1461 massX[dx][1] += wx * curl[qz][qy][qx][2];
1464 for (int dy = 0; dy < D1Dy; ++dy)
1466 const double wy = Bct(dy,qy);
1467 const double wDy = Gct(dy,qy);
1469 for (int dx = 0; dx < D1Dx; ++dx)
1471 gradXY21[dy][dx] += massX[dx][0] * wy;
1472 gradXY12[dy][dx] += massX[dx][1] * wDy;
1477 for (int dz = 0; dz < D1Dz; ++dz)
1479 const double wz = Bct(dz,qz);
1480 const double wDz = Gct(dz,qz);
1481 for (int dy = 0; dy < D1Dy; ++dy)
1483 for (int dx = 0; dx < D1Dx; ++dx)
1485 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
1486 // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
1487 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
1488 e) += (gradXY21[dy][dx] * wDz) - (gradXY12[dy][dx] * wz);
1494 osc += D1Dx * D1Dy * D1Dz;
1499 const int D1Dz = D1D;
1500 const int D1Dy = D1D - 1;
1501 const int D1Dx = D1D;
1503 for (int qz = 0; qz < Q1D; ++qz)
1505 double gradXY02[MAX_D1D][MAX_D1D];
1506 double gradXY20[MAX_D1D][MAX_D1D];
1508 for (int dy = 0; dy < D1Dy; ++dy)
1510 for (int dx = 0; dx < D1Dx; ++dx)
1512 gradXY02[dy][dx] = 0.0;
1513 gradXY20[dy][dx] = 0.0;
1516 for (int qx = 0; qx < Q1D; ++qx)
1518 double massY[MAX_D1D][2];
1519 for (int dy = 0; dy < D1Dy; ++dy)
1524 for (int qy = 0; qy < Q1D; ++qy)
1526 for (int dy = 0; dy < D1Dy; ++dy)
1528 const double wy = Bot(dy,qy);
1530 massY[dy][0] += wy * curl[qz][qy][qx][2];
1531 massY[dy][1] += wy * curl[qz][qy][qx][0];
1534 for (int dx = 0; dx < D1Dx; ++dx)
1536 const double wx = Bct(dx,qx);
1537 const double wDx = Gct(dx,qx);
1539 for (int dy = 0; dy < D1Dy; ++dy)
1541 gradXY02[dy][dx] += massY[dy][0] * wDx;
1542 gradXY20[dy][dx] += massY[dy][1] * wx;
1547 for (int dz = 0; dz < D1Dz; ++dz)
1549 const double wz = Bct(dz,qz);
1550 const double wDz = Gct(dz,qz);
1551 for (int dy = 0; dy < D1Dy; ++dy)
1553 for (int dx = 0; dx < D1Dx; ++dx)
1555 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
1556 // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
1557 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
1558 e) += (-gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
1564 osc += D1Dx * D1Dy * D1Dz;
1569 const int D1Dz = D1D - 1;
1570 const int D1Dy = D1D;
1571 const int D1Dx = D1D;
1573 for (int qx = 0; qx < Q1D; ++qx)
1575 double gradYZ01[MAX_D1D][MAX_D1D];
1576 double gradYZ10[MAX_D1D][MAX_D1D];
1578 for (int dy = 0; dy < D1Dy; ++dy)
1580 for (int dz = 0; dz < D1Dz; ++dz)
1582 gradYZ01[dz][dy] = 0.0;
1583 gradYZ10[dz][dy] = 0.0;
1586 for (int qy = 0; qy < Q1D; ++qy)
1588 double massZ[MAX_D1D][2];
1589 for (int dz = 0; dz < D1Dz; ++dz)
1591 for (int n = 0; n < 2; ++n)
1596 for (int qz = 0; qz < Q1D; ++qz)
1598 for (int dz = 0; dz < D1Dz; ++dz)
1600 const double wz = Bot(dz,qz);
1602 massZ[dz][0] += wz * curl[qz][qy][qx][0];
1603 massZ[dz][1] += wz * curl[qz][qy][qx][1];
1606 for (int dy = 0; dy < D1Dy; ++dy)
1608 const double wy = Bct(dy,qy);
1609 const double wDy = Gct(dy,qy);
1611 for (int dz = 0; dz < D1Dz; ++dz)
1613 gradYZ01[dz][dy] += wy * massZ[dz][1];
1614 gradYZ10[dz][dy] += wDy * massZ[dz][0];
1619 for (int dx = 0; dx < D1Dx; ++dx)
1621 const double wx = Bct(dx,qx);
1622 const double wDx = Gct(dx,qx);
1624 for (int dy = 0; dy < D1Dy; ++dy)
1626 for (int dz = 0; dz < D1Dz; ++dz)
1628 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
1629 // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
1630 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
1631 e) += (gradYZ10[dz][dy] * wx) - (gradYZ01[dz][dy] * wDx);
1637 }); // end of element loop
1640 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
1641 static void SmemPACurlCurlApply3D(const int D1D,
1643 const bool symmetric,
1645 const Array<double> &bo,
1646 const Array<double> &bc,
1647 const Array<double> &bot,
1648 const Array<double> &bct,
1649 const Array<double> &gc,
1650 const Array<double> &gct,
1651 const Vector &pa_data,
1655 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
1656 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
1657 // Using (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u} (p. 78 of Monk), we get
1658 // (\nabla\times u) \cdot (\nabla\times v) = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{\nabla}\times\hat{v}
1659 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
1660 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
1661 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
1663 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
1664 auto Bc = Reshape(bc.Read(), Q1D, D1D);
1665 auto Gc = Reshape(gc.Read(), Q1D, D1D);
1666 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
1667 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
1668 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
1670 const int s = symmetric ? 6 : 9;
1672 auto device_kernel = [=] MFEM_DEVICE (int e)
1674 constexpr int VDIM = 3;
1676 MFEM_SHARED double sBo[MAX_D1D][MAX_Q1D];
1677 MFEM_SHARED double sBc[MAX_D1D][MAX_Q1D];
1678 MFEM_SHARED double sGc[MAX_D1D][MAX_Q1D];
1681 MFEM_SHARED double sop[9][MAX_Q1D][MAX_Q1D];
1682 MFEM_SHARED double curl[MAX_Q1D][MAX_Q1D][3];
1684 MFEM_SHARED double sX[MAX_D1D][MAX_D1D][MAX_D1D];
1686 MFEM_FOREACH_THREAD(qx,x,Q1D)
1688 MFEM_FOREACH_THREAD(qy,y,Q1D)
1690 MFEM_FOREACH_THREAD(qz,z,Q1D)
1692 for (int i=0; i<s; ++i)
1694 ope[i] = op(qx,qy,qz,i,e);
1700 const int tidx = MFEM_THREAD_ID(x);
1701 const int tidy = MFEM_THREAD_ID(y);
1702 const int tidz = MFEM_THREAD_ID(z);
1706 MFEM_FOREACH_THREAD(d,y,D1D)
1708 MFEM_FOREACH_THREAD(q,x,Q1D)
1710 sBc[d][q] = Bc(q,d);
1711 sGc[d][q] = Gc(q,d);
1714 sBo[d][q] = Bo(q,d);
1721 for (int qz=0; qz < Q1D; ++qz)
1725 MFEM_FOREACH_THREAD(qy,y,Q1D)
1727 MFEM_FOREACH_THREAD(qx,x,Q1D)
1729 for (int i=0; i<3; ++i)
1731 curl[qy][qx][i] = 0.0;
1738 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1740 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
1741 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
1742 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
1744 MFEM_FOREACH_THREAD(dz,z,D1Dz)
1746 MFEM_FOREACH_THREAD(dy,y,D1Dy)
1748 MFEM_FOREACH_THREAD(dx,x,D1Dx)
1750 sX[dz][dy][dx] = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1760 for (int i=0; i<s; ++i)
1762 sop[i][tidx][tidy] = ope[i];
1766 MFEM_FOREACH_THREAD(qy,y,Q1D)
1768 MFEM_FOREACH_THREAD(qx,x,Q1D)
1773 // We treat x, y, z components separately for optimization specific to each.
1774 if (c == 0) // x component
1776 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
1778 for (int dz = 0; dz < D1Dz; ++dz)
1780 const double wz = sBc[dz][qz];
1781 const double wDz = sGc[dz][qz];
1783 for (int dy = 0; dy < D1Dy; ++dy)
1785 const double wy = sBc[dy][qy];
1786 const double wDy = sGc[dy][qy];
1788 for (int dx = 0; dx < D1Dx; ++dx)
1790 const double wx = sX[dz][dy][dx] * sBo[dx][qx];
1797 curl[qy][qx][1] += v; // (u_0)_{x_2}
1798 curl[qy][qx][2] -= u; // -(u_0)_{x_1}
1800 else if (c == 1) // y component
1802 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
1804 for (int dz = 0; dz < D1Dz; ++dz)
1806 const double wz = sBc[dz][qz];
1807 const double wDz = sGc[dz][qz];
1809 for (int dy = 0; dy < D1Dy; ++dy)
1811 const double wy = sBo[dy][qy];
1813 for (int dx = 0; dx < D1Dx; ++dx)
1815 const double t = sX[dz][dy][dx];
1816 const double wx = t * sBc[dx][qx];
1817 const double wDx = t * sGc[dx][qx];
1825 curl[qy][qx][0] -= v; // -(u_1)_{x_2}
1826 curl[qy][qx][2] += u; // (u_1)_{x_0}
1830 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
1832 for (int dz = 0; dz < D1Dz; ++dz)
1834 const double wz = sBo[dz][qz];
1836 for (int dy = 0; dy < D1Dy; ++dy)
1838 const double wy = sBc[dy][qy];
1839 const double wDy = sGc[dy][qy];
1841 for (int dx = 0; dx < D1Dx; ++dx)
1843 const double t = sX[dz][dy][dx];
1844 const double wx = t * sBc[dx][qx];
1845 const double wDx = t * sGc[dx][qx];
1853 curl[qy][qx][0] += v; // (u_2)_{x_1}
1854 curl[qy][qx][1] -= u; // -(u_2)_{x_0}
1860 osc += D1Dx * D1Dy * D1Dz;
1868 MFEM_FOREACH_THREAD(dz,z,D1D)
1870 const double wcz = sBc[dz][qz];
1871 const double wcDz = sGc[dz][qz];
1872 const double wz = (dz < D1D-1) ? sBo[dz][qz] : 0.0;
1874 MFEM_FOREACH_THREAD(dy,y,D1D)
1876 MFEM_FOREACH_THREAD(dx,x,D1D)
1878 for (int qy = 0; qy < Q1D; ++qy)
1880 const double wcy = sBc[dy][qy];
1881 const double wcDy = sGc[dy][qy];
1882 const double wy = (dy < D1D-1) ? sBo[dy][qy] : 0.0;
1884 for (int qx = 0; qx < Q1D; ++qx)
1886 const double O11 = sop[0][qx][qy];
1887 const double O12 = sop[1][qx][qy];
1888 const double O13 = sop[2][qx][qy];
1889 const double O21 = symmetric ? O12 : sop[3][qx][qy];
1890 const double O22 = symmetric ? sop[3][qx][qy] : sop[4][qx][qy];
1891 const double O23 = symmetric ? sop[4][qx][qy] : sop[5][qx][qy];
1892 const double O31 = symmetric ? O13 : sop[6][qx][qy];
1893 const double O32 = symmetric ? O23 : sop[7][qx][qy];
1894 const double O33 = symmetric ? sop[5][qx][qy] : sop[8][qx][qy];
1896 const double c1 = (O11 * curl[qy][qx][0]) + (O12 * curl[qy][qx][1]) +
1897 (O13 * curl[qy][qx][2]);
1898 const double c2 = (O21 * curl[qy][qx][0]) + (O22 * curl[qy][qx][1]) +
1899 (O23 * curl[qy][qx][2]);
1900 const double c3 = (O31 * curl[qy][qx][0]) + (O32 * curl[qy][qx][1]) +
1901 (O33 * curl[qy][qx][2]);
1903 const double wcx = sBc[dx][qx];
1904 const double wDx = sGc[dx][qx];
1908 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
1909 // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
1910 const double wx = sBo[dx][qx];
1911 dxyz1 += (wx * c2 * wcy * wcDz) - (wx * c3 * wcDy * wcz);
1914 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
1915 // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
1916 dxyz2 += (-wy * c1 * wcx * wcDz) + (wy * c3 * wDx * wcz);
1918 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
1919 // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
1920 dxyz3 += (wcDy * wz * c1 * wcx) - (wcy * wz * c2 * wDx);
1929 MFEM_FOREACH_THREAD(dz,z,D1D)
1931 MFEM_FOREACH_THREAD(dy,y,D1D)
1933 MFEM_FOREACH_THREAD(dx,x,D1D)
1937 Y(dx + ((dy + (dz * D1D)) * (D1D-1)), e) += dxyz1;
1941 Y(dx + ((dy + (dz * (D1D-1))) * D1D) + ((D1D-1)*D1D*D1D), e) += dxyz2;
1945 Y(dx + ((dy + (dz * D1D)) * D1D) + (2*(D1D-1)*D1D*D1D), e) += dxyz3;
1951 }; // end of element loop
1953 auto host_kernel = [&] MFEM_LAMBDA (int)
1955 MFEM_ABORT_KERNEL("This kernel should only be used on GPU.
");
1958 ForallWrap<3>(true, NE, device_kernel, host_kernel, Q1D, Q1D, Q1D);
1961 static void PACurlL2Apply2D(const int D1D,
1965 const Array<double> &bo,
1966 const Array<double> &bot,
1967 const Array<double> &bt,
1968 const Array<double> &gc,
1969 const Vector &pa_data,
1970 const Vector &x, // trial = H(curl)
1971 Vector &y) // test = L2 or H1
1973 constexpr static int VDIM = 2;
1974 constexpr static int MAX_D1D = HCURL_MAX_D1D;
1975 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
1976 const int H1 = (D1Dtest == D1D);
1978 MFEM_VERIFY(y.Size() == NE*D1Dtest*D1Dtest, "Test vector of wrong
dimension");
1980 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
1981 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
1982 auto Bt = Reshape(bt.Read(), D1D, Q1D);
1983 auto Gc = Reshape(gc.Read(), Q1D, D1D);
1984 auto op = Reshape(pa_data.Read(), Q1D, Q1D, NE);
1985 auto X = Reshape(x.Read(), 2*(D1D-1)*D1D, NE);
1986 auto Y = Reshape(y.ReadWrite(), D1Dtest, D1Dtest, NE);
1990 double curl[MAX_Q1D][MAX_Q1D];
1992 // curl[qy][qx] will be computed as du_y/dx - du_x/dy
1994 for (int qy = 0; qy < Q1D; ++qy)
1996 for (int qx = 0; qx < Q1D; ++qx)
2004 for (int c = 0; c < VDIM; ++c) // loop over x, y components
2006 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2007 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2009 for (int dy = 0; dy < D1Dy; ++dy)
2011 double gradX[MAX_Q1D];
2012 for (int qx = 0; qx < Q1D; ++qx)
2017 for (int dx = 0; dx < D1Dx; ++dx)
2019 const double t = X(dx + (dy * D1Dx) + osc, e);
2020 for (int qx = 0; qx < Q1D; ++qx)
2022 gradX[qx] += t * ((c == 0) ? Bo(qx,dx) : Gc(qx,dx));
2026 for (int qy = 0; qy < Q1D; ++qy)
2028 const double wy = (c == 0) ? -Gc(qy,dy) : Bo(qy,dy);
2029 for (int qx = 0; qx < Q1D; ++qx)
2031 curl[qy][qx] += gradX[qx] * wy;
2037 } // loop (c) over components
2039 // Apply D operator.
2040 for (int qy = 0; qy < Q1D; ++qy)
2042 for (int qx = 0; qx < Q1D; ++qx)
2044 curl[qy][qx] *= op(qx,qy,e);
2048 for (int qy = 0; qy < Q1D; ++qy)
2050 double sol_x[MAX_D1D];
2051 for (int dx = 0; dx < D1Dtest; ++dx)
2055 for (int qx = 0; qx < Q1D; ++qx)
2057 const double s = curl[qy][qx];
2058 for (int dx = 0; dx < D1Dtest; ++dx)
2060 sol_x[dx] += s * ((H1 == 1) ? Bt(dx,qx) : Bot(dx,qx));
2063 for (int dy = 0; dy < D1Dtest; ++dy)
2065 const double wy = (H1 == 1) ? Bt(dy,qy) : Bot(dy,qy);
2067 for (int dx = 0; dx < D1Dtest; ++dx)
2069 Y(dx,dy,e) += sol_x[dx] * wy;
2073 }); // end of element loop
2076 static void PACurlL2ApplyTranspose2D(const int D1D,
2080 const Array<double> &bo,
2081 const Array<double> &bot,
2082 const Array<double> &b,
2083 const Array<double> &gct,
2084 const Vector &pa_data,
2085 const Vector &x, // trial = H(curl)
2086 Vector &y) // test = L2 or H1
2088 constexpr static int VDIM = 2;
2089 constexpr static int MAX_D1D = HCURL_MAX_D1D;
2090 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
2091 const int H1 = (D1Dtest == D1D);
2093 MFEM_VERIFY(x.Size() == NE*D1Dtest*D1Dtest, "Test vector of wrong
dimension");
2095 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
2096 auto B = Reshape(b.Read(), Q1D, D1D);
2097 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
2098 auto Gct = Reshape(gct.Read(), D1D, Q1D);
2099 auto op = Reshape(pa_data.Read(), Q1D, Q1D, NE);
2100 auto X = Reshape(x.Read(), D1Dtest, D1Dtest, NE);
2101 auto Y = Reshape(y.ReadWrite(), 2*(D1D-1)*D1D, NE);
2105 double mass[MAX_Q1D][MAX_Q1D];
2107 // Zero-order term in L2 or H1 test space
2109 for (int qy = 0; qy < Q1D; ++qy)
2111 for (int qx = 0; qx < Q1D; ++qx)
2117 for (int dy = 0; dy < D1Dtest; ++dy)
2119 double sol_x[MAX_Q1D];
2120 for (int qy = 0; qy < Q1D; ++qy)
2124 for (int dx = 0; dx < D1Dtest; ++dx)
2126 const double s = X(dx,dy,e);
2127 for (int qx = 0; qx < Q1D; ++qx)
2129 sol_x[qx] += s * ((H1 == 1) ? B(qx,dx) : Bo(qx,dx));
2132 for (int qy = 0; qy < Q1D; ++qy)
2134 const double d2q = (H1 == 1) ? B(qy,dy) : Bo(qy,dy);
2135 for (int qx = 0; qx < Q1D; ++qx)
2137 mass[qy][qx] += d2q * sol_x[qx];
2142 // Apply D operator.
2143 for (int qy = 0; qy < Q1D; ++qy)
2145 for (int qx = 0; qx < Q1D; ++qx)
2147 mass[qy][qx] *= op(qx,qy,e);
2151 for (int qy = 0; qy < Q1D; ++qy)
2155 for (int c = 0; c < VDIM; ++c) // loop over x, y components
2157 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2158 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2160 double gradX[MAX_D1D];
2161 for (int dx = 0; dx < D1Dx; ++dx)
2165 for (int qx = 0; qx < Q1D; ++qx)
2167 for (int dx = 0; dx < D1Dx; ++dx)
2169 gradX[dx] += mass[qy][qx] * ((c == 0) ? Bot(dx,qx) : Gct(dx,qx));
2172 for (int dy = 0; dy < D1Dy; ++dy)
2174 const double wy = (c == 0) ? -Gct(dy,qy) : Bot(dy,qy);
2176 for (int dx = 0; dx < D1Dx; ++dx)
2178 Y(dx + (dy * D1Dx) + osc, e) += gradX[dx] * wy;
2185 }); // end of element loop
2188 void CurlCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
2192 if (Device::Allows(Backend::DEVICE_MASK))
2194 const int ID = (dofs1D << 4) | quad1D;
2197 case 0x23: return SmemPACurlCurlApply3D<2,3>(dofs1D, quad1D, symmetric, ne,
2198 mapsO->B, mapsC->B, mapsO->Bt,
2199 mapsC->Bt, mapsC->G, mapsC->Gt, pa_data, x, y);
2200 case 0x34: return SmemPACurlCurlApply3D<3,4>(dofs1D, quad1D, symmetric, ne,
2201 mapsO->B, mapsC->B, mapsO->Bt,
2202 mapsC->Bt, mapsC->G, mapsC->Gt, pa_data, x, y);
2203 case 0x45: return SmemPACurlCurlApply3D<4,5>(dofs1D, quad1D, symmetric, ne,
2205 mapsC->B, mapsO->Bt,
2206 mapsC->Bt, mapsC->G, mapsC->Gt, pa_data, x, y);
2207 case 0x56: return SmemPACurlCurlApply3D<5,6>(dofs1D, quad1D, symmetric, ne,
2208 mapsO->B, mapsC->B, mapsO->Bt,
2209 mapsC->Bt, mapsC->G, mapsC->Gt, pa_data, x, y);
2210 default: return SmemPACurlCurlApply3D(dofs1D, quad1D, symmetric, ne, mapsO->B,
2211 mapsC->B, mapsO->Bt, mapsC->Bt,
2212 mapsC->G, mapsC->Gt, pa_data, x, y);
2216 PACurlCurlApply3D(dofs1D, quad1D, symmetric, ne, mapsO->B, mapsC->B, mapsO->Bt,
2217 mapsC->Bt, mapsC->G, mapsC->Gt, pa_data, x, y);
2221 PACurlCurlApply2D(dofs1D, quad1D, ne, mapsO->B, mapsO->Bt,
2222 mapsC->G, mapsC->Gt, pa_data, x, y);
2230 static void PACurlCurlAssembleDiagonal2D(const int D1D,
2233 const Array<double> &bo,
2234 const Array<double> &gc,
2235 const Vector &pa_data,
2238 constexpr static int VDIM = 2;
2239 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
2241 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
2242 auto Gc = Reshape(gc.Read(), Q1D, D1D);
2243 auto op = Reshape(pa_data.Read(), Q1D, Q1D, NE);
2244 auto D = Reshape(diag.ReadWrite(), 2*(D1D-1)*D1D, NE);
2250 for (int c = 0; c < VDIM; ++c) // loop over x, y components
2252 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2253 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2257 for (int dy = 0; dy < D1Dy; ++dy)
2259 for (int qx = 0; qx < Q1D; ++qx)
2262 for (int qy = 0; qy < Q1D; ++qy)
2264 const double wy = (c == 1) ? Bo(qy,dy) : -Gc(qy,dy);
2265 t[qx] += wy * wy * op(qx,qy,e);
2269 for (int dx = 0; dx < D1Dx; ++dx)
2271 for (int qx = 0; qx < Q1D; ++qx)
2273 const double wx = ((c == 0) ? Bo(qx,dx) : Gc(qx,dx));
2274 D(dx + (dy * D1Dx) + osc, e) += t[qx] * wx * wx;
2281 }); // end of element loop
2284 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
2285 static void PACurlCurlAssembleDiagonal3D(const int D1D,
2287 const bool symmetric,
2289 const Array<double> &bo,
2290 const Array<double> &bc,
2291 const Array<double> &go,
2292 const Array<double> &gc,
2293 const Vector &pa_data,
2296 constexpr static int VDIM = 3;
2297 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
2298 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
2300 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
2301 auto Bc = Reshape(bc.Read(), Q1D, D1D);
2302 auto Go = Reshape(go.Read(), Q1D, D1D-1);
2303 auto Gc = Reshape(gc.Read(), Q1D, D1D);
2304 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, (symmetric ? 6 : 9), NE);
2305 auto D = Reshape(diag.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
2307 const int s = symmetric ? 6 : 9;
2311 const int i21 = symmetric ? i12 : 3;
2312 const int i22 = symmetric ? 3 : 4;
2313 const int i23 = symmetric ? 4 : 5;
2314 const int i31 = symmetric ? i13 : 6;
2315 const int i32 = symmetric ? i23 : 7;
2316 const int i33 = symmetric ? 5 : 8;
2320 // Using (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u} (p. 78 of Monk), we get
2321 // (\nabla\times u) \cdot (\nabla\times u) = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{\nabla}\times\hat{u}
2322 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
2323 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
2324 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
2326 // For each c, we will keep 9 arrays for derivatives multiplied by the 9 entries of the 3x3 matrix (dF^T C dF),
2327 // which may be non-symmetric depending on a possibly non-symmetric matrix coefficient.
2331 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
2333 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
2334 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2335 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2337 double zt[MAX_Q1D][MAX_Q1D][MAX_D1D][9][3];
2340 for (int qx = 0; qx < Q1D; ++qx)
2342 for (int qy = 0; qy < Q1D; ++qy)
2344 for (int dz = 0; dz < D1Dz; ++dz)
2346 for (int i=0; i<s; ++i)
2348 for (int d=0; d<3; ++d)
2350 zt[qx][qy][dz][i][d] = 0.0;
2354 for (int qz = 0; qz < Q1D; ++qz)
2356 const double wz = ((c == 2) ? Bo(qz,dz) : Bc(qz,dz));
2357 const double wDz = ((c == 2) ? Go(qz,dz) : Gc(qz,dz));
2359 for (int i=0; i<s; ++i)
2361 zt[qx][qy][dz][i][0] += wz * wz * op(qx,qy,qz,i,e);
2362 zt[qx][qy][dz][i][1] += wDz * wz * op(qx,qy,qz,i,e);
2363 zt[qx][qy][dz][i][2] += wDz * wDz * op(qx,qy,qz,i,e);
2368 } // end of z contraction
2370 double yt[MAX_Q1D][MAX_D1D][MAX_D1D][9][3][3];
2373 for (int qx = 0; qx < Q1D; ++qx)
2375 for (int dz = 0; dz < D1Dz; ++dz)
2377 for (int dy = 0; dy < D1Dy; ++dy)
2379 for (int i=0; i<s; ++i)
2381 for (int d=0; d<3; ++d)
2382 for (int j=0; j<3; ++j)
2384 yt[qx][dy][dz][i][d][j] = 0.0;
2388 for (int qy = 0; qy < Q1D; ++qy)
2390 const double wy = ((c == 1) ? Bo(qy,dy) : Bc(qy,dy));
2391 const double wDy = ((c == 1) ? Go(qy,dy) : Gc(qy,dy));
2393 for (int i=0; i<s; ++i)
2395 for (int d=0; d<3; ++d)
2397 yt[qx][dy][dz][i][d][0] += wy * wy * zt[qx][qy][dz][i][d];
2398 yt[qx][dy][dz][i][d][1] += wDy * wy * zt[qx][qy][dz][i][d];
2399 yt[qx][dy][dz][i][d][2] += wDy * wDy * zt[qx][qy][dz][i][d];
2405 } // end of y contraction
2408 for (int dz = 0; dz < D1Dz; ++dz)
2410 for (int dy = 0; dy < D1Dy; ++dy)
2412 for (int dx = 0; dx < D1Dx; ++dx)
2414 for (int qx = 0; qx < Q1D; ++qx)
2416 const double wx = ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
2417 const double wDx = ((c == 0) ? Go(qx,dx) : Gc(qx,dx));
2419 // Using (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u} (p. 78 of Monk), we get
2420 // (\nabla\times u) \cdot (\nabla\times u) = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{\nabla}\times\hat{u}
2421 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
2422 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
2423 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
2426 const double O11 = op(q,0,e);
2427 const double O12 = op(q,1,e);
2428 const double O13 = op(q,2,e);
2429 const double O22 = op(q,3,e);
2430 const double O23 = op(q,4,e);
2431 const double O33 = op(q,5,e);
2436 // (u_0)_{x_2} (O22 (u_0)_{x_2} - O23 (u_0)_{x_1}) - (u_0)_{x_1} (O32 (u_0)_{x_2} - O33 (u_0)_{x_1})
2437 const double sumy = yt[qx][dy][dz][i22][2][0] - yt[qx][dy][dz][i23][1][1]
2438 - yt[qx][dy][dz][i32][1][1] + yt[qx][dy][dz][i33][0][2];
2440 D(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += sumy * wx * wx;
2444 // (u_1)_{x_2} (O11 (u_1)_{x_2} - O13 (u_1)_{x_0}) + (u_1)_{x_0} (-O31 (u_1)_{x_2} + O33 (u_1)_{x_0})
2445 const double d = (yt[qx][dy][dz][i11][2][0] * wx * wx)
2446 - ((yt[qx][dy][dz][i13][1][0] + yt[qx][dy][dz][i31][1][0]) * wDx * wx)
2447 + (yt[qx][dy][dz][i33][0][0] * wDx * wDx);
2449 D(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += d;
2453 // (u_2)_{x_1} (O11 (u_2)_{x_1} - O12 (u_2)_{x_0}) - (u_2)_{x_0} (O21 (u_2)_{x_1} - O22 (u_2)_{x_0})
2454 const double d = (yt[qx][dy][dz][i11][0][2] * wx * wx)
2455 - ((yt[qx][dy][dz][i12][0][1] + yt[qx][dy][dz][i21][0][1]) * wDx * wx)
2456 + (yt[qx][dy][dz][i22][0][0] * wDx * wDx);
2458 D(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += d;
2463 } // end of x contraction
2465 osc += D1Dx * D1Dy * D1Dz;
2467 }); // end of element loop
2470 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
2471 static void SmemPACurlCurlAssembleDiagonal3D(const int D1D,
2473 const bool symmetric,
2475 const Array<double> &bo,
2476 const Array<double> &bc,
2477 const Array<double> &go,
2478 const Array<double> &gc,
2479 const Vector &pa_data,
2482 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
2483 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
2485 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
2486 auto Bc = Reshape(bc.Read(), Q1D, D1D);
2487 auto Go = Reshape(go.Read(), Q1D, D1D-1);
2488 auto Gc = Reshape(gc.Read(), Q1D, D1D);
2489 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, (symmetric ? 6 : 9), NE);
2490 auto D = Reshape(diag.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
2492 const int s = symmetric ? 6 : 9;
2496 const int i21 = symmetric ? i12 : 3;
2497 const int i22 = symmetric ? 3 : 4;
2498 const int i23 = symmetric ? 4 : 5;
2499 const int i31 = symmetric ? i13 : 6;
2500 const int i32 = symmetric ? i23 : 7;
2501 const int i33 = symmetric ? 5 : 8;
2503 MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
2505 // Using (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u} (p. 78 of Monk), we get
2506 // (\nabla\times u) \cdot (\nabla\times u) = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{\nabla}\times\hat{u}
2507 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
2508 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
2509 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
2511 constexpr int VDIM = 3;
2513 MFEM_SHARED double sBo[MAX_Q1D][MAX_D1D];
2514 MFEM_SHARED double sBc[MAX_Q1D][MAX_D1D];
2515 MFEM_SHARED double sGo[MAX_Q1D][MAX_D1D];
2516 MFEM_SHARED double sGc[MAX_Q1D][MAX_D1D];
2519 MFEM_SHARED double sop[9][MAX_Q1D][MAX_Q1D];
2521 MFEM_FOREACH_THREAD(qx,x,Q1D)
2523 MFEM_FOREACH_THREAD(qy,y,Q1D)
2525 MFEM_FOREACH_THREAD(qz,z,Q1D)
2527 for (int i=0; i<s; ++i)
2529 ope[i] = op(qx,qy,qz,i,e);
2535 const int tidx = MFEM_THREAD_ID(x);
2536 const int tidy = MFEM_THREAD_ID(y);
2537 const int tidz = MFEM_THREAD_ID(z);
2541 MFEM_FOREACH_THREAD(d,y,D1D)
2543 MFEM_FOREACH_THREAD(q,x,Q1D)
2545 sBc[q][d] = Bc(q,d);
2546 sGc[q][d] = Gc(q,d);
2549 sBo[q][d] = Bo(q,d);
2550 sGo[q][d] = Go(q,d);
2558 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
2560 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
2561 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2562 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2566 for (int qz=0; qz < Q1D; ++qz)
2570 for (int i=0; i<s; ++i)
2572 sop[i][tidx][tidy] = ope[i];
2578 MFEM_FOREACH_THREAD(dz,z,D1Dz)
2580 const double wz = ((c == 2) ? sBo[qz][dz] : sBc[qz][dz]);
2581 const double wDz = ((c == 2) ? sGo[qz][dz] : sGc[qz][dz]);
2583 MFEM_FOREACH_THREAD(dy,y,D1Dy)
2585 MFEM_FOREACH_THREAD(dx,x,D1Dx)
2587 for (int qy = 0; qy < Q1D; ++qy)
2589 const double wy = ((c == 1) ? sBo[qy][dy] : sBc[qy][dy]);
2590 const double wDy = ((c == 1) ? sGo[qy][dy] : sGc[qy][dy]);
2592 for (int qx = 0; qx < Q1D; ++qx)
2594 const double wx = ((c == 0) ? sBo[qx][dx] : sBc[qx][dx]);
2595 const double wDx = ((c == 0) ? sGo[qx][dx] : sGc[qx][dx]);
2599 // (u_0)_{x_2} (O22 (u_0)_{x_2} - O23 (u_0)_{x_1}) - (u_0)_{x_1} (O32 (u_0)_{x_2} - O33 (u_0)_{x_1})
2601 // (u_0)_{x_2} O22 (u_0)_{x_2}
2602 dxyz += sop[i22][qx][qy] * wx * wx * wy * wy * wDz * wDz;
2604 // -(u_0)_{x_2} O23 (u_0)_{x_1} - (u_0)_{x_1} O32 (u_0)_{x_2}
2605 dxyz += -(sop[i23][qx][qy] + sop[i32][qx][qy]) * wx * wx * wDy * wy * wDz * wz;
2607 // (u_0)_{x_1} O33 (u_0)_{x_1}
2608 dxyz += sop[i33][qx][qy] * wx * wx * wDy * wDy * wz * wz;
2612 // (u_1)_{x_2} (O11 (u_1)_{x_2} - O13 (u_1)_{x_0}) + (u_1)_{x_0} (-O31 (u_1)_{x_2} + O33 (u_1)_{x_0})
2614 // (u_1)_{x_2} O11 (u_1)_{x_2}
2615 dxyz += sop[i11][qx][qy] * wx * wx * wy * wy * wDz * wDz;
2617 // -(u_1)_{x_2} O13 (u_1)_{x_0} - (u_1)_{x_0} O31 (u_1)_{x_2}
2618 dxyz += -(sop[i13][qx][qy] + sop[i31][qx][qy]) * wDx * wx * wy * wy * wDz * wz;
2620 // (u_1)_{x_0} O33 (u_1)_{x_0})
2621 dxyz += sop[i33][qx][qy] * wDx * wDx * wy * wy * wz * wz;
2625 // (u_2)_{x_1} (O11 (u_2)_{x_1} - O12 (u_2)_{x_0}) - (u_2)_{x_0} (O21 (u_2)_{x_1} - O22 (u_2)_{x_0})
2627 // (u_2)_{x_1} O11 (u_2)_{x_1}
2628 dxyz += sop[i11][qx][qy] * wx * wx * wDy * wDy * wz * wz;
2630 // -(u_2)_{x_1} O12 (u_2)_{x_0} - (u_2)_{x_0} O21 (u_2)_{x_1}
2631 dxyz += -(sop[i12][qx][qy] + sop[i21][qx][qy]) * wDx * wx * wDy * wy * wz * wz;
2633 // (u_2)_{x_0} O22 (u_2)_{x_0}
2634 dxyz += sop[i22][qx][qy] * wDx * wDx * wy * wy * wz * wz;
2645 MFEM_FOREACH_THREAD(dz,z,D1Dz)
2647 MFEM_FOREACH_THREAD(dy,y,D1Dy)
2649 MFEM_FOREACH_THREAD(dx,x,D1Dx)
2651 D(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += dxyz;
2656 osc += D1Dx * D1Dy * D1Dz;
2658 }); // end of element loop
2661 void CurlCurlIntegrator::AssembleDiagonalPA(Vector& diag)
2665 if (Device::Allows(Backend::DEVICE_MASK))
2667 const int ID = (dofs1D << 4) | quad1D;
2670 case 0x23: return SmemPACurlCurlAssembleDiagonal3D<2,3>(dofs1D, quad1D,
2675 case 0x34: return SmemPACurlCurlAssembleDiagonal3D<3,4>(dofs1D, quad1D,
2680 case 0x45: return SmemPACurlCurlAssembleDiagonal3D<4,5>(dofs1D, quad1D,
2685 case 0x56: return SmemPACurlCurlAssembleDiagonal3D<5,6>(dofs1D, quad1D,
2690 default: return SmemPACurlCurlAssembleDiagonal3D(dofs1D, quad1D, symmetric, ne,
2697 PACurlCurlAssembleDiagonal3D(dofs1D, quad1D, symmetric, ne,
2704 PACurlCurlAssembleDiagonal2D(dofs1D, quad1D, ne,
2705 mapsO->B, mapsC->G, pa_data, diag);
2713 // Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
2714 // integrated against H(curl) test functions corresponding to y.
2715 void PAHcurlH1Apply3D(const int D1D,
2718 const Array<double> &bc,
2719 const Array<double> &gc,
2720 const Array<double> &bot,
2721 const Array<double> &bct,
2722 const Vector &pa_data,
2726 constexpr static int MAX_D1D = HCURL_MAX_D1D;
2727 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
2729 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
2730 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
2732 constexpr static int VDIM = 3;
2734 auto Bc = Reshape(bc.Read(), Q1D, D1D);
2735 auto Gc = Reshape(gc.Read(), Q1D, D1D);
2736 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
2737 auto Bct = Reshape(bct.Read(), D1D, Q1D);
2738 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
2739 auto X = Reshape(x.Read(), D1D, D1D, D1D, NE);
2740 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
2744 double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
2746 for (int qz = 0; qz < Q1D; ++qz)
2748 for (int qy = 0; qy < Q1D; ++qy)
2750 for (int qx = 0; qx < Q1D; ++qx)
2752 for (int c = 0; c < VDIM; ++c)
2754 mass[qz][qy][qx][c] = 0.0;
2760 for (int dz = 0; dz < D1D; ++dz)
2762 double gradXY[MAX_Q1D][MAX_Q1D][3];
2763 for (int qy = 0; qy < Q1D; ++qy)
2765 for (int qx = 0; qx < Q1D; ++qx)
2767 gradXY[qy][qx][0] = 0.0;
2768 gradXY[qy][qx][1] = 0.0;
2769 gradXY[qy][qx][2] = 0.0;
2772 for (int dy = 0; dy < D1D; ++dy)
2774 double gradX[MAX_Q1D][2];
2775 for (int qx = 0; qx < Q1D; ++qx)
2780 for (int dx = 0; dx < D1D; ++dx)
2782 const double s = X(dx,dy,dz,e);
2783 for (int qx = 0; qx < Q1D; ++qx)
2785 gradX[qx][0] += s * Bc(qx,dx);
2786 gradX[qx][1] += s * Gc(qx,dx);
2789 for (int qy = 0; qy < Q1D; ++qy)
2791 const double wy = Bc(qy,dy);
2792 const double wDy = Gc(qy,dy);
2793 for (int qx = 0; qx < Q1D; ++qx)
2795 const double wx = gradX[qx][0];
2796 const double wDx = gradX[qx][1];
2797 gradXY[qy][qx][0] += wDx * wy;
2798 gradXY[qy][qx][1] += wx * wDy;
2799 gradXY[qy][qx][2] += wx * wy;
2803 for (int qz = 0; qz < Q1D; ++qz)
2805 const double wz = Bc(qz,dz);
2806 const double wDz = Gc(qz,dz);
2807 for (int qy = 0; qy < Q1D; ++qy)
2809 for (int qx = 0; qx < Q1D; ++qx)
2811 mass[qz][qy][qx][0] += gradXY[qy][qx][0] * wz;
2812 mass[qz][qy][qx][1] += gradXY[qy][qx][1] * wz;
2813 mass[qz][qy][qx][2] += gradXY[qy][qx][2] * wDz;
2819 // Apply D operator.
2820 for (int qz = 0; qz < Q1D; ++qz)
2822 for (int qy = 0; qy < Q1D; ++qy)
2824 for (int qx = 0; qx < Q1D; ++qx)
2826 const double O11 = op(qx,qy,qz,0,e);
2827 const double O12 = op(qx,qy,qz,1,e);
2828 const double O13 = op(qx,qy,qz,2,e);
2829 const double O22 = op(qx,qy,qz,3,e);
2830 const double O23 = op(qx,qy,qz,4,e);
2831 const double O33 = op(qx,qy,qz,5,e);
2832 const double massX = mass[qz][qy][qx][0];
2833 const double massY = mass[qz][qy][qx][1];
2834 const double massZ = mass[qz][qy][qx][2];
2835 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
2836 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
2837 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
2842 for (int qz = 0; qz < Q1D; ++qz)
2844 double massXY[MAX_D1D][MAX_D1D];
2848 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
2850 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
2851 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2852 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2854 for (int dy = 0; dy < D1Dy; ++dy)
2856 for (int dx = 0; dx < D1Dx; ++dx)
2858 massXY[dy][dx] = 0.0;
2861 for (int qy = 0; qy < Q1D; ++qy)
2863 double massX[MAX_D1D];
2864 for (int dx = 0; dx < D1Dx; ++dx)
2868 for (int qx = 0; qx < Q1D; ++qx)
2870 for (int dx = 0; dx < D1Dx; ++dx)
2872 massX[dx] += mass[qz][qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
2875 for (int dy = 0; dy < D1Dy; ++dy)
2877 const double wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
2878 for (int dx = 0; dx < D1Dx; ++dx)
2880 massXY[dy][dx] += massX[dx] * wy;
2885 for (int dz = 0; dz < D1Dz; ++dz)
2887 const double wz = (c == 2) ? Bot(dz,qz) : Bct(dz,qz);
2888 for (int dy = 0; dy < D1Dy; ++dy)
2890 for (int dx = 0; dx < D1Dx; ++dx)
2892 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += massXY[dy][dx] * wz;
2897 osc += D1Dx * D1Dy * D1Dz;
2900 }); // end of element loop
2903 // Apply to x corresponding to DOFs in H(curl), integrated
2904 // against gradients of H^1 functions corresponding to y.
2905 void PAHcurlH1ApplyTranspose3D(const int D1D,
2908 const Array<double> &bc,
2909 const Array<double> &bo,
2910 const Array<double> &bct,
2911 const Array<double> &gct,
2912 const Vector &pa_data,
2916 constexpr static int MAX_D1D = HCURL_MAX_D1D;
2917 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
2919 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
2920 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
2922 constexpr static int VDIM = 3;
2924 auto Bc = Reshape(bc.Read(), Q1D, D1D);
2925 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
2926 auto Bt = Reshape(bct.Read(), D1D, Q1D);
2927 auto Gt = Reshape(gct.Read(), D1D, Q1D);
2928 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
2929 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
2930 auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, NE);
2934 double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
2936 for (int qz = 0; qz < Q1D; ++qz)
2938 for (int qy = 0; qy < Q1D; ++qy)
2940 for (int qx = 0; qx < Q1D; ++qx)
2942 for (int c = 0; c < VDIM; ++c)
2944 mass[qz][qy][qx][c] = 0.0;
2952 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
2954 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
2955 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
2956 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
2958 for (int dz = 0; dz < D1Dz; ++dz)
2960 double massXY[MAX_Q1D][MAX_Q1D];
2961 for (int qy = 0; qy < Q1D; ++qy)
2963 for (int qx = 0; qx < Q1D; ++qx)
2965 massXY[qy][qx] = 0.0;
2969 for (int dy = 0; dy < D1Dy; ++dy)
2971 double massX[MAX_Q1D];
2972 for (int qx = 0; qx < Q1D; ++qx)
2977 for (int dx = 0; dx < D1Dx; ++dx)
2979 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
2980 for (int qx = 0; qx < Q1D; ++qx)
2982 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
2986 for (int qy = 0; qy < Q1D; ++qy)
2988 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
2989 for (int qx = 0; qx < Q1D; ++qx)
2991 const double wx = massX[qx];
2992 massXY[qy][qx] += wx * wy;
2997 for (int qz = 0; qz < Q1D; ++qz)
2999 const double wz = (c == 2) ? Bo(qz,dz) : Bc(qz,dz);
3000 for (int qy = 0; qy < Q1D; ++qy)
3002 for (int qx = 0; qx < Q1D; ++qx)
3004 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
3010 osc += D1Dx * D1Dy * D1Dz;
3011 } // loop (c) over components
3013 // Apply D operator.
3014 for (int qz = 0; qz < Q1D; ++qz)
3016 for (int qy = 0; qy < Q1D; ++qy)
3018 for (int qx = 0; qx < Q1D; ++qx)
3020 const double O11 = op(qx,qy,qz,0,e);
3021 const double O12 = op(qx,qy,qz,1,e);
3022 const double O13 = op(qx,qy,qz,2,e);
3023 const double O22 = op(qx,qy,qz,3,e);
3024 const double O23 = op(qx,qy,qz,4,e);
3025 const double O33 = op(qx,qy,qz,5,e);
3026 const double massX = mass[qz][qy][qx][0];
3027 const double massY = mass[qz][qy][qx][1];
3028 const double massZ = mass[qz][qy][qx][2];
3029 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
3030 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
3031 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
3036 for (int qz = 0; qz < Q1D; ++qz)
3038 double gradXY[MAX_D1D][MAX_D1D][3];
3039 for (int dy = 0; dy < D1D; ++dy)
3041 for (int dx = 0; dx < D1D; ++dx)
3043 gradXY[dy][dx][0] = 0;
3044 gradXY[dy][dx][1] = 0;
3045 gradXY[dy][dx][2] = 0;
3048 for (int qy = 0; qy < Q1D; ++qy)
3050 double gradX[MAX_D1D][3];
3051 for (int dx = 0; dx < D1D; ++dx)
3057 for (int qx = 0; qx < Q1D; ++qx)
3059 const double gX = mass[qz][qy][qx][0];
3060 const double gY = mass[qz][qy][qx][1];
3061 const double gZ = mass[qz][qy][qx][2];
3062 for (int dx = 0; dx < D1D; ++dx)
3064 const double wx = Bt(dx,qx);
3065 const double wDx = Gt(dx,qx);
3066 gradX[dx][0] += gX * wDx;
3067 gradX[dx][1] += gY * wx;
3068 gradX[dx][2] += gZ * wx;
3071 for (int dy = 0; dy < D1D; ++dy)
3073 const double wy = Bt(dy,qy);
3074 const double wDy = Gt(dy,qy);
3075 for (int dx = 0; dx < D1D; ++dx)
3077 gradXY[dy][dx][0] += gradX[dx][0] * wy;
3078 gradXY[dy][dx][1] += gradX[dx][1] * wDy;
3079 gradXY[dy][dx][2] += gradX[dx][2] * wy;
3083 for (int dz = 0; dz < D1D; ++dz)
3085 const double wz = Bt(dz,qz);
3086 const double wDz = Gt(dz,qz);
3087 for (int dy = 0; dy < D1D; ++dy)
3089 for (int dx = 0; dx < D1D; ++dx)
3092 ((gradXY[dy][dx][0] * wz) +
3093 (gradXY[dy][dx][1] * wz) +
3094 (gradXY[dy][dx][2] * wDz));
3099 }); // end of element loop
3102 // Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
3103 // integrated against H(curl) test functions corresponding to y.
3104 void PAHcurlH1Apply2D(const int D1D,
3107 const Array<double> &bc,
3108 const Array<double> &gc,
3109 const Array<double> &bot,
3110 const Array<double> &bct,
3111 const Vector &pa_data,
3115 constexpr static int VDIM = 2;
3116 constexpr static int MAX_D1D = HCURL_MAX_D1D;
3117 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
3119 auto Bc = Reshape(bc.Read(), Q1D, D1D);
3120 auto Gc = Reshape(gc.Read(), Q1D, D1D);
3121 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
3122 auto Bct = Reshape(bct.Read(), D1D, Q1D);
3123 auto op = Reshape(pa_data.Read(), Q1D, Q1D, 3, NE);
3124 auto X = Reshape(x.Read(), D1D, D1D, NE);
3125 auto Y = Reshape(y.ReadWrite(), 2*(D1D-1)*D1D, NE);
3129 double mass[MAX_Q1D][MAX_Q1D][VDIM];
3131 for (int qy = 0; qy < Q1D; ++qy)
3133 for (int qx = 0; qx < Q1D; ++qx)
3135 for (int c = 0; c < VDIM; ++c)
3137 mass[qy][qx][c] = 0.0;
3142 for (int dy = 0; dy < D1D; ++dy)
3144 double gradX[MAX_Q1D][2];
3145 for (int qx = 0; qx < Q1D; ++qx)
3150 for (int dx = 0; dx < D1D; ++dx)
3152 const double s = X(dx,dy,e);
3153 for (int qx = 0; qx < Q1D; ++qx)
3155 gradX[qx][0] += s * Bc(qx,dx);
3156 gradX[qx][1] += s * Gc(qx,dx);
3159 for (int qy = 0; qy < Q1D; ++qy)
3161 const double wy = Bc(qy,dy);
3162 const double wDy = Gc(qy,dy);
3163 for (int qx = 0; qx < Q1D; ++qx)
3165 const double wx = gradX[qx][0];
3166 const double wDx = gradX[qx][1];
3167 mass[qy][qx][0] += wDx * wy;
3168 mass[qy][qx][1] += wx * wDy;
3173 // Apply D operator.
3174 for (int qy = 0; qy < Q1D; ++qy)
3176 for (int qx = 0; qx < Q1D; ++qx)
3178 const double O11 = op(qx,qy,0,e);
3179 const double O12 = op(qx,qy,1,e);
3180 const double O22 = op(qx,qy,2,e);
3181 const double massX = mass[qy][qx][0];
3182 const double massY = mass[qy][qx][1];
3183 mass[qy][qx][0] = (O11*massX)+(O12*massY);
3184 mass[qy][qx][1] = (O12*massX)+(O22*massY);
3188 for (int qy = 0; qy < Q1D; ++qy)
3192 for (int c = 0; c < VDIM; ++c) // loop over x, y components
3194 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
3195 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
3197 double massX[MAX_D1D];
3198 for (int dx = 0; dx < D1Dx; ++dx)
3202 for (int qx = 0; qx < Q1D; ++qx)
3204 for (int dx = 0; dx < D1Dx; ++dx)
3206 massX[dx] += mass[qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
3210 for (int dy = 0; dy < D1Dy; ++dy)
3212 const double wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
3214 for (int dx = 0; dx < D1Dx; ++dx)
3216 Y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
3223 }); // end of element loop
3226 // Apply to x corresponding to DOFs in H(curl), integrated
3227 // against gradients of H^1 functions corresponding to y.
3228 void PAHcurlH1ApplyTranspose2D(const int D1D,
3231 const Array<double> &bc,
3232 const Array<double> &bo,
3233 const Array<double> &bct,
3234 const Array<double> &gct,
3235 const Vector &pa_data,
3239 constexpr static int VDIM = 2;
3240 constexpr static int MAX_D1D = HCURL_MAX_D1D;
3241 constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
3243 auto Bc = Reshape(bc.Read(), Q1D, D1D);
3244 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
3245 auto Bt = Reshape(bct.Read(), D1D, Q1D);
3246 auto Gt = Reshape(gct.Read(), D1D, Q1D);
3247 auto op = Reshape(pa_data.Read(), Q1D, Q1D, 3, NE);
3248 auto X = Reshape(x.Read(), 2*(D1D-1)*D1D, NE);
3249 auto Y = Reshape(y.ReadWrite(), D1D, D1D, NE);
3253 double mass[MAX_Q1D][MAX_Q1D][VDIM];
3255 for (int qy = 0; qy < Q1D; ++qy)
3257 for (int qx = 0; qx < Q1D; ++qx)
3259 for (int c = 0; c < VDIM; ++c)
3261 mass[qy][qx][c] = 0.0;
3268 for (int c = 0; c < VDIM; ++c) // loop over x, y components
3270 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
3271 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
3273 for (int dy = 0; dy < D1Dy; ++dy)
3275 double massX[MAX_Q1D];
3276 for (int qx = 0; qx < Q1D; ++qx)
3281 for (int dx = 0; dx < D1Dx; ++dx)
3283 const double t = X(dx + (dy * D1Dx) + osc, e);
3284 for (int qx = 0; qx < Q1D; ++qx)
3286 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
3290 for (int qy = 0; qy < Q1D; ++qy)
3292 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
3293 for (int qx = 0; qx < Q1D; ++qx)
3295 mass[qy][qx][c] += massX[qx] * wy;
3301 } // loop (c) over components
3303 // Apply D operator.
3304 for (int qy = 0; qy < Q1D; ++qy)
3306 for (int qx = 0; qx < Q1D; ++qx)
3308 const double O11 = op(qx,qy,0,e);
3309 const double O12 = op(qx,qy,1,e);
3310 const double O22 = op(qx,qy,2,e);
3311 const double massX = mass[qy][qx][0];
3312 const double massY = mass[qy][qx][1];
3313 mass[qy][qx][0] = (O11*massX)+(O12*massY);
3314 mass[qy][qx][1] = (O12*massX)+(O22*massY);
3318 for (int qy = 0; qy < Q1D; ++qy)
3320 double gradX[MAX_D1D][2];
3321 for (int dx = 0; dx < D1D; ++dx)
3326 for (int qx = 0; qx < Q1D; ++qx)
3328 const double gX = mass[qy][qx][0];
3329 const double gY = mass[qy][qx][1];
3330 for (int dx = 0; dx < D1D; ++dx)
3332 const double wx = Bt(dx,qx);
3333 const double wDx = Gt(dx,qx);
3334 gradX[dx][0] += gX * wDx;
3335 gradX[dx][1] += gY * wx;
3338 for (int dy = 0; dy < D1D; ++dy)
3340 const double wy = Bt(dy,qy);
3341 const double wDy = Gt(dy,qy);
3342 for (int dx = 0; dx < D1D; ++dx)
3344 Y(dx,dy,e) += ((gradX[dx][0] * wy) + (gradX[dx][1] * wDy));
3348 }); // end of element loop
3351 // PA H(curl) Mass Assemble 3D kernel
3352 void PAHcurlL2Setup(const int NQ,
3355 const Array<double> &w,
3360 auto C = Reshape(coeff.Read(), coeffDim, NQ, NE);
3361 auto y = Reshape(op.Write(), coeffDim, NQ, NE);
3365 for (int q = 0; q < NQ; ++q)
3367 for (int c=0; c<coeffDim; ++c)
3369 y(c,q,e) = W[q] * C(c,q,e);
3375 void MixedScalarCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
3376 const FiniteElementSpace &test_fes)
3378 // Assumes tensor-product elements
3379 Mesh *mesh = trial_fes.GetMesh();
3380 const FiniteElement *fel = trial_fes.GetFE(0); // In H(curl)
3381 const FiniteElement *eltest = test_fes.GetFE(0); // In scalar space
3383 const VectorTensorFiniteElement *el =
3384 dynamic_cast<const VectorTensorFiniteElement*>(fel);
3387 if (el->GetDerivType() != mfem::FiniteElement::CURL)
3389 MFEM_ABORT("Unknown kernel.
");
3392 const IntegrationRule *ir
3393 = IntRule ? IntRule : &MassIntegrator::GetRule(*eltest, *eltest,
3394 *mesh->GetElementTransformation(0));
3396 const int dims = el->GetDim();
3397 MFEM_VERIFY(dims == 2, "");
3399 const int nq = ir->GetNPoints();
3400 dim = mesh->Dimension();
3401 MFEM_VERIFY(dim == 2, "");
3403 ne = test_fes.GetNE();
3404 mapsC = &el->GetDofToQuad(*ir, DofToQuad::TENSOR);
3405 mapsO = &el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
3406 dofs1D = mapsC->ndof;
3407 quad1D = mapsC->nqpt;
3409 MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
3411 if (el->GetOrder() == eltest->GetOrder())
3413 dofs1Dtest = dofs1D;
3417 dofs1Dtest = dofs1D - 1;
3420 pa_data.SetSize(nq * ne, Device::GetMemoryType());
3422 QuadratureSpace qs(*mesh, *ir);
3423 CoefficientVector coeff(Q, qs, CoefficientStorage::FULL);
3427 PACurlL2Setup2D(quad1D, ne, ir->GetWeights(), coeff, pa_data);
3435 void MixedScalarCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
3439 PACurlL2Apply2D(dofs1D, dofs1Dtest, quad1D, ne, mapsO->B, mapsO->Bt,
3440 mapsC->Bt, mapsC->G, pa_data, x, y);
3448 void MixedScalarCurlIntegrator::AddMultTransposePA(const Vector &x,
3453 PACurlL2ApplyTranspose2D(dofs1D, dofs1Dtest, quad1D, ne, mapsO->B, mapsO->Bt,
3454 mapsC->B, mapsC->Gt, pa_data, x, y);
3462 void MixedVectorCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
3463 const FiniteElementSpace &test_fes)
3465 // Assumes tensor-product elements, with vector test and trial spaces.
3466 Mesh *mesh = trial_fes.GetMesh();
3467 const FiniteElement *trial_fel = trial_fes.GetFE(0);
3468 const FiniteElement *test_fel = test_fes.GetFE(0);
3470 const VectorTensorFiniteElement *trial_el =
3471 dynamic_cast<const VectorTensorFiniteElement*>(trial_fel);
3474 const VectorTensorFiniteElement *test_el =
3475 dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
3478 const IntegrationRule *ir
3479 = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
3480 *mesh->GetElementTransformation(0));
3481 const int dims = trial_el->GetDim();
3482 MFEM_VERIFY(dims == 3, "");
3484 const int nq = ir->GetNPoints();
3485 dim = mesh->Dimension();
3486 MFEM_VERIFY(dim == 3, "");
3488 MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(), "");
3490 ne = trial_fes.GetNE();
3491 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
3492 mapsC = &trial_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
3493 mapsO = &trial_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
3494 mapsCtest = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
3495 mapsOtest = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
3496 dofs1D = mapsC->ndof;
3497 quad1D = mapsC->nqpt;
3498 dofs1Dtest = mapsCtest->ndof;
3500 MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
3502 testType = test_el->GetDerivType();
3503 trialType = trial_el->GetDerivType();
3505 const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
3506 coeffDim = (DQ ? 3 : 1);
3508 const bool curlSpaces = (testType == mfem::FiniteElement::CURL &&
3509 trialType == mfem::FiniteElement::CURL);
3511 const int ndata = curlSpaces ? (coeffDim == 1 ? 1 : 9) : symmDims;
3512 pa_data.SetSize(ndata * nq * ne, Device::GetMemoryType());
3514 QuadratureSpace qs(*mesh, *ir);
3515 CoefficientVector coeff(qs, CoefficientStorage::FULL);
3516 if (Q) { coeff.Project(*Q); }
3517 else if (DQ) { coeff.Project(*DQ); }
3518 else { coeff.SetConstant(1.0); }
3520 if (testType == mfem::FiniteElement::CURL &&
3521 trialType == mfem::FiniteElement::CURL && dim == 3)
3525 PAHcurlL2Setup(nq, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
3529 PAHcurlHdivSetup3D(quad1D, coeffDim, ne, false, ir->GetWeights(),
3530 geom->J, coeff, pa_data);
3533 else if (testType == mfem::FiniteElement::DIV &&
3534 trialType == mfem::FiniteElement::CURL && dim == 3 &&
3535 test_fel->GetOrder() == trial_fel->GetOrder())
3537 PACurlCurlSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J, coeff,
3542 MFEM_ABORT("Unknown kernel.
");
3546 // Apply to x corresponding to DOFs in H(curl) (trial), whose curl is
3547 // integrated against H(curl) test functions corresponding to y.
3548 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
3549 static void PAHcurlL2Apply3D(const int D1D,
3553 const Array<double> &bo,
3554 const Array<double> &bc,
3555 const Array<double> &bot,
3556 const Array<double> &bct,
3557 const Array<double> &gc,
3558 const Vector &pa_data,
3562 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
3563 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
3564 // Using u = dF^{-T} \hat{u} and (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u} (p. 78 of Monk), we get
3565 // (\nabla\times u) \cdot v = 1/det(dF) \hat{\nabla}\times\hat{u}^T dF^T dF^{-T} \hat{v}
3566 // = 1/det(dF) \hat{\nabla}\times\hat{u}^T \hat{v}
3567 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
3568 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
3569 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
3571 constexpr static int VDIM = 3;
3573 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
3574 auto Bc = Reshape(bc.Read(), Q1D, D1D);
3575 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
3576 auto Bct = Reshape(bct.Read(), D1D, Q1D);
3577 auto Gc = Reshape(gc.Read(), Q1D, D1D);
3578 auto op = Reshape(pa_data.Read(), coeffDim, Q1D, Q1D, Q1D, NE);
3579 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
3580 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
3584 double curl[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
3585 // curl[qz][qy][qx] will be computed as the vector curl at each quadrature point.
3587 for (int qz = 0; qz < Q1D; ++qz)
3589 for (int qy = 0; qy < Q1D; ++qy)
3591 for (int qx = 0; qx < Q1D; ++qx)
3593 for (int c = 0; c < VDIM; ++c)
3595 curl[qz][qy][qx][c] = 0.0;
3601 // We treat x, y, z components separately for optimization specific to each.
3607 const int D1Dz = D1D;
3608 const int D1Dy = D1D;
3609 const int D1Dx = D1D - 1;
3611 for (int dz = 0; dz < D1Dz; ++dz)
3613 double gradXY[MAX_Q1D][MAX_Q1D][2];
3614 for (int qy = 0; qy < Q1D; ++qy)
3616 for (int qx = 0; qx < Q1D; ++qx)
3618 for (int d = 0; d < 2; ++d)
3620 gradXY[qy][qx][d] = 0.0;
3625 for (int dy = 0; dy < D1Dy; ++dy)
3627 double massX[MAX_Q1D];
3628 for (int qx = 0; qx < Q1D; ++qx)
3633 for (int dx = 0; dx < D1Dx; ++dx)
3635 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
3636 for (int qx = 0; qx < Q1D; ++qx)
3638 massX[qx] += t * Bo(qx,dx);
3642 for (int qy = 0; qy < Q1D; ++qy)
3644 const double wy = Bc(qy,dy);
3645 const double wDy = Gc(qy,dy);
3646 for (int qx = 0; qx < Q1D; ++qx)
3648 const double wx = massX[qx];
3649 gradXY[qy][qx][0] += wx * wDy;
3650 gradXY[qy][qx][1] += wx * wy;
3655 for (int qz = 0; qz < Q1D; ++qz)
3657 const double wz = Bc(qz,dz);
3658 const double wDz = Gc(qz,dz);
3659 for (int qy = 0; qy < Q1D; ++qy)
3661 for (int qx = 0; qx < Q1D; ++qx)
3663 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
3664 curl[qz][qy][qx][1] += gradXY[qy][qx][1] * wDz; // (u_0)_{x_2}
3665 curl[qz][qy][qx][2] -= gradXY[qy][qx][0] * wz; // -(u_0)_{x_1}
3671 osc += D1Dx * D1Dy * D1Dz;
3676 const int D1Dz = D1D;
3677 const int D1Dy = D1D - 1;
3678 const int D1Dx = D1D;
3680 for (int dz = 0; dz < D1Dz; ++dz)
3682 double gradXY[MAX_Q1D][MAX_Q1D][2];
3683 for (int qy = 0; qy < Q1D; ++qy)
3685 for (int qx = 0; qx < Q1D; ++qx)
3687 for (int d = 0; d < 2; ++d)
3689 gradXY[qy][qx][d] = 0.0;
3694 for (int dx = 0; dx < D1Dx; ++dx)
3696 double massY[MAX_Q1D];
3697 for (int qy = 0; qy < Q1D; ++qy)
3702 for (int dy = 0; dy < D1Dy; ++dy)
3704 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
3705 for (int qy = 0; qy < Q1D; ++qy)
3707 massY[qy] += t * Bo(qy,dy);
3711 for (int qx = 0; qx < Q1D; ++qx)
3713 const double wx = Bc(qx,dx);
3714 const double wDx = Gc(qx,dx);
3715 for (int qy = 0; qy < Q1D; ++qy)
3717 const double wy = massY[qy];
3718 gradXY[qy][qx][0] += wDx * wy;
3719 gradXY[qy][qx][1] += wx * wy;
3724 for (int qz = 0; qz < Q1D; ++qz)
3726 const double wz = Bc(qz,dz);
3727 const double wDz = Gc(qz,dz);
3728 for (int qy = 0; qy < Q1D; ++qy)
3730 for (int qx = 0; qx < Q1D; ++qx)
3732 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
3733 curl[qz][qy][qx][0] -= gradXY[qy][qx][1] * wDz; // -(u_1)_{x_2}
3734 curl[qz][qy][qx][2] += gradXY[qy][qx][0] * wz; // (u_1)_{x_0}
3740 osc += D1Dx * D1Dy * D1Dz;
3745 const int D1Dz = D1D - 1;
3746 const int D1Dy = D1D;
3747 const int D1Dx = D1D;
3749 for (int dx = 0; dx < D1Dx; ++dx)
3751 double gradYZ[MAX_Q1D][MAX_Q1D][2];
3752 for (int qz = 0; qz < Q1D; ++qz)
3754 for (int qy = 0; qy < Q1D; ++qy)
3756 for (int d = 0; d < 2; ++d)
3758 gradYZ[qz][qy][d] = 0.0;
3763 for (int dy = 0; dy < D1Dy; ++dy)
3765 double massZ[MAX_Q1D];
3766 for (int qz = 0; qz < Q1D; ++qz)
3771 for (int dz = 0; dz < D1Dz; ++dz)
3773 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
3774 for (int qz = 0; qz < Q1D; ++qz)
3776 massZ[qz] += t * Bo(qz,dz);
3780 for (int qy = 0; qy < Q1D; ++qy)
3782 const double wy = Bc(qy,dy);
3783 const double wDy = Gc(qy,dy);
3784 for (int qz = 0; qz < Q1D; ++qz)
3786 const double wz = massZ[qz];
3787 gradYZ[qz][qy][0] += wz * wy;
3788 gradYZ[qz][qy][1] += wz * wDy;
3793 for (int qx = 0; qx < Q1D; ++qx)
3795 const double wx = Bc(qx,dx);
3796 const double wDx = Gc(qx,dx);
3798 for (int qy = 0; qy < Q1D; ++qy)
3800 for (int qz = 0; qz < Q1D; ++qz)
3802 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
3803 curl[qz][qy][qx][0] += gradYZ[qz][qy][1] * wx; // (u_2)_{x_1}
3804 curl[qz][qy][qx][1] -= gradYZ[qz][qy][0] * wDx; // -(u_2)_{x_0}
3811 // Apply D operator.
3812 for (int qz = 0; qz < Q1D; ++qz)
3814 for (int qy = 0; qy < Q1D; ++qy)
3816 for (int qx = 0; qx < Q1D; ++qx)
3818 const double O11 = op(0,qx,qy,qz,e);
3821 for (int c = 0; c < VDIM; ++c)
3823 curl[qz][qy][qx][c] *= O11;
3828 const double O21 = op(1,qx,qy,qz,e);
3829 const double O31 = op(2,qx,qy,qz,e);
3830 const double O12 = op(3,qx,qy,qz,e);
3831 const double O22 = op(4,qx,qy,qz,e);
3832 const double O32 = op(5,qx,qy,qz,e);
3833 const double O13 = op(6,qx,qy,qz,e);
3834 const double O23 = op(7,qx,qy,qz,e);
3835 const double O33 = op(8,qx,qy,qz,e);
3836 const double curlX = curl[qz][qy][qx][0];
3837 const double curlY = curl[qz][qy][qx][1];
3838 const double curlZ = curl[qz][qy][qx][2];
3839 curl[qz][qy][qx][0] = (O11*curlX)+(O12*curlY)+(O13*curlZ);
3840 curl[qz][qy][qx][1] = (O21*curlX)+(O22*curlY)+(O23*curlZ);
3841 curl[qz][qy][qx][2] = (O31*curlX)+(O32*curlY)+(O33*curlZ);
3847 for (int qz = 0; qz < Q1D; ++qz)
3849 double massXY[MAX_D1D][MAX_D1D];
3853 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
3855 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
3856 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
3857 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
3859 for (int dy = 0; dy < D1Dy; ++dy)
3861 for (int dx = 0; dx < D1Dx; ++dx)
3866 for (int qy = 0; qy < Q1D; ++qy)
3868 double massX[MAX_D1D];
3869 for (int dx = 0; dx < D1Dx; ++dx)
3873 for (int qx = 0; qx < Q1D; ++qx)
3875 for (int dx = 0; dx < D1Dx; ++dx)
3877 massX[dx] += curl[qz][qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
3881 for (int dy = 0; dy < D1Dy; ++dy)
3883 const double wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
3884 for (int dx = 0; dx < D1Dx; ++dx)
3886 massXY[dy][dx] += massX[dx] * wy;
3891 for (int dz = 0; dz < D1Dz; ++dz)
3893 const double wz = (c == 2) ? Bot(dz,qz) : Bct(dz,qz);
3894 for (int dy = 0; dy < D1Dy; ++dy)
3896 for (int dx = 0; dx < D1Dx; ++dx)
3898 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += massXY[dy][dx] * wz;
3903 osc += D1Dx * D1Dy * D1Dz;
3906 }); // end of element loop
3909 // Apply to x corresponding to DOFs in H(curl) (trial), whose curl is
3910 // integrated against H(curl) test functions corresponding to y.
3911 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
3912 static void SmemPAHcurlL2Apply3D(const int D1D,
3916 const Array<double> &bo,
3917 const Array<double> &bc,
3918 const Array<double> &gc,
3919 const Vector &pa_data,
3923 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
3924 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
3926 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
3927 auto Bc = Reshape(bc.Read(), Q1D, D1D);
3928 auto Gc = Reshape(gc.Read(), Q1D, D1D);
3929 auto op = Reshape(pa_data.Read(), coeffDim, Q1D, Q1D, Q1D, NE);
3930 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
3931 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
3933 auto device_kernel = [=] MFEM_DEVICE (int e)
3935 constexpr int VDIM = 3;
3936 constexpr int maxCoeffDim = 9;
3938 MFEM_SHARED double sBo[MAX_D1D][MAX_Q1D];
3939 MFEM_SHARED double sBc[MAX_D1D][MAX_Q1D];
3940 MFEM_SHARED double sGc[MAX_D1D][MAX_Q1D];
3942 double opc[maxCoeffDim];
3943 MFEM_SHARED double sop[maxCoeffDim][MAX_Q1D][MAX_Q1D];
3944 MFEM_SHARED double curl[MAX_Q1D][MAX_Q1D][3];
3946 MFEM_SHARED double sX[MAX_D1D][MAX_D1D][MAX_D1D];
3948 MFEM_FOREACH_THREAD(qx,x,Q1D)
3950 MFEM_FOREACH_THREAD(qy,y,Q1D)
3952 MFEM_FOREACH_THREAD(qz,z,Q1D)
3954 for (int i=0; i<coeffDim; ++i)
3956 opc[i] = op(i,qx,qy,qz,e);
3962 const int tidx = MFEM_THREAD_ID(x);
3963 const int tidy = MFEM_THREAD_ID(y);
3964 const int tidz = MFEM_THREAD_ID(z);
3968 MFEM_FOREACH_THREAD(d,y,D1D)
3970 MFEM_FOREACH_THREAD(q,x,Q1D)
3972 sBc[d][q] = Bc(q,d);
3973 sGc[d][q] = Gc(q,d);
3976 sBo[d][q] = Bo(q,d);
3983 for (int qz=0; qz < Q1D; ++qz)
3987 MFEM_FOREACH_THREAD(qy,y,Q1D)
3989 MFEM_FOREACH_THREAD(qx,x,Q1D)
3991 for (int i=0; i<3; ++i)
3993 curl[qy][qx][i] = 0.0;
4000 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
4002 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
4003 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
4004 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
4006 MFEM_FOREACH_THREAD(dz,z,D1Dz)
4008 MFEM_FOREACH_THREAD(dy,y,D1Dy)
4010 MFEM_FOREACH_THREAD(dx,x,D1Dx)
4012 sX[dz][dy][dx] = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
4022 for (int i=0; i<coeffDim; ++i)
4024 sop[i][tidx][tidy] = opc[i];
4028 MFEM_FOREACH_THREAD(qy,y,Q1D)
4030 MFEM_FOREACH_THREAD(qx,x,Q1D)
4035 // We treat x, y, z components separately for optimization specific to each.
4036 if (c == 0) // x component
4038 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
4040 for (int dz = 0; dz < D1Dz; ++dz)
4042 const double wz = sBc[dz][qz];
4043 const double wDz = sGc[dz][qz];
4045 for (int dy = 0; dy < D1Dy; ++dy)
4047 const double wy = sBc[dy][qy];
4048 const double wDy = sGc[dy][qy];
4050 for (int dx = 0; dx < D1Dx; ++dx)
4052 const double wx = sX[dz][dy][dx] * sBo[dx][qx];
4059 curl[qy][qx][1] += v; // (u_0)_{x_2}
4060 curl[qy][qx][2] -= u; // -(u_0)_{x_1}
4062 else if (c == 1) // y component
4064 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
4066 for (int dz = 0; dz < D1Dz; ++dz)
4068 const double wz = sBc[dz][qz];
4069 const double wDz = sGc[dz][qz];
4071 for (int dy = 0; dy < D1Dy; ++dy)
4073 const double wy = sBo[dy][qy];
4075 for (int dx = 0; dx < D1Dx; ++dx)
4077 const double t = sX[dz][dy][dx];
4078 const double wx = t * sBc[dx][qx];
4079 const double wDx = t * sGc[dx][qx];
4087 curl[qy][qx][0] -= v; // -(u_1)_{x_2}
4088 curl[qy][qx][2] += u; // (u_1)_{x_0}
4092 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
4094 for (int dz = 0; dz < D1Dz; ++dz)
4096 const double wz = sBo[dz][qz];
4098 for (int dy = 0; dy < D1Dy; ++dy)
4100 const double wy = sBc[dy][qy];
4101 const double wDy = sGc[dy][qy];
4103 for (int dx = 0; dx < D1Dx; ++dx)
4105 const double t = sX[dz][dy][dx];
4106 const double wx = t * sBc[dx][qx];
4107 const double wDx = t * sGc[dx][qx];
4115 curl[qy][qx][0] += v; // (u_2)_{x_1}
4116 curl[qy][qx][1] -= u; // -(u_2)_{x_0}
4122 osc += D1Dx * D1Dy * D1Dz;
4130 MFEM_FOREACH_THREAD(dz,z,D1D)
4132 const double wcz = sBc[dz][qz];
4133 const double wz = (dz < D1D-1) ? sBo[dz][qz] : 0.0;
4135 MFEM_FOREACH_THREAD(dy,y,D1D)
4137 MFEM_FOREACH_THREAD(dx,x,D1D)
4139 for (int qy = 0; qy < Q1D; ++qy)
4141 const double wcy = sBc[dy][qy];
4142 const double wy = (dy < D1D-1) ? sBo[dy][qy] : 0.0;
4144 for (int qx = 0; qx < Q1D; ++qx)
4146 const double O11 = sop[0][qx][qy];
4150 c1 = O11 * curl[qy][qx][0];
4151 c2 = O11 * curl[qy][qx][1];
4152 c3 = O11 * curl[qy][qx][2];
4156 const double O21 = sop[1][qx][qy];
4157 const double O31 = sop[2][qx][qy];
4158 const double O12 = sop[3][qx][qy];
4159 const double O22 = sop[4][qx][qy];
4160 const double O32 = sop[5][qx][qy];
4161 const double O13 = sop[6][qx][qy];
4162 const double O23 = sop[7][qx][qy];
4163 const double O33 = sop[8][qx][qy];
4164 c1 = (O11*curl[qy][qx][0])+(O12*curl[qy][qx][1])+(O13*curl[qy][qx][2]);
4165 c2 = (O21*curl[qy][qx][0])+(O22*curl[qy][qx][1])+(O23*curl[qy][qx][2]);
4166 c3 = (O31*curl[qy][qx][0])+(O32*curl[qy][qx][1])+(O33*curl[qy][qx][2]);
4169 const double wcx = sBc[dx][qx];
4173 const double wx = sBo[dx][qx];
4174 dxyz1 += c1 * wx * wcy * wcz;
4177 dxyz2 += c2 * wcx * wy * wcz;
4178 dxyz3 += c3 * wcx * wcy * wz;
4187 MFEM_FOREACH_THREAD(dz,z,D1D)
4189 MFEM_FOREACH_THREAD(dy,y,D1D)
4191 MFEM_FOREACH_THREAD(dx,x,D1D)
4195 Y(dx + ((dy + (dz * D1D)) * (D1D-1)), e) += dxyz1;
4199 Y(dx + ((dy + (dz * (D1D-1))) * D1D) + ((D1D-1)*D1D*D1D), e) += dxyz2;
4203 Y(dx + ((dy + (dz * D1D)) * D1D) + (2*(D1D-1)*D1D*D1D), e) += dxyz3;
4209 }; // end of element loop
4211 auto host_kernel = [&] MFEM_LAMBDA (int)
4213 MFEM_ABORT_KERNEL("This kernel should only be used on GPU.
");
4216 ForallWrap<3>(true, NE, device_kernel, host_kernel, Q1D, Q1D, Q1D);
4219 // Apply to x corresponding to DOFs in H(curl) (trial), whose curl is
4220 // integrated against H(div) test functions corresponding to y.
4221 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
4222 static void PAHcurlHdivApply3D(const int D1D,
4226 const Array<double> &bo,
4227 const Array<double> &bc,
4228 const Array<double> &bot,
4229 const Array<double> &bct,
4230 const Array<double> &gc,
4231 const Vector &pa_data,
4235 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
4236 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
4237 // Using Piola transformations (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u}
4238 // for u in H(curl) and w = (1 / det (dF)) dF \hat{w} for w in H(div), we get
4239 // (\nabla\times u) \cdot w = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{w}
4240 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
4241 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
4242 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
4244 constexpr static int VDIM = 3;
4246 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
4247 auto Bc = Reshape(bc.Read(), Q1D, D1D);
4248 auto Bot = Reshape(bot.Read(), D1Dtest-1, Q1D);
4249 auto Bct = Reshape(bct.Read(), D1Dtest, Q1D);
4250 auto Gc = Reshape(gc.Read(), Q1D, D1D);
4251 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
4252 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
4253 auto Y = Reshape(y.ReadWrite(), 3*(D1Dtest-1)*(D1Dtest-1)*D1D, NE);
4257 double curl[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
4258 // curl[qz][qy][qx] will be computed as the vector curl at each quadrature point.
4260 for (int qz = 0; qz < Q1D; ++qz)
4262 for (int qy = 0; qy < Q1D; ++qy)
4264 for (int qx = 0; qx < Q1D; ++qx)
4266 for (int c = 0; c < VDIM; ++c)
4268 curl[qz][qy][qx][c] = 0.0;
4274 // We treat x, y, z components separately for optimization specific to each.
4280 const int D1Dz = D1D;
4281 const int D1Dy = D1D;
4282 const int D1Dx = D1D - 1;
4284 for (int dz = 0; dz < D1Dz; ++dz)
4286 double gradXY[MAX_Q1D][MAX_Q1D][2];
4287 for (int qy = 0; qy < Q1D; ++qy)
4289 for (int qx = 0; qx < Q1D; ++qx)
4291 for (int d = 0; d < 2; ++d)
4293 gradXY[qy][qx][d] = 0.0;
4298 for (int dy = 0; dy < D1Dy; ++dy)
4300 double massX[MAX_Q1D];
4301 for (int qx = 0; qx < Q1D; ++qx)
4306 for (int dx = 0; dx < D1Dx; ++dx)
4308 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
4309 for (int qx = 0; qx < Q1D; ++qx)
4311 massX[qx] += t * Bo(qx,dx);
4315 for (int qy = 0; qy < Q1D; ++qy)
4317 const double wy = Bc(qy,dy);
4318 const double wDy = Gc(qy,dy);
4319 for (int qx = 0; qx < Q1D; ++qx)
4321 const double wx = massX[qx];
4322 gradXY[qy][qx][0] += wx * wDy;
4323 gradXY[qy][qx][1] += wx * wy;
4328 for (int qz = 0; qz < Q1D; ++qz)
4330 const double wz = Bc(qz,dz);
4331 const double wDz = Gc(qz,dz);
4332 for (int qy = 0; qy < Q1D; ++qy)
4334 for (int qx = 0; qx < Q1D; ++qx)
4336 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
4337 curl[qz][qy][qx][1] += gradXY[qy][qx][1] * wDz; // (u_0)_{x_2}
4338 curl[qz][qy][qx][2] -= gradXY[qy][qx][0] * wz; // -(u_0)_{x_1}
4344 osc += D1Dx * D1Dy * D1Dz;
4349 const int D1Dz = D1D;
4350 const int D1Dy = D1D - 1;
4351 const int D1Dx = D1D;
4353 for (int dz = 0; dz < D1Dz; ++dz)
4355 double gradXY[MAX_Q1D][MAX_Q1D][2];
4356 for (int qy = 0; qy < Q1D; ++qy)
4358 for (int qx = 0; qx < Q1D; ++qx)
4360 for (int d = 0; d < 2; ++d)
4362 gradXY[qy][qx][d] = 0.0;
4367 for (int dx = 0; dx < D1Dx; ++dx)
4369 double massY[MAX_Q1D];
4370 for (int qy = 0; qy < Q1D; ++qy)
4375 for (int dy = 0; dy < D1Dy; ++dy)
4377 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
4378 for (int qy = 0; qy < Q1D; ++qy)
4380 massY[qy] += t * Bo(qy,dy);
4384 for (int qx = 0; qx < Q1D; ++qx)
4386 const double wx = Bc(qx,dx);
4387 const double wDx = Gc(qx,dx);
4388 for (int qy = 0; qy < Q1D; ++qy)
4390 const double wy = massY[qy];
4391 gradXY[qy][qx][0] += wDx * wy;
4392 gradXY[qy][qx][1] += wx * wy;
4397 for (int qz = 0; qz < Q1D; ++qz)
4399 const double wz = Bc(qz,dz);
4400 const double wDz = Gc(qz,dz);
4401 for (int qy = 0; qy < Q1D; ++qy)
4403 for (int qx = 0; qx < Q1D; ++qx)
4405 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
4406 curl[qz][qy][qx][0] -= gradXY[qy][qx][1] * wDz; // -(u_1)_{x_2}
4407 curl[qz][qy][qx][2] += gradXY[qy][qx][0] * wz; // (u_1)_{x_0}
4413 osc += D1Dx * D1Dy * D1Dz;
4418 const int D1Dz = D1D - 1;
4419 const int D1Dy = D1D;
4420 const int D1Dx = D1D;
4422 for (int dx = 0; dx < D1Dx; ++dx)
4424 double gradYZ[MAX_Q1D][MAX_Q1D][2];
4425 for (int qz = 0; qz < Q1D; ++qz)
4427 for (int qy = 0; qy < Q1D; ++qy)
4429 for (int d = 0; d < 2; ++d)
4431 gradYZ[qz][qy][d] = 0.0;
4436 for (int dy = 0; dy < D1Dy; ++dy)
4438 double massZ[MAX_Q1D];
4439 for (int qz = 0; qz < Q1D; ++qz)
4444 for (int dz = 0; dz < D1Dz; ++dz)
4446 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
4447 for (int qz = 0; qz < Q1D; ++qz)
4449 massZ[qz] += t * Bo(qz,dz);
4453 for (int qy = 0; qy < Q1D; ++qy)
4455 const double wy = Bc(qy,dy);
4456 const double wDy = Gc(qy,dy);
4457 for (int qz = 0; qz < Q1D; ++qz)
4459 const double wz = massZ[qz];
4460 gradYZ[qz][qy][0] += wz * wy;
4461 gradYZ[qz][qy][1] += wz * wDy;
4466 for (int qx = 0; qx < Q1D; ++qx)
4468 const double wx = Bc(qx,dx);
4469 const double wDx = Gc(qx,dx);
4471 for (int qy = 0; qy < Q1D; ++qy)
4473 for (int qz = 0; qz < Q1D; ++qz)
4475 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
4476 curl[qz][qy][qx][0] += gradYZ[qz][qy][1] * wx; // (u_2)_{x_1}
4477 curl[qz][qy][qx][1] -= gradYZ[qz][qy][0] * wDx; // -(u_2)_{x_0}
4484 // Apply D operator.
4485 for (int qz = 0; qz < Q1D; ++qz)
4487 for (int qy = 0; qy < Q1D; ++qy)
4489 for (int qx = 0; qx < Q1D; ++qx)
4491 const double O11 = op(qx,qy,qz,0,e);
4492 const double O12 = op(qx,qy,qz,1,e);
4493 const double O13 = op(qx,qy,qz,2,e);
4494 const double O22 = op(qx,qy,qz,3,e);
4495 const double O23 = op(qx,qy,qz,4,e);
4496 const double O33 = op(qx,qy,qz,5,e);
4498 const double c1 = (O11 * curl[qz][qy][qx][0]) + (O12 * curl[qz][qy][qx][1]) +
4499 (O13 * curl[qz][qy][qx][2]);
4500 const double c2 = (O12 * curl[qz][qy][qx][0]) + (O22 * curl[qz][qy][qx][1]) +
4501 (O23 * curl[qz][qy][qx][2]);
4502 const double c3 = (O13 * curl[qz][qy][qx][0]) + (O23 * curl[qz][qy][qx][1]) +
4503 (O33 * curl[qz][qy][qx][2]);
4505 curl[qz][qy][qx][0] = c1;
4506 curl[qz][qy][qx][1] = c2;
4507 curl[qz][qy][qx][2] = c3;
4512 for (int qz = 0; qz < Q1D; ++qz)
4514 double massXY[HCURL_MAX_D1D][HCURL_MAX_D1D]; // Assuming HDIV_MAX_D1D <= HCURL_MAX_D1D
4518 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
4520 const int D1Dz = (c == 2) ? D1Dtest : D1Dtest - 1;
4521 const int D1Dy = (c == 1) ? D1Dtest : D1Dtest - 1;
4522 const int D1Dx = (c == 0) ? D1Dtest : D1Dtest - 1;
4524 for (int dy = 0; dy < D1Dy; ++dy)
4526 for (int dx = 0; dx < D1Dx; ++dx)
4531 for (int qy = 0; qy < Q1D; ++qy)
4533 double massX[HCURL_MAX_D1D];
4534 for (int dx = 0; dx < D1Dx; ++dx)
4538 for (int qx = 0; qx < Q1D; ++qx)
4540 for (int dx = 0; dx < D1Dx; ++dx)
4542 massX[dx] += curl[qz][qy][qx][c] *
4543 ((c == 0) ? Bct(dx,qx) : Bot(dx,qx));
4546 for (int dy = 0; dy < D1Dy; ++dy)
4548 const double wy = (c == 1) ? Bct(dy,qy) : Bot(dy,qy);
4549 for (int dx = 0; dx < D1Dx; ++dx)
4551 massXY[dy][dx] += massX[dx] * wy;
4556 for (int dz = 0; dz < D1Dz; ++dz)
4558 const double wz = (c == 2) ? Bct(dz,qz) : Bot(dz,qz);
4559 for (int dy = 0; dy < D1Dy; ++dy)
4561 for (int dx = 0; dx < D1Dx; ++dx)
4563 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) +=
4564 massXY[dy][dx] * wz;
4569 osc += D1Dx * D1Dy * D1Dz;
4572 }); // end of element loop
4575 // Apply to x corresponding to DOFs in H(div) (test), integrated against the
4576 // curl of H(curl) trial functions corresponding to y.
4577 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
4578 static void PAHcurlHdivApply3DTranspose(const int D1D,
4582 const Array<double> &bo,
4583 const Array<double> &bc,
4584 const Array<double> &bot,
4585 const Array<double> &bct,
4586 const Array<double> &gct,
4587 const Vector &pa_data,
4591 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
4592 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
4593 // Using Piola transformations (\nabla\times u) F = 1/det(dF) dF \hat{\nabla}\times\hat{u}
4594 // for u in H(curl) and w = (1 / det (dF)) dF \hat{w} for w in H(div), we get
4595 // (\nabla\times u) \cdot w = 1/det(dF)^2 \hat{\nabla}\times\hat{u}^T dF^T dF \hat{w}
4596 // If c = 0, \hat{\nabla}\times\hat{u} reduces to [0, (u_0)_{x_2}, -(u_0)_{x_1}]
4597 // If c = 1, \hat{\nabla}\times\hat{u} reduces to [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
4598 // If c = 2, \hat{\nabla}\times\hat{u} reduces to [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
4600 constexpr static int VDIM = 3;
4602 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
4603 auto Bc = Reshape(bc.Read(), Q1D, D1D);
4604 auto Bot = Reshape(bot.Read(), D1Dtest-1, Q1D);
4605 auto Bct = Reshape(bct.Read(), D1Dtest, Q1D);
4606 auto Gct = Reshape(gct.Read(), D1D, Q1D);
4607 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
4608 auto X = Reshape(x.Read(), 3*(D1Dtest-1)*(D1Dtest-1)*D1D, NE);
4609 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
4613 double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM]; // Assuming HDIV_MAX_D1D <= HCURL_MAX_D1D
4615 for (int qz = 0; qz < Q1D; ++qz)
4617 for (int qy = 0; qy < Q1D; ++qy)
4619 for (int qx = 0; qx < Q1D; ++qx)
4621 for (int c = 0; c < VDIM; ++c)
4623 mass[qz][qy][qx][c] = 0.0;
4631 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
4633 const int D1Dz = (c == 2) ? D1D : D1D - 1;
4634 const int D1Dy = (c == 1) ? D1D : D1D - 1;
4635 const int D1Dx = (c == 0) ? D1D : D1D - 1;
4637 for (int dz = 0; dz < D1Dz; ++dz)
4639 double massXY[HDIV_MAX_Q1D][HDIV_MAX_Q1D];
4640 for (int qy = 0; qy < Q1D; ++qy)
4642 for (int qx = 0; qx < Q1D; ++qx)
4644 massXY[qy][qx] = 0.0;
4648 for (int dy = 0; dy < D1Dy; ++dy)
4650 double massX[HDIV_MAX_Q1D];
4651 for (int qx = 0; qx < Q1D; ++qx)
4656 for (int dx = 0; dx < D1Dx; ++dx)
4658 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
4659 for (int qx = 0; qx < Q1D; ++qx)
4661 massX[qx] += t * ((c == 0) ? Bc(qx,dx) : Bo(qx,dx));
4665 for (int qy = 0; qy < Q1D; ++qy)
4667 const double wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
4668 for (int qx = 0; qx < Q1D; ++qx)
4670 const double wx = massX[qx];
4671 massXY[qy][qx] += wx * wy;
4676 for (int qz = 0; qz < Q1D; ++qz)
4678 const double wz = (c == 2) ? Bc(qz,dz) : Bo(qz,dz);
4679 for (int qy = 0; qy < Q1D; ++qy)
4681 for (int qx = 0; qx < Q1D; ++qx)
4683 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
4689 osc += D1Dx * D1Dy * D1Dz;
4690 } // loop (c) over components
4692 // Apply D operator.
4693 for (int qz = 0; qz < Q1D; ++qz)
4695 for (int qy = 0; qy < Q1D; ++qy)
4697 for (int qx = 0; qx < Q1D; ++qx)
4699 const double O11 = op(qx,qy,qz,0,e);
4700 const double O12 = op(qx,qy,qz,1,e);
4701 const double O13 = op(qx,qy,qz,2,e);
4702 const double O22 = op(qx,qy,qz,3,e);
4703 const double O23 = op(qx,qy,qz,4,e);
4704 const double O33 = op(qx,qy,qz,5,e);
4705 const double massX = mass[qz][qy][qx][0];
4706 const double massY = mass[qz][qy][qx][1];
4707 const double massZ = mass[qz][qy][qx][2];
4708 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
4709 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
4710 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
4718 const int D1Dz = D1D;
4719 const int D1Dy = D1D;
4720 const int D1Dx = D1D - 1;
4722 for (int qz = 0; qz < Q1D; ++qz)
4724 double gradXY12[MAX_D1D][MAX_D1D];
4725 double gradXY21[MAX_D1D][MAX_D1D];
4727 for (int dy = 0; dy < D1Dy; ++dy)
4729 for (int dx = 0; dx < D1Dx; ++dx)
4731 gradXY12[dy][dx] = 0.0;
4732 gradXY21[dy][dx] = 0.0;
4735 for (int qy = 0; qy < Q1D; ++qy)
4737 double massX[MAX_D1D][2];
4738 for (int dx = 0; dx < D1Dx; ++dx)
4740 for (int n = 0; n < 2; ++n)
4745 for (int qx = 0; qx < Q1D; ++qx)
4747 for (int dx = 0; dx < D1Dx; ++dx)
4749 const double wx = Bot(dx,qx);
4751 massX[dx][0] += wx * mass[qz][qy][qx][1];
4752 massX[dx][1] += wx * mass[qz][qy][qx][2];
4755 for (int dy = 0; dy < D1Dy; ++dy)
4757 const double wy = Bct(dy,qy);
4758 const double wDy = Gct(dy,qy);
4760 for (int dx = 0; dx < D1Dx; ++dx)
4762 gradXY21[dy][dx] += massX[dx][0] * wy;
4763 gradXY12[dy][dx] += massX[dx][1] * wDy;
4768 for (int dz = 0; dz < D1Dz; ++dz)
4770 const double wz = Bct(dz,qz);
4771 const double wDz = Gct(dz,qz);
4772 for (int dy = 0; dy < D1Dy; ++dy)
4774 for (int dx = 0; dx < D1Dx; ++dx)
4776 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
4777 // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
4778 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
4779 e) += (gradXY21[dy][dx] * wDz) - (gradXY12[dy][dx] * wz);
4785 osc += D1Dx * D1Dy * D1Dz;
4790 const int D1Dz = D1D;
4791 const int D1Dy = D1D - 1;
4792 const int D1Dx = D1D;
4794 for (int qz = 0; qz < Q1D; ++qz)
4796 double gradXY02[MAX_D1D][MAX_D1D];
4797 double gradXY20[MAX_D1D][MAX_D1D];
4799 for (int dy = 0; dy < D1Dy; ++dy)
4801 for (int dx = 0; dx < D1Dx; ++dx)
4803 gradXY02[dy][dx] = 0.0;
4804 gradXY20[dy][dx] = 0.0;
4807 for (int qx = 0; qx < Q1D; ++qx)
4809 double massY[MAX_D1D][2];
4810 for (int dy = 0; dy < D1Dy; ++dy)
4815 for (int qy = 0; qy < Q1D; ++qy)
4817 for (int dy = 0; dy < D1Dy; ++dy)
4819 const double wy = Bot(dy,qy);
4821 massY[dy][0] += wy * mass[qz][qy][qx][2];
4822 massY[dy][1] += wy * mass[qz][qy][qx][0];
4825 for (int dx = 0; dx < D1Dx; ++dx)
4827 const double wx = Bct(dx,qx);
4828 const double wDx = Gct(dx,qx);
4830 for (int dy = 0; dy < D1Dy; ++dy)
4832 gradXY02[dy][dx] += massY[dy][0] * wDx;
4833 gradXY20[dy][dx] += massY[dy][1] * wx;
4838 for (int dz = 0; dz < D1Dz; ++dz)
4840 const double wz = Bct(dz,qz);
4841 const double wDz = Gct(dz,qz);
4842 for (int dy = 0; dy < D1Dy; ++dy)
4844 for (int dx = 0; dx < D1Dx; ++dx)
4846 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
4847 // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
4848 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
4849 e) += (-gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
4855 osc += D1Dx * D1Dy * D1Dz;
4860 const int D1Dz = D1D - 1;
4861 const int D1Dy = D1D;
4862 const int D1Dx = D1D;
4864 for (int qx = 0; qx < Q1D; ++qx)
4866 double gradYZ01[MAX_D1D][MAX_D1D];
4867 double gradYZ10[MAX_D1D][MAX_D1D];
4869 for (int dy = 0; dy < D1Dy; ++dy)
4871 for (int dz = 0; dz < D1Dz; ++dz)
4873 gradYZ01[dz][dy] = 0.0;
4874 gradYZ10[dz][dy] = 0.0;
4877 for (int qy = 0; qy < Q1D; ++qy)
4879 double massZ[MAX_D1D][2];
4880 for (int dz = 0; dz < D1Dz; ++dz)
4882 for (int n = 0; n < 2; ++n)
4887 for (int qz = 0; qz < Q1D; ++qz)
4889 for (int dz = 0; dz < D1Dz; ++dz)
4891 const double wz = Bot(dz,qz);
4893 massZ[dz][0] += wz * mass[qz][qy][qx][0];
4894 massZ[dz][1] += wz * mass[qz][qy][qx][1];
4897 for (int dy = 0; dy < D1Dy; ++dy)
4899 const double wy = Bct(dy,qy);
4900 const double wDy = Gct(dy,qy);
4902 for (int dz = 0; dz < D1Dz; ++dz)
4904 gradYZ01[dz][dy] += wy * massZ[dz][1];
4905 gradYZ10[dz][dy] += wDy * massZ[dz][0];
4910 for (int dx = 0; dx < D1Dx; ++dx)
4912 const double wx = Bct(dx,qx);
4913 const double wDx = Gct(dx,qx);
4915 for (int dy = 0; dy < D1Dy; ++dy)
4917 for (int dz = 0; dz < D1Dz; ++dz)
4919 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
4920 // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
4921 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
4922 e) += (gradYZ10[dz][dy] * wx) - (gradYZ01[dz][dy] * wDx);
4928 }); // end of element loop
4931 void MixedVectorCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
4933 if (testType == mfem::FiniteElement::CURL &&
4934 trialType == mfem::FiniteElement::CURL && dim == 3)
4936 const int ndata = coeffDim == 1 ? 1 : 9;
4938 if (Device::Allows(Backend::DEVICE_MASK))
4940 const int ID = (dofs1D << 4) | quad1D;
4943 case 0x23: return SmemPAHcurlL2Apply3D<2,3>(dofs1D, quad1D, ndata, ne,
4945 mapsC->G, pa_data, x, y);
4946 case 0x34: return SmemPAHcurlL2Apply3D<3,4>(dofs1D, quad1D, ndata, ne,
4948 mapsC->G, pa_data, x, y);
4949 case 0x45: return SmemPAHcurlL2Apply3D<4,5>(dofs1D, quad1D, ndata, ne,
4951 mapsC->G, pa_data, x, y);
4952 case 0x56: return SmemPAHcurlL2Apply3D<5,6>(dofs1D, quad1D, ndata, ne,
4954 mapsC->G, pa_data, x, y);
4955 default: return SmemPAHcurlL2Apply3D(dofs1D, quad1D, ndata, ne,
4956 mapsO->B, mapsC->B, mapsC->G,
4961 PAHcurlL2Apply3D(dofs1D, quad1D, ndata, ne, mapsO->B, mapsC->B,
4962 mapsO->Bt, mapsC->Bt, mapsC->G, pa_data, x, y);
4964 else if (testType == mfem::FiniteElement::DIV &&
4965 trialType == mfem::FiniteElement::CURL && dim == 3)
4966 PAHcurlHdivApply3D(dofs1D, dofs1Dtest, quad1D, ne, mapsO->B,
4967 mapsC->B, mapsOtest->Bt, mapsCtest->Bt, mapsC->G,
4975 void MixedVectorCurlIntegrator::AddMultTransposePA(const Vector &x,
4978 if (testType == mfem::FiniteElement::DIV &&
4979 trialType == mfem::FiniteElement::CURL && dim == 3)
4980 PAHcurlHdivApply3DTranspose(dofs1D, dofs1Dtest, quad1D, ne, mapsO->B,
4981 mapsC->B, mapsOtest->Bt, mapsCtest->Bt,
4982 mapsC->Gt, pa_data, x, y);
4989 void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
4991 const FiniteElementSpace &test_fes)
4993 // Assumes tensor-product elements, with vector test and trial spaces.
4994 Mesh *mesh = trial_fes.GetMesh();
4995 const FiniteElement *trial_fel = trial_fes.GetFE(0);
4996 const FiniteElement *test_fel = test_fes.GetFE(0);
4998 const VectorTensorFiniteElement *trial_el =
4999 dynamic_cast<const VectorTensorFiniteElement*>(trial_fel);
5002 const VectorTensorFiniteElement *test_el =
5003 dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
5006 const IntegrationRule *ir
5007 = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
5008 *mesh->GetElementTransformation(0));
5009 const int dims = trial_el->GetDim();
5010 MFEM_VERIFY(dims == 3, "");
5012 const int nq = ir->GetNPoints();
5013 dim = mesh->Dimension();
5014 MFEM_VERIFY(dim == 3, "");
5016 MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(), "");
5018 ne = trial_fes.GetNE();
5019 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
5020 mapsC = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
5021 mapsO = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
5022 dofs1D = mapsC->ndof;
5023 quad1D = mapsC->nqpt;
5025 MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
5027 testType = test_el->GetDerivType();
5028 trialType = trial_el->GetDerivType();
5030 const bool curlSpaces = (testType == mfem::FiniteElement::CURL &&
5031 trialType == mfem::FiniteElement::CURL);
5033 const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
5035 coeffDim = DQ ? 3 : 1;
5036 const int ndata = curlSpaces ? (DQ ? 9 : 1) : symmDims;
5038 pa_data.SetSize(ndata * nq * ne, Device::GetMemoryType());
5040 QuadratureSpace qs(*mesh, *ir);
5041 CoefficientVector coeff(qs, CoefficientStorage::FULL);
5042 if (Q) { coeff.Project(*Q); }
5043 else if (DQ) { coeff.Project(*DQ); }
5044 else { coeff.SetConstant(1.0); }
5046 if (trialType == mfem::FiniteElement::CURL && dim == 3)
5050 PAHcurlL2Setup(nq, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
5054 PAHcurlHdivSetup3D(quad1D, coeffDim, ne, false, ir->GetWeights(),
5055 geom->J, coeff, pa_data);
5058 else if (trialType == mfem::FiniteElement::DIV && dim == 3 &&
5059 test_el->GetOrder() == trial_el->GetOrder())
5061 PACurlCurlSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J, coeff,
5066 MFEM_ABORT("Unknown kernel.
");
5070 // Apply to x corresponding to DOFs in H(curl) (trial), integrated against curl
5071 // of H(curl) test functions corresponding to y.
5072 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
5073 static void PAHcurlL2Apply3DTranspose(const int D1D,
5077 const Array<double> &bo,
5078 const Array<double> &bc,
5079 const Array<double> &bot,
5080 const Array<double> &bct,
5081 const Array<double> &gct,
5082 const Vector &pa_data,
5086 // See PAHcurlL2Apply3D for comments.
5088 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
5089 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
5091 constexpr static int VDIM = 3;
5093 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
5094 auto Bc = Reshape(bc.Read(), Q1D, D1D);
5095 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
5096 auto Bct = Reshape(bct.Read(), D1D, Q1D);
5097 auto Gct = Reshape(gct.Read(), D1D, Q1D);
5098 auto op = Reshape(pa_data.Read(), coeffDim, Q1D, Q1D, Q1D, NE);
5099 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
5100 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
5104 double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
5106 for (int qz = 0; qz < Q1D; ++qz)
5108 for (int qy = 0; qy < Q1D; ++qy)
5110 for (int qx = 0; qx < Q1D; ++qx)
5112 for (int c = 0; c < VDIM; ++c)
5114 mass[qz][qy][qx][c] = 0.0;
5122 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
5124 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
5125 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
5126 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
5128 for (int dz = 0; dz < D1Dz; ++dz)
5130 double massXY[MAX_Q1D][MAX_Q1D];
5131 for (int qy = 0; qy < Q1D; ++qy)
5133 for (int qx = 0; qx < Q1D; ++qx)
5135 massXY[qy][qx] = 0.0;
5139 for (int dy = 0; dy < D1Dy; ++dy)
5141 double massX[MAX_Q1D];
5142 for (int qx = 0; qx < Q1D; ++qx)
5147 for (int dx = 0; dx < D1Dx; ++dx)
5149 const double t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
5150 for (int qx = 0; qx < Q1D; ++qx)
5152 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
5156 for (int qy = 0; qy < Q1D; ++qy)
5158 const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
5159 for (int qx = 0; qx < Q1D; ++qx)
5161 const double wx = massX[qx];
5162 massXY[qy][qx] += wx * wy;
5167 for (int qz = 0; qz < Q1D; ++qz)
5169 const double wz = (c == 2) ? Bo(qz,dz) : Bc(qz,dz);
5170 for (int qy = 0; qy < Q1D; ++qy)
5172 for (int qx = 0; qx < Q1D; ++qx)
5174 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
5180 osc += D1Dx * D1Dy * D1Dz;
5181 } // loop (c) over components
5183 // Apply D operator.
5184 for (int qz = 0; qz < Q1D; ++qz)
5186 for (int qy = 0; qy < Q1D; ++qy)
5188 for (int qx = 0; qx < Q1D; ++qx)
5190 const double O11 = op(0,qx,qy,qz,e);
5193 for (int c = 0; c < VDIM; ++c)
5195 mass[qz][qy][qx][c] *= O11;
5200 const double O12 = op(1,qx,qy,qz,e);
5201 const double O13 = op(2,qx,qy,qz,e);
5202 const double O21 = op(3,qx,qy,qz,e);
5203 const double O22 = op(4,qx,qy,qz,e);
5204 const double O23 = op(5,qx,qy,qz,e);
5205 const double O31 = op(6,qx,qy,qz,e);
5206 const double O32 = op(7,qx,qy,qz,e);
5207 const double O33 = op(8,qx,qy,qz,e);
5208 const double massX = mass[qz][qy][qx][0];
5209 const double massY = mass[qz][qy][qx][1];
5210 const double massZ = mass[qz][qy][qx][2];
5211 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
5212 mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
5213 mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
5222 const int D1Dz = D1D;
5223 const int D1Dy = D1D;
5224 const int D1Dx = D1D - 1;
5226 for (int qz = 0; qz < Q1D; ++qz)
5228 double gradXY12[MAX_D1D][MAX_D1D];
5229 double gradXY21[MAX_D1D][MAX_D1D];
5231 for (int dy = 0; dy < D1Dy; ++dy)
5233 for (int dx = 0; dx < D1Dx; ++dx)
5235 gradXY12[dy][dx] = 0.0;
5236 gradXY21[dy][dx] = 0.0;
5239 for (int qy = 0; qy < Q1D; ++qy)
5241 double massX[MAX_D1D][2];
5242 for (int dx = 0; dx < D1Dx; ++dx)
5244 for (int n = 0; n < 2; ++n)
5249 for (int qx = 0; qx < Q1D; ++qx)
5251 for (int dx = 0; dx < D1Dx; ++dx)
5253 const double wx = Bot(dx,qx);
5255 massX[dx][0] += wx * mass[qz][qy][qx][1];
5256 massX[dx][1] += wx * mass[qz][qy][qx][2];
5259 for (int dy = 0; dy < D1Dy; ++dy)
5261 const double wy = Bct(dy,qy);
5262 const double wDy = Gct(dy,qy);
5264 for (int dx = 0; dx < D1Dx; ++dx)
5266 gradXY21[dy][dx] += massX[dx][0] * wy;
5267 gradXY12[dy][dx] += massX[dx][1] * wDy;
5272 for (int dz = 0; dz < D1Dz; ++dz)
5274 const double wz = Bct(dz,qz);
5275 const double wDz = Gct(dz,qz);
5276 for (int dy = 0; dy < D1Dy; ++dy)
5278 for (int dx = 0; dx < D1Dx; ++dx)
5280 // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
5281 // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
5282 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
5283 e) += (gradXY21[dy][dx] * wDz) - (gradXY12[dy][dx] * wz);
5289 osc += D1Dx * D1Dy * D1Dz;
5294 const int D1Dz = D1D;
5295 const int D1Dy = D1D - 1;
5296 const int D1Dx = D1D;
5298 for (int qz = 0; qz < Q1D; ++qz)
5300 double gradXY02[MAX_D1D][MAX_D1D];
5301 double gradXY20[MAX_D1D][MAX_D1D];
5303 for (int dy = 0; dy < D1Dy; ++dy)
5305 for (int dx = 0; dx < D1Dx; ++dx)
5307 gradXY02[dy][dx] = 0.0;
5308 gradXY20[dy][dx] = 0.0;
5311 for (int qx = 0; qx < Q1D; ++qx)
5313 double massY[MAX_D1D][2];
5314 for (int dy = 0; dy < D1Dy; ++dy)
5319 for (int qy = 0; qy < Q1D; ++qy)
5321 for (int dy = 0; dy < D1Dy; ++dy)
5323 const double wy = Bot(dy,qy);
5325 massY[dy][0] += wy * mass[qz][qy][qx][2];
5326 massY[dy][1] += wy * mass[qz][qy][qx][0];
5329 for (int dx = 0; dx < D1Dx; ++dx)
5331 const double wx = Bct(dx,qx);
5332 const double wDx = Gct(dx,qx);
5334 for (int dy = 0; dy < D1Dy; ++dy)
5336 gradXY02[dy][dx] += massY[dy][0] * wDx;
5337 gradXY20[dy][dx] += massY[dy][1] * wx;
5342 for (int dz = 0; dz < D1Dz; ++dz)
5344 const double wz = Bct(dz,qz);
5345 const double wDz = Gct(dz,qz);
5346 for (int dy = 0; dy < D1Dy; ++dy)
5348 for (int dx = 0; dx < D1Dx; ++dx)
5350 // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
5351 // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
5352 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
5353 e) += (-gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
5359 osc += D1Dx * D1Dy * D1Dz;
5364 const int D1Dz = D1D - 1;
5365 const int D1Dy = D1D;
5366 const int D1Dx = D1D;
5368 for (int qx = 0; qx < Q1D; ++qx)
5370 double gradYZ01[MAX_D1D][MAX_D1D];
5371 double gradYZ10[MAX_D1D][MAX_D1D];
5373 for (int dy = 0; dy < D1Dy; ++dy)
5375 for (int dz = 0; dz < D1Dz; ++dz)
5377 gradYZ01[dz][dy] = 0.0;
5378 gradYZ10[dz][dy] = 0.0;
5381 for (int qy = 0; qy < Q1D; ++qy)
5383 double massZ[MAX_D1D][2];
5384 for (int dz = 0; dz < D1Dz; ++dz)
5386 for (int n = 0; n < 2; ++n)
5391 for (int qz = 0; qz < Q1D; ++qz)
5393 for (int dz = 0; dz < D1Dz; ++dz)
5395 const double wz = Bot(dz,qz);
5397 massZ[dz][0] += wz * mass[qz][qy][qx][0];
5398 massZ[dz][1] += wz * mass[qz][qy][qx][1];
5401 for (int dy = 0; dy < D1Dy; ++dy)
5403 const double wy = Bct(dy,qy);
5404 const double wDy = Gct(dy,qy);
5406 for (int dz = 0; dz < D1Dz; ++dz)
5408 gradYZ01[dz][dy] += wy * massZ[dz][1];
5409 gradYZ10[dz][dy] += wDy * massZ[dz][0];
5414 for (int dx = 0; dx < D1Dx; ++dx)
5416 const double wx = Bct(dx,qx);
5417 const double wDx = Gct(dx,qx);
5419 for (int dy = 0; dy < D1Dy; ++dy)
5421 for (int dz = 0; dz < D1Dz; ++dz)
5423 // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
5424 // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
5425 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
5426 e) += (gradYZ10[dz][dy] * wx) - (gradYZ01[dz][dy] * wDx);
5435 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
5436 static void SmemPAHcurlL2Apply3DTranspose(const int D1D,
5440 const Array<double> &bo,
5441 const Array<double> &bc,
5442 const Array<double> &gc,
5443 const Vector &pa_data,
5447 MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D
");
5448 MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D
");
5450 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
5451 auto Bc = Reshape(bc.Read(), Q1D, D1D);
5452 auto Gc = Reshape(gc.Read(), Q1D, D1D);
5453 auto op = Reshape(pa_data.Read(), coeffDim, Q1D, Q1D, Q1D, NE);
5454 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
5455 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
5457 auto device_kernel = [=] MFEM_DEVICE (int e)
5459 constexpr int VDIM = 3;
5460 constexpr int maxCoeffDim = 9;
5462 MFEM_SHARED double sBo[MAX_D1D][MAX_Q1D];
5463 MFEM_SHARED double sBc[MAX_D1D][MAX_Q1D];
5464 MFEM_SHARED double sGc[MAX_D1D][MAX_Q1D];
5466 double opc[maxCoeffDim];
5467 MFEM_SHARED double sop[maxCoeffDim][MAX_Q1D][MAX_Q1D];
5468 MFEM_SHARED double mass[MAX_Q1D][MAX_Q1D][3];
5470 MFEM_SHARED double sX[MAX_D1D][MAX_D1D][MAX_D1D];
5472 MFEM_FOREACH_THREAD(qx,x,Q1D)
5474 MFEM_FOREACH_THREAD(qy,y,Q1D)
5476 MFEM_FOREACH_THREAD(qz,z,Q1D)
5478 for (int i=0; i<coeffDim; ++i)
5480 opc[i] = op(i,qx,qy,qz,e);
5486 const int tidx = MFEM_THREAD_ID(x);
5487 const int tidy = MFEM_THREAD_ID(y);
5488 const int tidz = MFEM_THREAD_ID(z);
5492 MFEM_FOREACH_THREAD(d,y,D1D)
5494 MFEM_FOREACH_THREAD(q,x,Q1D)
5496 sBc[d][q] = Bc(q,d);
5497 sGc[d][q] = Gc(q,d);
5500 sBo[d][q] = Bo(q,d);
5507 for (int qz=0; qz < Q1D; ++qz)
5511 MFEM_FOREACH_THREAD(qy,y,Q1D)
5513 MFEM_FOREACH_THREAD(qx,x,Q1D)
5515 for (int i=0; i<3; ++i)
5517 mass[qy][qx][i] = 0.0;
5524 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
5526 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
5527 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
5528 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
5530 MFEM_FOREACH_THREAD(dz,z,D1Dz)
5532 MFEM_FOREACH_THREAD(dy,y,D1Dy)
5534 MFEM_FOREACH_THREAD(dx,x,D1Dx)
5536 sX[dz][dy][dx] = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
5546 for (int i=0; i<coeffDim; ++i)
5548 sop[i][tidx][tidy] = opc[i];
5552 MFEM_FOREACH_THREAD(qy,y,Q1D)
5554 MFEM_FOREACH_THREAD(qx,x,Q1D)
5558 for (int dz = 0; dz < D1Dz; ++dz)
5560 const double wz = (c == 2) ? sBo[dz][qz] : sBc[dz][qz];
5562 for (int dy = 0; dy < D1Dy; ++dy)
5564 const double wy = (c == 1) ? sBo[dy][qy] : sBc[dy][qy];
5566 for (int dx = 0; dx < D1Dx; ++dx)
5568 const double wx = sX[dz][dy][dx] * ((c == 0) ? sBo[dx][qx] : sBc[dx][qx]);
5574 mass[qy][qx][c] += u;
5579 osc += D1Dx * D1Dy * D1Dz;
5587 MFEM_FOREACH_THREAD(dz,z,D1D)
5589 const double wcz = sBc[dz][qz];
5590 const double wcDz = sGc[dz][qz];
5591 const double wz = (dz < D1D-1) ? sBo[dz][qz] : 0.0;
5593 MFEM_FOREACH_THREAD(dy,y,D1D)
5595 MFEM_FOREACH_THREAD(dx,x,D1D)
5597 for (int qy = 0; qy < Q1D; ++qy)
5599 const double wcy = sBc[dy][qy];
5600 const double wcDy = sGc[dy][qy];
5601 const double wy = (dy < D1D-1) ? sBo[dy][qy] : 0.0;
5603 for (int qx = 0; qx < Q1D; ++qx)
5605 const double O11 = sop[0][qx][qy];
5609 c1 = O11 * mass[qy][qx][0];
5610 c2 = O11 * mass[qy][qx][1];
5611 c3 = O11 * mass[qy][qx][2];
5615 const double O12 = sop[1][qx][qy];
5616 const double O13 = sop[2][qx][qy];
5617 const double O21 = sop[3][qx][qy];
5618 const double O22 = sop[4][qx][qy];
5619 const double O23 = sop[5][qx][qy];
5620 const double O31 = sop[6][qx][qy];
5621 const double O32 = sop[7][qx][qy];
5622 const double O33 = sop[8][qx][qy];
5624 c1 = (O11*mass[qy][qx][0])+(O12*mass[qy][qx][1])+(O13*mass[qy][qx][2]);
5625 c2 = (O21*mass[qy][qx][0])+(O22*mass[qy][qx][1])+(O23*mass[qy][qx][2]);
5626 c3 = (O31*mass[qy][qx][0])+(O32*mass[qy][qx][1])+(O33*mass[qy][qx][2]);
5629 const double wcx = sBc[dx][qx];
5630 const double wDx = sGc[dx][qx];
5634 const double wx = sBo[dx][qx];
5635 dxyz1 += (wx * c2 * wcy * wcDz) - (wx * c3 * wcDy * wcz);
5638 dxyz2 += (-wy * c1 * wcx * wcDz) + (wy * c3 * wDx * wcz);
5640 dxyz3 += (wcDy * wz * c1 * wcx) - (wcy * wz * c2 * wDx);
5649 MFEM_FOREACH_THREAD(dz,z,D1D)
5651 MFEM_FOREACH_THREAD(dy,y,D1D)
5653 MFEM_FOREACH_THREAD(dx,x,D1D)
5657 Y(dx + ((dy + (dz * D1D)) * (D1D-1)), e) += dxyz1;
5661 Y(dx + ((dy + (dz * (D1D-1))) * D1D) + ((D1D-1)*D1D*D1D), e) += dxyz2;
5665 Y(dx + ((dy + (dz * D1D)) * D1D) + (2*(D1D-1)*D1D*D1D), e) += dxyz3;
5671 }; // end of element loop
5673 auto host_kernel = [&] MFEM_LAMBDA (int)
5675 MFEM_ABORT_KERNEL("This kernel should only be used on GPU.
");
5678 ForallWrap<3>(true, NE, device_kernel, host_kernel, Q1D, Q1D, Q1D);
5681 void MixedVectorWeakCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
5683 if (testType == mfem::FiniteElement::CURL &&
5684 trialType == mfem::FiniteElement::CURL && dim == 3)
5686 const int ndata = coeffDim == 1 ? 1 : 9;
5687 if (Device::Allows(Backend::DEVICE_MASK))
5689 const int ID = (dofs1D << 4) | quad1D;
5692 case 0x23: return SmemPAHcurlL2Apply3DTranspose<2,3>(dofs1D, quad1D, ndata,
5693 ne, mapsO->B, mapsC->B,
5694 mapsC->G, pa_data, x, y);
5695 case 0x34: return SmemPAHcurlL2Apply3DTranspose<3,4>(dofs1D, quad1D, ndata,
5696 ne, mapsO->B, mapsC->B,
5697 mapsC->G, pa_data, x, y);
5698 case 0x45: return SmemPAHcurlL2Apply3DTranspose<4,5>(dofs1D, quad1D, ndata,
5699 ne, mapsO->B, mapsC->B,
5700 mapsC->G, pa_data, x, y);
5701 case 0x56: return SmemPAHcurlL2Apply3DTranspose<5,6>(dofs1D, quad1D, ndata,
5702 ne, mapsO->B, mapsC->B,
5703 mapsC->G, pa_data, x, y);
5704 default: return SmemPAHcurlL2Apply3DTranspose(dofs1D, quad1D, ndata, ne,
5706 mapsC->G, pa_data, x, y);
5710 PAHcurlL2Apply3DTranspose(dofs1D, quad1D, ndata, ne, mapsO->B,
5711 mapsC->B, mapsO->Bt, mapsC->Bt, mapsC->Gt, pa_data, x, y);
5713 else if (testType == mfem::FiniteElement::CURL &&
5714 trialType == mfem::FiniteElement::DIV && dim == 3)
5716 PAHcurlHdivApply3DTranspose(dofs1D, dofs1D, quad1D, ne, mapsO->B,
5717 mapsC->B, mapsO->Bt, mapsC->Bt,
5718 mapsC->Gt, pa_data, x, y);
5726 void MixedVectorWeakCurlIntegrator::AddMultTransposePA(const Vector &x,
5729 if (testType == mfem::FiniteElement::CURL &&
5730 trialType == mfem::FiniteElement::DIV && dim == 3)
5732 PAHcurlHdivApply3D(dofs1D, dofs1D, quad1D, ne, mapsO->B,
5733 mapsC->B, mapsO->Bt, mapsC->Bt, mapsC->G,
5742 // Apply to x corresponding to DOFs in H^1 (domain) the (topological) gradient
5743 // to get a dof in H(curl) (range). You can think of the range as the "test
" space
5744 // and the domain as the "trial
" space, but there's no integration.
5745 static void PAHcurlApplyGradient2D(const int c_dofs1D,
5748 const Array<double> &B_,
5749 const Array<double> &G_,
5753 auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
5754 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
5756 auto x = Reshape(x_.Read(), c_dofs1D, c_dofs1D, NE);
5757 auto y = Reshape(y_.ReadWrite(), 2 * c_dofs1D * o_dofs1D, NE);
5759 constexpr static int MAX_D1D = HCURL_MAX_D1D;
5760 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
5764 double w[MAX_D1D][MAX_D1D];
5767 for (int dx = 0; dx < c_dofs1D; ++dx)
5769 for (int ey = 0; ey < c_dofs1D; ++ey)
5772 for (int dy = 0; dy < c_dofs1D; ++dy)
5774 w[dx][ey] += B(ey, dy) * x(dx, dy, e);
5779 for (int ey = 0; ey < c_dofs1D; ++ey)
5781 for (int ex = 0; ex < o_dofs1D; ++ex)
5784 for (int dx = 0; dx < c_dofs1D; ++dx)
5786 s += G(ex, dx) * w[dx][ey];
5788 const int local_index = ey*o_dofs1D + ex;
5789 y(local_index, e) += s;
5794 for (int dx = 0; dx < c_dofs1D; ++dx)
5796 for (int ey = 0; ey < o_dofs1D; ++ey)
5799 for (int dy = 0; dy < c_dofs1D; ++dy)
5801 w[dx][ey] += G(ey, dy) * x(dx, dy, e);
5806 for (int ey = 0; ey < o_dofs1D; ++ey)
5808 for (int ex = 0; ex < c_dofs1D; ++ex)
5811 for (int dx = 0; dx < c_dofs1D; ++dx)
5813 s += B(ex, dx) * w[dx][ey];
5815 const int local_index = c_dofs1D * o_dofs1D + ey*c_dofs1D + ex;
5816 y(local_index, e) += s;
5822 // Specialization of PAHcurlApplyGradient2D to the case where B is identity
5823 static void PAHcurlApplyGradient2DBId(const int c_dofs1D,
5826 const Array<double> &G_,
5830 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
5832 auto x = Reshape(x_.Read(), c_dofs1D, c_dofs1D, NE);
5833 auto y = Reshape(y_.ReadWrite(), 2 * c_dofs1D * o_dofs1D, NE);
5835 constexpr static int MAX_D1D = HCURL_MAX_D1D;
5836 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
5840 double w[MAX_D1D][MAX_D1D];
5843 for (int dx = 0; dx < c_dofs1D; ++dx)
5845 for (int ey = 0; ey < c_dofs1D; ++ey)
5848 w[dx][ey] = x(dx, dy, e);
5852 for (int ey = 0; ey < c_dofs1D; ++ey)
5854 for (int ex = 0; ex < o_dofs1D; ++ex)
5857 for (int dx = 0; dx < c_dofs1D; ++dx)
5859 s += G(ex, dx) * w[dx][ey];
5861 const int local_index = ey*o_dofs1D + ex;
5862 y(local_index, e) += s;
5867 for (int dx = 0; dx < c_dofs1D; ++dx)
5869 for (int ey = 0; ey < o_dofs1D; ++ey)
5872 for (int dy = 0; dy < c_dofs1D; ++dy)
5874 w[dx][ey] += G(ey, dy) * x(dx, dy, e);
5879 for (int ey = 0; ey < o_dofs1D; ++ey)
5881 for (int ex = 0; ex < c_dofs1D; ++ex)
5884 const double s = w[dx][ey];
5885 const int local_index = c_dofs1D * o_dofs1D + ey*c_dofs1D + ex;
5886 y(local_index, e) += s;
5892 static void PAHcurlApplyGradientTranspose2D(
5893 const int c_dofs1D, const int o_dofs1D, const int NE,
5894 const Array<double> &B_, const Array<double> &G_,
5895 const Vector &x_, Vector &y_)
5897 auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
5898 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
5900 auto x = Reshape(x_.Read(), 2 * c_dofs1D * o_dofs1D, NE);
5901 auto y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, NE);
5903 constexpr static int MAX_D1D = HCURL_MAX_D1D;
5904 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
5908 double w[MAX_D1D][MAX_D1D];
5910 // horizontal part (open x, closed y)
5911 for (int dy = 0; dy < c_dofs1D; ++dy)
5913 for (int ex = 0; ex < o_dofs1D; ++ex)
5916 for (int ey = 0; ey < c_dofs1D; ++ey)
5918 const int local_index = ey*o_dofs1D + ex;
5919 w[dy][ex] += B(ey, dy) * x(local_index, e);
5924 for (int dy = 0; dy < c_dofs1D; ++dy)
5926 for (int dx = 0; dx < c_dofs1D; ++dx)
5929 for (int ex = 0; ex < o_dofs1D; ++ex)
5931 s += G(ex, dx) * w[dy][ex];
5937 // vertical part (open y, closed x)
5938 for (int dy = 0; dy < c_dofs1D; ++dy)
5940 for (int ex = 0; ex < c_dofs1D; ++ex)
5943 for (int ey = 0; ey < o_dofs1D; ++ey)
5945 const int local_index = c_dofs1D * o_dofs1D + ey*c_dofs1D + ex;
5946 w[dy][ex] += G(ey, dy) * x(local_index, e);
5951 for (int dy = 0; dy < c_dofs1D; ++dy)
5953 for (int dx = 0; dx < c_dofs1D; ++dx)
5956 for (int ex = 0; ex < c_dofs1D; ++ex)
5958 s += B(ex, dx) * w[dy][ex];
5966 // Specialization of PAHcurlApplyGradientTranspose2D to the case where
5968 static void PAHcurlApplyGradientTranspose2DBId(
5969 const int c_dofs1D, const int o_dofs1D, const int NE,
5970 const Array<double> &G_,
5971 const Vector &x_, Vector &y_)
5973 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
5975 auto x = Reshape(x_.Read(), 2 * c_dofs1D * o_dofs1D, NE);
5976 auto y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, NE);
5978 constexpr static int MAX_D1D = HCURL_MAX_D1D;
5979 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
5983 double w[MAX_D1D][MAX_D1D];
5985 // horizontal part (open x, closed y)
5986 for (int dy = 0; dy < c_dofs1D; ++dy)
5988 for (int ex = 0; ex < o_dofs1D; ++ex)
5991 const int local_index = ey*o_dofs1D + ex;
5992 w[dy][ex] = x(local_index, e);
5996 for (int dy = 0; dy < c_dofs1D; ++dy)
5998 for (int dx = 0; dx < c_dofs1D; ++dx)
6001 for (int ex = 0; ex < o_dofs1D; ++ex)
6003 s += G(ex, dx) * w[dy][ex];
6009 // vertical part (open y, closed x)
6010 for (int dy = 0; dy < c_dofs1D; ++dy)
6012 for (int ex = 0; ex < c_dofs1D; ++ex)
6015 for (int ey = 0; ey < o_dofs1D; ++ey)
6017 const int local_index = c_dofs1D * o_dofs1D + ey*c_dofs1D + ex;
6018 w[dy][ex] += G(ey, dy) * x(local_index, e);
6023 for (int dy = 0; dy < c_dofs1D; ++dy)
6025 for (int dx = 0; dx < c_dofs1D; ++dx)
6028 const double s = w[dy][ex];
6035 static void PAHcurlApplyGradient3D(const int c_dofs1D,
6038 const Array<double> &B_,
6039 const Array<double> &G_,
6043 auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
6044 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
6046 auto x = Reshape(x_.Read(), c_dofs1D, c_dofs1D, c_dofs1D, NE);
6047 auto y = Reshape(y_.ReadWrite(), (3 * c_dofs1D * c_dofs1D * o_dofs1D), NE);
6049 constexpr static int MAX_D1D = HCURL_MAX_D1D;
6050 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
6054 double w1[MAX_D1D][MAX_D1D][MAX_D1D];
6055 double w2[MAX_D1D][MAX_D1D][MAX_D1D];
6058 // dofs that point parallel to x-axis (open in x, closed in y, z)
6062 for (int ez = 0; ez < c_dofs1D; ++ez)
6064 for (int dx = 0; dx < c_dofs1D; ++dx)
6066 for (int dy = 0; dy < c_dofs1D; ++dy)
6068 w1[dx][dy][ez] = 0.0;
6069 for (int dz = 0; dz < c_dofs1D; ++dz)
6071 w1[dx][dy][ez] += B(ez, dz) * x(dx, dy, dz, e);
6078 for (int ez = 0; ez < c_dofs1D; ++ez)
6080 for (int ey = 0; ey < c_dofs1D; ++ey)
6082 for (int dx = 0; dx < c_dofs1D; ++dx)
6084 w2[dx][ey][ez] = 0.0;
6085 for (int dy = 0; dy < c_dofs1D; ++dy)
6087 w2[dx][ey][ez] += B(ey, dy) * w1[dx][dy][ez];
6094 for (int ez = 0; ez < c_dofs1D; ++ez)
6096 for (int ey = 0; ey < c_dofs1D; ++ey)
6098 for (int ex = 0; ex < o_dofs1D; ++ex)
6101 for (int dx = 0; dx < c_dofs1D; ++dx)
6103 s += G(ex, dx) * w2[dx][ey][ez];
6105 const int local_index = ez*c_dofs1D*o_dofs1D + ey*o_dofs1D + ex;
6106 y(local_index, e) += s;
6112 // dofs that point parallel to y-axis (open in y, closed in x, z)
6116 for (int ez = 0; ez < c_dofs1D; ++ez)
6118 for (int dx = 0; dx < c_dofs1D; ++dx)
6120 for (int dy = 0; dy < c_dofs1D; ++dy)
6122 w1[dx][dy][ez] = 0.0;
6123 for (int dz = 0; dz < c_dofs1D; ++dz)
6125 w1[dx][dy][ez] += B(ez, dz) * x(dx, dy, dz, e);
6132 for (int ez = 0; ez < c_dofs1D; ++ez)
6134 for (int ey = 0; ey < o_dofs1D; ++ey)
6136 for (int dx = 0; dx < c_dofs1D; ++dx)
6138 w2[dx][ey][ez] = 0.0;
6139 for (int dy = 0; dy < c_dofs1D; ++dy)
6141 w2[dx][ey][ez] += G(ey, dy) * w1[dx][dy][ez];
6148 for (int ez = 0; ez < c_dofs1D; ++ez)
6150 for (int ey = 0; ey < o_dofs1D; ++ey)
6152 for (int ex = 0; ex < c_dofs1D; ++ex)
6155 for (int dx = 0; dx < c_dofs1D; ++dx)
6157 s += B(ex, dx) * w2[dx][ey][ez];
6159 const int local_index = c_dofs1D*c_dofs1D*o_dofs1D +
6160 ez*c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
6161 y(local_index, e) += s;
6167 // dofs that point parallel to z-axis (open in z, closed in x, y)
6171 for (int ez = 0; ez < o_dofs1D; ++ez)
6173 for (int dx = 0; dx < c_dofs1D; ++dx)
6175 for (int dy = 0; dy < c_dofs1D; ++dy)
6177 w1[dx][dy][ez] = 0.0;
6178 for (int dz = 0; dz < c_dofs1D; ++dz)
6180 w1[dx][dy][ez] += G(ez, dz) * x(dx, dy, dz, e);
6187 for (int ez = 0; ez < o_dofs1D; ++ez)
6189 for (int ey = 0; ey < c_dofs1D; ++ey)
6191 for (int dx = 0; dx < c_dofs1D; ++dx)
6193 w2[dx][ey][ez] = 0.0;
6194 for (int dy = 0; dy < c_dofs1D; ++dy)
6196 w2[dx][ey][ez] += B(ey, dy) * w1[dx][dy][ez];
6203 for (int ez = 0; ez < o_dofs1D; ++ez)
6205 for (int ey = 0; ey < c_dofs1D; ++ey)
6207 for (int ex = 0; ex < c_dofs1D; ++ex)
6210 for (int dx = 0; dx < c_dofs1D; ++dx)
6212 s += B(ex, dx) * w2[dx][ey][ez];
6214 const int local_index = 2*c_dofs1D*c_dofs1D*o_dofs1D +
6215 ez*c_dofs1D*c_dofs1D + ey*c_dofs1D + ex;
6216 y(local_index, e) += s;
6223 // Specialization of PAHcurlApplyGradient3D to the case where
6224 static void PAHcurlApplyGradient3DBId(const int c_dofs1D,
6227 const Array<double> &G_,
6231 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
6233 auto x = Reshape(x_.Read(), c_dofs1D, c_dofs1D, c_dofs1D, NE);
6234 auto y = Reshape(y_.ReadWrite(), (3 * c_dofs1D * c_dofs1D * o_dofs1D), NE);
6236 constexpr static int MAX_D1D = HCURL_MAX_D1D;
6237 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
6241 double w1[MAX_D1D][MAX_D1D][MAX_D1D];
6242 double w2[MAX_D1D][MAX_D1D][MAX_D1D];
6245 // dofs that point parallel to x-axis (open in x, closed in y, z)
6249 for (int ez = 0; ez < c_dofs1D; ++ez)
6251 for (int dx = 0; dx < c_dofs1D; ++dx)
6253 for (int dy = 0; dy < c_dofs1D; ++dy)
6256 w1[dx][dy][ez] = x(dx, dy, dz, e);
6262 for (int ez = 0; ez < c_dofs1D; ++ez)
6264 for (int ey = 0; ey < c_dofs1D; ++ey)
6266 for (int dx = 0; dx < c_dofs1D; ++dx)
6269 w2[dx][ey][ez] = w1[dx][dy][ez];
6275 for (int ez = 0; ez < c_dofs1D; ++ez)
6277 for (int ey = 0; ey < c_dofs1D; ++ey)
6279 for (int ex = 0; ex < o_dofs1D; ++ex)
6282 for (int dx = 0; dx < c_dofs1D; ++dx)
6284 s += G(ex, dx) * w2[dx][ey][ez];
6286 const int local_index = ez*c_dofs1D*o_dofs1D + ey*o_dofs1D + ex;
6287 y(local_index, e) += s;
6293 // dofs that point parallel to y-axis (open in y, closed in x, z)
6297 for (int ez = 0; ez < c_dofs1D; ++ez)
6299 for (int dx = 0; dx < c_dofs1D; ++dx)
6301 for (int dy = 0; dy < c_dofs1D; ++dy)
6304 w1[dx][dy][ez] = x(dx, dy, dz, e);
6310 for (int ez = 0; ez < c_dofs1D; ++ez)
6312 for (int ey = 0; ey < o_dofs1D; ++ey)
6314 for (int dx = 0; dx < c_dofs1D; ++dx)
6316 w2[dx][ey][ez] = 0.0;
6317 for (int dy = 0; dy < c_dofs1D; ++dy)
6319 w2[dx][ey][ez] += G(ey, dy) * w1[dx][dy][ez];
6326 for (int ez = 0; ez < c_dofs1D; ++ez)
6328 for (int ey = 0; ey < o_dofs1D; ++ey)
6330 for (int ex = 0; ex < c_dofs1D; ++ex)
6333 const double s = w2[dx][ey][ez];
6334 const int local_index = c_dofs1D*c_dofs1D*o_dofs1D +
6335 ez*c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
6336 y(local_index, e) += s;
6342 // dofs that point parallel to z-axis (open in z, closed in x, y)
6346 for (int ez = 0; ez < o_dofs1D; ++ez)
6348 for (int dx = 0; dx < c_dofs1D; ++dx)
6350 for (int dy = 0; dy < c_dofs1D; ++dy)
6352 w1[dx][dy][ez] = 0.0;
6353 for (int dz = 0; dz < c_dofs1D; ++dz)
6355 w1[dx][dy][ez] += G(ez, dz) * x(dx, dy, dz, e);
6362 for (int ez = 0; ez < o_dofs1D; ++ez)
6364 for (int ey = 0; ey < c_dofs1D; ++ey)
6366 for (int dx = 0; dx < c_dofs1D; ++dx)
6369 w2[dx][ey][ez] = w1[dx][dy][ez];
6375 for (int ez = 0; ez < o_dofs1D; ++ez)
6377 for (int ey = 0; ey < c_dofs1D; ++ey)
6379 for (int ex = 0; ex < c_dofs1D; ++ex)
6382 const double s = w2[dx][ey][ez];
6383 const int local_index = 2*c_dofs1D*c_dofs1D*o_dofs1D +
6384 ez*c_dofs1D*c_dofs1D + ey*c_dofs1D + ex;
6385 y(local_index, e) += s;
6392 static void PAHcurlApplyGradientTranspose3D(
6393 const int c_dofs1D, const int o_dofs1D, const int NE,
6394 const Array<double> &B_, const Array<double> &G_,
6395 const Vector &x_, Vector &y_)
6397 auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
6398 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
6400 auto x = Reshape(x_.Read(), (3 * c_dofs1D * c_dofs1D * o_dofs1D), NE);
6401 auto y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, c_dofs1D, NE);
6403 constexpr static int MAX_D1D = HCURL_MAX_D1D;
6404 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
6408 double w1[MAX_D1D][MAX_D1D][MAX_D1D];
6409 double w2[MAX_D1D][MAX_D1D][MAX_D1D];
6411 // dofs that point parallel to x-axis (open in x, closed in y, z)
6415 for (int dz = 0; dz < c_dofs1D; ++dz)
6417 for (int ex = 0; ex < o_dofs1D; ++ex)
6419 for (int ey = 0; ey < c_dofs1D; ++ey)
6421 w1[ex][ey][dz] = 0.0;
6422 for (int ez = 0; ez < c_dofs1D; ++ez)
6424 const int local_index = ez*c_dofs1D*o_dofs1D + ey*o_dofs1D + ex;
6425 w1[ex][ey][dz] += B(ez, dz) * x(local_index, e);
6432 for (int dz = 0; dz < c_dofs1D; ++dz)
6434 for (int dy = 0; dy < c_dofs1D; ++dy)
6436 for (int ex = 0; ex < o_dofs1D; ++ex)
6438 w2[ex][dy][dz] = 0.0;
6439 for (int ey = 0; ey < c_dofs1D; ++ey)
6441 w2[ex][dy][dz] += B(ey, dy) * w1[ex][ey][dz];
6448 for (int dz = 0; dz < c_dofs1D; ++dz)
6450 for (int dy = 0; dy < c_dofs1D; ++dy)
6452 for (int dx = 0; dx < c_dofs1D; ++dx)
6455 for (int ex = 0; ex < o_dofs1D; ++ex)
6457 s += G(ex, dx) * w2[ex][dy][dz];
6459 y(dx, dy, dz, e) += s;
6465 // dofs that point parallel to y-axis (open in y, closed in x, z)
6469 for (int dz = 0; dz < c_dofs1D; ++dz)
6471 for (int ex = 0; ex < c_dofs1D; ++ex)
6473 for (int ey = 0; ey < o_dofs1D; ++ey)
6475 w1[ex][ey][dz] = 0.0;
6476 for (int ez = 0; ez < c_dofs1D; ++ez)
6478 const int local_index = c_dofs1D*c_dofs1D*o_dofs1D +
6479 ez*c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
6480 w1[ex][ey][dz] += B(ez, dz) * x(local_index, e);
6487 for (int dz = 0; dz < c_dofs1D; ++dz)
6489 for (int dy = 0; dy < c_dofs1D; ++dy)
6491 for (int ex = 0; ex < c_dofs1D; ++ex)
6493 w2[ex][dy][dz] = 0.0;
6494 for (int ey = 0; ey < o_dofs1D; ++ey)
6496 w2[ex][dy][dz] += G(ey, dy) * w1[ex][ey][dz];
6503 for (int dz = 0; dz < c_dofs1D; ++dz)
6505 for (int dy = 0; dy < c_dofs1D; ++dy)
6507 for (int dx = 0; dx < c_dofs1D; ++dx)
6510 for (int ex = 0; ex < c_dofs1D; ++ex)
6512 s += B(ex, dx) * w2[ex][dy][dz];
6514 y(dx, dy, dz, e) += s;
6520 // dofs that point parallel to z-axis (open in z, closed in x, y)
6524 for (int dz = 0; dz < c_dofs1D; ++dz)
6526 for (int ex = 0; ex < c_dofs1D; ++ex)
6528 for (int ey = 0; ey < c_dofs1D; ++ey)
6530 w1[ex][ey][dz] = 0.0;
6531 for (int ez = 0; ez < o_dofs1D; ++ez)
6533 const int local_index = 2*c_dofs1D*c_dofs1D*o_dofs1D +
6534 ez*c_dofs1D*c_dofs1D + ey*c_dofs1D + ex;
6535 w1[ex][ey][dz] += G(ez, dz) * x(local_index, e);
6542 for (int dz = 0; dz < c_dofs1D; ++dz)
6544 for (int dy = 0; dy < c_dofs1D; ++dy)
6546 for (int ex = 0; ex < c_dofs1D; ++ex)
6548 w2[ex][dy][dz] = 0.0;
6549 for (int ey = 0; ey < c_dofs1D; ++ey)
6551 w2[ex][dy][dz] += B(ey, dy) * w1[ex][ey][dz];
6558 for (int dz = 0; dz < c_dofs1D; ++dz)
6560 for (int dy = 0; dy < c_dofs1D; ++dy)
6562 for (int dx = 0; dx < c_dofs1D; ++dx)
6565 for (int ex = 0; ex < c_dofs1D; ++ex)
6567 s += B(ex, dx) * w2[ex][dy][dz];
6569 y(dx, dy, dz, e) += s;
6576 // Specialization of PAHcurlApplyGradientTranspose3D to the case where
6577 static void PAHcurlApplyGradientTranspose3DBId(
6578 const int c_dofs1D, const int o_dofs1D, const int NE,
6579 const Array<double> &G_,
6580 const Vector &x_, Vector &y_)
6582 auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
6584 auto x = Reshape(x_.Read(), (3 * c_dofs1D * c_dofs1D * o_dofs1D), NE);
6585 auto y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, c_dofs1D, NE);
6587 constexpr static int MAX_D1D = HCURL_MAX_D1D;
6588 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
6592 double w1[MAX_D1D][MAX_D1D][MAX_D1D];
6593 double w2[MAX_D1D][MAX_D1D][MAX_D1D];
6595 // dofs that point parallel to x-axis (open in x, closed in y, z)
6599 for (int dz = 0; dz < c_dofs1D; ++dz)
6601 for (int ex = 0; ex < o_dofs1D; ++ex)
6603 for (int ey = 0; ey < c_dofs1D; ++ey)
6606 const int local_index = ez*c_dofs1D*o_dofs1D + ey*o_dofs1D + ex;
6607 w1[ex][ey][dz] = x(local_index, e);
6613 for (int dz = 0; dz < c_dofs1D; ++dz)
6615 for (int dy = 0; dy < c_dofs1D; ++dy)
6617 for (int ex = 0; ex < o_dofs1D; ++ex)
6620 w2[ex][dy][dz] = w1[ex][ey][dz];
6626 for (int dz = 0; dz < c_dofs1D; ++dz)
6628 for (int dy = 0; dy < c_dofs1D; ++dy)
6630 for (int dx = 0; dx < c_dofs1D; ++dx)
6633 for (int ex = 0; ex < o_dofs1D; ++ex)
6635 s += G(ex, dx) * w2[ex][dy][dz];
6637 y(dx, dy, dz, e) += s;
6643 // dofs that point parallel to y-axis (open in y, closed in x, z)
6647 for (int dz = 0; dz < c_dofs1D; ++dz)
6649 for (int ex = 0; ex < c_dofs1D; ++ex)
6651 for (int ey = 0; ey < o_dofs1D; ++ey)
6654 const int local_index = c_dofs1D*c_dofs1D*o_dofs1D +
6655 ez*c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
6656 w1[ex][ey][dz] = x(local_index, e);
6662 for (int dz = 0; dz < c_dofs1D; ++dz)
6664 for (int dy = 0; dy < c_dofs1D; ++dy)
6666 for (int ex = 0; ex < c_dofs1D; ++ex)
6668 w2[ex][dy][dz] = 0.0;
6669 for (int ey = 0; ey < o_dofs1D; ++ey)
6671 w2[ex][dy][dz] += G(ey, dy) * w1[ex][ey][dz];
6678 for (int dz = 0; dz < c_dofs1D; ++dz)
6680 for (int dy = 0; dy < c_dofs1D; ++dy)
6682 for (int dx = 0; dx < c_dofs1D; ++dx)
6685 double s = w2[ex][dy][dz];
6686 y(dx, dy, dz, e) += s;
6692 // dofs that point parallel to z-axis (open in z, closed in x, y)
6696 for (int dz = 0; dz < c_dofs1D; ++dz)
6698 for (int ex = 0; ex < c_dofs1D; ++ex)
6700 for (int ey = 0; ey < c_dofs1D; ++ey)
6702 w1[ex][ey][dz] = 0.0;
6703 for (int ez = 0; ez < o_dofs1D; ++ez)
6705 const int local_index = 2*c_dofs1D*c_dofs1D*o_dofs1D +
6706 ez*c_dofs1D*c_dofs1D + ey*c_dofs1D + ex;
6707 w1[ex][ey][dz] += G(ez, dz) * x(local_index, e);
6714 for (int dz = 0; dz < c_dofs1D; ++dz)
6716 for (int dy = 0; dy < c_dofs1D; ++dy)
6718 for (int ex = 0; ex < c_dofs1D; ++ex)
6721 w2[ex][dy][dz] = w1[ex][ey][dz];
6727 for (int dz = 0; dz < c_dofs1D; ++dz)
6729 for (int dy = 0; dy < c_dofs1D; ++dy)
6731 for (int dx = 0; dx < c_dofs1D; ++dx)
6734 double s = w2[ex][dy][dz];
6735 y(dx, dy, dz, e) += s;
6742 void GradientInterpolator::AssemblePA(const FiniteElementSpace &trial_fes,
6743 const FiniteElementSpace &test_fes)
6745 // Assumes tensor-product elements, with a vector test space and H^1 trial space.
6746 Mesh *mesh = trial_fes.GetMesh();
6747 const FiniteElement *trial_fel = trial_fes.GetFE(0);
6748 const FiniteElement *test_fel = test_fes.GetFE(0);
6750 const NodalTensorFiniteElement *trial_el =
6751 dynamic_cast<const NodalTensorFiniteElement*>(trial_fel);
6754 const VectorTensorFiniteElement *test_el =
6755 dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
6758 const int dims = trial_el->GetDim();
6759 MFEM_VERIFY(dims == 2 || dims == 3, "Bad
dimension!
");
6760 dim = mesh->Dimension();
6761 MFEM_VERIFY(dim == 2 || dim == 3, "Bad
dimension!
");
6762 MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(),
6763 "Orders
do not match!
");
6764 ne = trial_fes.GetNE();
6766 const int order = trial_el->GetOrder();
6767 dofquad_fe = new H1_SegmentElement(order, trial_el->GetBasisType());
6768 mfem::QuadratureFunctions1D qf1d;
6769 mfem::IntegrationRule closed_ir;
6770 closed_ir.SetSize(order + 1);
6771 qf1d.GaussLobatto(order + 1, &closed_ir);
6772 mfem::IntegrationRule open_ir;
6773 open_ir.SetSize(order);
6774 qf1d.GaussLegendre(order, &open_ir);
6776 maps_O_C = &dofquad_fe->GetDofToQuad(open_ir, DofToQuad::TENSOR);
6777 o_dofs1D = maps_O_C->nqpt;
6778 if (trial_el->GetBasisType() == BasisType::GaussLobatto)
6781 c_dofs1D = maps_O_C->ndof;
6786 maps_C_C = &dofquad_fe->GetDofToQuad(closed_ir, DofToQuad::TENSOR);
6787 c_dofs1D = maps_C_C->nqpt;
6791 void GradientInterpolator::AddMultPA(const Vector &x, Vector &y) const
6797 PAHcurlApplyGradient3DBId(c_dofs1D, o_dofs1D, ne,
6802 PAHcurlApplyGradient3D(c_dofs1D, o_dofs1D, ne, maps_C_C->B,
6810 PAHcurlApplyGradient2DBId(c_dofs1D, o_dofs1D, ne,
6815 PAHcurlApplyGradient2D(c_dofs1D, o_dofs1D, ne, maps_C_C->B, maps_O_C->G,
6825 void GradientInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
6831 PAHcurlApplyGradientTranspose3DBId(c_dofs1D, o_dofs1D, ne,
6836 PAHcurlApplyGradientTranspose3D(c_dofs1D, o_dofs1D, ne, maps_C_C->B,
6844 PAHcurlApplyGradientTranspose2DBId(c_dofs1D, o_dofs1D, ne,
6849 PAHcurlApplyGradientTranspose2D(c_dofs1D, o_dofs1D, ne, maps_C_C->B,
6859 static void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
6862 const Array<double> &Bclosed,
6863 const Array<double> &Bopen,
6864 const Vector &pa_data,
6868 auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
6869 auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
6871 auto x = Reshape(x_.Read(), c_dofs1D, c_dofs1D, c_dofs1D, 3, NE);
6872 auto y = Reshape(y_.ReadWrite(), (3 * c_dofs1D * c_dofs1D * o_dofs1D), NE);
6874 auto vk = Reshape(pa_data.Read(), 3, (3 * c_dofs1D * c_dofs1D * o_dofs1D),
6877 constexpr static int MAX_D1D = HCURL_MAX_D1D;
6878 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
6882 double w1[3][MAX_D1D][MAX_D1D][MAX_D1D];
6883 double w2[3][MAX_D1D][MAX_D1D][MAX_D1D];
6885 // dofs that point parallel to x-axis (open in x, closed in y, z)
6888 for (int ez = 0; ez < c_dofs1D; ++ez)
6890 for (int dx = 0; dx < c_dofs1D; ++dx)
6892 for (int dy = 0; dy < c_dofs1D; ++dy)
6894 for (int j=0; j<3; ++j)
6896 w1[j][dx][dy][ez] = 0.0;
6897 for (int dz = 0; dz < c_dofs1D; ++dz)
6899 w1[j][dx][dy][ez] += Bc(ez, dz) * x(dx, dy, dz, j, e);
6907 for (int ez = 0; ez < c_dofs1D; ++ez)
6909 for (int ey = 0; ey < c_dofs1D; ++ey)
6911 for (int dx = 0; dx < c_dofs1D; ++dx)
6913 for (int j=0; j<3; ++j)
6915 w2[j][dx][ey][ez] = 0.0;
6916 for (int dy = 0; dy < c_dofs1D; ++dy)
6918 w2[j][dx][ey][ez] += Bc(ey, dy) * w1[j][dx][dy][ez];
6926 for (int ez = 0; ez < c_dofs1D; ++ez)
6928 for (int ey = 0; ey < c_dofs1D; ++ey)
6930 for (int ex = 0; ex < o_dofs1D; ++ex)
6932 for (int j=0; j<3; ++j)
6935 for (int dx = 0; dx < c_dofs1D; ++dx)
6937 s += Bo(ex, dx) * w2[j][dx][ey][ez];
6939 const int local_index = ez*c_dofs1D*o_dofs1D + ey*o_dofs1D + ex;
6940 y(local_index, e) += s * vk(j, local_index, e);
6946 // dofs that point parallel to y-axis (open in y, closed in x, z)
6949 for (int ez = 0; ez < c_dofs1D; ++ez)
6951 for (int dx = 0; dx < c_dofs1D; ++dx)
6953 for (int dy = 0; dy < c_dofs1D; ++dy)
6955 for (int j=0; j<3; ++j)
6957 w1[j][dx][dy][ez] = 0.0;
6958 for (int dz = 0; dz < c_dofs1D; ++dz)
6960 w1[j][dx][dy][ez] += Bc(ez, dz) * x(dx, dy, dz, j, e);
6968 for (int ez = 0; ez < c_dofs1D; ++ez)
6970 for (int ey = 0; ey < o_dofs1D; ++ey)
6972 for (int dx = 0; dx < c_dofs1D; ++dx)
6974 for (int j=0; j<3; ++j)
6976 w2[j][dx][ey][ez] = 0.0;
6977 for (int dy = 0; dy < c_dofs1D; ++dy)
6979 w2[j][dx][ey][ez] += Bo(ey, dy) * w1[j][dx][dy][ez];
6987 for (int ez = 0; ez < c_dofs1D; ++ez)
6989 for (int ey = 0; ey < o_dofs1D; ++ey)
6991 for (int ex = 0; ex < c_dofs1D; ++ex)
6993 for (int j=0; j<3; ++j)
6996 for (int dx = 0; dx < c_dofs1D; ++dx)
6998 s += Bc(ex, dx) * w2[j][dx][ey][ez];
7000 const int local_index = c_dofs1D*c_dofs1D*o_dofs1D +
7001 ez*c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
7002 y(local_index, e) += s * vk(j, local_index, e);
7008 // dofs that point parallel to z-axis (open in z, closed in x, y)
7011 for (int ez = 0; ez < o_dofs1D; ++ez)
7013 for (int dx = 0; dx < c_dofs1D; ++dx)
7015 for (int dy = 0; dy < c_dofs1D; ++dy)
7017 for (int j=0; j<3; ++j)
7019 w1[j][dx][dy][ez] = 0.0;
7020 for (int dz = 0; dz < c_dofs1D; ++dz)
7022 w1[j][dx][dy][ez] += Bo(ez, dz) * x(dx, dy, dz, j, e);
7030 for (int ez = 0; ez < o_dofs1D; ++ez)
7032 for (int ey = 0; ey < c_dofs1D; ++ey)
7034 for (int dx = 0; dx < c_dofs1D; ++dx)
7036 for (int j=0; j<3; ++j)
7038 w2[j][dx][ey][ez] = 0.0;
7039 for (int dy = 0; dy < c_dofs1D; ++dy)
7041 w2[j][dx][ey][ez] += Bc(ey, dy) * w1[j][dx][dy][ez];
7049 for (int ez = 0; ez < o_dofs1D; ++ez)
7051 for (int ey = 0; ey < c_dofs1D; ++ey)
7053 for (int ex = 0; ex < c_dofs1D; ++ex)
7055 for (int j=0; j<3; ++j)
7058 for (int dx = 0; dx < c_dofs1D; ++dx)
7060 s += Bc(ex, dx) * w2[j][dx][ey][ez];
7062 const int local_index = 2*c_dofs1D*c_dofs1D*o_dofs1D +
7063 ez*c_dofs1D*c_dofs1D + ey*c_dofs1D + ex;
7064 y(local_index, e) += s * vk(j, local_index, e);
7072 static void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
7075 const Array<double> &Bclosed,
7076 const Array<double> &Bopen,
7077 const Vector &pa_data,
7081 auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
7082 auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
7084 auto x = Reshape(x_.Read(), (3 * c_dofs1D * c_dofs1D * o_dofs1D), NE);
7085 auto y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, c_dofs1D, 3, NE);
7087 auto vk = Reshape(pa_data.Read(), 3, (3 * c_dofs1D * c_dofs1D * o_dofs1D),
7090 constexpr static int MAX_D1D = HCURL_MAX_D1D;
7092 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
7096 double w1[3][MAX_D1D][MAX_D1D][MAX_D1D];
7097 double w2[3][MAX_D1D][MAX_D1D][MAX_D1D];
7099 // dofs that point parallel to x-axis (open in x, closed in y, z)
7102 for (int ez = 0; ez < c_dofs1D; ++ez)
7104 for (int ey = 0; ey < c_dofs1D; ++ey)
7106 for (int j=0; j<3; ++j)
7108 for (int dx = 0; dx < c_dofs1D; ++dx)
7110 w2[j][dx][ey][ez] = 0.0;
7112 for (int ex = 0; ex < o_dofs1D; ++ex)
7114 const int local_index = ez*c_dofs1D*o_dofs1D + ey*o_dofs1D + ex;
7115 const double xv = x(local_index, e) * vk(j, local_index, e);
7116 for (int dx = 0; dx < c_dofs1D; ++dx)
7118 w2[j][dx][ey][ez] += xv * Bo(ex, dx);
7126 for (int ez = 0; ez < c_dofs1D; ++ez)
7128 for (int dx = 0; dx < c_dofs1D; ++dx)
7130 for (int dy = 0; dy < c_dofs1D; ++dy)
7132 for (int j=0; j<3; ++j)
7134 w1[j][dx][dy][ez] = 0.0;
7135 for (int ey = 0; ey < c_dofs1D; ++ey)
7137 w1[j][dx][dy][ez] += w2[j][dx][ey][ez] * Bc(ey, dy);
7145 for (int dx = 0; dx < c_dofs1D; ++dx)
7147 for (int dy = 0; dy < c_dofs1D; ++dy)
7149 for (int dz = 0; dz < c_dofs1D; ++dz)
7151 for (int j=0; j<3; ++j)
7154 for (int ez = 0; ez < c_dofs1D; ++ez)
7156 s += w1[j][dx][dy][ez] * Bc(ez, dz);
7158 y(dx, dy, dz, j, e) += s;
7164 // dofs that point parallel to y-axis (open in y, closed in x, z)
7167 for (int ez = 0; ez < c_dofs1D; ++ez)
7169 for (int ey = 0; ey < o_dofs1D; ++ey)
7171 for (int j=0; j<3; ++j)
7173 for (int dx = 0; dx < c_dofs1D; ++dx)
7175 w2[j][dx][ey][ez] = 0.0;
7177 for (int ex = 0; ex < c_dofs1D; ++ex)
7179 const int local_index = c_dofs1D*c_dofs1D*o_dofs1D +
7180 ez*c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
7181 const double xv = x(local_index, e) * vk(j, local_index, e);
7182 for (int dx = 0; dx < c_dofs1D; ++dx)
7184 w2[j][dx][ey][ez] += xv * Bc(ex, dx);
7192 for (int ez = 0; ez < c_dofs1D; ++ez)
7194 for (int dx = 0; dx < c_dofs1D; ++dx)
7196 for (int dy = 0; dy < c_dofs1D; ++dy)
7198 for (int j=0; j<3; ++j)
7200 w1[j][dx][dy][ez] = 0.0;
7201 for (int ey = 0; ey < o_dofs1D; ++ey)
7203 w1[j][dx][dy][ez] += w2[j][dx][ey][ez] * Bo(ey, dy);
7211 for (int dx = 0; dx < c_dofs1D; ++dx)
7213 for (int dy = 0; dy < c_dofs1D; ++dy)
7215 for (int dz = 0; dz < c_dofs1D; ++dz)
7217 for (int j=0; j<3; ++j)
7220 for (int ez = 0; ez < c_dofs1D; ++ez)
7222 s += w1[j][dx][dy][ez] * Bc(ez, dz);
7224 y(dx, dy, dz, j, e) += s;
7230 // dofs that point parallel to z-axis (open in z, closed in x, y)
7233 for (int ez = 0; ez < o_dofs1D; ++ez)
7235 for (int ey = 0; ey < c_dofs1D; ++ey)
7237 for (int j=0; j<3; ++j)
7239 for (int dx = 0; dx < c_dofs1D; ++dx)
7241 w2[j][dx][ey][ez] = 0.0;
7243 for (int ex = 0; ex < c_dofs1D; ++ex)
7245 const int local_index = 2*c_dofs1D*c_dofs1D*o_dofs1D +
7246 ez*c_dofs1D*c_dofs1D + ey*c_dofs1D + ex;
7247 const double xv = x(local_index, e) * vk(j, local_index, e);
7248 for (int dx = 0; dx < c_dofs1D; ++dx)
7250 w2[j][dx][ey][ez] += xv * Bc(ex, dx);
7258 for (int ez = 0; ez < o_dofs1D; ++ez)
7260 for (int dx = 0; dx < c_dofs1D; ++dx)
7262 for (int dy = 0; dy < c_dofs1D; ++dy)
7264 for (int j=0; j<3; ++j)
7266 w1[j][dx][dy][ez] = 0.0;
7267 for (int ey = 0; ey < c_dofs1D; ++ey)
7269 w1[j][dx][dy][ez] += w2[j][dx][ey][ez] * Bc(ey, dy);
7277 for (int dx = 0; dx < c_dofs1D; ++dx)
7279 for (int dy = 0; dy < c_dofs1D; ++dy)
7281 for (int dz = 0; dz < c_dofs1D; ++dz)
7283 for (int j=0; j<3; ++j)
7286 for (int ez = 0; ez < o_dofs1D; ++ez)
7288 s += w1[j][dx][dy][ez] * Bo(ez, dz);
7290 y(dx, dy, dz, j, e) += s;
7298 static void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
7301 const Array<double> &Bclosed,
7302 const Array<double> &Bopen,
7303 const Vector &pa_data,
7307 auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
7308 auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
7310 auto x = Reshape(x_.Read(), c_dofs1D, c_dofs1D, 2, NE);
7311 auto y = Reshape(y_.ReadWrite(), (2 * c_dofs1D * o_dofs1D), NE);
7313 auto vk = Reshape(pa_data.Read(), 2, (2 * c_dofs1D * o_dofs1D), NE);
7315 constexpr static int MAX_D1D = HCURL_MAX_D1D;
7317 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
7321 double w[2][MAX_D1D][MAX_D1D];
7323 // dofs that point parallel to x-axis (open in x, closed in y)
7326 for (int ey = 0; ey < c_dofs1D; ++ey)
7328 for (int dx = 0; dx < c_dofs1D; ++dx)
7330 for (int j=0; j<2; ++j)
7333 for (int dy = 0; dy < c_dofs1D; ++dy)
7335 w[j][dx][ey] += Bc(ey, dy) * x(dx, dy, j, e);
7342 for (int ey = 0; ey < c_dofs1D; ++ey)
7344 for (int ex = 0; ex < o_dofs1D; ++ex)
7346 for (int j=0; j<2; ++j)
7349 for (int dx = 0; dx < c_dofs1D; ++dx)
7351 s += Bo(ex, dx) * w[j][dx][ey];
7353 const int local_index = ey*o_dofs1D + ex;
7354 y(local_index, e) += s * vk(j, local_index, e);
7359 // dofs that point parallel to y-axis (open in y, closed in x)
7362 for (int ey = 0; ey < o_dofs1D; ++ey)
7364 for (int dx = 0; dx < c_dofs1D; ++dx)
7366 for (int j=0; j<2; ++j)
7369 for (int dy = 0; dy < c_dofs1D; ++dy)
7371 w[j][dx][ey] += Bo(ey, dy) * x(dx, dy, j, e);
7378 for (int ey = 0; ey < o_dofs1D; ++ey)
7380 for (int ex = 0; ex < c_dofs1D; ++ex)
7382 for (int j=0; j<2; ++j)
7385 for (int dx = 0; dx < c_dofs1D; ++dx)
7387 s += Bc(ex, dx) * w[j][dx][ey];
7389 const int local_index = c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
7390 y(local_index, e) += s * vk(j, local_index, e);
7397 static void PAHcurlVecH1IdentityApplyTranspose2D(const int c_dofs1D,
7400 const Array<double> &Bclosed,
7401 const Array<double> &Bopen,
7402 const Vector &pa_data,
7406 auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
7407 auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
7409 auto x = Reshape(x_.Read(), (2 * c_dofs1D * o_dofs1D), NE);
7410 auto y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, 2, NE);
7412 auto vk = Reshape(pa_data.Read(), 2, (2 * c_dofs1D * o_dofs1D), NE);
7414 constexpr static int MAX_D1D = HCURL_MAX_D1D;
7415 //constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
7417 MFEM_VERIFY(c_dofs1D <= MAX_D1D && o_dofs1D <= c_dofs1D, "");
7421 double w[2][MAX_D1D][MAX_D1D];
7423 // dofs that point parallel to x-axis (open in x, closed in y)
7426 for (int ey = 0; ey < c_dofs1D; ++ey)
7428 for (int dx = 0; dx < c_dofs1D; ++dx)
7430 for (int j=0; j<2; ++j) { w[j][dx][ey] = 0.0; }
7432 for (int ex = 0; ex < o_dofs1D; ++ex)
7434 const int local_index = ey*o_dofs1D + ex;
7435 const double xd = x(local_index, e);
7437 for (int dx = 0; dx < c_dofs1D; ++dx)
7439 for (int j=0; j<2; ++j)
7441 w[j][dx][ey] += Bo(ex, dx) * xd * vk(j, local_index, e);
7448 for (int dx = 0; dx < c_dofs1D; ++dx)
7450 for (int dy = 0; dy < c_dofs1D; ++dy)
7452 for (int j=0; j<2; ++j)
7455 for (int ey = 0; ey < c_dofs1D; ++ey)
7457 s += w[j][dx][ey] * Bc(ey, dy);
7459 y(dx, dy, j, e) += s;
7464 // dofs that point parallel to y-axis (open in y, closed in x)
7467 for (int ey = 0; ey < o_dofs1D; ++ey)
7469 for (int dx = 0; dx < c_dofs1D; ++dx)
7471 for (int j=0; j<2; ++j) { w[j][dx][ey] = 0.0; }
7473 for (int ex = 0; ex < c_dofs1D; ++ex)
7475 const int local_index = c_dofs1D*o_dofs1D + ey*c_dofs1D + ex;
7476 const double xd = x(local_index, e);
7477 for (int dx = 0; dx < c_dofs1D; ++dx)
7479 for (int j=0; j<2; ++j)
7481 w[j][dx][ey] += Bc(ex, dx) * xd * vk(j, local_index, e);
7488 for (int dx = 0; dx < c_dofs1D; ++dx)
7490 for (int dy = 0; dy < c_dofs1D; ++dy)
7492 for (int j=0; j<2; ++j)
7495 for (int ey = 0; ey < o_dofs1D; ++ey)
7497 s += w[j][dx][ey] * Bo(ey, dy);
7499 y(dx, dy, j, e) += s;
7506 void IdentityInterpolator::AssemblePA(const FiniteElementSpace &trial_fes,
7507 const FiniteElementSpace &test_fes)
7509 // Assumes tensor-product elements, with a vector test space and H^1 trial space.
7510 Mesh *mesh = trial_fes.GetMesh();
7511 const FiniteElement *trial_fel = trial_fes.GetFE(0);
7512 const FiniteElement *test_fel = test_fes.GetFE(0);
7514 const NodalTensorFiniteElement *trial_el =
7515 dynamic_cast<const NodalTensorFiniteElement*>(trial_fel);
7518 const VectorTensorFiniteElement *test_el =
7519 dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
7522 const int dims = trial_el->GetDim();
7523 MFEM_VERIFY(dims == 2 || dims == 3, "");
7525 dim = mesh->Dimension();
7526 MFEM_VERIFY(dim == 2 || dim == 3, "");
7528 MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(), "");
7530 ne = trial_fes.GetNE();
7532 const int order = trial_el->GetOrder();
7533 dofquad_fe = new H1_SegmentElement(order);
7534 mfem::QuadratureFunctions1D qf1d;
7535 mfem::IntegrationRule closed_ir;
7536 closed_ir.SetSize(order + 1);
7537 qf1d.GaussLobatto(order + 1, &closed_ir);
7538 mfem::IntegrationRule open_ir;
7539 open_ir.SetSize(order);
7540 qf1d.GaussLegendre(order, &open_ir);
7542 maps_C_C = &dofquad_fe->GetDofToQuad(closed_ir, DofToQuad::TENSOR);
7543 maps_O_C = &dofquad_fe->GetDofToQuad(open_ir, DofToQuad::TENSOR);
7545 o_dofs1D = maps_O_C->nqpt;
7546 c_dofs1D = maps_C_C->nqpt;
7547 MFEM_VERIFY(maps_O_C->ndof == c_dofs1D &&
7548 maps_C_C->ndof == c_dofs1D, "Discrepancy in the number of DOFs
");
7550 const int ndof_test = (dim == 3) ? 3 * c_dofs1D * c_dofs1D * o_dofs1D
7551 : 2 * c_dofs1D * o_dofs1D;
7553 const IntegrationRule & Nodes = test_el->GetNodes();
7555 pa_data.SetSize(dim * ndof_test * ne, Device::GetMemoryType());
7556 auto op = Reshape(pa_data.HostWrite(), dim, ndof_test, ne);
7558 const Array<int> &dofmap = test_el->GetDofMap();
7562 // Note that ND_HexahedronElement uses 6 vectors in tk rather than 3, with
7563 // the last 3 having negative signs. Here the signs are all positive, as
7564 // signs are applied in ElementRestriction.
7566 const double tk[9] = { 1.,0.,0., 0.,1.,0., 0.,0.,1. };
7568 for (int c=0; c<3; ++c)
7570 for (int i=0; i<ndof_test/3; ++i)
7572 const int d = (c*ndof_test/3) + i;
7573 // ND_HexahedronElement sets dof2tk = (dofmap < 0) ? 3+c : c, but here
7574 // no signs should be applied due to ElementRestriction.
7575 const int dof2tk = c;
7576 const int id = (dofmap[d] >= 0) ? dofmap[d] : -1 - dofmap[d];
7578 for (int e=0; e<ne; ++e)
7581 ElementTransformation *tr = mesh->GetElementTransformation(e);
7582 tr->SetIntPoint(&Nodes.IntPoint(id));
7583 tr->Jacobian().Mult(tk + dof2tk*dim, v);
7585 for (int j=0; j<3; ++j)
7595 const double tk[4] = { 1.,0., 0.,1. };
7596 for (int c=0; c<2; ++c)
7598 for (int i=0; i<ndof_test/2; ++i)
7600 const int d = (c*ndof_test/2) + i;
7601 // ND_QuadrilateralElement sets dof2tk = (dofmap < 0) ? 2+c : c, but here
7602 // no signs should be applied due to ElementRestriction.
7603 const int dof2tk = c;
7604 const int id = (dofmap[d] >= 0) ? dofmap[d] : -1 - dofmap[d];
7606 for (int e=0; e<ne; ++e)
7609 ElementTransformation *tr = mesh->GetElementTransformation(e);
7610 tr->SetIntPoint(&Nodes.IntPoint(id));
7611 tr->Jacobian().Mult(tk + dof2tk*dim, v);
7613 for (int j=0; j<2; ++j)
7623 void IdentityInterpolator::AddMultPA(const Vector &x, Vector &y) const
7627 PAHcurlVecH1IdentityApply3D(c_dofs1D, o_dofs1D, ne, maps_C_C->B, maps_O_C->B,
7632 PAHcurlVecH1IdentityApply2D(c_dofs1D, o_dofs1D, ne, maps_C_C->B, maps_O_C->B,
7641 void IdentityInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
7645 PAHcurlVecH1IdentityApplyTranspose3D(c_dofs1D, o_dofs1D, ne, maps_C_C->B,
7646 maps_O_C->B, pa_data, x, y);
7650 PAHcurlVecH1IdentityApplyTranspose2D(c_dofs1D, o_dofs1D, ne, maps_C_C->B,
7651 maps_O_C->B, pa_data, x, y);
7659 template void SmemPAHcurlMassAssembleDiagonal3D<0,0>(const int D1D,
7662 const bool symmetric,
7663 const Array<double> &bo,
7664 const Array<double> &bc,
7665 const Vector &pa_data,
7668 template void SmemPAHcurlMassAssembleDiagonal3D<2,3>(const int D1D,
7671 const bool symmetric,
7672 const Array<double> &bo,
7673 const Array<double> &bc,
7674 const Vector &pa_data,
7677 template void SmemPAHcurlMassAssembleDiagonal3D<3,4>(const int D1D,
7680 const bool symmetric,
7681 const Array<double> &bo,
7682 const Array<double> &bc,
7683 const Vector &pa_data,
7686 template void SmemPAHcurlMassAssembleDiagonal3D<4,5>(const int D1D,
7689 const bool symmetric,
7690 const Array<double> &bo,
7691 const Array<double> &bc,
7692 const Vector &pa_data,
7695 template void SmemPAHcurlMassAssembleDiagonal3D<5,6>(const int D1D,
7698 const bool symmetric,
7699 const Array<double> &bo,
7700 const Array<double> &bc,
7701 const Vector &pa_data,
7704 template void SmemPAHcurlMassApply3D<0,0>(const int D1D,
7707 const bool symmetric,
7708 const Array<double> &bo,
7709 const Array<double> &bc,
7710 const Array<double> &bot,
7711 const Array<double> &bct,
7712 const Vector &pa_data,
7716 template void SmemPAHcurlMassApply3D<2,3>(const int D1D,
7719 const bool symmetric,
7720 const Array<double> &bo,
7721 const Array<double> &bc,
7722 const Array<double> &bot,
7723 const Array<double> &bct,
7724 const Vector &pa_data,
7728 template void SmemPAHcurlMassApply3D<3,4>(const int D1D,
7731 const bool symmetric,
7732 const Array<double> &bo,
7733 const Array<double> &bc,
7734 const Array<double> &bot,
7735 const Array<double> &bct,
7736 const Vector &pa_data,
7740 template void SmemPAHcurlMassApply3D<4,5>(const int D1D,
7743 const bool symmetric,
7744 const Array<double> &bo,
7745 const Array<double> &bc,
7746 const Array<double> &bot,
7747 const Array<double> &bct,
7748 const Vector &pa_data,
7752 template void SmemPAHcurlMassApply3D<5,6>(const int D1D,
7755 const bool symmetric,
7756 const Array<double> &bo,
7757 const Array<double> &bc,
7758 const Array<double> &bot,
7759 const Array<double> &bct,
7760 const Vector &pa_data,
void PAHcurlMassAssembleDiagonal2D(const int D1D, const int Q1D, const int NE, const bool symmetric, const Array< double > &bo, const Array< double > &bc, const Vector &pa_data, Vector &diag)
constexpr int HCURL_MAX_D1D
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
void PAHcurlMassAssembleDiagonal3D(const int D1D, const int Q1D, const int NE, const bool symmetric, const Array< double > &bo, const Array< double > &bc, const Vector &pa_data, Vector &diag)
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
void PAHcurlMassApply2D(const int D1D, const int Q1D, const int NE, const bool symmetric, const Array< double > &bo, const Array< double > &bc, const Array< double > &bot, const Array< double > &bct, const Vector &pa_data, const Vector &x, Vector &y)
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
void PAHcurlHdivSetup3D(const int Q1D, const int coeffDim, const int NE, const bool transpose, const Array< double > &w_, const Vector &j, Vector &coeff_, Vector &op)
constexpr int HCURL_MAX_Q1D
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.