MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
quadinterpolator_face.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_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->GetTypicalFE();
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 = q_layout == QVectorLayout::byNODES ? // only tangential der
123 Reshape(q_der.Write(), NQ1D, VDIM, NF):
124 Reshape(q_der.Write(), VDIM, NQ1D, NF);
125 auto det = Reshape(q_det.Write(), NQ1D, NF);
126 auto n = q_layout == QVectorLayout::byNODES ?
127 Reshape(q_nor.Write(), NQ1D, 2, NF):
128 Reshape(q_nor.Write(), 2, NQ1D, NF);
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 (eval_flags & DERIVATIVES)
178 {
179 for (int c = 0; c < VDIM; c++)
180 {
181 if (q_layout == QVectorLayout::byVDIM)
182 {
183 der(c,q,f) = D[c];
184 }
185 else // q_layout == QVectorLayout::byNODES
186 {
187 der(q,c,f) = D[c];
188 }
189 }
190 }
191 if (VDIM == 2 &&
192 ((eval_flags & NORMALS)
193 || (eval_flags & DETERMINANTS)))
194 {
195 const real_t norm = sqrt(D[0]*D[0]+D[1]*D[1]);
196 if (eval_flags & DETERMINANTS)
197 {
198 det(q,f) = norm;
199 }
200 if (eval_flags & NORMALS)
201 {
202 const real_t s = sign[f] ? -1.0 : 1.0;
203 if (q_layout == QVectorLayout::byVDIM)
204 {
205 n(0,q,f) = s*D[1]/norm;
206 n(1,q,f) = -s*D[0]/norm;
207 }
208 if (q_layout == QVectorLayout::byNODES)
209 {
210 n(q,0,f) = s*D[1]/norm;
211 n(q,1,f) = -s*D[0]/norm;
212 }
213 }
214 }
215 }
216 }
217 });
218}
219
220template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
221void FaceQuadratureInterpolator::Eval3D(
222 const int NF,
223 const int vdim,
224 const QVectorLayout q_layout,
225 const DofToQuad &maps,
226 const Array<bool> &signs,
227 const Vector &e_vec,
228 Vector &q_val,
229 Vector &q_der,
230 Vector &q_det,
231 Vector &q_nor,
232 const int eval_flags)
233{
234 const int nd1d = maps.ndof;
235 const int nq1d = maps.nqpt;
236 const int ND1D = T_ND1D ? T_ND1D : nd1d;
237 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
238 const int VDIM = T_VDIM ? T_VDIM : vdim;
239 MFEM_VERIFY(ND1D <= MAX_ND1D, "");
240 MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
241 MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
242 auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
243 auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
244 auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
245 auto sign = signs.Read();
246 auto val = q_layout == QVectorLayout::byNODES ?
247 Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF):
248 Reshape(q_val.Write(), VDIM, NQ1D, NQ1D, NF);
249 auto der = q_layout == QVectorLayout::byNODES ?
250 Reshape(q_der.Write(), NQ1D, NQ1D, VDIM, 2, NF):
251 Reshape(q_der.Write(), VDIM, 2, NQ1D, NQ1D, NF);
252 auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
253 auto nor = q_layout == QVectorLayout::byNODES ?
254 Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF):
255 Reshape(q_nor.Write(), 3, NQ1D, NQ1D, NF);
256 mfem::forall(NF, [=] MFEM_HOST_DEVICE (int f)
257 {
258 constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
259 constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
260 constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
261 real_t r_F[max_ND1D][max_ND1D][max_VDIM];
262 for (int d1 = 0; d1 < ND1D; d1++)
263 {
264 for (int d2 = 0; d2 < ND1D; d2++)
265 {
266 for (int c = 0; c < VDIM; c++)
267 {
268 r_F[d1][d2][c] = F(d1,d2,c,f);
269 }
270 }
271 }
272 if (eval_flags & VALUES)
273 {
274 real_t Bu[max_NQ1D][max_ND1D][max_VDIM];
275 for (int d2 = 0; d2 < ND1D; ++d2)
276 {
277 for (int q = 0; q < NQ1D; ++q)
278 {
279 for (int c = 0; c < VDIM; c++) { Bu[q][d2][c] = 0.0; }
280 for (int d1 = 0; d1 < ND1D; ++d1)
281 {
282 const real_t b = B(q,d1);
283 for (int c = 0; c < VDIM; c++)
284 {
285 Bu[q][d2][c] += b*r_F[d1][d2][c];
286 }
287 }
288 }
289 }
290 real_t BBu[max_NQ1D][max_NQ1D][max_VDIM];
291 for (int q2 = 0; q2 < NQ1D; ++q2)
292 {
293 for (int q1 = 0; q1 < NQ1D; ++q1)
294 {
295 for (int c = 0; c < VDIM; c++) { BBu[q2][q1][c] = 0.0; }
296 for (int d2 = 0; d2 < ND1D; ++d2)
297 {
298 const real_t b = B(q2,d2);
299 for (int c = 0; c < VDIM; c++)
300 {
301 BBu[q2][q1][c] += b*Bu[q1][d2][c];
302 }
303 }
304 for (int c = 0; c < VDIM; c++)
305 {
306 const real_t v = BBu[q2][q1][c];
307 if (q_layout == QVectorLayout::byVDIM) { val(c,q1,q2,f) = v; }
308 if (q_layout == QVectorLayout::byNODES) { val(q1,q2,c,f) = v; }
309 }
310 }
311 }
312 }
313 if ((eval_flags & DERIVATIVES)
314 || (eval_flags & DETERMINANTS)
315 || (eval_flags & NORMALS))
316 {
317 // We only compute the tangential derivatives
318 real_t Gu[max_NQ1D][max_ND1D][max_VDIM];
319 real_t Bu[max_NQ1D][max_ND1D][max_VDIM];
320 for (int d2 = 0; d2 < ND1D; ++d2)
321 {
322 for (int q = 0; q < NQ1D; ++q)
323 {
324 for (int c = 0; c < VDIM; c++)
325 {
326 Gu[q][d2][c] = 0.0;
327 Bu[q][d2][c] = 0.0;
328 }
329 for (int d1 = 0; d1 < ND1D; ++d1)
330 {
331 const real_t b = B(q,d1);
332 const real_t g = G(q,d1);
333 for (int c = 0; c < VDIM; c++)
334 {
335 const real_t u = r_F[d1][d2][c];
336 Gu[q][d2][c] += g*u;
337 Bu[q][d2][c] += b*u;
338 }
339 }
340 }
341 }
342 real_t BGu[max_NQ1D][max_NQ1D][max_VDIM];
343 real_t GBu[max_NQ1D][max_NQ1D][max_VDIM];
344 for (int q2 = 0; q2 < NQ1D; ++q2)
345 {
346 for (int q1 = 0; q1 < NQ1D; ++q1)
347 {
348 for (int c = 0; c < VDIM; c++)
349 {
350 BGu[q2][q1][c] = 0.0;
351 GBu[q2][q1][c] = 0.0;
352 }
353 for (int d2 = 0; d2 < ND1D; ++d2)
354 {
355 const real_t b = B(q2,d2);
356 const real_t g = G(q2,d2);
357 for (int c = 0; c < VDIM; c++)
358 {
359 BGu[q2][q1][c] += b*Gu[q1][d2][c];
360 GBu[q2][q1][c] += g*Bu[q1][d2][c];
361 }
362 }
363 }
364 }
365 if (eval_flags & DERIVATIVES)
366 {
367 for (int c = 0; c < VDIM; c++)
368 {
369 for (int q2 = 0; q2 < NQ1D; ++q2)
370 {
371 for (int q1 = 0; q1 < NQ1D; ++q1)
372 {
373 if (q_layout == QVectorLayout::byVDIM)
374 {
375 der(c,0,q1,q2,f) = BGu[q2][q1][c];
376 der(c,1,q1,q2,f) = GBu[q2][q1][c];
377 }
378 else // q_layout == QVectorLayout::byNODES
379 {
380 der(q1,q2,c,0,f) = BGu[q2][q1][c];
381 der(q1,q2,c,1,f) = GBu[q2][q1][c];
382 }
383 }
384 }
385 }
386 }
387 if (VDIM == 3 && ((eval_flags & NORMALS) ||
388 (eval_flags & DETERMINANTS)))
389 {
390 real_t n[3];
391 for (int q2 = 0; q2 < NQ1D; ++q2)
392 {
393 for (int q1 = 0; q1 < NQ1D; ++q1)
394 {
395 const real_t s = sign[f] ? -1.0 : 1.0;
396 n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
397 BGu[q2][q1][2] );
398 n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
399 BGu[q2][q1][2] );
400 n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
401 BGu[q2][q1][1] );
402 const real_t norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
403 if (eval_flags & DETERMINANTS) { det(q1,q2,f) = norm; }
404 if (eval_flags & NORMALS)
405 {
406 if (q_layout == QVectorLayout::byVDIM)
407 {
408 nor(0,q1,q2,f) = n[0]/norm;
409 nor(1,q1,q2,f) = n[1]/norm;
410 nor(2,q1,q2,f) = n[2]/norm;
411 }
412 if (q_layout == QVectorLayout::byNODES)
413 {
414 nor(q1,q2,0,f) = n[0]/norm;
415 nor(q1,q2,1,f) = n[1]/norm;
416 nor(q1,q2,2,f) = n[2]/norm;
417 }
418 }
419 }
420 }
421 }
422 }
423 });
424}
425
426template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
427void FaceQuadratureInterpolator::SmemEval3D(
428 const int NF,
429 const int vdim,
430 const QVectorLayout q_layout,
431 const DofToQuad &maps,
432 const Array<bool> &signs,
433 const Vector &e_vec,
434 Vector &q_val,
435 Vector &q_der,
436 Vector &q_det,
437 Vector &q_nor,
438 const int eval_flags)
439{
440 MFEM_PERF_SCOPE("FaceQuadInterpolator::SmemEval3D");
441 const int nd1d = maps.ndof;
442 const int nq1d = maps.nqpt;
443 const int ND1D = T_ND1D ? T_ND1D : nd1d;
444 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
445 const int VDIM = T_VDIM ? T_VDIM : vdim;
446 MFEM_VERIFY(ND1D <= MAX_ND1D, "");
447 MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
448 MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
449 auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
450 auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
451 auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
452 auto sign = signs.Read();
453 auto val = q_layout == QVectorLayout::byNODES ?
454 Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF):
455 Reshape(q_val.Write(), VDIM, NQ1D, NQ1D, NF);
456 auto der = q_layout == QVectorLayout::byNODES ?
457 Reshape(q_der.Write(), NQ1D, NQ1D, VDIM, 2, NF):
458 Reshape(q_der.Write(), VDIM, 2, NQ1D, NQ1D, NF);
459 auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
460 auto nor = q_layout == QVectorLayout::byNODES ?
461 Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF):
462 Reshape(q_nor.Write(), 3, NQ1D, NQ1D, NF);
463
464 mfem::forall_3D(NF, NQ1D, NQ1D, VDIM, [=] MFEM_HOST_DEVICE (int f)
465 {
466 constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
467 constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
468 constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
469
470 MFEM_SHARED real_t sm1[max_NQ1D*max_NQ1D*max_VDIM];
471 MFEM_SHARED real_t sm2[max_NQ1D*max_ND1D*max_VDIM];
472
473 auto s_F = (real_t(*)[max_ND1D][max_VDIM])sm1;
474 MFEM_FOREACH_THREAD(d1,x,ND1D)
475 {
476 MFEM_FOREACH_THREAD(d2,y,ND1D)
477 {
478 MFEM_FOREACH_THREAD(c,z,VDIM)
479 {
480 s_F[d1][d2][c] = F(d1,d2,c,f);
481 }
482 }
483 }
484 MFEM_SYNC_THREAD;
485
486 if (eval_flags & VALUES)
487 {
488 auto Bu = (real_t (*)[max_ND1D][max_VDIM])sm2;
489 MFEM_FOREACH_THREAD(d2,x,ND1D)
490 {
491 MFEM_FOREACH_THREAD(q1,y,NQ1D)
492 {
493 MFEM_FOREACH_THREAD(c,z,VDIM)
494 {
495 real_t thrdBu = 0.0;
496 for (int d1 = 0; d1 < ND1D; ++d1)
497 {
498 thrdBu += B(q1,d1)*s_F[d1][d2][c];
499 }
500 Bu[q1][d2][c] = thrdBu;
501 }
502 }
503 }
504 MFEM_SYNC_THREAD;
505
506 MFEM_FOREACH_THREAD(q2,x,NQ1D)
507 {
508 MFEM_FOREACH_THREAD(q1,y,NQ1D)
509 {
510 MFEM_FOREACH_THREAD(c,z,VDIM)
511 {
512 real_t v = 0.0;
513 for (int d2 = 0; d2 < ND1D; ++d2)
514 {
515 v += B(q2,d2)*Bu[q1][d2][c];
516 }
517 if (q_layout == QVectorLayout::byVDIM) { val(c,q1,q2,f) = v; }
518 if (q_layout == QVectorLayout::byNODES) { val(q1,q2,c,f) = v; }
519 }
520 }
521 }
522 }
523
524 if ((eval_flags & DERIVATIVES)
525 || (eval_flags & DETERMINANTS)
526 || (eval_flags & NORMALS))
527 {
528 // We only compute the tangential derivatives
529 auto Gu = (real_t (*)[max_ND1D][max_VDIM])sm2;
530 MFEM_SHARED real_t Bu[max_NQ1D][max_ND1D][max_VDIM];
531 MFEM_FOREACH_THREAD(d2,x,ND1D)
532 {
533 MFEM_FOREACH_THREAD(q1,y,NQ1D)
534 {
535 MFEM_FOREACH_THREAD(c,z,VDIM)
536 {
537 real_t thrdGu = 0;
538 real_t thrdBu = 0;
539 for (int d1 = 0; d1 < ND1D; ++d1)
540 {
541 const real_t u = s_F[d1][d2][c];
542 thrdBu += B(q1,d1)*u;
543 thrdGu += G(q1,d1)*u;
544 }
545 Gu[q1][d2][c] = thrdGu;
546 Bu[q1][d2][c] = thrdBu;
547 }
548 }
549 }
550 MFEM_SYNC_THREAD;
551
552 auto BGu = (real_t (*)[max_NQ1D][max_VDIM])sm1;
553 MFEM_SHARED real_t GBu[max_NQ1D][max_NQ1D][max_VDIM];
554 MFEM_FOREACH_THREAD(q2,x,NQ1D)
555 {
556 MFEM_FOREACH_THREAD(q1,y,NQ1D)
557 {
558 MFEM_FOREACH_THREAD(c,z,VDIM)
559 {
560 real_t thrdBGu = 0.0;
561 real_t thrdGBu = 0.0;
562 for (int d2 = 0; d2 < ND1D; ++d2)
563 {
564 thrdBGu += B(q2,d2)*Gu[q1][d2][c];
565 thrdGBu += G(q2,d2)*Bu[q1][d2][c];
566 }
567 BGu[q2][q1][c] = thrdBGu;
568 GBu[q2][q1][c] = thrdGBu;
569 }
570 }
571 }
572 MFEM_SYNC_THREAD;
573
574 if (eval_flags & DERIVATIVES)
575 {
576 MFEM_FOREACH_THREAD(q2,x,NQ1D)
577 {
578 MFEM_FOREACH_THREAD(q1,y,NQ1D)
579 {
580 MFEM_FOREACH_THREAD(c,z,VDIM)
581 {
582 if (q_layout == QVectorLayout::byVDIM)
583 {
584 der(c,0,q1,q2,f) = BGu[q2][q1][c];
585 der(c,1,q1,q2,f) = GBu[q2][q1][c];
586 }
587 else // q_layout == QVectorLayout::byNODES
588 {
589 der(q1,q2,c,0,f) = BGu[q2][q1][c];
590 der(q1,q2,c,1,f) = GBu[q2][q1][c];
591 }
592 }
593 }
594 }
595 }
596
597 if (VDIM == 3 && ((eval_flags & NORMALS) ||
598 (eval_flags & DETERMINANTS)))
599 {
600 real_t n[3];
601 MFEM_FOREACH_THREAD(q2,x,NQ1D)
602 {
603 MFEM_FOREACH_THREAD(q1,y,NQ1D)
604 {
605 if (MFEM_THREAD_ID(z) == 0)
606 {
607 const real_t s = sign[f] ? -1.0 : 1.0;
608 n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
609 BGu[q2][q1][2] );
610 n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
611 BGu[q2][q1][2] );
612 n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
613 BGu[q2][q1][1] );
614
615 const real_t norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
616
617 if (eval_flags & DETERMINANTS) { det(q1,q2,f) = norm; }
618
619 if (eval_flags & NORMALS)
620 {
621 if (q_layout == QVectorLayout::byVDIM)
622 {
623 nor(0,q1,q2,f) = n[0]/norm;
624 nor(1,q1,q2,f) = n[1]/norm;
625 nor(2,q1,q2,f) = n[2]/norm;
626 }
627 if (q_layout == QVectorLayout::byNODES)
628 {
629 nor(q1,q2,0,f) = n[0]/norm;
630 nor(q1,q2,1,f) = n[1]/norm;
631 nor(q1,q2,2,f) = n[2]/norm;
632 }
633 }
634 }
635 }
636 }
637 }
638 }
639 });
640}
641
642void FaceQuadratureInterpolator::Mult(
643 const Vector &e_vec, unsigned eval_flags,
644 Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
645{
646 if (nf == 0) { return; }
647 const int vdim = fespace->GetVDim();
648 const int dim = fespace->GetMesh()->Dimension();
649 const FiniteElement *fe = fespace->GetTypicalTraceElement();
650 const IntegrationRule *ir = IntRule;
651 const DofToQuad &maps = fe->GetDofToQuad(*ir, DofToQuad::TENSOR);
652 const int nd1d = maps.ndof;
653 const int nq1d = maps.nqpt;
654 void (*eval_func)(
655 const int NF,
656 const int vdim,
657 const QVectorLayout q_layout,
658 const DofToQuad &maps,
659 const Array<bool> &signs,
660 const Vector &e_vec,
661 Vector &q_val,
662 Vector &q_der,
663 Vector &q_det,
664 Vector &q_nor,
665 const int eval_flags) = NULL;
666 if (vdim == 1)
667 {
668 if (dim == 2)
669 {
670 switch (10*nd1d + nq1d)
671 {
672 // Q0
673 case 11: eval_func = &Eval2D<1,1,1>; break;
674 case 12: eval_func = &Eval2D<1,1,2>; break;
675 // Q1
676 case 22: eval_func = &Eval2D<1,2,2>; break;
677 case 23: eval_func = &Eval2D<1,2,3>; break;
678 // Q2
679 case 33: eval_func = &Eval2D<1,3,3>; break;
680 case 34: eval_func = &Eval2D<1,3,4>; break;
681 // Q3
682 case 44: eval_func = &Eval2D<1,4,4>; break;
683 case 45: eval_func = &Eval2D<1,4,5>; break;
684 case 46: eval_func = &Eval2D<1,4,6>; break;
685 // Q4
686 case 55: eval_func = &Eval2D<1,5,5>; break;
687 case 56: eval_func = &Eval2D<1,5,6>; break;
688 case 57: eval_func = &Eval2D<1,5,7>; break;
689 case 58: eval_func = &Eval2D<1,5,8>; break;
690 }
691 if (nq1d >= 10 || !eval_func)
692 {
693 eval_func = &Eval2D<1>;
694 }
695 }
696 else if (dim == 3)
697 {
698 switch (10*nd1d + nq1d)
699 {
700 // Q0
701 case 11: eval_func = &SmemEval3D<1,1,1>; break;
702 case 12: eval_func = &SmemEval3D<1,1,2>; break;
703 // Q1
704 case 22: eval_func = &SmemEval3D<1,2,2>; break;
705 case 23: eval_func = &SmemEval3D<1,2,3>; break;
706 case 24: eval_func = &SmemEval3D<1,2,4>; break;
707 // Q2
708 case 33: eval_func = &SmemEval3D<1,3,3>; break;
709 case 34: eval_func = &SmemEval3D<1,3,4>; break;
710 // Q3
711 case 44: eval_func = &SmemEval3D<1,4,4>; break;
712 case 45: eval_func = &SmemEval3D<1,4,5>; break;
713 case 46: eval_func = &SmemEval3D<1,4,6>; break;
714 // Q4
715 case 55: eval_func = &SmemEval3D<1,5,5>; break;
716 case 56: eval_func = &SmemEval3D<1,5,6>; break;
717 }
718 if (nq1d >= 10 || !eval_func)
719 {
720 eval_func = &Eval3D<1>;
721 }
722 }
723 }
724 else if (vdim == dim)
725 {
726 if (dim == 2)
727 {
728 switch (10*nd1d + nq1d)
729 {
730 // Q1
731 case 22: eval_func = &Eval2D<2,2,2>; break;
732 case 23: eval_func = &Eval2D<2,2,3>; break;
733 // Q2
734 case 33: eval_func = &Eval2D<2,3,3>; break;
735 case 34: eval_func = &Eval2D<2,3,4>; break;
736 // Q3
737 case 44: eval_func = &Eval2D<2,4,4>; break;
738 case 45: eval_func = &Eval2D<2,4,5>; break;
739 case 46: eval_func = &Eval2D<2,4,6>; break;
740 // Q4
741 case 55: eval_func = &Eval2D<2,5,5>; break;
742 case 56: eval_func = &Eval2D<2,5,6>; break;
743 case 57: eval_func = &Eval2D<2,5,7>; break;
744 case 58: eval_func = &Eval2D<2,5,8>; break;
745 }
746 if (nq1d >= 10 || !eval_func)
747 {
748 eval_func = &Eval2D<2>;
749 }
750 }
751 else if (dim == 3)
752 {
753 switch (10*nd1d + nq1d)
754 {
755 // Q1
756 case 22: eval_func = &SmemEval3D<3,2,2>; break;
757 case 23: eval_func = &SmemEval3D<3,2,3>; break;
758 case 24: eval_func = &SmemEval3D<3,2,4>; break;
759 // Q2
760 case 33: eval_func = &SmemEval3D<3,3,3>; break;
761 case 34: eval_func = &SmemEval3D<3,3,4>; break;
762 // Q3
763 case 44: eval_func = &SmemEval3D<3,4,4>; break;
764 case 45: eval_func = &SmemEval3D<3,4,5>; break;
765 case 46: eval_func = &SmemEval3D<3,4,6>; break;
766 // Q4
767 case 55: eval_func = &SmemEval3D<3,5,5>; break;
768 case 56: eval_func = &SmemEval3D<3,5,6>; break;
769 }
770 if (nq1d >= 10 || !eval_func)
771 {
772 eval_func = &Eval3D<3>;
773 }
774 }
775 }
776 if (eval_func)
777 {
778 eval_func(nf, vdim, q_layout, maps, signs, e_vec,
779 q_val, q_der, q_det, q_nor, eval_flags);
780 }
781 else
782 {
783 MFEM_ABORT("case not supported yet");
784 }
785}
786
787void FaceQuadratureInterpolator::Values(
788 const Vector &e_vec, Vector &q_val) const
789{
790 Vector q_der, q_det, q_nor;
791 Mult(e_vec, VALUES, q_val, q_der, q_det, q_nor);
792}
793
794} // 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:337
mfem::BasisType::GaussLobatto
@ GaussLobatto
Closed type.
Definition fe_base.hpp:36
mfem::BasisType::Positive
@ Positive
Bernstein polynomials.
Definition fe_base.hpp:37
mfem::DofToQuad
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:141
mfem::DofToQuad::G
Array< real_t > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition fe_base.hpp:214
mfem::DofToQuad::TENSOR
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
mfem::DofToQuad::B
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
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:178
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
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:787
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:221
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:642
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:427
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:244
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:891
mfem::FiniteElementSpace::GetTypicalTraceElement
const FiniteElement * GetTypicalTraceElement() const
Return a "typical" trace element.
Definition fespace.cpp:3959
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns the vector dimension of the finite element space.
Definition fespace.hpp:841
mfem::FiniteElementSpace::GetTypicalFE
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:3871
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:244
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::Dimension
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
mfem::ScalarFiniteElement
Class for finite elements with basis functions that return scalar values.
Definition fe_base.hpp:662
mfem::Vector
Vector data type.
Definition vector.hpp:82
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:494
mfem::Vector::Write
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:502
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:774
mfem::QVectorLayout
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition fespace.hpp:53
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:753
mfem::FaceType
FaceType
Definition mesh.hpp:48
quadinterpolator_face.hpp
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0