MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
quadinterpolator_face.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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_face.hpp"
13#include "../general/annotation.hpp"
14#include "../general/forall.hpp"
15#include "../linalg/dtensor.hpp"
16#include "../linalg/kernels.hpp"
17
18namespace mfem
19{
20
21/// Return the sign to apply to the normals on each face to point from e1 to e2.
22static void GetSigns(const FiniteElementSpace &fes, const FaceType type,
23 Array<bool> &signs)
24{
25 const Mesh &mesh = *fes.GetMesh();
26 const int dim = mesh.SpaceDimension();
27 int face_id;
28 int f_ind = 0;
29 for (int f = 0; f < mesh.GetNumFacesWithGhost(); ++f)
30 {
31 Mesh::FaceInformation face = mesh.GetFaceInformation(f);
32 face_id = face.element[0].local_face_id;
33 if (face.IsNonconformingCoarse())
34 {
35 // We skip nonconforming coarse-fine faces as they are treated
36 // by the corresponding nonconforming fine-coarse faces.
37 continue;
38 }
39 else if ( face.IsOfFaceType(type) )
40 {
41 if (dim==2)
42 {
43 if (face_id==2 || face_id==3)
44 {
45 signs[f_ind] = true;
46 }
47 else
48 {
49 signs[f_ind] = false;
50 }
51 }
52 else if (dim==3)
53 {
54 if (face_id==0 || face_id==3 || face_id==4)
55 {
56 signs[f_ind] = true;
57 }
58 else
59 {
60 signs[f_ind] = false;
61 }
62 }
63 f_ind++;
64 }
65 }
66}
67
68FaceQuadratureInterpolator::FaceQuadratureInterpolator(
69 const FiniteElementSpace &fes,
70 const IntegrationRule &ir, FaceType type_)
71 : type(type_), nf(fes.GetNFbyType(type)), signs(nf),
72 q_layout(QVectorLayout::byNODES)
73{
74 fespace = &fes;
75 IntRule = &ir;
76 use_tensor_products = true; // not implemented yet (not used)
77
78 if (fespace->GetNE() == 0) { return; }
79 GetSigns(*fespace, type, signs);
80 const FiniteElement *fe = fespace->GetFE(0);
81 const ScalarFiniteElement *sfe =
82 dynamic_cast<const ScalarFiniteElement*>(fe);
83 const TensorBasisElement *tfe =
84 dynamic_cast<const TensorBasisElement*>(fe);
85 MFEM_VERIFY(sfe != NULL, "Only scalar finite elements are supported");
86 MFEM_VERIFY(tfe != NULL &&
87 (tfe->GetBasisType()==BasisType::GaussLobatto ||
88 tfe->GetBasisType()==BasisType::Positive),
89 "Only Gauss-Lobatto and Bernstein basis are supported in "
90 "FaceQuadratureInterpolator.");
91}
92
93template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
94void FaceQuadratureInterpolator::Eval2D(
95 const int NF,
96 const int vdim,
97 const QVectorLayout q_layout,
98 const DofToQuad &maps,
99 const Array<bool> &signs,
100 const Vector &f_vec,
101 Vector &q_val,
102 Vector &q_der,
103 Vector &q_det,
104 Vector &q_nor,
105 const int eval_flags)
106{
107 const int nd1d = maps.ndof;
108 const int nq1d = maps.nqpt;
109 const int ND1D = T_ND1D ? T_ND1D : nd1d;
110 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
111 const int VDIM = T_VDIM ? T_VDIM : vdim;
112 MFEM_VERIFY(ND1D <= MAX_ND1D, "");
113 MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
114 MFEM_VERIFY(VDIM == 2 || !(eval_flags & DETERMINANTS), "");
115 auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
116 auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
117 auto F = Reshape(f_vec.Read(), ND1D, VDIM, NF);
118 auto sign = signs.Read();
119 auto val = q_layout == QVectorLayout::byNODES ?
120 Reshape(q_val.Write(), NQ1D, VDIM, NF):
121 Reshape(q_val.Write(), VDIM, NQ1D, NF);
122 // auto der = Reshape(q_der.Write(), NQ1D, VDIM, NF); // only tangential der
123 auto det = Reshape(q_det.Write(), NQ1D, NF);
124 auto n = q_layout == QVectorLayout::byNODES ?
125 Reshape(q_nor.Write(), NQ1D, 2, NF):
126 Reshape(q_nor.Write(), 2, NQ1D, NF);
127 MFEM_VERIFY(eval_flags | DERIVATIVES,
128 "Derivatives on the faces are not yet supported.");
129 // If Gauss-Lobatto
130 mfem::forall(NF, [=] MFEM_HOST_DEVICE (int f)
131 {
132 const int ND1D = T_ND1D ? T_ND1D : nd1d;
133 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
134 const int VDIM = T_VDIM ? T_VDIM : vdim;
135 constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
136 constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM2D;
137 real_t r_F[max_ND1D][max_VDIM];
138 for (int d = 0; d < ND1D; d++)
139 {
140 for (int c = 0; c < VDIM; c++)
141 {
142 r_F[d][c] = F(d,c,f);
143 }
144 }
145 for (int q = 0; q < NQ1D; ++q)
146 {
147 if (eval_flags & VALUES)
148 {
149 real_t ed[max_VDIM];
150 for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
151 for (int d = 0; d < ND1D; ++d)
152 {
153 const real_t b = B(q,d);
154 for (int c = 0; c < VDIM; c++) { ed[c] += b*r_F[d][c]; }
155 }
156 for (int c = 0; c < VDIM; c++)
157 {
158 if (q_layout == QVectorLayout::byVDIM) { val(c,q,f) = ed[c]; }
159 if (q_layout == QVectorLayout::byNODES) { val(q,c,f) = ed[c]; }
160 }
161 }
162 if ((eval_flags & DERIVATIVES)
163 || (eval_flags & DETERMINANTS)
164 || (eval_flags & NORMALS))
165 {
166 real_t D[max_VDIM];
167 for (int i = 0; i < VDIM; i++) { D[i] = 0.0; }
168 for (int d = 0; d < ND1D; ++d)
169 {
170 const real_t w = G(q,d);
171 for (int c = 0; c < VDIM; c++)
172 {
173 real_t s_e = r_F[d][c];
174 D[c] += s_e * w;
175 }
176 }
177 if (VDIM == 2 &&
178 ((eval_flags & NORMALS)
179 || (eval_flags & DETERMINANTS)))
180 {
181 const real_t norm = sqrt(D[0]*D[0]+D[1]*D[1]);
182 if (eval_flags & DETERMINANTS)
183 {
184 det(q,f) = norm;
185 }
186 if (eval_flags & NORMALS)
187 {
188 const real_t s = sign[f] ? -1.0 : 1.0;
189 if (q_layout == QVectorLayout::byVDIM)
190 {
191 n(0,q,f) = s*D[1]/norm;
192 n(1,q,f) = -s*D[0]/norm;
193 }
194 if (q_layout == QVectorLayout::byNODES)
195 {
196 n(q,0,f) = s*D[1]/norm;
197 n(q,1,f) = -s*D[0]/norm;
198 }
199 }
200 }
201 }
202 }
203 });
204}
205
206template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
207void FaceQuadratureInterpolator::Eval3D(
208 const int NF,
209 const int vdim,
210 const QVectorLayout q_layout,
211 const DofToQuad &maps,
212 const Array<bool> &signs,
213 const Vector &e_vec,
214 Vector &q_val,
215 Vector &q_der,
216 Vector &q_det,
217 Vector &q_nor,
218 const int eval_flags)
219{
220 const int nd1d = maps.ndof;
221 const int nq1d = maps.nqpt;
222 const int ND1D = T_ND1D ? T_ND1D : nd1d;
223 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
224 const int VDIM = T_VDIM ? T_VDIM : vdim;
225 MFEM_VERIFY(ND1D <= MAX_ND1D, "");
226 MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
227 MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
228 auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
229 auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
230 auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
231 auto sign = signs.Read();
232 auto val = q_layout == QVectorLayout::byNODES ?
233 Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF):
234 Reshape(q_val.Write(), VDIM, NQ1D, NQ1D, NF);
235 // auto der = Reshape(q_der.Write(), NQ1D, VDIM, 3, NF);
236 auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
237 auto nor = q_layout == QVectorLayout::byNODES ?
238 Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF):
239 Reshape(q_nor.Write(), 3, NQ1D, NQ1D, NF);
240 MFEM_VERIFY(eval_flags | DERIVATIVES,
241 "Derivatives on the faces are not yet supported.");
242 mfem::forall(NF, [=] MFEM_HOST_DEVICE (int f)
243 {
244 constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
245 constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
246 constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
247 real_t r_F[max_ND1D][max_ND1D][max_VDIM];
248 for (int d1 = 0; d1 < ND1D; d1++)
249 {
250 for (int d2 = 0; d2 < ND1D; d2++)
251 {
252 for (int c = 0; c < VDIM; c++)
253 {
254 r_F[d1][d2][c] = F(d1,d2,c,f);
255 }
256 }
257 }
258 if (eval_flags & VALUES)
259 {
260 real_t Bu[max_NQ1D][max_ND1D][max_VDIM];
261 for (int d2 = 0; d2 < ND1D; ++d2)
262 {
263 for (int q = 0; q < NQ1D; ++q)
264 {
265 for (int c = 0; c < VDIM; c++) { Bu[q][d2][c] = 0.0; }
266 for (int d1 = 0; d1 < ND1D; ++d1)
267 {
268 const real_t b = B(q,d1);
269 for (int c = 0; c < VDIM; c++)
270 {
271 Bu[q][d2][c] += b*r_F[d1][d2][c];
272 }
273 }
274 }
275 }
276 real_t BBu[max_NQ1D][max_NQ1D][max_VDIM];
277 for (int q2 = 0; q2 < NQ1D; ++q2)
278 {
279 for (int q1 = 0; q1 < NQ1D; ++q1)
280 {
281 for (int c = 0; c < VDIM; c++) { BBu[q2][q1][c] = 0.0; }
282 for (int d2 = 0; d2 < ND1D; ++d2)
283 {
284 const real_t b = B(q2,d2);
285 for (int c = 0; c < VDIM; c++)
286 {
287 BBu[q2][q1][c] += b*Bu[q1][d2][c];
288 }
289 }
290 for (int c = 0; c < VDIM; c++)
291 {
292 const real_t v = BBu[q2][q1][c];
293 if (q_layout == QVectorLayout::byVDIM) { val(c,q1,q2,f) = v; }
294 if (q_layout == QVectorLayout::byNODES) { val(q1,q2,c,f) = v; }
295 }
296 }
297 }
298 }
299 if ((eval_flags & DERIVATIVES)
300 || (eval_flags & DETERMINANTS)
301 || (eval_flags & NORMALS))
302 {
303 // We only compute the tangential derivatives
304 real_t Gu[max_NQ1D][max_ND1D][max_VDIM];
305 real_t Bu[max_NQ1D][max_ND1D][max_VDIM];
306 for (int d2 = 0; d2 < ND1D; ++d2)
307 {
308 for (int q = 0; q < NQ1D; ++q)
309 {
310 for (int c = 0; c < VDIM; c++)
311 {
312 Gu[q][d2][c] = 0.0;
313 Bu[q][d2][c] = 0.0;
314 }
315 for (int d1 = 0; d1 < ND1D; ++d1)
316 {
317 const real_t b = B(q,d1);
318 const real_t g = G(q,d1);
319 for (int c = 0; c < VDIM; c++)
320 {
321 const real_t u = r_F[d1][d2][c];
322 Gu[q][d2][c] += g*u;
323 Bu[q][d2][c] += b*u;
324 }
325 }
326 }
327 }
328 real_t BGu[max_NQ1D][max_NQ1D][max_VDIM];
329 real_t GBu[max_NQ1D][max_NQ1D][max_VDIM];
330 for (int q2 = 0; q2 < NQ1D; ++q2)
331 {
332 for (int q1 = 0; q1 < NQ1D; ++q1)
333 {
334 for (int c = 0; c < VDIM; c++)
335 {
336 BGu[q2][q1][c] = 0.0;
337 GBu[q2][q1][c] = 0.0;
338 }
339 for (int d2 = 0; d2 < ND1D; ++d2)
340 {
341 const real_t b = B(q2,d2);
342 const real_t g = G(q2,d2);
343 for (int c = 0; c < VDIM; c++)
344 {
345 BGu[q2][q1][c] += b*Gu[q1][d2][c];
346 GBu[q2][q1][c] += g*Bu[q1][d2][c];
347 }
348 }
349 }
350 }
351 if (VDIM == 3 && ((eval_flags & NORMALS) ||
352 (eval_flags & DETERMINANTS)))
353 {
354 real_t n[3];
355 for (int q2 = 0; q2 < NQ1D; ++q2)
356 {
357 for (int q1 = 0; q1 < NQ1D; ++q1)
358 {
359 const real_t s = sign[f] ? -1.0 : 1.0;
360 n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
361 BGu[q2][q1][2] );
362 n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
363 BGu[q2][q1][2] );
364 n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
365 BGu[q2][q1][1] );
366 const real_t norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
367 if (eval_flags & DETERMINANTS) { det(q1,q2,f) = norm; }
368 if (eval_flags & NORMALS)
369 {
370 if (q_layout == QVectorLayout::byVDIM)
371 {
372 nor(0,q1,q2,f) = n[0]/norm;
373 nor(1,q1,q2,f) = n[1]/norm;
374 nor(2,q1,q2,f) = n[2]/norm;
375 }
376 if (q_layout == QVectorLayout::byNODES)
377 {
378 nor(q1,q2,0,f) = n[0]/norm;
379 nor(q1,q2,1,f) = n[1]/norm;
380 nor(q1,q2,2,f) = n[2]/norm;
381 }
382 }
383 }
384 }
385 }
386 }
387 });
388}
389
390template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
391void FaceQuadratureInterpolator::SmemEval3D(
392 const int NF,
393 const int vdim,
394 const QVectorLayout q_layout,
395 const DofToQuad &maps,
396 const Array<bool> &signs,
397 const Vector &e_vec,
398 Vector &q_val,
399 Vector &q_der,
400 Vector &q_det,
401 Vector &q_nor,
402 const int eval_flags)
403{
404 MFEM_PERF_SCOPE("FaceQuadInterpolator::SmemEval3D");
405 const int nd1d = maps.ndof;
406 const int nq1d = maps.nqpt;
407 const int ND1D = T_ND1D ? T_ND1D : nd1d;
408 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
409 const int VDIM = T_VDIM ? T_VDIM : vdim;
410 MFEM_VERIFY(ND1D <= MAX_ND1D, "");
411 MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
412 MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
413 auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
414 auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
415 auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
416 auto sign = signs.Read();
417 auto val = q_layout == QVectorLayout::byNODES ?
418 Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF):
419 Reshape(q_val.Write(), VDIM, NQ1D, NQ1D, NF);
420 // auto der = Reshape(q_der.Write(), NQ1D, VDIM, 3, NF);
421 auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
422 auto nor = q_layout == QVectorLayout::byNODES ?
423 Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF):
424 Reshape(q_nor.Write(), 3, NQ1D, NQ1D, NF);
425 MFEM_VERIFY(eval_flags | DERIVATIVES,
426 "Derivatives on the faces are not yet supported.");
427
428 mfem::forall_3D(NF, NQ1D, NQ1D, VDIM, [=] MFEM_HOST_DEVICE (int f)
429 {
430 constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
431 constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
432 constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
433
434 MFEM_SHARED real_t sm1[max_NQ1D*max_NQ1D*max_VDIM];
435 MFEM_SHARED real_t sm2[max_NQ1D*max_ND1D*max_VDIM];
436
437 auto s_F = (real_t(*)[max_ND1D][max_VDIM])sm1;
438 MFEM_FOREACH_THREAD(d1,x,ND1D)
439 {
440 MFEM_FOREACH_THREAD(d2,y,ND1D)
441 {
442 MFEM_FOREACH_THREAD(c,z,VDIM)
443 {
444 s_F[d1][d2][c] = F(d1,d2,c,f);
445 }
446 }
447 }
448 MFEM_SYNC_THREAD;
449
450 if (eval_flags & VALUES)
451 {
452 auto Bu = (real_t (*)[max_ND1D][max_VDIM])sm2;
453 MFEM_FOREACH_THREAD(d2,x,ND1D)
454 {
455 MFEM_FOREACH_THREAD(q1,y,NQ1D)
456 {
457 MFEM_FOREACH_THREAD(c,z,VDIM)
458 {
459 real_t thrdBu = 0.0;
460 for (int d1 = 0; d1 < ND1D; ++d1)
461 {
462 thrdBu += B(q1,d1)*s_F[d1][d2][c];
463 }
464 Bu[q1][d2][c] = thrdBu;
465 }
466 }
467 }
468 MFEM_SYNC_THREAD;
469
470 MFEM_FOREACH_THREAD(q2,x,NQ1D)
471 {
472 MFEM_FOREACH_THREAD(q1,y,NQ1D)
473 {
474 MFEM_FOREACH_THREAD(c,z,VDIM)
475 {
476 real_t v = 0.0;
477 for (int d2 = 0; d2 < ND1D; ++d2)
478 {
479 v += B(q2,d2)*Bu[q1][d2][c];
480 }
481 if (q_layout == QVectorLayout::byVDIM) { val(c,q1,q2,f) = v; }
482 if (q_layout == QVectorLayout::byNODES) { val(q1,q2,c,f) = v; }
483 }
484 }
485 }
486 }
487
488 if ((eval_flags & DERIVATIVES)
489 || (eval_flags & DETERMINANTS)
490 || (eval_flags & NORMALS))
491 {
492 // We only compute the tangential derivatives
493 auto Gu = (real_t (*)[max_ND1D][max_VDIM])sm2;
494 MFEM_SHARED real_t Bu[max_NQ1D][max_ND1D][max_VDIM];
495 MFEM_FOREACH_THREAD(d2,x,ND1D)
496 {
497 MFEM_FOREACH_THREAD(q1,y,NQ1D)
498 {
499 MFEM_FOREACH_THREAD(c,z,VDIM)
500 {
501 real_t thrdGu = 0;
502 real_t thrdBu = 0;
503 for (int d1 = 0; d1 < ND1D; ++d1)
504 {
505 const real_t u = s_F[d1][d2][c];
506 thrdBu += B(q1,d1)*u;
507 thrdGu += G(q1,d1)*u;
508 }
509 Gu[q1][d2][c] = thrdGu;
510 Bu[q1][d2][c] = thrdBu;
511 }
512 }
513 }
514 MFEM_SYNC_THREAD;
515
516 auto BGu = (real_t (*)[max_NQ1D][max_VDIM])sm1;
517 MFEM_SHARED real_t GBu[max_NQ1D][max_NQ1D][max_VDIM];
518 MFEM_FOREACH_THREAD(q2,x,NQ1D)
519 {
520 MFEM_FOREACH_THREAD(q1,y,NQ1D)
521 {
522 MFEM_FOREACH_THREAD(c,z,VDIM)
523 {
524 real_t thrdBGu = 0.0;
525 real_t thrdGBu = 0.0;
526 for (int d2 = 0; d2 < ND1D; ++d2)
527 {
528 thrdBGu += B(q2,d2)*Gu[q1][d2][c];
529 thrdGBu += G(q2,d2)*Bu[q1][d2][c];
530 }
531 BGu[q2][q1][c] = thrdBGu;
532 GBu[q2][q1][c] = thrdGBu;
533 }
534 }
535 }
536 MFEM_SYNC_THREAD;
537
538 if (VDIM == 3 && ((eval_flags & NORMALS) ||
539 (eval_flags & DETERMINANTS)))
540 {
541 real_t n[3];
542 MFEM_FOREACH_THREAD(q2,x,NQ1D)
543 {
544 MFEM_FOREACH_THREAD(q1,y,NQ1D)
545 {
546 if (MFEM_THREAD_ID(z) == 0)
547 {
548 const real_t s = sign[f] ? -1.0 : 1.0;
549 n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
550 BGu[q2][q1][2] );
551 n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
552 BGu[q2][q1][2] );
553 n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
554 BGu[q2][q1][1] );
555
556 const real_t norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
557
558 if (eval_flags & DETERMINANTS) { det(q1,q2,f) = norm; }
559
560 if (eval_flags & NORMALS)
561 {
562 if (q_layout == QVectorLayout::byVDIM)
563 {
564 nor(0,q1,q2,f) = n[0]/norm;
565 nor(1,q1,q2,f) = n[1]/norm;
566 nor(2,q1,q2,f) = n[2]/norm;
567 }
568 if (q_layout == QVectorLayout::byNODES)
569 {
570 nor(q1,q2,0,f) = n[0]/norm;
571 nor(q1,q2,1,f) = n[1]/norm;
572 nor(q1,q2,2,f) = n[2]/norm;
573 }
574 }
575 }
576 }
577 }
578 }
579 }
580 });
581}
582
583void FaceQuadratureInterpolator::Mult(
584 const Vector &e_vec, unsigned eval_flags,
585 Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
586{
587 if (nf == 0) { return; }
588 const int vdim = fespace->GetVDim();
589 const int dim = fespace->GetMesh()->Dimension();
590 const FiniteElement *fe =
591 fespace->GetTraceElement(0, fespace->GetMesh()->GetFaceGeometry(0));
592 const IntegrationRule *ir = IntRule;
593 const DofToQuad &maps = fe->GetDofToQuad(*ir, DofToQuad::TENSOR);
594 const int nd1d = maps.ndof;
595 const int nq1d = maps.nqpt;
596 void (*eval_func)(
597 const int NF,
598 const int vdim,
599 const QVectorLayout q_layout,
600 const DofToQuad &maps,
601 const Array<bool> &signs,
602 const Vector &e_vec,
603 Vector &q_val,
604 Vector &q_der,
605 Vector &q_det,
606 Vector &q_nor,
607 const int eval_flags) = NULL;
608 if (vdim == 1)
609 {
610 if (dim == 2)
611 {
612 switch (10*nd1d + nq1d)
613 {
614 // Q0
615 case 11: eval_func = &Eval2D<1,1,1>; break;
616 case 12: eval_func = &Eval2D<1,1,2>; break;
617 // Q1
618 case 22: eval_func = &Eval2D<1,2,2>; break;
619 case 23: eval_func = &Eval2D<1,2,3>; break;
620 // Q2
621 case 33: eval_func = &Eval2D<1,3,3>; break;
622 case 34: eval_func = &Eval2D<1,3,4>; break;
623 // Q3
624 case 44: eval_func = &Eval2D<1,4,4>; break;
625 case 45: eval_func = &Eval2D<1,4,5>; break;
626 case 46: eval_func = &Eval2D<1,4,6>; break;
627 // Q4
628 case 55: eval_func = &Eval2D<1,5,5>; break;
629 case 56: eval_func = &Eval2D<1,5,6>; break;
630 case 57: eval_func = &Eval2D<1,5,7>; break;
631 case 58: eval_func = &Eval2D<1,5,8>; break;
632 }
633 if (nq1d >= 10 || !eval_func)
634 {
635 eval_func = &Eval2D<1>;
636 }
637 }
638 else if (dim == 3)
639 {
640 switch (10*nd1d + nq1d)
641 {
642 // Q0
643 case 11: eval_func = &SmemEval3D<1,1,1>; break;
644 case 12: eval_func = &SmemEval3D<1,1,2>; break;
645 // Q1
646 case 22: eval_func = &SmemEval3D<1,2,2>; break;
647 case 23: eval_func = &SmemEval3D<1,2,3>; break;
648 case 24: eval_func = &SmemEval3D<1,2,4>; break;
649 // Q2
650 case 33: eval_func = &SmemEval3D<1,3,3>; break;
651 case 34: eval_func = &SmemEval3D<1,3,4>; break;
652 // Q3
653 case 44: eval_func = &SmemEval3D<1,4,4>; break;
654 case 45: eval_func = &SmemEval3D<1,4,5>; break;
655 case 46: eval_func = &SmemEval3D<1,4,6>; break;
656 // Q4
657 case 55: eval_func = &SmemEval3D<1,5,5>; break;
658 case 56: eval_func = &SmemEval3D<1,5,6>; break;
659 }
660 if (nq1d >= 10 || !eval_func)
661 {
662 eval_func = &Eval3D<1>;
663 }
664 }
665 }
666 else if (vdim == dim)
667 {
668 if (dim == 2)
669 {
670 switch (10*nd1d + nq1d)
671 {
672 // Q1
673 case 22: eval_func = &Eval2D<2,2,2>; break;
674 case 23: eval_func = &Eval2D<2,2,3>; break;
675 // Q2
676 case 33: eval_func = &Eval2D<2,3,3>; break;
677 case 34: eval_func = &Eval2D<2,3,4>; break;
678 // Q3
679 case 44: eval_func = &Eval2D<2,4,4>; break;
680 case 45: eval_func = &Eval2D<2,4,5>; break;
681 case 46: eval_func = &Eval2D<2,4,6>; break;
682 // Q4
683 case 55: eval_func = &Eval2D<2,5,5>; break;
684 case 56: eval_func = &Eval2D<2,5,6>; break;
685 case 57: eval_func = &Eval2D<2,5,7>; break;
686 case 58: eval_func = &Eval2D<2,5,8>; break;
687 }
688 if (nq1d >= 10 || !eval_func)
689 {
690 eval_func = &Eval2D<2>;
691 }
692 }
693 else if (dim == 3)
694 {
695 switch (10*nd1d + nq1d)
696 {
697 // Q1
698 case 22: eval_func = &SmemEval3D<3,2,2>; break;
699 case 23: eval_func = &SmemEval3D<3,2,3>; break;
700 case 24: eval_func = &SmemEval3D<3,2,4>; break;
701 // Q2
702 case 33: eval_func = &SmemEval3D<3,3,3>; break;
703 case 34: eval_func = &SmemEval3D<3,3,4>; break;
704 // Q3
705 case 44: eval_func = &SmemEval3D<3,4,4>; break;
706 case 45: eval_func = &SmemEval3D<3,4,5>; break;
707 case 46: eval_func = &SmemEval3D<3,4,6>; break;
708 // Q4
709 case 55: eval_func = &SmemEval3D<3,5,5>; break;
710 case 56: eval_func = &SmemEval3D<3,5,6>; break;
711 }
712 if (nq1d >= 10 || !eval_func)
713 {
714 eval_func = &Eval3D<3>;
715 }
716 }
717 }
718 if (eval_func)
719 {
720 eval_func(nf, vdim, q_layout, maps, signs, e_vec,
721 q_val, q_der, q_det, q_nor, eval_flags);
722 }
723 else
724 {
725 MFEM_ABORT("case not supported yet");
726 }
727}
728
729void FaceQuadratureInterpolator::Values(
730 const Vector &e_vec, Vector &q_val) const
731{
732 Vector q_der, q_det, q_nor;
733 Mult(e_vec, VALUES, q_val, q_der, q_det, q_nor);
734}
735
736} // namespace mfem
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:317
mfem::BasisType::GaussLobatto
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
mfem::BasisType::Positive
@ Positive
Bernstein polynomials.
Definition fe_base.hpp:33
mfem::DofToQuad
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:137
mfem::DofToQuad::G
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:210
mfem::DofToQuad::TENSOR
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:161
mfem::DofToQuad::B
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:189
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:174
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
mfem::FaceQuadratureInterpolator::DERIVATIVES
@ DERIVATIVES
Evaluate the derivatives at quadrature points.
Definition quadinterpolator_face.hpp:56
mfem::FaceQuadratureInterpolator::DETERMINANTS
@ DETERMINANTS
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
Definition quadinterpolator_face.hpp:60
mfem::FaceQuadratureInterpolator::VALUES
@ VALUES
Evaluate the values at quadrature points.
Definition quadinterpolator_face.hpp:55
mfem::FaceQuadratureInterpolator::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_face.cpp:729
mfem::FaceQuadratureInterpolator::fespace
const FiniteElementSpace * fespace
Not owned.
Definition quadinterpolator_face.hpp:34
mfem::FaceQuadratureInterpolator::Eval3D
static void Eval3D(const int NF, const int vdim, const QVectorLayout q_layout, const DofToQuad &maps, const Array< bool > &signs, const Vector &e_vec, Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor, const int eval_flags)
Template compute kernel for 3D.
Definition quadinterpolator_face.cpp:207
mfem::FaceQuadratureInterpolator::FaceQuadratureInterpolator
FaceQuadratureInterpolator(const FiniteElementSpace &fes, const IntegrationRule &ir, FaceType type)
Definition quadinterpolator_face.cpp:68
mfem::FaceQuadratureInterpolator::q_layout
QVectorLayout q_layout
Output Q-vector layout.
Definition quadinterpolator_face.hpp:36
mfem::FaceQuadratureInterpolator::Mult
void Mult(const Vector &e_vec, unsigned eval_flags, Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
Interpolate the E-vector e_vec to quadrature points.
Definition quadinterpolator_face.cpp:583
mfem::FaceQuadratureInterpolator::IntRule
const IntegrationRule * IntRule
Not owned.
Definition quadinterpolator_face.hpp:35
mfem::FaceQuadratureInterpolator::SmemEval3D
static void SmemEval3D(const int NF, const int vdim, const QVectorLayout q_layout, const DofToQuad &maps, const Array< bool > &signs, const Vector &e_vec, Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor, const int eval_flags)
Definition quadinterpolator_face.cpp:391
mfem::FaceQuadratureInterpolator::Eval2D
static void Eval2D(const int NF, const int vdim, const QVectorLayout q_layout, const DofToQuad &maps, const Array< bool > &signs, const Vector &e_vec, Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor, const int eval_flags)
Template compute kernel for 2D.
Definition quadinterpolator_face.cpp:94
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
mfem::FiniteElementSpace::GetFE
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:3168
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:740
mfem::FiniteElementSpace::GetTraceElement
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
Definition fespace.cpp:3276
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:239
mfem::FiniteElement::GetDofToQuad
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition fe_base.cpp:365
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::Mesh::GetFaceGeometry
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
Definition mesh.cpp:1452
mfem::Mesh::Dimension
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
mfem::ScalarFiniteElement
Class for finite elements with basis functions that return scalar values.
Definition fe_base.hpp:656
mfem::Vector
Vector data type.
Definition vector.hpp:80
mfem::Vector::Read
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:474
mfem::Vector::Write
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:482
dim
int dim
Definition ex24.cpp:53
b
real_t b
Definition lissajous.cpp:42
mfem::u
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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::f_vec
std::function< void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
Definition lor_mms.hpp:68
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition forall.hpp:775
mfem::QVectorLayout
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition fespace.hpp:53
mfem::QVectorLayout::byNODES
@ byNODES
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
mfem::QVectorLayout::byVDIM
@ byVDIM
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
mfem::real_t
float real_t
Definition config.hpp:43
mfem::f
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
mfem::forall
void forall(int N, lambda &&body)
Definition forall.hpp:754
mfem::FaceType
FaceType
Definition mesh.hpp:47
s
RefCoord s[3]
Definition ncmesh_tables.hpp:182
quadinterpolator_face.hpp
doxygen 1.11.0