dune-pdelab  2.0.0
stokesdg_vecfem.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_STOKESDG_VECFEM_HH
3 #define DUNE_PDELAB_STOKESDG_VECFEM_HH
4 
5 #include <dune/common/exceptions.hh>
6 #include <dune/common/fvector.hh>
7 #include <dune/common/static_assert.hh>
8 #include <dune/common/parametertree.hh>
9 
10 #include <dune/geometry/quadraturerules.hh>
11 
12 #include <dune/localfunctions/common/interfaceswitch.hh>
14 
15 #include "defaultimp.hh"
16 #include "pattern.hh"
17 #include "flags.hh"
18 #include "stokesdg.hh"
19 
20 #ifndef VBLOCK
21 #define VBLOCK 0
22 #endif
23 #define PBLOCK (- VBLOCK + 1)
24 
25 namespace Dune {
26  namespace PDELab {
27 
28  template<class Basis, class Dummy = void>
31  typedef typename Basis::Traits::DomainLocal DomainLocal;
33  typedef typename Basis::Traits::RangeField RangeField;
35  static const std::size_t dimRange = Basis::Traits::dimRange;
36 
38  template<typename Geometry>
39  static void jacobian(const Basis& basis, const Geometry& geometry,
40  const DomainLocal& xl,
41  std::vector<FieldMatrix<RangeField, dimRange,
42  Geometry::coorddimension> >& jac)
43  {
44  jac.resize(basis.size());
45  basis.evaluateJacobian(xl, jac);
46  }
47  };
48 
50  template<class Basis>
52  Basis, typename enable_if<Basis::Traits::dimDomain>::type
53  >
54  {
56  typedef typename Basis::Traits::DomainType DomainLocal;
58  typedef typename Basis::Traits::RangeFieldType RangeField;
60  static const std::size_t dimRange = Basis::Traits::dimRange;
61 
63  template<typename Geometry>
64  static void jacobian(const Basis& basis, const Geometry& geometry,
65  const DomainLocal& xl,
66  std::vector<FieldMatrix<RangeField, dimRange,
67  Geometry::coorddimension> >& jac)
68  {
69  jac.resize(basis.size());
70 
71  std::vector<FieldMatrix<
72  RangeField, dimRange, Geometry::coorddimension> > ljac(basis.size());
73  basis.evaluateJacobian(xl, ljac);
74 
75  const typename Geometry::Jacobian& geojac =
76  geometry.jacobianInverseTransposed(xl);
77 
78  for(std::size_t i = 0; i < basis.size(); ++i)
79  for(std::size_t row=0; row < dimRange; ++row)
80  geojac.mv(ljac[i][row], jac[i][row]);
81  }
82  };
83 
94  template<typename F, typename B, typename V, typename P,
95  typename IP = DefaultInteriorPenalty<typename V::Traits::RangeFieldType> >
99  public Dune::PDELab::NumericalJacobianApplyVolume<StokesDGVectorFEM<F,B,V,P,IP> >,
100  public Dune::PDELab::NumericalJacobianVolume<StokesDGVectorFEM<F,B,V,P,IP> >,
101  public Dune::PDELab::NumericalJacobianApplySkeleton<StokesDGVectorFEM<F,B,V,P,IP> >,
102  public Dune::PDELab::NumericalJacobianSkeleton<StokesDGVectorFEM<F,B,V,P,IP> >,
103  public Dune::PDELab::NumericalJacobianApplyBoundary<StokesDGVectorFEM<F,B,V,P,IP> >,
104  public Dune::PDELab::NumericalJacobianBoundary<StokesDGVectorFEM<F,B,V,P,IP> >,
106  {
107  typedef StokesBoundaryCondition BC;
108  typedef typename V::Traits::RangeFieldType RF;
109 
110  public:
112 
113  // pattern assembly flags
114  enum { doPatternVolume = true };
115  enum { doPatternSkeleton = true };
116 
117  // residual assembly flags
118  enum { doAlphaVolume = true };
119  enum { doAlphaSkeleton = true };
120  enum { doAlphaBoundary = true };
121 
122  StokesDGVectorFEM (const Dune::ParameterTree & configuration,const IP & ip_factor_, const RF mu_,
123  const F & _f, const B & _b, const V & _v, const P & _p, int _qorder=4) :
124  f(_f), b(_b), v(_v), p(_p), qorder(_qorder), mu(mu_), ip_factor(ip_factor_)
125  {
126  epsilon = configuration.get<int>("epsilon");
127  }
128 
129  // volume integral depending on test and ansatz functions
130  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
131  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
132  {
133  // dimensions
134  static const unsigned int dim = EG::Geometry::dimension;
135  static const unsigned int dimw = EG::Geometry::dimensionworld;
136 
137  // subspaces
138  dune_static_assert
139  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for StokesDGVectorFEM");
140 
141  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
142  const LFSV_V& lfsv_v = lfsv.template child<VBLOCK>();
143  const LFSV_V& lfsu_v = lfsu.template child<VBLOCK>();
144 
145  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
146  const LFSV_P& lfsv_p = lfsv.template child<PBLOCK>();
147  const LFSV_P& lfsu_p = lfsu.template child<PBLOCK>();
148 
149  // domain and range field type
150  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
151  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
153  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
154  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
155  typedef typename BasisSwitch_V::DomainField DF;
156  typedef typename BasisSwitch_V::RangeField RF;
157  typedef typename BasisSwitch_V::Range Range_V;
158  typedef typename BasisSwitch_P::Range Range_P;
159  typedef typename LFSV::Traits::SizeType size_type;
160 
161  // select quadrature rule
162  Dune::GeometryType gt = eg.geometry().type();
163  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
164 
165  // loop over quadrature points
166  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
167  {
168  // Compute weight and coordinates
169  const Dune::FieldVector<DF,dim> local = it->position();
170  const Dune::FieldVector<DF,dimw> global = eg.geometry().global(local);
171  const RF weight = it->weight() * eg.geometry().integrationElement(it->position());
172 
173  // values of velocity shape functions
174  std::vector<Range_V> phi_v(lfsv_v.size());
175  FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
176 
177  // values of velocity shape functions
178  std::vector<Range_P> phi_p(lfsv_p.size());
179  FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
180 
181  // evaluate jacobian of basis functions on reference element
182  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
183  VectorBasisSwitch_V::jacobian
184  (FESwitch_V::basis(lfsv_v.finiteElement()), eg.geometry(), it->position(), jac_phi_v);
185 
186  // compute divergence of test functions
187  std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
188  for (size_type i=0; i<lfsv_v.size(); i++)
189  for (size_type d=0; d<dim; d++)
190  div_phi_v[i] += jac_phi_v[i][d][d];
191 
192  // compute velocity value, jacobian, and divergence
193  Range_V val_u(0.0);
194  Dune::FieldMatrix<RF,dim,dim> jac_u(0.0);
195  RF div_u(0.0);
196  for (size_type i=0; i<lfsu_v.size(); i++){
197  val_u.axpy(x(lfsu_v,i),phi_v[i]);
198  jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
199  div_u += x(lfsu_v,i) * div_phi_v[i];
200  }
201 
202  // compute pressure value
203  Range_P val_p(0.0);
204  for (size_type i=0; i<lfsu_p.size(); i++)
205  val_p.axpy(x(lfsu_p,i),phi_p[i]);
206 
207  // evaluate source term
208  typename F::Traits::RangeType fval;
209  f.evaluateGlobal(global,fval);
210 
211  {// Integrate (f*v)
212  const RF factor = weight;
213  for (size_type i=0; i<lfsv_v.size(); i++)
214  r.accumulate(lfsv_v, i, (fval * phi_v[i]) * factor );
215  }
216 
217  {// Integrate (mu * d_i u_j d_i v_j)
218  const RF factor = mu * weight;
219  for (size_type i=0; i<lfsv_v.size(); i++){
220  RF dvdu(0); contraction(jac_u,jac_phi_v[i],dvdu);
221  r.accumulate(lfsv_v, i, dvdu * factor);
222  }
223  }
224 
225  {// Integrate - p div v
226  for (size_type i=0; i<lfsv_v.size(); i++)
227  r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
228  }
229 
230  {// Integrate - q div u
231  for (size_type i=0; i<lfsv_p.size(); i++)
232  r.accumulate(lfsv_p, i, - div_u * phi_p[i] * weight);
233  }
234 
235  }
236  }
237 
238  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
239  void alpha_skeleton (const IG& ig,
240  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
241  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
242  R& r_s, R& r_n) const
243  {
244  // dimensions
245  static const unsigned int dim = IG::Geometry::dimension;
246  static const unsigned int dimw = IG::Geometry::dimensionworld;
247 
248  // subspaces
249  dune_static_assert
250  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for StokesDGVectorFEM");
251 
252  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
253  const LFSV_V& lfsv_v_s = lfsv_s.template child<VBLOCK>();
254  const LFSV_V& lfsu_v_s = lfsu_s.template child<VBLOCK>();
255  const LFSV_V& lfsv_v_n = lfsv_n.template child<VBLOCK>();
256  const LFSV_V& lfsu_v_n = lfsu_n.template child<VBLOCK>();
257 
258  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
259  const LFSV_P& lfsv_p_s = lfsv_s.template child<PBLOCK>();
260  const LFSV_P& lfsu_p_s = lfsu_s.template child<PBLOCK>();
261  const LFSV_P& lfsv_p_n = lfsv_n.template child<PBLOCK>();
262  const LFSV_P& lfsu_p_n = lfsu_n.template child<PBLOCK>();
263 
264  // domain and range field type
265  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
266  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
268  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
269  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
270  typedef typename BasisSwitch_V::DomainField DF;
271  typedef typename BasisSwitch_V::RangeField RF;
272  typedef typename BasisSwitch_V::Range Range_V;
273  typedef typename BasisSwitch_P::Range Range_P;
274  typedef typename LFSV::Traits::SizeType size_type;
275 
276  // select quadrature rule
277  Dune::GeometryType gt = ig.geometry().type();
278  const Dune::QuadratureRule<DF,dim-1>& rule
279  = Dune::QuadratureRules<DF,dim-1>::rule(gt,qorder);
280 
281  const typename IG::EntityPointer self = ig.inside();
282  const typename IG::EntityPointer neighbor = ig.outside();
283 
284  // loop over quadrature points
285  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it){
286 
287  // position of quadrature point in adjacent elements
288  const Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
289  const Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
290  const Dune::FieldVector<DF,dim> global = ig.geometry().global(it->position());
291  const RF weight = it->weight() * ig.geometry().integrationElement(it->position());
292 
293  // values of velocity shape functions
294  std::vector<Range_V> phi_v_s(lfsv_v_s.size());
295  std::vector<Range_V> phi_v_n(lfsv_v_n.size());
296  FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
297  FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
298 
299  // values of velocity shape functions
300  std::vector<Range_P> phi_p_s(lfsv_p_s.size());
301  std::vector<Range_P> phi_p_n(lfsv_p_n.size());
302  FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
303  FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
304 
305  // evaluate jacobian of basis functions on reference element
306  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
307  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
308  VectorBasisSwitch_V::jacobian
309  (FESwitch_V::basis(lfsv_v_s.finiteElement()), ig.inside()->geometry(), local_s, jac_phi_v_s);
310  VectorBasisSwitch_V::jacobian
311  (FESwitch_V::basis(lfsv_v_n.finiteElement()), ig.outside()->geometry(), local_n, jac_phi_v_n);
312 
313  // compute divergence of test functions
314  std::vector<RF> div_phi_v_s(lfsv_v_s.size(),0.0);
315  std::vector<RF> div_phi_v_n(lfsv_v_s.size(),0.0);
316  for (size_type d=0; d<dim; d++){
317  for (size_type i=0; i<lfsv_v_s.size(); i++)
318  div_phi_v_s[i] += jac_phi_v_s[i][d][d];
319  for (size_type i=0; i<lfsv_v_n.size(); i++)
320  div_phi_v_n[i] += jac_phi_v_n[i][d][d];
321  }
322 
323  // compute velocity value, jacobian, and divergence
324  Range_V val_u_s(0.0);
325  Range_V val_u_n(0.0);
326  Dune::FieldMatrix<RF,dim,dim> jac_u_s(0.0);
327  Dune::FieldMatrix<RF,dim,dim> jac_u_n(0.0);
328  RF div_u_s(0.0);
329  RF div_u_n(0.0);
330  for (size_type i=0; i<lfsu_v_s.size(); i++){
331  val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
332  jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
333  div_u_s += x_s(lfsu_v_s,i) * div_phi_v_s[i];
334  }
335  for (size_type i=0; i<lfsu_v_n.size(); i++){
336  val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
337  jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
338  div_u_n += x_n(lfsu_v_n,i) * div_phi_v_n[i];
339  }
340 
341  // compute pressure value
342  Range_P val_p_s(0.0);
343  Range_P val_p_n(0.0);
344  for (size_type i=0; i<lfsu_p_s.size(); i++)
345  val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
346  for (size_type i=0; i<lfsu_p_n.size(); i++)
347  val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
348 
349  // face normal vector
350  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
351 
352  // compute jump in solution
353  const Dune::FieldVector<DF,dimw> jump = val_u_s - val_u_n;
354 
355  {// update residual (flux term)
356 
357  // compute flux
358  Dune::FieldVector<DF,dimw> flux(0.0);
359  add_compute_flux(jac_u_s,normal,flux);
360  add_compute_flux(jac_u_n,normal,flux);
361  flux *= 0.5;
362 
363  const RF factor = weight * mu;
364  for (size_t i=0; i<lfsv_v_s.size(); i++)
365  r_s.accumulate(lfsv_v_s, i, -(flux * phi_v_s[i]) * factor);
366  for (size_t i=0; i<lfsv_v_n.size(); i++)
367  r_n.accumulate(lfsv_v_n, i, (flux * phi_v_n[i]) * factor);
368  }
369 
370  {// update residual (symmetry term)
371  const RF factor = weight * mu;
372 
373  for (size_t i=0; i<lfsv_v_s.size(); i++){
374  Dune::FieldVector<DF,dimw> flux(0.0);
375  add_compute_flux(jac_phi_v_s[i],normal,flux);
376  r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux * jump) * factor);
377  }
378  for (size_t i=0; i<lfsv_v_n.size(); i++){
379  Dune::FieldVector<DF,dimw> flux(0.0);
380  add_compute_flux(jac_phi_v_n[i],normal,flux);
381  r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux * jump) * factor);
382  }
383  }
384 
385  {// interior penalty
386  const RF factor = weight;
387  const RF gamma = ip_factor.getFaceIP(ig);
388  for (size_t i=0; i<lfsv_v_s.size(); i++)
389  r_s.accumulate(lfsv_v_s,i, gamma * (jump * phi_v_s[i]) * factor);
390  for (size_t i=0; i<lfsv_v_n.size(); i++)
391  r_n.accumulate(lfsv_v_n,i, - gamma * (jump * phi_v_n[i]) * factor);
392  }
393 
394  {// pressure and incompressibility
395  const RF factor = weight;
396  const RF val_p = 0.5 * (val_p_s + val_p_n);
397  for (size_t i=0; i<lfsv_v_s.size(); i++)
398  r_s.accumulate(lfsv_v_s, i, val_p * (phi_v_s[i] * normal) * factor);
399  for (size_t i=0; i<lfsv_v_n.size(); i++)
400  r_n.accumulate(lfsv_v_n, i, - val_p * (phi_v_n[i] * normal) * factor);
401 
402  for (size_t i=0; i<lfsv_p_s.size(); i++)
403  r_s.accumulate(lfsv_p_s, i, 0.5 * phi_p_s[i] * (jump*normal) * factor);
404  for (size_t i=0; i<lfsv_p_n.size(); i++)
405  r_n.accumulate(lfsv_p_n, i, 0.5 * phi_p_n[i] * (jump*normal) * factor);
406  }
407 
408  } // it - quadrature
409 
410  }
411 
412  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
413  void alpha_boundary (const IG& ig,
414  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
415  R& r_s) const
416  {
417  // dimensions
418  static const unsigned int dim = IG::Geometry::dimension;
419  static const unsigned int dimw = IG::Geometry::dimensionworld;
420 
421  // subspaces
422  dune_static_assert
423  ((LFSV::CHILDREN == 2), "You seem to use the wrong function space for StokesDGVectorFEM");
424 
425  typedef typename LFSV::template Child<VBLOCK>::Type LFSV_V;
426  const LFSV_V& lfsv_v_s = lfsv_s.template child<VBLOCK>();
427  const LFSV_V& lfsu_v_s = lfsu_s.template child<VBLOCK>();
428 
429  typedef typename LFSV::template Child<PBLOCK>::Type LFSV_P;
430  const LFSV_P& lfsv_p_s = lfsv_s.template child<PBLOCK>();
431  const LFSV_P& lfsu_p_s = lfsu_s.template child<PBLOCK>();
432 
433  // domain and range field type
434  typedef FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType > FESwitch_V;
435  typedef BasisInterfaceSwitch<typename FESwitch_V::Basis > BasisSwitch_V;
437  typedef FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType > FESwitch_P;
438  typedef BasisInterfaceSwitch<typename FESwitch_P::Basis > BasisSwitch_P;
439  typedef typename BasisSwitch_V::DomainField DF;
440  typedef typename BasisSwitch_V::RangeField RF;
441  typedef typename BasisSwitch_V::Range Range_V;
442  typedef typename BasisSwitch_P::Range Range_P;
443  typedef typename LFSV::Traits::SizeType size_type;
444 
445  // select quadrature rule
446  Dune::GeometryType gt = ig.geometry().type();
447  const Dune::QuadratureRule<DF,dim-1>& rule
448  = Dune::QuadratureRules<DF,dim-1>::rule(gt,qorder);
449 
450  const typename IG::EntityPointer self = ig.inside();
451 
452  // loop over quadrature points
453  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it){
454 
455  // position of quadrature point in adjacent elements
456  const Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
457  const Dune::FieldVector<DF,dim> global = ig.geometry().global(it->position());
458  const RF weight = it->weight() * ig.geometry().integrationElement(it->position());
459 
460  // evaluate boundary condition type
461  typename B::Traits::RangeType bctype;
462  b.evaluate(ig,it->position(),bctype);
463 
464  // values of velocity shape functions
465  std::vector<Range_V> phi_v_s(lfsv_v_s.size());
466  FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
467 
468  // values of velocity shape functions
469  std::vector<Range_P> phi_p_s(lfsv_p_s.size());
470  FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
471 
472  // evaluate jacobian of basis functions on reference element
473  std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
474  VectorBasisSwitch_V::jacobian
475  (FESwitch_V::basis(lfsv_v_s.finiteElement()), ig.inside()->geometry(), local_s, jac_phi_v_s);
476 
477  // compute divergence of test functions
478  std::vector<RF> div_phi_v_s(lfsv_v_s.size(),0.0);
479 
480  for (size_type d=0; d<dim; d++){
481  for (size_type i=0; i<lfsv_v_s.size(); i++)
482  div_phi_v_s[i] += jac_phi_v_s[i][d][d];
483  }
484 
485  // compute velocity value, jacobian, and divergence
486  Range_V val_u_s(0.0);
487  Dune::FieldMatrix<RF,dim,dim> jac_u_s(0.0);
488  RF div_u_s(0.0);
489  for (size_type i=0; i<lfsu_v_s.size(); i++){
490  val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
491  jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
492  div_u_s += x_s(lfsu_v_s,i) * div_phi_v_s[i];
493  }
494 
495  // compute pressure value
496  Range_P val_p_s(0.0);
497  for (size_type i=0; i<lfsu_p_s.size(); i++)
498  val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
499 
500  // face normal vector
501  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
502 
503  if (bctype == BC::VelocityDirichlet) {
504 
505  // Compute jump relative to Dirichlet value
506  typename V::Traits::RangeType u0;
507  v.evaluateGlobal(global,u0);
508  const Dune::FieldVector<DF,dimw> jump = val_u_s - u0;
509 
510  {// update residual (flux term)
511 
512  // compute flux
513  Dune::FieldVector<DF,dimw> flux(0.0);
514  add_compute_flux(jac_u_s,normal,flux);
515 
516  const RF factor = weight * mu;
517  for (size_t i=0; i<lfsv_v_s.size(); i++)
518  r_s.accumulate(lfsv_v_s, i, -(flux * phi_v_s[i]) * factor);
519  }
520 
521  {// update residual (symmetry term)
522  const RF factor = weight * mu;
523 
524  for (size_t i=0; i<lfsv_v_s.size(); i++){
525  Dune::FieldVector<DF,dimw> flux(0.0);
526  add_compute_flux(jac_phi_v_s[i],normal,flux);
527  r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux * jump) * factor);
528  }
529  }
530 
531  {// interior penalty
532  const RF factor = weight;
533  const RF gamma = ip_factor.getFaceIP(ig);
534  for (size_t i=0; i<lfsv_v_s.size(); i++)
535  r_s.accumulate(lfsv_v_s,i, gamma * (jump * phi_v_s[i]) * factor);
536  }
537 
538  {// pressure and incompressibility
539  const RF factor = weight;
540  const RF val_p = val_p_s;
541  for (size_t i=0; i<lfsv_v_s.size(); i++)
542  r_s.accumulate(lfsv_v_s, i, val_p * (phi_v_s[i] * normal) * factor);
543  for (size_t i=0; i<lfsv_p_s.size(); i++)
544  r_s.accumulate(lfsv_p_s, i, 0.5 * phi_p_s[i] * (jump*normal) * factor);
545  }
546 
547  } // DirichletVelocity
548 
549  if (bctype == BC::StressNeumann){
550  typename P::Traits::RangeType p0;
551  p.evaluateGlobal(global,p0);
552 
553  for (size_t i=0; i<lfsv_v_s.size(); i++){
554  const RF val = p0 * (normal*phi_v_s[i]) * weight;
555  r_s.accumulate(lfsv_v_s, i, val);
556  }
557  } // PressureDirichlet
558 
559  } // it - quadrature
560 
561  }
562 
563  private:
564 
565  template<class M, class RF>
566  static void contraction(const M & a, const M & b, RF & v)
567  {
568  v = 0;
569  const int an = a.N();
570  const int am = a.M();
571  for(int r=0; r<an; ++r)
572  for(int c=0; c<am; ++c)
573  v += a[r][c] * b[r][c];
574  }
575 
576  template<class DU, class R>
577  static void add_compute_flux(const DU & du, const R & n, R & result)
578  {
579  const int N = du.N();
580  const int M = du.M();
581  for(int r=0; r<N; ++r)
582  for(int c=0; c<M; ++c)
583  result[r] += du[r][c] * n[c];
584  }
585 
586  const F& f;
587  const B& b;
588  const V& v;
589  const P& p;
590  // values for NIPG / NIPG
591  int epsilon;
592  int qorder;
593  // physical parameters
594  double mu;
595  const IP & ip_factor;
596  };
597 
599  } // namespace PDELab
600 } // namespace Dune
601 
602 #endif
Definition: stokesdg_vecfem.hh:120
sparsity pattern generator
Definition: pattern.hh:30
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: stokesdg_vecfem.hh:413
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
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: stokesdg_vecfem.hh:39
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: stokesdg_vecfem.hh:239
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
StokesDGVectorFEM(const Dune::ParameterTree &configuration, const IP &ip_factor_, const RF mu_, const F &_f, const B &_b, const V &_v, const P &_p, int _qorder=4)
Definition: stokesdg_vecfem.hh:122
Basis::Traits::DomainLocal DomainLocal
export vector type of the local coordinates
Definition: stokesdg_vecfem.hh:31
static const int dim
Definition: adaptivity.hh:82
Definition: stokesparameter.hh:32
Definition: stokesdg_vecfem.hh:118
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
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: stokesdg_vecfem.hh:64
const IG & ig
Definition: common/constraints.hh:146
static const std::size_t dimRange
export dimension of the values
Definition: stokesdg_vecfem.hh:35
Basis::Traits::RangeField RangeField
export field type of the values
Definition: stokesdg_vecfem.hh:33
Definition: stokesdg_vecfem.hh:115
A local operator for solving the stokes equation using a DG discretization and a vector valued finite...
Definition: stokesdg_vecfem.hh:96
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: stokesdg_vecfem.hh:131
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: stokesdg_vecfem.hh:114
const EG & eg
Definition: common/constraints.hh:277
IP InteriorPenaltyFactor
Definition: stokesdg_vecfem.hh:111
Definition: stokesdg_vecfem.hh:29
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: stokesdg_vecfem.hh:56
Basis::Traits::RangeFieldType RangeField
export field type of the values
Definition: stokesdg_vecfem.hh:58
Definition: stokesdg_vecfem.hh:119