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