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