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