1 // Copyright (c) 2010-2021, 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
QuadratureInterpolator(const FiniteElementSpace & fes,const IntegrationRule & ir)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
QuadratureInterpolator(const FiniteElementSpace & fes,const QuadratureSpace & qs)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>
Eval2D(const int NE,const int vdim,const QVectorLayout q_layout,const GeometricFactors * geom,const DofToQuad & maps,const Vector & e_vec,Vector & q_val,Vector & q_der,Vector & q_det,const int eval_flags)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>
Eval3D(const int NE,const int vdim,const QVectorLayout q_layout,const GeometricFactors * geom,const DofToQuad & maps,const Vector & e_vec,Vector & q_val,Vector & q_der,Vector & q_det,const int eval_flags)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
Mult(const Vector & e_vec,unsigned eval_flags,Vector & q_val,Vector & q_der,Vector & q_det) const370 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
MultTranspose(unsigned eval_flags,const Vector & q_val,const Vector & q_der,Vector & e_vec) const598 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
Values(const Vector & e_vec,Vector & q_val) const610 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
Derivatives(const Vector & e_vec,Vector & q_der) const617 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
PhysDerivatives(const Vector & e_vec,Vector & q_der) const624 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
Determinants(const Vector & e_vec,Vector & q_det) const631 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
639