MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
quadinterpolator_face.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "quadinterpolator_face.hpp"
13 #include "../general/annotation.hpp"
14 #include "../general/forall.hpp"
15 #include "../linalg/dtensor.hpp"
16 #include "../linalg/kernels.hpp"
17 
18 namespace mfem
19 {
20 
21 /// Return the sign to apply to the normals on each face to point from e1 to e2.
22 static 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 
68 FaceQuadratureInterpolator::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 
93 template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
94 void 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(f, NF,
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  double 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  double 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 double 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  double 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 double w = G(q,d);
171  for (int c = 0; c < VDIM; c++)
172  {
173  double 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 double 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 double 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 
206 template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
207 void 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(f, NF,
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  double 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  double 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 double 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  double 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 double 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 double 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  double Gu[max_NQ1D][max_ND1D][max_VDIM];
305  double 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 double b = B(q,d1);
318  const double g = G(q,d1);
319  for (int c = 0; c < VDIM; c++)
320  {
321  const double 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  double BGu[max_NQ1D][max_NQ1D][max_VDIM];
329  double 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 double b = B(q2,d2);
342  const double 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  double n[3];
355  for (int q2 = 0; q2 < NQ1D; ++q2)
356  {
357  for (int q1 = 0; q1 < NQ1D; ++q1)
358  {
359  const double 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 double 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 
390 template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
391 void 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(f, NF, NQ1D, NQ1D, VDIM,
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 double sm1[max_NQ1D*max_NQ1D*max_VDIM];
435  MFEM_SHARED double sm2[max_NQ1D*max_ND1D*max_VDIM];
436 
437  auto s_F = (double(*)[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 = (double (*)[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  double 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  double 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 = (double (*)[max_ND1D][max_VDIM])sm2;
494  MFEM_SHARED double 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  double thrdGu = 0;
502  double thrdBu = 0;
503  for (int d1 = 0; d1 < ND1D; ++d1)
504  {
505  const double 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 = (double (*)[max_NQ1D][max_VDIM])sm1;
517  MFEM_SHARED double 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  double thrdBGu = 0.0;
525  double 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  double 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 double 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 double 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 
583 void 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()->GetFaceBaseGeometry(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 
729 void 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::FiniteElement
Abstract class for all finite elements.
Definition: fe_base.hpp:235
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::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::DofToQuad::TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:162
mfem::FaceQuadratureInterpolator::FaceQuadratureInterpolator
FaceQuadratureInterpolator(const FiniteElementSpace &fes, const IntegrationRule &ir, FaceType type)
Definition: quadinterpolator_face.cpp:68
mfem::FiniteElementSpace::GetTraceElement
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element &#39;i&#39; to the given &#39;geom_type&#39;.
Definition: fespace.cpp:3173
mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:174
mfem::DofToQuad::ndof
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:170
mfem::Mesh::GetFaceBaseGeometry
Geometry::Type GetFaceBaseGeometry(int i) const
Definition: mesh.hpp:1066
mfem::FaceQuadratureInterpolator::IntRule
const IntegrationRule * IntRule
Not owned.
Definition: quadinterpolator_face.hpp:35
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
mfem::ScalarFiniteElement
Class for finite elements with basis functions that return scalar values.
Definition: fe_base.hpp:629
mfem::TensorBasisElement::GetBasisType
int GetBasisType() const
Definition: fe_base.hpp:1174
mfem::FaceType
FaceType
Definition: mesh.hpp:45
mfem::f
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
mfem::Vector::Write
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:457
b
double b
Definition: lissajous.cpp:42
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::BasisType::Positive
Bernstein polynomials.
Definition: fe_base.hpp:32
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:1006
mfem::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
quadinterpolator_face.hpp
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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::DofToQuad
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:136
mfem::FaceQuadratureInterpolator::VALUES
Evaluate the values at quadrature points.
Definition: quadinterpolator_face.hpp:55
mfem::DofToQuad::B
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:185
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:366
mfem::QVectorLayout
QVectorLayout
Type describing possible layouts for Q-vectors.
Definition: fespace.hpp:52
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
dim
int dim
Definition: ex24.cpp:53
mfem::DofToQuad::G
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition: fe_base.hpp:206
mfem::FiniteElementSpace::GetFE
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2781
mfem::QVectorLayout::byNODES
NQPT x VDIM x NE (values) / NQPT x VDIM x DIM x NE (grads)
mfem::FaceQuadratureInterpolator::DERIVATIVES
Evaluate the derivatives at quadrature points.
Definition: quadinterpolator_face.hpp:56
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::f_vec
void f_vec(const Vector &xvec, Vector &f)
Definition: lor_mms.hpp:68
s
RefCoord s[3]
Definition: ncmesh_tables.hpp:158
mfem::u
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
mfem::FaceQuadratureInterpolator::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::fespace
const FiniteElementSpace * fespace
Not owned.
Definition: quadinterpolator_face.hpp:34
mfem::FaceQuadratureInterpolator::q_layout
QVectorLayout q_layout
Output Q-vector layout.
Definition: quadinterpolator_face.hpp:36
mfem::Vector::Read
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:449
mfem::Reshape
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
mfem::QVectorLayout::byVDIM
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
Generated on Wed Dec 7 2022 10:41:46 for MFEM by   doxygen 1.8.5