dune-pdelab  2.0.0
cg_stokes.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_CGSTOKES_HH
3 #define DUNE_PDELAB_CGSTOKES_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/common/static_assert.hh>
10 
11 #include<dune/geometry/type.hh>
12 #include<dune/geometry/quadraturerules.hh>
14 
15 #include"defaultimp.hh"
16 #include"pattern.hh"
17 #include"idefault.hh"
18 #include"flags.hh"
19 #include"l2.hh"
20 #include"stokesparameter.hh"
21 
22 namespace Dune {
23  namespace PDELab {
24 
28 
54  template<typename P>
56  public FullVolumePattern,
58  public InstationaryLocalOperatorDefaultMethods<typename P::Traits::RangeField>
59  {
60  public:
63 
64  static const bool navier = P::assemble_navier;
65  static const bool full_tensor = P::assemble_full_tensor;
66 
67  // pattern assembly flags
68  enum { doPatternVolume = true };
69 
70  // residual assembly flags
71  enum { doAlphaVolume = true };
72  enum { doLambdaVolume = true };
73  enum { doLambdaBoundary = true };
74 
75  typedef P PhysicalParameters;
76 
77  TaylorHoodNavierStokes (const PhysicalParameters & p, std::size_t quadrature_order = 4)
78 
79  : _p(p)
80  , _quadrature_order(quadrature_order)
81  {}
82 
83  // volume integral depending on test and ansatz functions
84  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
85  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
86  {
87  // dimensions
88  const int dim = EG::Geometry::dimension;
89 
90  // extract local function spaces
91  typedef typename LFSU::template Child<0>::Type LFSU_V_PFS;
92  const LFSU_V_PFS& lfsu_v_pfs = lfsu.template child<0>();
93 
94  typedef typename LFSU_V_PFS::template Child<0>::Type LFSU_V;
95  const unsigned int vsize = lfsu_v_pfs.child(0).size();
96 
97  // domain and range field type
98  typedef typename LFSU_V::Traits::FiniteElementType::
99  Traits::LocalBasisType::Traits::RangeFieldType RF;
100  typedef typename LFSU_V::Traits::FiniteElementType::
101  Traits::LocalBasisType::Traits::RangeType RT_V;
102  typedef typename LFSU_V::Traits::FiniteElementType::
103  Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
104 
105 
106  typedef typename LFSU::template Child<1>::Type LFSU_P;
107  const LFSU_P& lfsu_p = lfsu.template child<1>();
108  const unsigned int psize = lfsu_p.size();
109 
110  typedef typename LFSU_P::Traits::FiniteElementType::
111  Traits::LocalBasisType::Traits::DomainFieldType DF;
112  typedef typename LFSU_P::Traits::FiniteElementType::
113  Traits::LocalBasisType::Traits::RangeType RT_P;
114 
115  // select quadrature rule
116  Dune::GeometryType gt = eg.geometry().type();
117  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,_quadrature_order);
118 
119  // loop over quadrature points
120  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),
121  endit = rule.end();
122  it != endit;
123  ++it)
124  {
125  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
126  std::vector<JacobianType_V> js(vsize);
127  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
128 
129  // transform gradient to real element
130  const typename EG::Geometry::JacobianInverseTransposed jac =
131  eg.geometry().jacobianInverseTransposed(it->position());
132  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
133  for (size_t i=0; i<vsize; i++)
134  {
135  gradphi[i] = 0.0;
136  jac.umv(js[i][0],gradphi[i]);
137  }
138 
139  // evaluate basis functions
140  std::vector<RT_P> psi(psize);
141  lfsu_p.finiteElement().localBasis().evaluateFunction(it->position(),psi);
142 
143  // compute u (if Navier term enabled)
144  Dune::FieldVector<RF,dim> vu(0.0);
145 
146  std::vector<RT_V> phi(vsize);
147  if(navier)
148  {
149  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
150 
151  for(int d=0; d<dim; ++d)
152  {
153  const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
154  for (size_t i=0; i<lfsu_v.size(); i++)
155  vu[d] += x(lfsu_v,i) * phi[i];
156  }
157  }
158 
159  // Compute velocity jacobian
160  Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
161  for(int d=0; d<dim; ++d){
162  const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
163  for (size_t i=0; i<lfsu_v.size(); i++)
164  jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
165  }
166 
167  // compute pressure
168  RT_P func_p(0.0);
169  for (size_t i=0; i<lfsu_p.size(); i++)
170  func_p += x(lfsu_p,i) * psi[i];
171 
172  // Viscosity and density
173  const RF mu = _p.mu(eg,it->position());
174  const RF rho = _p.rho(eg,it->position());
175 
176  // geometric weight
177  const RF factor = it->weight() * eg.geometry().integrationElement(it->position());
178 
179  for(int d=0; d<dim; ++d){
180 
181  const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
182 
183  //compute u * grad u_d
184  const RF u_nabla_u = vu * jacu[d];
185 
186  for (size_t i=0; i<vsize; i++){
187 
188  // integrate grad u * grad phi_i
189  r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
190 
191  // integrate (grad u)^T * grad phi_i
192  if (full_tensor)
193  for(int dd=0; dd<dim; ++dd)
194  r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
195 
196  // integrate div phi_i * p
197  r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
198 
199  // integrate u * grad u * phi_i
200  if(navier)
201  r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
202  }
203 
204  }
205 
206  // compute divergence of u
207  RF divu(0.0);
208  for (size_t i=0; i<vsize; i++)
209  for(int d=0; d<dim; ++d)
210  divu += x(lfsu_v_pfs.child(d),i) * gradphi[i][d];
211 
212  // integrate div u * psi_i
213  for (size_t i=0; i<lfsu_p.size(); i++)
214  {
215  r.accumulate(lfsu_p,i, -1.0 * divu * psi[i] * factor);
216  }
217 
218  }
219  }
220 
221 
222  // volume integral depending on test functions
223  template<typename EG, typename LFSV, typename R>
224  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
225  {
226  // dimensions
227  const int dim = EG::Geometry::dimension;
228 
229  // extract local function spaces
230  typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
231  const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
232 
233  typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
234  const unsigned int vsize = lfsv_v_pfs.child(0).size();
235 
236  // domain and range field type
237  typedef typename LFSV_V::Traits::FiniteElementType::
238  Traits::LocalBasisType::Traits::RangeFieldType RF;
239  typedef typename LFSV_V::Traits::FiniteElementType::
240  Traits::LocalBasisType::Traits::RangeType RT_V;
241 
242  typedef typename LFSV::template Child<1>::Type LFSV_P;
243  const LFSV_P& lfsv_p = lfsv.template child<1>();
244  const unsigned int psize = lfsv_p.size();
245 
246  typedef typename LFSV_V::Traits::FiniteElementType::
247  Traits::LocalBasisType::Traits::DomainFieldType DF;
248  typedef typename LFSV_P::Traits::FiniteElementType::
249  Traits::LocalBasisType::Traits::RangeType RT_P;
250 
251  // select quadrature rule
252  Dune::GeometryType gt = eg.geometry().type();
253  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,_quadrature_order);
254 
255  // loop over quadrature points
256  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),
257  endit = rule.end();
258  it != endit;
259  ++it)
260  {
261  std::vector<RT_V> phi(vsize);
262  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
263 
264  std::vector<RT_P> psi(psize);
265  lfsv_p.finiteElement().localBasis().evaluateFunction(it->position(),psi);
266 
267  // forcing term
268  const Dune::FieldVector<RF,dim> f1 = _p.f(eg,it->position());
269 
270  // geometric weight
271  const RF factor = it->weight() * eg.geometry().integrationElement(it->position());
272 
273  for(int d=0; d<dim; ++d){
274 
275  const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
276 
277  for (size_t i=0; i<vsize; i++)
278  {
279  // integrate f1 * phi_i
280  r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
281  }
282 
283  }
284 
285  const RF g2 = _p.g2(eg,it->position());
286 
287  // integrate div u * psi_i
288  for (size_t i=0; i<lfsv_p.size(); i++)
289  {
290  r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
291  }
292 
293  }
294  }
295 
296 
297  // residual of boundary term
298  template<typename IG, typename LFSV, typename R>
299  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
300  {
301  // dimensions
302  static const unsigned int dim = IG::Geometry::dimension;
303  static const unsigned int dimw = IG::Geometry::dimensionworld;
304 
305  // extract local velocity function spaces
306  typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
307  const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
308 
309  typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
310  const unsigned int vsize = lfsv_v_pfs.child(0).size();
311 
312  typedef typename LFSV_V::Traits::FiniteElementType::
313  Traits::LocalBasisType::Traits::RangeType RT_V;
314 
315  // the range field type (equal for velocity and pressure)
316  typedef typename LFSV_V::Traits::FiniteElementType::
317  Traits::LocalBasisType::Traits::RangeFieldType RF;
318 
319  // the domain field type (equal for velocity and pressure)
320  typedef typename LFSV_V::Traits::FiniteElementType::
321  Traits::LocalBasisType::Traits::DomainFieldType DF;
322 
323  // select quadrature rule
324  Dune::GeometryType gtface = ig.geometryInInside().type();
325  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,_quadrature_order);
326 
327  // loop over quadrature points and integrate normal flux
328  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it = rule.begin(),
329  endit = rule.end();
330  it != endit;
331  ++it)
332  {
333  // evaluate boundary condition type
334  typename BC::Type bctype = _p.bctype(ig,it->position());
335 
336  // skip rest if we are on Dirichlet boundary
337  if (bctype != BC::StressNeumann)
338  continue;
339 
340  // position of quadrature point in local coordinates of element
341  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
342 
343  // evaluate basis functions
344  std::vector<RT_V> phi(vsize);
345  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
346 
347  const RF factor = it->weight() * ig.geometry().integrationElement(it->position());
348  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
349 
350  // evaluate flux boundary condition
351  const Dune::FieldVector<DF,dimw> neumann_stress = _p.j(ig,it->position(),normal);
352 
353  for(unsigned int d=0; d<dim; ++d)
354  {
355 
356  const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
357 
358  for (size_t i=0; i<vsize; i++)
359  {
360  r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
361  }
362 
363  }
364  }
365  }
366 
367 
368  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
369  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
370  M& mat) const
371  {
372  // dimensions
373  const int dim = EG::Geometry::dimension;
374 
375 
376  // extract local function spaces
377  typedef typename LFSU::template Child<0>::Type LFSU_V_PFS;
378  const LFSU_V_PFS& lfsu_v_pfs = lfsu.template child<0>();
379  const unsigned int vsize = lfsu_v_pfs.child(0).size();
380 
381  typedef typename LFSU_V_PFS::template Child<0>::Type LFSU_V;
382 
383  // domain and range field type
384  typedef typename LFSU_V::Traits::FiniteElementType::
385  Traits::LocalBasisType::Traits::RangeFieldType RF;
386  typedef typename LFSU_V::Traits::FiniteElementType::
387  Traits::LocalBasisType::Traits::RangeType RT_V;
388  typedef typename LFSU_V::Traits::FiniteElementType::
389  Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
390 
391 
392  typedef typename LFSU::template Child<1>::Type LFSU_P;
393  const LFSU_P& lfsu_p = lfsu.template child<1>();
394  const unsigned int psize = lfsu_p.size();
395 
396  typedef typename LFSU_P::Traits::FiniteElementType::
397  Traits::LocalBasisType::Traits::DomainFieldType DF;
398 
399  typedef typename LFSU_P::Traits::FiniteElementType::
400  Traits::LocalBasisType::Traits::RangeType RT_P;
401 
402  // select quadrature rule
403  Dune::GeometryType gt = eg.geometry().type();
404  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,_quadrature_order);
405 
406  // loop over quadrature points
407  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),
408  endit = rule.end();
409  it != endit;
410  ++it)
411  {
412  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
413  std::vector<JacobianType_V> js(vsize);
414  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
415 
416  // transform gradient to real element
417  const typename EG::Geometry::JacobianInverseTransposed jac =
418  eg.geometry().jacobianInverseTransposed(it->position());
419  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
420  for (size_t i=0; i<vsize; i++)
421  {
422  gradphi[i] = 0.0;
423  jac.umv(js[i][0],gradphi[i]);
424  }
425 
426  // evaluate basis functions
427  std::vector<RT_P> psi(psize);
428  lfsu_p.finiteElement().localBasis().evaluateFunction(it->position(),psi);
429 
430  // compute u (if Navier term enabled)
431  std::vector<RT_V> phi(vsize);
432  Dune::FieldVector<RF,dim> vu(0.0);
433  if(navier){
434  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
435 
436  for(int d = 0; d < dim; ++d){
437  const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
438  for(size_t l = 0; l < vsize; ++l)
439  vu[d] += x(lfsv_v,l) * phi[l];
440  }
441  }
442 
443  // Viscosity and density
444  const RF mu = _p.mu(eg,it->position());
445  const RF rho = _p.rho(eg,it->position());
446 
447  const RF factor = it->weight() * eg.geometry().integrationElement(it->position());
448 
449  for(int d=0; d<dim; ++d){
450 
451  const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
452  const LFSU_V & lfsu_v = lfsv_v;
453 
454  // Derivatives of d-th velocity component
455  Dune::FieldVector<RF,dim> gradu_d(0.0);
456  if(navier)
457  for(size_t l =0; l < vsize; ++l)
458  gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
459 
460  for (size_t i=0; i<lfsv_v.size(); i++){
461 
462  // integrate grad phi_u_i * grad phi_v_i (viscous force)
463  for (size_t j=0; j<lfsv_v.size(); j++){
464  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
465 
466  // integrate (grad phi_u_i)^T * grad phi_v_i (viscous force)
467  if(full_tensor)
468  for(int dd=0; dd<dim; ++dd){
469  const LFSU_V & lfsu_v = lfsu_v_pfs.child(dd);
470  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
471  }
472 
473  }
474 
475  // integrate grad_d phi_v_d * p_u (pressure force)
476  for (size_t j=0; j<lfsu_p.size(); j++)
477  mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
478 
479  if(navier){
480  for(int k =0; k < dim; ++k){
481  const LFSU_V & lfsu_v = lfsu_v_pfs.child(k);
482 
483  const RF pre_factor = factor * rho * gradu_d[k] * phi[i];
484 
485  for(size_t j=0; j< lfsu_v.size(); ++j)
486  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
487  } // k
488  }
489 
490  if(navier){
491  const RF pre_factor = factor * rho * phi[i];
492  for(size_t j=0; j< lfsu_v.size(); ++j)
493  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
494  }
495 
496  }
497 
498  for (size_t i=0; i<lfsu_p.size(); i++){
499  for (size_t j=0; j<lfsu_v.size(); j++)
500  mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
501  }
502 
503  } // d
504 
505  } // it
506  }
507 
508  private:
509  const P& _p;
510  const std::size_t _quadrature_order;
511  };
512 
513 
521  template< typename P >
523  : public NumericalJacobianApplyVolume<NavierStokesMass<P> >
524  , public FullVolumePattern
527  {
528  public:
529  // pattern assembly flags
530  enum { doPatternVolume = true };
531 
532  // residual assembly flags
533  enum { doAlphaVolume = true };
534 
535  NavierStokesMass (const P & p_, int intorder_=4)
536  : p(p_), intorder(intorder_)
537  {}
538 
539  // volume integral depending on test and ansatz functions
540  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
541  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
542  {
543  typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
544  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
545 
546  for(unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
547  {
548  scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
549  }
550  }
551 
552  // jacobian of volume term
553  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
554  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
555  M& mat) const
556  {
557  typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
558  const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
559 
560  for(unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
561  {
562  scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
563  }
564  }
565 
566  private:
567  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
568  void scalar_alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
569  R& r) const
570  {
571 
572  // Switches between local and global interface
573  typedef FiniteElementInterfaceSwitch<
574  typename LFSU::Traits::FiniteElementType
575  > FESwitch;
576  typedef BasisInterfaceSwitch<
577  typename FESwitch::Basis
578  > BasisSwitch;
579 
580  // domain and range field type
581  typedef typename BasisSwitch::DomainField DF;
582  typedef typename BasisSwitch::RangeField RF;
583  typedef typename BasisSwitch::Range RangeType;
584 
585  typedef typename LFSU::Traits::SizeType size_type;
586 
587  // dimensions
588  const int dim = EG::Geometry::dimension;
589 
590  // select quadrature rule
591  Dune::GeometryType gt = eg.geometry().type();
592  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
593 
594  // loop over quadrature points
595  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
596  {
597  // evaluate basis functions
598  std::vector<RangeType> phi(lfsu.size());
599  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(it->position(),phi);
600 
601  RF rho = p.rho(eg,it->position());
602  // evaluate u
603  RF u=0.0;
604  for (size_type i=0; i<lfsu.size(); i++)
605  u += x(lfsu,i)*phi[i];
606 
607  // u*phi_i
608  RF factor = it->weight() * rho * eg.geometry().integrationElement(it->position());
609 
610  for (size_type i=0; i<lfsu.size(); i++)
611  r.accumulate(lfsv,i, u*phi[i]*factor);
612  }
613  }
614 
615  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
616  void scalar_jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
617  M& mat) const
618  {
619 
620  // Switches between local and global interface
621  typedef FiniteElementInterfaceSwitch<
622  typename LFSU::Traits::FiniteElementType
623  > FESwitch;
624  typedef BasisInterfaceSwitch<
625  typename FESwitch::Basis
626  > BasisSwitch;
627 
628  // domain and range field type
629  typedef typename BasisSwitch::DomainField DF;
630  typedef typename BasisSwitch::RangeField RF;
631  typedef typename BasisSwitch::Range RangeType;
632  typedef typename LFSU::Traits::SizeType size_type;
633 
634  // dimensions
635  const int dim = EG::Geometry::dimension;
636 
637  // select quadrature rule
638  Dune::GeometryType gt = eg.geometry().type();
639  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
640 
641  // loop over quadrature points
642  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
643  {
644  // evaluate basis functions
645  std::vector<RangeType> phi(lfsu.size());
646  FESwitch::basis(lfsu.finiteElement()).evaluateFunction(it->position(),phi);
647 
648  // integrate phi_j*phi_i
649  RF rho = p.rho(eg,it->position());
650  RF factor = it->weight() * rho * eg.geometry().integrationElement(it->position());
651  for (size_type j=0; j<lfsu.size(); j++)
652  for (size_type i=0; i<lfsu.size(); i++)
653  mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
654  }
655  }
656 
657  const P & p;
658  int intorder;
659  };
660 
662  } // namespace PDELab
663 } // namespace Dune
664 
665 #endif
Type
Definition: stokesparameter.hh:33
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:85
Default flags for all local operators.
Definition: flags.hh:18
A local operator for the Navier-Stokes equations.
Definition: cg_stokes.hh:55
NavierStokesMass(const P &p_, int intorder_=4)
Definition: cg_stokes.hh:535
static const bool navier
Definition: cg_stokes.hh:64
P PhysicalParameters
Definition: cg_stokes.hh:75
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:541
TaylorHoodNavierStokes(const PhysicalParameters &p, std::size_t quadrature_order=4)
Definition: cg_stokes.hh:77
static const int dim
Definition: adaptivity.hh:82
StokesBoundaryCondition BC
Boundary condition indicator type.
Definition: cg_stokes.hh:62
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:299
Definition: stokesparameter.hh:32
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:224
const IG & ig
Definition: common/constraints.hh:146
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: cg_stokes.hh:554
Definition: cg_stokes.hh:522
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
static const bool full_tensor
Definition: cg_stokes.hh:65
const EG & eg
Definition: common/constraints.hh:277
R p(int i, D x)
Definition: qkdg.hh:62
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: cg_stokes.hh:369