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