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