13 #include "../general/forall.hpp"
14 #include "../linalg/dtensor.hpp"
15 #include "../linalg/kernels.hpp"
21 static void GetSigns(
const FiniteElementSpace &fes,
const FaceType type,
24 const int dim = fes.GetMesh()->SpaceDimension();
29 for (
int f = 0;
f < fes.GetNF(); ++
f)
31 fes.GetMesh()->GetFaceElements(
f, &e1, &e2);
32 fes.GetMesh()->GetFaceInfos(
f, &inf1, &inf2);
39 if (face_id==2 || face_id==3)
50 if (face_id==0 || face_id==3 || face_id==4)
67 : type(type_), nf(fes.GetNFbyType(type)), signs(nf)
74 GetSigns(*
fespace, type, signs);
80 MFEM_VERIFY(sfe != NULL,
"Only scalar finite elements are supported");
81 MFEM_VERIFY(tfe != NULL &&
84 "Only Gauss-Lobatto and Bernstein basis are supported in "
85 "FaceQuadratureInterpolator.");
88 template<const
int T_VDIM, const
int T_ND1D, const
int T_NQ1D>
101 const int nd1d = maps.
ndof;
102 const int nq1d = maps.
nqpt;
103 const int ND1D = T_ND1D ? T_ND1D : nd1d;
104 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
105 const int VDIM = T_VDIM ? T_VDIM : vdim;
108 MFEM_VERIFY(VDIM == 2 || !(eval_flags &
DETERMINANTS),
"");
112 auto sign = signs.
Read();
118 "Derivatives on the faces are not yet supported.");
122 const int ND1D = T_ND1D ? T_ND1D : nd1d;
123 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
124 const int VDIM = T_VDIM ? T_VDIM : vdim;
125 constexpr
int max_ND1D = T_ND1D ? T_ND1D :
MAX_ND1D;
126 constexpr
int max_VDIM = T_VDIM ? T_VDIM :
MAX_VDIM2D;
127 double r_F[max_ND1D][max_VDIM];
128 for (
int d = 0; d < ND1D; d++)
130 for (
int c = 0; c < VDIM; c++)
132 r_F[d][c] = F(d,c,
f);
135 for (
int q = 0; q < NQ1D; ++q)
140 for (
int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
141 for (
int d = 0; d < ND1D; ++d)
143 const double b = B(q,d);
144 for (
int c = 0; c < VDIM; c++) { ed[c] += b*r_F[d][c]; }
146 for (
int c = 0; c < VDIM; c++) { val(q,c,
f) = ed[c]; }
149 || (eval_flags & DETERMINANTS)
153 for (
int i = 0; i < VDIM; i++) { D[i] = 0.0; }
154 for (
int d = 0; d < ND1D; ++d)
156 const double w = G(q,d);
157 for (
int c = 0; c < VDIM; c++)
159 double s_e = r_F[d][c];
164 ((eval_flags & NORMALS)
165 || (eval_flags & DETERMINANTS)))
167 const double norm = sqrt(D[0]*D[0]+D[1]*D[1]);
168 if (eval_flags & DETERMINANTS)
172 if (eval_flags & NORMALS)
174 const double s = sign[
f] ? -1.0 : 1.0;
175 n(q,0,
f) = s*D[1]/norm;
176 n(q,1,
f) = -s*D[0]/norm;
184 template<const
int T_VDIM, const
int T_ND1D, const
int T_NQ1D>
195 const int eval_flags)
197 const int nd1d = maps.
ndof;
198 const int nq1d = maps.
nqpt;
199 const int ND1D = T_ND1D ? T_ND1D : nd1d;
200 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
201 const int VDIM = T_VDIM ? T_VDIM : vdim;
204 MFEM_VERIFY(VDIM == 3 || !(eval_flags &
DETERMINANTS),
"");
207 auto F =
Reshape(e_vec.
Read(), ND1D, ND1D, VDIM, NF);
208 auto sign = signs.
Read();
209 auto val =
Reshape(q_val.
Write(), NQ1D, NQ1D, VDIM, NF);
214 "Derivatives on the faces are not yet supported.");
217 const int ND1D = T_ND1D ? T_ND1D : nd1d;
218 const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
219 const int VDIM = T_VDIM ? T_VDIM : vdim;
220 constexpr
int max_ND1D = T_ND1D ? T_ND1D :
MAX_ND1D;
221 constexpr
int max_NQ1D = T_NQ1D ? T_NQ1D :
MAX_NQ1D;
222 constexpr
int max_VDIM = T_VDIM ? T_VDIM :
MAX_VDIM3D;
223 double r_F[max_ND1D][max_ND1D][max_VDIM];
224 for (
int d1 = 0; d1 < ND1D; d1++)
226 for (
int d2 = 0; d2 < ND1D; d2++)
228 for (
int c = 0; c < VDIM; c++)
230 r_F[d1][d2][c] = F(d1,d2,c,
f);
236 double Bu[max_NQ1D][max_ND1D][max_VDIM];
237 for (
int d2 = 0; d2 < ND1D; ++d2)
239 for (
int q = 0; q < NQ1D; ++q)
241 for (
int c = 0; c < VDIM; c++) { Bu[q][d2][c] = 0.0; }
242 for (
int d1 = 0; d1 < ND1D; ++d1)
244 const double b = B(q,d1);
245 for (
int c = 0; c < VDIM; c++)
247 Bu[q][d2][c] += b*r_F[d1][d2][c];
252 double BBu[max_NQ1D][max_NQ1D][max_VDIM];
253 for (
int q2 = 0; q2 < NQ1D; ++q2)
255 for (
int q1 = 0; q1 < NQ1D; ++q1)
257 for (
int c = 0; c < VDIM; c++) { BBu[q2][q1][c] = 0.0; }
258 for (
int d2 = 0; d2 < ND1D; ++d2)
260 const double b = B(q2,d2);
261 for (
int c = 0; c < VDIM; c++)
263 BBu[q2][q1][c] += b*Bu[q1][d2][c];
266 for (
int c = 0; c < VDIM; c++)
268 val(q1,q2,c,
f) = BBu[q2][q1][c];
274 || (eval_flags & DETERMINANTS)
278 double Gu[max_NQ1D][max_ND1D][max_VDIM];
279 double Bu[max_NQ1D][max_ND1D][max_VDIM];
280 for (
int d2 = 0; d2 < ND1D; ++d2)
282 for (
int q = 0; q < NQ1D; ++q)
284 for (
int c = 0; c < VDIM; c++)
289 for (
int d1 = 0; d1 < ND1D; ++d1)
291 const double b = B(q,d1);
292 const double g = G(q,d1);
293 for (
int c = 0; c < VDIM; c++)
295 const double u = r_F[d1][d2][c];
302 double BGu[max_NQ1D][max_NQ1D][max_VDIM];
303 double GBu[max_NQ1D][max_NQ1D][max_VDIM];
304 for (
int q2 = 0; q2 < NQ1D; ++q2)
306 for (
int q1 = 0; q1 < NQ1D; ++q1)
308 for (
int c = 0; c < VDIM; c++)
310 BGu[q2][q1][c] = 0.0;
311 GBu[q2][q1][c] = 0.0;
313 for (
int d2 = 0; d2 < ND1D; ++d2)
315 const double b = B(q2,d2);
316 const double g = G(q2,d2);
317 for (
int c = 0; c < VDIM; c++)
319 BGu[q2][q1][c] += b*Gu[q1][d2][c];
320 GBu[q2][q1][c] += g*Bu[q1][d2][c];
325 if (VDIM == 3 && ((eval_flags & NORMALS) ||
326 (eval_flags & DETERMINANTS)))
329 for (
int q2 = 0; q2 < NQ1D; ++q2)
331 for (
int q1 = 0; q1 < NQ1D; ++q1)
333 const double s = sign[
f] ? -1.0 : 1.0;
334 n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
336 n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
338 n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
340 const double norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
341 if (eval_flags & DETERMINANTS) { det(q1,q2,
f) = norm; }
342 if (eval_flags & NORMALS)
344 nor(q1,q2,0,
f) = n[0]/norm;
345 nor(q1,q2,1,
f) = n[1]/norm;
346 nor(q1,q2,2,
f) = n[2]/norm;
356 const Vector &e_vec,
unsigned eval_flags,
359 if (nf == 0) {
return; }
366 const int nd1d = maps.
ndof;
367 const int nq1d = maps.
nqpt;
378 const int eval_flags) = NULL;
383 switch (10*nd1d + nq1d)
386 case 11: eval_func = &Eval2D<1,1,1>;
break;
387 case 12: eval_func = &Eval2D<1,1,2>;
break;
389 case 22: eval_func = &Eval2D<1,2,2>;
break;
390 case 23: eval_func = &Eval2D<1,2,3>;
break;
392 case 33: eval_func = &Eval2D<1,3,3>;
break;
393 case 34: eval_func = &Eval2D<1,3,4>;
break;
395 case 44: eval_func = &Eval2D<1,4,4>;
break;
396 case 45: eval_func = &Eval2D<1,4,5>;
break;
397 case 46: eval_func = &Eval2D<1,4,6>;
break;
399 case 55: eval_func = &Eval2D<1,5,5>;
break;
400 case 56: eval_func = &Eval2D<1,5,6>;
break;
401 case 57: eval_func = &Eval2D<1,5,7>;
break;
402 case 58: eval_func = &Eval2D<1,5,8>;
break;
404 if (nq1d >= 10 || !eval_func)
406 eval_func = &Eval2D<1>;
411 switch (10*nd1d + nq1d)
414 case 11: eval_func = &Eval3D<1,1,1>;
break;
415 case 12: eval_func = &Eval3D<1,1,2>;
break;
417 case 22: eval_func = &Eval3D<1,2,2>;
break;
418 case 23: eval_func = &Eval3D<1,2,3>;
break;
420 case 33: eval_func = &Eval3D<1,3,3>;
break;
421 case 34: eval_func = &Eval3D<1,3,4>;
break;
423 case 44: eval_func = &Eval3D<1,4,4>;
break;
424 case 45: eval_func = &Eval3D<1,4,5>;
break;
425 case 46: eval_func = &Eval3D<1,4,6>;
break;
427 case 55: eval_func = &Eval3D<1,5,5>;
break;
428 case 56: eval_func = &Eval3D<1,5,6>;
break;
430 if (nq1d >= 10 || !eval_func)
432 eval_func = &Eval3D<1>;
436 else if (vdim == dim)
440 switch (10*nd1d + nq1d)
443 case 22: eval_func = &Eval2D<2,2,2>;
break;
444 case 23: eval_func = &Eval2D<2,2,3>;
break;
446 case 33: eval_func = &Eval2D<2,3,3>;
break;
447 case 34: eval_func = &Eval2D<2,3,4>;
break;
449 case 44: eval_func = &Eval2D<2,4,4>;
break;
450 case 45: eval_func = &Eval2D<2,4,5>;
break;
451 case 46: eval_func = &Eval2D<2,4,6>;
break;
453 case 55: eval_func = &Eval2D<2,5,5>;
break;
454 case 56: eval_func = &Eval2D<2,5,6>;
break;
455 case 57: eval_func = &Eval2D<2,5,7>;
break;
456 case 58: eval_func = &Eval2D<2,5,8>;
break;
458 if (nq1d >= 10 || !eval_func)
460 eval_func = &Eval2D<2>;
465 switch (10*nd1d + nq1d)
468 case 22: eval_func = &Eval3D<3,2,2>;
break;
469 case 23: eval_func = &Eval3D<3,2,3>;
break;
471 case 33: eval_func = &Eval3D<3,3,3>;
break;
472 case 34: eval_func = &Eval3D<3,3,4>;
break;
474 case 44: eval_func = &Eval3D<3,4,4>;
break;
475 case 45: eval_func = &Eval3D<3,4,5>;
break;
476 case 46: eval_func = &Eval3D<3,4,6>;
break;
478 case 55: eval_func = &Eval3D<3,5,5>;
break;
479 case 56: eval_func = &Eval3D<3,5,6>;
break;
481 if (nq1d >= 10 || !eval_func)
483 eval_func = &Eval3D<3>;
489 eval_func(nf, vdim, maps, signs, e_vec,
490 q_val, q_der, q_det, q_nor, eval_flags);
494 MFEM_ABORT(
"case not supported yet");
501 Vector q_der, q_det, q_nor;
502 Mult(e_vec,
VALUES, q_val, q_der, q_det, q_nor);
Abstract class for all finite elements.
Class for an integration rule - an Array of IntegrationPoint.
static void Eval3D(const int NF, const int vdim, 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.
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
FaceQuadratureInterpolator(const FiniteElementSpace &fes, const IntegrationRule &ir, FaceType type)
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element 'i' to the given 'geom_type'.
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Geometry::Type GetFaceBaseGeometry(int i) const
const IntegrationRule * IntRule
Not owned.
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
int GetNE() const
Returns number of elements in the mesh.
Class for finite elements with basis functions that return scalar values.
double f(const Vector &xvec)
static const int MAX_ND1D
Mesh * GetMesh() const
Returns the mesh.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
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.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
int GetVDim() const
Returns vector dimension.
void Values(const Vector &e_vec, Vector &q_val) const
Interpolate the values of the E-vector e_vec at quadrature points.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
static const int MAX_VDIM3D
Evaluate the values at quadrature points.
Array< double > B
Basis functions evaluated at quadrature points.
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...
static const int MAX_NQ1D
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
static void Eval2D(const int NF, const int vdim, 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.
Evaluate the derivatives at quadrature points.
void f_vec(const Vector &xvec, Vector &f)
double u(const Vector &xvec)
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
const FiniteElementSpace * fespace
Not owned.
static const int MAX_VDIM2D
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).