MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilininteg_dgtrace.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 "restriction.hpp"
16 
17 using namespace std;
18 
19 namespace mfem
20 {
21 // PA DG Trace Integrator
22 static void PADGTraceSetup2D(const int Q1D,
23  const int NF,
24  const Array<double> &w,
25  const Vector &det,
26  const Vector &nor,
27  const Vector &rho,
28  const Vector &vel,
29  const double alpha,
30  const double beta,
31  Vector &op)
32 {
33  const int VDIM = 2;
34 
35  auto d = Reshape(det.Read(), Q1D, NF);
36  auto n = Reshape(nor.Read(), Q1D, VDIM, NF);
37  const bool const_r = rho.Size() == 1;
38  auto R =
39  const_r ? Reshape(rho.Read(), 1,1) : Reshape(rho.Read(), Q1D,NF);
40  const bool const_v = vel.Size() == 2;
41  auto V =
42  const_v ? Reshape(vel.Read(), 2,1,1) : Reshape(vel.Read(), 2,Q1D,NF);
43  auto W = w.Read();
44  auto qd = Reshape(op.Write(), Q1D, 2, 2, NF);
45 
46  MFEM_FORALL(f, NF,//can be optimized with Q1D thread for NF blocks
47  {
48  for (int q = 0; q < Q1D; ++q)
49  {
50  const double r = const_r ? R(0,0) : R(q,f);
51  const double v0 = const_v ? V(0,0,0) : V(0,q,f);
52  const double v1 = const_v ? V(1,0,0) : V(1,q,f);
53  const double dot = n(q,0,f) * v0 + n(q,1,f) * v1;
54  const double abs = dot > 0.0 ? dot : -dot;
55  const double w = W[q]*r*d(q,f);
56  qd(q,0,0,f) = w*( alpha/2 * dot + beta * abs );
57  qd(q,1,0,f) = w*( alpha/2 * dot - beta * abs );
58  qd(q,0,1,f) = w*(-alpha/2 * dot - beta * abs );
59  qd(q,1,1,f) = w*(-alpha/2 * dot + beta * abs );
60  }
61  });
62 }
63 
64 static void PADGTraceSetup3D(const int Q1D,
65  const int NF,
66  const Array<double> &w,
67  const Vector &det,
68  const Vector &nor,
69  const Vector &rho,
70  const Vector &vel,
71  const double alpha,
72  const double beta,
73  Vector &op)
74 {
75  const int VDIM = 3;
76 
77  auto d = Reshape(det.Read(), Q1D, Q1D, NF);
78  auto n = Reshape(nor.Read(), Q1D, Q1D, VDIM, NF);
79  const bool const_r = rho.Size() == 1;
80  auto R =
81  const_r ? Reshape(rho.Read(), 1,1,1) : Reshape(rho.Read(), Q1D,Q1D,NF);
82  const bool const_v = vel.Size() == 3;
83  auto V =
84  const_v ? Reshape(vel.Read(), 3,1,1,1) : Reshape(vel.Read(), 3,Q1D,Q1D,NF);
85  auto W = w.Read();
86  auto qd = Reshape(op.Write(), Q1D, Q1D, 2, 2, NF);
87 
88  MFEM_FORALL(f, NF,//can be optimized with Q1D*Q1D threads for NF blocks
89  {
90  for (int q1 = 0; q1 < Q1D; ++q1)
91  {
92  for (int q2 = 0; q2 < Q1D; ++q2)
93  {
94  const double r = const_r ? R(0,0,0) : R(q1,q2,f);
95  const double v0 = const_v ? V(0,0,0,0) : V(0,q1,q2,f);
96  const double v1 = const_v ? V(1,0,0,0) : V(1,q1,q2,f);
97  const double v2 = const_v ? V(2,0,0,0) : V(2,q1,q2,f);
98  const double dot = n(q1,q2,0,f) * v0 + n(q1,q2,1,f) * v1 +
99  /* */ n(q1,q2,2,f) * v2;
100  const double abs = dot > 0.0 ? dot : -dot;
101  const double w = W[q1+q2*Q1D]*r*d(q1,q2,f);
102  qd(q1,q2,0,0,f) = w*( alpha/2 * dot + beta * abs );
103  qd(q1,q2,1,0,f) = w*( alpha/2 * dot - beta * abs );
104  qd(q1,q2,0,1,f) = w*(-alpha/2 * dot - beta * abs );
105  qd(q1,q2,1,1,f) = w*(-alpha/2 * dot + beta * abs );
106  }
107  }
108  });
109 }
110 
111 static void PADGTraceSetup(const int dim,
112  const int D1D,
113  const int Q1D,
114  const int NF,
115  const Array<double> &W,
116  const Vector &det,
117  const Vector &nor,
118  const Vector &rho,
119  const Vector &u,
120  const double alpha,
121  const double beta,
122  Vector &op)
123 {
124  if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADGTraceSetup"); }
125  if (dim == 2)
126  {
127  PADGTraceSetup2D(Q1D, NF, W, det, nor, rho, u, alpha, beta, op);
128  }
129  if (dim == 3)
130  {
131  PADGTraceSetup3D(Q1D, NF, W, det, nor, rho, u, alpha, beta, op);
132  }
133 }
134 
135 void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
136 {
137  nf = fes.GetNFbyType(type);
138  if (nf==0) { return; }
139  // Assumes tensor-product elements
140  Mesh *mesh = fes.GetMesh();
141  const FiniteElement &el =
142  *fes.GetTraceElement(0, fes.GetMesh()->GetFaceBaseGeometry(0));
143  FaceElementTransformations &T =
144  *fes.GetMesh()->GetFaceElementTransformations(0);
145  const IntegrationRule *ir = IntRule?
146  IntRule:
147  &GetRule(el.GetGeomType(), el.GetOrder(), T);
148  const int symmDims = 4;
149  const int nq = ir->GetNPoints();
150  dim = mesh->Dimension();
151  geom = mesh->GetFaceGeometricFactors(
152  *ir,
153  FaceGeometricFactors::DETERMINANTS |
154  FaceGeometricFactors::NORMALS, type);
155  maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
156  dofs1D = maps->ndof;
157  quad1D = maps->nqpt;
158  pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
159  Vector r;
160  if (rho==nullptr)
161  {
162  r.SetSize(1);
163  r(0) = 1.0;
164  }
165  else if (ConstantCoefficient *c_rho = dynamic_cast<ConstantCoefficient*>(rho))
166  {
167  r.SetSize(1);
168  r(0) = c_rho->constant;
169  }
170  else
171  {
172  r.SetSize(nq * nf);
173  auto C = Reshape(r.HostWrite(), nq, nf);
174  int f_ind = 0;
175  for (int f = 0; f < fes.GetNF(); ++f)
176  {
177  int e1, e2;
178  int inf1, inf2;
179  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
180  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
181  int face_id = inf1 / 64;
182  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
183  (type==FaceType::Boundary && e2<0 && inf2<0) )
184  {
185  ElementTransformation& T = *fes.GetMesh()->GetFaceTransformation(f);
186  for (int q = 0; q < nq; ++q)
187  {
188  // Convert to lexicographic ordering
189  int iq = ToLexOrdering(dim, face_id, quad1D, q);
190  C(iq,f_ind) = rho->Eval(T, ir->IntPoint(q));
191  }
192  f_ind++;
193  }
194  }
195  MFEM_VERIFY(f_ind==nf, "Incorrect number of faces.");
196  }
197  Vector vel;
198  if (VectorConstantCoefficient *c_u = dynamic_cast<VectorConstantCoefficient*>
199  (u))
200  {
201  vel = c_u->GetVec();
202  }
203  else
204  {
205  vel.SetSize(dim * nq * nf);
206  auto C = Reshape(vel.HostWrite(), dim, nq, nf);
207  Vector Vq(dim);
208  int f_ind = 0;
209  for (int f = 0; f < fes.GetNF(); ++f)
210  {
211  int e1, e2;
212  int inf1, inf2;
213  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
214  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
215  int face_id = inf1 / 64;
216  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
217  (type==FaceType::Boundary && e2<0 && inf2<0) )
218  {
219  ElementTransformation& T = *fes.GetMesh()->GetFaceTransformation(f);
220  for (int q = 0; q < nq; ++q)
221  {
222  // Convert to lexicographic ordering
223  int iq = ToLexOrdering(dim, face_id, quad1D, q);
224  u->Eval(Vq, T, ir->IntPoint(q));
225  for (int i = 0; i < dim; ++i)
226  {
227  C(i,iq,f_ind) = Vq(i);
228  }
229  }
230  f_ind++;
231  }
232  }
233  MFEM_VERIFY(f_ind==nf, "Incorrect number of faces.");
234  }
235  PADGTraceSetup(dim, dofs1D, quad1D, nf, ir->GetWeights(),
236  geom->detJ, geom->normal, r, vel,
237  alpha, beta, pa_data);
238 }
239 
240 void DGTraceIntegrator::AssemblePAInteriorFaces(const FiniteElementSpace& fes)
241 {
242  SetupPA(fes, FaceType::Interior);
243 }
244 
245 void DGTraceIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace& fes)
246 {
247  SetupPA(fes, FaceType::Boundary);
248 }
249 
250 // PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
251 template<int T_D1D = 0, int T_Q1D = 0> static
252 void PADGTraceApply2D(const int NF,
253  const Array<double> &b,
254  const Array<double> &bt,
255  const Vector &_op,
256  const Vector &_x,
257  Vector &_y,
258  const int d1d = 0,
259  const int q1d = 0)
260 {
261  const int VDIM = 1;
262  const int D1D = T_D1D ? T_D1D : d1d;
263  const int Q1D = T_Q1D ? T_Q1D : q1d;
264  MFEM_VERIFY(D1D <= MAX_D1D, "");
265  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
266  auto B = Reshape(b.Read(), Q1D, D1D);
267  auto Bt = Reshape(bt.Read(), D1D, Q1D);
268  auto op = Reshape(_op.Read(), Q1D, 2, 2, NF);
269  auto x = Reshape(_x.Read(), D1D, VDIM, 2, NF);
270  auto y = Reshape(_y.ReadWrite(), D1D, VDIM, 2, NF);
271 
272  MFEM_FORALL(f, NF,
273  {
274  const int VDIM = 1;
275  const int D1D = T_D1D ? T_D1D : d1d;
276  const int Q1D = T_Q1D ? T_Q1D : q1d;
277  // the following variables are evaluated at compile time
278  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
279  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
280  double u0[max_D1D][VDIM];
281  double u1[max_D1D][VDIM];
282  for (int d = 0; d < D1D; d++)
283  {
284  for (int c = 0; c < VDIM; c++)
285  {
286  u0[d][c] = x(d,c,0,f);
287  u1[d][c] = x(d,c,1,f);
288  }
289  }
290  double Bu0[max_Q1D][VDIM];
291  double Bu1[max_Q1D][VDIM];
292  for (int q = 0; q < Q1D; ++q)
293  {
294  for (int c = 0; c < VDIM; c++)
295  {
296  Bu0[q][c] = 0.0;
297  Bu1[q][c] = 0.0;
298  }
299  for (int d = 0; d < D1D; ++d)
300  {
301  const double b = B(q,d);
302  for (int c = 0; c < VDIM; c++)
303  {
304  Bu0[q][c] += b*u0[d][c];
305  Bu1[q][c] += b*u1[d][c];
306  }
307  }
308  }
309  double DBu[max_Q1D][VDIM];
310  for (int q = 0; q < Q1D; ++q)
311  {
312  for (int c = 0; c < VDIM; c++)
313  {
314  DBu[q][c] = op(q,0,0,f)*Bu0[q][c] + op(q,1,0,f)*Bu1[q][c];
315  }
316  }
317  double BDBu[max_D1D][VDIM];
318  for (int d = 0; d < D1D; ++d)
319  {
320  for (int c = 0; c < VDIM; c++)
321  {
322  BDBu[d][c] = 0.0;
323  }
324  for (int q = 0; q < Q1D; ++q)
325  {
326  const double b = Bt(d,q);
327  for (int c = 0; c < VDIM; c++)
328  {
329  BDBu[d][c] += b*DBu[q][c];
330  }
331  }
332  for (int c = 0; c < VDIM; c++)
333  {
334  y(d,c,0,f) += BDBu[d][c];
335  y(d,c,1,f) += -BDBu[d][c];
336  }
337  }
338  });
339 }
340 
341 // PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
342 template<int T_D1D = 0, int T_Q1D = 0> static
343 void PADGTraceApply3D(const int NF,
344  const Array<double> &b,
345  const Array<double> &bt,
346  const Vector &_op,
347  const Vector &_x,
348  Vector &_y,
349  const int d1d = 0,
350  const int q1d = 0)
351 {
352  const int VDIM = 1;
353  const int D1D = T_D1D ? T_D1D : d1d;
354  const int Q1D = T_Q1D ? T_Q1D : q1d;
355  MFEM_VERIFY(D1D <= MAX_D1D, "");
356  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
357  auto B = Reshape(b.Read(), Q1D, D1D);
358  auto Bt = Reshape(bt.Read(), D1D, Q1D);
359  auto op = Reshape(_op.Read(), Q1D, Q1D, 2, 2, NF);
360  auto x = Reshape(_x.Read(), D1D, D1D, VDIM, 2, NF);
361  auto y = Reshape(_y.ReadWrite(), D1D, D1D, VDIM, 2, NF);
362 
363  MFEM_FORALL(f, NF,
364  {
365  const int VDIM = 1;
366  const int D1D = T_D1D ? T_D1D : d1d;
367  const int Q1D = T_Q1D ? T_Q1D : q1d;
368  // the following variables are evaluated at compile time
369  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
370  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
371  double u0[max_D1D][max_D1D][VDIM];
372  double u1[max_D1D][max_D1D][VDIM];
373  for (int d1 = 0; d1 < D1D; d1++)
374  {
375  for (int d2 = 0; d2 < D1D; d2++)
376  {
377  for (int c = 0; c < VDIM; c++)
378  {
379  u0[d1][d2][c] = x(d1,d2,c,0,f);
380  u1[d1][d2][c] = x(d1,d2,c,1,f);
381  }
382  }
383  }
384  double Bu0[max_Q1D][max_D1D][VDIM];
385  double Bu1[max_Q1D][max_D1D][VDIM];
386  for (int q = 0; q < Q1D; ++q)
387  {
388  for (int d2 = 0; d2 < D1D; d2++)
389  {
390  for (int c = 0; c < VDIM; c++)
391  {
392  Bu0[q][d2][c] = 0.0;
393  Bu1[q][d2][c] = 0.0;
394  }
395  for (int d1 = 0; d1 < D1D; ++d1)
396  {
397  const double b = B(q,d1);
398  for (int c = 0; c < VDIM; c++)
399  {
400  Bu0[q][d2][c] += b*u0[d1][d2][c];
401  Bu1[q][d2][c] += b*u1[d1][d2][c];
402  }
403  }
404  }
405  }
406  double BBu0[max_Q1D][max_Q1D][VDIM];
407  double BBu1[max_Q1D][max_Q1D][VDIM];
408  for (int q1 = 0; q1 < Q1D; ++q1)
409  {
410  for (int q2 = 0; q2 < Q1D; q2++)
411  {
412  for (int c = 0; c < VDIM; c++)
413  {
414  BBu0[q1][q2][c] = 0.0;
415  BBu1[q1][q2][c] = 0.0;
416  }
417  for (int d2 = 0; d2 < D1D; ++d2)
418  {
419  const double b = B(q2,d2);
420  for (int c = 0; c < VDIM; c++)
421  {
422  BBu0[q1][q2][c] += b*Bu0[q1][d2][c];
423  BBu1[q1][q2][c] += b*Bu1[q1][d2][c];
424  }
425  }
426  }
427  }
428  double DBBu[max_Q1D][max_Q1D][VDIM];
429  for (int q1 = 0; q1 < Q1D; ++q1)
430  {
431  for (int q2 = 0; q2 < Q1D; q2++)
432  {
433  for (int c = 0; c < VDIM; c++)
434  {
435  DBBu[q1][q2][c] = op(q1,q2,0,0,f)*BBu0[q1][q2][c] +
436  op(q1,q2,1,0,f)*BBu1[q1][q2][c];
437  }
438  }
439  }
440  double BDBBu[max_Q1D][max_D1D][VDIM];
441  for (int q1 = 0; q1 < Q1D; ++q1)
442  {
443  for (int d2 = 0; d2 < D1D; d2++)
444  {
445  for (int c = 0; c < VDIM; c++)
446  {
447  BDBBu[q1][d2][c] = 0.0;
448  }
449  for (int q2 = 0; q2 < Q1D; ++q2)
450  {
451  const double b = Bt(d2,q2);
452  for (int c = 0; c < VDIM; c++)
453  {
454  BDBBu[q1][d2][c] += b*DBBu[q1][q2][c];
455  }
456  }
457  }
458  }
459  double BBDBBu[max_D1D][max_D1D][VDIM];
460  for (int d1 = 0; d1 < D1D; ++d1)
461  {
462  for (int d2 = 0; d2 < D1D; d2++)
463  {
464  for (int c = 0; c < VDIM; c++)
465  {
466  BBDBBu[d1][d2][c] = 0.0;
467  }
468  for (int q1 = 0; q1 < Q1D; ++q1)
469  {
470  const double b = Bt(d1,q1);
471  for (int c = 0; c < VDIM; c++)
472  {
473  BBDBBu[d1][d2][c] += b*BDBBu[q1][d2][c];
474  }
475  }
476  for (int c = 0; c < VDIM; c++)
477  {
478  y(d1,d2,c,0,f) += BBDBBu[d1][d2][c];
479  y(d1,d2,c,1,f) += -BBDBBu[d1][d2][c];
480  }
481  }
482  }
483  });
484 }
485 
486 // Optimized PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
487 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
488 void SmemPADGTraceApply3D(const int NF,
489  const Array<double> &b,
490  const Array<double> &bt,
491  const Vector &_op,
492  const Vector &_x,
493  Vector &_y,
494  const int d1d = 0,
495  const int q1d = 0)
496 {
497  const int D1D = T_D1D ? T_D1D : d1d;
498  const int Q1D = T_Q1D ? T_Q1D : q1d;
499  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
500  MFEM_VERIFY(D1D <= MAX_D1D, "");
501  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
502  auto B = Reshape(b.Read(), Q1D, D1D);
503  auto Bt = Reshape(bt.Read(), D1D, Q1D);
504  auto op = Reshape(_op.Read(), Q1D, Q1D, 2, 2, NF);
505  auto x = Reshape(_x.Read(), D1D, D1D, 2, NF);
506  auto y = Reshape(_y.ReadWrite(), D1D, D1D, 2, NF);
507 
508  MFEM_FORALL_2D(f, NF, Q1D, Q1D, NBZ,
509  {
510  const int tidz = MFEM_THREAD_ID(z);
511  const int D1D = T_D1D ? T_D1D : d1d;
512  const int Q1D = T_Q1D ? T_Q1D : q1d;
513  // the following variables are evaluated at compile time
514  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
515  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
516  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
517  MFEM_SHARED double u0[NBZ][max_D1D][max_D1D];
518  MFEM_SHARED double u1[NBZ][max_D1D][max_D1D];
519  MFEM_FOREACH_THREAD(d1,x,D1D)
520  {
521  MFEM_FOREACH_THREAD(d2,y,D1D)
522  {
523  u0[tidz][d1][d2] = x(d1,d2,0,f+tidz);
524  u1[tidz][d1][d2] = x(d1,d2,1,f+tidz);
525  }
526  }
527  MFEM_SYNC_THREAD;
528  MFEM_SHARED double Bu0[NBZ][max_Q1D][max_D1D];
529  MFEM_SHARED double Bu1[NBZ][max_Q1D][max_D1D];
530  MFEM_FOREACH_THREAD(q1,x,Q1D)
531  {
532  MFEM_FOREACH_THREAD(d2,y,D1D)
533  {
534  double Bu0_ = 0.0;
535  double Bu1_ = 0.0;
536  for (int d1 = 0; d1 < D1D; ++d1)
537  {
538  const double b = B(q1,d1);
539  Bu0_ += b*u0[tidz][d1][d2];
540  Bu1_ += b*u1[tidz][d1][d2];
541  }
542  Bu0[tidz][q1][d2] = Bu0_;
543  Bu1[tidz][q1][d2] = Bu1_;
544  }
545  }
546  MFEM_SYNC_THREAD;
547  MFEM_SHARED double BBu0[NBZ][max_Q1D][max_Q1D];
548  MFEM_SHARED double BBu1[NBZ][max_Q1D][max_Q1D];
549  MFEM_FOREACH_THREAD(q1,x,Q1D)
550  {
551  MFEM_FOREACH_THREAD(q2,y,Q1D)
552  {
553  double BBu0_ = 0.0;
554  double BBu1_ = 0.0;
555  for (int d2 = 0; d2 < D1D; ++d2)
556  {
557  const double b = B(q2,d2);
558  BBu0_ += b*Bu0[tidz][q1][d2];
559  BBu1_ += b*Bu1[tidz][q1][d2];
560  }
561  BBu0[tidz][q1][q2] = BBu0_;
562  BBu1[tidz][q1][q2] = BBu1_;
563  }
564  }
565  MFEM_SYNC_THREAD;
566  MFEM_SHARED double DBBu[NBZ][max_Q1D][max_Q1D];
567  MFEM_FOREACH_THREAD(q1,x,Q1D)
568  {
569  MFEM_FOREACH_THREAD(q2,y,Q1D)
570  {
571  DBBu[tidz][q1][q2] = op(q1,q2,0,0,f+tidz)*BBu0[tidz][q1][q2] +
572  op(q1,q2,1,0,f+tidz)*BBu1[tidz][q1][q2];
573  }
574  }
575  MFEM_SYNC_THREAD;
576  MFEM_SHARED double BDBBu[NBZ][max_Q1D][max_D1D];
577  MFEM_FOREACH_THREAD(q1,x,Q1D)
578  {
579  MFEM_FOREACH_THREAD(d2,y,D1D)
580  {
581  double BDBBu_ = 0.0;
582  for (int q2 = 0; q2 < Q1D; ++q2)
583  {
584  const double b = Bt(d2,q2);
585  BDBBu_ += b*DBBu[tidz][q1][q2];
586  }
587  BDBBu[tidz][q1][d2] = BDBBu_;
588  }
589  }
590  MFEM_SYNC_THREAD;
591  MFEM_FOREACH_THREAD(d1,x,D1D)
592  {
593  MFEM_FOREACH_THREAD(d2,y,D1D)
594  {
595  double BBDBBu_ = 0.0;
596  for (int q1 = 0; q1 < Q1D; ++q1)
597  {
598  const double b = Bt(d1,q1);
599  BBDBBu_ += b*BDBBu[tidz][q1][d2];
600  }
601  y(d1,d2,0,f+tidz) += BBDBBu_;
602  y(d1,d2,1,f+tidz) += -BBDBBu_;
603  }
604  }
605  });
606 }
607 
608 static void PADGTraceApply(const int dim,
609  const int D1D,
610  const int Q1D,
611  const int NF,
612  const Array<double> &B,
613  const Array<double> &Bt,
614  const Vector &op,
615  const Vector &x,
616  Vector &y)
617 {
618  if (dim == 2)
619  {
620  switch ((D1D << 4 ) | Q1D)
621  {
622  case 0x22: return PADGTraceApply2D<2,2>(NF,B,Bt,op,x,y);
623  case 0x33: return PADGTraceApply2D<3,3>(NF,B,Bt,op,x,y);
624  case 0x44: return PADGTraceApply2D<4,4>(NF,B,Bt,op,x,y);
625  case 0x55: return PADGTraceApply2D<5,5>(NF,B,Bt,op,x,y);
626  case 0x66: return PADGTraceApply2D<6,6>(NF,B,Bt,op,x,y);
627  case 0x77: return PADGTraceApply2D<7,7>(NF,B,Bt,op,x,y);
628  case 0x88: return PADGTraceApply2D<8,8>(NF,B,Bt,op,x,y);
629  case 0x99: return PADGTraceApply2D<9,9>(NF,B,Bt,op,x,y);
630  default: return PADGTraceApply2D(NF,B,Bt,op,x,y,D1D,Q1D);
631  }
632  }
633  else if (dim == 3)
634  {
635  switch ((D1D << 4 ) | Q1D)
636  {
637  case 0x23: return SmemPADGTraceApply3D<2,3,1>(NF,B,Bt,op,x,y);
638  case 0x34: return SmemPADGTraceApply3D<3,4,2>(NF,B,Bt,op,x,y);
639  case 0x45: return SmemPADGTraceApply3D<4,5,2>(NF,B,Bt,op,x,y);
640  case 0x56: return SmemPADGTraceApply3D<5,6,1>(NF,B,Bt,op,x,y);
641  case 0x67: return SmemPADGTraceApply3D<6,7,1>(NF,B,Bt,op,x,y);
642  case 0x78: return SmemPADGTraceApply3D<7,8,1>(NF,B,Bt,op,x,y);
643  case 0x89: return SmemPADGTraceApply3D<8,9,1>(NF,B,Bt,op,x,y);
644  default: return PADGTraceApply3D(NF,B,Bt,op,x,y,D1D,Q1D);
645  }
646  }
647  MFEM_ABORT("Unknown kernel.");
648 }
649 
650 // PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
651 template<int T_D1D = 0, int T_Q1D = 0> static
652 void PADGTraceApplyTranspose2D(const int NF,
653  const Array<double> &b,
654  const Array<double> &bt,
655  const Vector &_op,
656  const Vector &_x,
657  Vector &_y,
658  const int d1d = 0,
659  const int q1d = 0)
660 {
661  const int VDIM = 1;
662  const int D1D = T_D1D ? T_D1D : d1d;
663  const int Q1D = T_Q1D ? T_Q1D : q1d;
664  MFEM_VERIFY(D1D <= MAX_D1D, "");
665  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
666  auto B = Reshape(b.Read(), Q1D, D1D);
667  auto Bt = Reshape(bt.Read(), D1D, Q1D);
668  auto op = Reshape(_op.Read(), Q1D, 2, 2, NF);
669  auto x = Reshape(_x.Read(), D1D, VDIM, 2, NF);
670  auto y = Reshape(_y.ReadWrite(), D1D, VDIM, 2, NF);
671 
672  MFEM_FORALL(f, NF,
673  {
674  const int VDIM = 1;
675  const int D1D = T_D1D ? T_D1D : d1d;
676  const int Q1D = T_Q1D ? T_Q1D : q1d;
677  // the following variables are evaluated at compile time
678  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
679  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
680  double u0[max_D1D][VDIM];
681  double u1[max_D1D][VDIM];
682  for (int d = 0; d < D1D; d++)
683  {
684  for (int c = 0; c < VDIM; c++)
685  {
686  u0[d][c] = x(d,c,0,f);
687  u1[d][c] = x(d,c,1,f);
688  }
689  }
690  double Bu0[max_Q1D][VDIM];
691  double Bu1[max_Q1D][VDIM];
692  for (int q = 0; q < Q1D; ++q)
693  {
694  for (int c = 0; c < VDIM; c++)
695  {
696  Bu0[q][c] = 0.0;
697  Bu1[q][c] = 0.0;
698  }
699  for (int d = 0; d < D1D; ++d)
700  {
701  const double b = B(q,d);
702  for (int c = 0; c < VDIM; c++)
703  {
704  Bu0[q][c] += b*u0[d][c];
705  Bu1[q][c] += b*u1[d][c];
706  }
707  }
708  }
709  double DBu0[max_Q1D][VDIM];
710  double DBu1[max_Q1D][VDIM];
711  for (int q = 0; q < Q1D; ++q)
712  {
713  for (int c = 0; c < VDIM; c++)
714  {
715  DBu0[q][c] = op(q,0,0,f)*Bu0[q][c] + op(q,0,1,f)*Bu1[q][c];
716  DBu1[q][c] = op(q,1,0,f)*Bu0[q][c] + op(q,1,1,f)*Bu1[q][c];
717  }
718  }
719  double BDBu0[max_D1D][VDIM];
720  double BDBu1[max_D1D][VDIM];
721  for (int d = 0; d < D1D; ++d)
722  {
723  for (int c = 0; c < VDIM; c++)
724  {
725  BDBu0[d][c] = 0.0;
726  BDBu1[d][c] = 0.0;
727  }
728  for (int q = 0; q < Q1D; ++q)
729  {
730  const double b = Bt(d,q);
731  for (int c = 0; c < VDIM; c++)
732  {
733  BDBu0[d][c] += b*DBu0[q][c];
734  BDBu1[d][c] += b*DBu1[q][c];
735  }
736  }
737  for (int c = 0; c < VDIM; c++)
738  {
739  y(d,c,0,f) += BDBu0[d][c];
740  y(d,c,1,f) += BDBu1[d][c];
741  }
742  }
743  });
744 }
745 
746 // PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
747 template<int T_D1D = 0, int T_Q1D = 0> static
748 void PADGTraceApplyTranspose3D(const int NF,
749  const Array<double> &b,
750  const Array<double> &bt,
751  const Vector &_op,
752  const Vector &_x,
753  Vector &_y,
754  const int d1d = 0,
755  const int q1d = 0)
756 {
757  const int VDIM = 1;
758  const int D1D = T_D1D ? T_D1D : d1d;
759  const int Q1D = T_Q1D ? T_Q1D : q1d;
760  MFEM_VERIFY(D1D <= MAX_D1D, "");
761  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
762  auto B = Reshape(b.Read(), Q1D, D1D);
763  auto Bt = Reshape(bt.Read(), D1D, Q1D);
764  auto op = Reshape(_op.Read(), Q1D, Q1D, 2, 2, NF);
765  auto x = Reshape(_x.Read(), D1D, D1D, VDIM, 2, NF);
766  auto y = Reshape(_y.ReadWrite(), D1D, D1D, VDIM, 2, NF);
767 
768  MFEM_FORALL(f, NF,
769  {
770  const int VDIM = 1;
771  const int D1D = T_D1D ? T_D1D : d1d;
772  const int Q1D = T_Q1D ? T_Q1D : q1d;
773  // the following variables are evaluated at compile time
774  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
775  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
776  double u0[max_D1D][max_D1D][VDIM];
777  double u1[max_D1D][max_D1D][VDIM];
778  for (int d1 = 0; d1 < D1D; d1++)
779  {
780  for (int d2 = 0; d2 < D1D; d2++)
781  {
782  for (int c = 0; c < VDIM; c++)
783  {
784  u0[d1][d2][c] = x(d1,d2,c,0,f);
785  u1[d1][d2][c] = x(d1,d2,c,1,f);
786  }
787  }
788  }
789  double Bu0[max_Q1D][max_D1D][VDIM];
790  double Bu1[max_Q1D][max_D1D][VDIM];
791  for (int q1 = 0; q1 < Q1D; ++q1)
792  {
793  for (int d2 = 0; d2 < D1D; ++d2)
794  {
795  for (int c = 0; c < VDIM; c++)
796  {
797  Bu0[q1][d2][c] = 0.0;
798  Bu1[q1][d2][c] = 0.0;
799  }
800  for (int d1 = 0; d1 < D1D; ++d1)
801  {
802  const double b = B(q1,d1);
803  for (int c = 0; c < VDIM; c++)
804  {
805  Bu0[q1][d2][c] += b*u0[d1][d2][c];
806  Bu1[q1][d2][c] += b*u1[d1][d2][c];
807  }
808  }
809  }
810  }
811  double BBu0[max_Q1D][max_Q1D][VDIM];
812  double BBu1[max_Q1D][max_Q1D][VDIM];
813  for (int q1 = 0; q1 < Q1D; ++q1)
814  {
815  for (int q2 = 0; q2 < Q1D; ++q2)
816  {
817  for (int c = 0; c < VDIM; c++)
818  {
819  BBu0[q1][q2][c] = 0.0;
820  BBu1[q1][q2][c] = 0.0;
821  }
822  for (int d2 = 0; d2 < D1D; ++d2)
823  {
824  const double b = B(q2,d2);
825  for (int c = 0; c < VDIM; c++)
826  {
827  BBu0[q1][q2][c] += b*Bu0[q1][d2][c];
828  BBu1[q1][q2][c] += b*Bu1[q1][d2][c];
829  }
830  }
831  }
832  }
833  double DBu0[max_Q1D][max_Q1D][VDIM];
834  double DBu1[max_Q1D][max_Q1D][VDIM];
835  for (int q1 = 0; q1 < Q1D; ++q1)
836  {
837  for (int q2 = 0; q2 < Q1D; ++q2)
838  {
839  const double D00 = op(q1,q2,0,0,f);
840  const double D01 = op(q1,q2,0,1,f);
841  const double D10 = op(q1,q2,1,0,f);
842  const double D11 = op(q1,q2,1,1,f);
843  for (int c = 0; c < VDIM; c++)
844  {
845  DBu0[q1][q2][c] = D00*BBu0[q1][q2][c] + D01*BBu1[q1][q2][c];
846  DBu1[q1][q2][c] = D10*BBu0[q1][q2][c] + D11*BBu1[q1][q2][c];
847  }
848  }
849  }
850  double BDBu0[max_D1D][max_Q1D][VDIM];
851  double BDBu1[max_D1D][max_Q1D][VDIM];
852  for (int d1 = 0; d1 < D1D; ++d1)
853  {
854  for (int q2 = 0; q2 < Q1D; ++q2)
855  {
856  for (int c = 0; c < VDIM; c++)
857  {
858  BDBu0[d1][q2][c] = 0.0;
859  BDBu1[d1][q2][c] = 0.0;
860  }
861  for (int q1 = 0; q1 < Q1D; ++q1)
862  {
863  const double b = Bt(d1,q1);
864  for (int c = 0; c < VDIM; c++)
865  {
866  BDBu0[d1][q2][c] += b*DBu0[q1][q2][c];
867  BDBu1[d1][q2][c] += b*DBu1[q1][q2][c];
868  }
869  }
870  }
871  }
872  double BBDBu0[max_D1D][max_D1D][VDIM];
873  double BBDBu1[max_D1D][max_D1D][VDIM];
874  for (int d1 = 0; d1 < D1D; ++d1)
875  {
876  for (int d2 = 0; d2 < D1D; ++d2)
877  {
878  for (int c = 0; c < VDIM; c++)
879  {
880  BBDBu0[d1][d2][c] = 0.0;
881  BBDBu1[d1][d2][c] = 0.0;
882  }
883  for (int q2 = 0; q2 < Q1D; ++q2)
884  {
885  const double b = Bt(d2,q2);
886  for (int c = 0; c < VDIM; c++)
887  {
888  BBDBu0[d1][d2][c] += b*BDBu0[d1][q2][c];
889  BBDBu1[d1][d2][c] += b*BDBu1[d1][q2][c];
890  }
891  }
892  for (int c = 0; c < VDIM; c++)
893  {
894  y(d1,d2,c,0,f) += BBDBu0[d1][d2][c];
895  y(d1,d2,c,1,f) += BBDBu1[d1][d2][c];
896  }
897  }
898  }
899  });
900 }
901 
902 // Optimized PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
903 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
904 void SmemPADGTraceApplyTranspose3D(const int NF,
905  const Array<double> &b,
906  const Array<double> &bt,
907  const Vector &_op,
908  const Vector &_x,
909  Vector &_y,
910  const int d1d = 0,
911  const int q1d = 0)
912 {
913  const int D1D = T_D1D ? T_D1D : d1d;
914  const int Q1D = T_Q1D ? T_Q1D : q1d;
915  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
916  MFEM_VERIFY(D1D <= MAX_D1D, "");
917  MFEM_VERIFY(Q1D <= MAX_Q1D, "");
918  auto B = Reshape(b.Read(), Q1D, D1D);
919  auto Bt = Reshape(bt.Read(), D1D, Q1D);
920  auto op = Reshape(_op.Read(), Q1D, Q1D, 2, 2, NF);
921  auto x = Reshape(_x.Read(), D1D, D1D, 2, NF);
922  auto y = Reshape(_y.ReadWrite(), D1D, D1D, 2, NF);
923 
924  MFEM_FORALL_2D(f, NF, Q1D, Q1D, NBZ,
925  {
926  const int tidz = MFEM_THREAD_ID(z);
927  const int D1D = T_D1D ? T_D1D : d1d;
928  const int Q1D = T_Q1D ? T_Q1D : q1d;
929  // the following variables are evaluated at compile time
930  constexpr int NBZ = T_NBZ ? T_NBZ : 1;
931  constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
932  constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;
933  MFEM_SHARED double u0[NBZ][max_D1D][max_D1D];
934  MFEM_SHARED double u1[NBZ][max_D1D][max_D1D];
935  MFEM_FOREACH_THREAD(d1,x,D1D)
936  {
937  MFEM_FOREACH_THREAD(d2,y,D1D)
938  {
939  u0[tidz][d1][d2] = x(d1,d2,0,f+tidz);
940  u1[tidz][d1][d2] = x(d1,d2,1,f+tidz);
941  }
942  }
943  MFEM_SYNC_THREAD;
944  MFEM_SHARED double Bu0[NBZ][max_Q1D][max_D1D];
945  MFEM_SHARED double Bu1[NBZ][max_Q1D][max_D1D];
946  MFEM_FOREACH_THREAD(q1,x,Q1D)
947  {
948  MFEM_FOREACH_THREAD(d2,y,D1D)
949  {
950  double Bu0_ = 0.0;
951  double Bu1_ = 0.0;
952  for (int d1 = 0; d1 < D1D; ++d1)
953  {
954  const double b = B(q1,d1);
955  Bu0_ += b*u0[tidz][d1][d2];
956  Bu1_ += b*u1[tidz][d1][d2];
957  }
958  Bu0[tidz][q1][d2] = Bu0_;
959  Bu1[tidz][q1][d2] = Bu1_;
960  }
961  }
962  MFEM_SYNC_THREAD;
963  MFEM_SHARED double BBu0[NBZ][max_Q1D][max_Q1D];
964  MFEM_SHARED double BBu1[NBZ][max_Q1D][max_Q1D];
965  MFEM_FOREACH_THREAD(q1,x,Q1D)
966  {
967  MFEM_FOREACH_THREAD(q2,y,Q1D)
968  {
969  double BBu0_ = 0.0;
970  double BBu1_ = 0.0;
971  for (int d2 = 0; d2 < D1D; ++d2)
972  {
973  const double b = B(q2,d2);
974  BBu0_ += b*Bu0[tidz][q1][d2];
975  BBu1_ += b*Bu1[tidz][q1][d2];
976  }
977  BBu0[tidz][q1][q2] = BBu0_;
978  BBu1[tidz][q1][q2] = BBu1_;
979  }
980  }
981  MFEM_SYNC_THREAD;
982  MFEM_SHARED double DBBu0[NBZ][max_Q1D][max_Q1D];
983  MFEM_SHARED double DBBu1[NBZ][max_Q1D][max_Q1D];
984  MFEM_FOREACH_THREAD(q1,x,Q1D)
985  {
986  MFEM_FOREACH_THREAD(q2,y,Q1D)
987  {
988  const double D00 = op(q1,q2,0,0,f+tidz);
989  const double D01 = op(q1,q2,0,1,f+tidz);
990  const double D10 = op(q1,q2,1,0,f+tidz);
991  const double D11 = op(q1,q2,1,1,f+tidz);
992  const double u0 = BBu0[tidz][q1][q2];
993  const double u1 = BBu1[tidz][q1][q2];
994  DBBu0[tidz][q1][q2] = D00*u0 + D01*u1;
995  DBBu1[tidz][q1][q2] = D10*u0 + D11*u1;
996  }
997  }
998  MFEM_SYNC_THREAD;
999  MFEM_SHARED double BDBBu0[NBZ][max_Q1D][max_D1D];
1000  MFEM_SHARED double BDBBu1[NBZ][max_Q1D][max_D1D];
1001  MFEM_FOREACH_THREAD(q1,x,Q1D)
1002  {
1003  MFEM_FOREACH_THREAD(d2,y,D1D)
1004  {
1005  double BDBBu0_ = 0.0;
1006  double BDBBu1_ = 0.0;
1007  for (int q2 = 0; q2 < Q1D; ++q2)
1008  {
1009  const double b = Bt(d2,q2);
1010  BDBBu0_ += b*DBBu0[tidz][q1][q2];
1011  BDBBu1_ += b*DBBu1[tidz][q1][q2];
1012  }
1013  BDBBu0[tidz][q1][d2] = BDBBu0_;
1014  BDBBu1[tidz][q1][d2] = BDBBu1_;
1015  }
1016  }
1017  MFEM_SYNC_THREAD;
1018  MFEM_FOREACH_THREAD(d1,x,D1D)
1019  {
1020  MFEM_FOREACH_THREAD(d2,y,D1D)
1021  {
1022  double BBDBBu0_ = 0.0;
1023  double BBDBBu1_ = 0.0;
1024  for (int q1 = 0; q1 < Q1D; ++q1)
1025  {
1026  const double b = Bt(d1,q1);
1027  BBDBBu0_ += b*BDBBu0[tidz][q1][d2];
1028  BBDBBu1_ += b*BDBBu1[tidz][q1][d2];
1029  }
1030  y(d1,d2,0,f+tidz) += BBDBBu0_;
1031  y(d1,d2,1,f+tidz) += BBDBBu1_;
1032  }
1033  }
1034  });
1035 }
1036 
1037 static void PADGTraceApplyTranspose(const int dim,
1038  const int D1D,
1039  const int Q1D,
1040  const int NF,
1041  const Array<double> &B,
1042  const Array<double> &Bt,
1043  const Vector &op,
1044  const Vector &x,
1045  Vector &y)
1046 {
1047  if (dim == 2)
1048  {
1049  switch ((D1D << 4 ) | Q1D)
1050  {
1051  case 0x22: return PADGTraceApplyTranspose2D<2,2>(NF,B,Bt,op,x,y);
1052  case 0x33: return PADGTraceApplyTranspose2D<3,3>(NF,B,Bt,op,x,y);
1053  case 0x44: return PADGTraceApplyTranspose2D<4,4>(NF,B,Bt,op,x,y);
1054  case 0x55: return PADGTraceApplyTranspose2D<5,5>(NF,B,Bt,op,x,y);
1055  case 0x66: return PADGTraceApplyTranspose2D<6,6>(NF,B,Bt,op,x,y);
1056  case 0x77: return PADGTraceApplyTranspose2D<7,7>(NF,B,Bt,op,x,y);
1057  case 0x88: return PADGTraceApplyTranspose2D<8,8>(NF,B,Bt,op,x,y);
1058  case 0x99: return PADGTraceApplyTranspose2D<9,9>(NF,B,Bt,op,x,y);
1059  default: return PADGTraceApplyTranspose2D(NF,B,Bt,op,x,y,D1D,Q1D);
1060  }
1061  }
1062  else if (dim == 3)
1063  {
1064  switch ((D1D << 4 ) | Q1D)
1065  {
1066  case 0x23: return SmemPADGTraceApplyTranspose3D<2,3>(NF,B,Bt,op,x,y);
1067  case 0x34: return SmemPADGTraceApplyTranspose3D<3,4>(NF,B,Bt,op,x,y);
1068  case 0x45: return SmemPADGTraceApplyTranspose3D<4,5>(NF,B,Bt,op,x,y);
1069  case 0x56: return SmemPADGTraceApplyTranspose3D<5,6>(NF,B,Bt,op,x,y);
1070  case 0x67: return SmemPADGTraceApplyTranspose3D<6,7>(NF,B,Bt,op,x,y);
1071  case 0x78: return SmemPADGTraceApplyTranspose3D<7,8>(NF,B,Bt,op,x,y);
1072  case 0x89: return SmemPADGTraceApplyTranspose3D<8,9>(NF,B,Bt,op,x,y);
1073  default: return PADGTraceApplyTranspose3D(NF,B,Bt,op,x,y,D1D,Q1D);
1074  }
1075  }
1076  MFEM_ABORT("Unknown kernel.");
1077 }
1078 
1079 // PA DGTraceIntegrator Apply kernel
1080 void DGTraceIntegrator::AddMultPA(const Vector &x, Vector &y) const
1081 {
1082  PADGTraceApply(dim, dofs1D, quad1D, nf,
1083  maps->B, maps->Bt,
1084  pa_data, x, y);
1085 }
1086 
1087 void DGTraceIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
1088 {
1089  PADGTraceApplyTranspose(dim, dofs1D, quad1D, nf,
1090  maps->B, maps->Bt,
1091  pa_data, x, y);
1092 }
1093 
1094 } // namespace mfem
geom
const Geometry::Type geom
Definition: ex1.cpp:40
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::FaceType
FaceType
Definition: mesh.hpp:42
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
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::ToLexOrdering
int ToLexOrdering(const int dim, const int face_id, const int size1d, const int index)
Definition: restriction.cpp:955
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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::GetMemoryType
MemoryType GetMemoryType(MemoryClass mc)
Return a suitable MemoryType for a given MemoryClass.
Definition: mem_manager.cpp:49
mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:35
alpha
const double alpha
Definition: ex15.cpp:336
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