dune-pdelab  2.4-dev
taylorhoodnavierstokes.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/quadraturerules.hh>
13 
14 #include"defaultimp.hh"
15 #include"pattern.hh"
16 #include"idefault.hh"
17 #include"flags.hh"
18 #include"l2.hh"
19 #include"stokesparameter.hh"
20 #include"navierstokesmass.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, int superintegration_order_ = 0)
78 
79  : _p(p)
80  , superintegration_order(superintegration_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 int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
118  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
119  const int jac_order = gt.isSimplex() ? 0 : 1;
120  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
121  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
122 
123  // loop over quadrature points
124  for (const auto& ip : rule)
125  {
126  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
127  std::vector<JacobianType_V> js(vsize);
128  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
129 
130  // transform gradient to real element
131  const typename EG::Geometry::JacobianInverseTransposed jac =
132  eg.geometry().jacobianInverseTransposed(ip.position());
133  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
134  for (size_t i=0; i<vsize; i++)
135  {
136  gradphi[i] = 0.0;
137  jac.umv(js[i][0],gradphi[i]);
138  }
139 
140  // evaluate basis functions
141  std::vector<RT_P> psi(psize);
142  lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
143 
144  // compute u (if Navier term enabled)
145  Dune::FieldVector<RF,dim> vu(0.0);
146 
147  std::vector<RT_V> phi(vsize);
148  if(navier)
149  {
150  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
151 
152  for(int d=0; d<dim; ++d)
153  {
154  const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
155  for (size_t i=0; i<lfsu_v.size(); i++)
156  vu[d] += x(lfsu_v,i) * phi[i];
157  }
158  }
159 
160  // Compute velocity jacobian
161  Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
162  for(int d=0; d<dim; ++d){
163  const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
164  for (size_t i=0; i<lfsu_v.size(); i++)
165  jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
166  }
167 
168  // compute pressure
169  RT_P func_p(0.0);
170  for (size_t i=0; i<lfsu_p.size(); i++)
171  func_p += x(lfsu_p,i) * psi[i];
172 
173  // Viscosity and density
174  const RF mu = _p.mu(eg,ip.position());
175  const RF rho = _p.rho(eg,ip.position());
176 
177  // geometric weight
178  const RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
179 
180  for(int d=0; d<dim; ++d){
181 
182  const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
183 
184  //compute u * grad u_d
185  const RF u_nabla_u = vu * jacu[d];
186 
187  for (size_t i=0; i<vsize; i++){
188 
189  // integrate grad u * grad phi_i
190  r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
191 
192  // integrate (grad u)^T * grad phi_i
193  if (full_tensor)
194  for(int dd=0; dd<dim; ++dd)
195  r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
196 
197  // integrate div phi_i * p
198  r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
199 
200  // integrate u * grad u * phi_i
201  if(navier)
202  r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
203  }
204 
205  }
206 
207  // integrate div u * psi_i
208  for (size_t i=0; i<psize; i++)
209  for(int d=0; d<dim; ++d)
210  // divergence of u is the trace of the velocity jacobian
211  r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
212 
213  }
214  }
215 
216 
217  // volume integral depending on test functions
218  template<typename EG, typename LFSV, typename R>
219  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
220  {
221  // dimensions
222  const int dim = EG::Geometry::dimension;
223 
224  // extract local function spaces
225  typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
226  const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
227 
228  typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
229  const unsigned int vsize = lfsv_v_pfs.child(0).size();
230 
231  // domain and range field type
232  typedef typename LFSV_V::Traits::FiniteElementType::
233  Traits::LocalBasisType::Traits::RangeFieldType RF;
234  typedef typename LFSV_V::Traits::FiniteElementType::
235  Traits::LocalBasisType::Traits::RangeType RT_V;
236 
237  typedef typename LFSV::template Child<1>::Type LFSV_P;
238  const LFSV_P& lfsv_p = lfsv.template child<1>();
239  const unsigned int psize = lfsv_p.size();
240 
241  typedef typename LFSV_V::Traits::FiniteElementType::
242  Traits::LocalBasisType::Traits::DomainFieldType DF;
243  typedef typename LFSV_P::Traits::FiniteElementType::
244  Traits::LocalBasisType::Traits::RangeType RT_P;
245 
246  // select quadrature rule
247  Dune::GeometryType gt = eg.geometry().type();
248  const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
249  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
250  const int qorder = 2*v_order + det_jac_order + superintegration_order;
251  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
252 
253  // loop over quadrature points
254  for (const auto& ip : rule)
255  {
256  std::vector<RT_V> phi(vsize);
257  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
258 
259  std::vector<RT_P> psi(psize);
260  lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
261 
262  // forcing term
263  const Dune::FieldVector<RF,dim> f1 = _p.f(eg,ip.position());
264 
265  // geometric weight
266  const RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
267 
268  for(int d=0; d<dim; ++d){
269 
270  const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
271 
272  for (size_t i=0; i<vsize; i++)
273  {
274  // integrate f1 * phi_i
275  r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
276  }
277 
278  }
279 
280  const RF g2 = _p.g2(eg,ip.position());
281 
282  // integrate div u * psi_i
283  for (size_t i=0; i<psize; i++)
284  {
285  r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
286  }
287 
288  }
289  }
290 
291 
292  // residual of boundary term
293  template<typename IG, typename LFSV, typename R>
294  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
295  {
296  // dimensions
297  static const unsigned int dim = IG::Geometry::dimension;
298  static const unsigned int dimw = IG::Geometry::dimensionworld;
299 
300  // extract local velocity function spaces
301  typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
302  const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
303 
304  typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
305  const unsigned int vsize = lfsv_v_pfs.child(0).size();
306 
307  typedef typename LFSV_V::Traits::FiniteElementType::
308  Traits::LocalBasisType::Traits::RangeType RT_V;
309 
310  // the range field type (equal for velocity and pressure)
311  typedef typename LFSV_V::Traits::FiniteElementType::
312  Traits::LocalBasisType::Traits::RangeFieldType RF;
313 
314  // the domain field type (equal for velocity and pressure)
315  typedef typename LFSV_V::Traits::FiniteElementType::
316  Traits::LocalBasisType::Traits::DomainFieldType DF;
317 
318  // select quadrature rule
319  Dune::GeometryType gtface = ig.geometryInInside().type();
320  const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
321  const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
322  const int jac_order = gtface.isSimplex() ? 0 : 1;
323  const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
324  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
325 
326  // loop over quadrature points and integrate normal flux
327  for (const auto& ip : rule)
328  {
329  // evaluate boundary condition type
330  typename BC::Type bctype = _p.bctype(ig,ip.position());
331 
332  // skip rest if we are on Dirichlet boundary
333  if (bctype != BC::StressNeumann)
334  continue;
335 
336  // position of quadrature point in local coordinates of element
337  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
338 
339  // evaluate basis functions
340  std::vector<RT_V> phi(vsize);
341  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
342 
343  const RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
344  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
345 
346  // evaluate flux boundary condition
347  const Dune::FieldVector<DF,dimw> neumann_stress = _p.j(ig,ip.position(),normal);
348 
349  for(unsigned int d=0; d<dim; ++d)
350  {
351 
352  const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
353 
354  for (size_t i=0; i<vsize; i++)
355  {
356  r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
357  }
358 
359  }
360  }
361  }
362 
363 
364  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
365  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
366  M& mat) const
367  {
368  // dimensions
369  const int dim = EG::Geometry::dimension;
370 
371 
372  // extract local function spaces
373  typedef typename LFSU::template Child<0>::Type LFSU_V_PFS;
374  const LFSU_V_PFS& lfsu_v_pfs = lfsu.template child<0>();
375  const unsigned int vsize = lfsu_v_pfs.child(0).size();
376 
377  typedef typename LFSU_V_PFS::template Child<0>::Type LFSU_V;
378 
379  // domain and range field type
380  typedef typename LFSU_V::Traits::FiniteElementType::
381  Traits::LocalBasisType::Traits::RangeFieldType RF;
382  typedef typename LFSU_V::Traits::FiniteElementType::
383  Traits::LocalBasisType::Traits::RangeType RT_V;
384  typedef typename LFSU_V::Traits::FiniteElementType::
385  Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
386 
387 
388  typedef typename LFSU::template Child<1>::Type LFSU_P;
389  const LFSU_P& lfsu_p = lfsu.template child<1>();
390  const unsigned int psize = lfsu_p.size();
391 
392  typedef typename LFSU_P::Traits::FiniteElementType::
393  Traits::LocalBasisType::Traits::DomainFieldType DF;
394 
395  typedef typename LFSU_P::Traits::FiniteElementType::
396  Traits::LocalBasisType::Traits::RangeType RT_P;
397 
398  // select quadrature rule
399  Dune::GeometryType gt = eg.geometry().type();
400  const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
401  const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
402  const int jac_order = gt.isSimplex() ? 0 : 1;
403  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
404  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
405 
406  // loop over quadrature points
407  for (const auto& ip : rule)
408  {
409  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
410  std::vector<JacobianType_V> js(vsize);
411  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
412 
413  // transform gradient to real element
414  const typename EG::Geometry::JacobianInverseTransposed jac =
415  eg.geometry().jacobianInverseTransposed(ip.position());
416  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
417  for (size_t i=0; i<vsize; i++)
418  {
419  gradphi[i] = 0.0;
420  jac.umv(js[i][0],gradphi[i]);
421  }
422 
423  // evaluate basis functions
424  std::vector<RT_P> psi(psize);
425  lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
426 
427  // compute u (if Navier term enabled)
428  std::vector<RT_V> phi(vsize);
429  Dune::FieldVector<RF,dim> vu(0.0);
430  if(navier){
431  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
432 
433  for(int d = 0; d < dim; ++d){
434  const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
435  for(size_t l = 0; l < vsize; ++l)
436  vu[d] += x(lfsv_v,l) * phi[l];
437  }
438  }
439 
440  // Viscosity and density
441  const RF mu = _p.mu(eg,ip.position());
442  const RF rho = _p.rho(eg,ip.position());
443 
444  const RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
445 
446  for(int d=0; d<dim; ++d){
447 
448  const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
449  const LFSU_V & lfsu_v = lfsv_v;
450 
451  // Derivatives of d-th velocity component
452  Dune::FieldVector<RF,dim> gradu_d(0.0);
453  if(navier)
454  for(size_t l =0; l < vsize; ++l)
455  gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
456 
457  for (size_t i=0; i<lfsv_v.size(); i++){
458 
459  // integrate grad phi_u_i * grad phi_v_i (viscous force)
460  for (size_t j=0; j<lfsv_v.size(); j++){
461  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
462 
463  // integrate (grad phi_u_i)^T * grad phi_v_i (viscous force)
464  if(full_tensor)
465  for(int dd=0; dd<dim; ++dd){
466  const LFSU_V & lfsu_v = lfsu_v_pfs.child(dd);
467  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
468  }
469 
470  }
471 
472  // integrate grad_d phi_v_d * p_u (pressure force)
473  for (size_t j=0; j<psize; j++)
474  mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
475 
476  if(navier){
477  for(int k =0; k < dim; ++k){
478  const LFSU_V & lfsu_v = lfsu_v_pfs.child(k);
479 
480  const RF pre_factor = factor * rho * gradu_d[k] * phi[i];
481 
482  for(size_t j=0; j< lfsu_v.size(); ++j)
483  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
484  } // k
485  }
486 
487  if(navier){
488  const RF pre_factor = factor * rho * phi[i];
489  for(size_t j=0; j< lfsu_v.size(); ++j)
490  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
491  }
492 
493  }
494 
495  for (size_t i=0; i<psize; i++){
496  for (size_t j=0; j<lfsu_v.size(); j++)
497  mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
498  }
499 
500  } // d
501 
502  } // it
503  }
504 
505  private:
506  const P& _p;
507  const int superintegration_order;
508  };
509 
511  } // namespace PDELab
512 } // namespace Dune
513 
514 #endif
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: taylorhoodnavierstokes.hh:365
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:85
TaylorHoodNavierStokes(const PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:77
static const bool navier
Definition: taylorhoodnavierstokes.hh:64
Definition: taylorhoodnavierstokes.hh:73
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:75
StokesBoundaryCondition BC
Boundary condition indicator type.
Definition: taylorhoodnavierstokes.hh:62
Definition: taylorhoodnavierstokes.hh:71
const IG & ig
Definition: constraints.hh:147
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:294
Definition: taylorhoodnavierstokes.hh:68
Definition: adaptivity.hh:27
Type
Definition: stokesparameter.hh:33
static const int dim
Definition: adaptivity.hh:83
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:55
const P & p
Definition: constraints.hh:146
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:219
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:65
Definition: taylorhoodnavierstokes.hh:72
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89