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