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, const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
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
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)
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
constexpr int HCURL_MAX_Q1D