dune-pdelab  2.5-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>
12 
16 
17 #include"defaultimp.hh"
18 #include"pattern.hh"
19 #include"idefault.hh"
20 #include"flags.hh"
21 #include"l2.hh"
22 #include"stokesparameter.hh"
23 #include"navierstokesmass.hh"
24 
25 namespace Dune {
26  namespace PDELab {
27 
31 
57  template<typename P>
59  public FullVolumePattern,
61  public InstationaryLocalOperatorDefaultMethods<typename P::Traits::RangeField>
62  {
63  public:
66 
67  static const bool navier = P::assemble_navier;
68  static const bool full_tensor = P::assemble_full_tensor;
69 
70  // pattern assembly flags
71  enum { doPatternVolume = true };
72 
73  // residual assembly flags
74  enum { doAlphaVolume = true };
75  enum { doLambdaVolume = true };
76  enum { doLambdaBoundary = true };
77 
78  using PhysicalParameters = P;
79 
80  TaylorHoodNavierStokes (const PhysicalParameters & p, int superintegration_order_ = 0)
81 
82  : _p(p)
83  , superintegration_order(superintegration_order_)
84  {}
85 
86  // volume integral depending on test and ansatz functions
87  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
88  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
89  {
90  // define types
91  using namespace Indices;
92  using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
93  using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
94  using LFSU_P = TypeTree::Child<LFSU,_1>;
95  using RF = typename LFSU_V::Traits::FiniteElementType::
96  Traits::LocalBasisType::Traits::RangeFieldType;
97  using RT_V = typename LFSU_V::Traits::FiniteElementType::
98  Traits::LocalBasisType::Traits::RangeType;
99  using JacobianType_V = typename LFSU_V::Traits::FiniteElementType::
100  Traits::LocalBasisType::Traits::JacobianType;
101  using RT_P = typename LFSU_P::Traits::FiniteElementType::
102  Traits::LocalBasisType::Traits::RangeType;
103 
104  // extract local function spaces
105  const auto& lfsu_v_pfs = child(lfsu,_0);
106  const unsigned int vsize = lfsu_v_pfs.child(0).size();
107  const auto& lfsu_p = child(lfsu,_1);
108  const unsigned int psize = lfsu_p.size();
109 
110  // dimensions
111  const int dim = EG::Geometry::mydimension;
112 
113  // get geometry
114  auto geo = eg.geometry();
115 
116  // determine quadrature order
117  const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
118  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
119  const int jac_order = geo.type().isSimplex() ? 0 : 1;
120  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
121 
122  // Initialize vectors outside for loop
123  typename EG::Geometry::JacobianInverseTransposed jac;
124  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
125  std::vector<RT_P> psi(psize);
126  Dune::FieldVector<RF,dim> vu(0.0);
127  std::vector<RT_V> phi(vsize);
128  Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
129 
130  // loop over quadrature points
131  for (const auto& ip : quadratureRule(geo,qorder))
132  {
133  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
134  std::vector<JacobianType_V> js(vsize);
135  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
136 
137  // transform gradient to real element
138  jac = geo.jacobianInverseTransposed(ip.position());
139  for (size_t i=0; i<vsize; i++)
140  {
141  gradphi[i] = 0.0;
142  jac.umv(js[i][0],gradphi[i]);
143  }
144 
145  // evaluate basis functions
146  lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
147 
148  // compute u (if Navier term enabled)
149  if(navier)
150  {
151  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
152 
153  for(int d=0; d<dim; ++d)
154  {
155  vu[d] = 0.0;
156  const auto& lfsu_v = lfsu_v_pfs.child(d);
157  for (size_t i=0; i<lfsu_v.size(); i++)
158  vu[d] += x(lfsu_v,i) * phi[i];
159  }
160  }
161 
162  // Compute velocity jacobian
163  for(int d=0; d<dim; ++d){
164  jacu[d] = 0.0;
165  const auto& lfsu_v = lfsu_v_pfs.child(d);
166  for (size_t i=0; i<lfsu_v.size(); i++)
167  jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
168  }
169 
170  // compute pressure
171  RT_P func_p(0.0);
172  for (size_t i=0; i<lfsu_p.size(); i++)
173  func_p += x(lfsu_p,i) * psi[i];
174 
175  // Viscosity and density
176  const auto mu = _p.mu(eg,ip.position());
177  const auto rho = _p.rho(eg,ip.position());
178 
179  // geometric weight
180  const auto factor = ip.weight() * geo.integrationElement(ip.position());
181 
182  for(int d=0; d<dim; ++d){
183 
184  const auto& lfsu_v = lfsu_v_pfs.child(d);
185 
186  //compute u * grad u_d
187  const auto u_nabla_u = vu * jacu[d];
188 
189  for (size_t i=0; i<vsize; i++){
190 
191  // integrate grad u * grad phi_i
192  r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
193 
194  // integrate (grad u)^T * grad phi_i
195  if (full_tensor)
196  for(int dd=0; dd<dim; ++dd)
197  r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
198 
199  // integrate div phi_i * p
200  r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
201 
202  // integrate u * grad u * phi_i
203  if(navier)
204  r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
205  }
206 
207  }
208 
209  // integrate div u * psi_i
210  for (size_t i=0; i<psize; i++)
211  for(int d=0; d<dim; ++d)
212  // divergence of u is the trace of the velocity jacobian
213  r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
214 
215  }
216  }
217 
218 
219  // volume integral depending on test functions
220  template<typename EG, typename LFSV, typename R>
221  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
222  {
223  // define types
224  using namespace Indices;
225  using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
226  using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
227  using LFSV_P = TypeTree::Child<LFSV,_1>;
228  using RT_V = typename LFSV_V::Traits::FiniteElementType::
229  Traits::LocalBasisType::Traits::RangeType;
230  using RT_P = typename LFSV_P::Traits::FiniteElementType::
231  Traits::LocalBasisType::Traits::RangeType;
232 
233  // extract local function spaces
234  const auto& lfsv_v_pfs = child(lfsv,_0);
235  const unsigned int vsize = lfsv_v_pfs.child(0).size();
236  const auto& lfsv_p = child(lfsv,_1);
237  const unsigned int psize = lfsv_p.size();
238 
239  // dimensions
240  const int dim = EG::Geometry::mydimension;
241 
242  // get cell
243  const auto& cell = eg.entity();
244 
245  // get geometry
246  auto geo = eg.geometry();
247 
248  // determine quadrature order
249  const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
250  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
251  const int qorder = 2*v_order + det_jac_order + superintegration_order;
252 
253  // Initialize vectors outside for loop
254  std::vector<RT_V> phi(vsize);
255  std::vector<RT_P> psi(psize);
256 
257  // loop over quadrature points
258  for (const auto& ip : quadratureRule(geo,qorder))
259  {
260  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
261 
262  lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
263 
264  // forcing term
265  const auto f1 = _p.f(cell,ip.position());
266 
267  // geometric weight
268  const auto factor = ip.weight() * geo.integrationElement(ip.position());
269 
270  for(int d=0; d<dim; ++d){
271 
272  const auto& lfsv_v = lfsv_v_pfs.child(d);
273 
274  for (size_t i=0; i<vsize; i++)
275  {
276  // integrate f1 * phi_i
277  r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
278  }
279 
280  }
281 
282  const auto g2 = _p.g2(eg,ip.position());
283 
284  // integrate div u * psi_i
285  for (size_t i=0; i<psize; i++)
286  {
287  r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
288  }
289 
290  }
291  }
292 
293 
294  // residual of boundary term
295  template<typename IG, typename LFSV, typename R>
296  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
297  {
298  // define types
299  using namespace Indices;
300  using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
301  using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
302  using RT_V = typename LFSV_V::Traits::FiniteElementType::
303  Traits::LocalBasisType::Traits::RangeType;
304 
305  // extract local velocity function spaces
306  const auto& lfsv_v_pfs = child(lfsv,_0);
307  const unsigned int vsize = lfsv_v_pfs.child(0).size();
308 
309  // dimensions
310  static const unsigned int dim = IG::dimension;
311 
312  // get geometry
313  auto geo = ig.geometry();
314 
315  // Get geometry of intersection in local coordinates of inside cell
316  auto geo_in_inside = ig.geometryInInside();
317 
318  // determine quadrature order
319  const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
320  const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
321  const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
322  const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
323 
324  // Initialize vectors outside for loop
325  std::vector<RT_V> phi(vsize);
326 
327  // loop over quadrature points and integrate normal flux
328  for (const auto& ip : quadratureRule(geo,qorder))
329  {
330  // evaluate boundary condition type
331  auto bctype = _p.bctype(ig,ip.position());
332 
333  // skip rest if we are on Dirichlet boundary
334  if (bctype != BC::StressNeumann)
335  continue;
336 
337  // position of quadrature point in local coordinates of element
338  auto local = geo_in_inside.global(ip.position());
339 
340  // evaluate basis functions
341  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
342 
343  const auto factor = ip.weight() * geo.integrationElement(ip.position());
344  const auto normal = ig.unitOuterNormal(ip.position());
345 
346  // evaluate flux boundary condition
347  const auto neumann_stress = _p.j(ig,ip.position(),normal);
348 
349  for(unsigned int d=0; d<dim; ++d)
350  {
351 
352  const auto& 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  // define types
369  using namespace Indices;
370  using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
371  using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
372  using LFSU_P = TypeTree::Child<LFSU,_1>;
373  using RF = typename LFSU_V::Traits::FiniteElementType::
374  Traits::LocalBasisType::Traits::RangeFieldType;
375  using RT_V = typename LFSU_V::Traits::FiniteElementType::
376  Traits::LocalBasisType::Traits::RangeType;
377  using JacobianType_V = typename LFSU_V::Traits::FiniteElementType::
378  Traits::LocalBasisType::Traits::JacobianType;
379  using RT_P = typename LFSU_P::Traits::FiniteElementType::
380  Traits::LocalBasisType::Traits::RangeType;
381 
382  // extract local function spaces
383  const auto& lfsu_v_pfs = child(lfsu,_0);
384  const unsigned int vsize = lfsu_v_pfs.child(0).size();
385  const auto& lfsu_p = child(lfsu,_1);
386  const unsigned int psize = lfsu_p.size();
387 
388  // dimensions
389  const int dim = EG::Geometry::mydimension;
390 
391  // get geometry
392  auto geo = eg.geometry();
393 
394  // determine quadrature order
395  const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
396  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
397  const int jac_order = geo.type().isSimplex() ? 0 : 1;
398  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
399 
400  // Initialize vectors outside for loop
401  typename EG::Geometry::JacobianInverseTransposed jac;
402  std::vector<JacobianType_V> js(vsize);
403  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
404  std::vector<RT_P> psi(psize);
405  std::vector<RT_V> phi(vsize);
406  Dune::FieldVector<RF,dim> vu(0.0);
407  Dune::FieldVector<RF,dim> gradu_d(0.0);
408 
409  // loop over quadrature points
410  for (const auto& ip : quadratureRule(geo,qorder))
411  {
412  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
413  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
414 
415  // transform gradient to real element
416  jac = geo.jacobianInverseTransposed(ip.position());
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  lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
425 
426  // compute u (if Navier term enabled)
427  if(navier){
428  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
429 
430  for(int d = 0; d < dim; ++d){
431  vu[d] = 0.0;
432  const auto& lfsv_v = lfsu_v_pfs.child(d);
433  for(size_t l = 0; l < vsize; ++l)
434  vu[d] += x(lfsv_v,l) * phi[l];
435  }
436  }
437 
438  // Viscosity and density
439  const auto mu = _p.mu(eg,ip.position());
440  const auto rho = _p.rho(eg,ip.position());
441 
442  const auto factor = ip.weight() * geo.integrationElement(ip.position());
443 
444  for(int d=0; d<dim; ++d){
445 
446  const auto& lfsv_v = lfsu_v_pfs.child(d);
447  const auto& lfsu_v = lfsv_v;
448 
449  // Derivatives of d-th velocity component
450  gradu_d = 0.0;
451  if(navier)
452  for(size_t l =0; l < vsize; ++l)
453  gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
454 
455  for (size_t i=0; i<lfsv_v.size(); i++){
456 
457  // integrate grad phi_u_i * grad phi_v_i (viscous force)
458  for (size_t j=0; j<lfsv_v.size(); j++){
459  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
460 
461  // integrate (grad phi_u_i)^T * grad phi_v_i (viscous force)
462  if(full_tensor)
463  for(int dd=0; dd<dim; ++dd){
464  const auto& lfsu_v = lfsu_v_pfs.child(dd);
465  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
466  }
467 
468  }
469 
470  // integrate grad_d phi_v_d * p_u (pressure force)
471  for (size_t j=0; j<psize; j++)
472  mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
473 
474  if(navier){
475  for(int k =0; k < dim; ++k){
476  const auto& lfsu_v = lfsu_v_pfs.child(k);
477 
478  const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
479 
480  for(size_t j=0; j< lfsu_v.size(); ++j)
481  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
482  } // k
483  }
484 
485  if(navier){
486  const auto pre_factor = factor * rho * phi[i];
487  for(size_t j=0; j< lfsu_v.size(); ++j)
488  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
489  }
490 
491  }
492 
493  for (size_t i=0; i<psize; i++){
494  for (size_t j=0; j<lfsu_v.size(); j++)
495  mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
496  }
497  } // d
498  } // it
499  }
500 
501  private:
502  const P& _p;
503  const int superintegration_order;
504  };
505 
507  } // namespace PDELab
508 } // namespace Dune
509 
510 #endif // DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
static const int dim
Definition: adaptivity.hh:83
const IG & ig
Definition: constraints.hh:148
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:88
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: stokesparameter.hh:32
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Definition: taylorhoodnavierstokes.hh:71
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:221
sparsity pattern generator
Definition: pattern.hh:13
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
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: taylorhoodnavierstokes.hh:365
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:58
Definition: taylorhoodnavierstokes.hh:75
Default flags for all local operators.
Definition: flags.hh:18
Definition: taylorhoodnavierstokes.hh:76
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:68
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:78
TaylorHoodNavierStokes(const PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:80
const P & p
Definition: constraints.hh:147
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:296
Definition: taylorhoodnavierstokes.hh:74
static const bool navier
Definition: taylorhoodnavierstokes.hh:67