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