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");
Abstract class for 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.
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.
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
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 associated with i'th element.
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