dune-pdelab  2.4-dev
dgnavierstokesvelvecfem.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
6 
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/parametertree.hh>
10 
11 #include <dune/geometry/quadraturerules.hh>
12 
13 #include <dune/localfunctions/common/interfaceswitch.hh>
15 
21 
22 #ifndef VBLOCK
23 #define VBLOCK 0
24 #endif
25 #define PBLOCK (- VBLOCK + 1)
26 
27 namespace Dune {
28  namespace PDELab {
29 
30  template<class Basis, class Dummy = void>
33  typedef typename Basis::Traits::DomainLocal DomainLocal;
35  typedef typename Basis::Traits::RangeField RangeField;
37  static const std::size_t dimRange = Basis::Traits::dimRange;
38 
40  template<typename Geometry>
41  static void jacobian(const Basis& basis, const Geometry& geometry,
42  const DomainLocal& xl,
43  std::vector<FieldMatrix<RangeField, dimRange,
44  Geometry::coorddimension> >& jac)
45  {
46  jac.resize(basis.size());
47  basis.evaluateJacobian(xl, jac);
48  }
49  };
50 
52  template<class Basis>
54  Basis, typename enable_if<
55  Std::to_true_type<
56  integral_constant<
57  std::size_t,
58  Basis::Traits::dimDomain
59  >
60  >::value
61  >::type
62  >
63  {
65  typedef typename Basis::Traits::DomainType DomainLocal;
67  typedef typename Basis::Traits::RangeFieldType RangeField;
69  static const std::size_t dimRange = Basis::Traits::dimRange;
70 
72  template<typename Geometry>
73  static void jacobian(const Basis& basis, const Geometry& geometry,
74  const DomainLocal& xl,
75  std::vector<FieldMatrix<RangeField, dimRange,
76  Geometry::coorddimension> >& jac)
77  {
78  jac.resize(basis.size());
79 
80  std::vector<FieldMatrix<
81  RangeField, dimRange, Geometry::coorddimension> > ljac(basis.size());
82  basis.evaluateJacobian(xl, ljac);
83 
84  const typename Geometry::JacobianInverseTransposed& geojac =
85  geometry.jacobianInverseTransposed(xl);
86 
87  for(std::size_t i = 0; i < basis.size(); ++i)
88  for(std::size_t row=0; row < dimRange; ++row)
89  geojac.mv(ljac[i][row], jac[i][row]);
90  }
91  };
92 
100  template<typename PRM>
103  public FullSkeletonPattern, public FullVolumePattern,
105  {
106  typedef StokesBoundaryCondition BC;
107  typedef typename PRM::Traits::RangeField RF;
108 
110  typedef typename InstatBase::RealType Real;
111 
112  static const bool navier = PRM::assemble_navier;
113  static const bool full_tensor = PRM::assemble_full_tensor;
114 
115  public :
116 
117  // pattern assembly flags
118  enum { doPatternVolume = true };
119  enum { doPatternSkeleton = true };
120 
121  // call the assembler for each face only once
122  enum { doSkeletonTwoSided = false };
123 
124  // residual assembly flags
125  enum { doAlphaVolume = true };
126  enum { doAlphaSkeleton = true };
127  enum { doAlphaBoundary = true };
128  enum { doLambdaVolume = true };
129 
142  DGNavierStokesVelVecFEM (PRM& _prm, int _superintegration_order=0) :
143  prm(_prm), superintegration_order(_superintegration_order),
144  current_dt(1.0)
145  {}
146 
147  // Store current dt
148  void preStep (RealType , RealType dt, int )
149  {
150  current_dt = dt;
151  }
152 
153  // set time in parameter class
154  void setTime(Real t)
155  {
157  prm.setTime(t);
158  }
159 
160  // volume integral depending on test and ansatz functions
161  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
162  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
163  {
164  const unsigned int dim = EG::Geometry::dimension;
165 
166  // subspaces
167  static_assert
168  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
169 
170  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
171  const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
172  const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
173 
174  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
175  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
176  const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
177 
178  // domain and range field type
179  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
180  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
182  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
183  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
184  typedef typename BasisSwitch_V::DomainField DF;
185  typedef typename BasisSwitch_V::RangeField RF;
186  typedef typename BasisSwitch_V::Range Range_V;
187  typedef typename BasisSwitch_P::Range Range_P;
188  typedef typename LFSV::Traits::SizeType size_type;
189 
190  // select quadrature rule
191  Dune::GeometryType gt = eg.geometry().type();
192  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
193  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
194  const int jac_order = gt.isSimplex() ? 0 : 1;
195  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
196  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
197 
198  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
199 
200  // loop over quadrature points
201  for (const auto& ip : rule)
202  {
203  const Dune::FieldVector<DF,dim> local = ip.position();
204  const RF mu = prm.mu(eg,local);
205  const RF rho = prm.rho(eg,local);
206 
207  // compute u (if Navier term enabled)
208  std::vector<Range_V> phi_v(lfsv_v.size());
209  Range_V val_u(0.0);
210  if(navier) {
211  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
212  for (size_type i=0; i<lfsu_v.size(); i++)
213  val_u.axpy(x(lfsu_v,i),phi_v[i]);
214  }
215 
216  // values of pressure shape functions
217  std::vector<Range_P> phi_p(lfsv_p.size());
218  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
219 
220  // compute pressure value
221  Range_P val_p(0.0);
222  for (size_type i=0; i<lfsu_p.size(); i++)
223  val_p.axpy(x(lfsu_p,i),phi_p[i]);
224 
225  // evaluate jacobian of velocity shape functions on reference element
226  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
227  VectorBasisSwitch_V::jacobian
228  (FESwitch_V::basis(lfsv_v.finiteElement()), eg.geometry(), local, jac_phi_v);
229 
230  // compute divergence of test functions
231  std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
232  for (size_type i=0; i<lfsv_v.size(); i++)
233  for (size_type d=0; d<dim; d++)
234  div_phi_v[i] += jac_phi_v[i][d][d];
235 
236  // compute velocity jacobian and divergence
237  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
238  RF div_u(0.0);
239  for (size_type i=0; i<lfsu_v.size(); i++){
240  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
241  div_u += x(lfsu_v,i) * div_phi_v[i];
242  }
243 
244  const RF detj = eg.geometry().integrationElement(ip.position());
245  const RF weight = ip.weight() * detj;
246 
247  for (size_type i=0; i<lfsv_v.size(); i++) {
248  //================================================//
249  // \int (mu*grad_u*grad_v)
250  //================================================//
251  RF dvdu(0); contraction(jac_u,jac_phi_v[i],dvdu);
252  r.accumulate(lfsv_v, i, dvdu * mu * weight);
253 
254  //================================================//
255  // \int -p \nabla\cdot v
256  //================================================//
257  r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
258 
259  //================================================//
260  // \int \rho ((u\cdot\nabla ) u )\cdot v
261  //================================================//
262  if(navier) {
263  // compute (grad u) u (matrix-vector product)
264  Range_V nabla_u_u(0.0);
265  jac_u.mv(val_u,nabla_u_u);
266  r.accumulate(lfsv_v, i, rho * (nabla_u_u*phi_v[i]) * weight);
267  } // end navier
268 
269  } // end i
270 
271  for (size_type i=0; i<lfsv_p.size(); i++) {
272  //================================================//
273  // \int -q \nabla\cdot u
274  //================================================//
275  r.accumulate(lfsv_p, i, - div_u * phi_p[i] * incomp_scaling * weight);
276  }
277 
278  } // end loop quadrature points
279  } // end alpha_volume
280 
281  // jacobian of volume term
282  template<typename EG, typename LFSU, typename X, typename LFSV,
283  typename LocalMatrix>
284  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
285  LocalMatrix& mat) const
286  {
287  const unsigned int dim = EG::Geometry::dimension;
288 
289  // subspaces
290  static_assert
291  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
292 
293  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
294  const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
295  const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
296 
297  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
298  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
299  const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
300 
301  // domain and range field type
302  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
303  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
305  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
306  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
307  typedef typename BasisSwitch_V::DomainField DF;
308  typedef typename BasisSwitch_V::RangeField RF;
309  typedef typename BasisSwitch_V::Range Range_V;
310  typedef typename BasisSwitch_P::Range Range_P;
311  typedef typename LFSV::Traits::SizeType size_type;
312 
313  // select quadrature rule
314  Dune::GeometryType gt = eg.geometry().type();
315  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
316  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
317  const int jac_order = gt.isSimplex() ? 0 : 1;
318  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
319  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
320 
321  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
322 
323  // loop over quadrature points
324  for (const auto& ip : rule)
325  {
326  const Dune::FieldVector<DF,dim> local = ip.position();
327  const RF mu = prm.mu(eg,local);
328  const RF rho = prm.rho(eg,local);
329 
330  // compute u (if Navier term enabled)
331  std::vector<Range_V> phi_v(lfsv_v.size());
332  Range_V val_u(0.0);
333  if(navier) {
334  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
335  for (size_type i=0; i<lfsu_v.size(); i++)
336  val_u.axpy(x(lfsu_v,i),phi_v[i]);
337  }
338 
339  // values of pressure shape functions
340  std::vector<Range_P> phi_p(lfsv_p.size());
341  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
342 
343  // evaluate jacobian of velocity shape functions on reference element
344  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
345  VectorBasisSwitch_V::jacobian
346  (FESwitch_V::basis(lfsv_v.finiteElement()), eg.geometry(), local, jac_phi_v);
347 
348  assert(lfsu_v.size() == lfsv_v.size());
349  // compute divergence of velocity shape functions
350  std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
351  for (size_type i=0; i<lfsv_v.size(); i++)
352  for (size_type d=0; d<dim; d++)
353  div_phi_v[i] += jac_phi_v[i][d][d];
354 
355  // compute velocity jacobian (if Navier term enabled)
356  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
357  if(navier) {
358  for (size_type i=0; i<lfsu_v.size(); i++){
359  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
360  }
361  }
362 
363  const RF detj = eg.geometry().integrationElement(ip.position());
364  const RF weight = ip.weight() * detj;
365 
366  for(size_type i=0; i<lfsv_v.size(); i++) {
367 
368  for(size_type j=0; j<lfsu_v.size(); j++) {
369  //================================================//
370  // \int (mu*grad_u*grad_v)
371  //================================================//
372  RF dvdu(0.0); contraction(jac_phi_v[j],jac_phi_v[i],dvdu);
373  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * dvdu * weight);
374 
375  //================================================//
376  // \int \rho ((u\cdot\nabla ) u )\cdot v
377  //================================================//
378  if(navier) {
379  // compute (grad u) phi_v_j (matrix-vector product)
380  Range_V nabla_u_phi_v_j(0.0);
381  jac_u.mv(phi_v[j],nabla_u_phi_v_j);
382  // compute (grad phi_v_j) u (matrix-vector product)
383  Range_V nabla_phi_v_j_u(0.0);
384  jac_phi_v[j].mv(val_u,nabla_phi_v_j_u);
385  mat.accumulate(lfsv_v,i,lfsu_v,j, rho * ((nabla_u_phi_v_j*phi_v[i]) + (nabla_phi_v_j_u*phi_v[i])) * weight);
386  } // end navier
387 
388  } // end j
389 
390  for(size_type j=0; j<lfsv_p.size(); j++) {
391  //================================================//
392  // - p * div v
393  // - q * div u
394  //================================================//
395  mat.accumulate(lfsv_v,i,lfsu_p,j, -phi_p[j] * div_phi_v[i] * weight);
396  mat.accumulate(lfsv_p,j,lfsu_v,i, -phi_p[j] * div_phi_v[i] * incomp_scaling * weight);
397  }
398  } // end i
399 
400  } // end loop quadrature points
401  } // end jacobian_volume
402 
403  // skeleton term, each face is only visited ONCE
404  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
405  void alpha_skeleton (const IG& ig,
406  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
407  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
408  R& r_s, R& r_n) const
409  {
410  // dimensions
411  const unsigned int dim = IG::Geometry::dimension;
412  const unsigned int dimw = IG::Geometry::dimensionworld;
413 
414  // subspaces
415  static_assert
416  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
417 
418  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
419  const LFSV_V& lfsv_v_s = lfsv_s.template child<VBLOCK>();
420  const LFSV_V& lfsu_v_s = lfsu_s.template child<VBLOCK>();
421  const LFSV_V& lfsv_v_n = lfsv_n.template child<VBLOCK>();
422  const LFSV_V& lfsu_v_n = lfsu_n.template child<VBLOCK>();
423 
424  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
425  const LFSV_P& lfsv_p_s = lfsv_s.template child<PBLOCK>();
426  const LFSV_P& lfsu_p_s = lfsu_s.template child<PBLOCK>();
427  const LFSV_P& lfsv_p_n = lfsv_n.template child<PBLOCK>();
428  const LFSV_P& lfsu_p_n = lfsu_n.template child<PBLOCK>();
429 
430  // domain and range field type
431  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
432  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
434  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
435  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
436  typedef typename BasisSwitch_V::DomainField DF;
437  typedef typename BasisSwitch_V::RangeField RF;
438  typedef typename BasisSwitch_V::Range Range_V;
439  typedef typename BasisSwitch_P::Range Range_P;
440  typedef typename LFSV::Traits::SizeType size_type;
441 
442  // make copy of inside and outside cell w.r.t. the intersection
443  auto inside_cell = ig.inside();
444  auto outside_cell = ig.outside();
445 
446  // select quadrature rule
447  Dune::GeometryType gtface = ig.geometry().type();
448  const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
449  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
450  const int qorder = 2*v_order + det_jac_order + superintegration_order;
451  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
452 
453  const int epsilon = prm.epsilonIPSymmetryFactor();
454  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
455 
456  // loop over quadrature points and integrate normal flux
457  for (const auto& ip : rule)
458  {
459 
460  // position of quadrature point in local coordinates of element
461  Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
462  Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
463 
464  const RF penalty_factor = prm.getFaceIP(ig,ip.position());
465 
466  // values of velocity shape functions
467  std::vector<Range_V> phi_v_s(lfsv_v_s.size());
468  std::vector<Range_V> phi_v_n(lfsv_v_n.size());
469  FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
470  FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
471 
472  // values of pressure shape functions
473  std::vector<Range_P> phi_p_s(lfsv_p_s.size());
474  std::vector<Range_P> phi_p_n(lfsv_p_n.size());
475  FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
476  FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
477 
478  // compute pressure value
479  Range_P val_p_s(0.0);
480  Range_P val_p_n(0.0);
481  for (size_type i=0; i<lfsu_p_s.size(); i++)
482  val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
483  for (size_type i=0; i<lfsu_p_n.size(); i++)
484  val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
485 
486  // evaluate jacobian of velocity shape functions on reference element
487  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
488  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
489  VectorBasisSwitch_V::jacobian
490  (FESwitch_V::basis(lfsv_v_s.finiteElement()), inside_cell.geometry(), local_s, jac_phi_v_s);
491  VectorBasisSwitch_V::jacobian
492  (FESwitch_V::basis(lfsv_v_n.finiteElement()), outside_cell.geometry(), local_n, jac_phi_v_n);
493 
494  // compute velocity value, jacobian, and divergence
495  Range_V val_u_s(0.0);
496  Range_V val_u_n(0.0);
497  Dune::FieldMatrix<RF,dim,dim> jac_u_s(0.0);
498  Dune::FieldMatrix<RF,dim,dim> jac_u_n(0.0);
499  for (size_type i=0; i<lfsu_v_s.size(); i++){
500  val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
501  jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
502  }
503  for (size_type i=0; i<lfsu_v_n.size(); i++){
504  val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
505  jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
506  }
507 
508  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
509  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
510  const RF mu = prm.mu(ig,ip.position());
511 
512  const RF factor = mu * weight;
513 
514  // compute jump in velocity
515  const Dune::FieldVector<DF,dimw> jump = val_u_s - val_u_n;
516 
517  // compute mean in pressure
518  const RF mean_p = 0.5*(val_p_s + val_p_n);
519 
520  // compute flux of velocity jacobian
521  Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
522  add_compute_flux(jac_u_s,normal,flux_jac_u);
523  add_compute_flux(jac_u_n,normal,flux_jac_u);
524  flux_jac_u *= 0.5;
525 
526  // loop over test functions, same element
527  for (size_t i=0; i<lfsv_v_s.size(); i++) {
528  //================================================//
529  // diffusion term
530  //================================================//
531  r_s.accumulate(lfsv_v_s, i, -(flux_jac_u * phi_v_s[i]) * factor);
532 
533  //================================================//
534  // (non-)symmetric IP term
535  //================================================//
536  Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
537  add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi);
538  r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
539 
540  //================================================//
541  // standard IP term integral
542  //================================================//
543  r_s.accumulate(lfsv_v_s,i, penalty_factor * (jump*phi_v_s[i]) * weight);
544 
545  //================================================//
546  // pressure-velocity-coupling in momentum equation
547  //================================================//
548  r_s.accumulate(lfsv_v_s,i, mean_p * (phi_v_s[i]*normal) * weight);
549  }
550 
551  // loop over test functions, neighbour element
552  for (size_t i=0; i<lfsv_v_n.size(); i++) {
553  //================================================//
554  // diffusion term
555  //================================================//
556  r_n.accumulate(lfsv_v_n, i, (flux_jac_u * phi_v_n[i]) * factor);
557 
558  //================================================//
559  // (non-)symmetric IP term
560  //================================================//
561  Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
562  add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi);
563  r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
564 
565  //================================================//
566  // standard IP term integral
567  //================================================//
568  r_n.accumulate(lfsv_v_n,i, -penalty_factor * (jump*phi_v_n[i]) * weight);
569 
570  //================================================//
571  // pressure-velocity-coupling in momentum equation
572  //================================================//
573  r_n.accumulate(lfsv_v_n,i, -mean_p * (phi_v_n[i]*normal) * weight);
574  }
575 
576  //================================================//
577  // incompressibility constraint
578  //================================================//
579  for (size_t i=0; i<lfsv_p_s.size(); i++)
580  r_s.accumulate(lfsv_p_s,i, 0.5*phi_p_s[i] * (jump*normal) * incomp_scaling * weight);
581  for (size_t i=0; i<lfsv_p_n.size(); i++)
582  r_n.accumulate(lfsv_p_n,i, 0.5*phi_p_n[i] * (jump*normal) * incomp_scaling * weight);
583 
584  } // end loop quadrature points
585  } // end alpha_skeleton
586 
587  // jacobian of skeleton term, each face is only visited ONCE
588  template<typename IG, typename LFSU, typename X, typename LFSV,
589  typename LocalMatrix>
590  void jacobian_skeleton (const IG& ig,
591  const LFSU& lfsu_s, const X&, const LFSV& lfsv_s,
592  const LFSU& lfsu_n, const X&, const LFSV& lfsv_n,
593  LocalMatrix& mat_ss, LocalMatrix& mat_sn,
594  LocalMatrix& mat_ns, LocalMatrix& mat_nn) const
595  {
596  // dimensions
597  const unsigned int dim = IG::Geometry::dimension;
598  const unsigned int dimw = IG::Geometry::dimensionworld;
599 
600  // subspaces
601  static_assert
602  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
603 
604  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
605  const LFSV_V& lfsv_v_s = lfsv_s.template child<VBLOCK>();
606  const LFSV_V& lfsu_v_s = lfsu_s.template child<VBLOCK>();
607  const LFSV_V& lfsv_v_n = lfsv_n.template child<VBLOCK>();
608  const LFSV_V& lfsu_v_n = lfsu_n.template child<VBLOCK>();
609 
610  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
611  const LFSV_P& lfsv_p_s = lfsv_s.template child<PBLOCK>();
612  const LFSV_P& lfsu_p_s = lfsu_s.template child<PBLOCK>();
613  const LFSV_P& lfsv_p_n = lfsv_n.template child<PBLOCK>();
614  const LFSV_P& lfsu_p_n = lfsu_n.template child<PBLOCK>();
615 
616  // domain and range field type
617  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
618  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
620  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
621  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
622  typedef typename BasisSwitch_V::DomainField DF;
623  typedef typename BasisSwitch_V::RangeField RF;
624  typedef typename BasisSwitch_V::Range Range_V;
625  typedef typename BasisSwitch_P::Range Range_P;
626  typedef typename LFSV::Traits::SizeType size_type;
627 
628  // make copy of inside and outside cell w.r.t. the intersection
629  auto inside_cell = ig.inside();
630  auto outside_cell = ig.outside();
631 
632  // select quadrature rule
633  Dune::GeometryType gtface = ig.geometry().type();
634  const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
635  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
636  const int qorder = 2*v_order + det_jac_order + superintegration_order;
637  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
638 
639  const int epsilon = prm.epsilonIPSymmetryFactor();
640  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
641 
642  // loop over quadrature points and integrate normal flux
643  for (const auto& ip : rule)
644  {
645 
646  // position of quadrature point in local coordinates of element
647  Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
648  Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
649 
650  const RF penalty_factor = prm.getFaceIP(ig,ip.position());
651 
652  // values of velocity shape functions
653  std::vector<Range_V> phi_v_s(lfsv_v_s.size());
654  std::vector<Range_V> phi_v_n(lfsv_v_n.size());
655  FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
656  FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
657 
658  // values of pressure shape functions
659  std::vector<Range_P> phi_p_s(lfsv_p_s.size());
660  std::vector<Range_P> phi_p_n(lfsv_p_n.size());
661  FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
662  FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
663 
664  // evaluate jacobian of velocity shape functions on reference element
665  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
666  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
667  VectorBasisSwitch_V::jacobian
668  (FESwitch_V::basis(lfsv_v_s.finiteElement()), inside_cell.geometry(), local_s, jac_phi_v_s);
669  VectorBasisSwitch_V::jacobian
670  (FESwitch_V::basis(lfsv_v_n.finiteElement()), outside_cell.geometry(), local_n, jac_phi_v_n);
671 
672  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
673  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
674  const RF mu = prm.mu(ig,ip.position());
675 
676  const RF factor = mu * weight;
677 
678  //============================================
679  // loop over test functions, same element
680  //============================================
681  for(size_type i=0; i<lfsv_v_s.size(); i++) {
682 
683  // compute flux
684  Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
685  add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi_i);
686 
687  //============================================
688  // diffusion
689  // (non-)symmetric IP-Term
690  // standard IP integral
691  //============================================
692  for(size_type j=0; j<lfsu_v_s.size(); j++) {
693  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
694  add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
695 
696  mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
697  mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
698  mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, penalty_factor * (phi_v_s[j]*phi_v_s[i]) * weight);
699  }
700 
701  for(size_type j=0; j<lfsu_v_n.size(); j++) {
702  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
703  add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
704 
705  mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
706  mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
707  mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -penalty_factor * (phi_v_n[j]*phi_v_s[i]) * weight);
708  }
709 
710  //============================================
711  // pressure-velocity coupling in momentum equation
712  //============================================
713  for(size_type j=0; j<lfsu_p_s.size(); j++) {
714  mat_ss.accumulate(lfsv_v_s,i,lfsu_p_s,j, 0.5*phi_p_s[j] * (phi_v_s[i]*normal) * weight);
715  }
716 
717  for(size_type j=0; j<lfsu_p_n.size(); j++) {
718  mat_sn.accumulate(lfsv_v_s,i,lfsu_p_n,j, 0.5*phi_p_n[j] * (phi_v_s[i]*normal) * weight);
719  }
720  } // end i (same)
721 
722  //============================================
723  // loop over test functions, neighbour element
724  //============================================
725  for(size_type i=0; i<lfsv_v_n.size(); i++) {
726 
727  // compute flux
728  Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
729  add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi_i);
730 
731  //============================================
732  // diffusion
733  // (non-)symmetric IP-Term
734  // standard IP integral
735  //============================================
736  for(size_type j=0; j<lfsu_v_s.size(); j++) {
737  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
738  add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
739 
740  mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
741  mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
742  mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, -penalty_factor * (phi_v_s[j]*phi_v_n[i]) * weight);
743  }
744 
745  for(size_type j=0; j<lfsu_v_n.size(); j++) {
746  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
747  add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
748 
749  mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
750  mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
751  mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, penalty_factor * (phi_v_n[j]*phi_v_n[i]) * weight);
752  }
753 
754  //============================================
755  // pressure-velocity coupling in momentum equation
756  //============================================
757  for(size_type j=0; j<lfsu_p_s.size(); j++) {
758  mat_ns.accumulate(lfsv_v_n,i,lfsu_p_s,j, -0.5*phi_p_s[j] * (phi_v_n[i]*normal) * weight);
759  }
760 
761  for(size_type j=0; j<lfsu_p_n.size(); j++) {
762  mat_nn.accumulate(lfsv_v_n,i,lfsu_p_n,j, -0.5*phi_p_n[j] * (phi_v_n[i]*normal) * weight);
763  }
764  } // end i (neighbour)
765 
766  //================================================//
767  // \int <q> [u] n
768  //================================================//
769  for(size_type i=0; i<lfsv_p_s.size(); i++) {
770  for(size_type j=0; j<lfsu_v_s.size(); j++)
771  mat_ss.accumulate(lfsv_p_s,i,lfsu_v_s,j, 0.5*phi_p_s[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
772 
773  for(size_type j=0; j<lfsu_v_n.size(); j++)
774  mat_sn.accumulate(lfsv_p_s,i,lfsu_v_n,j, -0.5*phi_p_s[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
775  }
776 
777  for(size_type i=0; i<lfsv_p_n.size(); i++) {
778  for(size_type j=0; j<lfsu_v_s.size(); j++)
779  mat_ns.accumulate(lfsv_p_n,i,lfsu_v_s,j, 0.5*phi_p_n[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
780 
781  for(size_type j=0; j<lfsu_v_n.size(); j++)
782  mat_nn.accumulate(lfsv_p_n,i,lfsu_v_n,j, -0.5*phi_p_n[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
783  }
784 
785  } // end loop quadrature points
786  } // end jacobian_skeleton
787 
788  // boundary term
789  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
790  void alpha_boundary (const IG& ig,
791  const LFSU& lfsu, const X& x, const LFSV& lfsv,
792  R& r) const
793  {
794  // dimensions
795  const unsigned int dim = IG::Geometry::dimension;
796  const unsigned int dimw = IG::Geometry::dimensionworld;
797 
798  // subspaces
799  static_assert
800  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
801 
802  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
803  const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
804  const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
805 
806  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
807  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
808  const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
809 
810  // domain and range field type
811  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
812  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
814  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
815  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
816  typedef typename BasisSwitch_V::DomainField DF;
817  typedef typename BasisSwitch_V::RangeField RF;
818  typedef typename BasisSwitch_V::Range Range_V;
819  typedef typename BasisSwitch_P::Range Range_P;
820  typedef typename LFSV::Traits::SizeType size_type;
821 
822  // make copy of inside cell w.r.t. the boundary
823  auto inside_cell = ig.inside();
824 
825  // select quadrature rule
826  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
827  Dune::GeometryType gtface = ig.geometry().type();
828  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
829  const int qorder = 2*v_order + det_jac_order + superintegration_order;
830  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
831 
832  const int epsilon = prm.epsilonIPSymmetryFactor();
833  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
834 
835  // loop over quadrature points and integrate normal flux
836  for (const auto& ip : rule)
837  {
838  // position of quadrature point in local coordinates of element
839  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
840 
841  const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
842 
843  // values of velocity shape functions
844  std::vector<Range_V> phi_v(lfsv_v.size());
845  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
846 
847  // values of pressure shape functions
848  std::vector<Range_P> phi_p(lfsv_p.size());
849  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
850 
851  // evaluate jacobian of basis functions on reference element
852  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
853  VectorBasisSwitch_V::jacobian
854  (FESwitch_V::basis(lfsv_v.finiteElement()), inside_cell.geometry(), local, jac_phi_v);
855 
856  // compute pressure value
857  Range_P val_p(0.0);
858  for (size_type i=0; i<lfsu_p.size(); i++)
859  val_p.axpy(x(lfsu_p,i),phi_p[i]);
860 
861  // compute u and velocity jacobian
862  Range_V val_u(0.0);
863  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
864  for (size_type i=0; i<lfsu_v.size(); i++){
865  val_u.axpy(x(lfsu_v,i),phi_v[i]);
866  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
867  }
868 
869  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
870  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
871  const RF mu = prm.mu(ig,ip.position());
872 
873  // evaluate boundary condition type
874  typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
875 
876  if (bctype == BC::VelocityDirichlet) {
877  // compute jump relative to Dirichlet value
878  typename PRM::Traits::VelocityRange u0(prm.g(inside_cell,local));
879  const Dune::FieldVector<DF,dimw> jump = val_u - u0;
880 
881  // compute flux of velocity jacobian
882  Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
883  add_compute_flux(jac_u,normal,flux_jac_u);
884 
885  for (size_t i=0; i<lfsv_v.size(); i++) {
886  //================================================//
887  // diffusion term
888  //================================================//
889  r.accumulate(lfsv_v,i, -mu * (flux_jac_u * phi_v[i]) * weight);
890 
891  //================================================//
892  // (non-)symmetric IP term
893  //================================================//
894  Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
895  add_compute_flux(jac_phi_v[i],normal,flux_jac_phi);
896  r.accumulate(lfsv_v,i, mu * epsilon * (flux_jac_phi*jump) * weight);
897 
898  //================================================//
899  // standard IP term integral
900  //================================================//
901  r.accumulate(lfsv_v,i, (jump*phi_v[i]) * penalty_factor * weight);
902 
903  //================================================//
904  // pressure-velocity-coupling in momentum equation
905  //================================================//
906  r.accumulate(lfsv_v,i, val_p * (phi_v[i]*normal) * weight);
907  } // end i
908 
909  //================================================//
910  // incompressibility constraint
911  //================================================//
912  for(size_type i=0; i<lfsv_p.size(); i++) {
913  r.accumulate(lfsv_p,i, phi_p[i] * (jump*normal) * incomp_scaling * weight);
914  }
915  } // Velocity Dirichlet
916 
917  if (bctype == BC::StressNeumann) {
918  typename PRM::Traits::VelocityRange stress(prm.j(ig,ip.position(),normal));
919 
920  for(size_type i=0; i<lfsv_v.size(); i++) {
921  r.accumulate(lfsv_v,i, (stress*phi_v[i]) * weight);
922  }
923  } // Pressure Dirichlet
924 
925  } // end loop quadrature points
926  } // end alpha_boundary
927 
928  // jacobian of boundary term
929  template<typename IG, typename LFSU, typename X, typename LFSV,
930  typename LocalMatrix>
931  void jacobian_boundary (const IG& ig,
932  const LFSU& lfsu, const X& x, const LFSV& lfsv,
933  LocalMatrix& mat) const
934  {
935  // dimensions
936  const unsigned int dim = IG::Geometry::dimension;
937  const unsigned int dimw = IG::Geometry::dimensionworld;
938 
939  // subspaces
940  static_assert
941  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
942 
943  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
944  const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
945  const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
946 
947  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
948  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
949  const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
950 
951  // domain and range field type
952  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
953  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
955  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
956  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
957  typedef typename BasisSwitch_V::DomainField DF;
958  typedef typename BasisSwitch_V::RangeField RF;
959  typedef typename BasisSwitch_V::Range Range_V;
960  typedef typename BasisSwitch_P::Range Range_P;
961  typedef typename LFSV::Traits::SizeType size_type;
962 
963  // make copy of inside cell w.r.t. the boundary
964  auto inside_cell = ig.inside();
965 
966  // select quadrature rule
967  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
968  Dune::GeometryType gtface = ig.geometry().type();
969  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
970  const int qorder = 2*v_order + det_jac_order + superintegration_order;
971  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
972 
973  const int epsilon = prm.epsilonIPSymmetryFactor();
974  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
975 
976  // loop over quadrature points and integrate normal flux
977  for (const auto& ip : rule)
978  {
979  // position of quadrature point in local coordinates of element
980  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
981 
982  const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
983 
984  // values of velocity shape functions
985  std::vector<Range_V> phi_v(lfsv_v.size());
986  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
987 
988  // values of pressure shape functions
989  std::vector<Range_P> phi_p(lfsv_p.size());
990  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
991 
992  // evaluate jacobian of basis functions on reference element
993  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
994  VectorBasisSwitch_V::jacobian
995  (FESwitch_V::basis(lfsv_v.finiteElement()), inside_cell.geometry(), local, jac_phi_v);
996 
997  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
998  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
999  const RF mu = prm.mu(ig,ip.position());
1000 
1001  // evaluate boundary condition type
1002  typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
1003 
1004  if (bctype == BC::VelocityDirichlet) {
1005 
1006  for(size_type i=0; i<lfsv_v.size(); i++) {
1007  // compute flux
1008  Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
1009  add_compute_flux(jac_phi_v[i],normal,flux_jac_phi_i);
1010 
1011  for(size_type j=0; j<lfsu_v.size(); j++) {
1012  //================================================//
1013  // diffusion term
1014  // (non-)symmetric IP term
1015  //================================================//
1016  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
1017  add_compute_flux(jac_phi_v[j],normal,flux_jac_phi_j);
1018 
1019  mat.accumulate(lfsv_v,i,lfsu_v,j, -mu * (flux_jac_phi_j*phi_v[i]) * weight);
1020  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * epsilon * (flux_jac_phi_i*phi_v[j]) *weight);
1021 
1022  //================================================//
1023  // standard IP term integral
1024  //================================================//
1025  mat.accumulate(lfsv_v,i,lfsu_v,j, (phi_v[j]*phi_v[i]) * penalty_factor * weight);
1026  }
1027 
1028  //================================================//
1029  // pressure-velocity-coupling in momentum equation
1030  //================================================//
1031  for(size_type j=0; j<lfsu_p.size(); j++) {
1032  mat.accumulate(lfsv_v,i,lfsu_p,j, phi_p[j] * (phi_v[i]*normal) * weight);
1033  }
1034  } // end i
1035 
1036  //================================================//
1037  // incompressibility constraint
1038  //================================================//
1039  for(size_type i=0; i<lfsv_p.size(); i++) {
1040  for(size_type j=0; j<lfsu_v.size(); j++) {
1041  mat.accumulate(lfsv_p,i,lfsu_v,j, phi_p[i] * (phi_v[j]*normal) * incomp_scaling * weight);
1042  }
1043  }
1044 
1045  } // Velocity Dirichlet
1046 
1047  } // end loop quadrature points
1048  } // end jacobian_boundary
1049 
1050  // volume integral depending only on test functions,
1051  // contains f on the right hand side
1052  template<typename EG, typename LFSV, typename R>
1053  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
1054  {
1055  const unsigned int dim = EG::Geometry::dimension;
1056 
1057  // subspaces
1058  static_assert
1059  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokesVelVecFEM");
1060 
1061  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
1062  const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
1063 
1064  // domain and range field type
1065  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1066  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1067  typedef VectorBasisInterfaceSwitch<typename FESwitch_V::Basis > VectorBasisSwitch_V;
1068  typedef typename BasisSwitch_V::DomainField DF;
1069  typedef typename BasisSwitch_V::RangeField RF;
1070  typedef typename BasisSwitch_V::Range Range_V;
1071  typedef typename LFSV::Traits::SizeType size_type;
1072 
1073  // select quadrature rule
1074  Dune::GeometryType gt = eg.geometry().type();
1075  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1076  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1077  const int qorder = 2*v_order + det_jac_order + superintegration_order;
1078 
1079  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1080 
1081  // loop over quadrature points
1082  for (const auto& ip : rule)
1083  {
1084  const Dune::FieldVector<DF,dim> local = ip.position();
1085  //const Dune::FieldVector<DF,dimw> global = eg.geometry().global(local);
1086 
1087  // values of velocity shape functions
1088  std::vector<Range_V> phi_v(lfsv_v.size());
1089  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1090 
1091  const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1092 
1093  // evaluate source term
1094  typename PRM::Traits::VelocityRange fval(prm.f(eg,local));
1095 
1096  //================================================//
1097  // \int (f*v)
1098  //================================================//
1099  for(size_type i=0; i<lfsv_v.size(); i++)
1100  r.accumulate(lfsv_v,i, -(fval*phi_v[i]) * weight);
1101 
1102  } // end loop quadrature points
1103  } // end lambda_volume
1104 
1105  private :
1106 
1107  template<class M, class RF>
1108  static void contraction(const M & a, const M & b, RF & v)
1109  {
1110  v = 0;
1111  const int an = a.N();
1112  const int am = a.M();
1113  for(int r=0; r<an; ++r)
1114  for(int c=0; c<am; ++c)
1115  v += a[r][c] * b[r][c];
1116  }
1117 
1118  template<class DU, class R>
1119  static void add_compute_flux(const DU & du, const R & n, R & result)
1120  {
1121  const int N = du.N();
1122  const int M = du.M();
1123  for(int r=0; r<N; ++r)
1124  for(int c=0; c<M; ++c)
1125  result[r] += du[r][c] * n[c];
1126  }
1127 
1128  PRM& prm;
1129  const int superintegration_order;
1130  Real current_dt;
1131  }; // end class DGNavierStokesVelVecFEM
1132 
1133  } // end namespace PDELab
1134 } // end namespace Dune
1135 #endif
Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:65
Definition: dgnavierstokesvelvecfem.hh:31
void setTime(doublet_)
set time for subsequent evaluation
Definition: idefault.hh:104
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:162
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokesvelvecfem.hh:284
static const std::size_t dimRange
export dimension of the values
Definition: dgnavierstokesvelvecfem.hh:37
A local operator for solving the stokes equation using a DG discretization and a vector-valued finite...
Definition: dgnavierstokesvelvecfem.hh:101
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:1053
const IG & ig
Definition: constraints.hh:147
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Definition: dgnavierstokesvelvecfem.hh:122
static void jacobian(const Basis &basis, const Geometry &geometry, const DomainLocal &xl, std::vector< FieldMatrix< RangeField, dimRange, Geometry::coorddimension > > &jac)
Compute global jacobian matrix for vector valued bases.
Definition: dgnavierstokesvelvecfem.hh:41
Basis::Traits::RangeField RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:35
Default flags for all local operators.
Definition: flags.hh:18
DGNavierStokesVelVecFEM(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokesvelvecfem.hh:142
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
void preStep(RealType, RealType dt, int)
Definition: dgnavierstokesvelvecfem.hh:148
void setTime(Real t)
Definition: dgnavierstokesvelvecfem.hh:154
Definition: dgnavierstokesvelvecfem.hh:128
Definition: dgnavierstokesvelvecfem.hh:127
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: dgnavierstokesvelvecfem.hh:405
static void jacobian(const Basis &basis, const Geometry &geometry, const DomainLocal &xl, std::vector< FieldMatrix< RangeField, dimRange, Geometry::coorddimension > > &jac)
Compute global jacobian matrix for vector valued bases.
Definition: dgnavierstokesvelvecfem.hh:73
Definition: adaptivity.hh:27
Definition: dgnavierstokesvelvecfem.hh:125
static const int dim
Definition: adaptivity.hh:83
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:790
sparsity pattern generator
Definition: pattern.hh:29
Definition: dgnavierstokesvelvecfem.hh:119
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &, const LFSV &lfsv_n, LocalMatrix &mat_ss, LocalMatrix &mat_sn, LocalMatrix &mat_ns, LocalMatrix &mat_nn) const
Definition: dgnavierstokesvelvecfem.hh:590
Basis::Traits::DomainLocal DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:33
Definition: dgnavierstokesvelvecfem.hh:118
const EG & eg
Definition: constraints.hh:280
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokesvelvecfem.hh:931
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: dgnavierstokesvelvecfem.hh:126