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