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