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