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