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