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