MFEM  v4.5.2
Finite element discretization library
bilininteg_vectorfe.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 "qspace.hpp"
15 #include "gridfunc.hpp"
16 
17 namespace mfem
18 {
19 
20 void PADiffusionSetup3D(const int Q1D,
21  const int coeffDim,
22  const int NE,
23  const Array<double> &w,
24  const Vector &j,
25  const Vector &coeff_,
26  Vector &op);
27 
28 void PAHcurlMassAssembleDiagonal2D(const int D1D,
29  const int Q1D,
30  const int NE,
31  const bool symmetric,
32  const Array<double> &bo,
33  const Array<double> &bc,
34  const Vector &pa_data,
35  Vector &diag);
36 
37 void PAHcurlMassAssembleDiagonal3D(const int D1D,
38  const int Q1D,
39  const int NE,
40  const bool symmetric,
41  const Array<double> &bo,
42  const Array<double> &bc,
43  const Vector &pa_data,
44  Vector &diag);
45 
46 template<int T_D1D = 0, int T_Q1D = 0>
47 void SmemPAHcurlMassAssembleDiagonal3D(const int D1D,
48  const int Q1D,
49  const int NE,
50  const bool symmetric,
51  const Array<double> &bo,
52  const Array<double> &bc,
53  const Vector &pa_data,
54  Vector &diag);
55 
56 void PAHcurlMassApply2D(const int D1D,
57  const int Q1D,
58  const int NE,
59  const bool symmetric,
60  const Array<double> &bo,
61  const Array<double> &bc,
62  const Array<double> &bot,
63  const Array<double> &bct,
64  const Vector &pa_data,
65  const Vector &x,
66  Vector &y);
67 
68 void PAHcurlMassApply3D(const int D1D,
69  const int Q1D,
70  const int NE,
71  const bool symmetric,
72  const Array<double> &bo,
73  const Array<double> &bc,
74  const Array<double> &bot,
75  const Array<double> &bct,
76  const Vector &pa_data,
77  const Vector &x,
78  Vector &y);
79 
80 template<int T_D1D = 0, int T_Q1D = 0>
81 void SmemPAHcurlMassApply3D(const int D1D,
82  const int Q1D,
83  const int NE,
84  const bool symmetric,
85  const Array<double> &bo,
86  const Array<double> &bc,
87  const Array<double> &bot,
88  const Array<double> &bct,
89  const Vector &pa_data,
90  const Vector &x,
91  Vector &y);
92 
93 void PAHdivSetup2D(const int Q1D,
94  const int coeffDim,
95  const int NE,
96  const Array<double> &w,
97  const Vector &j,
98  Vector &coeff_,
99  Vector &op);
100 
101 void PAHdivSetup3D(const int Q1D,
102  const int coeffDim,
103  const int NE,
104  const Array<double> &w,
105  const Vector &j,
106  Vector &coeff_,
107  Vector &op);
108 
109 void PAHcurlH1Apply2D(const int D1D,
110  const int Q1D,
111  const int NE,
112  const Array<double> &bc,
113  const Array<double> &gc,
114  const Array<double> &bot,
115  const Array<double> &bct,
116  const Vector &pa_data,
117  const Vector &x,
118  Vector &y);
119 
120 void PAHcurlH1ApplyTranspose2D(const int D1D,
121  const int Q1D,
122  const int NE,
123  const Array<double> &bc,
124  const Array<double> &bo,
125  const Array<double> &bct,
126  const Array<double> &gct,
127  const Vector &pa_data,
128  const Vector &x,
129  Vector &y);
130 
131 void PAHcurlH1Apply3D(const int D1D,
132  const int Q1D,
133  const int NE,
134  const Array<double> &bc,
135  const Array<double> &gc,
136  const Array<double> &bot,
137  const Array<double> &bct,
138  const Vector &pa_data,
139  const Vector &x,
140  Vector &y);
141 
142 void PAHcurlH1ApplyTranspose3D(const int D1D,
143  const int Q1D,
144  const int NE,
145  const Array<double> &bc,
146  const Array<double> &bo,
147  const Array<double> &bct,
148  const Array<double> &gct,
149  const Vector &pa_data,
150  const Vector &x,
151  Vector &y);
152 
153 void PAHdivMassAssembleDiagonal2D(const int D1D,
154  const int Q1D,
155  const int NE,
156  const bool symmetric,
157  const Array<double> &Bo_,
158  const Array<double> &Bc_,
159  const Vector &op_,
160  Vector &diag_);
161 
162 void PAHdivMassAssembleDiagonal3D(const int D1D,
163  const int Q1D,
164  const int NE,
165  const bool symmetric,
166  const Array<double> &Bo_,
167  const Array<double> &Bc_,
168  const Vector &op_,
169  Vector &diag_);
170 
171 void PAHdivMassApply(const int dim,
172  const int D1D,
173  const int Q1D,
174  const int NE,
175  const bool symmetric,
176  const Array<double> &Bo,
177  const Array<double> &Bc,
178  const Array<double> &Bot,
179  const Array<double> &Bct,
180  const Vector &op,
181  const Vector &x,
182  Vector &y);
183 
184 void PAHcurlL2Setup(const int NQ,
185  const int coeffDim,
186  const int NE,
187  const Array<double> &w,
188  Vector &coeff_,
189  Vector &op);
190 
191 // PA H(curl) x H(div) mass assemble 3D kernel, with factor
192 // dF^{-1} C dF for a vector or matrix coefficient C.
193 // If transpose, use dF^T C dF^{-T} for H(div) x H(curl).
194 void PAHcurlHdivSetup3D(const int Q1D,
195  const int coeffDim,
196  const int NE,
197  const bool transpose,
198  const Array<double> &w_,
199  const Vector &j,
200  Vector &coeff_,
201  Vector &op)
202 {
203  const bool symmetric = (coeffDim != 9);
204  auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
205  auto J = Reshape(j.Read(), Q1D, Q1D, Q1D, 3, 3, NE);
206  auto coeff = Reshape(coeff_.Read(), coeffDim, Q1D, Q1D, Q1D, NE);
207  auto y = Reshape(op.Write(), 9, Q1D, Q1D, Q1D, NE);
208 
209  const int i11 = 0;
210  const int i12 = transpose ? 3 : 1;
211  const int i13 = transpose ? 6 : 2;
212  const int i21 = transpose ? 1 : 3;
213  const int i22 = 4;
214  const int i23 = transpose ? 7 : 5;
215  const int i31 = transpose ? 2 : 6;
216  const int i32 = transpose ? 5 : 7;
217  const int i33 = 8;
218 
219  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
220  {
221  MFEM_FOREACH_THREAD(qx,x,Q1D)
222  {
223  MFEM_FOREACH_THREAD(qy,y,Q1D)
224  {
225  MFEM_FOREACH_THREAD(qz,z,Q1D)
226  {
227  const double J11 = J(qx,qy,qz,0,0,e);
228  const double J21 = J(qx,qy,qz,1,0,e);
229  const double J31 = J(qx,qy,qz,2,0,e);
230  const double J12 = J(qx,qy,qz,0,1,e);
231  const double J22 = J(qx,qy,qz,1,1,e);
232  const double J32 = J(qx,qy,qz,2,1,e);
233  const double J13 = J(qx,qy,qz,0,2,e);
234  const double J23 = J(qx,qy,qz,1,2,e);
235  const double J33 = J(qx,qy,qz,2,2,e);
236  const double detJ = J11 * (J22 * J33 - J32 * J23) -
237  /* */ J21 * (J12 * J33 - J32 * J13) +
238  /* */ J31 * (J12 * J23 - J22 * J13);
239  const double w_detJ = W(qx,qy,qz) / detJ;
240  // adj(J)
241  const double A11 = (J22 * J33) - (J23 * J32);
242  const double A12 = (J32 * J13) - (J12 * J33);
243  const double A13 = (J12 * J23) - (J22 * J13);
244  const double A21 = (J31 * J23) - (J21 * J33);
245  const double A22 = (J11 * J33) - (J13 * J31);
246  const double A23 = (J21 * J13) - (J11 * J23);
247  const double A31 = (J21 * J32) - (J31 * J22);
248  const double A32 = (J31 * J12) - (J11 * J32);
249  const double A33 = (J11 * J22) - (J12 * J21);
250 
251  if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
252  {
253  // First compute entries of R = M^T J
254  const double M11 = (!symmetric) ? coeff(i11,qx,qy,qz,e) : coeff(0,qx,qy,qz,e);
255  const double M12 = (!symmetric) ? coeff(i12,qx,qy,qz,e) : coeff(1,qx,qy,qz,e);
256  const double M13 = (!symmetric) ? coeff(i13,qx,qy,qz,e) : coeff(2,qx,qy,qz,e);
257  const double M21 = (!symmetric) ? coeff(i21,qx,qy,qz,e) : M12;
258  const double M22 = (!symmetric) ? coeff(i22,qx,qy,qz,e) : coeff(3,qx,qy,qz,e);
259  const double M23 = (!symmetric) ? coeff(i23,qx,qy,qz,e) : coeff(4,qx,qy,qz,e);
260  const double M31 = (!symmetric) ? coeff(i31,qx,qy,qz,e) : M13;
261  const double M32 = (!symmetric) ? coeff(i32,qx,qy,qz,e) : M23;
262  const double M33 = (!symmetric) ? coeff(i33,qx,qy,qz,e) : coeff(5,qx,qy,qz,e);
263 
264  const double R11 = M11*J11 + M21*J21 + M31*J31;
265  const double R12 = M11*J12 + M21*J22 + M31*J32;
266  const double R13 = M11*J13 + M21*J23 + M31*J33;
267  const double R21 = M12*J11 + M22*J21 + M32*J31;
268  const double R22 = M12*J12 + M22*J22 + M32*J32;
269  const double R23 = M12*J13 + M22*J23 + M32*J33;
270  const double R31 = M13*J11 + M23*J21 + M33*J31;
271  const double R32 = M13*J12 + M23*J22 + M33*J32;
272  const double R33 = M13*J13 + M23*J23 + M33*J33;
273 
274  // y = (J^{-1} M^T J)^T
275  y(i11,qx,qy,qz,e) = w_detJ * (A11*R11 + A12*R21 + A13*R31); // 1,1
276  y(i21,qx,qy,qz,e) = w_detJ * (A11*R12 + A12*R22 + A13*R32); // 1,2
277  y(i31,qx,qy,qz,e) = w_detJ * (A11*R13 + A12*R23 + A13*R33); // 1,3
278  y(i12,qx,qy,qz,e) = w_detJ * (A21*R11 + A22*R21 + A23*R31); // 2,1
279  y(i22,qx,qy,qz,e) = w_detJ * (A21*R12 + A22*R22 + A23*R32); // 2,2
280  y(i32,qx,qy,qz,e) = w_detJ * (A21*R13 + A22*R23 + A23*R33); // 2,3
281  y(i13,qx,qy,qz,e) = w_detJ * (A31*R11 + A32*R21 + A33*R31); // 3,1
282  y(i23,qx,qy,qz,e) = w_detJ * (A31*R12 + A32*R22 + A33*R32); // 3,2
283  y(i33,qx,qy,qz,e) = w_detJ * (A31*R13 + A32*R23 + A33*R33); // 3,3
284  }
285  else if (coeffDim == 3) // Vector coefficient version
286  {
287  const double D1 = coeff(0,qx,qy,qz,e);
288  const double D2 = coeff(1,qx,qy,qz,e);
289  const double D3 = coeff(2,qx,qy,qz,e);
290  // detJ J^{-1} DJ = adj(J) DJ
291  // transpose
292  y(i11,qx,qy,qz,e) = w_detJ * (D1*A11*J11 + D2*A12*J21 + D3*A13*J31); // 1,1
293  y(i21,qx,qy,qz,e) = w_detJ * (D1*A11*J12 + D2*A12*J22 + D3*A13*J32); // 1,2
294  y(i31,qx,qy,qz,e) = w_detJ * (D1*A11*J13 + D2*A12*J23 + D3*A13*J33); // 1,3
295  y(i12,qx,qy,qz,e) = w_detJ * (D1*A21*J11 + D2*A22*J21 + D3*A23*J31); // 2,1
296  y(i22,qx,qy,qz,e) = w_detJ * (D1*A21*J12 + D2*A22*J22 + D3*A23*J32); // 2,2
297  y(i32,qx,qy,qz,e) = w_detJ * (D1*A21*J13 + D2*A22*J23 + D3*A23*J33); // 2,3
298  y(i13,qx,qy,qz,e) = w_detJ * (D1*A31*J11 + D2*A32*J21 + D3*A33*J31); // 3,1
299  y(i23,qx,qy,qz,e) = w_detJ * (D1*A31*J12 + D2*A32*J22 + D3*A33*J32); // 3,2
300  y(i33,qx,qy,qz,e) = w_detJ * (D1*A31*J13 + D2*A32*J23 + D3*A33*J33); // 3,3
301  }
302  }
303  }
304  }
305  });
306 }
307 
308 // PA H(curl) x H(div) mass assemble 2D kernel, with factor
309 // dF^{-1} C dF for a vector or matrix coefficient C.
310 // If transpose, use dF^T C dF^{-T} for H(div) x H(curl).
311 void PAHcurlHdivSetup2D(const int Q1D,
312  const int coeffDim,
313  const int NE,
314  const bool transpose,
315  const Array<double> &w_,
316  const Vector &j,
317  Vector &coeff_,
318  Vector &op)
319 {
320  const bool symmetric = (coeffDim != 4);
321  auto W = Reshape(w_.Read(), Q1D, Q1D);
322  auto J = Reshape(j.Read(), Q1D, Q1D, 2, 2, NE);
323  auto coeff = Reshape(coeff_.Read(), coeffDim, Q1D, Q1D, NE);
324  auto y = Reshape(op.Write(), 4, Q1D, Q1D, NE);
325 
326  const int i11 = 0;
327  const int i12 = transpose ? 2 : 1;
328  const int i21 = transpose ? 1 : 2;
329  const int i22 = 3;
330 
331  MFEM_FORALL_2D(e, NE, Q1D, Q1D, 1,
332  {
333  MFEM_FOREACH_THREAD(qx,x,Q1D)
334  {
335  MFEM_FOREACH_THREAD(qy,y,Q1D)
336  {
337  const double J11 = J(qx,qy,0,0,e);
338  const double J21 = J(qx,qy,1,0,e);
339  const double J12 = J(qx,qy,0,1,e);
340  const double J22 = J(qx,qy,1,1,e);
341  const double w_detJ = W(qx,qy) / ((J11*J22) - (J21*J12));
342 
343  if (coeffDim == 3 || coeffDim == 4) // Matrix coefficient version
344  {
345  // First compute entries of R = MJ
346  const double M11 = coeff(i11,qx,qy,e);
347  const double M12 = (!symmetric) ? coeff(i12,qx,qy,e) : coeff(1,qx,qy,e);
348  const double M21 = (!symmetric) ? coeff(i21,qx,qy,e) : M12;
349  const double M22 = (!symmetric) ? coeff(i22,qx,qy,e) : coeff(2,qx,qy,e);
350 
351  // J^{-1} M^T
352  const double R11 = ( J22*M11 - J12*M12); // 1,1
353  const double R12 = ( J22*M21 - J12*M22); // 1,2
354  const double R21 = (-J21*M11 + J11*M12); // 2,1
355  const double R22 = (-J21*M21 + J11*M22); // 2,2
356 
357  // (RJ)^T
358  y(i11,qx,qy,e) = w_detJ * (R11*J11 + R12*J21); // 1,1
359  y(i21,qx,qy,e) = w_detJ * (R11*J12 + R12*J22); // 1,2 (transpose)
360  y(i12,qx,qy,e) = w_detJ * (R21*J11 + R22*J21); // 2,1 (transpose)
361  y(i22,qx,qy,e) = w_detJ * (R21*J12 + R22*J22); // 2,2
362  }
363  else if (coeffDim == 2) // Vector coefficient version
364  {
365  const double D1 = coeff(0,qx,qy,e);
366  const double D2 = coeff(1,qx,qy,e);
367  const double R11 = D1*J11;
368  const double R12 = D1*J12;
369  const double R21 = D2*J21;
370  const double R22 = D2*J22;
371  y(i11,qx,qy,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
372  y(i21,qx,qy,e) = w_detJ * ( J22*R12 - J12*R22); // 1,2 (transpose)
373  y(i12,qx,qy,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1 (transpose)
374  y(i22,qx,qy,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
375  }
376  }
377  }
378  });
379 }
380 
381 // Mass operator for H(curl) and H(div) functions, using Piola transformations
382 // u = dF^{-T} \hat{u} in H(curl), v = (1 / det dF) dF \hat{v} in H(div).
383 void PAHcurlHdivMassApply3D(const int D1D,
384  const int D1Dtest,
385  const int Q1D,
386  const int NE,
387  const bool scalarCoeff,
388  const bool trialHcurl,
389  const bool transpose,
390  const Array<double> &Bo_,
391  const Array<double> &Bc_,
392  const Array<double> &Bot_,
393  const Array<double> &Bct_,
394  const Vector &op_,
395  const Vector &x_,
396  Vector &y_)
397 {
398  constexpr static int MAX_D1D = HCURL_MAX_D1D;
399  constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
400 
401  MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
402  MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
403  constexpr static int VDIM = 3;
404 
405  auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
406  auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
407  auto Bot = Reshape(Bot_.Read(), D1Dtest-1, Q1D);
408  auto Bct = Reshape(Bct_.Read(), D1Dtest, Q1D);
409  auto op = Reshape(op_.Read(), scalarCoeff ? 1 : 9, Q1D, Q1D, Q1D, NE);
410  auto x = Reshape(x_.Read(), 3*(D1D-1)*D1D*(trialHcurl ? D1D : D1D-1), NE);
411  auto y = Reshape(y_.ReadWrite(), 3*(D1Dtest-1)*D1Dtest*
412  (trialHcurl ? D1Dtest-1 : D1Dtest), NE);
413 
414  const int i12 = transpose ? 3 : 1;
415  const int i13 = transpose ? 6 : 2;
416  const int i21 = transpose ? 1 : 3;
417  const int i23 = transpose ? 7 : 5;
418  const int i31 = transpose ? 2 : 6;
419  const int i32 = transpose ? 5 : 7;
420 
421  MFEM_FORALL(e, NE,
422  {
423  double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
424 
425  for (int qz = 0; qz < Q1D; ++qz)
426  {
427  for (int qy = 0; qy < Q1D; ++qy)
428  {
429  for (int qx = 0; qx < Q1D; ++qx)
430  {
431  for (int c = 0; c < VDIM; ++c)
432  {
433  mass[qz][qy][qx][c] = 0.0;
434  }
435  }
436  }
437  }
438 
439  int osc = 0;
440  for (int c = 0; c < VDIM; ++c) // loop over x, y, z trial components
441  {
442  const int D1Dz = trialHcurl ? ((c == 2) ? D1D - 1 : D1D) :
443  ((c == 2) ? D1D : D1D - 1);
444  const int D1Dy = trialHcurl ? ((c == 1) ? D1D - 1 : D1D) :
445  ((c == 1) ? D1D : D1D - 1);
446  const int D1Dx = trialHcurl ? ((c == 0) ? D1D - 1 : D1D) :
447  ((c == 0) ? D1D : D1D - 1);
448 
449  for (int dz = 0; dz < D1Dz; ++dz)
450  {
451  double massXY[MAX_Q1D][MAX_Q1D];
452  for (int qy = 0; qy < Q1D; ++qy)
453  {
454  for (int qx = 0; qx < Q1D; ++qx)
455  {
456  massXY[qy][qx] = 0.0;
457  }
458  }
459 
460  for (int dy = 0; dy < D1Dy; ++dy)
461  {
462  double massX[MAX_Q1D];
463  for (int qx = 0; qx < Q1D; ++qx)
464  {
465  massX[qx] = 0.0;
466  }
467 
468  for (int dx = 0; dx < D1Dx; ++dx)
469  {
470  const double t = x(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
471  for (int qx = 0; qx < Q1D; ++qx)
472  {
473  massX[qx] += t * (trialHcurl ? ((c == 0) ? Bo(qx,dx) : Bc(qx,dx)) :
474  ((c == 0) ? Bc(qx,dx) : Bo(qx,dx)));
475  }
476  }
477 
478  for (int qy = 0; qy < Q1D; ++qy)
479  {
480  const double wy = trialHcurl ? ((c == 1) ? Bo(qy,dy) : Bc(qy,dy)) :
481  ((c == 1) ? Bc(qy,dy) : Bo(qy,dy));
482  for (int qx = 0; qx < Q1D; ++qx)
483  {
484  const double wx = massX[qx];
485  massXY[qy][qx] += wx * wy;
486  }
487  }
488  }
489 
490  for (int qz = 0; qz < Q1D; ++qz)
491  {
492  const double wz = trialHcurl ? ((c == 2) ? Bo(qz,dz) : Bc(qz,dz)) :
493  ((c == 2) ? Bc(qz,dz) : Bo(qz,dz));
494  for (int qy = 0; qy < Q1D; ++qy)
495  {
496  for (int qx = 0; qx < Q1D; ++qx)
497  {
498  mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
499  }
500  }
501  }
502  }
503 
504  osc += D1Dx * D1Dy * D1Dz;
505  } // loop (c) over components
506 
507  // Apply D operator.
508  for (int qz = 0; qz < Q1D; ++qz)
509  {
510  for (int qy = 0; qy < Q1D; ++qy)
511  {
512  for (int qx = 0; qx < Q1D; ++qx)
513  {
514  const double O11 = op(0,qx,qy,qz,e);
515  const double O12 = scalarCoeff ? 0.0 : op(i12,qx,qy,qz,e);
516  const double O13 = scalarCoeff ? 0.0 : op(i13,qx,qy,qz,e);
517  const double O21 = scalarCoeff ? 0.0 : op(i21,qx,qy,qz,e);
518  const double O22 = scalarCoeff ? O11 : op(4,qx,qy,qz,e);
519  const double O23 = scalarCoeff ? 0.0 : op(i23,qx,qy,qz,e);
520  const double O31 = scalarCoeff ? 0.0 : op(i31,qx,qy,qz,e);
521  const double O32 = scalarCoeff ? 0.0 : op(i32,qx,qy,qz,e);
522  const double O33 = scalarCoeff ? O11 : op(8,qx,qy,qz,e);
523  const double massX = mass[qz][qy][qx][0];
524  const double massY = mass[qz][qy][qx][1];
525  const double massZ = mass[qz][qy][qx][2];
526  mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
527  mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
528  mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
529  }
530  }
531  }
532 
533  for (int qz = 0; qz < Q1D; ++qz)
534  {
535  double massXY[HDIV_MAX_D1D][HDIV_MAX_D1D];
536 
537  osc = 0;
538  for (int c = 0; c < VDIM; ++c) // loop over x, y, z test components
539  {
540  const int D1Dz = trialHcurl ? ((c == 2) ? D1Dtest : D1Dtest - 1) :
541  ((c == 2) ? D1Dtest - 1 : D1Dtest);
542  const int D1Dy = trialHcurl ? ((c == 1) ? D1Dtest : D1Dtest - 1) :
543  ((c == 1) ? D1Dtest - 1 : D1Dtest);
544  const int D1Dx = trialHcurl ? ((c == 0) ? D1Dtest : D1Dtest - 1) :
545  ((c == 0) ? D1Dtest - 1 : D1Dtest);
546 
547  for (int dy = 0; dy < D1Dy; ++dy)
548  {
549  for (int dx = 0; dx < D1Dx; ++dx)
550  {
551  massXY[dy][dx] = 0.0;
552  }
553  }
554  for (int qy = 0; qy < Q1D; ++qy)
555  {
556  double massX[HDIV_MAX_D1D];
557  for (int dx = 0; dx < D1Dx; ++dx)
558  {
559  massX[dx] = 0.0;
560  }
561  for (int qx = 0; qx < Q1D; ++qx)
562  {
563  for (int dx = 0; dx < D1Dx; ++dx)
564  {
565  massX[dx] += mass[qz][qy][qx][c] * (trialHcurl ?
566  ((c == 0) ? Bct(dx,qx) : Bot(dx,qx)) :
567  ((c == 0) ? Bot(dx,qx) : Bct(dx,qx)));
568  }
569  }
570  for (int dy = 0; dy < D1Dy; ++dy)
571  {
572  const double wy = trialHcurl ? ((c == 1) ? Bct(dy,qy) : Bot(dy,qy)) :
573  ((c == 1) ? Bot(dy,qy) : Bct(dy,qy));
574  for (int dx = 0; dx < D1Dx; ++dx)
575  {
576  massXY[dy][dx] += massX[dx] * wy;
577  }
578  }
579  }
580 
581  for (int dz = 0; dz < D1Dz; ++dz)
582  {
583  const double wz = trialHcurl ? ((c == 2) ? Bct(dz,qz) : Bot(dz,qz)) :
584  ((c == 2) ? Bot(dz,qz) : Bct(dz,qz));
585  for (int dy = 0; dy < D1Dy; ++dy)
586  {
587  for (int dx = 0; dx < D1Dx; ++dx)
588  {
589  y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) +=
590  massXY[dy][dx] * wz;
591  }
592  }
593  }
594 
595  osc += D1Dx * D1Dy * D1Dz;
596  } // loop c
597  } // loop qz
598  }); // end of element loop
599 }
600 
601 // Mass operator for H(curl) and H(div) functions, using Piola transformations
602 // u = dF^{-T} \hat{u} in H(curl), v = (1 / det dF) dF \hat{v} in H(div).
603 void PAHcurlHdivMassApply2D(const int D1D,
604  const int D1Dtest,
605  const int Q1D,
606  const int NE,
607  const bool scalarCoeff,
608  const bool trialHcurl,
609  const bool transpose,
610  const Array<double> &Bo_,
611  const Array<double> &Bc_,
612  const Array<double> &Bot_,
613  const Array<double> &Bct_,
614  const Vector &op_,
615  const Vector &x_,
616  Vector &y_)
617 {
618  constexpr static int MAX_D1D = HCURL_MAX_D1D;
619  constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
620 
621  MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
622  MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
623  constexpr static int VDIM = 2;
624 
625  auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
626  auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
627  auto Bot = Reshape(Bot_.Read(), D1Dtest-1, Q1D);
628  auto Bct = Reshape(Bct_.Read(), D1Dtest, Q1D);
629  auto op = Reshape(op_.Read(), scalarCoeff ? 1 : 4, Q1D, Q1D, NE);
630  auto x = Reshape(x_.Read(), 2*(D1D-1)*D1D, NE);
631  auto y = Reshape(y_.ReadWrite(), 2*(D1Dtest-1)*D1Dtest, NE);
632 
633  const int i12 = transpose ? 2 : 1;
634  const int i21 = transpose ? 1 : 2;
635 
636  MFEM_FORALL(e, NE,
637  {
638  double mass[MAX_Q1D][MAX_Q1D][VDIM];
639 
640  for (int qy = 0; qy < Q1D; ++qy)
641  {
642  for (int qx = 0; qx < Q1D; ++qx)
643  {
644  for (int c = 0; c < VDIM; ++c)
645  {
646  mass[qy][qx][c] = 0.0;
647  }
648  }
649  }
650 
651  int osc = 0;
652  for (int c = 0; c < VDIM; ++c) // loop over x, y trial components
653  {
654  const int D1Dy = trialHcurl ? ((c == 1) ? D1D - 1 : D1D) :
655  ((c == 1) ? D1D : D1D - 1);
656  const int D1Dx = trialHcurl ? ((c == 0) ? D1D - 1 : D1D) :
657  ((c == 0) ? D1D : D1D - 1);
658 
659  for (int dy = 0; dy < D1Dy; ++dy)
660  {
661  double massX[MAX_Q1D];
662  for (int qx = 0; qx < Q1D; ++qx)
663  {
664  massX[qx] = 0.0;
665  }
666 
667  for (int dx = 0; dx < D1Dx; ++dx)
668  {
669  const double t = x(dx + (dy * D1Dx) + osc, e);
670  for (int qx = 0; qx < Q1D; ++qx)
671  {
672  massX[qx] += t * (trialHcurl ? ((c == 0) ? Bo(qx,dx) : Bc(qx,dx)) :
673  ((c == 0) ? Bc(qx,dx) : Bo(qx,dx)));
674  }
675  }
676 
677  for (int qy = 0; qy < Q1D; ++qy)
678  {
679  const double wy = trialHcurl ? ((c == 1) ? Bo(qy,dy) : Bc(qy,dy)) :
680  ((c == 1) ? Bc(qy,dy) : Bo(qy,dy));
681  for (int qx = 0; qx < Q1D; ++qx)
682  {
683  mass[qy][qx][c] += massX[qx] * wy;
684  }
685  }
686  }
687 
688  osc += D1Dx * D1Dy;
689  } // loop (c) over components
690 
691  // Apply D operator.
692  for (int qy = 0; qy < Q1D; ++qy)
693  {
694  for (int qx = 0; qx < Q1D; ++qx)
695  {
696  const double O11 = op(0,qx,qy,e);
697  const double O12 = scalarCoeff ? 0.0 : op(i12,qx,qy,e);
698  const double O21 = scalarCoeff ? 0.0 : op(i21,qx,qy,e);
699  const double O22 = scalarCoeff ? O11 : op(3,qx,qy,e);
700  const double massX = mass[qy][qx][0];
701  const double massY = mass[qy][qx][1];
702  mass[qy][qx][0] = (O11*massX)+(O12*massY);
703  mass[qy][qx][1] = (O21*massX)+(O22*massY);
704  }
705  }
706 
707  osc = 0;
708  for (int c = 0; c < VDIM; ++c) // loop over x, y test components
709  {
710  const int D1Dy = trialHcurl ? ((c == 1) ? D1Dtest : D1Dtest - 1) :
711  ((c == 1) ? D1Dtest - 1 : D1Dtest);
712  const int D1Dx = trialHcurl ? ((c == 0) ? D1Dtest : D1Dtest - 1) :
713  ((c == 0) ? D1Dtest - 1 : D1Dtest);
714 
715  for (int qy = 0; qy < Q1D; ++qy)
716  {
717  double massX[HDIV_MAX_D1D];
718  for (int dx = 0; dx < D1Dx; ++dx)
719  {
720  massX[dx] = 0.0;
721  }
722  for (int qx = 0; qx < Q1D; ++qx)
723  {
724  for (int dx = 0; dx < D1Dx; ++dx)
725  {
726  massX[dx] += mass[qy][qx][c] * (trialHcurl ?
727  ((c == 0) ? Bct(dx,qx) : Bot(dx,qx)) :
728  ((c == 0) ? Bot(dx,qx) : Bct(dx,qx)));
729  }
730  }
731  for (int dy = 0; dy < D1Dy; ++dy)
732  {
733  const double wy = trialHcurl ? ((c == 1) ? Bct(dy,qy) : Bot(dy,qy)) :
734  ((c == 1) ? Bot(dy,qy) : Bct(dy,qy));
735  for (int dx = 0; dx < D1Dx; ++dx)
736  {
737  y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
738  }
739  }
740  }
741 
742  osc += D1Dx * D1Dy;
743  } // loop c
744  }); // end of element loop
745 }
746 
747 void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
748 {
749  AssemblePA(fes, fes);
750 }
751 
752 void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
753  const FiniteElementSpace &test_fes)
754 {
755  // Assumes tensor-product elements
756  Mesh *mesh = trial_fes.GetMesh();
757 
758  const FiniteElement *trial_fel = trial_fes.GetFE(0);
759  const VectorTensorFiniteElement *trial_el =
760  dynamic_cast<const VectorTensorFiniteElement*>(trial_fel);
761  MFEM_VERIFY(trial_el != NULL, "Only VectorTensorFiniteElement is supported!");
762 
763  const FiniteElement *test_fel = test_fes.GetFE(0);
764  const VectorTensorFiniteElement *test_el =
765  dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
766  MFEM_VERIFY(test_el != NULL, "Only VectorTensorFiniteElement is supported!");
767 
768  const IntegrationRule *ir
769  = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
770  *mesh->GetElementTransformation(0));
771  const int dims = trial_el->GetDim();
772  MFEM_VERIFY(dims == 2 || dims == 3, "");
773 
774  const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
775  nq = ir->GetNPoints();
776  dim = mesh->Dimension();
777  MFEM_VERIFY(dim == 2 || dim == 3, "");
778 
779  ne = trial_fes.GetNE();
780  MFEM_VERIFY(ne == test_fes.GetNE(),
781  "Different meshes for test and trial spaces");
782  geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
783  mapsC = &trial_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
784  mapsO = &trial_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
785  dofs1D = mapsC->ndof;
786  quad1D = mapsC->nqpt;
787 
788  mapsCtest = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
789  mapsOtest = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
790  dofs1Dtest = mapsCtest->ndof;
791 
792  MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
793 
794  trial_fetype = trial_el->GetDerivType();
795  test_fetype = test_el->GetDerivType();
796 
797  const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
798  const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
799  const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
800  const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
801 
802  QuadratureSpace qs(*mesh, *ir);
803  CoefficientVector coeff(qs, CoefficientStorage::SYMMETRIC);
804  if (Q) { coeff.Project(*Q); }
805  else if (MQ) { coeff.ProjectTranspose(*MQ); }
806  else if (DQ) { coeff.Project(*DQ); }
807  else { coeff.SetConstant(1.0); }
808 
809  const int coeff_dim = coeff.GetVDim();
810  symmetric = (coeff_dim != dim*dim);
811 
812  if ((trial_curl && test_div) || (trial_div && test_curl))
813  pa_data.SetSize((coeff_dim == 1 ? 1 : dim*dim) * nq * ne,
814  Device::GetMemoryType());
815  else
816  pa_data.SetSize((symmetric ? symmDims : dims*dims) * nq * ne,
817  Device::GetMemoryType());
818 
819  if (trial_curl && test_curl && dim == 3)
820  {
821  PADiffusionSetup3D(quad1D, coeff_dim, ne, ir->GetWeights(), geom->J,
822  coeff, pa_data);
823  }
824  else if (trial_curl && test_curl && dim == 2)
825  {
826  PADiffusionSetup2D<2>(quad1D, coeff_dim, ne, ir->GetWeights(), geom->J,
827  coeff, pa_data);
828  }
829  else if (trial_div && test_div && dim == 3)
830  {
831  PAHdivSetup3D(quad1D, coeff_dim, ne, ir->GetWeights(), geom->J,
832  coeff, pa_data);
833  }
834  else if (trial_div && test_div && dim == 2)
835  {
836  PAHdivSetup2D(quad1D, coeff_dim, ne, ir->GetWeights(), geom->J,
837  coeff, pa_data);
838  }
839  else if (((trial_curl && test_div) || (trial_div && test_curl)) &&
840  test_fel->GetOrder() == trial_fel->GetOrder())
841  {
842  if (coeff_dim == 1)
843  {
844  PAHcurlL2Setup(nq, coeff_dim, ne, ir->GetWeights(), coeff, pa_data);
845  }
846  else
847  {
848  const bool tr = (trial_div && test_curl);
849  if (dim == 3)
850  PAHcurlHdivSetup3D(quad1D, coeff_dim, ne, tr, ir->GetWeights(),
851  geom->J, coeff, pa_data);
852  else
853  PAHcurlHdivSetup2D(quad1D, coeff_dim, ne, tr, ir->GetWeights(),
854  geom->J, coeff, pa_data);
855  }
856  }
857  else
858  {
859  MFEM_ABORT("Unknown kernel.");
860  }
861 }
862 
863 void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
864 {
865  if (dim == 3)
866  {
867  if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
868  {
869  if (Device::Allows(Backend::DEVICE_MASK))
870  {
871  const int ID = (dofs1D << 4) | quad1D;
872  switch (ID)
873  {
874  case 0x23: return SmemPAHcurlMassAssembleDiagonal3D<2,3>(dofs1D, quad1D, ne,
875  symmetric,
876  mapsO->B, mapsC->B, pa_data, diag);
877  case 0x34: return SmemPAHcurlMassAssembleDiagonal3D<3,4>(dofs1D, quad1D, ne,
878  symmetric,
879  mapsO->B, mapsC->B, pa_data, diag);
880  case 0x45: return SmemPAHcurlMassAssembleDiagonal3D<4,5>(dofs1D, quad1D, ne,
881  symmetric,
882  mapsO->B, mapsC->B, pa_data, diag);
883  case 0x56: return SmemPAHcurlMassAssembleDiagonal3D<5,6>(dofs1D, quad1D, ne,
884  symmetric,
885  mapsO->B, mapsC->B, pa_data, diag);
886  default: return SmemPAHcurlMassAssembleDiagonal3D(dofs1D, quad1D, ne, symmetric,
887  mapsO->B, mapsC->B, pa_data, diag);
888  }
889  }
890  else
891  PAHcurlMassAssembleDiagonal3D(dofs1D, quad1D, ne, symmetric,
892  mapsO->B, mapsC->B, pa_data, diag);
893  }
894  else if (trial_fetype == mfem::FiniteElement::DIV &&
895  test_fetype == trial_fetype)
896  {
897  PAHdivMassAssembleDiagonal3D(dofs1D, quad1D, ne, symmetric,
898  mapsO->B, mapsC->B, pa_data, diag);
899  }
900  else
901  {
902  MFEM_ABORT("Unknown kernel.");
903  }
904  }
905  else // 2D
906  {
907  if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
908  {
909  PAHcurlMassAssembleDiagonal2D(dofs1D, quad1D, ne, symmetric,
910  mapsO->B, mapsC->B, pa_data, diag);
911  }
912  else if (trial_fetype == mfem::FiniteElement::DIV &&
913  test_fetype == trial_fetype)
914  {
915  PAHdivMassAssembleDiagonal2D(dofs1D, quad1D, ne, symmetric,
916  mapsO->B, mapsC->B, pa_data, diag);
917  }
918  else
919  {
920  MFEM_ABORT("Unknown kernel.");
921  }
922  }
923 }
924 
925 void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
926 {
927  const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
928  const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
929  const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
930  const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
931 
932  if (dim == 3)
933  {
934  if (trial_curl && test_curl)
935  {
936  if (Device::Allows(Backend::DEVICE_MASK))
937  {
938  const int ID = (dofs1D << 4) | quad1D;
939  switch (ID)
940  {
941  case 0x23: return SmemPAHcurlMassApply3D<2,3>(dofs1D, quad1D, ne, symmetric,
942  mapsO->B,
943  mapsC->B, mapsO->Bt,
944  mapsC->Bt, pa_data, x, y);
945  case 0x34: return SmemPAHcurlMassApply3D<3,4>(dofs1D, quad1D, ne, symmetric,
946  mapsO->B,
947  mapsC->B, mapsO->Bt,
948  mapsC->Bt, pa_data, x, y);
949  case 0x45: return SmemPAHcurlMassApply3D<4,5>(dofs1D, quad1D, ne, symmetric,
950  mapsO->B,
951  mapsC->B, mapsO->Bt,
952  mapsC->Bt, pa_data, x, y);
953  case 0x56: return SmemPAHcurlMassApply3D<5,6>(dofs1D, quad1D, ne, symmetric,
954  mapsO->B,
955  mapsC->B, mapsO->Bt,
956  mapsC->Bt, pa_data, x, y);
957  default: return SmemPAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric, mapsO->B,
958  mapsC->B,
959  mapsO->Bt, mapsC->Bt, pa_data, x, y);
960  }
961  }
962  else
963  PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B, mapsO->Bt,
964  mapsC->Bt, pa_data, x, y);
965  }
966  else if (trial_div && test_div)
967  {
968  PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B, mapsO->Bt,
969  mapsC->Bt, pa_data, x, y);
970  }
971  else if (trial_curl && test_div)
972  {
973  const bool scalarCoeff = !(DQ || MQ);
974  PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
975  true, false, mapsO->B, mapsC->B, mapsOtest->Bt,
976  mapsCtest->Bt, pa_data, x, y);
977  }
978  else if (trial_div && test_curl)
979  {
980  const bool scalarCoeff = !(DQ || MQ);
981  PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
982  false, false, mapsO->B, mapsC->B, mapsOtest->Bt,
983  mapsCtest->Bt, pa_data, x, y);
984  }
985  else
986  {
987  MFEM_ABORT("Unknown kernel.");
988  }
989  }
990  else // 2D
991  {
992  if (trial_curl && test_curl)
993  {
994  PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
995  mapsO->Bt, mapsC->Bt, pa_data, x, y);
996  }
997  else if (trial_div && test_div)
998  {
999  PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B, mapsO->Bt,
1000  mapsC->Bt, pa_data, x, y);
1001  }
1002  else if ((trial_curl && test_div) || (trial_div && test_curl))
1003  {
1004  const bool scalarCoeff = !(DQ || MQ);
1005  PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
1006  trial_curl, false, mapsO->B, mapsC->B,
1007  mapsOtest->Bt, mapsCtest->Bt, pa_data, x, y);
1008  }
1009  else
1010  {
1011  MFEM_ABORT("Unknown kernel.");
1012  }
1013  }
1014 }
1015 
1016 void VectorFEMassIntegrator::AddMultTransposePA(const Vector &x,
1017  Vector &y) const
1018 {
1019  const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
1020  const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
1021  const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
1022  const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
1023 
1024  bool symmetricSpaces = true;
1025 
1026  if (dim == 3 && ((trial_div && test_curl) || (trial_curl && test_div)))
1027  {
1028  const bool scalarCoeff = !(DQ || MQ);
1029  PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
1030  trial_div, true, mapsO->B, mapsC->B, mapsOtest->Bt,
1031  mapsCtest->Bt, pa_data, x, y);
1032  symmetricSpaces = false;
1033  }
1034  else if (dim == 2 && ((trial_curl && test_div) || (trial_div && test_curl)))
1035  {
1036  const bool scalarCoeff = !(DQ || MQ);
1037  PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
1038  !trial_curl, true, mapsO->B, mapsC->B, mapsOtest->Bt,
1039  mapsCtest->Bt, pa_data, x, y);
1040  symmetricSpaces = false;
1041  }
1042 
1043  if (symmetricSpaces)
1044  {
1045  if (MQ && dynamic_cast<SymmetricMatrixCoefficient*>(MQ) == NULL)
1046  {
1047  MFEM_ABORT("VectorFEMassIntegrator transpose not implemented for asymmetric MatrixCoefficient");
1048  }
1049 
1050  this->AddMultPA(x, y);
1051  }
1052 }
1053 
1054 void MixedVectorGradientIntegrator::AssemblePA(const FiniteElementSpace
1055  &trial_fes,
1056  const FiniteElementSpace &test_fes)
1057 {
1058  // Assumes tensor-product elements, with a vector test space and H^1 trial space.
1059  Mesh *mesh = trial_fes.GetMesh();
1060  const FiniteElement *trial_fel = trial_fes.GetFE(0);
1061  const FiniteElement *test_fel = test_fes.GetFE(0);
1062 
1063  const NodalTensorFiniteElement *trial_el =
1064  dynamic_cast<const NodalTensorFiniteElement*>(trial_fel);
1065  MFEM_VERIFY(trial_el != NULL, "Only NodalTensorFiniteElement is supported!");
1066 
1067  const VectorTensorFiniteElement *test_el =
1068  dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
1069  MFEM_VERIFY(test_el != NULL, "Only VectorTensorFiniteElement is supported!");
1070 
1071  const IntegrationRule *ir
1072  = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
1073  *mesh->GetElementTransformation(0));
1074  const int dims = trial_el->GetDim();
1075  MFEM_VERIFY(dims == 2 || dims == 3, "");
1076 
1077  const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
1078  const int nq = ir->GetNPoints();
1079  dim = mesh->Dimension();
1080  MFEM_VERIFY(dim == 2 || dim == 3, "");
1081 
1082  MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(), "");
1083 
1084  ne = trial_fes.GetNE();
1085  geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
1086  mapsC = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
1087  mapsO = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
1088  dofs1D = mapsC->ndof;
1089  quad1D = mapsC->nqpt;
1090 
1091  MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
1092 
1093  pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
1094 
1095  QuadratureSpace qs(*mesh, *ir);
1096  CoefficientVector coeff(Q, qs, CoefficientStorage::FULL);
1097 
1098  // Use the same setup functions as VectorFEMassIntegrator.
1099  if (test_el->GetDerivType() == mfem::FiniteElement::CURL && dim == 3)
1100  {
1101  PADiffusionSetup3D(quad1D, 1, ne, ir->GetWeights(), geom->J,
1102  coeff, pa_data);
1103  }
1104  else if (test_el->GetDerivType() == mfem::FiniteElement::CURL && dim == 2)
1105  {
1106  PADiffusionSetup2D<2>(quad1D, 1, ne, ir->GetWeights(), geom->J,
1107  coeff, pa_data);
1108  }
1109  else
1110  {
1111  MFEM_ABORT("Unknown kernel.");
1112  }
1113 }
1114 
1115 void MixedVectorGradientIntegrator::AddMultPA(const Vector &x, Vector &y) const
1116 {
1117  if (dim == 3)
1118  PAHcurlH1Apply3D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
1119  mapsO->Bt, mapsC->Bt, pa_data, x, y);
1120  else if (dim == 2)
1121  PAHcurlH1Apply2D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
1122  mapsO->Bt, mapsC->Bt, pa_data, x, y);
1123  else
1124  {
1125  MFEM_ABORT("Unsupported dimension!");
1126  }
1127 }
1128 
1129 void MixedVectorGradientIntegrator::AddMultTransposePA(const Vector &x,
1130  Vector &y) const
1131 {
1132  if (dim == 3)
1133  PAHcurlH1ApplyTranspose3D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
1134  mapsC->Bt, mapsC->Gt, pa_data, x, y);
1135  else if (dim == 2)
1136  PAHcurlH1ApplyTranspose2D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
1137  mapsC->Bt, mapsC->Gt, pa_data, x, y);
1138  else
1139  {
1140  MFEM_ABORT("Unsupported dimension!");
1141  }
1142 }
1143 
1144 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
void PAHdivMassApply(const int dim, 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 &op, const Vector &x, Vector &y)
void PAHcurlH1Apply3D(const int D1D, const int Q1D, const int NE, const Array< double > &bc, const Array< double > &gc, const Array< double > &bot, const Array< double > &bct, const Vector &pa_data, const Vector &x, Vector &y)
void PAHcurlMassApply3D(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)
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:448
void PAHcurlHdivMassApply3D(const int D1D, const int D1Dtest, const int Q1D, const int NE, const bool scalarCoeff, const bool trialHcurl, const bool transpose, const Array< double > &Bo_, const Array< double > &Bc_, const Array< double > &Bot_, const Array< double > &Bct_, const Vector &op_, const Vector &x_, Vector &y_)
void PAHcurlL2Setup(const int NQ, const int coeffDim, const int NE, const Array< double > &w, Vector &coeff, Vector &op)
void SmemPAHcurlMassApply3D(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)
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)
void PAHdivMassAssembleDiagonal3D(const int D1D, const int Q1D, const int NE, const bool symmetric, const Array< double > &Bo_, const Array< double > &Bc_, const Vector &op_, Vector &diag_)
const int MAX_Q1D
Definition: forall.hpp:29
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:456
constexpr int HCURL_MAX_D1D
Definition: bilininteg.hpp:24
void PAHdivSetup2D(const int Q1D, const int coeffDim, const int NE, const Array< double > &w, const Vector &j, Vector &coeff_, Vector &op)
void PAHcurlH1ApplyTranspose2D(const int D1D, const int Q1D, const int NE, const Array< double > &bc, const Array< double > &bo, const Array< double > &bct, const Array< double > &gct, const Vector &pa_data, const Vector &x, Vector &y)
void PAHdivMassAssembleDiagonal2D(const int D1D, const int Q1D, const int NE, const bool symmetric, const Array< double > &Bo_, const Array< double > &Bc_, const Vector &op_, Vector &diag_)
Base class for Matrix Coefficients that optionally depend on time and space.
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)
void PADiffusionSetup3D(const int Q1D, const int coeffDim, const int NE, const Array< double > &w, const Vector &j, const Vector &c, Vector &d)
int dim
Definition: ex24.cpp:53
const int MAX_D1D
Definition: forall.hpp:28
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.
Definition: hooke.cpp:45
void SmemPAHcurlMassAssembleDiagonal3D(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)
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)
Vector data type.
Definition: vector.hpp:60
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
void PAHcurlHdivSetup2D(const int Q1D, const int coeffDim, const int NE, const bool transpose, const Array< double > &w_, const Vector &j, Vector &coeff_, Vector &op)
void PAHcurlH1Apply2D(const int D1D, const int Q1D, const int NE, const Array< double > &bc, const Array< double > &gc, const Array< double > &bot, const Array< double > &bct, const Vector &pa_data, const Vector &x, Vector &y)
constexpr int HCURL_MAX_Q1D
Definition: bilininteg.hpp:26
void PAHdivSetup3D(const int Q1D, const int coeffDim, const int NE, const Array< double > &w, const Vector &j, Vector &coeff_, Vector &op)
void PAHcurlH1ApplyTranspose3D(const int D1D, const int Q1D, const int NE, const Array< double > &bc, const Array< double > &bo, const Array< double > &bct, const Array< double > &gct, const Vector &pa_data, const Vector &x, Vector &y)