dune-pdelab  2.4-dev
dgnavierstokes.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_DGNAVIERSTOKES_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKES_HH
6 
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/parametertreeparser.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 
35  template<typename PRM>
40  {
42  typedef typename PRM::Traits::RangeField RF;
43 
45  typedef typename InstatBase::RealType Real;
46 
47  static const bool navier = PRM::assemble_navier;
48  static const bool full_tensor = PRM::assemble_full_tensor;
49 
50  public:
51 
52  // pattern assembly flags
53  enum { doPatternVolume = true };
54  enum { doPatternSkeleton = true };
55 
56  // call the assembler for each face only once
57  enum { doSkeletonTwoSided = false };
58 
59  // residual assembly flags
60  enum { doAlphaVolume = true };
61  enum { doAlphaSkeleton = true };
62  enum { doAlphaBoundary = true };
63  enum { doLambdaVolume = true };
64  enum { doLambdaBoundary = true };
65 
78  DGNavierStokes (PRM& _prm, int _superintegration_order=0) :
79  prm(_prm), superintegration_order(_superintegration_order),
80  current_dt(1.0)
81  {}
82 
83  // Store current dt
84  void preStep (RealType , RealType dt, int )
85  {
86  current_dt = dt;
87  }
88 
89  // set time in parameter class
90  void setTime(Real t)
91  {
93  prm.setTime(t);
94  }
95 
96  // volume integral depending on test and ansatz functions
97  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
98  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
99  {
100  // dimensions
101  const unsigned int dim = EG::Geometry::dimension;
102 
103  // subspaces
104  static_assert
105  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
106  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
107  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
108  static_assert
109  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
110 
111  // ... we assume all velocity components are the same
112  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
113  const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
114  const unsigned int vsize = lfsv_v.size();
115  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
116  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
117  const unsigned int psize = lfsv_p.size();
118 
119  // domain and range field type
120  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
121  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
122  typedef typename BasisSwitch_V::DomainField DF;
123  typedef typename BasisSwitch_V::Range RT;
124  typedef typename BasisSwitch_V::RangeField RF;
125  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
126  typedef typename LFSV::Traits::SizeType size_type;
127 
128  // select quadrature rule
129  Dune::GeometryType gt = eg.geometry().type();
130  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
131  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
132  const int jac_order = gt.isSimplex() ? 0 : 1;
133  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
134  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
135 
136  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
137 
138  // loop over quadrature points
139  for (const auto& ip : rule)
140  {
141  const Dune::FieldVector<DF,dim> local = ip.position();
142  const RF mu = prm.mu(eg,local);
143  const RF rho = prm.rho(eg,local);
144 
145  // compute u (if Navier term enabled)
146  std::vector<RT> phi_v(vsize);
147  Dune::FieldVector<RF,dim> vu(0.0);
148  if(navier) {
149  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
150 
151  for(unsigned int d=0; d<dim; ++d) {
152  const LFSV_V& lfsu_v = lfsv_pfs_v.child(d);
153  for(size_type i=0; i<vsize; i++)
154  vu[d] += x(lfsu_v,i) * phi_v[i];
155  }
156  } // end navier
157 
158  // and value of pressure shape functions
159  std::vector<RT> phi_p(psize);
160  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
161 
162  // compute pressure
163  RF p(0.0);
164  for(size_type i=0; i<psize; i++)
165  p += x(lfsv_p,i) * phi_p[i];
166 
167  // compute gradients
168  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
169  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
170  eg.geometry(), local, grad_phi_v);
171 
172  // compute velocity jacobian
173  Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
174  for(unsigned int d = 0; d<dim; ++d) {
175  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
176  for(size_type i=0; i<vsize; i++)
177  jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
178  }
179 
180  const RF detj = eg.geometry().integrationElement(ip.position());
181  const RF weight = ip.weight() * detj;
182 
183  for(unsigned int d = 0; d<dim; ++d) {
184  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
185 
186  for(size_type i=0; i<vsize; i++) {
187  //================================================//
188  // \int (mu*grad_u*grad_v)
189  //================================================//
190  r.accumulate(lfsv_v,i, mu * (jacu[d]*grad_phi_v[i][0]) * weight);
191 
192  // Assemble symmetric part for (grad u)^T
193  if(full_tensor)
194  for(unsigned int dd = 0; dd<dim; ++dd)
195  r.accumulate(lfsv_v,i, mu * jacu[dd][d] * grad_phi_v[i][0][dd] * weight);
196 
197  //================================================//
198  // \int -p \nabla\cdot v
199  //================================================//
200  r.accumulate(lfsv_v,i, -p * grad_phi_v[i][0][d] * weight);
201 
202  //================================================//
203  // \int \rho ((u\cdot\nabla ) u )\cdot v
204  //================================================//
205  if(navier) {
206  // compute u * grad u_d
207  const RF u_nabla_u = vu * jacu[d];
208 
209  r.accumulate(lfsv_v,i, rho * u_nabla_u * phi_v[i] * weight);
210  } // end navier
211 
212  } // end i
213  } // end d
214 
215  //================================================//
216  // \int -q \nabla\cdot u
217  //================================================//
218  for(size_type i=0; i<psize; i++)
219  for(unsigned int d = 0; d < dim; ++d)
220  // divergence of u is the trace of the velocity jacobian
221  r.accumulate(lfsv_p,i, -jacu[d][d] * phi_p[i] * incomp_scaling * weight);
222 
223  } // end loop quadrature points
224  } // end alpha_volume
225 
226  // jacobian of volume term
227  template<typename EG, typename LFSU, typename X, typename LFSV,
228  typename LocalMatrix>
229  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
230  LocalMatrix& mat) const
231  {
232  // dimensions
233  const unsigned int dim = EG::Geometry::dimension;
234 
235  // subspaces
236  static_assert
237  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
238  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
239  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
240  static_assert
241  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
242 
243  // ... we assume all velocity components are the same
244  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
245  const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
246  const unsigned int vsize = lfsv_v.size();
247  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
248  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
249  const unsigned int psize = lfsv_p.size();
250 
251  // domain and range field type
252  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
253  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
254  typedef typename BasisSwitch_V::DomainField DF;
255  typedef typename BasisSwitch_V::Range RT;
256  typedef typename BasisSwitch_V::RangeField RF;
257  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
258  typedef typename LFSV::Traits::SizeType size_type;
259 
260  // select quadrature rule
261  Dune::GeometryType gt = eg.geometry().type();
262  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
263  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
264  const int jac_order = gt.isSimplex() ? 0 : 1;
265  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
266  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
267 
268  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
269 
270  // loop over quadrature points
271  for (const auto& ip : rule)
272  {
273  const Dune::FieldVector<DF,dim> local = ip.position();
274  const RF mu = prm.mu(eg,local);
275  const RF rho = prm.rho(eg,local);
276 
277  // compute u (if Navier term enabled)
278  std::vector<RT> phi_v(vsize);
279  Dune::FieldVector<RF,dim> vu(0.0);
280  if(navier) {
281  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
282 
283  for(unsigned int d=0; d<dim; ++d) {
284  const LFSV_V& lfsu_v = lfsv_pfs_v.child(d);
285  for(size_type i=0; i<vsize; i++)
286  vu[d] += x(lfsu_v,i) * phi_v[i];
287  }
288  } // end navier
289 
290  // and value of pressure shape functions
291  std::vector<RT> phi_p(psize);
292  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
293 
294  // compute gradients
295  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
296  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
297  eg.geometry(), local, grad_phi_v);
298 
299  const RF detj = eg.geometry().integrationElement(ip.position());
300  const RF weight = ip.weight() * detj;
301 
302  for(unsigned int dv = 0; dv<dim; ++dv) {
303  const LFSV_V& lfsv_v = lfsv_pfs_v.child(dv);
304 
305  // gradient of dv-th velocity component
306  Dune::FieldVector<RF,dim> gradu_dv(0.0);
307  if(navier)
308  for(size_type l=0; l<vsize; ++l)
309  gradu_dv.axpy(x(lfsv_v,l), grad_phi_v[l][0]);
310 
311  for(size_type i=0; i<vsize; i++) {
312 
313  for(size_type j=0; j<vsize; j++) {
314  //================================================//
315  // \int (mu*grad_u*grad_v)
316  //================================================//
317  mat.accumulate(lfsv_v,i,lfsv_v,j, mu * (grad_phi_v[j][0]*grad_phi_v[i][0]) * weight);
318 
319  // Assemble symmetric part for (grad u)^T
320  if(full_tensor)
321  for(unsigned int du = 0; du<dim; ++du) {
322  const LFSV_V& lfsu_v = lfsv_pfs_v.child(du);
323  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (grad_phi_v[j][0][dv]*grad_phi_v[i][0][du]) * weight);
324  }
325  }
326 
327  //================================================//
328  // - q * div u
329  // - p * div v
330  //================================================//
331  for(size_type j=0; j<psize; j++) {
332  mat.accumulate(lfsv_p,j,lfsv_v,i, -phi_p[j] * grad_phi_v[i][0][dv] * incomp_scaling * weight);
333  mat.accumulate(lfsv_v,i,lfsv_p,j, -phi_p[j] * grad_phi_v[i][0][dv] * weight);
334  }
335 
336  //================================================//
337  // \int \rho ((u\cdot\nabla ) u )\cdot v
338  //================================================//
339  if(navier) {
340 
341  // block diagonal contribution
342  for(size_type j=0; j<vsize; j++)
343  mat.accumulate(lfsv_v,i,lfsv_v,j, rho * (vu * grad_phi_v[j][0]) * phi_v[i] * weight);
344 
345  // remaining contribution
346  for(unsigned int du = 0; du < dim; ++du) {
347  const LFSV_V& lfsu_v = lfsv_pfs_v.child(du);
348  for(size_type j=0; j<vsize; j++)
349  mat.accumulate(lfsv_v,i,lfsu_v,j, rho * phi_v[j] * gradu_dv[du] * phi_v[i] * weight);
350  }
351 
352  } // end navier
353 
354  } // end i
355  } // end dv
356 
357  } // end loop quadrature points
358  } // end jacobian_volume
359 
360  // skeleton integral depending on test and ansatz functions
361  // each face is only visited ONCE!
362  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
363  void alpha_skeleton (const IG& ig,
364  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
365  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
366  R& r_s, R& r_n) const
367  {
368  // dimensions
369  const unsigned int dim = IG::Geometry::dimension;
370  const unsigned int dimw = IG::Geometry::dimensionworld;
371 
372  // subspaces
373  static_assert
374  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
375 
376  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
377  const LFSV_PFS_V& lfsv_s_pfs_v = lfsv_s.template child<VBLOCK>();
378  const LFSV_PFS_V& lfsv_n_pfs_v = lfsv_n.template child<VBLOCK>();
379  static_assert
380  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
381 
382  // ... we assume all velocity components are the same
383  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
384  const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.template child<0>();
385  const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.template child<0>();
386  const unsigned int vsize_s = lfsv_s_v.size();
387  const unsigned int vsize_n = lfsv_n_v.size();
388  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
389  const LFSV_P& lfsv_s_p = lfsv_s.template child<PBLOCK>();
390  const LFSV_P& lfsv_n_p = lfsv_n.template child<PBLOCK>();
391  const unsigned int psize_s = lfsv_s_p.size();
392  const unsigned int psize_n = lfsv_n_p.size();
393 
394  // domain and range field type
395  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
396  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
397  typedef typename BasisSwitch_V::DomainField DF;
398  typedef typename BasisSwitch_V::Range RT;
399  typedef typename BasisSwitch_V::RangeField RF;
400  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
401 
402  // make copy of inside and outside cell w.r.t. the intersection
403  auto inside_cell = ig.inside();
404  auto outside_cell = ig.outside();
405 
406  // select quadrature rule
407  Dune::GeometryType gtface = ig.geometry().type();
408  const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
409  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
410  const int qorder = 2*v_order + det_jac_order + superintegration_order;
411  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
412 
413  const int epsilon = prm.epsilonIPSymmetryFactor();
414  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
415 
416  // loop over quadrature points and integrate normal flux
417  for (const auto& ip : rule)
418  {
419 
420  // position of quadrature point in local coordinates of element
421  Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
422  Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
423 
424  const RF penalty_factor = prm.getFaceIP(ig,ip.position());
425 
426  // value of velocity shape functions
427  std::vector<RT> phi_v_s(vsize_s);
428  std::vector<RT> phi_v_n(vsize_n);
429  FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
430  FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
431 
432  // evaluate u
433  assert(vsize_s == vsize_n);
434  Dune::FieldVector<RF,dim> u_s(0.0), u_n(0.0);
435  for(unsigned int d=0; d<dim; ++d) {
436  const LFSV_V & lfsv_s_v = lfsv_s_pfs_v.child(d);
437  const LFSV_V & lfsv_n_v = lfsv_n_pfs_v.child(d);
438  for(unsigned int i=0; i<vsize_s; i++) {
439  u_s[d] += x_s(lfsv_s_v,i) * phi_v_s[i];
440  u_n[d] += x_n(lfsv_n_v,i) * phi_v_n[i];
441  }
442  }
443 
444  // value of pressure shape functions
445  std::vector<RT> phi_p_s(psize_s);
446  std::vector<RT> phi_p_n(psize_n);
447  FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
448  FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
449 
450  // evaluate pressure
451  assert(psize_s == psize_n);
452  RF p_s(0.0), p_n(0.0);
453  for(unsigned int i=0; i<psize_s; i++) {
454  p_s += x_s(lfsv_s_p,i) * phi_p_s[i];
455  p_n += x_n(lfsv_n_p,i) * phi_p_n[i];
456  }
457 
458  // compute gradients
459  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
460  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
461  inside_cell.geometry(), local_s, grad_phi_v_s);
462 
463  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
464  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
465  outside_cell.geometry(), local_n, grad_phi_v_n);
466 
467  // evaluate velocity jacobian
468  Dune::FieldMatrix<RF,dim,dim> jacu_s(0.0), jacu_n(0.0);
469  for(unsigned int d=0; d<dim; ++d) {
470  const LFSV_V & lfsv_s_v = lfsv_s_pfs_v.child(d);
471  const LFSV_V & lfsv_n_v = lfsv_n_pfs_v.child(d);
472  for(unsigned int i=0; i<vsize_s; i++) {
473  jacu_s[d].axpy(x_s(lfsv_s_v,i), grad_phi_v_s[i][0]);
474  jacu_n[d].axpy(x_n(lfsv_n_v,i), grad_phi_v_n[i][0]);
475  }
476  }
477 
478  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
479  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
480  const RF mu = prm.mu(ig,ip.position());
481 
482  const RF factor = mu * weight;
483 
484  for(unsigned int d=0; d<dim; ++d) {
485  const LFSV_V & lfsv_s_v = lfsv_s_pfs_v.child(d);
486  const LFSV_V & lfsv_n_v = lfsv_n_pfs_v.child(d);
487 
488  //================================================//
489  // diffusion term
490  //================================================//
491  RF val = 0.5 * ((jacu_s[d] * normal) + (jacu_n[d] * normal)) * factor;
492  for(unsigned int i=0; i<vsize_s; i++) {
493  r_s.accumulate(lfsv_s_v,i, -val * phi_v_s[i]);
494  r_n.accumulate(lfsv_n_v,i, val * phi_v_n[i]);
495 
496  if(full_tensor) {
497  for(unsigned int dd=0; dd<dim; ++dd) {
498  RF Tval = 0.5 * (jacu_s[dd][d] + jacu_n[dd][d]) * normal[dd] * factor;
499  r_s.accumulate(lfsv_s_v,i, -Tval * phi_v_s[i]);
500  r_n.accumulate(lfsv_n_v,i, Tval * phi_v_n[i]);
501  }
502  }
503  } // end i
504 
505  //================================================//
506  // (non-)symmetric IP term
507  //================================================//
508  RF jumpu_d = u_s[d] - u_n[d];
509  for(unsigned int i=0; i<vsize_s; i++) {
510  r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * (grad_phi_v_s[i][0] * normal) * jumpu_d * factor);
511  r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * (grad_phi_v_n[i][0] * normal) * jumpu_d * factor);
512 
513  if(full_tensor) {
514  for(unsigned int dd=0; dd<dim; ++dd) {
515  r_s.accumulate(lfsv_s_v,i, epsilon * 0.5 * grad_phi_v_s[i][0][dd] * (u_s[dd] - u_n[dd]) * normal[d] * factor);
516  r_n.accumulate(lfsv_n_v,i, epsilon * 0.5 * grad_phi_v_n[i][0][dd] * (u_s[dd] - u_n[dd]) * normal[d] * factor);
517  }
518  }
519  } // end i
520 
521  //================================================//
522  // standard IP term integral
523  //================================================//
524  for(unsigned int i=0; i<vsize_s; i++) {
525  r_s.accumulate(lfsv_s_v,i, penalty_factor * jumpu_d * phi_v_s[i] * weight);
526  r_n.accumulate(lfsv_n_v,i, -penalty_factor * jumpu_d * phi_v_n[i] * weight);
527  } // end i
528 
529  //================================================//
530  // pressure-velocity-coupling in momentum equation
531  //================================================//
532  RF mean_p = 0.5*(p_s + p_n);
533  for(unsigned int i=0; i<vsize_s; i++) {
534  r_s.accumulate(lfsv_s_v,i, mean_p * phi_v_s[i] * normal[d] * weight);
535  r_n.accumulate(lfsv_n_v,i, -mean_p * phi_v_n[i] * normal[d] * weight);
536  } // end i
537  } // end d
538 
539  //================================================//
540  // incompressibility constraint
541  //================================================//
542  RF jumpu_n = (u_s*normal) - (u_n*normal);
543  for(unsigned int i=0; i<psize_s; i++) {
544  r_s.accumulate(lfsv_s_p,i, 0.5 * phi_p_s[i] * jumpu_n * incomp_scaling * weight);
545  r_n.accumulate(lfsv_n_p,i, 0.5 * phi_p_n[i] * jumpu_n * incomp_scaling * weight);
546  } // end i
547 
548  } // end loop quadrature points
549  } // end alpha_skeleton
550 
551  // jacobian of skeleton term
552  template<typename IG, typename LFSU, typename X, typename LFSV,
553  typename LocalMatrix>
554  void jacobian_skeleton (const IG& ig,
555  const LFSU& lfsu_s, const X&, const LFSV& lfsv_s,
556  const LFSU& lfsu_n, const X&, const LFSV& lfsv_n,
557  LocalMatrix& mat_ss, LocalMatrix& mat_sn,
558  LocalMatrix& mat_ns, LocalMatrix& mat_nn) const
559  {
560  // dimensions
561  const unsigned int dim = IG::Geometry::dimension;
562  const unsigned int dimw = IG::Geometry::dimensionworld;
563 
564  // subspaces
565  static_assert
566  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
567 
568  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
569  const LFSV_PFS_V& lfsv_s_pfs_v = lfsv_s.template child<VBLOCK>();
570  const LFSV_PFS_V& lfsv_n_pfs_v = lfsv_n.template child<VBLOCK>();
571  static_assert
572  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
573 
574  // ... we assume all velocity components are the same
575  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
576  const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.template child<0>();
577  const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.template child<0>();
578  const unsigned int vsize_s = lfsv_s_v.size();
579  const unsigned int vsize_n = lfsv_n_v.size();
580  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
581  const LFSV_P& lfsv_s_p = lfsv_s.template child<PBLOCK>();
582  const LFSV_P& lfsv_n_p = lfsv_n.template child<PBLOCK>();
583  const unsigned int psize_s = lfsv_s_p.size();
584  const unsigned int psize_n = lfsv_n_p.size();
585 
586  // domain and range field type
587  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
588  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
589  typedef typename BasisSwitch_V::DomainField DF;
590  typedef typename BasisSwitch_V::Range RT;
591  typedef typename BasisSwitch_V::RangeField RF;
592  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
593 
594  // make copy of inside and outside cell w.r.t. the intersection
595  auto inside_cell = ig.inside();
596  auto outside_cell = ig.outside();
597 
598  // select quadrature rule
599  Dune::GeometryType gtface = ig.geometry().type();
600  const int v_order = FESwitch_V::basis(lfsv_s_v.finiteElement()).order();
601  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
602  const int qorder = 2*v_order + det_jac_order + superintegration_order;
603  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
604 
605  const int epsilon = prm.epsilonIPSymmetryFactor();
606  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
607 
608  // loop over quadrature points and integrate normal flux
609  for (const auto& ip : rule)
610  {
611 
612  // position of quadrature point in local coordinates of element
613  Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(ip.position());
614  Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(ip.position());
615 
616  const RF penalty_factor = prm.getFaceIP(ig,ip.position());
617 
618  // value of velocity shape functions
619  std::vector<RT> phi_v_s(vsize_s);
620  std::vector<RT> phi_v_n(vsize_n);
621  FESwitch_V::basis(lfsv_s_v.finiteElement()).evaluateFunction(local_s,phi_v_s);
622  FESwitch_V::basis(lfsv_n_v.finiteElement()).evaluateFunction(local_n,phi_v_n);
623  // and value of pressure shape functions
624  std::vector<RT> phi_p_s(psize_s);
625  std::vector<RT> phi_p_n(psize_n);
626  FESwitch_P::basis(lfsv_s_p.finiteElement()).evaluateFunction(local_s,phi_p_s);
627  FESwitch_P::basis(lfsv_n_p.finiteElement()).evaluateFunction(local_n,phi_p_n);
628 
629  // compute gradients
630  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_s(vsize_s);
631  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_s_v.finiteElement()),
632  inside_cell.geometry(), local_s, grad_phi_v_s);
633 
634  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v_n(vsize_n);
635  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_n_v.finiteElement()),
636  outside_cell.geometry(), local_n, grad_phi_v_n);
637 
638  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
639  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
640  const RF mu = prm.mu(ig,ip.position());
641 
642  assert(vsize_s == vsize_n);
643  const RF factor = mu * weight;
644 
645  for(unsigned int d = 0; d < dim; ++d) {
646  const LFSV_V& lfsv_s_v = lfsv_s_pfs_v.child(d);
647  const LFSV_V& lfsv_n_v = lfsv_n_pfs_v.child(d);
648 
649  //================================================//
650  // - (\mu \int < \nabla u > . normal . [v])
651  // \mu \int \frac{\sigma}{|\gamma|^\beta} [u] \cdot [v]
652  //================================================//
653  for(unsigned int i=0; i<vsize_s; ++i) {
654 
655  for(unsigned int j=0; j<vsize_s; ++j) {
656  RF val = (0.5*(grad_phi_v_s[i][0]*normal)*phi_v_s[j]) * factor;
657  mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v,i, -val);
658  mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, epsilon * val);
659  mat_ss.accumulate(lfsv_s_v,i,lfsv_s_v,j, phi_v_s[i] * phi_v_s[j] * penalty_factor * weight);
660 
661  // Assemble symmetric part for (grad u)^T
662  if(full_tensor) {
663  for(unsigned int dd = 0; dd < dim; ++dd) {
664  RF Tval = (0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
665  const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
666  mat_ss.accumulate(lfsv_s_v,j,lfsv_s_v_dd,i, - Tval);
667  mat_ss.accumulate(lfsv_s_v_dd,i,lfsv_s_v,j, epsilon*Tval );
668  }
669  }
670  }
671 
672  for(unsigned int j=0; j<vsize_n; ++j) {
673  // the normal vector flipped, thus the sign flips
674  RF val = (-0.5*(grad_phi_v_s[i][0]*normal)*phi_v_n[j]) * factor;
675  mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i,- val);
676  mat_sn.accumulate(lfsv_s_v,i,lfsv_n_v,j, epsilon*val);
677  mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v,i, -phi_v_s[i] * phi_v_n[j] * penalty_factor * weight);
678 
679  // Assemble symmetric part for (grad u)^T
680  if(full_tensor) {
681  for (unsigned int dd=0;dd<dim;++dd) {
682  RF Tval = (-0.5*(grad_phi_v_s[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
683  const LFSV_V& lfsv_s_v_dd = lfsv_s_pfs_v.child(dd);
684  mat_ns.accumulate(lfsv_n_v,j,lfsv_s_v_dd,i,- Tval);
685  mat_sn.accumulate(lfsv_s_v_dd,i,lfsv_n_v,j, epsilon*Tval);
686  }
687  }
688  }
689 
690  for(unsigned int j=0; j<vsize_s; ++j) {
691  RF val = (0.5*(grad_phi_v_n[i][0]*normal)*phi_v_s[j]) * factor;
692  mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, - val);
693  mat_ns.accumulate(lfsv_n_v,i,lfsv_s_v,j, epsilon*val );
694  mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v,i, -phi_v_n[i] * phi_v_s[j] * penalty_factor * weight);
695 
696  // Assemble symmetric part for (grad u)^T
697  if(full_tensor) {
698  for (unsigned int dd=0;dd<dim;++dd) {
699  RF Tval = (0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_s[j]) * factor;
700  const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
701  mat_sn.accumulate(lfsv_s_v,j,lfsv_n_v_dd,i, - Tval);
702  mat_ns.accumulate(lfsv_n_v_dd,i,lfsv_s_v,j, epsilon*Tval );
703  }
704  }
705  }
706 
707  for(unsigned int j=0; j<vsize_n; ++j) {
708  // the normal vector flipped, thus the sign flips
709  RF val = (-0.5*(grad_phi_v_n[i][0]*normal)*phi_v_n[j]) * factor;
710  mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, - val);
711  mat_nn.accumulate(lfsv_n_v,i,lfsv_n_v,j, epsilon*val);
712  mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v,i, phi_v_n[i] * phi_v_n[j] * penalty_factor * weight);
713 
714  // Assemble symmetric part for (grad u)^T
715  if(full_tensor) {
716  for (unsigned int dd=0;dd<dim;++dd) {
717  RF Tval = (-0.5*(grad_phi_v_n[i][0][d]*normal[dd])*phi_v_n[j]) * factor;
718  const LFSV_V& lfsv_n_v_dd = lfsv_n_pfs_v.child(dd);
719  mat_nn.accumulate(lfsv_n_v,j,lfsv_n_v_dd,i,- Tval);
720  mat_nn.accumulate(lfsv_n_v_dd,i,lfsv_n_v,j, epsilon*Tval);
721  }
722  }
723  }
724 
725  //================================================//
726  // \int <q> [u] n
727  // \int <p> [v] n
728  //================================================//
729  for(unsigned int j=0; j<psize_s; ++j) {
730  RF val = 0.5*(phi_p_s[j]*normal[d]*phi_v_s[i]) * weight;
731  mat_ss.accumulate(lfsv_s_v,i,lfsv_s_p,j, val);
732  mat_ss.accumulate(lfsv_s_p,j,lfsv_s_v,i, val * incomp_scaling);
733  }
734 
735  for(unsigned int j=0; j<psize_n; ++j) {
736  RF val = 0.5*(phi_p_n[j]*normal[d]*phi_v_s[i]) * weight;
737  mat_sn.accumulate(lfsv_s_v,i,lfsv_n_p,j, val);
738  mat_ns.accumulate(lfsv_n_p,j,lfsv_s_v,i, val * incomp_scaling);
739  }
740 
741  for (unsigned int j=0; j<psize_s;++j) {
742  // the normal vector flipped, thus the sign flips
743  RF val = -0.5*(phi_p_s[j]*normal[d]*phi_v_n[i]) * weight;
744  mat_ns.accumulate(lfsv_n_v,i,lfsv_s_p,j, val);
745  mat_sn.accumulate(lfsv_s_p,j,lfsv_n_v,i, val * incomp_scaling);
746  }
747 
748  for (unsigned int j=0; j<psize_n;++j) {
749  // the normal vector flipped, thus the sign flips
750  RF val = -0.5*(phi_p_n[j]*normal[d]*phi_v_n[i]) * weight;
751  mat_nn.accumulate(lfsv_n_v,i,lfsv_n_p,j, val);
752  mat_nn.accumulate(lfsv_n_p,j,lfsv_n_v,i, val * incomp_scaling);
753  }
754  } // end i
755  } // end d
756 
757  } // end loop quadrature points
758  } // end jacobian_skeleton
759 
760  // boundary term
761  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
762  void alpha_boundary (const IG& ig,
763  const LFSU& lfsu, const X& x, const LFSV& lfsv,
764  R& r) const
765  {
766  // dimensions
767  const unsigned int dim = IG::Geometry::dimension;
768  const unsigned int dimw = IG::Geometry::dimensionworld;
769 
770  // subspaces
771  static_assert
772  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
773 
774  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
775  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
776  static_assert
777  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
778 
779  // ... we assume all velocity components are the same
780  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
781  const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
782  const unsigned int vsize = lfsv_v.size();
783  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
784  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
785  const unsigned int psize = lfsv_p.size();
786 
787  // domain and range field type
788  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
789  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
790  typedef typename BasisSwitch_V::DomainField DF;
791  typedef typename BasisSwitch_V::Range RT;
792  typedef typename BasisSwitch_V::RangeField RF;
793  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
794 
795  // make copy of inside cell w.r.t. the boundary
796  auto inside_cell = ig.inside();
797 
798  // select quadrature rule
799  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
800  Dune::GeometryType gtface = ig.geometry().type();
801  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
802  const int qorder = 2*v_order + det_jac_order + superintegration_order;
803  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
804 
805  const int epsilon = prm.epsilonIPSymmetryFactor();
806  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
807 
808  // loop over quadrature points and integrate normal flux
809  for (const auto& ip : rule)
810  {
811  // position of quadrature point in local coordinates of element
812  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
813 
814  const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
815 
816  // value of velocity shape functions
817  std::vector<RT> phi_v(vsize);
818  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
819 
820  // evaluate u
821  Dune::FieldVector<RF,dim> u(0.0);
822  for(unsigned int d=0; d<dim; ++d) {
823  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
824  for(unsigned int i=0; i<vsize; i++)
825  u[d] += x(lfsv_v,i) * phi_v[i];
826  }
827 
828  // value of pressure shape functions
829  std::vector<RT> phi_p(psize);
830  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
831 
832  // evaluate pressure
833  RF p(0.0);
834  for(unsigned int i=0; i<psize; i++)
835  p += x(lfsv_p,i) * phi_p[i];
836 
837  // compute gradients
838  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
839  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
840  inside_cell.geometry(), local, grad_phi_v);
841 
842  // evaluate velocity jacobian
843  Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
844  for(unsigned int d=0; d<dim; ++d) {
845  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
846  for(unsigned int i=0; i<vsize; i++)
847  jacu[d].axpy(x(lfsv_v,i), grad_phi_v[i][0]);
848  }
849 
850  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
851  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
852  const RF mu = prm.mu(ig,ip.position());
853 
854  // evaluate boundary condition type
855  typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
856 
857  // Slip factor smoothly switching between slip and no slip conditions.
858  RF slip_factor = 0.0;
859  typedef NavierStokesDGImp::VariableBoundarySlipSwitch<PRM> BoundarySlipSwitch;
860  if (bctype == BC::SlipVelocity)
861  // Calls boundarySlip(..) function of parameter
862  // class if available, i.e. if
863  // enable_variable_slip is defined. Otherwise
864  // returns 1.0;
865  slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
866 
867  // velocity boundary condition
868  if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
869  {
870  // on BC::VelocityDirichlet: 1.0 - slip_factor = 1.0
871  const RF factor = weight * (1.0 - slip_factor);
872 
873  for(unsigned int d = 0; d < dim; ++d) {
874  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
875 
876  RF val = (jacu[d] * normal) * factor * mu;
877  for(unsigned int i=0; i<vsize; i++) {
878  //================================================//
879  // - (\mu \int \nabla u. normal . v)
880  //================================================//
881  r.accumulate(lfsv_v,i, -val * phi_v[i]);
882  r.accumulate(lfsv_v,i, epsilon * mu * (grad_phi_v[i][0] * normal) * u[d] * factor);
883 
884  if(full_tensor) {
885  for(unsigned int dd=0; dd<dim; ++dd) {
886  r.accumulate(lfsv_v,i, -mu * jacu[dd][d] * normal[dd] * phi_v[i] * factor);
887  r.accumulate(lfsv_v,i, epsilon * mu * grad_phi_v[i][0][dd] * u[dd] * normal[d] * factor);
888  }
889  }
890 
891  //================================================//
892  // \mu \int \sigma / |\gamma|^\beta v u
893  //================================================//
894  r.accumulate(lfsv_v,i, u[d] * phi_v[i] * penalty_factor * factor);
895 
896  //================================================//
897  // \int p v n
898  //================================================//
899  r.accumulate(lfsv_v,i, p * phi_v[i] * normal[d] * weight);
900  } // end i
901  } // end d
902 
903  //================================================//
904  // \int q u n
905  //================================================//
906  for(unsigned int i=0; i<psize; i++) {
907  r.accumulate(lfsv_p,i, phi_p[i] * (u * normal) * incomp_scaling * weight);
908  }
909  } // Velocity Dirichlet
910 
911  if(bctype == BC::SlipVelocity)
912  {
913  const RF factor = weight * (slip_factor);
914 
915  RF ten_sum = 1.0;
916  if(full_tensor)
917  ten_sum = 2.0;
918 
919  for(unsigned int d = 0; d < dim; ++d) {
920  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
921 
922  for(unsigned int i=0; i<vsize; i++) {
923  //================================================//
924  // - (\mu \int \nabla u. normal . v)
925  //================================================//
926  for(unsigned int dd = 0; dd < dim; ++dd)
927  r.accumulate(lfsv_v,i, -ten_sum * mu * (jacu[dd] * normal) * normal[dd] *phi_v[i] * normal[d] * factor);
928  r.accumulate(lfsv_v,i, epsilon * ten_sum * mu * (u * normal) * (grad_phi_v[i][0] * normal) * normal[d] * factor);
929 
930  //================================================//
931  // \mu \int \sigma / |\gamma|^\beta v u
932  //================================================//
933  r.accumulate(lfsv_v,i, penalty_factor * (u * normal) * phi_v[i] * normal[d] * factor);
934  } // end i
935  } // end d
936 
937  } // Slip Velocity
938  } // end loop quadrature points
939  } // end alpha_boundary
940 
941  // jacobian of boundary term
942  template<typename IG, typename LFSU, typename X, typename LFSV,
943  typename LocalMatrix>
944  void jacobian_boundary (const IG& ig,
945  const LFSU& lfsu, const X& x, const LFSV& lfsv,
946  LocalMatrix& mat) const
947  {
948  // dimensions
949  const unsigned int dim = IG::Geometry::dimension;
950  const unsigned int dimw = IG::Geometry::dimensionworld;
951 
952  // subspaces
953  static_assert
954  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
955 
956  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
957  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
958  static_assert
959  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
960 
961  // ... we assume all velocity components are the same
962  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
963  const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
964  const unsigned int vsize = lfsv_v.size();
965  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
966  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
967  const unsigned int psize = lfsv_p.size();
968 
969  // domain and range field type
970  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
971  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
972  typedef typename BasisSwitch_V::DomainField DF;
973  typedef typename BasisSwitch_V::Range RT;
974  typedef typename BasisSwitch_V::RangeField RF;
975  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
976 
977  // make copy of inside cell w.r.t. the boundary
978  auto inside_cell = ig.inside();
979 
980  // select quadrature rule
981  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
982  Dune::GeometryType gtface = ig.geometry().type();
983  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-1);
984  const int qorder = 2*v_order + det_jac_order + superintegration_order;
985  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
986 
987  const int epsilon = prm.epsilonIPSymmetryFactor();
988  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
989 
990  // loop over quadrature points and integrate normal flux
991  for (const auto& ip : rule)
992  {
993  // position of quadrature point in local coordinates of element
994  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
995 
996  const RF penalty_factor = prm.getFaceIP(ig,ip.position() );
997 
998  // value of velocity shape functions
999  std::vector<RT> phi_v(vsize);
1000  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1001  // and value of pressure shape functions
1002  std::vector<RT> phi_p(psize);
1003  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1004 
1005  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1006  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1007  inside_cell.geometry(), local, grad_phi_v);
1008 
1009  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
1010  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
1011  const RF mu = prm.mu(ig,ip.position());
1012 
1013  // evaluate boundary condition type
1014  typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,ip.position()));
1015 
1016  // Slip factor smoothly switching between slip and no slip conditions.
1017  RF slip_factor = 0.0;
1018  typedef NavierStokesDGImp::VariableBoundarySlipSwitch<PRM> BoundarySlipSwitch;
1019  if (bctype == BC::SlipVelocity)
1020  // Calls boundarySlip(..) function of parameter
1021  // class if available, i.e. if
1022  // enable_variable_slip is defined. Otherwise
1023  // returns 1.0;
1024  slip_factor = BoundarySlipSwitch::boundarySlip(prm,ig,ip.position());
1025 
1026  // velocity boundary condition
1027  if (bctype == BC::VelocityDirichlet || bctype == BC::SlipVelocity)
1028  {
1029  // on BC::VelocityDirichlet: 1.0 - slip_factor = 1.0
1030  const RF factor = weight * (1.0 - slip_factor);
1031 
1032  for(unsigned int d = 0; d < dim; ++d) {
1033  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1034 
1035  for(unsigned int i=0; i<vsize; i++) {
1036 
1037  for(unsigned int j=0; j<vsize; j++) {
1038  //================================================//
1039  // - (\mu \int \nabla u. normal . v)
1040  //================================================//
1041  RF val = ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1042  mat.accumulate(lfsv_v,i,lfsv_v,j, - val);
1043  mat.accumulate(lfsv_v,j,lfsv_v,i, epsilon*val);
1044 
1045  // Assemble symmetric part for (grad u)^T
1046  if(full_tensor) {
1047  for(unsigned int dd = 0; dd < dim; ++dd) {
1048  const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1049  RF Tval = ((grad_phi_v[j][0][d]*normal[dd])*phi_v[i]) * factor * mu;
1050  mat.accumulate(lfsv_v,i,lfsv_v_dd,j, - Tval);
1051  mat.accumulate(lfsv_v_dd,j,lfsv_v,i, epsilon*Tval);
1052  }
1053  }
1054  //================================================//
1055  // \mu \int \sigma / |\gamma|^\beta v u
1056  //================================================//
1057  mat.accumulate(lfsv_v,j,lfsv_v,i, phi_v[i] * phi_v[j] * penalty_factor * factor);
1058  }
1059 
1060  //================================================//
1061  // \int q u n
1062  // \int p v n
1063  //================================================//
1064  for(unsigned int j=0; j<psize; j++) {
1065  mat.accumulate(lfsv_p,j,lfsv_v,i, (phi_p[j]*normal[d]*phi_v[i]) * weight * incomp_scaling);
1066  mat.accumulate(lfsv_v,i,lfsv_p,j, (phi_p[j]*normal[d]*phi_v[i]) * weight);
1067  }
1068  } // end i
1069  } // end d
1070 
1071  } // Velocity Dirichlet
1072 
1073  if (bctype == BC::SlipVelocity)
1074  {
1075  const RF factor = weight * (slip_factor);
1076 
1077  //================================================//
1078  // - (\mu \int \nabla u. normal . v)
1079  //================================================//
1080 
1081  for (unsigned int i=0;i<vsize;++i) // ansatz
1082  {
1083  for (unsigned int j=0;j<vsize;++j) // test
1084  {
1085  RF ten_sum = 1.0;
1086 
1087  // Assemble symmetric part for (grad u)^T
1088  if(full_tensor)
1089  ten_sum = 2.0;
1090 
1091  RF val = ten_sum * ((grad_phi_v[j][0]*normal)*phi_v[i]) * factor * mu;
1092  for (unsigned int d=0;d<dim;++d)
1093  {
1094  const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
1095 
1096  for (unsigned int dd=0;dd<dim;++dd)
1097  {
1098  const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1099 
1100  mat.accumulate(lfsv_v_dd,i,lfsv_v_d,j, -val*normal[d]*normal[dd]);
1101  mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, epsilon*val*normal[d]*normal[dd]);
1102  }
1103  }
1104  }
1105  }
1106 
1107  //================================================//
1108  // \mu \int \sigma / |\gamma|^\beta v u
1109  //================================================//
1110  const RF p_factor = penalty_factor * factor;
1111  for (unsigned int i=0;i<vsize;++i)
1112  {
1113  for (unsigned int j=0;j<vsize;++j)
1114  {
1115  RF val = phi_v[i]*phi_v[j] * p_factor;
1116  for (unsigned int d=0;d<dim;++d)
1117  {
1118  const LFSV_V& lfsv_v_d = lfsv_pfs_v.child(d);
1119  for (unsigned int dd=0;dd<dim;++dd)
1120  {
1121  const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1122  mat.accumulate(lfsv_v_d,j,lfsv_v_dd,i, val*normal[d]*normal[dd]);
1123  }
1124  }
1125  }
1126  }
1127 
1128  } // Slip Velocity
1129  } // end loop quadrature points
1130  } // end jacobian_boundary
1131 
1132  // volume integral depending only on test functions,
1133  // contains f on the right hand side
1134  template<typename EG, typename LFSV, typename R>
1135  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
1136  {
1137  // dimensions
1138  static const unsigned int dim = EG::Geometry::dimension;
1139 
1140  // subspaces
1141  static_assert
1142  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
1143 
1144  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1145  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1146 
1147  static_assert
1148  ((LFSV_PFS_V::CHILDREN == dim),"You seem to use the wrong function space for DGNavierStokes");
1149 
1150  // we assume all velocity components are the same type
1151  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1152  const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1153  const unsigned int vsize = lfsv_v.size();
1154  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
1155  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
1156  const unsigned int psize = lfsv_p.size();
1157 
1158  // domain and range field type
1159  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1160  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1161  typedef typename BasisSwitch_V::DomainField DF;
1162  typedef typename BasisSwitch_V::Range RT;
1163  typedef typename BasisSwitch_V::RangeField RF;
1164  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
1165  typedef typename LFSV::Traits::SizeType size_type;
1166 
1167  // select quadrature rule
1168  Dune::GeometryType gt = eg.geometry().type();
1169  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1170  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
1171  const int qorder = 2*v_order + det_jac_order + superintegration_order;
1172 
1173  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
1174 
1175  // loop over quadrature points
1176  for (const auto& ip : rule)
1177  {
1178  const Dune::FieldVector<DF,dim> local = ip.position();
1179  //const Dune::FieldVector<DF,dimw> global = eg.geometry().global(local);
1180 
1181  // values of velocity shape functions
1182  std::vector<RT> phi_v(vsize);
1183  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1184 
1185  // values of pressure shape functions
1186  std::vector<RT> phi_p(psize);
1187  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1188 
1189  const RF weight = ip.weight() * eg.geometry().integrationElement(ip.position());
1190 
1191  // evaluate source term
1192  typename PRM::Traits::VelocityRange fval(prm.f(eg,local));
1193 
1194  //================================================//
1195  // \int (f*v)
1196  //================================================//
1197  const RF factor = weight;
1198  for (unsigned int d=0; d<dim; d++) {
1199  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1200 
1201  // and store for each velocity component
1202  for (size_type i=0; i<vsize; i++) {
1203  RF val = phi_v[i]*factor;
1204  r.accumulate(lfsv_v,i, -fval[d] * val);
1205  }
1206  }
1207 
1208  const RF g2 = prm.g2(eg,ip.position());
1209 
1210  // integrate div u * psi_i
1211  for (size_t i=0; i<lfsv_p.size(); i++) {
1212  r.accumulate(lfsv_p,i, g2 * phi_p[i] * factor);
1213  }
1214 
1215  } // end loop quadrature points
1216  } // end lambda_volume
1217 
1218  // boundary integral independent of ansatz functions,
1219  // Neumann and Dirichlet boundary conditions, DG penalty term's right hand side
1220  template<typename IG, typename LFSV, typename R>
1221  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
1222  {
1223  // dimensions
1224  static const unsigned int dim = IG::Geometry::dimension;
1225 
1226  // subspaces
1227  static_assert
1228  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for DGNavierStokes");
1229 
1230  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_PFS_V;
1231  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<VBLOCK>();
1232 
1233  static_assert
1234  ((LFSV_PFS_V::CHILDREN == dim), "You seem to use the wrong function space for DGNavierStokes");
1235 
1236  // ... we assume all velocity components are the same
1237  typedef typename LFSV_PFS_V::template Child<0>::Type LFSV_V;
1238  const LFSV_V& lfsv_v = lfsv_pfs_v.template child<0>();
1239  const unsigned int vsize = lfsv_v.size();
1240  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
1241  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
1242  const unsigned int psize = lfsv_p.size();
1243 
1244  // domain and range field type
1245  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
1246  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
1247  typedef typename BasisSwitch_V::DomainField DF;
1248  typedef typename BasisSwitch_V::Range RT;
1249  typedef typename BasisSwitch_V::RangeField RF;
1250  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
1251 
1252  // make copy of inside cell w.r.t. the boundary
1253  auto inside_cell = ig.inside();
1254 
1255  // select quadrature rule
1256  Dune::GeometryType gtface = ig.geometry().type();
1257  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1258  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
1259  const int jac_order = gtface.isSimplex() ? 0 : 1;
1260  const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
1261  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
1262 
1263  const int epsilon = prm.epsilonIPSymmetryFactor();
1264  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
1265 
1266  // loop over quadrature points and integrate normal flux
1267  for (const auto& ip : rule)
1268  {
1269  // position of quadrature point in local coordinates of element
1270  Dune::FieldVector<DF,dim-1> flocal = ip.position();
1271  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(flocal);
1272  //Dune::FieldVector<DF,dimw> global = ig.geometry().global(flocal);
1273 
1274  const RF penalty_factor = prm.getFaceIP(ig,flocal);
1275 
1276  // value of velocity shape functions
1277  std::vector<RT> phi_v(vsize);
1278  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1279  // and value of pressure shape functions
1280  std::vector<RT> phi_p(psize);
1281  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1282 
1283  std::vector<Dune::FieldMatrix<RF,1,dim> > grad_phi_v(vsize);
1284  BasisSwitch_V::gradient(FESwitch_V::basis(lfsv_v.finiteElement()),
1285  inside_cell.geometry(), local, grad_phi_v);
1286 
1287  const Dune::FieldVector<DF,dim> normal = ig.unitOuterNormal(ip.position());
1288  const RF weight = ip.weight()*ig.geometry().integrationElement(ip.position());
1289  const RF mu = prm.mu(ig,flocal);
1290 
1291  // evaluate boundary condition type
1292  typename PRM::Traits::BoundaryCondition::Type bctype(prm.bctype(ig,flocal));
1293 
1294  if (bctype == BC::VelocityDirichlet)
1295  {
1296  typename PRM::Traits::VelocityRange u0(prm.g(inside_cell,local));
1297 
1298  RF factor = mu * weight;
1299  for(unsigned int d = 0; d < dim; ++d) {
1300  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1301 
1302  for(unsigned int i=0; i<vsize; i++) {
1303  //================================================//
1304  // \mu \int \nabla v \cdot u_0 \cdot n
1305  //================================================//
1306  r.accumulate(lfsv_v,i, -epsilon * (grad_phi_v[i][0] * normal) * factor * u0[d]);
1307 
1308  // Assemble symmetric part for (grad u)^T
1309  if(full_tensor) {
1310  for(unsigned int dd = 0; dd < dim; ++dd) {
1311  const LFSV_V& lfsv_v_dd = lfsv_pfs_v.child(dd);
1312  RF Tval = (grad_phi_v[i][0][d]*normal[dd]) * factor;
1313  r.accumulate(lfsv_v_dd,i, -epsilon * Tval * u0[d]);
1314  }
1315  }
1316  //================================================//
1317  // \int \sigma / |\gamma|^\beta v u_0
1318  //================================================//
1319  r.accumulate(lfsv_v,i, -phi_v[i] * penalty_factor * u0[d] * weight);
1320 
1321  } // end i
1322  } // end d
1323 
1324  //================================================//
1325  // \int q u_0 n
1326  //================================================//
1327  for (unsigned int i=0;i<psize;++i) // test
1328  {
1329  RF val = phi_p[i]*(u0 * normal) * weight;
1330  r.accumulate(lfsv_p,i, - val * incomp_scaling);
1331  }
1332  } // end BC velocity
1333  if (bctype == BC::StressNeumann)
1334  {
1335  typename PRM::Traits::VelocityRange stress(prm.j(ig,flocal,normal));
1336 
1337  //std::cout << "Pdirichlet\n";
1338  //================================================//
1339  // \int p u n
1340  //================================================//
1341  for(unsigned int d = 0; d < dim; ++d) {
1342  const LFSV_V& lfsv_v = lfsv_pfs_v.child(d);
1343 
1344  for(unsigned int i=0; i<vsize; i++)
1345  r.accumulate(lfsv_v,i, stress[d] * phi_v[i] * weight);
1346  }
1347  }
1348  } // end loop quadrature points
1349  } // end lambda_boundary
1350 
1351  private :
1352  PRM& prm;
1353  const int superintegration_order;
1354  Real current_dt;
1355  }; // end class DGNavierStokes
1356 
1357  } // end namespace PDELab
1358 } // end namespace Dune
1359 #endif
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:1135
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:762
Definition: dgnavierstokes.hh:61
void setTime(doublet_)
set time for subsequent evaluation
Definition: idefault.hh:104
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
DGNavierStokes(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokes.hh:78
Definition: dgnavierstokes.hh:57
void setTime(Real t)
Definition: dgnavierstokes.hh:90
const IG & ig
Definition: constraints.hh:147
Definition: dgnavierstokes.hh:60
Definition: dgnavierstokes.hh:62
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:1221
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: dgnavierstokes.hh:554
Definition: dgnavierstokes.hh:54
Definition: adaptivity.hh:27
static const int dim
Definition: adaptivity.hh:83
Definition: dgnavierstokes.hh:63
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: dgnavierstokes.hh:363
const P & p
Definition: constraints.hh:146
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokes.hh:944
sparsity pattern generator
Definition: pattern.hh:29
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokes.hh:229
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokes.hh:98
Definition: dgnavierstokes.hh:64
A local operator for solving the Navier-Stokes equations using a DG discretization.
Definition: dgnavierstokes.hh:36
Definition: dgnavierstokes.hh:53
Compile-time switch for the boundary slip factor.
Definition: dgnavierstokesparameter.hh:159
void preStep(RealType, RealType dt, int)
Definition: dgnavierstokes.hh:84
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89