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