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