MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_mixedvecgrad_pa.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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
21// Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
22// integrated against H(curl) test functions corresponding to y.
23static void PAHcurlH1Apply2D(const int D1D,
24 const int Q1D,
25 const int NE,
26 const Array<real_t> &bc,
27 const Array<real_t> &gc,
28 const Array<real_t> &bot,
29 const Array<real_t> &bct,
30 const Vector &pa_data,
31 const Vector &x,
32 Vector &y)
33{
34 auto Bc = Reshape(bc.Read(), Q1D, D1D);
35 auto Gc = Reshape(gc.Read(), Q1D, D1D);
36 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
37 auto Bct = Reshape(bct.Read(), D1D, Q1D);
38 auto op = Reshape(pa_data.Read(), Q1D, Q1D, 3, NE);
39 auto X = Reshape(x.Read(), D1D, D1D, NE);
40 auto Y = Reshape(y.ReadWrite(), 2*(D1D-1)*D1D, NE);
41
42 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
43 {
44 constexpr static int VDIM = 2;
45 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
46 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
47
48 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
49
50 for (int qy = 0; qy < Q1D; ++qy)
51 {
52 for (int qx = 0; qx < Q1D; ++qx)
53 {
54 for (int c = 0; c < VDIM; ++c)
55 {
56 mass[qy][qx][c] = 0.0;
57 }
58 }
59 }
60
61 for (int dy = 0; dy < D1D; ++dy)
62 {
63 real_t gradX[MAX_Q1D][2];
64 for (int qx = 0; qx < Q1D; ++qx)
65 {
66 gradX[qx][0] = 0.0;
67 gradX[qx][1] = 0.0;
68 }
69 for (int dx = 0; dx < D1D; ++dx)
70 {
71 const real_t s = X(dx,dy,e);
72 for (int qx = 0; qx < Q1D; ++qx)
73 {
74 gradX[qx][0] += s * Bc(qx,dx);
75 gradX[qx][1] += s * Gc(qx,dx);
76 }
77 }
78 for (int qy = 0; qy < Q1D; ++qy)
79 {
80 const real_t wy = Bc(qy,dy);
81 const real_t wDy = Gc(qy,dy);
82 for (int qx = 0; qx < Q1D; ++qx)
83 {
84 const real_t wx = gradX[qx][0];
85 const real_t wDx = gradX[qx][1];
86 mass[qy][qx][0] += wDx * wy;
87 mass[qy][qx][1] += wx * wDy;
88 }
89 }
90 }
91
92 // Apply D operator.
93 for (int qy = 0; qy < Q1D; ++qy)
94 {
95 for (int qx = 0; qx < Q1D; ++qx)
96 {
97 const real_t O11 = op(qx,qy,0,e);
98 const real_t O12 = op(qx,qy,1,e);
99 const real_t O22 = op(qx,qy,2,e);
100 const real_t massX = mass[qy][qx][0];
101 const real_t massY = mass[qy][qx][1];
102 mass[qy][qx][0] = (O11*massX)+(O12*massY);
103 mass[qy][qx][1] = (O12*massX)+(O22*massY);
104 }
105 }
106
107 for (int qy = 0; qy < Q1D; ++qy)
108 {
109 int osc = 0;
110
111 for (int c = 0; c < VDIM; ++c) // loop over x, y components
112 {
113 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
114 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
115
116 real_t massX[MAX_D1D];
117 for (int dx = 0; dx < D1Dx; ++dx)
118 {
119 massX[dx] = 0;
120 }
121 for (int qx = 0; qx < Q1D; ++qx)
122 {
123 for (int dx = 0; dx < D1Dx; ++dx)
124 {
125 massX[dx] += mass[qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
126 }
127 }
128
129 for (int dy = 0; dy < D1Dy; ++dy)
130 {
131 const real_t wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
132
133 for (int dx = 0; dx < D1Dx; ++dx)
134 {
135 Y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
136 }
137 }
138
139 osc += D1Dx * D1Dy;
140 } // loop c
141 }
142 }); // end of element loop
143}
144
145// Apply to x corresponding to DOFs in H(curl), integrated
146// against gradients of H^1 functions corresponding to y.
147static void PAHcurlH1ApplyTranspose2D(const int D1D,
148 const int Q1D,
149 const int NE,
150 const Array<real_t> &bc,
151 const Array<real_t> &bo,
152 const Array<real_t> &bct,
153 const Array<real_t> &gct,
154 const Vector &pa_data,
155 const Vector &x,
156 Vector &y)
157{
158 auto Bc = Reshape(bc.Read(), Q1D, D1D);
159 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
160 auto Bt = Reshape(bct.Read(), D1D, Q1D);
161 auto Gt = Reshape(gct.Read(), D1D, Q1D);
162 auto op = Reshape(pa_data.Read(), Q1D, Q1D, 3, NE);
163 auto X = Reshape(x.Read(), 2*(D1D-1)*D1D, NE);
164 auto Y = Reshape(y.ReadWrite(), D1D, D1D, NE);
165
166 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
167 {
168 constexpr static int VDIM = 2;
169 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
170 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
171
172 real_t mass[MAX_Q1D][MAX_Q1D][VDIM];
173
174 for (int qy = 0; qy < Q1D; ++qy)
175 {
176 for (int qx = 0; qx < Q1D; ++qx)
177 {
178 for (int c = 0; c < VDIM; ++c)
179 {
180 mass[qy][qx][c] = 0.0;
181 }
182 }
183 }
184
185 int osc = 0;
186
187 for (int c = 0; c < VDIM; ++c) // loop over x, y components
188 {
189 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
190 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
191
192 for (int dy = 0; dy < D1Dy; ++dy)
193 {
194 real_t massX[MAX_Q1D];
195 for (int qx = 0; qx < Q1D; ++qx)
196 {
197 massX[qx] = 0.0;
198 }
199
200 for (int dx = 0; dx < D1Dx; ++dx)
201 {
202 const real_t t = X(dx + (dy * D1Dx) + osc, e);
203 for (int qx = 0; qx < Q1D; ++qx)
204 {
205 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
206 }
207 }
208
209 for (int qy = 0; qy < Q1D; ++qy)
210 {
211 const real_t wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
212 for (int qx = 0; qx < Q1D; ++qx)
213 {
214 mass[qy][qx][c] += massX[qx] * wy;
215 }
216 }
217 }
218
219 osc += D1Dx * D1Dy;
220 } // loop (c) over components
221
222 // Apply D operator.
223 for (int qy = 0; qy < Q1D; ++qy)
224 {
225 for (int qx = 0; qx < Q1D; ++qx)
226 {
227 const real_t O11 = op(qx,qy,0,e);
228 const real_t O12 = op(qx,qy,1,e);
229 const real_t O22 = op(qx,qy,2,e);
230 const real_t massX = mass[qy][qx][0];
231 const real_t massY = mass[qy][qx][1];
232 mass[qy][qx][0] = (O11*massX)+(O12*massY);
233 mass[qy][qx][1] = (O12*massX)+(O22*massY);
234 }
235 }
236
237 for (int qy = 0; qy < Q1D; ++qy)
238 {
239 real_t gradX[MAX_D1D][2];
240 for (int dx = 0; dx < D1D; ++dx)
241 {
242 gradX[dx][0] = 0;
243 gradX[dx][1] = 0;
244 }
245 for (int qx = 0; qx < Q1D; ++qx)
246 {
247 const real_t gX = mass[qy][qx][0];
248 const real_t gY = mass[qy][qx][1];
249 for (int dx = 0; dx < D1D; ++dx)
250 {
251 const real_t wx = Bt(dx,qx);
252 const real_t wDx = Gt(dx,qx);
253 gradX[dx][0] += gX * wDx;
254 gradX[dx][1] += gY * wx;
255 }
256 }
257 for (int dy = 0; dy < D1D; ++dy)
258 {
259 const real_t wy = Bt(dy,qy);
260 const real_t wDy = Gt(dy,qy);
261 for (int dx = 0; dx < D1D; ++dx)
262 {
263 Y(dx,dy,e) += ((gradX[dx][0] * wy) + (gradX[dx][1] * wDy));
264 }
265 }
266 }
267 }); // end of element loop
268}
269
270// Apply to x corresponding to DOFs in H^1 (trial), whose gradients are
271// integrated against H(curl) test functions corresponding to y.
272static void PAHcurlH1Apply3D(const int D1D,
273 const int Q1D,
274 const int NE,
275 const Array<real_t> &bc,
276 const Array<real_t> &gc,
277 const Array<real_t> &bot,
278 const Array<real_t> &bct,
279 const Vector &pa_data,
280 const Vector &x,
281 Vector &y)
282{
283 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
284 "Error: D1D > MAX_D1D");
285 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
286 "Error: Q1D > MAX_Q1D");
287
288 constexpr static int VDIM = 3;
289
290 auto Bc = Reshape(bc.Read(), Q1D, D1D);
291 auto Gc = Reshape(gc.Read(), Q1D, D1D);
292 auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
293 auto Bct = Reshape(bct.Read(), D1D, Q1D);
294 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
295 auto X = Reshape(x.Read(), D1D, D1D, D1D, NE);
296 auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
297
298 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
299 {
300 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
301 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
302
303 real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
304
305 for (int qz = 0; qz < Q1D; ++qz)
306 {
307 for (int qy = 0; qy < Q1D; ++qy)
308 {
309 for (int qx = 0; qx < Q1D; ++qx)
310 {
311 for (int c = 0; c < VDIM; ++c)
312 {
313 mass[qz][qy][qx][c] = 0.0;
314 }
315 }
316 }
317 }
318
319 for (int dz = 0; dz < D1D; ++dz)
320 {
321 real_t gradXY[MAX_Q1D][MAX_Q1D][3];
322 for (int qy = 0; qy < Q1D; ++qy)
323 {
324 for (int qx = 0; qx < Q1D; ++qx)
325 {
326 gradXY[qy][qx][0] = 0.0;
327 gradXY[qy][qx][1] = 0.0;
328 gradXY[qy][qx][2] = 0.0;
329 }
330 }
331 for (int dy = 0; dy < D1D; ++dy)
332 {
333 real_t gradX[MAX_Q1D][2];
334 for (int qx = 0; qx < Q1D; ++qx)
335 {
336 gradX[qx][0] = 0.0;
337 gradX[qx][1] = 0.0;
338 }
339 for (int dx = 0; dx < D1D; ++dx)
340 {
341 const real_t s = X(dx,dy,dz,e);
342 for (int qx = 0; qx < Q1D; ++qx)
343 {
344 gradX[qx][0] += s * Bc(qx,dx);
345 gradX[qx][1] += s * Gc(qx,dx);
346 }
347 }
348 for (int qy = 0; qy < Q1D; ++qy)
349 {
350 const real_t wy = Bc(qy,dy);
351 const real_t wDy = Gc(qy,dy);
352 for (int qx = 0; qx < Q1D; ++qx)
353 {
354 const real_t wx = gradX[qx][0];
355 const real_t wDx = gradX[qx][1];
356 gradXY[qy][qx][0] += wDx * wy;
357 gradXY[qy][qx][1] += wx * wDy;
358 gradXY[qy][qx][2] += wx * wy;
359 }
360 }
361 }
362 for (int qz = 0; qz < Q1D; ++qz)
363 {
364 const real_t wz = Bc(qz,dz);
365 const real_t wDz = Gc(qz,dz);
366 for (int qy = 0; qy < Q1D; ++qy)
367 {
368 for (int qx = 0; qx < Q1D; ++qx)
369 {
370 mass[qz][qy][qx][0] += gradXY[qy][qx][0] * wz;
371 mass[qz][qy][qx][1] += gradXY[qy][qx][1] * wz;
372 mass[qz][qy][qx][2] += gradXY[qy][qx][2] * wDz;
373 }
374 }
375 }
376 }
377
378 // Apply D operator.
379 for (int qz = 0; qz < Q1D; ++qz)
380 {
381 for (int qy = 0; qy < Q1D; ++qy)
382 {
383 for (int qx = 0; qx < Q1D; ++qx)
384 {
385 const real_t O11 = op(qx,qy,qz,0,e);
386 const real_t O12 = op(qx,qy,qz,1,e);
387 const real_t O13 = op(qx,qy,qz,2,e);
388 const real_t O22 = op(qx,qy,qz,3,e);
389 const real_t O23 = op(qx,qy,qz,4,e);
390 const real_t O33 = op(qx,qy,qz,5,e);
391 const real_t massX = mass[qz][qy][qx][0];
392 const real_t massY = mass[qz][qy][qx][1];
393 const real_t massZ = mass[qz][qy][qx][2];
394 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
395 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
396 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
397 }
398 }
399 }
400
401 for (int qz = 0; qz < Q1D; ++qz)
402 {
403 real_t massXY[MAX_D1D][MAX_D1D];
404
405 int osc = 0;
406
407 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
408 {
409 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
410 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
411 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
412
413 for (int dy = 0; dy < D1Dy; ++dy)
414 {
415 for (int dx = 0; dx < D1Dx; ++dx)
416 {
417 massXY[dy][dx] = 0.0;
418 }
419 }
420 for (int qy = 0; qy < Q1D; ++qy)
421 {
422 real_t massX[MAX_D1D];
423 for (int dx = 0; dx < D1Dx; ++dx)
424 {
425 massX[dx] = 0;
426 }
427 for (int qx = 0; qx < Q1D; ++qx)
428 {
429 for (int dx = 0; dx < D1Dx; ++dx)
430 {
431 massX[dx] += mass[qz][qy][qx][c] * ((c == 0) ? Bot(dx,qx) : Bct(dx,qx));
432 }
433 }
434 for (int dy = 0; dy < D1Dy; ++dy)
435 {
436 const real_t wy = (c == 1) ? Bot(dy,qy) : Bct(dy,qy);
437 for (int dx = 0; dx < D1Dx; ++dx)
438 {
439 massXY[dy][dx] += massX[dx] * wy;
440 }
441 }
442 }
443
444 for (int dz = 0; dz < D1Dz; ++dz)
445 {
446 const real_t wz = (c == 2) ? Bot(dz,qz) : Bct(dz,qz);
447 for (int dy = 0; dy < D1Dy; ++dy)
448 {
449 for (int dx = 0; dx < D1Dx; ++dx)
450 {
451 Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += massXY[dy][dx] * wz;
452 }
453 }
454 }
455
456 osc += D1Dx * D1Dy * D1Dz;
457 } // loop c
458 } // loop qz
459 }); // end of element loop
460}
461
462// Apply to x corresponding to DOFs in H(curl), integrated
463// against gradients of H^1 functions corresponding to y.
464static void PAHcurlH1ApplyTranspose3D(const int D1D,
465 const int Q1D,
466 const int NE,
467 const Array<real_t> &bc,
468 const Array<real_t> &bo,
469 const Array<real_t> &bct,
470 const Array<real_t> &gct,
471 const Vector &pa_data,
472 const Vector &x,
473 Vector &y)
474{
475 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
476 "Error: D1D > MAX_D1D");
477 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
478 "Error: Q1D > MAX_Q1D");
479
480 constexpr static int VDIM = 3;
481
482 auto Bc = Reshape(bc.Read(), Q1D, D1D);
483 auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
484 auto Bt = Reshape(bct.Read(), D1D, Q1D);
485 auto Gt = Reshape(gct.Read(), D1D, Q1D);
486 auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, 6, NE);
487 auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
488 auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, NE);
489
490 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
491 {
492 constexpr static int MAX_D1D = DofQuadLimits::HCURL_MAX_D1D;
493 constexpr static int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
494
495 real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
496
497 for (int qz = 0; qz < Q1D; ++qz)
498 {
499 for (int qy = 0; qy < Q1D; ++qy)
500 {
501 for (int qx = 0; qx < Q1D; ++qx)
502 {
503 for (int c = 0; c < VDIM; ++c)
504 {
505 mass[qz][qy][qx][c] = 0.0;
506 }
507 }
508 }
509 }
510
511 int osc = 0;
512
513 for (int c = 0; c < VDIM; ++c) // loop over x, y, z components
514 {
515 const int D1Dz = (c == 2) ? D1D - 1 : D1D;
516 const int D1Dy = (c == 1) ? D1D - 1 : D1D;
517 const int D1Dx = (c == 0) ? D1D - 1 : D1D;
518
519 for (int dz = 0; dz < D1Dz; ++dz)
520 {
521 real_t massXY[MAX_Q1D][MAX_Q1D];
522 for (int qy = 0; qy < Q1D; ++qy)
523 {
524 for (int qx = 0; qx < Q1D; ++qx)
525 {
526 massXY[qy][qx] = 0.0;
527 }
528 }
529
530 for (int dy = 0; dy < D1Dy; ++dy)
531 {
532 real_t massX[MAX_Q1D];
533 for (int qx = 0; qx < Q1D; ++qx)
534 {
535 massX[qx] = 0.0;
536 }
537
538 for (int dx = 0; dx < D1Dx; ++dx)
539 {
540 const real_t t = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
541 for (int qx = 0; qx < Q1D; ++qx)
542 {
543 massX[qx] += t * ((c == 0) ? Bo(qx,dx) : Bc(qx,dx));
544 }
545 }
546
547 for (int qy = 0; qy < Q1D; ++qy)
548 {
549 const real_t wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
550 for (int qx = 0; qx < Q1D; ++qx)
551 {
552 const real_t wx = massX[qx];
553 massXY[qy][qx] += wx * wy;
554 }
555 }
556 }
557
558 for (int qz = 0; qz < Q1D; ++qz)
559 {
560 const real_t wz = (c == 2) ? Bo(qz,dz) : Bc(qz,dz);
561 for (int qy = 0; qy < Q1D; ++qy)
562 {
563 for (int qx = 0; qx < Q1D; ++qx)
564 {
565 mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
566 }
567 }
568 }
569 }
570
571 osc += D1Dx * D1Dy * D1Dz;
572 } // loop (c) over components
573
574 // Apply D operator.
575 for (int qz = 0; qz < Q1D; ++qz)
576 {
577 for (int qy = 0; qy < Q1D; ++qy)
578 {
579 for (int qx = 0; qx < Q1D; ++qx)
580 {
581 const real_t O11 = op(qx,qy,qz,0,e);
582 const real_t O12 = op(qx,qy,qz,1,e);
583 const real_t O13 = op(qx,qy,qz,2,e);
584 const real_t O22 = op(qx,qy,qz,3,e);
585 const real_t O23 = op(qx,qy,qz,4,e);
586 const real_t O33 = op(qx,qy,qz,5,e);
587 const real_t massX = mass[qz][qy][qx][0];
588 const real_t massY = mass[qz][qy][qx][1];
589 const real_t massZ = mass[qz][qy][qx][2];
590 mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
591 mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
592 mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
593 }
594 }
595 }
596
597 for (int qz = 0; qz < Q1D; ++qz)
598 {
599 real_t gradXY[MAX_D1D][MAX_D1D][3];
600 for (int dy = 0; dy < D1D; ++dy)
601 {
602 for (int dx = 0; dx < D1D; ++dx)
603 {
604 gradXY[dy][dx][0] = 0;
605 gradXY[dy][dx][1] = 0;
606 gradXY[dy][dx][2] = 0;
607 }
608 }
609 for (int qy = 0; qy < Q1D; ++qy)
610 {
611 real_t gradX[MAX_D1D][3];
612 for (int dx = 0; dx < D1D; ++dx)
613 {
614 gradX[dx][0] = 0;
615 gradX[dx][1] = 0;
616 gradX[dx][2] = 0;
617 }
618 for (int qx = 0; qx < Q1D; ++qx)
619 {
620 const real_t gX = mass[qz][qy][qx][0];
621 const real_t gY = mass[qz][qy][qx][1];
622 const real_t gZ = mass[qz][qy][qx][2];
623 for (int dx = 0; dx < D1D; ++dx)
624 {
625 const real_t wx = Bt(dx,qx);
626 const real_t wDx = Gt(dx,qx);
627 gradX[dx][0] += gX * wDx;
628 gradX[dx][1] += gY * wx;
629 gradX[dx][2] += gZ * wx;
630 }
631 }
632 for (int dy = 0; dy < D1D; ++dy)
633 {
634 const real_t wy = Bt(dy,qy);
635 const real_t wDy = Gt(dy,qy);
636 for (int dx = 0; dx < D1D; ++dx)
637 {
638 gradXY[dy][dx][0] += gradX[dx][0] * wy;
639 gradXY[dy][dx][1] += gradX[dx][1] * wDy;
640 gradXY[dy][dx][2] += gradX[dx][2] * wy;
641 }
642 }
643 }
644 for (int dz = 0; dz < D1D; ++dz)
645 {
646 const real_t wz = Bt(dz,qz);
647 const real_t wDz = Gt(dz,qz);
648 for (int dy = 0; dy < D1D; ++dy)
649 {
650 for (int dx = 0; dx < D1D; ++dx)
651 {
652 Y(dx,dy,dz,e) +=
653 ((gradXY[dy][dx][0] * wz) +
654 (gradXY[dy][dx][1] * wz) +
655 (gradXY[dy][dx][2] * wDz));
656 }
657 }
658 }
659 } // loop qz
660 }); // end of element loop
661}
662
663void MixedVectorGradientIntegrator::AssemblePA(const FiniteElementSpace
664 &trial_fes,
665 const FiniteElementSpace &test_fes)
666{
667 // Assumes tensor-product elements, with a vector test space and H^1 trial space.
668 Mesh *mesh = trial_fes.GetMesh();
669 const FiniteElement *trial_fel = trial_fes.GetTypicalFE();
670 const FiniteElement *test_fel = test_fes.GetTypicalFE();
671
672 const NodalTensorFiniteElement *trial_el =
673 dynamic_cast<const NodalTensorFiniteElement*>(trial_fel);
674 MFEM_VERIFY(trial_el != NULL, "Only NodalTensorFiniteElement is supported!");
675
676 const VectorTensorFiniteElement *test_el =
677 dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
678 MFEM_VERIFY(test_el != NULL, "Only VectorTensorFiniteElement is supported!");
679
680 const IntegrationRule *ir
681 = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
682 *mesh->GetTypicalElementTransformation());
683 const int dims = trial_el->GetDim();
684 MFEM_VERIFY(dims == 2 || dims == 3, "");
685
686 const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
687 const int nq = ir->GetNPoints();
688 dim = mesh->Dimension();
689 MFEM_VERIFY(dim == 2 || dim == 3, "");
690
691 MFEM_VERIFY(trial_el->GetOrder() == test_el->GetOrder(), "");
692
693 ne = trial_fes.GetNE();
694 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
695 mapsC = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
696 mapsO = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
697 dofs1D = mapsC->ndof;
698 quad1D = mapsC->nqpt;
699
700 MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
701
702 pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
703
704 QuadratureSpace qs(*mesh, *ir);
705 CoefficientVector coeff(Q, qs, CoefficientStorage::FULL);
706
707 // Use the same setup functions as VectorFEMassIntegrator.
708 if (test_el->GetDerivType() == mfem::FiniteElement::CURL && dim == 3)
709 {
710 internal::PADiffusionSetup3D(quad1D, 1, ne, ir->GetWeights(), geom->J,
711 coeff, pa_data);
712 }
713 else if (test_el->GetDerivType() == mfem::FiniteElement::CURL && dim == 2)
714 {
715 internal::PADiffusionSetup2D<2>(quad1D, 1, ne, ir->GetWeights(), geom->J,
716 coeff, pa_data);
717 }
718 else
719 {
720 MFEM_ABORT("Unknown kernel.");
721 }
722}
723
724void MixedVectorGradientIntegrator::AddMultPA(const Vector &x, Vector &y) const
725{
726 if (dim == 3)
727 {
728 PAHcurlH1Apply3D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
729 mapsO->Bt, mapsC->Bt, pa_data, x, y);
730 }
731 else if (dim == 2)
732 {
733 PAHcurlH1Apply2D(dofs1D, quad1D, ne, mapsC->B, mapsC->G,
734 mapsO->Bt, mapsC->Bt, pa_data, x, y);
735 }
736 else
737 {
738 MFEM_ABORT("Unsupported dimension!");
739 }
740}
741
742void MixedVectorGradientIntegrator::AddMultTransposePA(const Vector &x,
743 Vector &y) const
744{
745 if (dim == 3)
746 {
747 PAHcurlH1ApplyTranspose3D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
748 mapsC->Bt, mapsC->Gt, pa_data, x, y);
749 }
750 else if (dim == 2)
751 {
752 PAHcurlH1ApplyTranspose2D(dofs1D, quad1D, ne, mapsC->B, mapsO->B,
753 mapsC->Bt, mapsC->Gt, pa_data, x, y);
754 }
755 else
756 {
757 MFEM_ABORT("Unsupported dimension!");
758 }
759}
760
761} // namespace mfem
bilininteg_diffusion_kernels.hpp
mfem::CoefficientVector
Class to represent a coefficient evaluated at quadrature points.
Definition coefficient.hpp:2429
mfem::Device::GetMemoryType
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:278
mfem::DofToQuad::G
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:214
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:193
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:178
mfem::DofToQuad::Gt
Array< real_t > Gt
Transpose of G.
Definition fe_base.hpp:221
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
mfem::DofToQuad::Bt
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:199
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:891
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
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:3871
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:244
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:338
mfem::FiniteElement::GetDerivType
int GetDerivType() const
Returns the FiniteElement::DerivType of the element describing the spatial derivative method implemen...
Definition fe_base.hpp:365
mfem::FiniteElement::GetDim
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:321
mfem::FiniteElement::CURL
@ CURL
Implements CalcCurlShape methods.
Definition fe_base.hpp:307
mfem::GeometricFactors::J
Vector J
Jacobians of the element transformations at all quadrature points.
Definition mesh.hpp:2976
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
mfem::IntegrationRule::GetWeights
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition intrules.cpp:86
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)
Definition bilininteg.cpp:1421
mfem::Mesh
Mesh data type.
Definition mesh.hpp:64
mfem::Mesh::Dimension
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
mfem::Mesh::GetTypicalElementTransformation
ElementTransformation * GetTypicalElementTransformation()
If the local mesh is not empty return GetElementTransformation(0); otherwise, return the identity tra...
Definition mesh.cpp:390
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:880
mfem::MixedVectorGradientIntegrator::AddMultPA
void AddMultPA(const Vector &, Vector &) const override
Method for partially assembled action.
Definition bilininteg_mixedvecgrad_pa.cpp:724
mfem::MixedVectorGradientIntegrator::AssemblePA
void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes) override
Definition bilininteg_mixedvecgrad_pa.cpp:663
mfem::MixedVectorGradientIntegrator::AddMultTransposePA
void AddMultTransposePA(const Vector &, Vector &) const override
Method for partially assembled transposed action.
Definition bilininteg_mixedvecgrad_pa.cpp:742
mfem::QuadratureSpace
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:120
mfem::VectorTensorFiniteElement::GetDofToQuadOpen
const DofToQuad & GetDofToQuadOpen(const IntegrationRule &ir, DofToQuad::Mode mode) const
Definition fe_base.hpp:1362
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:1354
mfem::Vector
Vector data type.
Definition vector.hpp:82
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:131
mfem::CoefficientStorage::FULL
@ FULL
Store the coefficient as a full QuadratureFunction.
mfem::real_t
float real_t
Definition config.hpp:43
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:753
mfem::DeviceDofQuadLimits::Get
static const DeviceDofQuadLimits & Get()
Return a const reference to the DeviceDofQuadLimits singleton.
Definition forall.hpp:128
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0