dune-pdelab  2.5-dev
dgnavierstokesvelvecfem.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_DGNAVIERSTOKESVELVECFEM_HH
5 #define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
6 
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fvector.hh>
9 #include <dune/common/parametertree.hh>
10 
11 #include <dune/geometry/quadraturerules.hh>
12 
13 #include <dune/localfunctions/common/interfaceswitch.hh>
15 
21 
22 #ifndef VBLOCK
23 #define VBLOCK 0
24 #endif
25 #define PBLOCK (- VBLOCK + 1)
26 
27 namespace Dune {
28  namespace PDELab {
29 
30  template<class Basis, class Dummy = void>
33  using DomainLocal = typename Basis::Traits::DomainLocal;
35  using RangeField = typename Basis::Traits::RangeField;
37  static const std::size_t dimRange = Basis::Traits::dimRange;
38 
40  template<typename Geometry>
41  static void jacobian(const Basis& basis, const Geometry& geometry,
42  const DomainLocal& xl,
43  std::vector<FieldMatrix<RangeField, dimRange,
44  Geometry::coorddimension> >& jac)
45  {
46  jac.resize(basis.size());
47  basis.evaluateJacobian(xl, jac);
48  }
49  };
50 
52  template<class Basis>
54  Basis, typename std::enable_if<
55  Std::to_true_type<
56  std::integral_constant<
57  std::size_t,
58  Basis::Traits::dimDomain
59  >
60  >::value
61  >::type
62  >
63  {
65  using DomainLocal = typename Basis::Traits::DomainType;
67  using RangeField = typename Basis::Traits::RangeFieldType;
69  static const std::size_t dimRange = Basis::Traits::dimRange;
70 
72  template<typename Geometry>
73  static void jacobian(const Basis& basis, const Geometry& geometry,
74  const DomainLocal& xl,
75  std::vector<FieldMatrix<RangeField, dimRange,
76  Geometry::coorddimension> >& jac)
77  {
78  jac.resize(basis.size());
79 
80  std::vector<FieldMatrix<
81  RangeField, dimRange, Geometry::coorddimension> > ljac(basis.size());
82  basis.evaluateJacobian(xl, ljac);
83 
84  const typename Geometry::JacobianInverseTransposed& geojac =
85  geometry.jacobianInverseTransposed(xl);
86 
87  for(std::size_t i = 0; i < basis.size(); ++i)
88  for(std::size_t row=0; row < dimRange; ++row)
89  geojac.mv(ljac[i][row], jac[i][row]);
90  }
91  };
92 
100  template<typename PRM>
103  public FullSkeletonPattern, public FullVolumePattern,
104  public InstationaryLocalOperatorDefaultMethods<typename PRM::Traits::RangeField>
105  {
106  using BC = StokesBoundaryCondition;
107  using RF = typename PRM::Traits::RangeField;
108 
110  using Real = typename InstatBase::RealType;
111 
112  static const bool navier = PRM::assemble_navier;
113  static const bool full_tensor = PRM::assemble_full_tensor;
114 
115  public :
116 
117  // pattern assembly flags
118  enum { doPatternVolume = true };
119  enum { doPatternSkeleton = true };
120 
121  // call the assembler for each face only once
122  enum { doSkeletonTwoSided = false };
123 
124  // residual assembly flags
125  enum { doAlphaVolume = true };
126  enum { doAlphaSkeleton = true };
127  enum { doAlphaBoundary = true };
128  enum { doLambdaVolume = true };
129 
142  DGNavierStokesVelVecFEM (PRM& _prm, int _superintegration_order=0) :
143  prm(_prm), superintegration_order(_superintegration_order),
144  current_dt(1.0)
145  {}
146 
147  // Store current dt
148  void preStep (Real , Real dt, int )
149  {
150  current_dt = dt;
151  }
152 
153  // set time in parameter class
154  void setTime(Real t)
155  {
157  prm.setTime(t);
158  }
159 
160  // volume integral depending on test and ansatz functions
161  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
162  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
163  {
164  const unsigned int dim = EG::Geometry::mydimension;
165 
166  // subspaces
167  using namespace Indices;
168  using LFSV_V = TypeTree::Child<LFSV,_0>;
169  const auto& lfsv_v = child(lfsv,_0);
170  const auto& lfsu_v = child(lfsu,_0);
171 
172  using LFSV_P = TypeTree::Child<LFSV,_1>;
173  const auto& lfsv_p = child(lfsv,_1);
174  const auto& lfsu_p = child(lfsu,_1);
175 
176  // domain and range field type
177  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
178  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
180  using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
181  using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
182  using RF = typename BasisSwitch_V::RangeField;
183  using Range_V = typename BasisSwitch_V::Range;
184  using Range_P = typename BasisSwitch_P::Range;
185  using size_type = typename LFSV::Traits::SizeType;
186 
187  // Get geometry
188  auto geo = eg.geometry();
189 
190  // Determine quadrature order
191  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
192  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
193  const int jac_order = geo.type().isSimplex() ? 0 : 1;
194  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
195 
196  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
197 
198  // loop over quadrature points
199  for (const auto& ip : quadratureRule(geo,qorder))
200  {
201  auto local = ip.position();
202  auto mu = prm.mu(eg,local);
203  auto rho = prm.rho(eg,local);
204 
205  // compute u (if Navier term enabled)
206  std::vector<Range_V> phi_v(lfsv_v.size());
207  Range_V val_u(0.0);
208  if(navier) {
209  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
210  for (size_type i=0; i<lfsu_v.size(); i++)
211  val_u.axpy(x(lfsu_v,i),phi_v[i]);
212  }
213 
214  // values of pressure shape functions
215  std::vector<Range_P> phi_p(lfsv_p.size());
216  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
217 
218  // compute pressure value
219  Range_P val_p(0.0);
220  for (size_type i=0; i<lfsu_p.size(); i++)
221  val_p.axpy(x(lfsu_p,i),phi_p[i]);
222 
223  // evaluate jacobian of velocity shape functions on reference element
224  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
225  VectorBasisSwitch_V::jacobian
226  (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
227 
228  // compute divergence of test functions
229  std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
230  for (size_type i=0; i<lfsv_v.size(); i++)
231  for (size_type d=0; d<dim; d++)
232  div_phi_v[i] += jac_phi_v[i][d][d];
233 
234  // compute velocity jacobian and divergence
235  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
236  RF div_u(0.0);
237  for (size_type i=0; i<lfsu_v.size(); i++){
238  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
239  div_u += x(lfsu_v,i) * div_phi_v[i];
240  }
241 
242  auto detj = geo.integrationElement(ip.position());
243  auto weight = ip.weight() * detj;
244 
245  for (size_type i=0; i<lfsv_v.size(); i++) {
246  //================================================//
247  // \int (mu*grad_u*grad_v)
248  //================================================//
249  RF dvdu(0);
250  contraction(jac_u,jac_phi_v[i],dvdu);
251  r.accumulate(lfsv_v, i, dvdu * mu * weight);
252 
253  //================================================//
254  // \int -p \nabla\cdot v
255  //================================================//
256  r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
257 
258  //================================================//
259  // \int \rho ((u\cdot\nabla ) u )\cdot v
260  //================================================//
261  if(navier) {
262  // compute (grad u) u (matrix-vector product)
263  Range_V nabla_u_u(0.0);
264  jac_u.mv(val_u,nabla_u_u);
265  r.accumulate(lfsv_v, i, rho * (nabla_u_u*phi_v[i]) * weight);
266  } // end navier
267 
268  } // end i
269 
270  for (size_type i=0; i<lfsv_p.size(); i++) {
271  //================================================//
272  // \int -q \nabla\cdot u
273  //================================================//
274  r.accumulate(lfsv_p, i, - div_u * phi_p[i] * incomp_scaling * weight);
275  }
276 
277  } // end loop quadrature points
278  } // end alpha_volume
279 
280  // jacobian of volume term
281  template<typename EG, typename LFSU, typename X, typename LFSV,
282  typename LocalMatrix>
283  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
284  LocalMatrix& mat) const
285  {
286  const unsigned int dim = EG::Geometry::mydimension;
287 
288  // subspaces
289  using namespace Indices;
290  using LFSV_V = TypeTree::Child<LFSV,_0>;
291  const auto& lfsv_v = child(lfsv,_0);
292  const auto& lfsu_v = child(lfsu,_0);
293 
294  using LFSV_P = TypeTree::Child<LFSV,_1>;
295  const auto& lfsv_p = child(lfsv,_1);
296  const auto& lfsu_p = child(lfsu,_1);
297 
298  // domain and range field type
299  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
300  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
302  using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
303  using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
304  using RF = typename BasisSwitch_V::RangeField;
305  using Range_V = typename BasisSwitch_V::Range;
306  using Range_P = typename BasisSwitch_P::Range;
307  using size_type = typename LFSV::Traits::SizeType;
308 
309  // Get geometry
310  auto geo = eg.geometry();
311 
312  // Determine quadrature order
313  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
314  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
315  const int jac_order = geo.type().isSimplex() ? 0 : 1;
316  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
317 
318  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
319 
320  // loop over quadrature points
321  for (const auto& ip : quadratureRule(geo,qorder))
322  {
323  auto local = ip.position();
324  auto mu = prm.mu(eg,local);
325  auto rho = prm.rho(eg,local);
326 
327  // compute u (if Navier term enabled)
328  std::vector<Range_V> phi_v(lfsv_v.size());
329  Range_V val_u(0.0);
330  if(navier) {
331  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
332  for (size_type i=0; i<lfsu_v.size(); i++)
333  val_u.axpy(x(lfsu_v,i),phi_v[i]);
334  }
335 
336  // values of pressure shape functions
337  std::vector<Range_P> phi_p(lfsv_p.size());
338  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
339 
340  // evaluate jacobian of velocity shape functions on reference element
341  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
342  VectorBasisSwitch_V::jacobian
343  (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
344 
345  assert(lfsu_v.size() == lfsv_v.size());
346  // compute divergence of velocity shape functions
347  std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
348  for (size_type i=0; i<lfsv_v.size(); i++)
349  for (size_type d=0; d<dim; d++)
350  div_phi_v[i] += jac_phi_v[i][d][d];
351 
352  // compute velocity jacobian (if Navier term enabled)
353  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
354  if(navier) {
355  for (size_type i=0; i<lfsu_v.size(); i++){
356  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
357  }
358  }
359 
360  auto detj = geo.integrationElement(ip.position());
361  auto weight = ip.weight() * detj;
362 
363  for(size_type i=0; i<lfsv_v.size(); i++) {
364 
365  for(size_type j=0; j<lfsu_v.size(); j++) {
366  //================================================//
367  // \int (mu*grad_u*grad_v)
368  //================================================//
369  RF dvdu(0.0);
370  contraction(jac_phi_v[j],jac_phi_v[i],dvdu);
371  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * dvdu * weight);
372 
373  //================================================//
374  // \int \rho ((u\cdot\nabla ) u )\cdot v
375  //================================================//
376  if(navier) {
377  // compute (grad u) phi_v_j (matrix-vector product)
378  Range_V nabla_u_phi_v_j(0.0);
379  jac_u.mv(phi_v[j],nabla_u_phi_v_j);
380  // compute (grad phi_v_j) u (matrix-vector product)
381  Range_V nabla_phi_v_j_u(0.0);
382  jac_phi_v[j].mv(val_u,nabla_phi_v_j_u);
383  mat.accumulate(lfsv_v,i,lfsu_v,j, rho * ((nabla_u_phi_v_j*phi_v[i]) + (nabla_phi_v_j_u*phi_v[i])) * weight);
384  } // end navier
385 
386  } // end j
387 
388  for(size_type j=0; j<lfsv_p.size(); j++) {
389  //================================================//
390  // - p * div v
391  // - q * div u
392  //================================================//
393  mat.accumulate(lfsv_v,i,lfsu_p,j, -phi_p[j] * div_phi_v[i] * weight);
394  mat.accumulate(lfsv_p,j,lfsu_v,i, -phi_p[j] * div_phi_v[i] * incomp_scaling * weight);
395  }
396  } // end i
397 
398  } // end loop quadrature points
399  } // end jacobian_volume
400 
401  // skeleton term, each face is only visited ONCE
402  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
403  void alpha_skeleton (const IG& ig,
404  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
405  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
406  R& r_s, R& r_n) const
407  {
408  // dimensions
409  const unsigned int dim = IG::Entity::dimension;
410  const unsigned int dimw = IG::coorddimension;
411 
412  // subspaces
413  using namespace Indices;
414  using LFSV_V = TypeTree::Child<LFSV,_0>;
415  const auto& lfsv_v_s = child(lfsv_s,_0);
416  const auto& lfsu_v_s = child(lfsu_s,_0);
417  const auto& lfsv_v_n = child(lfsv_n,_0);
418  const auto& lfsu_v_n = child(lfsu_n,_0);
419 
420  using LFSV_P = TypeTree::Child<LFSV,_1>;
421  const auto& lfsv_p_s = child(lfsv_s,_1);
422  const auto& lfsu_p_s = child(lfsu_s,_1);
423  const auto& lfsv_p_n = child(lfsv_n,_1);
424  const auto& lfsu_p_n = child(lfsu_n,_1);
425 
426  // domain and range field type
427  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
428  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
430  using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
431  using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
432  using DF = typename BasisSwitch_V::DomainField;
433  using RF = typename BasisSwitch_V::RangeField;
434  using Range_V = typename BasisSwitch_V::Range;
435  using Range_P = typename BasisSwitch_P::Range;
436  using size_type = typename LFSV::Traits::SizeType;
437 
438  // References to inside and outside cells
439  const auto& cell_inside = ig.inside();
440  const auto& cell_outside = ig.outside();
441 
442  // Get geometries
443  auto geo = ig.geometry();
444  auto geo_inside = cell_inside.geometry();
445  auto geo_outside = cell_outside.geometry();
446 
447  // Get geometry of intersection in local coordinates of cell_inside and cell_outside
448  auto geo_in_inside = ig.geometryInInside();
449  auto geo_in_outside = ig.geometryInOutside();
450 
451  // Determine quadrature order
452  const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
453  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
454  const int qorder = 2*v_order + det_jac_order + superintegration_order;
455 
456  const int epsilon = prm.epsilonIPSymmetryFactor();
457  const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
458 
459  auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
460 
461  // loop over quadrature points and integrate normal flux
462  for (const auto& ip : quadratureRule(geo,qorder))
463  {
464 
465  // position of quadrature point in local coordinates of element
466  auto local_s = geo_in_inside.global(ip.position());
467  auto local_n = geo_in_outside.global(ip.position());
468 
469  // values of velocity shape functions
470  std::vector<Range_V> phi_v_s(lfsv_v_s.size());
471  std::vector<Range_V> phi_v_n(lfsv_v_n.size());
472  FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
473  FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
474 
475  // values of pressure shape functions
476  std::vector<Range_P> phi_p_s(lfsv_p_s.size());
477  std::vector<Range_P> phi_p_n(lfsv_p_n.size());
478  FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
479  FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
480 
481  // compute pressure value
482  Range_P val_p_s(0.0);
483  Range_P val_p_n(0.0);
484  for (size_type i=0; i<lfsu_p_s.size(); i++)
485  val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
486  for (size_type i=0; i<lfsu_p_n.size(); i++)
487  val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
488 
489  // evaluate jacobian of velocity shape functions on reference element
490  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
491  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
492  VectorBasisSwitch_V::jacobian
493  (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
494  VectorBasisSwitch_V::jacobian
495  (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
496 
497  // compute velocity value, jacobian, and divergence
498  Range_V val_u_s(0.0);
499  Range_V val_u_n(0.0);
500  Dune::FieldMatrix<RF,dim,dim> jac_u_s(0.0);
501  Dune::FieldMatrix<RF,dim,dim> jac_u_n(0.0);
502  for (size_type i=0; i<lfsu_v_s.size(); i++){
503  val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
504  jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
505  }
506  for (size_type i=0; i<lfsu_v_n.size(); i++){
507  val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
508  jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
509  }
510 
511  auto normal = ig.unitOuterNormal(ip.position());
512  auto weight = ip.weight()*geo.integrationElement(ip.position());
513  auto mu = prm.mu(ig,ip.position());
514 
515  auto factor = mu * weight;
516 
517  // compute jump in velocity
518  auto jump = val_u_s - val_u_n;
519 
520  // compute mean in pressure
521  auto mean_p = 0.5*(val_p_s + val_p_n);
522 
523  // compute flux of velocity jacobian
524  Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
525  add_compute_flux(jac_u_s,normal,flux_jac_u);
526  add_compute_flux(jac_u_n,normal,flux_jac_u);
527  flux_jac_u *= 0.5;
528 
529  // loop over test functions, same element
530  for (size_t i=0; i<lfsv_v_s.size(); i++) {
531  //================================================//
532  // diffusion term
533  //================================================//
534  r_s.accumulate(lfsv_v_s, i, -(flux_jac_u * phi_v_s[i]) * factor);
535 
536  //================================================//
537  // (non-)symmetric IP term
538  //================================================//
539  Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
540  add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi);
541  r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
542 
543  //================================================//
544  // standard IP term integral
545  //================================================//
546  r_s.accumulate(lfsv_v_s,i, penalty_factor * (jump*phi_v_s[i]) * factor);
547 
548  //================================================//
549  // pressure-velocity-coupling in momentum equation
550  //================================================//
551  r_s.accumulate(lfsv_v_s,i, mean_p * (phi_v_s[i]*normal) * weight);
552  }
553 
554  // loop over test functions, neighbour element
555  for (size_t i=0; i<lfsv_v_n.size(); i++) {
556  //================================================//
557  // diffusion term
558  //================================================//
559  r_n.accumulate(lfsv_v_n, i, (flux_jac_u * phi_v_n[i]) * factor);
560 
561  //================================================//
562  // (non-)symmetric IP term
563  //================================================//
564  Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
565  add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi);
566  r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
567 
568  //================================================//
569  // standard IP term integral
570  //================================================//
571  r_n.accumulate(lfsv_v_n,i, -penalty_factor * (jump*phi_v_n[i]) * factor);
572 
573  //================================================//
574  // pressure-velocity-coupling in momentum equation
575  //================================================//
576  r_n.accumulate(lfsv_v_n,i, -mean_p * (phi_v_n[i]*normal) * weight);
577  }
578 
579  //================================================//
580  // incompressibility constraint
581  //================================================//
582  for (size_t i=0; i<lfsv_p_s.size(); i++)
583  r_s.accumulate(lfsv_p_s,i, 0.5*phi_p_s[i] * (jump*normal) * incomp_scaling * weight);
584  for (size_t i=0; i<lfsv_p_n.size(); i++)
585  r_n.accumulate(lfsv_p_n,i, 0.5*phi_p_n[i] * (jump*normal) * incomp_scaling * weight);
586 
587  } // end loop quadrature points
588  } // end alpha_skeleton
589 
590  // jacobian of skeleton term, each face is only visited ONCE
591  template<typename IG, typename LFSU, typename X, typename LFSV,
592  typename LocalMatrix>
593  void jacobian_skeleton (const IG& ig,
594  const LFSU& lfsu_s, const X&, const LFSV& lfsv_s,
595  const LFSU& lfsu_n, const X&, const LFSV& lfsv_n,
596  LocalMatrix& mat_ss, LocalMatrix& mat_sn,
597  LocalMatrix& mat_ns, LocalMatrix& mat_nn) const
598  {
599  // dimensions
600  const unsigned int dim = IG::Entity::dimension;
601  const unsigned int dimw = IG::coorddimension;
602 
603  // subspaces
604  using namespace Indices;
605  using LFSV_V = TypeTree::Child<LFSV,_0>;
606  const auto& lfsv_v_s = child(lfsv_s,_0);
607  const auto& lfsu_v_s = child(lfsu_s,_0);
608  const auto& lfsv_v_n = child(lfsv_n,_0);
609  const auto& lfsu_v_n = child(lfsu_n,_0);
610 
611  using LFSV_P = TypeTree::Child<LFSV,_1>;
612  const auto& lfsv_p_s = child(lfsv_s,_1);
613  const auto& lfsu_p_s = child(lfsu_s,_1);
614  const auto& lfsv_p_n = child(lfsv_n,_1);
615  const auto& lfsu_p_n = child(lfsu_n,_1);
616 
617  // domain and range field type
618  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
619  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
621  using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
622  using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
623  using DF = typename BasisSwitch_V::DomainField;
624  using RF = typename BasisSwitch_V::RangeField;
625  using Range_V = typename BasisSwitch_V::Range;
626  using Range_P = typename BasisSwitch_P::Range;
627  using size_type = typename LFSV::Traits::SizeType;
628 
629  // References to inside and outside cells
630  auto const& cell_inside = ig.inside();
631  auto const& cell_outside = ig.outside();
632 
633  // Get geometries
634  auto geo = ig.geometry();
635  auto geo_inside = cell_inside.geometry();
636  auto geo_outside = cell_outside.geometry();
637 
638  // Get geometry of intersection in local coordinates of cell_inside and cell_outside
639  auto geo_in_inside = ig.geometryInInside();
640  auto geo_in_outside = ig.geometryInOutside();
641 
642  // Determine quadrature order
643  const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
644  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
645  const int qorder = 2*v_order + det_jac_order + superintegration_order;
646 
647  auto epsilon = prm.epsilonIPSymmetryFactor();
648  auto incomp_scaling = prm.incompressibilityScaling(current_dt);
649 
650  auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
651 
652  // loop over quadrature points and integrate normal flux
653  for (const auto& ip : quadratureRule(geo,qorder))
654  {
655 
656  // position of quadrature point in local coordinates of element
657  auto local_s = geo_in_inside.global(ip.position());
658  auto local_n = geo_in_outside.global(ip.position());
659 
660  // values of velocity shape functions
661  std::vector<Range_V> phi_v_s(lfsv_v_s.size());
662  std::vector<Range_V> phi_v_n(lfsv_v_n.size());
663  FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
664  FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
665 
666  // values of pressure shape functions
667  std::vector<Range_P> phi_p_s(lfsv_p_s.size());
668  std::vector<Range_P> phi_p_n(lfsv_p_n.size());
669  FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
670  FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
671 
672  // evaluate jacobian of velocity shape functions on reference element
673  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
674  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
675  VectorBasisSwitch_V::jacobian
676  (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
677  VectorBasisSwitch_V::jacobian
678  (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
679 
680  auto normal = ig.unitOuterNormal(ip.position());
681  auto weight = ip.weight()*geo.integrationElement(ip.position());
682  auto mu = prm.mu(ig,ip.position());
683 
684  auto factor = mu * weight;
685 
686  //============================================
687  // loop over test functions, same element
688  //============================================
689  for(size_type i=0; i<lfsv_v_s.size(); i++) {
690 
691  // compute flux
692  Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
693  add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi_i);
694 
695  //============================================
696  // diffusion
697  // (non-)symmetric IP-Term
698  // standard IP integral
699  //============================================
700  for(size_type j=0; j<lfsu_v_s.size(); j++) {
701  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
702  add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
703 
704  mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
705  mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
706  mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, penalty_factor * (phi_v_s[j]*phi_v_s[i]) * factor);
707  }
708 
709  for(size_type j=0; j<lfsu_v_n.size(); j++) {
710  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
711  add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
712 
713  mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
714  mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
715  mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -penalty_factor * (phi_v_n[j]*phi_v_s[i]) * factor);
716  }
717 
718  //============================================
719  // pressure-velocity coupling in momentum equation
720  //============================================
721  for(size_type j=0; j<lfsu_p_s.size(); j++) {
722  mat_ss.accumulate(lfsv_v_s,i,lfsu_p_s,j, 0.5*phi_p_s[j] * (phi_v_s[i]*normal) * weight);
723  }
724 
725  for(size_type j=0; j<lfsu_p_n.size(); j++) {
726  mat_sn.accumulate(lfsv_v_s,i,lfsu_p_n,j, 0.5*phi_p_n[j] * (phi_v_s[i]*normal) * weight);
727  }
728  } // end i (same)
729 
730  //============================================
731  // loop over test functions, neighbour element
732  //============================================
733  for(size_type i=0; i<lfsv_v_n.size(); i++) {
734 
735  // compute flux
736  Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
737  add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi_i);
738 
739  //============================================
740  // diffusion
741  // (non-)symmetric IP-Term
742  // standard IP integral
743  //============================================
744  for(size_type j=0; j<lfsu_v_s.size(); j++) {
745  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
746  add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
747 
748  mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
749  mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
750  mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, -penalty_factor * (phi_v_s[j]*phi_v_n[i]) * factor);
751  }
752 
753  for(size_type j=0; j<lfsu_v_n.size(); j++) {
754  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
755  add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
756 
757  mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
758  mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
759  mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, penalty_factor * (phi_v_n[j]*phi_v_n[i]) * factor);
760  }
761 
762  //============================================
763  // pressure-velocity coupling in momentum equation
764  //============================================
765  for(size_type j=0; j<lfsu_p_s.size(); j++) {
766  mat_ns.accumulate(lfsv_v_n,i,lfsu_p_s,j, -0.5*phi_p_s[j] * (phi_v_n[i]*normal) * weight);
767  }
768 
769  for(size_type j=0; j<lfsu_p_n.size(); j++) {
770  mat_nn.accumulate(lfsv_v_n,i,lfsu_p_n,j, -0.5*phi_p_n[j] * (phi_v_n[i]*normal) * weight);
771  }
772  } // end i (neighbour)
773 
774  //================================================//
775  // \int <q> [u] n
776  //================================================//
777  for(size_type i=0; i<lfsv_p_s.size(); i++) {
778  for(size_type j=0; j<lfsu_v_s.size(); j++)
779  mat_ss.accumulate(lfsv_p_s,i,lfsu_v_s,j, 0.5*phi_p_s[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
780 
781  for(size_type j=0; j<lfsu_v_n.size(); j++)
782  mat_sn.accumulate(lfsv_p_s,i,lfsu_v_n,j, -0.5*phi_p_s[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
783  }
784 
785  for(size_type i=0; i<lfsv_p_n.size(); i++) {
786  for(size_type j=0; j<lfsu_v_s.size(); j++)
787  mat_ns.accumulate(lfsv_p_n,i,lfsu_v_s,j, 0.5*phi_p_n[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
788 
789  for(size_type j=0; j<lfsu_v_n.size(); j++)
790  mat_nn.accumulate(lfsv_p_n,i,lfsu_v_n,j, -0.5*phi_p_n[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
791  }
792 
793  } // end loop quadrature points
794  } // end jacobian_skeleton
795 
796  // boundary term
797  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
798  void alpha_boundary (const IG& ig,
799  const LFSU& lfsu, const X& x, const LFSV& lfsv,
800  R& r) const
801  {
802  // dimensions
803  const unsigned int dim = IG::Entity::dimension;
804  const unsigned int dimw = IG::coorddimension;
805 
806  // subspaces
807  using namespace Indices;
808  using LFSV_V = TypeTree::Child<LFSV,_0>;
809  const auto& lfsv_v = child(lfsv,_0);
810  const auto& lfsu_v = child(lfsu,_0);
811 
812  using LFSV_P = TypeTree::Child<LFSV,_1>;
813  const auto& lfsv_p = child(lfsv,_1);
814  const auto& lfsu_p = child(lfsu,_1);
815 
816  // domain and range field type
817  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
818  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
820  using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
821  using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
822  using DF = typename BasisSwitch_V::DomainField;
823  using RF = typename BasisSwitch_V::RangeField;
824  using Range_V = typename BasisSwitch_V::Range;
825  using Range_P = typename BasisSwitch_P::Range;
826  using size_type = typename LFSV::Traits::SizeType;
827 
828  // References to inside cell
829  const auto& cell_inside = ig.inside();
830 
831  // Get geometries
832  auto geo = ig.geometry();
833  auto geo_inside = cell_inside.geometry();
834 
835  // Get geometry of intersection in local coordinates of cell_inside
836  auto geo_in_inside = ig.geometryInInside();
837 
838  // Determine quadrature order
839  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
840  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
841  const int qorder = 2*v_order + det_jac_order + superintegration_order;
842 
843  auto epsilon = prm.epsilonIPSymmetryFactor();
844  auto incomp_scaling = prm.incompressibilityScaling(current_dt);
845 
846  auto penalty_factor = prm.getFaceIP(geo,geo_inside);
847 
848  // loop over quadrature points and integrate normal flux
849  for (const auto& ip : quadratureRule(geo,qorder))
850  {
851  // position of quadrature point in local coordinates of element
852  auto local = geo_in_inside.global(ip.position());
853 
854  // values of velocity shape functions
855  std::vector<Range_V> phi_v(lfsv_v.size());
856  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
857 
858  // values of pressure shape functions
859  std::vector<Range_P> phi_p(lfsv_p.size());
860  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
861 
862  // evaluate jacobian of basis functions on reference element
863  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
864  VectorBasisSwitch_V::jacobian
865  (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
866 
867  // compute pressure value
868  Range_P val_p(0.0);
869  for (size_type i=0; i<lfsu_p.size(); i++)
870  val_p.axpy(x(lfsu_p,i),phi_p[i]);
871 
872  // compute u and velocity jacobian
873  Range_V val_u(0.0);
874  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
875  for (size_type i=0; i<lfsu_v.size(); i++){
876  val_u.axpy(x(lfsu_v,i),phi_v[i]);
877  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
878  }
879 
880  auto normal = ig.unitOuterNormal(ip.position());
881  auto weight = ip.weight()*geo.integrationElement(ip.position());
882  auto mu = prm.mu(ig,ip.position());
883 
884  // evaluate boundary condition type
885  auto bctype(prm.bctype(ig,ip.position()));
886 
887  if (bctype == BC::VelocityDirichlet) {
888  // compute jump relative to Dirichlet value
889  auto u0(prm.g(cell_inside,local));
890  auto jump = val_u - u0;
891 
892  // compute flux of velocity jacobian
893  Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
894  add_compute_flux(jac_u,normal,flux_jac_u);
895 
896  for (size_t i=0; i<lfsv_v.size(); i++) {
897  //================================================//
898  // diffusion term
899  //================================================//
900  r.accumulate(lfsv_v,i, -mu * (flux_jac_u * phi_v[i]) * weight);
901 
902  //================================================//
903  // (non-)symmetric IP term
904  //================================================//
905  Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
906  add_compute_flux(jac_phi_v[i],normal,flux_jac_phi);
907  r.accumulate(lfsv_v,i, mu * epsilon * (flux_jac_phi*jump) * weight);
908 
909  //================================================//
910  // standard IP term integral
911  //================================================//
912  r.accumulate(lfsv_v,i, mu * (jump*phi_v[i]) * penalty_factor * weight);
913 
914  //================================================//
915  // pressure-velocity-coupling in momentum equation
916  //================================================//
917  r.accumulate(lfsv_v,i, val_p * (phi_v[i]*normal) * weight);
918  } // end i
919 
920  //================================================//
921  // incompressibility constraint
922  //================================================//
923  for(size_type i=0; i<lfsv_p.size(); i++) {
924  r.accumulate(lfsv_p,i, phi_p[i] * (jump*normal) * incomp_scaling * weight);
925  }
926  } // Velocity Dirichlet
927 
928  if (bctype == BC::StressNeumann) {
929  auto stress(prm.j(ig,ip.position(),normal));
930 
931  for(size_type i=0; i<lfsv_v.size(); i++) {
932  r.accumulate(lfsv_v,i, (stress*phi_v[i]) * weight);
933  }
934  } // Pressure Dirichlet
935 
936  } // end loop quadrature points
937  } // end alpha_boundary
938 
939  // jacobian of boundary term
940  template<typename IG, typename LFSU, typename X, typename LFSV,
941  typename LocalMatrix>
942  void jacobian_boundary (const IG& ig,
943  const LFSU& lfsu, const X& x, const LFSV& lfsv,
944  LocalMatrix& mat) const
945  {
946  // dimensions
947  const unsigned int dim = IG::Entity::dimension;
948  const unsigned int dimw = IG::coorddimension;
949 
950  // subspaces
951  using namespace Indices;
952  using LFSV_V = TypeTree::Child<LFSV,_0>;
953  const auto& lfsv_v = child(lfsv,_0);
954  const auto& lfsu_v = child(lfsu,_0);
955 
956  using LFSV_P = TypeTree::Child<LFSV,_1>;
957  const auto& lfsv_p = child(lfsv,_1);
958  const auto& lfsu_p = child(lfsu,_1);
959 
960  // domain and range field type
961  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
962  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
964  using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
965  using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
966  using DF = typename BasisSwitch_V::DomainField;
967  using RF = typename BasisSwitch_V::RangeField;
968  using Range_V = typename BasisSwitch_V::Range;
969  using Range_P = typename BasisSwitch_P::Range;
970  using size_type = typename LFSV::Traits::SizeType;
971 
972  // References to inside cell
973  const auto& cell_inside = ig.inside();
974 
975  // Get geometries
976  auto geo = ig.geometry();
977  auto geo_inside = cell_inside.geometry();
978 
979  // Get geometry of intersection in local coordinates of cell_inside
980  auto geo_in_inside = ig.geometryInInside();
981 
982  // Determine quadrature order
983  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
984  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
985  const int qorder = 2*v_order + det_jac_order + superintegration_order;
986 
987  auto epsilon = prm.epsilonIPSymmetryFactor();
988  auto incomp_scaling = prm.incompressibilityScaling(current_dt);
989 
990  auto penalty_factor = prm.getFaceIP(geo,geo_inside);
991 
992  // loop over quadrature points and integrate normal flux
993  for (const auto& ip : quadratureRule(geo,qorder))
994  {
995  // position of quadrature point in local coordinates of element
996  auto local = geo_in_inside.global(ip.position());
997 
998  // values of velocity shape functions
999  std::vector<Range_V> phi_v(lfsv_v.size());
1000  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1001 
1002  // values of pressure shape functions
1003  std::vector<Range_P> phi_p(lfsv_p.size());
1004  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1005 
1006  // evaluate jacobian of basis functions on reference element
1007  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
1008  VectorBasisSwitch_V::jacobian
1009  (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
1010 
1011  auto normal = ig.unitOuterNormal(ip.position());
1012  auto weight = ip.weight()*geo.integrationElement(ip.position());
1013  auto mu = prm.mu(ig,ip.position());
1014 
1015  // evaluate boundary condition type
1016  auto bctype(prm.bctype(ig,ip.position()));
1017 
1018  if (bctype == BC::VelocityDirichlet) {
1019 
1020  for(size_type i=0; i<lfsv_v.size(); i++) {
1021  // compute flux
1022  Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
1023  add_compute_flux(jac_phi_v[i],normal,flux_jac_phi_i);
1024 
1025  for(size_type j=0; j<lfsu_v.size(); j++) {
1026  //================================================//
1027  // diffusion term
1028  // (non-)symmetric IP term
1029  //================================================//
1030  Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
1031  add_compute_flux(jac_phi_v[j],normal,flux_jac_phi_j);
1032 
1033  mat.accumulate(lfsv_v,i,lfsu_v,j, -mu * (flux_jac_phi_j*phi_v[i]) * weight);
1034  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * epsilon * (flux_jac_phi_i*phi_v[j]) *weight);
1035 
1036  //================================================//
1037  // standard IP term integral
1038  //================================================//
1039  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (phi_v[j]*phi_v[i]) * penalty_factor * weight);
1040  }
1041 
1042  //================================================//
1043  // pressure-velocity-coupling in momentum equation
1044  //================================================//
1045  for(size_type j=0; j<lfsu_p.size(); j++) {
1046  mat.accumulate(lfsv_v,i,lfsu_p,j, phi_p[j] * (phi_v[i]*normal) * weight);
1047  }
1048  } // end i
1049 
1050  //================================================//
1051  // incompressibility constraint
1052  //================================================//
1053  for(size_type i=0; i<lfsv_p.size(); i++) {
1054  for(size_type j=0; j<lfsu_v.size(); j++) {
1055  mat.accumulate(lfsv_p,i,lfsu_v,j, phi_p[i] * (phi_v[j]*normal) * incomp_scaling * weight);
1056  }
1057  }
1058 
1059  } // Velocity Dirichlet
1060 
1061  } // end loop quadrature points
1062  } // end jacobian_boundary
1063 
1064  // volume integral depending only on test functions,
1065  // contains f on the right hand side
1066  template<typename EG, typename LFSV, typename R>
1067  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
1068  {
1069  const unsigned int dim = EG::Geometry::mydimension;
1070 
1071  // subspaces
1072  using namespace Indices;
1073  using LFSV_V = TypeTree::Child<LFSV,_0>;
1074  const auto& lfsv_v = child(lfsv,_0);
1075 
1076  // domain and range field type
1077  using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1078  using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1079  using Range_V = typename BasisSwitch_V::Range;
1080  using size_type = typename LFSV::Traits::SizeType;
1081 
1082  // Get cell
1083  const auto& cell = eg.entity();
1084 
1085  // Get geometries
1086  auto geo = eg.geometry();
1087 
1088  // Determine quadrature order
1089  const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1090  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
1091  const int qorder = 2*v_order + det_jac_order + superintegration_order;
1092 
1093  // loop over quadrature points
1094  for (const auto& ip : quadratureRule(geo,qorder))
1095  {
1096  auto local = ip.position();
1097  //const Dune::FieldVector<DF,dimw> global = eg.geometry().global(local);
1098 
1099  // values of velocity shape functions
1100  std::vector<Range_V> phi_v(lfsv_v.size());
1101  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1102 
1103  auto weight = ip.weight() * geo.integrationElement(ip.position());
1104 
1105  // evaluate source term
1106  auto fval(prm.f(cell,local));
1107 
1108  //================================================//
1109  // \int (f*v)
1110  //================================================//
1111  for(size_type i=0; i<lfsv_v.size(); i++)
1112  r.accumulate(lfsv_v,i, -(fval*phi_v[i]) * weight);
1113 
1114  } // end loop quadrature points
1115  } // end lambda_volume
1116 
1117  private :
1118 
1119  template<class M, class RF>
1120  static void contraction(const M & a, const M & b, RF & v)
1121  {
1122  v = 0;
1123  const int an = a.N();
1124  const int am = a.M();
1125  for(int r=0; r<an; ++r)
1126  for(int c=0; c<am; ++c)
1127  v += a[r][c] * b[r][c];
1128  }
1129 
1130  template<class DU, class R>
1131  static void add_compute_flux(const DU & du, const R & n, R & result)
1132  {
1133  const int N = du.N();
1134  const int M = du.M();
1135  for(int r=0; r<N; ++r)
1136  for(int c=0; c<M; ++c)
1137  result[r] += du[r][c] * n[c];
1138  }
1139 
1140  PRM& prm;
1141  const int superintegration_order;
1142  Real current_dt;
1143  }; // end class DGNavierStokesVelVecFEM
1144 
1145  } // end namespace PDELab
1146 } // end namespace Dune
1147 #endif // DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
Dune::PDELab::DGNavierStokesVelVecFEM::doAlphaBoundary
@ doAlphaBoundary
Definition: dgnavierstokesvelvecfem.hh:127
flags.hh
Dune::PDELab::VectorBasisInterfaceSwitch< Basis, typename std::enable_if< Std::to_true_type< std::integral_constant< std::size_t, Basis::Traits::dimDomain > >::value >::type >::DomainLocal
typename Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:65
Dune::PDELab::StokesBoundaryCondition
Definition: stokesparameter.hh:32
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
dim
static const int dim
Definition: adaptivity.hh:84
Dune::PDELab::DGNavierStokesVelVecFEM::jacobian_volume
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokesvelvecfem.hh:283
dgnavierstokesparameter.hh
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:18
Dune::PDELab::FullSkeletonPattern
sparsity pattern generator
Definition: pattern.hh:29
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::PDELab::DGNavierStokesVelVecFEM::doAlphaSkeleton
@ doAlphaSkeleton
Definition: dgnavierstokesvelvecfem.hh:126
Dune::PDELab::DGNavierStokesVelVecFEM::doPatternVolume
@ doPatternVolume
Definition: dgnavierstokesvelvecfem.hh:118
Dune::PDELab::DGNavierStokesVelVecFEM::DGNavierStokesVelVecFEM
DGNavierStokesVelVecFEM(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokesvelvecfem.hh:142
navierstokesmass.hh
defaultimp.hh
idefault.hh
Dune::PDELab::VectorBasisInterfaceSwitch
Definition: dgnavierstokesvelvecfem.hh:31
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:13
Dune::PDELab::InstationaryLocalOperatorDefaultMethods::RealType
R RealType
Definition: idefault.hh:92
Dune::PDELab::VectorBasisInterfaceSwitch::jacobian
static void jacobian(const Basis &basis, const Geometry &geometry, const DomainLocal &xl, std::vector< FieldMatrix< RangeField, dimRange, Geometry::coorddimension > > &jac)
Compute global jacobian matrix for vector valued bases.
Definition: dgnavierstokesvelvecfem.hh:41
Dune::PDELab::DGNavierStokesVelVecFEM::doPatternSkeleton
@ doPatternSkeleton
Definition: dgnavierstokesvelvecfem.hh:119
Dune::PDELab::DGNavierStokesVelVecFEM::alpha_skeleton
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: dgnavierstokesvelvecfem.hh:403
Dune::PDELab::DGNavierStokesVelVecFEM::alpha_volume
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:162
Dune::PDELab::DGNavierStokesVelVecFEM::jacobian_skeleton
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: dgnavierstokesvelvecfem.hh:593
Dune::PDELab::DGNavierStokesVelVecFEM::alpha_boundary
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:798
ig
const IG & ig
Definition: constraints.hh:148
Dune::PDELab::DGNavierStokesVelVecFEM::doSkeletonTwoSided
@ doSkeletonTwoSided
Definition: dgnavierstokesvelvecfem.hh:122
Dune::PDELab::VectorBasisInterfaceSwitch::RangeField
typename Basis::Traits::RangeField RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:35
Dune::PDELab::DGNavierStokesVelVecFEM::setTime
void setTime(Real t)
Definition: dgnavierstokesvelvecfem.hh:154
Dune::PDELab::DGNavierStokesVelVecFEM::preStep
void preStep(Real, Real dt, int)
Definition: dgnavierstokesvelvecfem.hh:148
Dune::PDELab::DGNavierStokesVelVecFEM
A local operator for solving the Navier-Stokes equations using a DG discretization and a vector-value...
Definition: dgnavierstokesvelvecfem.hh:101
Dune::PDELab::DGNavierStokesVelVecFEM::doLambdaVolume
@ doLambdaVolume
Definition: dgnavierstokesvelvecfem.hh:128
Dune::PDELab::DGNavierStokesVelVecFEM::jacobian_boundary
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix &mat) const
Definition: dgnavierstokesvelvecfem.hh:942
Dune::PDELab::quadratureRule
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
Dune::PDELab::LocalMatrix
A dense matrix for storing data associated with the degrees of freedom of a pair of LocalFunctionSpac...
Definition: localmatrix.hh:184
Dune::PDELab::StokesBoundaryCondition::VelocityDirichlet
@ VelocityDirichlet
Definition: stokesparameter.hh:35
pattern.hh
Dune::PDELab::VectorBasisInterfaceSwitch::DomainLocal
typename Basis::Traits::DomainLocal DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:33
Dune::PDELab::VectorBasisInterfaceSwitch< Basis, typename std::enable_if< Std::to_true_type< std::integral_constant< std::size_t, Basis::Traits::dimDomain > >::value >::type >::jacobian
static void jacobian(const Basis &basis, const Geometry &geometry, const DomainLocal &xl, std::vector< FieldMatrix< RangeField, dimRange, Geometry::coorddimension > > &jac)
Compute global jacobian matrix for vector valued bases.
Definition: dgnavierstokesvelvecfem.hh:73
Dune::PDELab::DGNavierStokesVelVecFEM::doAlphaVolume
@ doAlphaVolume
Definition: dgnavierstokesvelvecfem.hh:125
Dune::PDELab::VectorBasisInterfaceSwitch< Basis, typename std::enable_if< Std::to_true_type< std::integral_constant< std::size_t, Basis::Traits::dimDomain > >::value >::type >::RangeField
typename Basis::Traits::RangeFieldType RangeField
export field type of the values
Definition: dgnavierstokesvelvecfem.hh:67
Dune::PDELab::DGNavierStokesVelVecFEM::lambda_volume
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: dgnavierstokesvelvecfem.hh:1067
Dune::PDELab::InstationaryLocalOperatorDefaultMethods::setTime
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Dune::PDELab::StokesBoundaryCondition::StressNeumann
@ StressNeumann
Definition: stokesparameter.hh:36
Dune::PDELab::VectorBasisInterfaceSwitch::dimRange
static const std::size_t dimRange
export dimension of the values
Definition: dgnavierstokesvelvecfem.hh:37