dune-pdelab  2.4.1
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 TypeTree::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 TypeTree::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 geometry
243  auto geo = eg.geometry();
244 
245  // determine quadrature order
246  const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
247  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
248  const int qorder = 2*v_order + det_jac_order + superintegration_order;
249 
250  // Initialize vectors outside for loop
251  std::vector<RT_V> phi(vsize);
252  std::vector<RT_P> psi(psize);
253 
254  // loop over quadrature points
255  for (const auto& ip : quadratureRule(geo,qorder))
256  {
257  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
258 
259  lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
260 
261  // forcing term
262  const auto f1 = _p.f(eg,ip.position());
263 
264  // geometric weight
265  const auto factor = ip.weight() * geo.integrationElement(ip.position());
266 
267  for(int d=0; d<dim; ++d){
268 
269  const auto& lfsv_v = lfsv_v_pfs.child(d);
270 
271  for (size_t i=0; i<vsize; i++)
272  {
273  // integrate f1 * phi_i
274  r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
275  }
276 
277  }
278 
279  const auto g2 = _p.g2(eg,ip.position());
280 
281  // integrate div u * psi_i
282  for (size_t i=0; i<psize; i++)
283  {
284  r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
285  }
286 
287  }
288  }
289 
290 
291  // residual of boundary term
292  template<typename IG, typename LFSV, typename R>
293  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
294  {
295  // define types
296  using namespace TypeTree::Indices;
297  using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
298  using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
299  using RT_V = typename LFSV_V::Traits::FiniteElementType::
300  Traits::LocalBasisType::Traits::RangeType;
301 
302  // extract local velocity function spaces
303  const auto& lfsv_v_pfs = child(lfsv,_0);
304  const unsigned int vsize = lfsv_v_pfs.child(0).size();
305 
306  // dimensions
307  static const unsigned int dim = IG::dimension;
308 
309  // get geometry
310  auto geo = ig.geometry();
311 
312  // Get geometry of intersection in local coordinates of inside cell
313  auto geo_in_inside = ig.geometryInInside();
314 
315  // determine quadrature order
316  const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
317  const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
318  const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
319  const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
320 
321  // Initialize vectors outside for loop
322  std::vector<RT_V> phi(vsize);
323 
324  // loop over quadrature points and integrate normal flux
325  for (const auto& ip : quadratureRule(geo,qorder))
326  {
327  // evaluate boundary condition type
328  auto bctype = _p.bctype(ig,ip.position());
329 
330  // skip rest if we are on Dirichlet boundary
331  if (bctype != BC::StressNeumann)
332  continue;
333 
334  // position of quadrature point in local coordinates of element
335  auto local = geo_in_inside.global(ip.position());
336 
337  // evaluate basis functions
338  lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
339 
340  const auto factor = ip.weight() * geo.integrationElement(ip.position());
341  const auto normal = ig.unitOuterNormal(ip.position());
342 
343  // evaluate flux boundary condition
344  const auto neumann_stress = _p.j(ig,ip.position(),normal);
345 
346  for(unsigned int d=0; d<dim; ++d)
347  {
348 
349  const auto& lfsv_v = lfsv_v_pfs.child(d);
350 
351  for (size_t i=0; i<vsize; i++)
352  {
353  r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
354  }
355 
356  }
357  }
358  }
359 
360 
361  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
362  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
363  M& mat) const
364  {
365  // define types
366  using namespace TypeTree::Indices;
367  using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
368  using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
369  using LFSU_P = TypeTree::Child<LFSU,_1>;
370  using RF = typename LFSU_V::Traits::FiniteElementType::
371  Traits::LocalBasisType::Traits::RangeFieldType;
372  using RT_V = typename LFSU_V::Traits::FiniteElementType::
373  Traits::LocalBasisType::Traits::RangeType;
374  using JacobianType_V = typename LFSU_V::Traits::FiniteElementType::
375  Traits::LocalBasisType::Traits::JacobianType;
376  using RT_P = typename LFSU_P::Traits::FiniteElementType::
377  Traits::LocalBasisType::Traits::RangeType;
378 
379  // extract local function spaces
380  const auto& lfsu_v_pfs = child(lfsu,_0);
381  const unsigned int vsize = lfsu_v_pfs.child(0).size();
382  const auto& lfsu_p = child(lfsu,_1);
383  const unsigned int psize = lfsu_p.size();
384 
385  // dimensions
386  const int dim = EG::Geometry::mydimension;
387 
388  // get geometry
389  auto geo = eg.geometry();
390 
391  // determine quadrature order
392  const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
393  const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
394  const int jac_order = geo.type().isSimplex() ? 0 : 1;
395  const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
396 
397  // Initialize vectors outside for loop
398  typename EG::Geometry::JacobianInverseTransposed jac;
399  std::vector<JacobianType_V> js(vsize);
400  std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
401  std::vector<RT_P> psi(psize);
402  std::vector<RT_V> phi(vsize);
403  Dune::FieldVector<RF,dim> vu(0.0);
404  Dune::FieldVector<RF,dim> gradu_d(0.0);
405 
406  // loop over quadrature points
407  for (const auto& ip : quadratureRule(geo,qorder))
408  {
409  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
410  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
411 
412  // transform gradient to real element
413  jac = geo.jacobianInverseTransposed(ip.position());
414  for (size_t i=0; i<vsize; i++)
415  {
416  gradphi[i] = 0.0;
417  jac.umv(js[i][0],gradphi[i]);
418  }
419 
420  // evaluate basis functions
421  lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
422 
423  // compute u (if Navier term enabled)
424  if(navier){
425  lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
426 
427  for(int d = 0; d < dim; ++d){
428  vu[d] = 0.0;
429  const auto& lfsv_v = lfsu_v_pfs.child(d);
430  for(size_t l = 0; l < vsize; ++l)
431  vu[d] += x(lfsv_v,l) * phi[l];
432  }
433  }
434 
435  // Viscosity and density
436  const auto mu = _p.mu(eg,ip.position());
437  const auto rho = _p.rho(eg,ip.position());
438 
439  const auto factor = ip.weight() * geo.integrationElement(ip.position());
440 
441  for(int d=0; d<dim; ++d){
442 
443  const auto& lfsv_v = lfsu_v_pfs.child(d);
444  const auto& lfsu_v = lfsv_v;
445 
446  // Derivatives of d-th velocity component
447  gradu_d = 0.0;
448  if(navier)
449  for(size_t l =0; l < vsize; ++l)
450  gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
451 
452  for (size_t i=0; i<lfsv_v.size(); i++){
453 
454  // integrate grad phi_u_i * grad phi_v_i (viscous force)
455  for (size_t j=0; j<lfsv_v.size(); j++){
456  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
457 
458  // integrate (grad phi_u_i)^T * grad phi_v_i (viscous force)
459  if(full_tensor)
460  for(int dd=0; dd<dim; ++dd){
461  const auto& lfsu_v = lfsu_v_pfs.child(dd);
462  mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
463  }
464 
465  }
466 
467  // integrate grad_d phi_v_d * p_u (pressure force)
468  for (size_t j=0; j<psize; j++)
469  mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
470 
471  if(navier){
472  for(int k =0; k < dim; ++k){
473  const auto& lfsu_v = lfsu_v_pfs.child(k);
474 
475  const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
476 
477  for(size_t j=0; j< lfsu_v.size(); ++j)
478  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
479  } // k
480  }
481 
482  if(navier){
483  const auto pre_factor = factor * rho * phi[i];
484  for(size_t j=0; j< lfsu_v.size(); ++j)
485  mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
486  }
487 
488  }
489 
490  for (size_t i=0; i<psize; i++){
491  for (size_t j=0; j<lfsu_v.size(); j++)
492  mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
493  }
494  } // d
495  } // it
496  }
497 
498  private:
499  const P& _p;
500  const int superintegration_order;
501  };
502 
504  } // namespace PDELab
505 } // namespace Dune
506 
507 #endif
Definition: taylorhoodnavierstokes.hh:71
sparsity pattern generator
Definition: pattern.hh:13
static const bool navier
Definition: taylorhoodnavierstokes.hh:67
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
Definition: taylorhoodnavierstokes.hh:76
Definition: adaptivity.hh:27
Definition: stokesparameter.hh:32
TaylorHoodNavierStokes(const PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:80
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:88
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: taylorhoodnavierstokes.hh:362
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
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:68
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:58
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:78
Definition: taylorhoodnavierstokes.hh:74
const P & p
Definition: constraints.hh:147
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:293
Definition: taylorhoodnavierstokes.hh:75
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:221