MFEM  v4.6.0
Finite element discretization library
quadinterpolator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "quadinterpolator.hpp"
13 #include "qinterp/dispatch.hpp"
14 #include "qspace.hpp"
15 #include "../general/forall.hpp"
16 #include "../linalg/dtensor.hpp"
17 #include "../linalg/kernels.hpp"
18 
19 namespace mfem
20 {
21 
23  const IntegrationRule &ir):
24 
25  fespace(&fes),
26  qspace(nullptr),
27  IntRule(&ir),
28  q_layout(QVectorLayout::byNODES),
29  use_tensor_products(UsesTensorBasis(fes))
30 {
31  d_buffer.UseDevice(true);
32  if (fespace->GetNE() == 0) { return; }
33  const FiniteElement *fe = fespace->GetFE(0);
34  MFEM_VERIFY(dynamic_cast<const ScalarFiniteElement*>(fe) != NULL,
35  "Only scalar finite elements are supported");
36 }
37 
39  const QuadratureSpace &qs):
40 
41  fespace(&fes),
42  qspace(&qs),
43  IntRule(nullptr),
44  q_layout(QVectorLayout::byNODES),
45  use_tensor_products(UsesTensorBasis(fes))
46 {
47  d_buffer.UseDevice(true);
48  if (fespace->GetNE() == 0) { return; }
49  const FiniteElement *fe = fespace->GetFE(0);
50  MFEM_VERIFY(dynamic_cast<const ScalarFiniteElement*>(fe) != NULL,
51  "Only scalar finite elements are supported");
52 }
53 
54 namespace internal
55 {
56 
57 namespace quadrature_interpolator
58 {
59 
60 // Compute kernel for 1D quadrature interpolation:
61 // * non-tensor product version,
62 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
63 // * assumes 'maps.mode == FULL'.
64 static void Eval1D(const int NE,
65  const int vdim,
66  const QVectorLayout q_layout,
67  const GeometricFactors *geom,
68  const DofToQuad &maps,
69  const Vector &e_vec,
70  Vector &q_val,
71  Vector &q_der,
72  Vector &q_det,
73  const int eval_flags)
74 {
75  using QI = QuadratureInterpolator;
76 
77  const int nd = maps.ndof;
78  const int nq = maps.nqpt;
79  MFEM_ASSERT(maps.mode == DofToQuad::FULL, "internal error");
80  MFEM_ASSERT(!geom || geom->mesh->SpaceDimension() == 1, "");
81  MFEM_VERIFY(vdim == 1 || !(eval_flags & QI::DETERMINANTS), "");
82  MFEM_VERIFY(bool(geom) == bool(eval_flags & QI::PHYSICAL_DERIVATIVES),
83  "'geom' must be given (non-null) only when evaluating physical"
84  " derivatives");
85  const auto B = Reshape(maps.B.Read(), nq, nd);
86  const auto G = Reshape(maps.G.Read(), nq, nd);
87  const auto J = Reshape(geom ? geom->J.Read() : nullptr, nq, NE);
88  const auto E = Reshape(e_vec.Read(), nd, vdim, NE);
89  auto val = q_layout == QVectorLayout::byNODES ?
90  Reshape(q_val.Write(), nq, vdim, NE):
91  Reshape(q_val.Write(), vdim, nq, NE);
92  auto der = q_layout == QVectorLayout::byNODES ?
93  Reshape(q_der.Write(), nq, vdim, NE):
94  Reshape(q_der.Write(), vdim, nq, NE);
95  auto det = Reshape(q_det.Write(), nq, NE);
96  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
97  {
98  for (int q = 0; q < nq; ++q)
99  {
100  if (eval_flags & QI::VALUES)
101  {
102  for (int c = 0; c < vdim; c++)
103  {
104  double q_val = 0.0;
105  for (int d = 0; d < nd; ++d)
106  {
107  q_val += B(q,d)*E(d,c,e);
108  }
109  if (q_layout == QVectorLayout::byVDIM) { val(c,q,e) = q_val; }
110  if (q_layout == QVectorLayout::byNODES) { val(q,c,e) = q_val; }
111  }
112  }
113  if ((eval_flags & QI::DERIVATIVES) ||
114  (eval_flags & QI::PHYSICAL_DERIVATIVES) ||
115  (eval_flags & QI::DETERMINANTS))
116  {
117  for (int c = 0; c < vdim; c++)
118  {
119  double q_d = 0.0;
120  for (int d = 0; d < nd; ++d)
121  {
122  q_d += G(q,d)*E(d,c,e);
123  }
124  if (eval_flags & QI::PHYSICAL_DERIVATIVES)
125  {
126  q_d /= J(q,e);
127  }
128  if (eval_flags & QI::DERIVATIVES || eval_flags & QI::PHYSICAL_DERIVATIVES)
129  {
130  if (q_layout == QVectorLayout::byVDIM) { der(c,q,e) = q_d; }
131  if (q_layout == QVectorLayout::byNODES) { der(q,c,e) = q_d; }
132  }
133  if (vdim == 1 && (eval_flags & QI::DETERMINANTS))
134  {
135  det(q,e) = q_d;
136  }
137  }
138  }
139  }
140  });
141 }
142 
143 // Template compute kernel for 2D quadrature interpolation:
144 // * non-tensor product version,
145 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
146 // * assumes 'maps.mode == FULL'.
147 template<const int T_VDIM, const int T_ND, const int T_NQ>
148 static void Eval2D(const int NE,
149  const int vdim,
150  const QVectorLayout q_layout,
151  const GeometricFactors *geom,
152  const DofToQuad &maps,
153  const Vector &e_vec,
154  Vector &q_val,
155  Vector &q_der,
156  Vector &q_det,
157  const int eval_flags)
158 {
159  using QI = QuadratureInterpolator;
160 
161  const int nd = maps.ndof;
162  const int nq = maps.nqpt;
163  const int ND = T_ND ? T_ND : nd;
164  const int NQ = T_NQ ? T_NQ : nq;
165  const int NMAX = NQ > ND ? NQ : ND;
166  const int VDIM = T_VDIM ? T_VDIM : vdim;
167  MFEM_ASSERT(maps.mode == DofToQuad::FULL, "internal error");
168  MFEM_ASSERT(!geom || geom->mesh->SpaceDimension() == 2, "");
169  MFEM_VERIFY(ND <= QI::MAX_ND2D, "");
170  MFEM_VERIFY(NQ <= QI::MAX_NQ2D, "");
171  MFEM_VERIFY(VDIM == 2 || !(eval_flags & QI::DETERMINANTS), "");
172  MFEM_VERIFY(bool(geom) == bool(eval_flags & QI::PHYSICAL_DERIVATIVES),
173  "'geom' must be given (non-null) only when evaluating physical"
174  " derivatives");
175  const auto B = Reshape(maps.B.Read(), NQ, ND);
176  const auto G = Reshape(maps.G.Read(), NQ, 2, ND);
177  const auto J = Reshape(geom ? geom->J.Read() : nullptr, NQ, 2, 2, NE);
178  const auto E = Reshape(e_vec.Read(), ND, VDIM, NE);
179  auto val = q_layout == QVectorLayout::byNODES ?
180  Reshape(q_val.Write(), NQ, VDIM, NE):
181  Reshape(q_val.Write(), VDIM, NQ, NE);
182  auto der = q_layout == QVectorLayout::byNODES ?
183  Reshape(q_der.Write(), NQ, VDIM, 2, NE):
184  Reshape(q_der.Write(), VDIM, 2, NQ, NE);
185  auto det = Reshape(q_det.Write(), NQ, NE);
186  mfem::forall_2D(NE, NMAX, 1, [=] MFEM_HOST_DEVICE (int e)
187  {
188  const int ND = T_ND ? T_ND : nd;
189  const int NQ = T_NQ ? T_NQ : nq;
190  const int VDIM = T_VDIM ? T_VDIM : vdim;
191  constexpr int max_ND = T_ND ? T_ND : QI::MAX_ND2D;
192  constexpr int max_VDIM = T_VDIM ? T_VDIM : QI::MAX_VDIM2D;
193  MFEM_SHARED double s_E[max_VDIM*max_ND];
194  MFEM_FOREACH_THREAD(d, x, ND)
195  {
196  for (int c = 0; c < VDIM; c++)
197  {
198  s_E[c+d*VDIM] = E(d,c,e);
199  }
200  }
201  MFEM_SYNC_THREAD;
202 
203  MFEM_FOREACH_THREAD(q, x, NQ)
204  {
205  if (eval_flags & QI::VALUES)
206  {
207  double ed[max_VDIM];
208  for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
209  for (int d = 0; d < ND; ++d)
210  {
211  const double b = B(q,d);
212  for (int c = 0; c < VDIM; c++) { ed[c] += b*s_E[c+d*VDIM]; }
213  }
214  for (int c = 0; c < VDIM; c++)
215  {
216  if (q_layout == QVectorLayout::byVDIM) { val(c,q,e) = ed[c]; }
217  if (q_layout == QVectorLayout::byNODES) { val(q,c,e) = ed[c]; }
218  }
219  }
220  if ((eval_flags & QI::DERIVATIVES) ||
221  (eval_flags & QI::PHYSICAL_DERIVATIVES) ||
222  (eval_flags & QI::DETERMINANTS))
223  {
224  // use MAX_VDIM2D to avoid "subscript out of range" warnings
225  double D[QI::MAX_VDIM2D*2];
226  for (int i = 0; i < 2*VDIM; i++) { D[i] = 0.0; }
227  for (int d = 0; d < ND; ++d)
228  {
229  const double wx = G(q,0,d);
230  const double wy = G(q,1,d);
231  for (int c = 0; c < VDIM; c++)
232  {
233  double s_e = s_E[c+d*VDIM];
234  D[c+VDIM*0] += s_e * wx;
235  D[c+VDIM*1] += s_e * wy;
236  }
237  }
238  if (eval_flags & QI::DERIVATIVES)
239  {
240  for (int c = 0; c < VDIM; c++)
241  {
242  if (q_layout == QVectorLayout::byVDIM)
243  {
244  der(c,0,q,e) = D[c+VDIM*0];
245  der(c,1,q,e) = D[c+VDIM*1];
246  }
247  if (q_layout == QVectorLayout::byNODES)
248  {
249  der(q,c,0,e) = D[c+VDIM*0];
250  der(q,c,1,e) = D[c+VDIM*1];
251  }
252  }
253  }
254  if (eval_flags & QI::PHYSICAL_DERIVATIVES)
255  {
256  double Jloc[4], Jinv[4];
257  Jloc[0] = J(q,0,0,e);
258  Jloc[1] = J(q,1,0,e);
259  Jloc[2] = J(q,0,1,e);
260  Jloc[3] = J(q,1,1,e);
261  kernels::CalcInverse<2>(Jloc, Jinv);
262  for (int c = 0; c < VDIM; c++)
263  {
264  const double u = D[c+VDIM*0];
265  const double v = D[c+VDIM*1];
266  const double JiU = Jinv[0]*u + Jinv[1]*v;
267  const double JiV = Jinv[2]*u + Jinv[3]*v;
268  if (q_layout == QVectorLayout::byVDIM)
269  {
270  der(c,0,q,e) = JiU;
271  der(c,1,q,e) = JiV;
272  }
273  if (q_layout == QVectorLayout::byNODES)
274  {
275  der(q,c,0,e) = JiU;
276  der(q,c,1,e) = JiV;
277  }
278  }
279  }
280  if (VDIM == 2 && (eval_flags & QI::DETERMINANTS))
281  {
282  // The check (VDIM == 2) should eliminate this block when VDIM is
283  // known at compile time and (VDIM != 2).
284  det(q,e) = kernels::Det<2>(D);
285  }
286  }
287  }
288  });
289 }
290 
291 // Template compute kernel for 3D quadrature interpolation:
292 // * non-tensor product version,
293 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
294 // * assumes 'maps.mode == FULL'.
295 template<const int T_VDIM, const int T_ND, const int T_NQ>
296 static void Eval3D(const int NE,
297  const int vdim,
298  const QVectorLayout q_layout,
299  const GeometricFactors *geom,
300  const DofToQuad &maps,
301  const Vector &e_vec,
302  Vector &q_val,
303  Vector &q_der,
304  Vector &q_det,
305  const int eval_flags)
306 {
307  using QI = QuadratureInterpolator;
308 
309  const int nd = maps.ndof;
310  const int nq = maps.nqpt;
311  const int ND = T_ND ? T_ND : nd;
312  const int NQ = T_NQ ? T_NQ : nq;
313  const int NMAX = NQ > ND ? NQ : ND;
314  const int VDIM = T_VDIM ? T_VDIM : vdim;
315  MFEM_ASSERT(maps.mode == DofToQuad::FULL, "internal error");
316  MFEM_ASSERT(!geom || geom->mesh->SpaceDimension() == 3, "");
317  MFEM_VERIFY(ND <= QI::MAX_ND3D, "");
318  MFEM_VERIFY(NQ <= QI::MAX_NQ3D, "");
319  MFEM_VERIFY(VDIM == 3 || !(eval_flags & QI::DETERMINANTS), "");
320  MFEM_VERIFY(bool(geom) == bool(eval_flags & QI::PHYSICAL_DERIVATIVES),
321  "'geom' must be given (non-null) only when evaluating physical"
322  " derivatives");
323  const auto B = Reshape(maps.B.Read(), NQ, ND);
324  const auto G = Reshape(maps.G.Read(), NQ, 3, ND);
325  const auto J = Reshape(geom ? geom->J.Read() : nullptr, NQ, 3, 3, NE);
326  const auto E = Reshape(e_vec.Read(), ND, VDIM, NE);
327  auto val = q_layout == QVectorLayout::byNODES ?
328  Reshape(q_val.Write(), NQ, VDIM, NE):
329  Reshape(q_val.Write(), VDIM, NQ, NE);
330  auto der = q_layout == QVectorLayout::byNODES ?
331  Reshape(q_der.Write(), NQ, VDIM, 3, NE):
332  Reshape(q_der.Write(), VDIM, 3, NQ, NE);
333  auto det = Reshape(q_det.Write(), NQ, NE);
334  mfem::forall_2D(NE, NMAX, 1, [=] MFEM_HOST_DEVICE (int e)
335  {
336  const int ND = T_ND ? T_ND : nd;
337  const int NQ = T_NQ ? T_NQ : nq;
338  const int VDIM = T_VDIM ? T_VDIM : vdim;
339  constexpr int max_ND = T_ND ? T_ND : QI::MAX_ND3D;
340  constexpr int max_VDIM = T_VDIM ? T_VDIM : QI::MAX_VDIM3D;
341  MFEM_SHARED double s_E[max_VDIM*max_ND];
342  MFEM_FOREACH_THREAD(d, x, ND)
343  {
344  for (int c = 0; c < VDIM; c++)
345  {
346  s_E[c+d*VDIM] = E(d,c,e);
347  }
348  }
349  MFEM_SYNC_THREAD;
350 
351  MFEM_FOREACH_THREAD(q, x, NQ)
352  {
353  if (eval_flags & QI::VALUES)
354  {
355  double ed[max_VDIM];
356  for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
357  for (int d = 0; d < ND; ++d)
358  {
359  const double b = B(q,d);
360  for (int c = 0; c < VDIM; c++) { ed[c] += b*s_E[c+d*VDIM]; }
361  }
362  for (int c = 0; c < VDIM; c++)
363  {
364  if (q_layout == QVectorLayout::byVDIM) { val(c,q,e) = ed[c]; }
365  if (q_layout == QVectorLayout::byNODES) { val(q,c,e) = ed[c]; }
366  }
367  }
368  if ((eval_flags & QI::DERIVATIVES) ||
369  (eval_flags & QI::PHYSICAL_DERIVATIVES) ||
370  (eval_flags & QI::DETERMINANTS))
371  {
372  // use MAX_VDIM3D to avoid "subscript out of range" warnings
373  double D[QI::MAX_VDIM3D*3];
374  for (int i = 0; i < 3*VDIM; i++) { D[i] = 0.0; }
375  for (int d = 0; d < ND; ++d)
376  {
377  const double wx = G(q,0,d);
378  const double wy = G(q,1,d);
379  const double wz = G(q,2,d);
380  for (int c = 0; c < VDIM; c++)
381  {
382  double s_e = s_E[c+d*VDIM];
383  D[c+VDIM*0] += s_e * wx;
384  D[c+VDIM*1] += s_e * wy;
385  D[c+VDIM*2] += s_e * wz;
386  }
387  }
388  if (eval_flags & QI::DERIVATIVES)
389  {
390  for (int c = 0; c < VDIM; c++)
391  {
392  if (q_layout == QVectorLayout::byVDIM)
393  {
394  der(c,0,q,e) = D[c+VDIM*0];
395  der(c,1,q,e) = D[c+VDIM*1];
396  der(c,2,q,e) = D[c+VDIM*2];
397  }
398  if (q_layout == QVectorLayout::byNODES)
399  {
400  der(q,c,0,e) = D[c+VDIM*0];
401  der(q,c,1,e) = D[c+VDIM*1];
402  der(q,c,2,e) = D[c+VDIM*2];
403  }
404  }
405  }
406  if (eval_flags & QI::PHYSICAL_DERIVATIVES)
407  {
408  double Jloc[9], Jinv[9];
409  for (int col = 0; col < 3; col++)
410  {
411  for (int row = 0; row < 3; row++)
412  {
413  Jloc[row+3*col] = J(q,row,col,e);
414  }
415  }
416  kernels::CalcInverse<3>(Jloc, Jinv);
417  for (int c = 0; c < VDIM; c++)
418  {
419  const double u = D[c+VDIM*0];
420  const double v = D[c+VDIM*1];
421  const double w = D[c+VDIM*2];
422  const double JiU = Jinv[0]*u + Jinv[1]*v + Jinv[2]*w;
423  const double JiV = Jinv[3]*u + Jinv[4]*v + Jinv[5]*w;
424  const double JiW = Jinv[6]*u + Jinv[7]*v + Jinv[8]*w;
425  if (q_layout == QVectorLayout::byVDIM)
426  {
427  der(c,0,q,e) = JiU;
428  der(c,1,q,e) = JiV;
429  der(c,2,q,e) = JiW;
430  }
431  if (q_layout == QVectorLayout::byNODES)
432  {
433  der(q,c,0,e) = JiU;
434  der(q,c,1,e) = JiV;
435  der(q,c,2,e) = JiW;
436  }
437  }
438  }
439  if (VDIM == 3 && (eval_flags & QI::DETERMINANTS))
440  {
441  // The check (VDIM == 3) should eliminate this block when VDIM is
442  // known at compile time and (VDIM != 3).
443  det(q,e) = kernels::Det<3>(D);
444  }
445  }
446  }
447  });
448 }
449 
450 } // namespace quadrature_interpolator
451 
452 } // namespace internal
453 
455  unsigned eval_flags,
456  Vector &q_val,
457  Vector &q_der,
458  Vector &q_det) const
459 {
460  using namespace internal::quadrature_interpolator;
461 
462  const int ne = fespace->GetNE();
463  if (ne == 0) { return; }
464  const int vdim = fespace->GetVDim();
465  const FiniteElement *fe = fespace->GetFE(0);
466  const bool use_tensor_eval =
468  dynamic_cast<const TensorBasisElement*>(fe) != nullptr;
469  const IntegrationRule *ir =
471  const DofToQuad::Mode mode =
472  use_tensor_eval ? DofToQuad::TENSOR : DofToQuad::FULL;
473  const DofToQuad &maps = fe->GetDofToQuad(*ir, mode);
474  const GeometricFactors *geom = nullptr;
475  if (eval_flags & PHYSICAL_DERIVATIVES)
476  {
477  const int jacobians = GeometricFactors::JACOBIANS;
478  geom = fespace->GetMesh()->GetGeometricFactors(*ir, jacobians);
479  }
480 
481  MFEM_ASSERT(fespace->GetMesh()->GetNumGeometries(
482  fespace->GetMesh()->Dimension()) == 1,
483  "mixed meshes are not supported");
484 
485  if (use_tensor_eval)
486  {
487  // TODO: use fused kernels
489  {
490  if (eval_flags & VALUES)
491  {
492  TensorValues<QVectorLayout::byNODES>(ne, vdim, maps, e_vec, q_val);
493  }
494  if (eval_flags & DERIVATIVES)
495  {
496  TensorDerivatives<QVectorLayout::byNODES>(
497  ne, vdim, maps, e_vec, q_der);
498  }
499  if (eval_flags & PHYSICAL_DERIVATIVES)
500  {
501  TensorPhysDerivatives<QVectorLayout::byNODES>(
502  ne, vdim, maps, *geom, e_vec, q_der);
503  }
504  }
505 
507  {
508  if (eval_flags & VALUES)
509  {
510  TensorValues<QVectorLayout::byVDIM>(ne, vdim, maps, e_vec, q_val);
511  }
512  if (eval_flags & DERIVATIVES)
513  {
514  TensorDerivatives<QVectorLayout::byVDIM>(
515  ne, vdim, maps, e_vec, q_der);
516  }
517  if (eval_flags & PHYSICAL_DERIVATIVES)
518  {
519  TensorPhysDerivatives<QVectorLayout::byVDIM>(
520  ne, vdim, maps, *geom, e_vec, q_der);
521  }
522  }
523  if (eval_flags & DETERMINANTS)
524  {
525  TensorDeterminants(ne, vdim, maps, e_vec, q_det, d_buffer);
526  }
527  }
528  else // use_tensor_eval == false
529  {
530  const int nd = maps.ndof;
531  const int nq = maps.nqpt;
532  const int dim = maps.FE->GetDim();
533 
534  void (*mult)(const int NE,
535  const int vdim,
536  const QVectorLayout q_layout,
537  const GeometricFactors *geom,
538  const DofToQuad &maps,
539  const Vector &e_vec,
540  Vector &q_val,
541  Vector &q_der,
542  Vector &q_det,
543  const int eval_flags) = NULL;
544 
545  if (dim == 1)
546  {
547  mult = &Eval1D;
548  }
549  else if (vdim == 1) // dim == 2 || dim == 3
550  {
551  if (dim == 2)
552  {
553  switch (100*nd + nq)
554  {
555  // Q0
556  case 101: mult = &Eval2D<1,1,1>; break;
557  case 104: mult = &Eval2D<1,1,4>; break;
558  // Q1
559  case 404: mult = &Eval2D<1,4,4>; break;
560  case 409: mult = &Eval2D<1,4,9>; break;
561  // Q2
562  case 909: mult = &Eval2D<1,9,9>; break;
563  case 916: mult = &Eval2D<1,9,16>; break;
564  // Q3
565  case 1616: mult = &Eval2D<1,16,16>; break;
566  case 1625: mult = &Eval2D<1,16,25>; break;
567  case 1636: mult = &Eval2D<1,16,36>; break;
568  // Q4
569  case 2525: mult = &Eval2D<1,25,25>; break;
570  case 2536: mult = &Eval2D<1,25,36>; break;
571  case 2549: mult = &Eval2D<1,25,49>; break;
572  case 2564: mult = &Eval2D<1,25,64>; break;
573  }
574  if (nq >= 100 || !mult)
575  {
576  mult = &Eval2D<1,0,0>;
577  }
578  }
579  else if (dim == 3)
580  {
581  switch (1000*nd + nq)
582  {
583  // Q0
584  case 1001: mult = &Eval3D<1,1,1>; break;
585  case 1008: mult = &Eval3D<1,1,8>; break;
586  // Q1
587  case 8008: mult = &Eval3D<1,8,8>; break;
588  case 8027: mult = &Eval3D<1,8,27>; break;
589  // Q2
590  case 27027: mult = &Eval3D<1,27,27>; break;
591  case 27064: mult = &Eval3D<1,27,64>; break;
592  // Q3
593  case 64064: mult = &Eval3D<1,64,64>; break;
594  case 64125: mult = &Eval3D<1,64,125>; break;
595  case 64216: mult = &Eval3D<1,64,216>; break;
596  // Q4
597  case 125125: mult = &Eval3D<1,125,125>; break;
598  case 125216: mult = &Eval3D<1,125,216>; break;
599  }
600  if (nq >= 1000 || !mult)
601  {
602  mult = &Eval3D<1,0,0>;
603  }
604  }
605  }
606  else if (vdim == 3 && dim == 2)
607  {
608  switch (100*nd + nq)
609  {
610  // Q0
611  case 101: mult = &Eval2D<3,1,1>; break;
612  case 104: mult = &Eval2D<3,1,4>; break;
613  // Q1
614  case 404: mult = &Eval2D<3,4,4>; break;
615  case 409: mult = &Eval2D<3,4,9>; break;
616  // Q2
617  case 904: mult = &Eval2D<3,9,4>; break;
618  case 909: mult = &Eval2D<3,9,9>; break;
619  case 916: mult = &Eval2D<3,9,16>; break;
620  case 925: mult = &Eval2D<3,9,25>; break;
621  // Q3
622  case 1616: mult = &Eval2D<3,16,16>; break;
623  case 1625: mult = &Eval2D<3,16,25>; break;
624  case 1636: mult = &Eval2D<3,16,36>; break;
625  // Q4
626  case 2525: mult = &Eval2D<3,25,25>; break;
627  case 2536: mult = &Eval2D<3,25,36>; break;
628  case 2549: mult = &Eval2D<3,25,49>; break;
629  case 2564: mult = &Eval2D<3,25,64>; break;
630  default: mult = &Eval2D<3,0,0>;
631  }
632  }
633  else if (vdim == dim)
634  {
635  if (dim == 2)
636  {
637  switch (100*nd + nq)
638  {
639  // Q1
640  case 404: mult = &Eval2D<2,4,4>; break;
641  case 409: mult = &Eval2D<2,4,9>; break;
642  // Q2
643  case 909: mult = &Eval2D<2,9,9>; break;
644  case 916: mult = &Eval2D<2,9,16>; break;
645  // Q3
646  case 1616: mult = &Eval2D<2,16,16>; break;
647  case 1625: mult = &Eval2D<2,16,25>; break;
648  case 1636: mult = &Eval2D<2,16,36>; break;
649  // Q4
650  case 2525: mult = &Eval2D<2,25,25>; break;
651  case 2536: mult = &Eval2D<2,25,36>; break;
652  case 2549: mult = &Eval2D<2,25,49>; break;
653  case 2564: mult = &Eval2D<2,25,64>; break;
654  }
655  if (nq >= 100 || !mult) { mult = &Eval2D<2,0,0>; }
656  }
657  else if (dim == 3)
658  {
659  switch (1000*nd + nq)
660  {
661  // Q1
662  case 8008: mult = &Eval3D<3,8,8>; break;
663  case 8027: mult = &Eval3D<3,8,27>; break;
664  // Q2
665  case 27027: mult = &Eval3D<3,27,27>; break;
666  case 27064: mult = &Eval3D<3,27,64>; break;
667  case 27125: mult = &Eval3D<3,27,125>; break;
668  // Q3
669  case 64064: mult = &Eval3D<3,64,64>; break;
670  case 64125: mult = &Eval3D<3,64,125>; break;
671  case 64216: mult = &Eval3D<3,64,216>; break;
672  // Q4
673  case 125125: mult = &Eval3D<3,125,125>; break;
674  case 125216: mult = &Eval3D<3,125,216>; break;
675  }
676  if (nq >= 1000 || !mult) { mult = &Eval3D<3,0,0>; }
677  }
678  }
679  if (mult)
680  {
681  mult(ne,vdim,q_layout,geom,maps,e_vec,q_val,q_der,q_det,eval_flags);
682  }
683  else { MFEM_ABORT("case not supported yet"); }
684  }
685 }
686 
687 void QuadratureInterpolator::MultTranspose(unsigned eval_flags,
688  const Vector &q_val,
689  const Vector &q_der,
690  Vector &e_vec) const
691 {
692  MFEM_CONTRACT_VAR(eval_flags);
693  MFEM_CONTRACT_VAR(q_val);
694  MFEM_CONTRACT_VAR(q_der);
695  MFEM_CONTRACT_VAR(e_vec);
696  MFEM_ABORT("this method is not implemented yet");
697 }
698 
700  Vector &q_val) const
701 {
702  Vector empty;
703  Mult(e_vec, VALUES, q_val, empty, empty);
704 }
705 
707  Vector &q_der) const
708 {
709  Vector empty;
710  Mult(e_vec, DERIVATIVES, empty, q_der, empty);
711 }
712 
714  Vector &q_der) const
715 {
716  Vector empty;
717  Mult(e_vec, PHYSICAL_DERIVATIVES, empty, q_der, empty);
718 }
719 
721  Vector &q_det) const
722 {
723  Vector empty;
724  Mult(e_vec, DETERMINANTS, empty, empty, q_det);
725 }
726 
727 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
Abstract class for all finite elements.
Definition: fe_base.hpp:233
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition: fe_base.hpp:141
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
void forall_2D(int N, int X, int Y, lambda &&body)
Definition: forall.hpp:751
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:161
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:847
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
const FiniteElementSpace * fespace
Not owned.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
const IntegrationRule & GetElementIntRule(int idx) const
Get the IntegrationRule associated with mesh element idx.
Definition: qspace.hpp:130
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:173
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:169
Evaluate the derivatives at quadrature points.
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition: mesh.hpp:2183
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition: fespace.hpp:1306
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:6274
void Determinants(const Vector &e_vec, Vector &q_det) const
Compute the determinants of the derivatives (with respect to reference coordinates) of the E-vector e...
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:365
Vector J
Jacobians of the element transformations at all quadrature points.
Definition: mesh.hpp:2224
QVectorLayout q_layout
Output Q-vector layout.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
double b
Definition: lissajous.cpp:42
const IntegrationRule * IntRule
Not owned.
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition: fe_base.hpp:165
bool use_tensor_products
Tensor product evaluation mode.
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
Vector d_buffer
Auxiliary device buffer.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
QuadratureInterpolator(const FiniteElementSpace &fes, const IntegrationRule &ir)
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
void Values(const Vector &e_vec, Vector &q_val) const
Interpolate the values of the E-vector e_vec at quadrature points.
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:136
void MultTranspose(unsigned eval_flags, const Vector &q_val, const Vector &q_der, Vector &e_vec) const
Perform the transpose operation of Mult(). (TODO)
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:184
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition: fespace.hpp:52
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition: fe_base.hpp:149
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition: fe_base.hpp:154
void Derivatives(const Vector &e_vec, Vector &q_der) const
Interpolate the derivatives (with respect to reference coordinates) of the E-vector e_vec at quadratu...
int dim
Definition: ex24.cpp:53
void Mult(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_der, Vector &q_det) const
Interpolate the E-vector e_vec to quadrature points.
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition: fe_base.hpp:205
void PhysDerivatives(const Vector &e_vec, Vector &q_der) const
Interpolate the derivatives in physical space of the E-vector e_vec at quadrature points...
const Mesh * mesh
Definition: mesh.hpp:2191
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
Vector data type.
Definition: vector.hpp:58
Class representing the storage layout of a QuadratureFunction.
Definition: qspace.hpp:101
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
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
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
Evaluate the values at quadrature points.
const QuadratureSpace * qspace
Not owned.