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