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 nd = maps.
ndof;
102 const int nq = maps.
nqpt;
103 const int ND1D = T_ND1D ? T_ND1D : nd;
104 const int NQ1D = T_NQ1D ? T_NQ1D : nq;
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 : nd;
123 const int NQ1D = T_NQ1D ? T_NQ1D : nq;
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 nd = maps.
ndof;
198 const int nq = maps.
nqpt;
199 const int ND1D = T_ND1D ? T_ND1D : nd;
200 const int NQ1D = T_NQ1D ? T_NQ1D : nq;
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 : nd;
218 const int NQ1D = T_NQ1D ? T_NQ1D : nq;
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][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][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][VDIM];
279 double Bu[max_NQ1D][max_ND1D][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][VDIM];
303 double GBu[max_NQ1D][max_NQ1D][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 nd = maps.
ndof;
367 const int nq = maps.
nqpt;
378 const int eval_flags) = NULL;
386 case 101: eval_func = &Eval2D<1,1,1>;
break;
387 case 104: eval_func = &Eval2D<1,1,4>;
break;
389 case 404: eval_func = &Eval2D<1,4,4>;
break;
390 case 409: eval_func = &Eval2D<1,4,9>;
break;
392 case 909: eval_func = &Eval2D<1,9,9>;
break;
393 case 916: eval_func = &Eval2D<1,9,16>;
break;
395 case 1616: eval_func = &Eval2D<1,16,16>;
break;
396 case 1625: eval_func = &Eval2D<1,16,25>;
break;
397 case 1636: eval_func = &Eval2D<1,16,36>;
break;
399 case 2525: eval_func = &Eval2D<1,25,25>;
break;
400 case 2536: eval_func = &Eval2D<1,25,36>;
break;
401 case 2549: eval_func = &Eval2D<1,25,49>;
break;
402 case 2564: eval_func = &Eval2D<1,25,64>;
break;
404 if (nq >= 100 || !eval_func)
406 eval_func = &Eval2D<1>;
411 switch (1000*nd + nq)
414 case 1001: eval_func = &Eval3D<1,1,1>;
break;
415 case 1008: eval_func = &Eval3D<1,1,8>;
break;
417 case 8008: eval_func = &Eval3D<1,8,8>;
break;
418 case 8027: eval_func = &Eval3D<1,8,27>;
break;
420 case 27027: eval_func = &Eval3D<1,27,27>;
break;
421 case 27064: eval_func = &Eval3D<1,27,64>;
break;
423 case 64064: eval_func = &Eval3D<1,64,64>;
break;
424 case 64125: eval_func = &Eval3D<1,64,125>;
break;
425 case 64216: eval_func = &Eval3D<1,64,216>;
break;
427 case 125125: eval_func = &Eval3D<1,125,125>;
break;
428 case 125216: eval_func = &Eval3D<1,125,216>;
break;
430 if (nq >= 1000 || !eval_func)
432 eval_func = &Eval3D<1>;
436 else if (vdim == dim)
443 case 404: eval_func = &Eval2D<2,4,4>;
break;
444 case 409: eval_func = &Eval2D<2,4,9>;
break;
446 case 909: eval_func = &Eval2D<2,9,9>;
break;
447 case 916: eval_func = &Eval2D<2,9,16>;
break;
449 case 1616: eval_func = &Eval2D<2,16,16>;
break;
450 case 1625: eval_func = &Eval2D<2,16,25>;
break;
451 case 1636: eval_func = &Eval2D<2,16,36>;
break;
453 case 2525: eval_func = &Eval2D<2,25,25>;
break;
454 case 2536: eval_func = &Eval2D<2,25,36>;
break;
455 case 2549: eval_func = &Eval2D<2,25,49>;
break;
456 case 2564: eval_func = &Eval2D<2,25,64>;
break;
458 if (nq >= 100 || !eval_func)
460 eval_func = &Eval2D<2>;
465 switch (1000*nd + nq)
468 case 8008: eval_func = &Eval3D<3,8,8>;
break;
469 case 8027: eval_func = &Eval3D<3,8,27>;
break;
471 case 27027: eval_func = &Eval3D<3,27,27>;
break;
472 case 27064: eval_func = &Eval3D<3,27,64>;
break;
474 case 64064: eval_func = &Eval3D<3,64,64>;
break;
475 case 64125: eval_func = &Eval3D<3,64,125>;
break;
476 case 64216: eval_func = &Eval3D<3,64,216>;
break;
478 case 125125: eval_func = &Eval3D<3,125,125>;
break;
479 case 125216: eval_func = &Eval3D<3,125,216>;
break;
481 if (nq >= 1000 || !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
double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
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.
static const int MAX_ND1D
Mesh * GetMesh() const
Returns the mesh.
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...
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
static const int MAX_NQ1D
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
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.
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
double f(const Vector &p)