MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_mixedvecgrad_pa.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2026, 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#include "../qfunction.hpp"
16#include "bilininteg_diffusion_kernels.hpp"
17
18namespace mfem
19{
20
21namespace internal
22{
23
24// Setup for MixedVectorGradientIntegrator with an H(div) test space.
25// The stored operator represents (w/detJ) * J^T * Q * adj(J)^T, where Q is the
26// (optional) coefficient in physical space and adj(J) is the adjugate of the
27// physical Jacobian.
28static void PAMixedVectorGradientSetupHdiv2D(const int Q1D,
29 const int coeffDim,
30 const int NE,
31 const Array<real_t> &w,
32 const Vector &j,
33 const Vector &c,
34 Vector &op)
35{
36 MFEM_VERIFY(coeffDim == 1 || coeffDim == 2 || coeffDim == 4,
37 "Unsupported coefficient dimension for 2D MixedVectorGradient PA setup.");
38 const bool const_c = c.Size() == coeffDim;
39 const auto W = Reshape(w.Read(), Q1D, Q1D);
40 const auto J = Reshape(j.Read(), Q1D, Q1D, 2, 2, NE);
41 const auto C = const_c ? Reshape(c.Read(), coeffDim, 1, 1, 1)
42 : Reshape(c.Read(), coeffDim, Q1D, Q1D, NE);
43 auto O = Reshape(op.Write(), Q1D, Q1D, 4, NE);
44
45 auto get_coeff = [const_c] MFEM_HOST_DEVICE
46 (const decltype(C) &C, int i, int qx, int qy, int e)
47 {
48 return const_c ? C(i,0,0,0) : C(i,qx,qy,e);
49 };
50
51 mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
52 {
53 MFEM_FOREACH_THREAD(qx, x, Q1D)
54 {
55 MFEM_FOREACH_THREAD(qy, y, Q1D)
56 {
57 const real_t J11 = J(qx,qy,0,0,e);
58 const real_t J21 = J(qx,qy,1,0,e);
59 const real_t J12 = J(qx,qy,0,1,e);
60 const real_t J22 = J(qx,qy,1,1,e);
61 const real_t detJ = (J11*J22) - (J21*J12);
62 const real_t w_detJ = W(qx,qy) / detJ;
63
64 real_t M11, M12, M21, M22;
65 if (coeffDim == 4) // matrix coefficient
66 {
67 // NOTE: values are interpreted in row-major ordering.
68 M11 = get_coeff(C,0,qx,qy,e);
69 M12 = get_coeff(C,1,qx,qy,e);
70 M21 = get_coeff(C,2,qx,qy,e);
71 M22 = get_coeff(C,3,qx,qy,e);
72 }
73 else if (coeffDim == 2) // diagonal coefficient
74 {
75 M11 = get_coeff(C,0,qx,qy,e);
76 M22 = get_coeff(C,1,qx,qy,e);
77 M12 = 0.0;
78 M21 = 0.0;
79 }
80 else // scalar coefficient
81 {
82 M11 = get_coeff(C,0,qx,qy,e);
83 M22 = M11;
84 M12 = 0.0;
85 M21 = 0.0;
86 }
87
88 // R = Q * adj(J)^T, without detJ.
89 const real_t R11 = M11*J22 - M12*J12;
90 const real_t R12 = -M11*J21 + M12*J11;
91 const real_t R21 = M21*J22 - M22*J12;
92 const real_t R22 = -M21*J21 + M22*J11;
93
94 // O = (w/detJ) * J^T * R.
95 const real_t O11 = w_detJ * (J11*R11 + J21*R21);
96 const real_t O12 = w_detJ * (J11*R12 + J21*R22);
97 const real_t O21 = w_detJ * (J12*R11 + J22*R21);
98 const real_t O22 = w_detJ * (J12*R12 + J22*R22);
99
100 // Store as (1,1), (2,1), (1,2), (2,2).
101 O(qx,qy,0,e) = O11;
102 O(qx,qy,1,e) = O21;
103 O(qx,qy,2,e) = O12;
104 O(qx,qy,3,e) = O22;
105 }
106 }
107 });
108}
109
110static void PAMixedVectorGradientSetupHdiv3D(const int Q1D,
111 const int coeffDim,
112 const int NE,
113 const Array<real_t> &w,
114 const Vector &j,
115 const Vector &c,
116 Vector &op)
117{
118 MFEM_VERIFY(coeffDim == 1 || coeffDim == 3 || coeffDim == 9,
119 "Unsupported coefficient dimension for 3D MixedVectorGradient PA setup.");
120 const bool const_c = c.Size() == coeffDim;
121 const auto W = Reshape(w.Read(), Q1D, Q1D, Q1D);
122 const auto J = Reshape(j.Read(), Q1D, Q1D, Q1D, 3, 3, NE);
123 const auto C = const_c ? Reshape(c.Read(), coeffDim, 1, 1, 1, 1)
124 : Reshape(c.Read(), coeffDim, Q1D, Q1D, Q1D, NE);
125 auto O = Reshape(op.Write(), Q1D, Q1D, Q1D, 9, NE);
126
127 auto get_coeff = [const_c] MFEM_HOST_DEVICE
128 (const decltype(C) &C, int i, int qx, int qy, int qz, int e)
129 {
130 return const_c ? C(i,0,0,0,0) : C(i,qx,qy,qz,e);
131 };
132
133 mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
134 {
135 MFEM_FOREACH_THREAD(qx, x, Q1D)
136 {
137 MFEM_FOREACH_THREAD(qy, y, Q1D)
138 {
139 MFEM_FOREACH_THREAD(qz, z, Q1D)
140 {
141 const real_t J11 = J(qx,qy,qz,0,0,e);
142 const real_t J21 = J(qx,qy,qz,1,0,e);
143 const real_t J31 = J(qx,qy,qz,2,0,e);
144 const real_t J12 = J(qx,qy,qz,0,1,e);
145 const real_t J22 = J(qx,qy,qz,1,1,e);
146 const real_t J32 = J(qx,qy,qz,2,1,e);
147 const real_t J13 = J(qx,qy,qz,0,2,e);
148 const real_t J23 = J(qx,qy,qz,1,2,e);
149 const real_t J33 = J(qx,qy,qz,2,2,e);
150
151 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
152 J21 * (J12 * J33 - J32 * J13) +
153 J31 * (J12 * J23 - J22 * J13);
154 const real_t w_detJ = W(qx,qy,qz) / detJ;
155
156 // adj(J)
157 const real_t A11 = (J22 * J33) - (J23 * J32);
158 const real_t A12 = (J32 * J13) - (J12 * J33);
159 const real_t A13 = (J12 * J23) - (J22 * J13);
160 const real_t A21 = (J31 * J23) - (J21 * J33);
161 const real_t A22 = (J11 * J33) - (J13 * J31);
162 const real_t A23 = (J21 * J13) - (J11 * J23);
163 const real_t A31 = (J21 * J32) - (J31 * J22);
164 const real_t A32 = (J31 * J12) - (J11 * J32);
165 const real_t A33 = (J11 * J22) - (J12 * J21);
166
167 real_t M11, M12, M13, M21, M22, M23, M31, M32, M33;
168 if (coeffDim == 9)
169 {
170 // NOTE: values are interpreted in row-major ordering.
171 M11 = get_coeff(C,0,qx,qy,qz,e);
172 M12 = get_coeff(C,1,qx,qy,qz,e);
173 M13 = get_coeff(C,2,qx,qy,qz,e);
174 M21 = get_coeff(C,3,qx,qy,qz,e);
175 M22 = get_coeff(C,4,qx,qy,qz,e);
176 M23 = get_coeff(C,5,qx,qy,qz,e);
177 M31 = get_coeff(C,6,qx,qy,qz,e);
178 M32 = get_coeff(C,7,qx,qy,qz,e);
179 M33 = get_coeff(C,8,qx,qy,qz,e);
180 }
181 else if (coeffDim == 3)
182 {
183 M11 = get_coeff(C,0,qx,qy,qz,e);
184 M22 = get_coeff(C,1,qx,qy,qz,e);
185 M33 = get_coeff(C,2,qx,qy,qz,e);
186 M12 = M13 = M21 = M23 = M31 = M32 = 0.0;
187 }
188 else
189 {
190 M11 = get_coeff(C,0,qx,qy,qz,e);
191 M22 = M11;
192 M33 = M11;
193 M12 = M13 = M21 = M23 = M31 = M32 = 0.0;
194 }
195
196 // R = Q * adj(J)^T, without detJ.
197 const real_t R11 = M11*A11 + M12*A12 + M13*A13;
198 const real_t R12 = M11*A21 + M12*A22 + M13*A23;
199 const real_t R13 = M11*A31 + M12*A32 + M13*A33;
200 const real_t R21 = M21*A11 + M22*A12 + M23*A13;
201 const real_t R22 = M21*A21 + M22*A22 + M23*A23;
202 const real_t R23 = M21*A31 + M22*A32 + M23*A33;
203 const real_t R31 = M31*A11 + M32*A12 + M33*A13;
204 const real_t R32 = M31*A21 + M32*A22 + M33*A23;
205 const real_t R33 = M31*A31 + M32*A32 + M33*A33;
206
207 // O = (w/detJ) * J^T * R.
208 const real_t O11 = w_detJ * (J11*R11 + J21*R21 + J31*R31);
209 const real_t O12 = w_detJ * (J11*R12 + J21*R22 + J31*R32);
210 const real_t O13 = w_detJ * (J11*R13 + J21*R23 + J31*R33);
211 const real_t O21 = w_detJ * (J12*R11 + J22*R21 + J32*R31);
212 const real_t O22 = w_detJ * (J12*R12 + J22*R22 + J32*R32);
213 const real_t O23 = w_detJ * (J12*R13 + J22*R23 + J32*R33);
214 const real_t O31 = w_detJ * (J13*R11 + J23*R21 + J33*R31);
215 const real_t O32 = w_detJ * (J13*R12 + J23*R22 + J33*R32);
216 const real_t O33 = w_detJ * (J13*R13 + J23*R23 + J33*R33);
217
218 // Store row-major.
219 O(qx,qy,qz,0,e) = O11;
220 O(qx,qy,qz,1,e) = O12;
221 O(qx,qy,qz,2,e) = O13;
222 O(qx,qy,qz,3,e) = O21;
223 O(qx,qy,qz,4,e) = O22;
224 O(qx,qy,qz,5,e) = O23;
225 O(qx,qy,qz,6,e) = O31;
226 O(qx,qy,qz,7,e) = O32;
227 O(qx,qy,qz,8,e) = O33;
228 }
229 }
230 }
231 });
232}
233
234} // namespace internal
235
236// Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
237// integrated against H(curl) test functions corresponding to y.
238static void PAHcurlH1Apply2D(const int D1D,
239 const int Q1D,
240 const int NE,
241 const Array<real_t> &bc,
242 const Array<real_t> &gc,
243 const Array<real_t> &bot,
244 const Array<real_t> &bct,
245 const int op_entries,
246 const Vector &pa_data,
247 const Vector &x,
248 Vector &y)
249{
250 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
251 "Error: D1D > MAX_D1D");
252 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
253 "Error: Q1D > MAX_Q1D");
254
255 auto Bc = Reshape(bc.Read(), Q1D, D1D);
256 auto Gc = Reshape(gc.Read(), Q1D, D1D);
257 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
258 auto Bct = Reshape(bct.Read(), D1D, Q1D);
259 auto op = Reshape(pa_data.Read(), Q1D, Q1D, op_entries, NE);
260 auto X = Reshape(x.Read(), D1D, D1D, NE);
261 auto Y = Reshape(y.ReadWrite(), 2*(D1D-1)*D1D, NE);
262
263 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
264 {
265 constexpr static int VDIM = 2;
266 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
267 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
268
269 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
270
271 for (int qy = 0; qy < Q1D; ++qy)
272 {
273 for (int qx = 0; qx < Q1D; ++qx)
274 {
275 for (int c = 0; c < VDIM; ++c)
276 {
277 mass[qy][qx][c] = 0.0;
278 }
279 }
280 }
281
282 for (int dy = 0; dy < D1D; ++dy)
283 {
284 real_t gradX[MAX_Q1D][2];
285 for (int qx = 0; qx < Q1D; ++qx)
286 {
287 gradX[qx][0] = 0.0;
288 gradX[qx][1] = 0.0;
289 }
290 for (int dx = 0; dx < D1D; ++dx)
291 {
292 const real_t s = X(dx,dy,e);
293 for (int qx = 0; qx < Q1D; ++qx)
294 {
295 gradX[qx][0] += s * Bc(qx,dx);
296 gradX[qx][1] += s * Gc(qx,dx);
297 }
298 }
299 for (int qy = 0; qy < Q1D; ++qy)
300 {
301 const real_t wy = Bc(qy,dy);
302 const real_t wDy = Gc(qy,dy);
303 for (int qx = 0; qx < Q1D; ++qx)
304 {
305 const real_t wx = gradX[qx][0];
306 const real_t wDx = gradX[qx][1];
307 mass[qy][qx][0] += wDx * wy;
308 mass[qy][qx][1] += wx * wDy;
309 }
310 }
311 }
312
313 // Apply D operator.
314 for (int qy = 0; qy < Q1D; ++qy)
315 {
316 for (int qx = 0; qx < Q1D; ++qx)
317 {
318 const real_t massX = mass[qy][qx][0];
319 const real_t massY = mass[qy][qx][1];
320 if (op_entries == 3)
321 {
322 const real_t O11 = op(qx,qy,0,e);
323 const real_t O12 = op(qx,qy,1,e);
324 const real_t O22 = op(qx,qy,2,e);
325 mass[qy][qx][0] = (O11*massX)+(O12*massY);
326 mass[qy][qx][1] = (O12*massX)+(O22*massY);
327 }
328 else
329 {
330 // Non-symmetric operator stored as (1,1), (2,1), (1,2), (2,2).
331 const real_t O11 = op(qx,qy,0,e);
332 const real_t O21 = op(qx,qy,1,e);
333 const real_t O12 = op(qx,qy,2,e);
334 const real_t O22 = op(qx,qy,3,e);
335 mass[qy][qx][0] = (O11*massX)+(O12*massY);
336 mass[qy][qx][1] = (O21*massX)+(O22*massY);
337 }
338 }
339 }
340
341 for (int qy = 0; qy < Q1D; ++qy)
342 {
343 int osc = 0;
344
345 for (int c = 0; c < VDIM; ++c) // loop over x, y components
346 {
347 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
348 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
349
350 real_t massX[MAX_D1D];
351 for (int dx = 0; dx < D1Dx; ++dx)
352 {
353 massX[dx] = 0;
354 }
355 for (int qx = 0; qx < Q1D; ++qx)
356 {
357 for (int dx = 0; dx < D1Dx; ++dx)
358 {
359 massX[dx] += mass[qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
360 }
361 }
362
363 for (int dy = 0; dy < D1Dy; ++dy)
364 {
365 const real_t wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
366
367 for (int dx = 0; dx < D1Dx; ++dx)
368 {
369 Y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
370 }
371 }
372
373 osc += D1Dx * D1Dy;
374 } // loop c
375 }
376 }); // end of element loop
377}
378
379// Apply to x corresponding to DOFs in H(curl), integrated
380// against gradients of H^1 functions corresponding to y.
381static void PAHcurlH1ApplyTranspose2D(const int D1D,
382 const int Q1D,
383 const int NE,
384 const Array<real_t> &bc,
385 const Array<real_t> &bo,
386 const Array<real_t> &bct,
387 const Array<real_t> &gct,
388 const int op_entries,
389 const Vector &pa_data,
390 const Vector &x,
391 Vector &y)
392{
393 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
394 "Error: D1D > MAX_D1D");
395 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
396 "Error: Q1D > MAX_Q1D");
397 auto Bc = Reshape(bc.Read(), Q1D, D1D);
398 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
399 auto Bt = Reshape(bct.Read(), D1D, Q1D);
400 auto Gt = Reshape(gct.Read(), D1D, Q1D);
401 auto op = Reshape(pa_data.Read(), Q1D, Q1D, op_entries, NE);
402 auto X = Reshape(x.Read(), 2*(D1D-1)*D1D, NE);
403 auto Y = Reshape(y.ReadWrite(), D1D, D1D, NE);
404
405 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
406 {
407 constexpr static int VDIM = 2;
408 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
409 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
410
411 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
412
413 for (int qy = 0; qy < Q1D; ++qy)
414 {
415 for (int qx = 0; qx < Q1D; ++qx)
416 {
417 for (int c = 0; c < VDIM; ++c)
418 {
419 mass[qy][qx][c] = 0.0;
420 }
421 }
422 }
423
424 int osc = 0;
425
426 for (int c = 0; c < VDIM; ++c) // loop over x, y components
427 {
428 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
429 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
430
431 for (int dy = 0; dy < D1Dy; ++dy)
432 {
433 real_t 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 real_t t = X(dx + (dy * D1Dx) + osc, e);
442 for (int qx = 0; qx < Q1D; ++qx)
443 {
444 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
445 }
446 }
447
448 for (int qy = 0; qy < Q1D; ++qy)
449 {
450 const real_t wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
451 for (int qx = 0; qx < Q1D; ++qx)
452 {
453 mass[qy][qx][c] += massX[qx] * wy;
454 }
455 }
456 }
457
458 osc += D1Dx * D1Dy;
459 } // loop (c) over components
460
461 // Apply D operator.
462 for (int qy = 0; qy < Q1D; ++qy)
463 {
464 for (int qx = 0; qx < Q1D; ++qx)
465 {
466 const real_t massX = mass[qy][qx][0];
467 const real_t massY = mass[qy][qx][1];
468 if (op_entries == 3)
469 {
470 const real_t O11 = op(qx,qy,0,e);
471 const real_t O12 = op(qx,qy,1,e);
472 const real_t O22 = op(qx,qy,2,e);
473 mass[qy][qx][0] = (O11*massX)+(O12*massY);
474 mass[qy][qx][1] = (O12*massX)+(O22*massY);
475 }
476 else
477 {
478 // Non-symmetric operator stored as (1,1), (2,1), (1,2), (2,2).
479 // For transpose, apply D^T.
480 const real_t O11 = op(qx,qy,0,e);
481 const real_t O21 = op(qx,qy,1,e);
482 const real_t O12 = op(qx,qy,2,e);
483 const real_t O22 = op(qx,qy,3,e);
484 mass[qy][qx][0] = (O11*massX)+(O21*massY);
485 mass[qy][qx][1] = (O12*massX)+(O22*massY);
486 }
487 }
488 }
489
490 for (int qy = 0; qy < Q1D; ++qy)
491 {
492 real_t gradX[MAX_D1D][2];
493 for (int dx = 0; dx < D1D; ++dx)
494 {
495 gradX[dx][0] = 0;
496 gradX[dx][1] = 0;
497 }
498 for (int qx = 0; qx < Q1D; ++qx)
499 {
500 const real_t gX = mass[qy][qx][0];
501 const real_t gY = mass[qy][qx][1];
502 for (int dx = 0; dx < D1D; ++dx)
503 {
504 const real_t wx = Bt(dx,qx);
505 const real_t wDx = Gt(dx,qx);
506 gradX[dx][0] += gX * wDx;
507 gradX[dx][1] += gY * wx;
508 }
509 }
510 for (int dy = 0; dy < D1D; ++dy)
511 {
512 const real_t wy = Bt(dy,qy);
513 const real_t wDy = Gt(dy,qy);
514 for (int dx = 0; dx < D1D; ++dx)
515 {
516 Y(dx,dy,e) += ((gradX[dx][0] * wy) + (gradX[dx][1] * wDy));
517 }
518 }
519 }
520 }); // end of element loop
521}
522
523// Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
524// integrated against H(div) test functions corresponding to y.
525static void PAHdivH1Apply2D(const int D1D,
526 const int Q1D,
527 const int NE,
528 const Array<real_t> &bc,
529 const Array<real_t> &gc,
530 const Array<real_t> &bot,
531 const Array<real_t> &bct,
532 const int op_entries,
533 const Vector &pa_data,
534 const Vector &x,
535 Vector &y)
536{
537 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
538 "Error: D1D > MAX_D1D");
539 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
540 "Error: Q1D > MAX_Q1D");
541
542 auto Bc = Reshape(bc.Read(), Q1D, D1D);
543 auto Gc = Reshape(gc.Read(), Q1D, D1D);
544 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
545 auto Bct = Reshape(bct.Read(), D1D, Q1D);
546 auto op = Reshape(pa_data.Read(), Q1D, Q1D, op_entries, NE);
547 auto X = Reshape(x.Read(), D1D, D1D, NE);
548 auto Y = Reshape(y.ReadWrite(), 2*(D1D-1)*D1D, NE);
549
550 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
551 {
552 constexpr static int VDIM = 2;
553 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
554 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
555
556 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
557 for (int qy = 0; qy < Q1D; ++qy)
558 {
559 for (int qx = 0; qx < Q1D; ++qx)
560 {
561 for (int c = 0; c < VDIM; ++c) { mass[qy][qx][c] = 0.0; }
562 }
563 }
564
565 for (int dy = 0; dy < D1D; ++dy)
566 {
567 real_t gradX[MAX_Q1D][2];
568 for (int qx = 0; qx < Q1D; ++qx) { gradX[qx][0] = 0.0; gradX[qx][1] = 0.0; }
569 for (int dx = 0; dx < D1D; ++dx)
570 {
571 const real_t s = X(dx,dy,e);
572 for (int qx = 0; qx < Q1D; ++qx)
573 {
574 gradX[qx][0] += s * Bc(qx,dx);
575 gradX[qx][1] += s * Gc(qx,dx);
576 }
577 }
578 for (int qy = 0; qy < Q1D; ++qy)
579 {
580 const real_t wy = Bc(qy,dy);
581 const real_t wDy = Gc(qy,dy);
582 for (int qx = 0; qx < Q1D; ++qx)
583 {
584 const real_t wx = gradX[qx][0];
585 const real_t wDx = gradX[qx][1];
586 mass[qy][qx][0] += wDx * wy;
587 mass[qy][qx][1] += wx * wDy;
588 }
589 }
590 }
591
592 // Apply D operator.
593 for (int qy = 0; qy < Q1D; ++qy)
594 {
595 for (int qx = 0; qx < Q1D; ++qx)
596 {
597 const real_t massX = mass[qy][qx][0];
598 const real_t massY = mass[qy][qx][1];
599 if (op_entries == 3)
600 {
601 const real_t O11 = op(qx,qy,0,e);
602 const real_t O12 = op(qx,qy,1,e);
603 const real_t O22 = op(qx,qy,2,e);
604 mass[qy][qx][0] = (O11*massX)+(O12*massY);
605 mass[qy][qx][1] = (O12*massX)+(O22*massY);
606 }
607 else
608 {
609 // Non-symmetric operator stored as (1,1), (2,1), (1,2), (2,2).
610 const real_t O11 = op(qx,qy,0,e);
611 const real_t O21 = op(qx,qy,1,e);
612 const real_t O12 = op(qx,qy,2,e);
613 const real_t O22 = op(qx,qy,3,e);
614 mass[qy][qx][0] = (O11*massX)+(O12*massY);
615 mass[qy][qx][1] = (O21*massX)+(O22*massY);
616 }
617 }
618 }
619
620 for (int qy = 0; qy < Q1D; ++qy)
621 {
622 int osc = 0;
623 for (int c = 0; c < VDIM; ++c) // x, y components
624 {
625 const int D1Dx = (c == 0) ? D1D : D1D - 1;
626 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
627
628 real_t massX[MAX_D1D];
629 for (int dx = 0; dx < D1Dx; ++dx) { massX[dx] = 0.0; }
630
631 for (int qx = 0; qx < Q1D; ++qx)
632 {
633 for (int dx = 0; dx < D1Dx; ++dx)
634 {
635 massX[dx] += mass[qy][qx][c] * ((c == 0) ? Bct(dx,qx) : Bot(dx,qx));
636 }
637 }
638
639 for (int dy = 0; dy < D1Dy; ++dy)
640 {
641 const real_t wy = (c == 0) ? Bot(dy,qy) : Bct(dy,qy);
642 for (int dx = 0; dx < D1Dx; ++dx)
643 {
644 Y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
645 }
646 }
647
648 osc += D1Dx * D1Dy;
649 }
650 }
651 });
652}
653
654// Apply to x corresponding to DOFs in H(div), integrated
655// against gradients of H^1 functions corresponding to y.
656static void PAHdivH1ApplyTranspose2D(const int D1D,
657 const int Q1D,
658 const int NE,
659 const Array<real_t> &bc,
660 const Array<real_t> &bo,
661 const Array<real_t> &bct,
662 const Array<real_t> &gct,
663 const int op_entries,
664 const Vector &pa_data,
665 const Vector &x,
666 Vector &y)
667{
668 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
669 "Error: D1D > MAX_D1D");
670 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
671 "Error: Q1D > MAX_Q1D");
672
673 auto Bc = Reshape(bc.Read(), Q1D, D1D);
674 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
675 auto Bt = Reshape(bct.Read(), D1D, Q1D);
676 auto Gt = Reshape(gct.Read(), D1D, Q1D);
677 auto op = Reshape(pa_data.Read(), Q1D, Q1D, op_entries, NE);
678 auto X = Reshape(x.Read(), 2*(D1D-1)*D1D, NE);
679 auto Y = Reshape(y.ReadWrite(), D1D, D1D, NE);
680
681 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
682 {
683 constexpr static int VDIM = 2;
684 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
685 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
686
687 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
688 for (int qy = 0; qy < Q1D; ++qy)
689 {
690 for (int qx = 0; qx < Q1D; ++qx)
691 {
692 for (int c = 0; c < VDIM; ++c) { mass[qy][qx][c] = 0.0; }
693 }
694 }
695
696 int osc = 0;
697 for (int c = 0; c < VDIM; ++c)
698 {
699 const int D1Dx = (c == 0) ? D1D : D1D - 1;
700 const int D1Dy = (c == 0) ? D1D - 1 : D1D;
701
702 for (int dy = 0; dy < D1Dy; ++dy)
703 {
704 real_t massX[MAX_Q1D];
705 for (int qx = 0; qx < Q1D; ++qx) { massX[qx] = 0.0; }
706
707 for (int dx = 0; dx < D1Dx; ++dx)
708 {
709 const real_t t = X(dx + (dy * D1Dx) + osc, e);
710 for (int qx = 0; qx < Q1D; ++qx)
711 {
712 massX[qx] += t * ((c == 0) ? Bc(qx,dx) : Bo(qx,dx));
713 }
714 }
715
716 for (int qy = 0; qy < Q1D; ++qy)
717 {
718 const real_t wy = (c == 0) ? Bo(qy,dy) : Bc(qy,dy);
719 for (int qx = 0; qx < Q1D; ++qx)
720 {
721 mass[qy][qx][c] += massX[qx] * wy;
722 }
723 }
724 }
725
726 osc += D1Dx * D1Dy;
727 }
728
729 // Apply D operator.
730 for (int qy = 0; qy < Q1D; ++qy)
731 {
732 for (int qx = 0; qx < Q1D; ++qx)
733 {
734 const real_t massX = mass[qy][qx][0];
735 const real_t massY = mass[qy][qx][1];
736 if (op_entries == 3)
737 {
738 const real_t O11 = op(qx,qy,0,e);
739 const real_t O12 = op(qx,qy,1,e);
740 const real_t O22 = op(qx,qy,2,e);
741 mass[qy][qx][0] = (O11*massX)+(O12*massY);
742 mass[qy][qx][1] = (O12*massX)+(O22*massY);
743 }
744 else
745 {
746 // Non-symmetric operator stored as (1,1), (2,1), (1,2), (2,2).
747 // For transpose, apply D^T.
748 const real_t O11 = op(qx,qy,0,e);
749 const real_t O21 = op(qx,qy,1,e);
750 const real_t O12 = op(qx,qy,2,e);
751 const real_t O22 = op(qx,qy,3,e);
752 mass[qy][qx][0] = (O11*massX)+(O21*massY);
753 mass[qy][qx][1] = (O12*massX)+(O22*massY);
754 }
755 }
756 }
757
758 for (int qy = 0; qy < Q1D; ++qy)
759 {
760 real_t gradX[MAX_D1D][2];
761 for (int dx = 0; dx < D1D; ++dx) { gradX[dx][0] = 0.0; gradX[dx][1] = 0.0; }
762
763 for (int qx = 0; qx < Q1D; ++qx)
764 {
765 const real_t gX = mass[qy][qx][0];
766 const real_t gY = mass[qy][qx][1];
767 for (int dx = 0; dx < D1D; ++dx)
768 {
769 const real_t wx = Bt(dx,qx);
770 const real_t wDx = Gt(dx,qx);
771 gradX[dx][0] += gX * wDx;
772 gradX[dx][1] += gY * wx;
773 }
774 }
775
776 for (int dy = 0; dy < D1D; ++dy)
777 {
778 const real_t wy = Bt(dy,qy);
779 const real_t wDy = Gt(dy,qy);
780 for (int dx = 0; dx < D1D; ++dx)
781 {
782 Y(dx,dy,e) += ((gradX[dx][0] * wy) + (gradX[dx][1] * wDy));
783 }
784 }
785 }
786 });
787}
788
789// Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
790// integrated against H(curl) test functions corresponding to y.
791static void PAHcurlH1Apply3D(const int D1D,
792 const int Q1D,
793 const int NE,
794 const Array<real_t> &bc,
795 const Array<real_t> &gc,
796 const Array<real_t> &bot,
797 const Array<real_t> &bct,
798 const int op_entries,
799 const Vector &pa_data,
800 const Vector &x,
801 Vector &y)
802{
803 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
804 "Error: D1D > MAX_D1D");
805 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
806 "Error: Q1D > MAX_Q1D");
807
808 constexpr static int VDIM = 3;
809
810 auto Bc = Reshape(bc.Read(), Q1D, D1D);
811 auto Gc = Reshape(gc.Read(), Q1D, D1D);
812 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
813 auto Bct = Reshape(bct.Read(), D1D, Q1D);
814 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, op_entries, NE);
815 auto X = Reshape(x.Read(), D1D, D1D, D1D, NE);
816 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
817
818 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
819 {
820 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
821 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
822
823 real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
824
825 for (int qz = 0; qz < Q1D; ++qz)
826 {
827 for (int qy = 0; qy < Q1D; ++qy)
828 {
829 for (int qx = 0; qx < Q1D; ++qx)
830 {
831 for (int c = 0; c < VDIM; ++c)
832 {
833 mass[qz][qy][qx][c] = 0.0;
834 }
835 }
836 }
837 }
838
839 for (int dz = 0; dz < D1D; ++dz)
840 {
841 real_t gradXY[MAX_Q1D][MAX_Q1D][3];
842 for (int qy = 0; qy < Q1D; ++qy)
843 {
844 for (int qx = 0; qx < Q1D; ++qx)
845 {
846 gradXY[qy][qx][0] = 0.0;
847 gradXY[qy][qx][1] = 0.0;
848 gradXY[qy][qx][2] = 0.0;
849 }
850 }
851 for (int dy = 0; dy < D1D; ++dy)
852 {
853 real_t gradX[MAX_Q1D][2];
854 for (int qx = 0; qx < Q1D; ++qx)
855 {
856 gradX[qx][0] = 0.0;
857 gradX[qx][1] = 0.0;
858 }
859 for (int dx = 0; dx < D1D; ++dx)
860 {
861 const real_t s = X(dx,dy,dz,e);
862 for (int qx = 0; qx < Q1D; ++qx)
863 {
864 gradX[qx][0] += s * Bc(qx,dx);
865 gradX[qx][1] += s * Gc(qx,dx);
866 }
867 }
868 for (int qy = 0; qy < Q1D; ++qy)
869 {
870 const real_t wy = Bc(qy,dy);
871 const real_t wDy = Gc(qy,dy);
872 for (int qx = 0; qx < Q1D; ++qx)
873 {
874 const real_t wx = gradX[qx][0];
875 const real_t wDx = gradX[qx][1];
876 gradXY[qy][qx][0] += wDx * wy;
877 gradXY[qy][qx][1] += wx * wDy;
878 gradXY[qy][qx][2] += wx * wy;
879 }
880 }
881 }
882 for (int qz = 0; qz < Q1D; ++qz)
883 {
884 const real_t wz = Bc(qz,dz);
885 const real_t wDz = Gc(qz,dz);
886 for (int qy = 0; qy < Q1D; ++qy)
887 {
888 for (int qx = 0; qx < Q1D; ++qx)
889 {
890 mass[qz][qy][qx][0] += gradXY[qy][qx][0] * wz;
891 mass[qz][qy][qx][1] += gradXY[qy][qx][1] * wz;
892 mass[qz][qy][qx][2] += gradXY[qy][qx][2] * wDz;
893 }
894 }
895 }
896 }
897
898 // Apply D operator.
899 for (int qz = 0; qz < Q1D; ++qz)
900 {
901 for (int qy = 0; qy < Q1D; ++qy)
902 {
903 for (int qx = 0; qx < Q1D; ++qx)
904 {
905 const real_t massX = mass[qz][qy][qx][0];
906 const real_t massY = mass[qz][qy][qx][1];
907 const real_t massZ = mass[qz][qy][qx][2];
908 if (op_entries == 6)
909 {
910 const real_t O11 = op(qx,qy,qz,0,e);
911 const real_t O12 = op(qx,qy,qz,1,e);
912 const real_t O13 = op(qx,qy,qz,2,e);
913 const real_t O22 = op(qx,qy,qz,3,e);
914 const real_t O23 = op(qx,qy,qz,4,e);
915 const real_t O33 = op(qx,qy,qz,5,e);
916 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
917 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
918 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
919 }
920 else
921 {
922 // Non-symmetric operator stored row-major.
923 const real_t O11 = op(qx,qy,qz,0,e);
924 const real_t O12 = op(qx,qy,qz,1,e);
925 const real_t O13 = op(qx,qy,qz,2,e);
926 const real_t O21 = op(qx,qy,qz,3,e);
927 const real_t O22 = op(qx,qy,qz,4,e);
928 const real_t O23 = op(qx,qy,qz,5,e);
929 const real_t O31 = op(qx,qy,qz,6,e);
930 const real_t O32 = op(qx,qy,qz,7,e);
931 const real_t O33 = op(qx,qy,qz,8,e);
932 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
933 mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
934 mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
935 }
936 }
937 }
938 }
939
940 for (int qz = 0; qz < Q1D; ++qz)
941 {
942 real_t massXY[MAX_D1D][MAX_D1D];
943
944 int osc = 0;
945
946 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
947 {
948 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
949 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
950 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
951
952 for (int dy = 0; dy < D1Dy; ++dy)
953 {
954 for (int dx = 0; dx < D1Dx; ++dx)
955 {
956 massXY[dy][dx] = 0.0;
957 }
958 }
959 for (int qy = 0; qy < Q1D; ++qy)
960 {
961 real_t massX[MAX_D1D];
962 for (int dx = 0; dx < D1Dx; ++dx)
963 {
964 massX[dx] = 0;
965 }
966 for (int qx = 0; qx < Q1D; ++qx)
967 {
968 for (int dx = 0; dx < D1Dx; ++dx)
969 {
970 massX[dx] += mass[qz][qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
971 }
972 }
973 for (int dy = 0; dy < D1Dy; ++dy)
974 {
975 const real_t wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
976 for (int dx = 0; dx < D1Dx; ++dx)
977 {
978 massXY[dy][dx] += massX[dx] * wy;
979 }
980 }
981 }
982
983 for (int dz = 0; dz < D1Dz; ++dz)
984 {
985 const real_t wz = (c == 2) ? Bot(dz,qz) : Bct(dz,qz);
986 for (int dy = 0; dy < D1Dy; ++dy)
987 {
988 for (int dx = 0; dx < D1Dx; ++dx)
989 {
990 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += massXY[dy][dx] * wz;
991 }
992 }
993 }
994
995 osc += D1Dx * D1Dy * D1Dz;
996 } // loop c
997 } // loop qz
998 }); // end of element loop
999}
1000
1001// Apply to x corresponding to DOFs in H(curl), integrated
1002// against gradients of H^1 functions corresponding to y.
1003static void PAHcurlH1ApplyTranspose3D(const int D1D,
1004 const int Q1D,
1005 const int NE,
1006 const Array<real_t> &bc,
1007 const Array<real_t> &bo,
1008 const Array<real_t> &bct,
1009 const Array<real_t> &gct,
1010 const int op_entries,
1011 const Vector &pa_data,
1012 const Vector &x,
1013 Vector &y)
1014{
1015 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
1016 "Error: D1D > MAX_D1D");
1017 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
1018 "Error: Q1D > MAX_Q1D");
1019
1020 constexpr static int VDIM = 3;
1021
1022 auto Bc = Reshape(bc.Read(), Q1D, D1D);
1023 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
1024 auto Bt = Reshape(bct.Read(), D1D, Q1D);
1025 auto Gt = Reshape(gct.Read(), D1D, Q1D);
1026 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, op_entries, NE);
1027 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
1028 auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, NE);
1029
1030 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1031 {
1032 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
1033 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
1034
1035 real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
1036
1037 for (int qz = 0; qz < Q1D; ++qz)
1038 {
1039 for (int qy = 0; qy < Q1D; ++qy)
1040 {
1041 for (int qx = 0; qx < Q1D; ++qx)
1042 {
1043 for (int c = 0; c < VDIM; ++c)
1044 {
1045 mass[qz][qy][qx][c] = 0.0;
1046 }
1047 }
1048 }
1049 }
1050
1051 int osc = 0;
1052
1053 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
1054 {
1055 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
1056 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
1057 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
1058
1059 for (int dz = 0; dz < D1Dz; ++dz)
1060 {
1061 real_t massXY[MAX_Q1D][MAX_Q1D];
1062 for (int qy = 0; qy < Q1D; ++qy)
1063 {
1064 for (int qx = 0; qx < Q1D; ++qx)
1065 {
1066 massXY[qy][qx] = 0.0;
1067 }
1068 }
1069
1070 for (int dy = 0; dy < D1Dy; ++dy)
1071 {
1072 real_t massX[MAX_Q1D];
1073 for (int qx = 0; qx < Q1D; ++qx)
1074 {
1075 massX[qx] = 0.0;
1076 }
1077
1078 for (int dx = 0; dx < D1Dx; ++dx)
1079 {
1080 const real_t t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1081 for (int qx = 0; qx < Q1D; ++qx)
1082 {
1083 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
1084 }
1085 }
1086
1087 for (int qy = 0; qy < Q1D; ++qy)
1088 {
1089 const real_t wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
1090 for (int qx = 0; qx < Q1D; ++qx)
1091 {
1092 const real_t wx = massX[qx];
1093 massXY[qy][qx] += wx * wy;
1094 }
1095 }
1096 }
1097
1098 for (int qz = 0; qz < Q1D; ++qz)
1099 {
1100 const real_t wz = (c == 2) ? Bo(qz,dz) : Bc(qz,dz);
1101 for (int qy = 0; qy < Q1D; ++qy)
1102 {
1103 for (int qx = 0; qx < Q1D; ++qx)
1104 {
1105 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
1106 }
1107 }
1108 }
1109 }
1110
1111 osc += D1Dx * D1Dy * D1Dz;
1112 } // loop (c) over components
1113
1114 // Apply D operator.
1115 for (int qz = 0; qz < Q1D; ++qz)
1116 {
1117 for (int qy = 0; qy < Q1D; ++qy)
1118 {
1119 for (int qx = 0; qx < Q1D; ++qx)
1120 {
1121 const real_t massX = mass[qz][qy][qx][0];
1122 const real_t massY = mass[qz][qy][qx][1];
1123 const real_t massZ = mass[qz][qy][qx][2];
1124 if (op_entries == 6)
1125 {
1126 const real_t O11 = op(qx,qy,qz,0,e);
1127 const real_t O12 = op(qx,qy,qz,1,e);
1128 const real_t O13 = op(qx,qy,qz,2,e);
1129 const real_t O22 = op(qx,qy,qz,3,e);
1130 const real_t O23 = op(qx,qy,qz,4,e);
1131 const real_t O33 = op(qx,qy,qz,5,e);
1132 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
1133 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
1134 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
1135 }
1136 else
1137 {
1138 // Non-symmetric operator stored row-major. For transpose, apply D^T.
1139 const real_t O11 = op(qx,qy,qz,0,e);
1140 const real_t O12 = op(qx,qy,qz,1,e);
1141 const real_t O13 = op(qx,qy,qz,2,e);
1142 const real_t O21 = op(qx,qy,qz,3,e);
1143 const real_t O22 = op(qx,qy,qz,4,e);
1144 const real_t O23 = op(qx,qy,qz,5,e);
1145 const real_t O31 = op(qx,qy,qz,6,e);
1146 const real_t O32 = op(qx,qy,qz,7,e);
1147 const real_t O33 = op(qx,qy,qz,8,e);
1148 mass[qz][qy][qx][0] = (O11*massX)+(O21*massY)+(O31*massZ);
1149 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O32*massZ);
1150 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
1151 }
1152 }
1153 }
1154 }
1155
1156 for (int qz = 0; qz < Q1D; ++qz)
1157 {
1158 real_t gradXY[MAX_D1D][MAX_D1D][3];
1159 for (int dy = 0; dy < D1D; ++dy)
1160 {
1161 for (int dx = 0; dx < D1D; ++dx)
1162 {
1163 gradXY[dy][dx][0] = 0;
1164 gradXY[dy][dx][1] = 0;
1165 gradXY[dy][dx][2] = 0;
1166 }
1167 }
1168 for (int qy = 0; qy < Q1D; ++qy)
1169 {
1170 real_t gradX[MAX_D1D][3];
1171 for (int dx = 0; dx < D1D; ++dx)
1172 {
1173 gradX[dx][0] = 0;
1174 gradX[dx][1] = 0;
1175 gradX[dx][2] = 0;
1176 }
1177 for (int qx = 0; qx < Q1D; ++qx)
1178 {
1179 const real_t gX = mass[qz][qy][qx][0];
1180 const real_t gY = mass[qz][qy][qx][1];
1181 const real_t gZ = mass[qz][qy][qx][2];
1182 for (int dx = 0; dx < D1D; ++dx)
1183 {
1184 const real_t wx = Bt(dx,qx);
1185 const real_t wDx = Gt(dx,qx);
1186 gradX[dx][0] += gX * wDx;
1187 gradX[dx][1] += gY * wx;
1188 gradX[dx][2] += gZ * wx;
1189 }
1190 }
1191 for (int dy = 0; dy < D1D; ++dy)
1192 {
1193 const real_t wy = Bt(dy,qy);
1194 const real_t wDy = Gt(dy,qy);
1195 for (int dx = 0; dx < D1D; ++dx)
1196 {
1197 gradXY[dy][dx][0] += gradX[dx][0] * wy;
1198 gradXY[dy][dx][1] += gradX[dx][1] * wDy;
1199 gradXY[dy][dx][2] += gradX[dx][2] * wy;
1200 }
1201 }
1202 }
1203 for (int dz = 0; dz < D1D; ++dz)
1204 {
1205 const real_t wz = Bt(dz,qz);
1206 const real_t wDz = Gt(dz,qz);
1207 for (int dy = 0; dy < D1D; ++dy)
1208 {
1209 for (int dx = 0; dx < D1D; ++dx)
1210 {
1211 Y(dx,dy,dz,e) +=
1212 ((gradXY[dy][dx][0] * wz) +
1213 (gradXY[dy][dx][1] * wz) +
1214 (gradXY[dy][dx][2] * wDz));
1215 }
1216 }
1217 }
1218 } // loop qz
1219 }); // end of element loop
1220}
1221
1222// Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
1223// integrated against H(div) test functions corresponding to y.
1224static void PAHdivH1Apply3D(const int D1D,
1225 const int Q1D,
1226 const int NE,
1227 const Array<real_t> &bc,
1228 const Array<real_t> &gc,
1229 const Array<real_t> &bot,
1230 const Array<real_t> &bct,
1231 const int op_entries,
1232 const Vector &pa_data,
1233 const Vector &x,
1234 Vector &y)
1235{
1236 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
1237 "Error: D1D > MAX_D1D");
1238 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
1239 "Error: Q1D > MAX_Q1D");
1240
1241 constexpr static int VDIM = 3;
1242
1243 auto Bc = Reshape(bc.Read(), Q1D, D1D);
1244 auto Gc = Reshape(gc.Read(), Q1D, D1D);
1245 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
1246 auto Bct = Reshape(bct.Read(), D1D, Q1D);
1247 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, op_entries, NE);
1248 auto X = Reshape(x.Read(), D1D, D1D, D1D, NE);
1249 auto Y = Reshape(y.ReadWrite(), 3*D1D*(D1D-1)*(D1D-1), NE);
1250
1251 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1252 {
1253 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
1254 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
1255
1256 real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
1257
1258 for (int qz = 0; qz < Q1D; ++qz)
1259 {
1260 for (int qy = 0; qy < Q1D; ++qy)
1261 {
1262 for (int qx = 0; qx < Q1D; ++qx)
1263 {
1264 for (int c = 0; c < VDIM; ++c) { mass[qz][qy][qx][c] = 0.0; }
1265 }
1266 }
1267 }
1268
1269 for (int dz = 0; dz < D1D; ++dz)
1270 {
1271 real_t gradXY[MAX_Q1D][MAX_Q1D][3];
1272 for (int qy = 0; qy < Q1D; ++qy)
1273 {
1274 for (int qx = 0; qx < Q1D; ++qx)
1275 {
1276 gradXY[qy][qx][0] = 0.0;
1277 gradXY[qy][qx][1] = 0.0;
1278 gradXY[qy][qx][2] = 0.0;
1279 }
1280 }
1281 for (int dy = 0; dy < D1D; ++dy)
1282 {
1283 real_t gradX[MAX_Q1D][2];
1284 for (int qx = 0; qx < Q1D; ++qx)
1285 {
1286 gradX[qx][0] = 0.0;
1287 gradX[qx][1] = 0.0;
1288 }
1289 for (int dx = 0; dx < D1D; ++dx)
1290 {
1291 const real_t s = X(dx,dy,dz,e);
1292 for (int qx = 0; qx < Q1D; ++qx)
1293 {
1294 gradX[qx][0] += s * Bc(qx,dx);
1295 gradX[qx][1] += s * Gc(qx,dx);
1296 }
1297 }
1298 for (int qy = 0; qy < Q1D; ++qy)
1299 {
1300 const real_t wy = Bc(qy,dy);
1301 const real_t wDy = Gc(qy,dy);
1302 for (int qx = 0; qx < Q1D; ++qx)
1303 {
1304 const real_t wx = gradX[qx][0];
1305 const real_t wDx = gradX[qx][1];
1306 gradXY[qy][qx][0] += wDx * wy;
1307 gradXY[qy][qx][1] += wx * wDy;
1308 gradXY[qy][qx][2] += wx * wy;
1309 }
1310 }
1311 }
1312 for (int qz = 0; qz < Q1D; ++qz)
1313 {
1314 const real_t wz = Bc(qz,dz);
1315 const real_t wDz = Gc(qz,dz);
1316 for (int qy = 0; qy < Q1D; ++qy)
1317 {
1318 for (int qx = 0; qx < Q1D; ++qx)
1319 {
1320 mass[qz][qy][qx][0] += gradXY[qy][qx][0] * wz;
1321 mass[qz][qy][qx][1] += gradXY[qy][qx][1] * wz;
1322 mass[qz][qy][qx][2] += gradXY[qy][qx][2] * wDz;
1323 }
1324 }
1325 }
1326 }
1327
1328 // Apply D operator.
1329 for (int qz = 0; qz < Q1D; ++qz)
1330 {
1331 for (int qy = 0; qy < Q1D; ++qy)
1332 {
1333 for (int qx = 0; qx < Q1D; ++qx)
1334 {
1335 const real_t massX = mass[qz][qy][qx][0];
1336 const real_t massY = mass[qz][qy][qx][1];
1337 const real_t massZ = mass[qz][qy][qx][2];
1338 if (op_entries == 6)
1339 {
1340 const real_t O11 = op(qx,qy,qz,0,e);
1341 const real_t O12 = op(qx,qy,qz,1,e);
1342 const real_t O13 = op(qx,qy,qz,2,e);
1343 const real_t O22 = op(qx,qy,qz,3,e);
1344 const real_t O23 = op(qx,qy,qz,4,e);
1345 const real_t O33 = op(qx,qy,qz,5,e);
1346 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
1347 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
1348 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
1349 }
1350 else
1351 {
1352 // Non-symmetric operator stored row-major.
1353 const real_t O11 = op(qx,qy,qz,0,e);
1354 const real_t O12 = op(qx,qy,qz,1,e);
1355 const real_t O13 = op(qx,qy,qz,2,e);
1356 const real_t O21 = op(qx,qy,qz,3,e);
1357 const real_t O22 = op(qx,qy,qz,4,e);
1358 const real_t O23 = op(qx,qy,qz,5,e);
1359 const real_t O31 = op(qx,qy,qz,6,e);
1360 const real_t O32 = op(qx,qy,qz,7,e);
1361 const real_t O33 = op(qx,qy,qz,8,e);
1362 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
1363 mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
1364 mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
1365 }
1366 }
1367 }
1368 }
1369
1370 // Project to H(div) DOFs.
1371 for (int qz = 0; qz < Q1D; ++qz)
1372 {
1373 real_t massXY[MAX_D1D][MAX_D1D];
1374 int osc = 0;
1375
1376 for (int c = 0; c < VDIM; ++c)
1377 {
1378 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1379 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1380 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1381
1382 for (int dy = 0; dy < D1Dy; ++dy)
1383 {
1384 for (int dx = 0; dx < D1Dx; ++dx) { massXY[dy][dx] = 0.0; }
1385 }
1386
1387 for (int qy = 0; qy < Q1D; ++qy)
1388 {
1389 real_t massX[MAX_D1D];
1390 for (int dx = 0; dx < D1Dx; ++dx) { massX[dx] = 0.0; }
1391
1392 for (int qx = 0; qx < Q1D; ++qx)
1393 {
1394 for (int dx = 0; dx < D1Dx; ++dx)
1395 {
1396 massX[dx] += mass[qz][qy][qx][c] * ((c == 0) ? Bct(dx,qx) : Bot(dx,qx));
1397 }
1398 }
1399 for (int dy = 0; dy < D1Dy; ++dy)
1400 {
1401 const real_t wy = (c == 1) ? Bct(dy,qy) : Bot(dy,qy);
1402 for (int dx = 0; dx < D1Dx; ++dx) { massXY[dy][dx] += massX[dx] * wy; }
1403 }
1404 }
1405
1406 for (int dz = 0; dz < D1Dz; ++dz)
1407 {
1408 const real_t wz = (c == 2) ? Bct(dz,qz) : Bot(dz,qz);
1409 for (int dy = 0; dy < D1Dy; ++dy)
1410 {
1411 for (int dx = 0; dx < D1Dx; ++dx)
1412 {
1413 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += massXY[dy][dx] * wz;
1414 }
1415 }
1416 }
1417
1418 osc += D1Dx * D1Dy * D1Dz;
1419 }
1420 }
1421 });
1422}
1423
1424// Apply to x corresponding to DOFs in H(div), integrated
1425// against gradients of H^1 functions corresponding to y.
1426static void PAHdivH1ApplyTranspose3D(const int D1D,
1427 const int Q1D,
1428 const int NE,
1429 const Array<real_t> &bc,
1430 const Array<real_t> &bo,
1431 const Array<real_t> &bct,
1432 const Array<real_t> &gct,
1433 const int op_entries,
1434 const Vector &pa_data,
1435 const Vector &x,
1436 Vector &y)
1437{
1438 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
1439 "Error: D1D > MAX_D1D");
1440 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
1441 "Error: Q1D > MAX_Q1D");
1442
1443 constexpr static int VDIM = 3;
1444
1445 auto Bc = Reshape(bc.Read(), Q1D, D1D);
1446 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
1447 auto Bt = Reshape(bct.Read(), D1D, Q1D);
1448 auto Gt = Reshape(gct.Read(), D1D, Q1D);
1449 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, op_entries, NE);
1450 auto X = Reshape(x.Read(), 3*(D1D-1)*(D1D-1)*D1D, NE);
1451 auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, NE);
1452
1453 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1454 {
1455 constexpr static int MAX_D1D = DofQuadLimits::HDIV_MAX_D1D;
1456 constexpr static int MAX_Q1D = DofQuadLimits::HDIV_MAX_Q1D;
1457
1458 real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
1459 for (int qz = 0; qz < Q1D; ++qz)
1460 {
1461 for (int qy = 0; qy < Q1D; ++qy)
1462 {
1463 for (int qx = 0; qx < Q1D; ++qx)
1464 {
1465 for (int c = 0; c < VDIM; ++c) { mass[qz][qy][qx][c] = 0.0; }
1466 }
1467 }
1468 }
1469
1470 int osc = 0;
1471 for (int c = 0; c < VDIM; ++c)
1472 {
1473 const int D1Dx = (c == 0) ? D1D : D1D - 1;
1474 const int D1Dy = (c == 1) ? D1D : D1D - 1;
1475 const int D1Dz = (c == 2) ? D1D : D1D - 1;
1476
1477 for (int dz = 0; dz < D1Dz; ++dz)
1478 {
1479 real_t massXY[MAX_Q1D][MAX_Q1D];
1480 for (int qy = 0; qy < Q1D; ++qy)
1481 {
1482 for (int qx = 0; qx < Q1D; ++qx) { massXY[qy][qx] = 0.0; }
1483 }
1484
1485 for (int dy = 0; dy < D1Dy; ++dy)
1486 {
1487 real_t massX[MAX_Q1D];
1488 for (int qx = 0; qx < Q1D; ++qx) { massX[qx] = 0.0; }
1489
1490 for (int dx = 0; dx < D1Dx; ++dx)
1491 {
1492 const real_t t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
1493 for (int qx = 0; qx < Q1D; ++qx)
1494 {
1495 massX[qx] += t * ((c == 0) ? Bc(qx,dx) : Bo(qx,dx));
1496 }
1497 }
1498
1499 for (int qy = 0; qy < Q1D; ++qy)
1500 {
1501 const real_t wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
1502 for (int qx = 0; qx < Q1D; ++qx) { massXY[qy][qx] += massX[qx] * wy; }
1503 }
1504 }
1505
1506 for (int qz = 0; qz < Q1D; ++qz)
1507 {
1508 const real_t wz = (c == 2) ? Bc(qz,dz) : Bo(qz,dz);
1509 for (int qy = 0; qy < Q1D; ++qy)
1510 {
1511 for (int qx = 0; qx < Q1D; ++qx)
1512 {
1513 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
1514 }
1515 }
1516 }
1517 }
1518
1519 osc += D1Dx * D1Dy * D1Dz;
1520 }
1521
1522 // Apply D operator.
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 const real_t massX = mass[qz][qy][qx][0];
1530 const real_t massY = mass[qz][qy][qx][1];
1531 const real_t massZ = mass[qz][qy][qx][2];
1532 if (op_entries == 6)
1533 {
1534 const real_t O11 = op(qx,qy,qz,0,e);
1535 const real_t O12 = op(qx,qy,qz,1,e);
1536 const real_t O13 = op(qx,qy,qz,2,e);
1537 const real_t O22 = op(qx,qy,qz,3,e);
1538 const real_t O23 = op(qx,qy,qz,4,e);
1539 const real_t O33 = op(qx,qy,qz,5,e);
1540 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
1541 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
1542 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
1543 }
1544 else
1545 {
1546 // Non-symmetric operator stored row-major. For transpose, apply D^T.
1547 const real_t O11 = op(qx,qy,qz,0,e);
1548 const real_t O12 = op(qx,qy,qz,1,e);
1549 const real_t O13 = op(qx,qy,qz,2,e);
1550 const real_t O21 = op(qx,qy,qz,3,e);
1551 const real_t O22 = op(qx,qy,qz,4,e);
1552 const real_t O23 = op(qx,qy,qz,5,e);
1553 const real_t O31 = op(qx,qy,qz,6,e);
1554 const real_t O32 = op(qx,qy,qz,7,e);
1555 const real_t O33 = op(qx,qy,qz,8,e);
1556 mass[qz][qy][qx][0] = (O11*massX)+(O21*massY)+(O31*massZ);
1557 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O32*massZ);
1558 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
1559 }
1560 }
1561 }
1562 }
1563
1564 // Contract with test gradients for H1 output.
1565 for (int qz = 0; qz < Q1D; ++qz)
1566 {
1567 real_t gradXY[MAX_D1D][MAX_D1D][3];
1568 for (int dy = 0; dy < D1D; ++dy)
1569 {
1570 for (int dx = 0; dx < D1D; ++dx)
1571 {
1572 gradXY[dy][dx][0] = 0.0;
1573 gradXY[dy][dx][1] = 0.0;
1574 gradXY[dy][dx][2] = 0.0;
1575 }
1576 }
1577
1578 for (int qy = 0; qy < Q1D; ++qy)
1579 {
1580 real_t gradX[MAX_D1D][3];
1581 for (int dx = 0; dx < D1D; ++dx)
1582 {
1583 gradX[dx][0] = 0.0;
1584 gradX[dx][1] = 0.0;
1585 gradX[dx][2] = 0.0;
1586 }
1587 for (int qx = 0; qx < Q1D; ++qx)
1588 {
1589 const real_t gX = mass[qz][qy][qx][0];
1590 const real_t gY = mass[qz][qy][qx][1];
1591 const real_t gZ = mass[qz][qy][qx][2];
1592 for (int dx = 0; dx < D1D; ++dx)
1593 {
1594 const real_t wx = Bt(dx,qx);
1595 const real_t wDx = Gt(dx,qx);
1596 gradX[dx][0] += gX * wDx;
1597 gradX[dx][1] += gY * wx;
1598 gradX[dx][2] += gZ * wx;
1599 }
1600 }
1601 for (int dy = 0; dy < D1D; ++dy)
1602 {
1603 const real_t wy = Bt(dy,qy);
1604 const real_t wDy = Gt(dy,qy);
1605 for (int dx = 0; dx < D1D; ++dx)
1606 {
1607 gradXY[dy][dx][0] += gradX[dx][0] * wy;
1608 gradXY[dy][dx][1] += gradX[dx][1] * wDy;
1609 gradXY[dy][dx][2] += gradX[dx][2] * wy;
1610 }
1611 }
1612 }
1613
1614 for (int dz = 0; dz < D1D; ++dz)
1615 {
1616 const real_t wz = Bt(dz,qz);
1617 const real_t wDz = Gt(dz,qz);
1618 for (int dy = 0; dy < D1D; ++dy)
1619 {
1620 for (int dx = 0; dx < D1D; ++dx)
1621 {
1622 Y(dx,dy,dz,e) +=
1623 ((gradXY[dy][dx][0] * wz) +
1624 (gradXY[dy][dx][1] * wz) +
1625 (gradXY[dy][dx][2] * wDz));
1626 }
1627 }
1628 }
1629 }
1630 });
1631}
1632
1633void MixedVectorGradientIntegrator::AssemblePA(const FiniteElementSpace
1634 &trial_fes,
1635 const FiniteElementSpace &test_fes)
1636{
1637 // Assumes tensor-product elements, with a vector test space and H^1 trial space.
1638 Mesh *mesh = trial_fes.GetMesh();
1639 const FiniteElement *trial_fel = trial_fes.GetTypicalFE();
1640 const FiniteElement *test_fel = test_fes.GetTypicalFE();
1641
1642 const NodalTensorFiniteElement *trial_el =
1643 dynamic_cast<const NodalTensorFiniteElement*>(trial_fel);
1644 MFEM_VERIFY(trial_el != NULL, "Only NodalTensorFiniteElement is supported!");
1645 MFEM_VERIFY(trial_el->GetMapType() == FiniteElement::VALUE,
1646 "Only value map type is supported!");
1647
1648 const VectorTensorFiniteElement *test_el =
1649 dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
1650 MFEM_VERIFY(test_el != NULL, "Only VectorTensorFiniteElement is supported!");
1651
1652 const IntegrationRule *ir
1653 = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
1654 *mesh->GetTypicalElementTransformation());
1655 const int dims = trial_el->GetDim();
1656 MFEM_VERIFY(dims == 2 || dims == 3, "");
1657
1658 const int nq = ir->GetNPoints();
1659 dim = mesh->Dimension();
1660 MFEM_VERIFY(dim == 2 || dim == 3, "");
1661
1662 MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(), "");
1663
1664 ne = trial_fes.GetNE();
1665 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
1666 mapsC = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
1667 mapsO = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
1668 dofs1D = mapsC->ndof;
1669 quad1D = mapsC->nqpt;
1670 test_fetype = static_cast<FiniteElement::DerivType>(test_el->GetDerivType());
1671
1672 MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
1673
1674 QuadratureSpace qs(*mesh, *ir);
1675 CoefficientVector coeff(qs, CoefficientStorage::FULL);
1676 // NOTE: MFEM MatrixCoefficient values are stored in column-major ordering
1677 // (DenseMatrix layout). The PA diffusion setup kernels interpret matrix
1678 // entries in row-major order, i.e. they effectively see the transpose.
1679 // Projecting the transpose here ensures the kernels operate on the
1680 // intended matrix coefficient.
1681 if (MQ) { coeff.ProjectTranspose(*MQ); }
1682 else if (DQ) { coeff.Project(*DQ); }
1683 else if (Q) { coeff.Project(*Q); }
1684 else { coeff.SetConstant(1.0); }
1685
1686 const int coeffDim = coeff.GetVDim();
1687 op_entries = 0;
1688 if (test_fetype == mfem::FiniteElement::CURL)
1689 {
1690 op_entries = (dim == 2 ? (coeffDim == 4 ? 4 : 3) : (coeffDim == 9 ? 9 : 6));
1691 }
1692 else if (test_fetype == mfem::FiniteElement::DIV)
1693 {
1694 // MixedVectorGradient with an H(div) test space is generally non-symmetric,
1695 // so store the full operator.
1696 op_entries = dim * dim;
1697 }
1698 else
1699 {
1700 MFEM_ABORT("Unsupported test space derivative type.");
1701 }
1702
1703 pa_data.SetSize(op_entries * nq * ne, Device::GetMemoryType());
1704
1705 if (test_fetype == mfem::FiniteElement::CURL)
1706 {
1707 // Use the same setup functions as VectorFEMassIntegrator (H(curl)).
1708 if (dim == 3)
1709 {
1710 internal::PADiffusionSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J,
1711 coeff, pa_data);
1712 }
1713 else if (dim == 2)
1714 {
1715 internal::PADiffusionSetup2D<2>(quad1D, coeffDim, ne, ir->GetWeights(), geom->J,
1716 coeff, pa_data);
1717 }
1718 else
1719 {
1720 MFEM_ABORT("Unsupported dimension.");
1721 }
1722 }
1723 else if (test_fetype == mfem::FiniteElement::DIV)
1724 {
1725 if (dim == 3)
1726 {
1727 internal::PAMixedVectorGradientSetupHdiv3D(quad1D, coeffDim, ne,
1728 ir->GetWeights(),
1729 geom->J, coeff, pa_data);
1730 }
1731 else if (dim == 2)
1732 {
1733 internal::PAMixedVectorGradientSetupHdiv2D(quad1D, coeffDim, ne,
1734 ir->GetWeights(),
1735 geom->J, coeff, pa_data);
1736 }
1737 else
1738 {
1739 MFEM_ABORT("Unsupported dimension.");
1740 }
1741 }
1742}
1743
1744void MixedVectorGradientIntegrator::AddMultPA(const Vector &x, Vector &y) const
1745{
1746 if (test_fetype == mfem::FiniteElement::CURL)
1747 {
1748 if (dim == 3)
1749 {
1750 PAHcurlH1Apply3D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
1751 mapsO->Bt, mapsC->Bt, op_entries, pa_data, x, y);
1752 }
1753 else if (dim == 2)
1754 {
1755 PAHcurlH1Apply2D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
1756 mapsO->Bt, mapsC->Bt, op_entries, pa_data, x, y);
1757 }
1758 else
1759 {
1760 MFEM_ABORT("Unsupported dimension!");
1761 }
1762 }
1763 else if (test_fetype == mfem::FiniteElement::DIV)
1764 {
1765 if (dim == 3)
1766 {
1767 PAHdivH1Apply3D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
1768 mapsO->Bt, mapsC->Bt, op_entries, pa_data, x, y);
1769 }
1770 else if (dim == 2)
1771 {
1772 PAHdivH1Apply2D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
1773 mapsO->Bt, mapsC->Bt, op_entries, pa_data, x, y);
1774 }
1775 else
1776 {
1777 MFEM_ABORT("Unsupported dimension!");
1778 }
1779 }
1780 else
1781 {
1782 MFEM_ABORT("Unsupported test space derivative type!");
1783 }
1784}
1785
1786void MixedVectorGradientIntegrator::AddMultTransposePA(const Vector &x,
1787 Vector &y) const
1788{
1789 if (test_fetype == mfem::FiniteElement::CURL)
1790 {
1791 if (dim == 3)
1792 {
1793 PAHcurlH1ApplyTranspose3D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
1794 mapsC->Bt, mapsC->Gt, op_entries, pa_data, x, y);
1795 }
1796 else if (dim == 2)
1797 {
1798 PAHcurlH1ApplyTranspose2D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
1799 mapsC->Bt, mapsC->Gt, op_entries, pa_data, x, y);
1800 }
1801 else
1802 {
1803 MFEM_ABORT("Unsupported dimension!");
1804 }
1805 }
1806 else if (test_fetype == mfem::FiniteElement::DIV)
1807 {
1808 if (dim == 3)
1809 {
1810 PAHdivH1ApplyTranspose3D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
1811 mapsC->Bt, mapsC->Gt, op_entries, pa_data, x, y);
1812 }
1813 else if (dim == 2)
1814 {
1815 PAHdivH1ApplyTranspose2D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
1816 mapsC->Bt, mapsC->Gt, op_entries, pa_data, x, y);
1817 }
1818 else
1819 {
1820 MFEM_ABORT("Unsupported dimension!");
1821 }
1822 }
1823 else
1824 {
1825 MFEM_ABORT("Unsupported test space derivative type!");
1826 }
1827}
1828
1829} // namespace mfem
bilininteg_diffusion_kernels.hpp
mfem::CoefficientVector
Class to represent a coefficient evaluated at quadrature points.
Definition coefficient.hpp:2539
mfem::CoefficientVector::SetConstant
void SetConstant(real_t constant)
Set this vector to the given constant.
Definition coefficient.cpp:2142
mfem::CoefficientVector::Project
void Project(Coefficient &coeff)
Evaluate the given Coefficient at the quadrature points defined by qs.
Definition coefficient.cpp:2058
mfem::CoefficientVector::GetVDim
int GetVDim() const
Return the number of values per quadrature point.
Definition coefficient.cpp:2202
mfem::CoefficientVector::ProjectTranspose
void ProjectTranspose(MatrixCoefficient &coeff)
Project the transpose of coeff.
Definition coefficient.cpp:2126
mfem::Device::GetMemoryType
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:302
mfem::DofToQuad::G
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:222
mfem::DofToQuad::TENSOR
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
mfem::DofToQuad::B
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:201
mfem::DofToQuad::ndof
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:186
mfem::DofToQuad::Gt
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:229
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:190
mfem::DofToQuad::Bt
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:207
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:210
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:867
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:639
mfem::FiniteElementSpace::GetTypicalFE
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3896
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:294
mfem::FiniteElement::GetOrder
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:414
mfem::FiniteElement::GetDerivType
int GetDerivType() const
Returns the FiniteElement::DerivType of the element describing the spatial derivative method implemen...
Definition fe_base.hpp:441
mfem::FiniteElement::GetDim
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:381
mfem::FiniteElement::GetMapType
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:436
mfem::FiniteElement::DerivType
DerivType
Enumeration for DerivType: defines which derivative method is implemented.
Definition fe_base.hpp:363
mfem::FiniteElement::DIV
@ DIV
Implements CalcDivShape methods.
Definition fe_base.hpp:366
mfem::FiniteElement::CURL
@ CURL
Implements CalcCurlShape methods.
Definition fe_base.hpp:367
mfem::GeometricFactors::J
Vector J
Jacobians of the element transformations at all quadrature points.
Definition mesh.hpp:3158
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:96
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:255
mfem::IntegrationRule::GetWeights
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition intrules.cpp:98
mfem::Integrator::IntRule
const IntegrationRule * IntRule
Definition integrator.hpp:70
mfem::MassIntegrator::GetRule
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, const ElementTransformation &Trans, const bool stroud=false)
Definition bilininteg.cpp:1459
mfem::Mesh
Mesh data type.
Definition mesh.hpp:67
mfem::Mesh::Dimension
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1314
mfem::Mesh::GetTypicalElementTransformation
ElementTransformation * GetTypicalElementTransformation()
If the local mesh is not empty return GetElementTransformation(0); otherwise, return the identity tra...
Definition mesh.cpp:394
mfem::Mesh::GetGeometricFactors
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition mesh.cpp:958
mfem::MixedVectorGradientIntegrator::AddMultPA
void AddMultPA(const Vector &, Vector &) const override
Method for partially assembled action.
Definition bilininteg_mixedvecgrad_pa.cpp:1744
mfem::MixedVectorGradientIntegrator::AssemblePA
void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes) override
Definition bilininteg_mixedvecgrad_pa.cpp:1633
mfem::MixedVectorGradientIntegrator::AddMultTransposePA
void AddMultTransposePA(const Vector &, Vector &) const override
Method for partially assembled transposed action.
Definition bilininteg_mixedvecgrad_pa.cpp:1786
mfem::MixedVectorIntegrator::MQ
MatrixCoefficient * MQ
Definition bilininteg.hpp:674
mfem::MixedVectorIntegrator::DQ
DiagonalMatrixCoefficient * DQ
Definition bilininteg.hpp:673
mfem::QuadratureSpace
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:164
mfem::VectorTensorFiniteElement::GetDofToQuadOpen
const DofToQuad & GetDofToQuadOpen(const IntegrationRule &ir, DofToQuad::Mode mode) const
Definition fe_base.hpp:1442
mfem::VectorTensorFiniteElement::GetDofToQuad
const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const override
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.hpp:1434
mfem::Vector
Vector data type.
Definition vector.hpp:82
real_t
mfem::real_t real_t
Definition mtop_solvers.hpp:19
ds_common::mass
@ mass
Definition ds-common.hpp:46
mfem::Read
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition device.hpp:369
mfem::Reshape
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:138
mfem::CoefficientStorage::FULL
@ FULL
Store the coefficient as a full QuadratureFunction.
mfem::forall_2D
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:1220
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:1244
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:1134
mfem::DeviceDofQuadLimits::Get
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:138
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0