MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_vecdiffusion_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 "../ceed/integrators/diffusion/diffusion.hpp"
17
18namespace mfem
19{
20
21// PA Diffusion Assemble 2D kernel
22static void PAVectorDiffusionSetup2D(const int Q1D,
23 const int NE,
24 const Array<real_t> &w,
25 const Vector &j,
26 const Vector &c,
27 Vector &op)
28{
29 const int NQ = Q1D*Q1D;
30 auto W = w.Read();
31
32 auto J = Reshape(j.Read(), NQ, 2, 2, NE);
33 auto y = Reshape(op.Write(), NQ, 3, NE);
34
35 const bool const_c = c.Size() == 1;
36 const auto C = const_c ? Reshape(c.Read(), 1,1) :
37 Reshape(c.Read(), NQ, NE);
38
39
40 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
41 {
42 for (int q = 0; q < NQ; ++q)
43 {
44 const real_t J11 = J(q,0,0,e);
45 const real_t J21 = J(q,1,0,e);
46 const real_t J12 = J(q,0,1,e);
47 const real_t J22 = J(q,1,1,e);
48
49 const real_t C1 = const_c ? C(0,0) : C(q,e);
50 const real_t c_detJ = W[q] * C1 / ((J11*J22)-(J21*J12));
51 y(q,0,e) = c_detJ * (J12*J12 + J22*J22); // 1,1
52 y(q,1,e) = -c_detJ * (J12*J11 + J22*J21); // 1,2
53 y(q,2,e) = c_detJ * (J11*J11 + J21*J21); // 2,2
54 }
55 });
56}
57
58// PA Diffusion Assemble 3D kernel
59static void PAVectorDiffusionSetup3D(const int Q1D,
60 const int NE,
61 const Array<real_t> &w,
62 const Vector &j,
63 const Vector &c,
64 Vector &op)
65{
66 const int NQ = Q1D*Q1D*Q1D;
67 auto W = w.Read();
68 auto J = Reshape(j.Read(), NQ, 3, 3, NE);
69 auto y = Reshape(op.Write(), NQ, 6, NE);
70
71 const bool const_c = c.Size() == 1;
72 const auto C = const_c ? Reshape(c.Read(), 1,1) :
73 Reshape(c.Read(), NQ,NE);
74
75
76 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
77 {
78 for (int q = 0; q < NQ; ++q)
79 {
80 const real_t J11 = J(q,0,0,e);
81 const real_t J21 = J(q,1,0,e);
82 const real_t J31 = J(q,2,0,e);
83 const real_t J12 = J(q,0,1,e);
84 const real_t J22 = J(q,1,1,e);
85 const real_t J32 = J(q,2,1,e);
86 const real_t J13 = J(q,0,2,e);
87 const real_t J23 = J(q,1,2,e);
88 const real_t J33 = J(q,2,2,e);
89 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
90 J21 * (J12 * J33 - J32 * J13) +
91 J31 * (J12 * J23 - J22 * J13);
92
93 const real_t C1 = const_c ? C(0,0) : C(q,e);
94
95 const real_t c_detJ = W[q] * C1 / detJ;
96 // adj(J)
97 const real_t A11 = (J22 * J33) - (J23 * J32);
98 const real_t A12 = (J32 * J13) - (J12 * J33);
99 const real_t A13 = (J12 * J23) - (J22 * J13);
100 const real_t A21 = (J31 * J23) - (J21 * J33);
101 const real_t A22 = (J11 * J33) - (J13 * J31);
102 const real_t A23 = (J21 * J13) - (J11 * J23);
103 const real_t A31 = (J21 * J32) - (J31 * J22);
104 const real_t A32 = (J31 * J12) - (J11 * J32);
105 const real_t A33 = (J11 * J22) - (J12 * J21);
106 // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
107 y(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
108 y(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
109 y(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
110 y(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
111 y(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
112 y(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
113 }
114 });
115}
116
117static void PAVectorDiffusionSetup(const int dim,
118 const int Q1D,
119 const int NE,
120 const Array<real_t> &W,
121 const Vector &J,
122 const Vector &C,
123 Vector &op)
124{
125 if (!(dim == 2 || dim == 3))
126 {
127 MFEM_ABORT("Dimension not supported.");
128 }
129 if (dim == 2)
130 {
131 PAVectorDiffusionSetup2D(Q1D, NE, W, J, C, op);
132 }
133 if (dim == 3)
134 {
135 PAVectorDiffusionSetup3D(Q1D, NE, W, J, C, op);
136 }
137}
138
139void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
140{
141 // Assumes tensor-product elements
142 Mesh *mesh = fes.GetMesh();
143 const FiniteElement &el = *fes.GetTypicalFE();
144 const IntegrationRule *ir
145 = IntRule ? IntRule : &DiffusionIntegrator::GetRule(el, el);
146 if (DeviceCanUseCeed())
147 {
148 delete ceedOp;
149 const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
150 fes.IsVariableOrder();
151 if (mixed)
152 {
153 ceedOp = new ceed::MixedPADiffusionIntegrator(*this, fes, Q);
154 }
155 else
156 {
157 ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
158 }
159 return;
160 }
161 const int dims = el.GetDim();
162 const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
163 const int nq = ir->GetNPoints();
164 dim = mesh->Dimension();
165 sdim = mesh->SpaceDimension();
166 ne = fes.GetNE();
167 geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
168 maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
169 dofs1D = maps->ndof;
170 quad1D = maps->nqpt;
171 pa_data.SetSize(symmDims * nq * ne, Device::GetDeviceMemoryType());
172
173 MFEM_VERIFY(!VQ && !MQ,
174 "Only scalar coefficient supported for partial assembly for VectorDiffusionIntegrator");
175
176 QuadratureSpace qs(*mesh, *ir);
177 CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);
178
179 const Array<real_t> &w = ir->GetWeights();
180 const Vector &j = geom->J;
181 Vector &d = pa_data;
182 if (dim == 1) { MFEM_ABORT("dim==1 not supported in PAVectorDiffusionSetup"); }
183 if (dim == 2 && sdim == 3)
184 {
185 constexpr int DIM = 2;
186 constexpr int SDIM = 3;
187 const int NQ = quad1D*quad1D;
188 auto W = w.Read();
189 auto J = Reshape(j.Read(), NQ, SDIM, DIM, ne);
190 auto D = Reshape(d.Write(), NQ, SDIM, ne);
191
192 const bool const_c = coeff.Size() == 1;
193 const auto C = const_c ? Reshape(coeff.Read(), 1,1) :
194 Reshape(coeff.Read(), NQ,ne);
195
196 mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
197 {
198 for (int q = 0; q < NQ; ++q)
199 {
200 const real_t wq = W[q];
201 const real_t J11 = J(q,0,0,e);
202 const real_t J21 = J(q,1,0,e);
203 const real_t J31 = J(q,2,0,e);
204 const real_t J12 = J(q,0,1,e);
205 const real_t J22 = J(q,1,1,e);
206 const real_t J32 = J(q,2,1,e);
207 const real_t E = J11*J11 + J21*J21 + J31*J31;
208 const real_t G = J12*J12 + J22*J22 + J32*J32;
209 const real_t F = J11*J12 + J21*J22 + J31*J32;
210 const real_t iw = 1.0 / sqrt(E*G - F*F);
211 const real_t C1 = const_c ? C(0,0) : C(q,e);
212 const real_t alpha = wq * C1 * iw;
213 D(q,0,e) = alpha * G; // 1,1
214 D(q,1,e) = -alpha * F; // 1,2
215 D(q,2,e) = alpha * E; // 2,2
216 }
217 });
218 }
219 else
220 {
221 PAVectorDiffusionSetup(dim, quad1D, ne, w, j, coeff, d);
222 }
223}
224
225template<int T_D1D = 0, int T_Q1D = 0>
226static void PAVectorDiffusionDiagonal2D(const int NE,
227 const Array<real_t> &b,
228 const Array<real_t> &g,
229 const Vector &d,
230 Vector &y,
231 const int d1d = 0,
232 const int q1d = 0)
233{
234 const int D1D = T_D1D ? T_D1D : d1d;
235 const int Q1D = T_Q1D ? T_Q1D : q1d;
236 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
237 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
238 auto B = Reshape(b.Read(), Q1D, D1D);
239 auto G = Reshape(g.Read(), Q1D, D1D);
240 // note the different shape for D, this is a (symmetric) matrix so we only
241 // store necessary entries
242 auto D = Reshape(d.Read(), Q1D*Q1D, 3, NE);
243 auto Y = Reshape(y.ReadWrite(), D1D, D1D, 2, NE);
244 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
245 {
246 const int D1D = T_D1D ? T_D1D : d1d;
247 const int Q1D = T_Q1D ? T_Q1D : q1d;
248 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
249 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
250 // gradphi \cdot Q \gradphi has four terms
251 real_t QD0[MQ1][MD1];
252 real_t QD1[MQ1][MD1];
253 real_t QD2[MQ1][MD1];
254 for (int qx = 0; qx < Q1D; ++qx)
255 {
256 for (int dy = 0; dy < D1D; ++dy)
257 {
258 QD0[qx][dy] = 0.0;
259 QD1[qx][dy] = 0.0;
260 QD2[qx][dy] = 0.0;
261 for (int qy = 0; qy < Q1D; ++qy)
262 {
263 const int q = qx + qy * Q1D;
264 const real_t D0 = D(q,0,e);
265 const real_t D1 = D(q,1,e);
266 const real_t D2 = D(q,2,e);
267 QD0[qx][dy] += B(qy, dy) * B(qy, dy) * D0;
268 QD1[qx][dy] += B(qy, dy) * G(qy, dy) * D1;
269 QD2[qx][dy] += G(qy, dy) * G(qy, dy) * D2;
270 }
271 }
272 }
273 for (int dy = 0; dy < D1D; ++dy)
274 {
275 for (int dx = 0; dx < D1D; ++dx)
276 {
277 real_t temp = 0.0;
278 for (int qx = 0; qx < Q1D; ++qx)
279 {
280 temp += G(qx, dx) * G(qx, dx) * QD0[qx][dy];
281 temp += G(qx, dx) * B(qx, dx) * QD1[qx][dy];
282 temp += B(qx, dx) * G(qx, dx) * QD1[qx][dy];
283 temp += B(qx, dx) * B(qx, dx) * QD2[qx][dy];
284 }
285 Y(dx,dy,0,e) += temp;
286 Y(dx,dy,1,e) += temp;
287 }
288 }
289 });
290}
291
292template<int T_D1D = 0, int T_Q1D = 0>
293static void PAVectorDiffusionDiagonal3D(const int NE,
294 const Array<real_t> &b,
295 const Array<real_t> &g,
296 const Vector &d,
297 Vector &y,
298 const int d1d = 0,
299 const int q1d = 0)
300{
301 constexpr int DIM = 3;
302 const int D1D = T_D1D ? T_D1D : d1d;
303 const int Q1D = T_Q1D ? T_Q1D : q1d;
304 const int max_q1d = T_Q1D ? T_Q1D : DeviceDofQuadLimits::Get().MAX_Q1D;
305 const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
306 MFEM_VERIFY(D1D <= max_d1d, "");
307 MFEM_VERIFY(Q1D <= max_q1d, "");
308 auto B = Reshape(b.Read(), Q1D, D1D);
309 auto G = Reshape(g.Read(), Q1D, D1D);
310 auto Q = Reshape(d.Read(), Q1D*Q1D*Q1D, 6, NE);
311 auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, 3, NE);
312 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
313 {
314 const int D1D = T_D1D ? T_D1D : d1d;
315 const int Q1D = T_Q1D ? T_Q1D : q1d;
316 constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
317 constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
318 real_t QQD[MQ1][MQ1][MD1];
319 real_t QDD[MQ1][MD1][MD1];
320 for (int i = 0; i < DIM; ++i)
321 {
322 for (int j = 0; j < DIM; ++j)
323 {
324 // first tensor contraction, along z direction
325 for (int qx = 0; qx < Q1D; ++qx)
326 {
327 for (int qy = 0; qy < Q1D; ++qy)
328 {
329 for (int dz = 0; dz < D1D; ++dz)
330 {
331 QQD[qx][qy][dz] = 0.0;
332 for (int qz = 0; qz < Q1D; ++qz)
333 {
334 const int q = qx + (qy + qz * Q1D) * Q1D;
335 const int k = j >= i ?
336 3 - (3-i)*(2-i)/2 + j:
337 3 - (3-j)*(2-j)/2 + i;
338 const real_t O = Q(q,k,e);
339 const real_t Bz = B(qz,dz);
340 const real_t Gz = G(qz,dz);
341 const real_t L = i==2 ? Gz : Bz;
342 const real_t R = j==2 ? Gz : Bz;
343 QQD[qx][qy][dz] += L * O * R;
344 }
345 }
346 }
347 }
348 // second tensor contraction, along y direction
349 for (int qx = 0; qx < Q1D; ++qx)
350 {
351 for (int dz = 0; dz < D1D; ++dz)
352 {
353 for (int dy = 0; dy < D1D; ++dy)
354 {
355 QDD[qx][dy][dz] = 0.0;
356 for (int qy = 0; qy < Q1D; ++qy)
357 {
358 const real_t By = B(qy,dy);
359 const real_t Gy = G(qy,dy);
360 const real_t L = i==1 ? Gy : By;
361 const real_t R = j==1 ? Gy : By;
362 QDD[qx][dy][dz] += L * QQD[qx][qy][dz] * R;
363 }
364 }
365 }
366 }
367 // third tensor contraction, along x direction
368 for (int dz = 0; dz < D1D; ++dz)
369 {
370 for (int dy = 0; dy < D1D; ++dy)
371 {
372 for (int dx = 0; dx < D1D; ++dx)
373 {
374 real_t temp = 0.0;
375 for (int qx = 0; qx < Q1D; ++qx)
376 {
377 const real_t Bx = B(qx,dx);
378 const real_t Gx = G(qx,dx);
379 const real_t L = i==0 ? Gx : Bx;
380 const real_t R = j==0 ? Gx : Bx;
381 temp += L * QDD[qx][dy][dz] * R;
382 }
383 Y(dx, dy, dz, 0, e) += temp;
384 Y(dx, dy, dz, 1, e) += temp;
385 Y(dx, dy, dz, 2, e) += temp;
386 }
387 }
388 }
389 }
390 }
391 });
392}
393
394static void PAVectorDiffusionAssembleDiagonal(const int dim,
395 const int D1D,
396 const int Q1D,
397 const int NE,
398 const Array<real_t> &B,
399 const Array<real_t> &G,
400 const Vector &op,
401 Vector &y)
402{
403 if (dim == 2)
404 {
405 return PAVectorDiffusionDiagonal2D(NE, B, G, op, y, D1D, Q1D);
406 }
407 else if (dim == 3)
408 {
409 return PAVectorDiffusionDiagonal3D(NE, B, G, op, y, D1D, Q1D);
410 }
411 MFEM_ABORT("Dimension not implemented.");
412}
413
414void VectorDiffusionIntegrator::AssembleDiagonalPA(Vector &diag)
415{
416 if (DeviceCanUseCeed())
417 {
418 ceedOp->GetDiagonal(diag);
419 }
420 else
421 {
422 PAVectorDiffusionAssembleDiagonal(dim, dofs1D, quad1D, ne,
423 maps->B, maps->G,
424 pa_data, diag);
425 }
426}
427
428// PA Diffusion Apply 2D kernel
429template<int T_D1D = 0, int T_Q1D = 0, int T_VDIM = 0> static
430void PAVectorDiffusionApply2D(const int NE,
431 const Array<real_t> &b,
432 const Array<real_t> &g,
433 const Array<real_t> &bt,
434 const Array<real_t> &gt,
435 const Vector &d_,
436 const Vector &x_,
437 Vector &y_,
438 const int d1d = 0,
439 const int q1d = 0,
440 const int vdim = 0)
441{
442 const int D1D = T_D1D ? T_D1D : d1d;
443 const int Q1D = T_Q1D ? T_Q1D : q1d;
444 const int VDIM = T_VDIM ? T_VDIM : vdim;
445 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
446 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
447 auto B = Reshape(b.Read(), Q1D, D1D);
448 auto G = Reshape(g.Read(), Q1D, D1D);
449 auto Bt = Reshape(bt.Read(), D1D, Q1D);
450 auto Gt = Reshape(gt.Read(), D1D, Q1D);
451 auto D = Reshape(d_.Read(), Q1D*Q1D, 3, NE);
452 auto x = Reshape(x_.Read(), D1D, D1D, VDIM, NE);
453 auto y = Reshape(y_.ReadWrite(), D1D, D1D, VDIM, NE);
454 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
455 {
456 const int D1D = T_D1D ? T_D1D : d1d;
457 const int Q1D = T_Q1D ? T_Q1D : q1d;
458 const int VDIM = T_VDIM ? T_VDIM : vdim;
459 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
460 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
461
462 real_t grad[max_Q1D][max_Q1D][2];
463 for (int c = 0; c < VDIM; c++)
464 {
465 for (int qy = 0; qy < Q1D; ++qy)
466 {
467 for (int qx = 0; qx < Q1D; ++qx)
468 {
469 grad[qy][qx][0] = 0.0;
470 grad[qy][qx][1] = 0.0;
471 }
472 }
473 for (int dy = 0; dy < D1D; ++dy)
474 {
475 real_t gradX[max_Q1D][2];
476 for (int qx = 0; qx < Q1D; ++qx)
477 {
478 gradX[qx][0] = 0.0;
479 gradX[qx][1] = 0.0;
480 }
481 for (int dx = 0; dx < D1D; ++dx)
482 {
483 const real_t s = x(dx,dy,c,e);
484 for (int qx = 0; qx < Q1D; ++qx)
485 {
486 gradX[qx][0] += s * B(qx,dx);
487 gradX[qx][1] += s * G(qx,dx);
488 }
489 }
490 for (int qy = 0; qy < Q1D; ++qy)
491 {
492 const real_t wy = B(qy,dy);
493 const real_t wDy = G(qy,dy);
494 for (int qx = 0; qx < Q1D; ++qx)
495 {
496 grad[qy][qx][0] += gradX[qx][1] * wy;
497 grad[qy][qx][1] += gradX[qx][0] * wDy;
498 }
499 }
500 }
501 // Calculate Dxy, xDy in plane
502 for (int qy = 0; qy < Q1D; ++qy)
503 {
504 for (int qx = 0; qx < Q1D; ++qx)
505 {
506 const int q = qx + qy * Q1D;
507 const real_t O11 = D(q,0,e);
508 const real_t O12 = D(q,1,e);
509 const real_t O22 = D(q,2,e);
510 const real_t gradX = grad[qy][qx][0];
511 const real_t gradY = grad[qy][qx][1];
512 grad[qy][qx][0] = (O11 * gradX) + (O12 * gradY);
513 grad[qy][qx][1] = (O12 * gradX) + (O22 * gradY);
514 }
515 }
516 for (int qy = 0; qy < Q1D; ++qy)
517 {
518 real_t gradX[max_D1D][2];
519 for (int dx = 0; dx < D1D; ++dx)
520 {
521 gradX[dx][0] = 0.0;
522 gradX[dx][1] = 0.0;
523 }
524 for (int qx = 0; qx < Q1D; ++qx)
525 {
526 const real_t gX = grad[qy][qx][0];
527 const real_t gY = grad[qy][qx][1];
528 for (int dx = 0; dx < D1D; ++dx)
529 {
530 const real_t wx = Bt(dx,qx);
531 const real_t wDx = Gt(dx,qx);
532 gradX[dx][0] += gX * wDx;
533 gradX[dx][1] += gY * wx;
534 }
535 }
536 for (int dy = 0; dy < D1D; ++dy)
537 {
538 const real_t wy = Bt(dy,qy);
539 const real_t wDy = Gt(dy,qy);
540 for (int dx = 0; dx < D1D; ++dx)
541 {
542 y(dx,dy,c,e) += ((gradX[dx][0] * wy) + (gradX[dx][1] * wDy));
543 }
544 }
545 }
546 }
547 });
548}
549
550// PA Diffusion Apply 3D kernel
551template<const int T_D1D = 0,
552 const int T_Q1D = 0> static
553void PAVectorDiffusionApply3D(const int NE,
554 const Array<real_t> &b,
555 const Array<real_t> &g,
556 const Array<real_t> &bt,
557 const Array<real_t> &gt,
558 const Vector &op_,
559 const Vector &x_,
560 Vector &y_,
561 int d1d = 0, int q1d = 0)
562{
563 const int D1D = T_D1D ? T_D1D : d1d;
564 const int Q1D = T_Q1D ? T_Q1D : q1d;
565 constexpr int VDIM = 3;
566 MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
567 MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
568 auto B = Reshape(b.Read(), Q1D, D1D);
569 auto G = Reshape(g.Read(), Q1D, D1D);
570 auto Bt = Reshape(bt.Read(), D1D, Q1D);
571 auto Gt = Reshape(gt.Read(), D1D, Q1D);
572 auto op = Reshape(op_.Read(), Q1D*Q1D*Q1D, 6, NE);
573 auto x = Reshape(x_.Read(), D1D, D1D, D1D, VDIM, NE);
574 auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
575 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
576 {
577 const int D1D = T_D1D ? T_D1D : d1d;
578 const int Q1D = T_Q1D ? T_Q1D : q1d;
579 constexpr int max_D1D = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
580 constexpr int max_Q1D = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
581 for (int c = 0; c < VDIM; ++ c)
582 {
583 real_t grad[max_Q1D][max_Q1D][max_Q1D][3];
584 for (int qz = 0; qz < Q1D; ++qz)
585 {
586 for (int qy = 0; qy < Q1D; ++qy)
587 {
588 for (int qx = 0; qx < Q1D; ++qx)
589 {
590 grad[qz][qy][qx][0] = 0.0;
591 grad[qz][qy][qx][1] = 0.0;
592 grad[qz][qy][qx][2] = 0.0;
593 }
594 }
595 }
596 for (int dz = 0; dz < D1D; ++dz)
597 {
598 real_t gradXY[max_Q1D][max_Q1D][3];
599 for (int qy = 0; qy < Q1D; ++qy)
600 {
601 for (int qx = 0; qx < Q1D; ++qx)
602 {
603 gradXY[qy][qx][0] = 0.0;
604 gradXY[qy][qx][1] = 0.0;
605 gradXY[qy][qx][2] = 0.0;
606 }
607 }
608 for (int dy = 0; dy < D1D; ++dy)
609 {
610 real_t gradX[max_Q1D][2];
611 for (int qx = 0; qx < Q1D; ++qx)
612 {
613 gradX[qx][0] = 0.0;
614 gradX[qx][1] = 0.0;
615 }
616 for (int dx = 0; dx < D1D; ++dx)
617 {
618 const real_t s = x(dx,dy,dz,c,e);
619 for (int qx = 0; qx < Q1D; ++qx)
620 {
621 gradX[qx][0] += s * B(qx,dx);
622 gradX[qx][1] += s * G(qx,dx);
623 }
624 }
625 for (int qy = 0; qy < Q1D; ++qy)
626 {
627 const real_t wy = B(qy,dy);
628 const real_t wDy = G(qy,dy);
629 for (int qx = 0; qx < Q1D; ++qx)
630 {
631 const real_t wx = gradX[qx][0];
632 const real_t wDx = gradX[qx][1];
633 gradXY[qy][qx][0] += wDx * wy;
634 gradXY[qy][qx][1] += wx * wDy;
635 gradXY[qy][qx][2] += wx * wy;
636 }
637 }
638 }
639 for (int qz = 0; qz < Q1D; ++qz)
640 {
641 const real_t wz = B(qz,dz);
642 const real_t wDz = G(qz,dz);
643 for (int qy = 0; qy < Q1D; ++qy)
644 {
645 for (int qx = 0; qx < Q1D; ++qx)
646 {
647 grad[qz][qy][qx][0] += gradXY[qy][qx][0] * wz;
648 grad[qz][qy][qx][1] += gradXY[qy][qx][1] * wz;
649 grad[qz][qy][qx][2] += gradXY[qy][qx][2] * wDz;
650 }
651 }
652 }
653 }
654 // Calculate Dxyz, xDyz, xyDz in plane
655 for (int qz = 0; qz < Q1D; ++qz)
656 {
657 for (int qy = 0; qy < Q1D; ++qy)
658 {
659 for (int qx = 0; qx < Q1D; ++qx)
660 {
661 const int q = qx + (qy + qz * Q1D) * Q1D;
662 const real_t O11 = op(q,0,e);
663 const real_t O12 = op(q,1,e);
664 const real_t O13 = op(q,2,e);
665 const real_t O22 = op(q,3,e);
666 const real_t O23 = op(q,4,e);
667 const real_t O33 = op(q,5,e);
668 const real_t gradX = grad[qz][qy][qx][0];
669 const real_t gradY = grad[qz][qy][qx][1];
670 const real_t gradZ = grad[qz][qy][qx][2];
671 grad[qz][qy][qx][0] = (O11*gradX)+(O12*gradY)+(O13*gradZ);
672 grad[qz][qy][qx][1] = (O12*gradX)+(O22*gradY)+(O23*gradZ);
673 grad[qz][qy][qx][2] = (O13*gradX)+(O23*gradY)+(O33*gradZ);
674 }
675 }
676 }
677 for (int qz = 0; qz < Q1D; ++qz)
678 {
679 real_t gradXY[max_D1D][max_D1D][3];
680 for (int dy = 0; dy < D1D; ++dy)
681 {
682 for (int dx = 0; dx < D1D; ++dx)
683 {
684 gradXY[dy][dx][0] = 0;
685 gradXY[dy][dx][1] = 0;
686 gradXY[dy][dx][2] = 0;
687 }
688 }
689 for (int qy = 0; qy < Q1D; ++qy)
690 {
691 real_t gradX[max_D1D][3];
692 for (int dx = 0; dx < D1D; ++dx)
693 {
694 gradX[dx][0] = 0;
695 gradX[dx][1] = 0;
696 gradX[dx][2] = 0;
697 }
698 for (int qx = 0; qx < Q1D; ++qx)
699 {
700 const real_t gX = grad[qz][qy][qx][0];
701 const real_t gY = grad[qz][qy][qx][1];
702 const real_t gZ = grad[qz][qy][qx][2];
703 for (int dx = 0; dx < D1D; ++dx)
704 {
705 const real_t wx = Bt(dx,qx);
706 const real_t wDx = Gt(dx,qx);
707 gradX[dx][0] += gX * wDx;
708 gradX[dx][1] += gY * wx;
709 gradX[dx][2] += gZ * wx;
710 }
711 }
712 for (int dy = 0; dy < D1D; ++dy)
713 {
714 const real_t wy = Bt(dy,qy);
715 const real_t wDy = Gt(dy,qy);
716 for (int dx = 0; dx < D1D; ++dx)
717 {
718 gradXY[dy][dx][0] += gradX[dx][0] * wy;
719 gradXY[dy][dx][1] += gradX[dx][1] * wDy;
720 gradXY[dy][dx][2] += gradX[dx][2] * wy;
721 }
722 }
723 }
724 for (int dz = 0; dz < D1D; ++dz)
725 {
726 const real_t wz = Bt(dz,qz);
727 const real_t wDz = Gt(dz,qz);
728 for (int dy = 0; dy < D1D; ++dy)
729 {
730 for (int dx = 0; dx < D1D; ++dx)
731 {
732 y(dx,dy,dz,c,e) +=
733 ((gradXY[dy][dx][0] * wz) +
734 (gradXY[dy][dx][1] * wz) +
735 (gradXY[dy][dx][2] * wDz));
736 }
737 }
738 }
739 }
740 }
741 });
742}
743
744// PA Diffusion Apply kernel
745void VectorDiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
746{
747 if (DeviceCanUseCeed())
748 {
749 ceedOp->AddMult(x, y);
750 }
751 else
752 {
753 const int D1D = dofs1D;
754 const int Q1D = quad1D;
755 const Array<real_t> &B = maps->B;
756 const Array<real_t> &G = maps->G;
757 const Array<real_t> &Bt = maps->Bt;
758 const Array<real_t> &Gt = maps->Gt;
759 const Vector &D = pa_data;
760
761 if (dim == 2 && sdim == 3)
762 {
763 switch ((dofs1D << 4 ) | quad1D)
764 {
765 case 0x22: return PAVectorDiffusionApply2D<2,2,3>(ne,B,G,Bt,Gt,D,x,y);
766 case 0x33: return PAVectorDiffusionApply2D<3,3,3>(ne,B,G,Bt,Gt,D,x,y);
767 case 0x44: return PAVectorDiffusionApply2D<4,4,3>(ne,B,G,Bt,Gt,D,x,y);
768 case 0x55: return PAVectorDiffusionApply2D<5,5,3>(ne,B,G,Bt,Gt,D,x,y);
769 default:
770 return PAVectorDiffusionApply2D(ne,B,G,Bt,Gt,D,x,y,D1D,Q1D,sdim);
771 }
772 }
773 if (dim == 2 && sdim == 2)
774 { return PAVectorDiffusionApply2D(ne,B,G,Bt,Gt,D,x,y,D1D,Q1D,sdim); }
775
776 if (dim == 3 && sdim == 3)
777 { return PAVectorDiffusionApply3D(ne,B,G,Bt,Gt,D,x,y,D1D,Q1D); }
778
779 MFEM_ABORT("Unknown kernel.");
780 }
781}
782
783} // namespace mfem
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:337
mfem::CoefficientVector
Class to represent a coefficient evaluated at quadrature points.
Definition coefficient.hpp:2429
mfem::Device::GetDeviceMemoryType
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition device.hpp:274
mfem::DiffusionIntegrator::GetRule
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe)
Definition bilininteg.cpp:1318
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::IsVariableOrder
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:709
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::GetDofToQuad
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:365
mfem::FiniteElement::GetDim
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:321
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::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::SpaceDimension
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
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::Mesh::GetNumGeometries
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:7243
mfem::QuadratureSpace
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:120
mfem::VectorDiffusionIntegrator::AssemblePA
void AssemblePA(const FiniteElementSpace &fes) override
Method defining partial assembly.
Definition bilininteg_vecdiffusion_pa.cpp:139
mfem::VectorDiffusionIntegrator::maps
const DofToQuad * maps
Not owned.
Definition bilininteg.hpp:3062
mfem::VectorDiffusionIntegrator::AddMultPA
void AddMultPA(const Vector &x, Vector &y) const override
Method for partially assembled action.
Definition bilininteg_vecdiffusion_pa.cpp:745
mfem::VectorDiffusionIntegrator::geom
const GeometricFactors * geom
Not owned.
Definition bilininteg.hpp:3063
mfem::VectorDiffusionIntegrator::AssembleDiagonalPA
void AssembleDiagonalPA(Vector &diag) override
Assemble diagonal and add it to Vector diag.
Definition bilininteg_vecdiffusion_pa.cpp:414
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::Read
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:494
mfem::Vector::ReadWrite
virtual real_t * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:510
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
mfem::Vector::Write
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:502
mfem::ceed::Operator::GetDiagonal
void GetDiagonal(mfem::Vector &diag) const
Definition operator.cpp:104
mfem::ceed::Operator::AddMult
void AddMult(const mfem::Vector &x, mfem::Vector &y, const real_t a=1.0) const override
Operator application: y+=A(x) (default) or y+=a*A(x).
Definition operator.cpp:72
mfem::ceed::PADiffusionIntegrator
Represent a DiffusionIntegrator with AssemblyLevel::Partial using libCEED.
Definition diffusion.hpp:27
alpha
const real_t alpha
Definition ex15.cpp:369
dim
int dim
Definition ex24.cpp:53
b
real_t b
Definition lissajous.cpp:42
SDIM
constexpr int SDIM
Definition minimal-surface.cpp:73
DIM
constexpr int DIM
Definition minimal-surface.cpp:72
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:341
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::DeviceCanUseCeed
bool DeviceCanUseCeed()
Function that determines if a CEED kernel should be used, based on the current mfem::Device configura...
Definition util.cpp:33
mfem::CoefficientStorage::COMPRESSED
@ COMPRESSED
Enable all above compressions.
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
mfem::DeviceDofQuadLimits::MAX_D1D
int MAX_D1D
Maximum number of 1D nodal points.
Definition forall.hpp:118
mfem::DeviceDofQuadLimits::MAX_Q1D
int MAX_Q1D
Maximum number of 1D quadrature points.
Definition forall.hpp:119
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0