MFEM  v4.1.0 Finite element discretization library
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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
13 #include "../general/forall.hpp"
14 #include "../linalg/dtensor.hpp"
15 #include "../linalg/kernels.hpp"
16
17 namespace mfem
18 {
19
20 /// Return the sign to apply to the normals on each face to point from e1 to e2.
21 static void GetSigns(const FiniteElementSpace &fes, const FaceType type,
22  Array<bool> &signs)
23 {
24  const int dim = fes.GetMesh()->SpaceDimension();
25  int e1, e2;
26  int inf1, inf2;
27  int face_id;
28  int f_ind = 0;
29  for (int f = 0; f < fes.GetNF(); ++f)
30  {
31  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
32  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
33  face_id = inf1 / 64;
34  if ( (type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
35  (type==FaceType::Boundary && e2<0 && inf2<0) )
36  {
37  if (dim==2)
38  {
39  if (face_id==2 || face_id==3)
40  {
41  signs[f_ind] = true;
42  }
43  else
44  {
45  signs[f_ind] = false;
46  }
47  }
48  else if (dim==3)
49  {
50  if (face_id==0 || face_id==3 || face_id==4)
51  {
52  signs[f_ind] = true;
53  }
54  else
55  {
56  signs[f_ind] = false;
57  }
58  }
59  f_ind++;
60  }
61  }
62 }
63
65  const FiniteElementSpace &fes,
66  const IntegrationRule &ir, FaceType type_)
67  : type(type_), nf(fes.GetNFbyType(type)), signs(nf)
68 {
69  fespace = &fes;
70  IntRule = &ir;
71  use_tensor_products = true; // not implemented yet (not used)
72
73  if (fespace->GetNE() == 0) { return; }
74  GetSigns(*fespace, type, signs);
75  const FiniteElement *fe = fespace->GetFE(0);
76  const ScalarFiniteElement *sfe =
77  dynamic_cast<const ScalarFiniteElement*>(fe);
78  const TensorBasisElement *tfe =
79  dynamic_cast<const TensorBasisElement*>(fe);
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 "
86 }
87
88 template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
90  const int NF,
91  const int vdim,
93  const Array<bool> &signs,
94  const Vector &f_vec,
95  Vector &q_val,
96  Vector &q_der,
97  Vector &q_det,
98  Vector &q_nor,
99  const int eval_flags)
100 {
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;
106  MFEM_VERIFY(ND1D <= MAX_ND1D, "");
107  MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
108  MFEM_VERIFY(VDIM == 2 || !(eval_flags & DETERMINANTS), "");
109  auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
110  auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
111  auto F = Reshape(f_vec.Read(), ND1D, VDIM, NF);
113  auto val = Reshape(q_val.Write(), NQ1D, VDIM, NF);
114  // auto der = Reshape(q_der.Write(), NQ1D, VDIM, NF); // only tangential der
115  auto det = Reshape(q_det.Write(), NQ1D, NF);
116  auto n = Reshape(q_nor.Write(), NQ1D, VDIM, NF);
117  MFEM_VERIFY(eval_flags | DERIVATIVES,
118  "Derivatives on the faces are not yet supported.");
119  // If Gauss-Lobatto
120  MFEM_FORALL(f, NF,
121  {
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++)
129  {
130  for (int c = 0; c < VDIM; c++)
131  {
132  r_F[d][c] = F(d,c,f);
133  }
134  }
135  for (int q = 0; q < NQ1D; ++q)
136  {
137  if (eval_flags & VALUES)
138  {
139  double ed[max_VDIM];
140  for (int c = 0; c < VDIM; c++) { ed[c] = 0.0; }
141  for (int d = 0; d < ND1D; ++d)
142  {
143  const double b = B(q,d);
144  for (int c = 0; c < VDIM; c++) { ed[c] += b*r_F[d][c]; }
145  }
146  for (int c = 0; c < VDIM; c++) { val(q,c,f) = ed[c]; }
147  }
148  if ((eval_flags & DERIVATIVES)
149  || (eval_flags & DETERMINANTS)
150  || (eval_flags & NORMALS))
151  {
152  double D[max_VDIM];
153  for (int i = 0; i < VDIM; i++) { D[i] = 0.0; }
154  for (int d = 0; d < ND1D; ++d)
155  {
156  const double w = G(q,d);
157  for (int c = 0; c < VDIM; c++)
158  {
159  double s_e = r_F[d][c];
160  D[c] += s_e * w;
161  }
162  }
163  if (VDIM == 2 &&
164  ((eval_flags & NORMALS)
165  || (eval_flags & DETERMINANTS)))
166  {
167  const double norm = sqrt(D[0]*D[0]+D[1]*D[1]);
168  if (eval_flags & DETERMINANTS)
169  {
170  det(q,f) = norm;
171  }
172  if (eval_flags & NORMALS)
173  {
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;
177  }
178  }
179  }
180  }
181  });
182 }
183
184 template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
186  const int NF,
187  const int vdim,
189  const Array<bool> &signs,
190  const Vector &e_vec,
191  Vector &q_val,
192  Vector &q_der,
193  Vector &q_det,
194  Vector &q_nor,
195  const int eval_flags)
196 {
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;
202  MFEM_VERIFY(ND1D <= MAX_ND1D, "");
203  MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
204  MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
205  auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
206  auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
207  auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
209  auto val = Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF);
210  // auto der = Reshape(q_der.Write(), NQ1D, VDIM, 3, NF);
211  auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
212  auto nor = Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF);
213  MFEM_VERIFY(eval_flags | DERIVATIVES,
214  "Derivatives on the faces are not yet supported.");
215  MFEM_FORALL(f, NF,
216  {
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++)
225  {
226  for (int d2 = 0; d2 < ND1D; d2++)
227  {
228  for (int c = 0; c < VDIM; c++)
229  {
230  r_F[d1][d2][c] = F(d1,d2,c,f);
231  }
232  }
233  }
234  if (eval_flags & VALUES)
235  {
236  double Bu[max_NQ1D][max_ND1D][VDIM];
237  for (int d2 = 0; d2 < ND1D; ++d2)
238  {
239  for (int q = 0; q < NQ1D; ++q)
240  {
241  for (int c = 0; c < VDIM; c++) { Bu[q][d2][c] = 0.0; }
242  for (int d1 = 0; d1 < ND1D; ++d1)
243  {
244  const double b = B(q,d1);
245  for (int c = 0; c < VDIM; c++)
246  {
247  Bu[q][d2][c] += b*r_F[d1][d2][c];
248  }
249  }
250  }
251  }
252  double BBu[max_NQ1D][max_NQ1D][VDIM];
253  for (int q2 = 0; q2 < NQ1D; ++q2)
254  {
255  for (int q1 = 0; q1 < NQ1D; ++q1)
256  {
257  for (int c = 0; c < VDIM; c++) { BBu[q2][q1][c] = 0.0; }
258  for (int d2 = 0; d2 < ND1D; ++d2)
259  {
260  const double b = B(q2,d2);
261  for (int c = 0; c < VDIM; c++)
262  {
263  BBu[q2][q1][c] += b*Bu[q1][d2][c];
264  }
265  }
266  for (int c = 0; c < VDIM; c++)
267  {
268  val(q1,q2,c,f) = BBu[q2][q1][c];
269  }
270  }
271  }
272  }
273  if ((eval_flags & DERIVATIVES)
274  || (eval_flags & DETERMINANTS)
275  || (eval_flags & NORMALS))
276  {
277  // We only compute the tangential derivatives
278  double Gu[max_NQ1D][max_ND1D][VDIM];
279  double Bu[max_NQ1D][max_ND1D][VDIM];
280  for (int d2 = 0; d2 < ND1D; ++d2)
281  {
282  for (int q = 0; q < NQ1D; ++q)
283  {
284  for (int c = 0; c < VDIM; c++)
285  {
286  Gu[q][d2][c] = 0.0;
287  Bu[q][d2][c] = 0.0;
288  }
289  for (int d1 = 0; d1 < ND1D; ++d1)
290  {
291  const double b = B(q,d1);
292  const double g = G(q,d1);
293  for (int c = 0; c < VDIM; c++)
294  {
295  const double u = r_F[d1][d2][c];
296  Gu[q][d2][c] += g*u;
297  Bu[q][d2][c] += b*u;
298  }
299  }
300  }
301  }
302  double BGu[max_NQ1D][max_NQ1D][VDIM];
303  double GBu[max_NQ1D][max_NQ1D][VDIM];
304  for (int q2 = 0; q2 < NQ1D; ++q2)
305  {
306  for (int q1 = 0; q1 < NQ1D; ++q1)
307  {
308  for (int c = 0; c < VDIM; c++)
309  {
310  BGu[q2][q1][c] = 0.0;
311  GBu[q2][q1][c] = 0.0;
312  }
313  for (int d2 = 0; d2 < ND1D; ++d2)
314  {
315  const double b = B(q2,d2);
316  const double g = G(q2,d2);
317  for (int c = 0; c < VDIM; c++)
318  {
319  BGu[q2][q1][c] += b*Gu[q1][d2][c];
320  GBu[q2][q1][c] += g*Bu[q1][d2][c];
321  }
322  }
323  }
324  }
325  if (VDIM == 3 && ((eval_flags & NORMALS) ||
326  (eval_flags & DETERMINANTS)))
327  {
328  double n[3];
329  for (int q2 = 0; q2 < NQ1D; ++q2)
330  {
331  for (int q1 = 0; q1 < NQ1D; ++q1)
332  {
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]*
335  BGu[q2][q1][2] );
336  n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
337  BGu[q2][q1][2] );
338  n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
339  BGu[q2][q1][1] );
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)
343  {
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;
347  }
348  }
349  }
350  }
351  }
352  });
353 }
354
356  const Vector &e_vec, unsigned eval_flags,
357  Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
358 {
359  if (nf == 0) { return; }
360  const int vdim = fespace->GetVDim();
361  const int dim = fespace->GetMesh()->Dimension();
362  const FiniteElement *fe =
364  const IntegrationRule *ir = IntRule;
366  const int nd = maps.ndof;
367  const int nq = maps.nqpt;
368  void (*eval_func)(
369  const int NF,
370  const int vdim,
372  const Array<bool> &signs,
373  const Vector &e_vec,
374  Vector &q_val,
375  Vector &q_der,
376  Vector &q_det,
377  Vector &q_nor,
378  const int eval_flags) = NULL;
379  if (vdim == 1)
380  {
381  if (dim == 2)
382  {
383  switch (100*nd + nq)
384  {
385  // Q0
386  case 101: eval_func = &Eval2D<1,1,1>; break;
387  case 104: eval_func = &Eval2D<1,1,4>; break;
388  // Q1
389  case 404: eval_func = &Eval2D<1,4,4>; break;
390  case 409: eval_func = &Eval2D<1,4,9>; break;
391  // Q2
392  case 909: eval_func = &Eval2D<1,9,9>; break;
393  case 916: eval_func = &Eval2D<1,9,16>; break;
394  // Q3
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;
398  // Q4
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;
403  }
404  if (nq >= 100 || !eval_func)
405  {
406  eval_func = &Eval2D<1>;
407  }
408  }
409  else if (dim == 3)
410  {
411  switch (1000*nd + nq)
412  {
413  // Q0
414  case 1001: eval_func = &Eval3D<1,1,1>; break;
415  case 1008: eval_func = &Eval3D<1,1,8>; break;
416  // Q1
417  case 8008: eval_func = &Eval3D<1,8,8>; break;
418  case 8027: eval_func = &Eval3D<1,8,27>; break;
419  // Q2
420  case 27027: eval_func = &Eval3D<1,27,27>; break;
421  case 27064: eval_func = &Eval3D<1,27,64>; break;
422  // Q3
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;
426  // Q4
427  case 125125: eval_func = &Eval3D<1,125,125>; break;
428  case 125216: eval_func = &Eval3D<1,125,216>; break;
429  }
430  if (nq >= 1000 || !eval_func)
431  {
432  eval_func = &Eval3D<1>;
433  }
434  }
435  }
436  else if (vdim == dim)
437  {
438  if (dim == 2)
439  {
440  switch (100*nd + nq)
441  {
442  // Q1
443  case 404: eval_func = &Eval2D<2,4,4>; break;
444  case 409: eval_func = &Eval2D<2,4,9>; break;
445  // Q2
446  case 909: eval_func = &Eval2D<2,9,9>; break;
447  case 916: eval_func = &Eval2D<2,9,16>; break;
448  // Q3
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;
452  // Q4
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;
457  }
458  if (nq >= 100 || !eval_func)
459  {
460  eval_func = &Eval2D<2>;
461  }
462  }
463  else if (dim == 3)
464  {
465  switch (1000*nd + nq)
466  {
467  // Q1
468  case 8008: eval_func = &Eval3D<3,8,8>; break;
469  case 8027: eval_func = &Eval3D<3,8,27>; break;
470  // Q2
471  case 27027: eval_func = &Eval3D<3,27,27>; break;
472  case 27064: eval_func = &Eval3D<3,27,64>; break;
473  // Q3
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;
477  // Q4
478  case 125125: eval_func = &Eval3D<3,125,125>; break;
479  case 125216: eval_func = &Eval3D<3,125,216>; break;
480  }
481  if (nq >= 1000 || !eval_func)
482  {
483  eval_func = &Eval3D<3>;
484  }
485  }
486  }
487  if (eval_func)
488  {
489  eval_func(nf, vdim, maps, signs, e_vec,
490  q_val, q_der, q_det, q_nor, eval_flags);
491  }
492  else
493  {
494  MFEM_ABORT("case not supported yet");
495  }
496 }
497
498 } // namespace mfem
Abstract class for Finite Elements.
Definition: fe.hpp:232
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
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...
Definition: fe.hpp:154
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 &#39;i&#39; to the given &#39;geom_type&#39;.
Definition: fespace.cpp:1936
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe.hpp:166
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe.hpp:162
Geometry::Type GetFaceBaseGeometry(int i) const
Definition: mesh.hpp:785
double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:341
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.
Definition: dtensor.hpp:134
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
int GetBasisType() const
Definition: fe.hpp:1842
FaceType
Definition: mesh.hpp:42
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
double b
Definition: lissajous.cpp:42
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
Definition: array.hpp:273
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
int Dimension() const
Definition: mesh.hpp:744
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe.hpp:128
Evaluate the values at quadrature points.
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe.hpp:177
Definition: fe.cpp:370
const double * Read(bool on_dev=true) const
Definition: vector.hpp:333
int dim
Definition: ex24.cpp:43
Array< double > G
Definition: fe.hpp:198
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
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.
Vector data type.
Definition: vector.hpp:48
Bernstein polynomials.
Definition: fe.hpp:35
Assuming the derivative at quadrature points form a matrix, this flag can be used to compute and stor...
const FiniteElementSpace * fespace
Not owned.