MFEM v4.9.0
Finite element discretization library
Loading...
Searching...
No Matches
quadinterpolator.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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/grad.hpp"
14#include "qinterp/eval.hpp"
15#include "qspace.hpp"
16#include "../general/forall.hpp"
17#include "../linalg/dtensor.hpp"
18#include "../linalg/kernels.hpp"
19
20namespace mfem
21{
22
23namespace internal
24{
25namespace quadrature_interpolator
26{
27void InitEvalByNodesKernels();
28void InitEvalByVDimKernels();
29void InitEvalKernels();
30void InitDetKernels();
31template <bool P> void InitGradByNodesKernels();
32template <bool P> void InitGradByVDimKernels();
33void InitTensorEvalHDivKernels();
34struct Kernels
35{
36 Kernels()
37 {
38 using namespace internal::quadrature_interpolator;
39
40 InitEvalByNodesKernels();
41 InitEvalByVDimKernels();
42 // Non-phys grad kernels
43 InitGradByNodesKernels<false>();
44 InitGradByVDimKernels<false>();
45 // Phys grad kernels
46 InitGradByNodesKernels<true>();
47 InitGradByVDimKernels<true>();
48 // Determinants
49 InitDetKernels();
50 // Non-tensor
51 InitEvalKernels();
52 // Tensor (quad,hex) H(div)
53 InitTensorEvalHDivKernels();
54 }
55};
56}
57}
58
60 const IntegrationRule &ir):
61
62 fespace(&fes),
63 qspace(nullptr),
64 IntRule(&ir),
65 q_layout(QVectorLayout::byNODES),
66 use_tensor_products(UsesTensorBasis(fes))
67{
68 static internal::quadrature_interpolator::Kernels kernels;
69
70 d_buffer.UseDevice(true);
71 if (fespace->GetNE() == 0) { return; }
72 MFEM_VERIFY(SupportsFESpace(fes),
73 "Only elements with MapType VALUE and H_DIV are supported!");
74}
75
77 const QuadratureSpace &qs):
78
79 fespace(&fes),
80 qspace(&qs),
81 IntRule(nullptr),
82 q_layout(QVectorLayout::byNODES),
83 use_tensor_products(UsesTensorBasis(fes))
84{
85 d_buffer.UseDevice(true);
86 if (fespace->GetNE() == 0) { return; }
87 MFEM_VERIFY(SupportsFESpace(fes),
88 "Only elements with MapType VALUE and H_DIV are supported!");
89}
90
100
101namespace internal
102{
103
104namespace quadrature_interpolator
105{
106
107// Compute kernel for 1D quadrature interpolation:
108// * non-tensor product version,
109// * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
110// * assumes 'maps.mode == FULL'.
111static void Eval1D(const int NE,
112 const int vdim,
113 const QVectorLayout q_layout,
114 const GeometricFactors *geom,
115 const DofToQuad &maps,
116 const Vector &e_vec,
117 Vector &q_val,
118 Vector &q_der,
119 Vector &q_det,
120 const int eval_flags)
121{
122 using QI = QuadratureInterpolator;
123
124 const int nd = maps.ndof;
125 const int nq = maps.nqpt;
126 MFEM_ASSERT(maps.mode == DofToQuad::FULL, "internal error");
127 MFEM_ASSERT(!geom || geom->mesh->SpaceDimension() == 1, "");
128 MFEM_VERIFY(vdim == 1 || !(eval_flags & QI::DETERMINANTS), "");
129 MFEM_VERIFY(bool(geom) == bool(eval_flags & QI::PHYSICAL_DERIVATIVES),
130 "'geom' must be given (non-null) only when evaluating physical"
131 " derivatives");
132 const auto B = Reshape(maps.B.Read(), nq, nd);
133 const auto G = Reshape(maps.G.Read(), nq, nd);
134 const auto J = Reshape(geom ? geom->J.Read() : nullptr, nq, NE);
135 const auto E = Reshape(e_vec.Read(), nd, vdim, NE);
136 auto val = q_layout == QVectorLayout::byNODES ?
137 Reshape(q_val.Write(), nq, vdim, NE):
138 Reshape(q_val.Write(), vdim, nq, NE);
139 auto der = q_layout == QVectorLayout::byNODES ?
140 Reshape(q_der.Write(), nq, vdim, NE):
141 Reshape(q_der.Write(), vdim, nq, NE);
142 auto det = Reshape(q_det.Write(), nq, NE);
143 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
144 {
145 for (int q = 0; q < nq; ++q)
146 {
147 if (eval_flags & (QI::VALUES | QI::PHYSICAL_VALUES))
148 {
149 for (int c = 0; c < vdim; c++)
150 {
151 real_t q_val = 0.0;
152 for (int d = 0; d < nd; ++d)
153 {
154 q_val += B(q,d)*E(d,c,e);
155 }
156 if (q_layout == QVectorLayout::byVDIM) { val(c,q,e) = q_val; }
157 if (q_layout == QVectorLayout::byNODES) { val(q,c,e) = q_val; }
158 }
159 }
160 if ((eval_flags & QI::DERIVATIVES) ||
161 (eval_flags & QI::PHYSICAL_DERIVATIVES) ||
162 (eval_flags & QI::DETERMINANTS))
163 {
164 for (int c = 0; c < vdim; c++)
165 {
166 real_t q_d = 0.0;
167 for (int d = 0; d < nd; ++d)
168 {
169 q_d += G(q,d)*E(d,c,e);
170 }
171 if (eval_flags & QI::PHYSICAL_DERIVATIVES)
172 {
173 q_d /= J(q,e);
174 }
175 if (eval_flags & QI::DERIVATIVES || eval_flags & QI::PHYSICAL_DERIVATIVES)
176 {
177 if (q_layout == QVectorLayout::byVDIM) { der(c,q,e) = q_d; }
178 if (q_layout == QVectorLayout::byNODES) { der(q,c,e) = q_d; }
179 }
180 if (vdim == 1 && (eval_flags & QI::DETERMINANTS))
181 {
182 det(q,e) = q_d;
183 }
184 }
185 }
186 }
187 });
188}
189
190// Template compute kernel for 2D quadrature interpolation:
191// * non-tensor product version,
192// * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
193// * assumes 'maps.mode == FULL'.
194template<const int T_VDIM, const int T_ND, const int T_NQ>
195static void Eval2D(const int NE,
196 const int vdim,
197 const QVectorLayout q_layout,
198 const GeometricFactors *geom,
199 const DofToQuad &maps,
200 const Vector &e_vec,
201 Vector &q_val,
202 Vector &q_der,
203 Vector &q_det,
204 const int eval_flags)
205{
206 using QI = QuadratureInterpolator;
207
208 const int nd = maps.ndof;
209 const int nq = maps.nqpt;
210 const int ND = T_ND ? T_ND : nd;
211 const int NQ = T_NQ ? T_NQ : nq;
212 const int NMAX = NQ > ND ? NQ : ND;
213 const int VDIM = T_VDIM ? T_VDIM : vdim;
214 MFEM_ASSERT(maps.mode == DofToQuad::FULL, "internal error");
215 MFEM_ASSERT(!geom || geom->mesh->SpaceDimension() == 2, "");
216 MFEM_VERIFY(ND <= QI::MAX_ND2D, "");
217 MFEM_VERIFY(NQ <= QI::MAX_NQ2D, "");
218 MFEM_VERIFY(bool(geom) == bool(eval_flags & QI::PHYSICAL_DERIVATIVES),
219 "'geom' must be given (non-null) only when evaluating physical"
220 " derivatives");
221 const auto B = Reshape(maps.B.Read(), NQ, ND);
222 const auto G = Reshape(maps.G.Read(), NQ, 2, ND);
223 const auto J = Reshape(geom ? geom->J.Read() : nullptr, NQ, 2, 2, NE);
224 const auto E = Reshape(e_vec.Read(), ND, VDIM, NE);
225 auto val = q_layout == QVectorLayout::byNODES ?
226 Reshape(q_val.Write(), NQ, VDIM, NE):
227 Reshape(q_val.Write(), VDIM, NQ, NE);
228 auto der = q_layout == QVectorLayout::byNODES ?
229 Reshape(q_der.Write(), NQ, VDIM, 2, NE):
230 Reshape(q_der.Write(), VDIM, 2, NQ, NE);
231 auto det = Reshape(q_det.Write(), NQ, NE);
232 mfem::forall_2D(NE, NMAX, 1, [=] MFEM_HOST_DEVICE (int e)
233 {
234 const int ND = T_ND ? T_ND : nd;
235 const int NQ = T_NQ ? T_NQ : nq;
236 const int VDIM = T_VDIM ? T_VDIM : vdim;
237 constexpr int max_ND = T_ND ? T_ND : QI::MAX_ND2D;
238 constexpr int max_VDIM = T_VDIM ? T_VDIM : QI::MAX_VDIM2D;
239 MFEM_SHARED real_t s_E[max_VDIM*max_ND];
240 MFEM_FOREACH_THREAD(d, x, ND)
241 {
242 for (int c = 0; c < VDIM; c++)
243 {
244 s_E[c+d*VDIM] = E(d,c,e);
245 }
246 }
247 MFEM_SYNC_THREAD;
248
249 MFEM_FOREACH_THREAD(q, x, NQ)
250 {
251 if (eval_flags & (QI::VALUES | QI::PHYSICAL_VALUES))
252 {
253 real_t ed[max_VDIM];
254 for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
255 for (int d = 0; d < ND; ++d)
256 {
257 const real_t b = B(q,d);
258 for (int c = 0; c < VDIM; c++) { ed[c] += b*s_E[c+d*VDIM]; }
259 }
260 for (int c = 0; c < VDIM; c++)
261 {
262 if (q_layout == QVectorLayout::byVDIM) { val(c,q,e) = ed[c]; }
263 if (q_layout == QVectorLayout::byNODES) { val(q,c,e) = ed[c]; }
264 }
265 }
266 if ((eval_flags & QI::DERIVATIVES) ||
267 (eval_flags & QI::PHYSICAL_DERIVATIVES) ||
268 (eval_flags & QI::DETERMINANTS))
269 {
270 // use MAX_VDIM2D to avoid "subscript out of range" warnings
271 real_t D[QI::MAX_VDIM2D*2];
272 for (int i = 0; i < 2*VDIM; i++) { D[i] = 0.0; }
273 for (int d = 0; d < ND; ++d)
274 {
275 const real_t wx = G(q,0,d);
276 const real_t wy = G(q,1,d);
277 for (int c = 0; c < VDIM; c++)
278 {
279 real_t s_e = s_E[c+d*VDIM];
280 D[c+VDIM*0] += s_e * wx;
281 D[c+VDIM*1] += s_e * wy;
282 }
283 }
284 if (eval_flags & QI::DERIVATIVES)
285 {
286 for (int c = 0; c < VDIM; c++)
287 {
288 if (q_layout == QVectorLayout::byVDIM)
289 {
290 der(c,0,q,e) = D[c+VDIM*0];
291 der(c,1,q,e) = D[c+VDIM*1];
292 }
293 if (q_layout == QVectorLayout::byNODES)
294 {
295 der(q,c,0,e) = D[c+VDIM*0];
296 der(q,c,1,e) = D[c+VDIM*1];
297 }
298 }
299 }
300 if (eval_flags & QI::PHYSICAL_DERIVATIVES)
301 {
302 real_t Jloc[4], Jinv[4];
303 Jloc[0] = J(q,0,0,e);
304 Jloc[1] = J(q,1,0,e);
305 Jloc[2] = J(q,0,1,e);
306 Jloc[3] = J(q,1,1,e);
307 kernels::CalcInverse<2>(Jloc, Jinv);
308 for (int c = 0; c < VDIM; c++)
309 {
310 const real_t u = D[c+VDIM*0];
311 const real_t v = D[c+VDIM*1];
312 const real_t JiU = Jinv[0]*u + Jinv[1]*v;
313 const real_t JiV = Jinv[2]*u + Jinv[3]*v;
314 if (q_layout == QVectorLayout::byVDIM)
315 {
316 der(c,0,q,e) = JiU;
317 der(c,1,q,e) = JiV;
318 }
319 if (q_layout == QVectorLayout::byNODES)
320 {
321 der(q,c,0,e) = JiU;
322 der(q,c,1,e) = JiV;
323 }
324 }
325 }
326 if (eval_flags & QI::DETERMINANTS)
327 {
328 if (VDIM == 2) { det(q,e) = kernels::Det<2>(D); }
329 else
330 {
331 DeviceTensor<2> j(D, 3, 2);
332 const double E = j(0,0)*j(0,0) + j(1,0)*j(1,0) + j(2,0)*j(2,0);
333 const double F = j(0,0)*j(0,1) + j(1,0)*j(1,1) + j(2,0)*j(2,1);
334 const double G = j(0,1)*j(0,1) + j(1,1)*j(1,1) + j(2,1)*j(2,1);
335 det(q,e) = std::sqrt(E*G - F*F);
336 }
337 }
338 }
339 }
340 });
341}
342
343// Template compute kernel for 3D quadrature interpolation:
344// * non-tensor product version,
345// * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
346// * assumes 'maps.mode == FULL'.
347template<const int T_VDIM, const int T_ND, const int T_NQ>
348static void Eval3D(const int NE,
349 const int vdim,
350 const QVectorLayout q_layout,
351 const GeometricFactors *geom,
352 const DofToQuad &maps,
353 const Vector &e_vec,
354 Vector &q_val,
355 Vector &q_der,
356 Vector &q_det,
357 const int eval_flags)
358{
359 using QI = QuadratureInterpolator;
360
361 const int nd = maps.ndof;
362 const int nq = maps.nqpt;
363 const int ND = T_ND ? T_ND : nd;
364 const int NQ = T_NQ ? T_NQ : nq;
365 const int NMAX = NQ > ND ? NQ : ND;
366 const int VDIM = T_VDIM ? T_VDIM : vdim;
367 MFEM_ASSERT(maps.mode == DofToQuad::FULL, "internal error");
368 MFEM_ASSERT(!geom || geom->mesh->SpaceDimension() == 3, "");
369 MFEM_VERIFY(ND <= QI::MAX_ND3D, "");
370 MFEM_VERIFY(NQ <= QI::MAX_NQ3D, "");
371 MFEM_VERIFY(VDIM == 3 || !(eval_flags & QI::DETERMINANTS), "");
372 MFEM_VERIFY(bool(geom) == bool(eval_flags & QI::PHYSICAL_DERIVATIVES),
373 "'geom' must be given (non-null) only when evaluating physical"
374 " derivatives");
375 const auto B = Reshape(maps.B.Read(), NQ, ND);
376 const auto G = Reshape(maps.G.Read(), NQ, 3, ND);
377 const auto J = Reshape(geom ? geom->J.Read() : nullptr, NQ, 3, 3, NE);
378 const auto E = Reshape(e_vec.Read(), ND, VDIM, NE);
379 auto val = q_layout == QVectorLayout::byNODES ?
380 Reshape(q_val.Write(), NQ, VDIM, NE):
381 Reshape(q_val.Write(), VDIM, NQ, NE);
382 auto der = q_layout == QVectorLayout::byNODES ?
383 Reshape(q_der.Write(), NQ, VDIM, 3, NE):
384 Reshape(q_der.Write(), VDIM, 3, NQ, NE);
385 auto det = Reshape(q_det.Write(), NQ, NE);
386 mfem::forall_2D(NE, NMAX, 1, [=] MFEM_HOST_DEVICE (int e)
387 {
388 const int ND = T_ND ? T_ND : nd;
389 const int NQ = T_NQ ? T_NQ : nq;
390 const int VDIM = T_VDIM ? T_VDIM : vdim;
391 constexpr int max_ND = T_ND ? T_ND : QI::MAX_ND3D;
392 constexpr int max_VDIM = T_VDIM ? T_VDIM : QI::MAX_VDIM3D;
393 MFEM_SHARED real_t s_E[max_VDIM*max_ND];
394 MFEM_FOREACH_THREAD(d, x, ND)
395 {
396 for (int c = 0; c < VDIM; c++)
397 {
398 s_E[c+d*VDIM] = E(d,c,e);
399 }
400 }
401 MFEM_SYNC_THREAD;
402
403 MFEM_FOREACH_THREAD(q, x, NQ)
404 {
405 if (eval_flags & (QI::VALUES | QI::PHYSICAL_VALUES))
406 {
407 real_t ed[max_VDIM];
408 for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
409 for (int d = 0; d < ND; ++d)
410 {
411 const real_t b = B(q,d);
412 for (int c = 0; c < VDIM; c++) { ed[c] += b*s_E[c+d*VDIM]; }
413 }
414 for (int c = 0; c < VDIM; c++)
415 {
416 if (q_layout == QVectorLayout::byVDIM) { val(c,q,e) = ed[c]; }
417 if (q_layout == QVectorLayout::byNODES) { val(q,c,e) = ed[c]; }
418 }
419 }
420 if ((eval_flags & QI::DERIVATIVES) ||
421 (eval_flags & QI::PHYSICAL_DERIVATIVES) ||
422 (eval_flags & QI::DETERMINANTS))
423 {
424 // use MAX_VDIM3D to avoid "subscript out of range" warnings
425 real_t D[QI::MAX_VDIM3D*3];
426 for (int i = 0; i < 3*VDIM; i++) { D[i] = 0.0; }
427 for (int d = 0; d < ND; ++d)
428 {
429 const real_t wx = G(q,0,d);
430 const real_t wy = G(q,1,d);
431 const real_t wz = G(q,2,d);
432 for (int c = 0; c < VDIM; c++)
433 {
434 real_t s_e = s_E[c+d*VDIM];
435 D[c+VDIM*0] += s_e * wx;
436 D[c+VDIM*1] += s_e * wy;
437 D[c+VDIM*2] += s_e * wz;
438 }
439 }
440 if (eval_flags & QI::DERIVATIVES)
441 {
442 for (int c = 0; c < VDIM; c++)
443 {
444 if (q_layout == QVectorLayout::byVDIM)
445 {
446 der(c,0,q,e) = D[c+VDIM*0];
447 der(c,1,q,e) = D[c+VDIM*1];
448 der(c,2,q,e) = D[c+VDIM*2];
449 }
450 if (q_layout == QVectorLayout::byNODES)
451 {
452 der(q,c,0,e) = D[c+VDIM*0];
453 der(q,c,1,e) = D[c+VDIM*1];
454 der(q,c,2,e) = D[c+VDIM*2];
455 }
456 }
457 }
458 if (eval_flags & QI::PHYSICAL_DERIVATIVES)
459 {
460 real_t Jloc[9], Jinv[9];
461 for (int col = 0; col < 3; col++)
462 {
463 for (int row = 0; row < 3; row++)
464 {
465 Jloc[row+3*col] = J(q,row,col,e);
466 }
467 }
468 kernels::CalcInverse<3>(Jloc, Jinv);
469 for (int c = 0; c < VDIM; c++)
470 {
471 const real_t u = D[c+VDIM*0];
472 const real_t v = D[c+VDIM*1];
473 const real_t w = D[c+VDIM*2];
474 const real_t JiU = Jinv[0]*u + Jinv[1]*v + Jinv[2]*w;
475 const real_t JiV = Jinv[3]*u + Jinv[4]*v + Jinv[5]*w;
476 const real_t JiW = Jinv[6]*u + Jinv[7]*v + Jinv[8]*w;
477 if (q_layout == QVectorLayout::byVDIM)
478 {
479 der(c,0,q,e) = JiU;
480 der(c,1,q,e) = JiV;
481 der(c,2,q,e) = JiW;
482 }
483 if (q_layout == QVectorLayout::byNODES)
484 {
485 der(q,c,0,e) = JiU;
486 der(q,c,1,e) = JiV;
487 der(q,c,2,e) = JiW;
488 }
489 }
490 }
491 if (VDIM == 3 && (eval_flags & QI::DETERMINANTS))
492 {
493 // The check (VDIM == 3) should eliminate this block when VDIM is
494 // known at compile time and (VDIM != 3).
495 det(q,e) = kernels::Det<3>(D);
496 }
497 }
498 }
499 });
500}
501
502} // namespace quadrature_interpolator
503
504} // namespace internal
505
507 unsigned eval_flags,
508 Vector &q_val,
509 Vector &q_der,
510 Vector &q_det) const
511{
512 using namespace internal::quadrature_interpolator;
513
514 const int ne = fespace->GetNE();
515 if (ne == 0) { return; }
516 const FiniteElement *fe = fespace->GetFE(0);
517
519 {
520 // q_der == q_div
521 return MultHDiv(e_vec, eval_flags, q_val, q_der);
522 }
523
524 const int vdim = fespace->GetVDim();
525 const int sdim = fespace->GetMesh()->SpaceDimension();
526 const bool use_tensor_eval =
528 dynamic_cast<const TensorBasisElement*>(fe) != nullptr;
529 const IntegrationRule *ir =
531 const DofToQuad::Mode mode =
532 use_tensor_eval ? DofToQuad::TENSOR : DofToQuad::FULL;
533 const DofToQuad &maps = fe->GetDofToQuad(*ir, mode);
534 const int dim = maps.FE->GetDim();
535 const int nd = maps.ndof;
536 const int nq = maps.nqpt;
537 const GeometricFactors *geom = nullptr;
538 if (eval_flags & PHYSICAL_DERIVATIVES)
539 {
540 const int jacobians = GeometricFactors::JACOBIANS;
541 geom = fespace->GetMesh()->GetGeometricFactors(*ir, jacobians);
542 }
543
544 MFEM_ASSERT(!(eval_flags & DETERMINANTS) || dim == vdim ||
545 (dim == 2 && vdim == 3), "Invalid dimensions for determinants.");
546 MFEM_ASSERT(fespace->GetMesh()->GetNumGeometries(
547 fespace->GetMesh()->Dimension()) == 1,
548 "mixed meshes are not supported");
549
550 if (use_tensor_eval)
551 {
552 if (eval_flags & (VALUES | PHYSICAL_VALUES))
553 {
554 TensorEvalKernels::Run(dim, q_layout, vdim, nd, nq, ne, maps.B.Read(),
555 e_vec.Read(), q_val.Write(), vdim, nd, nq);
556 }
557 if (eval_flags & (DERIVATIVES | PHYSICAL_DERIVATIVES))
558 {
559 const bool phys = (eval_flags & PHYSICAL_DERIVATIVES);
560 const real_t *J = phys ? geom->J.Read() : nullptr;
561 const int s_dim = phys ? sdim : dim;
562 GradKernels::Run(dim, q_layout, phys, vdim, nd, nq, ne,
563 maps.B.Read(), maps.G.Read(), J, e_vec.Read(),
564 q_der.Write(), s_dim, vdim, nd, nq);
565 }
566 if (eval_flags & DETERMINANTS)
567 {
568 DetKernels::Run(dim, vdim, nd, nq, ne, maps.B.Read(),
569 maps.G.Read(), e_vec.Read(), q_det.Write(), nd,
570 nq, &d_buffer);
571 }
572 }
573 else // use_tensor_eval == false
574 {
575 EvalKernels::Run(dim, vdim, maps.ndof, maps.nqpt, ne,vdim, q_layout,
576 geom, maps, e_vec, q_val, q_der, q_det, eval_flags);
577 }
578}
579
581 unsigned eval_flags,
582 Vector &q_val,
583 Vector &q_div) const
584{
585 const int ne = fespace->GetNE();
586 if (ne == 0) { return; }
587 MFEM_VERIFY(fespace->IsVariableOrder() == false,
588 "variable order spaces are not supported yet!");
589 const FiniteElement *fe = fespace->GetFE(0);
590 MFEM_VERIFY(fe->GetMapType() == FiniteElement::MapType::H_DIV,
591 "this method can be used only for H(div) spaces");
592 MFEM_VERIFY((eval_flags &
594 "only VALUES, PHYSICAL_VALUES, and PHYSICAL_MAGNITUDES"
595 " evaluations are implemented!");
596 const int dim = fespace->GetMesh()->Dimension();
597 const int sdim = fespace->GetMesh()->SpaceDimension();
598 MFEM_VERIFY((dim == 2 || dim == 3) && dim == sdim,
599 "dim = " << dim << ", sdim = " << sdim
600 << " is not supported yet!");
601 MFEM_VERIFY(fespace->GetMesh()->GetNumGeometries(dim) <= 1,
602 "mixed meshes are not supported yet!");
603 const int vdim = fespace->GetVDim();
604 MFEM_VERIFY(vdim == 1, "vdim != 1 is not supported yet!");
605 auto tfe = dynamic_cast<const VectorTensorFiniteElement *>(fe);
606 MFEM_VERIFY(tfe != nullptr, "only quad and hex elements are supported!");
607 MFEM_VERIFY(use_tensor_products,
608 "non-tensor-product evaluation are not supported yet!");
609 const bool use_tensor_eval = use_tensor_products && (tfe != nullptr);
610 const IntegrationRule *ir =
612 const DofToQuad::Mode mode =
613 use_tensor_eval ? DofToQuad::TENSOR : DofToQuad::FULL;
614 const DofToQuad &maps_c = tfe->GetDofToQuad(*ir, mode);
615 const DofToQuad &maps_o = tfe->GetDofToQuadOpen(*ir, mode);
616 const int nd = maps_c.ndof;
617 const int nq = maps_c.nqpt;
618 const GeometricFactors *geom = nullptr;
619 if (eval_flags & (PHYSICAL_VALUES | PHYSICAL_MAGNITUDES))
620 {
621 const int jacobians = GeometricFactors::JACOBIANS;
622 geom = fespace->GetMesh()->GetGeometricFactors(*ir, jacobians);
623 }
624 // Check that at most one of VALUES, PHYSICAL_VALUES, and PHYSICAL_MAGNITUDES
625 // is specified:
626 MFEM_VERIFY(bool(eval_flags & VALUES) + bool(eval_flags & PHYSICAL_VALUES) +
627 bool(eval_flags & PHYSICAL_MAGNITUDES) <= 1,
628 "only one of VALUES, PHYSICAL_VALUES, and PHYSICAL_MAGNITUDES"
629 " can be requested at a time!");
630 const unsigned value_eval_mode =
631 eval_flags & (VALUES | PHYSICAL_VALUES | PHYSICAL_MAGNITUDES);
632 if (value_eval_mode)
633 {
634 // For PHYSICAL_MAGNITUDES the QVectorLayouts are the same and we
635 // instantiate only QVectorLayout::byNODES:
636 const auto q_l = (eval_flags & PHYSICAL_MAGNITUDES) ?
638 TensorEvalHDivKernels::Run(
639 // dispatch params: dim + the template params of EvalHDiv2D/3D:
640 dim, q_l, value_eval_mode, nd, nq,
641 // runtime params, see the arguments of EvalHDiv2D/3D:
642 ne, maps_o.B.Read(), maps_c.B.Read(), geom ? geom->J.Read() : nullptr,
643 e_vec.Read(), q_val.Write(), nd, nq);
644 }
645 MFEM_CONTRACT_VAR(q_div);
646}
647
649 const Vector &q_val,
650 const Vector &q_der,
651 Vector &e_vec) const
652{
653 MFEM_CONTRACT_VAR(eval_flags);
654 MFEM_CONTRACT_VAR(q_val);
655 MFEM_CONTRACT_VAR(q_der);
656 MFEM_CONTRACT_VAR(e_vec);
657 MFEM_ABORT("this method is not implemented yet");
658}
659
661 Vector &q_val) const
662{
663 Vector empty;
664 Mult(e_vec, VALUES, q_val, empty, empty);
665}
666
668 Vector &q_val) const
669{
670 Vector empty;
671 Mult(e_vec, PHYSICAL_VALUES, q_val, empty, empty);
672}
673
675 Vector &q_der) const
676{
677 Vector empty;
678 Mult(e_vec, DERIVATIVES, empty, q_der, empty);
679}
680
682 Vector &q_der) const
683{
684 Vector empty;
685 Mult(e_vec, PHYSICAL_DERIVATIVES, empty, q_der, empty);
686}
687
689 Vector &q_det) const
690{
691 Vector empty;
692 Mult(e_vec, DETERMINANTS, empty, empty, q_det);
693}
694
695/// @cond Suppress_Doxygen_warnings
696
697namespace
698{
699
700using namespace internal::quadrature_interpolator;
701
703using TensorEvalKernel = QuadratureInterpolator::TensorEvalKernelType;
705using CollocatedGradKernel = QuadratureInterpolator::CollocatedGradKernelType;
706
707template <QVectorLayout Q_LAYOUT>
708TensorEvalKernel FallbackTensorEvalKernel(int DIM)
709{
710 if (DIM == 1) { return Values1D<Q_LAYOUT>; }
711 else if (DIM == 2) { return Values2D<Q_LAYOUT>; }
712 else if (DIM == 3) { return Values3D<Q_LAYOUT>; }
713 else { MFEM_ABORT(""); }
714}
715
716template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS>
717GradKernel GetGradKernel(int DIM)
718{
719 if (DIM == 1) { return Derivatives1D<Q_LAYOUT, GRAD_PHYS>; }
720 else if (DIM == 2) { return Derivatives2D<Q_LAYOUT, GRAD_PHYS>; }
721 else if (DIM == 3) { return Derivatives3D<Q_LAYOUT, GRAD_PHYS>; }
722 else { MFEM_ABORT(""); }
723}
724
725
726template<QVectorLayout Q_LAYOUT>
727GradKernel GetGradKernel(int DIM, bool GRAD_PHYS)
728{
729 if (GRAD_PHYS) { return GetGradKernel<Q_LAYOUT, true>(DIM); }
730 else { return GetGradKernel<Q_LAYOUT, false>(DIM); }
731}
732
733template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS>
734CollocatedGradKernel GetCollocatedGradKernel(int DIM)
735{
736 if (DIM == 1) { return CollocatedDerivatives1D<Q_LAYOUT, GRAD_PHYS>; }
737 else if (DIM == 2) { return CollocatedDerivatives2D<Q_LAYOUT, GRAD_PHYS>; }
738 else if (DIM == 3) { return CollocatedDerivatives3D<Q_LAYOUT, GRAD_PHYS>; }
739 else { MFEM_ABORT(""); }
740}
741
742template<QVectorLayout Q_LAYOUT>
743CollocatedGradKernel GetCollocatedGradKernel(int DIM, bool GRAD_PHYS)
744{
745 if (GRAD_PHYS) { return GetCollocatedGradKernel<Q_LAYOUT, true>(DIM); }
746 else { return GetCollocatedGradKernel<Q_LAYOUT, false>(DIM); }
747}
748} // namespace
749
750template <int DIM, int VDIM, int ND, int NQ>
751EvalKernel QuadratureInterpolator::EvalKernels::Kernel()
752{
753 using namespace internal::quadrature_interpolator;
754 if (DIM == 1) { return Eval1D; }
755 else if (DIM == 2) { return Eval2D<VDIM,ND,NQ>; }
756 else if (DIM == 3) { return Eval3D<VDIM,ND,NQ>; }
757 else { MFEM_ABORT(""); }
758}
759
760template <int DIM>
761EvalKernel GetEvalKernelVDimFallback(int VDIM)
762{
763 using EvalKernels = QuadratureInterpolator::EvalKernels;
764 if (VDIM == 1) { return EvalKernels::Kernel<DIM,1,0,0>(); }
765 else if (VDIM == 2) { return EvalKernels::Kernel<DIM,2,0,0>(); }
766 else if (VDIM == 3) { return EvalKernels::Kernel<DIM,3,0,0>(); }
767 else { MFEM_ABORT(""); }
768}
769
770EvalKernel QuadratureInterpolator::EvalKernels::Fallback(
771 int DIM, int VDIM, int ND, int NQ)
772{
773 if (DIM == 1) { return GetEvalKernelVDimFallback<1>(VDIM); }
774 else if (DIM == 2) { return GetEvalKernelVDimFallback<2>(VDIM); }
775 else if (DIM == 3) { return GetEvalKernelVDimFallback<3>(VDIM); }
776 else { MFEM_ABORT(""); }
777}
778
779TensorEvalKernel QuadratureInterpolator::TensorEvalKernels::Fallback(
780 int DIM, QVectorLayout Q_LAYOUT, int, int, int)
781{
782 if (Q_LAYOUT == QVectorLayout::byNODES) { return FallbackTensorEvalKernel<QVectorLayout::byNODES>(DIM); }
783 else { return FallbackTensorEvalKernel<QVectorLayout::byVDIM>(DIM); }
784}
785
786GradKernel QuadratureInterpolator::GradKernels::Fallback(
787 int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int, int, int)
788{
789 if (Q_LAYOUT == QVectorLayout::byNODES) { return GetGradKernel<QVectorLayout::byNODES>(DIM, GRAD_PHYS); }
790 else { return GetGradKernel<QVectorLayout::byVDIM>(DIM, GRAD_PHYS); }
791}
792
793CollocatedGradKernel QuadratureInterpolator::CollocatedGradKernels::Fallback(
794 int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int, int)
795{
796 if (Q_LAYOUT == QVectorLayout::byNODES) { return GetCollocatedGradKernel<QVectorLayout::byNODES>(DIM, GRAD_PHYS); }
797 else { return GetCollocatedGradKernel<QVectorLayout::byVDIM>(DIM, GRAD_PHYS); }
798}
799
800/// @endcond
801
802namespace internal
803{
804namespace quadrature_interpolator
805{
806void InitEvalKernels()
807{
808 using k = QuadratureInterpolator::EvalKernels;
809 // 2D, VDIM = 1
810 k::Specialization<2,1,1,1>::Add();
811 k::Specialization<2,1,1,4>::Add();
812 // Q1
813 k::Specialization<2,1,4,4>::Add();
814 k::Specialization<2,1,4,9>::Add();
815 // Q2
816 k::Specialization<2,1,9,9>::Add();
817 k::Specialization<2,1,9,16>::Add();
818 // Q3
819 k::Specialization<2,1,16,16>::Add();
820 k::Specialization<2,1,16,25>::Add();
821 k::Specialization<2,1,16,36>::Add();
822 // Q4
823 k::Specialization<2,1,25,25>::Add();
824 k::Specialization<2,1,25,36>::Add();
825 k::Specialization<2,1,25,49>::Add();
826 k::Specialization<2,1,25,64>::Add();
827
828 // 3D, VDIM = 1
829 // Q0
830 k::Specialization<3,1,1,1>::Add();
831 k::Specialization<3,1,1,8>::Add();
832 // Q1
833 k::Specialization<3,1,8,8>::Add();
834 k::Specialization<3,1,8,27>::Add();
835 // Q2
836 k::Specialization<3,1,27,27>::Add();
837 k::Specialization<3,1,27,64>::Add();
838 // Q3
839 k::Specialization<3,1,64,64>::Add();
840 k::Specialization<3,1,64,125>::Add();
841 k::Specialization<3,1,64,216>::Add();
842 // Q4
843 k::Specialization<3,1,125,125>::Add();
844 k::Specialization<3,1,125,216>::Add();
845
846 // 2D, VDIM = 3
847 // Q0
848 k::Specialization<2,3,1,1>::Add();
849 k::Specialization<2,3,1,4>::Add();
850 // Q1
851 k::Specialization<2,3,4,4>::Add();
852 k::Specialization<2,3,4,9>::Add();
853 // Q2
854 k::Specialization<2,3,9,4>::Add();
855 k::Specialization<2,3,9,9>::Add();
856 k::Specialization<2,3,9,16>::Add();
857 k::Specialization<2,3,9,25>::Add();
858 // Q3
859 k::Specialization<2,3,16,16>::Add();
860 k::Specialization<2,3,16,25>::Add();
861 k::Specialization<2,3,16,36>::Add();
862 // Q4
863 k::Specialization<2,3,25,25>::Add();
864 k::Specialization<2,3,25,36>::Add();
865 k::Specialization<2,3,25,49>::Add();
866 k::Specialization<2,3,25,64>::Add();
867
868 // 2D, VDIM = 2
869 // Q1
870 k::Specialization<2,2,4,4>::Add();
871 k::Specialization<2,2,4,9>::Add();
872 // Q2
873 k::Specialization<2,2,9,9>::Add();
874 k::Specialization<2,2,9,16>::Add();
875 // Q3
876 k::Specialization<2,2,16,16>::Add();
877 k::Specialization<2,2,16,25>::Add();
878 k::Specialization<2,2,16,36>::Add();
879 // Q4
880 k::Specialization<2,2,25,25>::Add();
881 k::Specialization<2,2,25,36>::Add();
882 k::Specialization<2,2,25,49>::Add();
883 k::Specialization<2,2,25,64>::Add();
884
885 // 3D, VDIM = 3
886 // Q1
887 k::Specialization<3,3,8,8>::Add();
888 k::Specialization<3,3,8,27>::Add();
889 // Q2
890 k::Specialization<3,3,27,27>::Add();
891 k::Specialization<3,3,27,64>::Add();
892 k::Specialization<3,3,27,125>::Add();
893 // Q3
894 k::Specialization<3,3,64,64>::Add();
895 k::Specialization<3,3,64,125>::Add();
896 k::Specialization<3,3,64,216>::Add();
897 // Q4
898 k::Specialization<3,3,125,125>::Add();
899 k::Specialization<3,3,125,216>::Add();
900}
901
902} // namespace quadrature_Interpolator
903} // namespace internal
904
905} // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:381
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:141
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition fe_base.hpp:174
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:214
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition fe_base.hpp:154
@ FULL
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition fe_base.hpp:158
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition fe_base.hpp:145
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:208
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:678
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3824
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:869
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:644
int GetVDim() const
Returns the vector dimension of the finite element space.
Definition fespace.hpp:819
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3860
Abstract class for all finite elements.
Definition fe_base.hpp:247
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:375
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:324
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:363
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition mesh.hpp:3076
const Mesh * mesh
Definition mesh.hpp:3082
Vector J
Jacobians of the element transformations at all quadrature points.
Definition mesh.hpp:3115
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
Mesh data type.
Definition mesh.hpp:65
bool IsMixedMesh() const
Returns true if the mesh is a mixed mesh, false otherwise.
Definition mesh.cpp:7582
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1306
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1309
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:883
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:7558
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
void(*)(const int, const real_t *, const real_t *, const real_t *, real_t *, const int, const int, const int) CollocatedGradKernelType
bool use_tensor_products
Tensor product evaluation mode.
@ VALUES
Evaluate the values at quadrature points.
@ DERIVATIVES
Evaluate the derivatives at quadrature points.
@ PHYSICAL_DERIVATIVES
Evaluate the physical derivatives.
@ DETERMINANTS
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
void(*)(const int, const int, const QVectorLayout, const GeometricFactors *, const DofToQuad &, const Vector &, Vector &, Vector &, Vector &, const int) EvalKernelType
QuadratureInterpolator(const FiniteElementSpace &fes, const IntegrationRule &ir)
void(*)(const int, const real_t *, const real_t *, const real_t *, const real_t *, real_t *, const int, const int, const int, const int) GradKernelType
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.
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...
void MultTranspose(unsigned eval_flags, const Vector &q_val, const Vector &q_der, Vector &e_vec) const
Perform the transpose operation of Mult(). (TODO)
void(*)(const int, const real_t *, const real_t *, real_t *, const int, const int, const int) TensorEvalKernelType
void Values(const Vector &e_vec, Vector &q_val) const
Interpolate the values of the E-vector e_vec at quadrature points.
void PhysValues(const Vector &e_vec, Vector &q_val) const
Interpolate the physical values of the E-vector e_vec at quadrature points.
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...
QVectorLayout q_layout
Output Q-vector layout.
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.
static bool SupportsFESpace(const FiniteElementSpace &fespace)
Returns true if the given finite element space is supported by QuadratureInterpolator.
const IntegrationRule * IntRule
Not owned.
Vector d_buffer
Auxiliary device buffer.
const FiniteElementSpace * fespace
Not owned.
const QuadratureSpace * qspace
Not owned.
void MultHDiv(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_div) const
Auxiliary method called by Mult() when using H(div)-conforming space.
Class representing the storage layout of a QuadratureFunction.
Definition qspace.hpp:164
const IntegrationRule & GetElementIntRule(int idx) const
Get the IntegrationRule associated with mesh element idx.
Definition qspace.hpp:193
Vector data type.
Definition vector.hpp:82
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:520
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:145
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:528
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
constexpr int DIM
mfem::real_t real_t
MFEM_HOST_DEVICE T det(const tensor< T, 1, 1 > &A)
Returns the determinant of a matrix.
Definition tensor.hpp:1406
MFEM_HOST_DEVICE void CalcInverse(const T *data, T *inv_data)
Return the inverse of a matrix with given size and data into the matrix with data inv_data.
Definition kernels.hpp:306
MFEM_HOST_DEVICE T Det(const T *data)
Compute the determinant of a square matrix of size dim with given data.
Definition kernels.hpp:297
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition device.hpp:365
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:138
SchrodingerBaseKernels< ParMesh, ParFiniteElementSpace, ParComplexGridFunction, ParGridFunction, ParBilinearForm, ParMixedBilinearForm, ParLinearForm > Kernels
void forall_2D(int N, int X, int Y, lambda &&body)
Definition forall.hpp:925
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition fespace.hpp:1548
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition fespace.hpp:31
float real_t
Definition config.hpp:46
void forall(int N, lambda &&body)
Definition forall.hpp:839