dune-pdelab  2.0.0
convectiondiffusion.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSION_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSION_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/common/static_assert.hh>
10 #include<dune/geometry/type.hh>
11 
12 #include<dune/geometry/referenceelements.hh>
13 #include<dune/geometry/quadraturerules.hh>
14 
17 
18 #include"defaultimp.hh"
19 #include"pattern.hh"
20 #include"flags.hh"
21 #include"idefault.hh"
22 
23 
24 namespace Dune {
25  namespace PDELab {
29 
39  template<typename GV, typename RF>
42  {
44  typedef GV GridViewType;
45 
47  enum {
49  dimDomain = GV::dimension
50  };
51 
53  typedef typename GV::Grid::ctype DomainFieldType;
54 
56  typedef Dune::FieldVector<DomainFieldType,dimDomain> DomainType;
57 
59  typedef Dune::FieldVector<DomainFieldType,dimDomain-1> IntersectionDomainType;
60 
62  typedef RF RangeFieldType;
63 
65  typedef Dune::FieldVector<RF,GV::dimensionworld> RangeType;
66 
68  typedef Dune::FieldMatrix<RangeFieldType,dimDomain,dimDomain> PermTensorType;
69 
71  typedef typename GV::Traits::template Codim<0>::Entity ElementType;
72  typedef typename GV::Intersection IntersectionType;
73  };
74 
76  template<class T, class Imp>
78  {
79  public:
80  typedef T Traits;
81 
83  typename Traits::RangeFieldType
84  f (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
85  typename Traits::RangeFieldType u) const
86  {
87  return asImp().f(e,x,u);
88  }
89 
91  typename Traits::RangeFieldType
92  w (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
93  typename Traits::RangeFieldType u) const
94  {
95  return asImp().w(e,x,u);
96  }
97 
99  typename Traits::RangeFieldType
100  v (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
101  typename Traits::RangeFieldType u) const
102  {
103  return asImp().v(e,x,u);
104  }
105 
107  typename Traits::PermTensorType
108  D (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
109  {
110  return asImp().D(e,x);
111  }
112 
114  typename Traits::RangeType
115  q (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
116  typename Traits::RangeFieldType u) const
117  {
118  return asImp().q(e,x,u);
119  }
120 
121  template<typename I>
123  const I & intersection, /*@\label{bcp:name}@*/
124  const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
125  ) const
126  {
127  return asImp().isDirichlet( intersection, coord );
128  }
129 
131  typename Traits::RangeFieldType
132  g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
133  {
134  return asImp().g(e,x);
135  }
136 
138  // Good: The dependence on u allows us to implement Robin type boundary conditions.
139  // Bad: This interface cannot be used for mixed finite elements where the flux is the essential b.c.
140  typename Traits::RangeFieldType
141  j (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
142  typename Traits::RangeFieldType u) const
143  {
144  return asImp().j(e,x);
145  }
146 
147  private:
148  Imp& asImp () {return static_cast<Imp &> (*this);}
149  const Imp& asImp () const {return static_cast<const Imp &>(*this);}
150  };
151 
152 
157  template<typename T>
159  : public Dune::PDELab::DirichletConstraintsParameters /*@\label{bcp:base}@*/
160  {
161  const typename T::Traits::GridViewType gv;
162  const T& t;
163 
164  public:
165  BCTypeParam_CD( const typename T::Traits::GridViewType& gv_, const T& t_ )
166  : gv( gv_ ), t( t_ )
167  {
168  }
169 
170  template<typename I>
172  const I & intersection, /*@\label{bcp:name}@*/
173  const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
174  ) const
175  {
176  return t.isDirichlet( intersection, coord );
177  }
178  };
179 
180 
185  template<typename T>
187  : public GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
188  typename T::Traits::RangeFieldType,
189  1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >
190  ,DirichletBoundaryCondition_CD<T> >
191  {
192  public:
193  typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
194  typename T::Traits::RangeFieldType,
195  1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> > Traits;
196 
198  DirichletBoundaryCondition_CD (const typename Traits::GridViewType& g_, const T& t_) : g(g_), t(t_) {}
199 
201  inline void evaluate (const typename Traits::ElementType& e,
202  const typename Traits::DomainType& x,
203  typename Traits::RangeType& y) const
204  {
205  y = t.g(e,x);
206  }
207 
208  inline const typename Traits::GridViewType& getGridView () const
209  {
210  return g;
211  }
212 
213  private:
214  typename Traits::GridViewType g;
215  const T& t;
216  };
217 
218 
224  template<typename T>
226  public NumericalJacobianApplyVolume<ConvectionDiffusion<T> >,
227  public NumericalJacobianApplyBoundary<ConvectionDiffusion<T> >,
228  public NumericalJacobianVolume<ConvectionDiffusion<T> >,
229  public NumericalJacobianBoundary<ConvectionDiffusion<T> >,
230  public FullVolumePattern,
233  {
234  public:
235  // pattern assembly flags
236  enum { doPatternVolume = true };
237 
238  // residual assembly flags
239  enum { doAlphaVolume = true };
240  enum { doAlphaBoundary = true };
241 
242  ConvectionDiffusion (T& param_, int intorder_=2)
243  : param(param_), intorder(intorder_)
244  {}
245 
246  // volume integral depending on test and ansatz functions
247  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
248  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
249  {
250  // domain and range field type
251  typedef typename LFSU::Traits::FiniteElementType::
252  Traits::LocalBasisType::Traits::DomainFieldType DF;
253  typedef typename LFSU::Traits::FiniteElementType::
254  Traits::LocalBasisType::Traits::RangeFieldType RF;
255  typedef typename LFSU::Traits::FiniteElementType::
256  Traits::LocalBasisType::Traits::JacobianType JacobianType;
257  typedef typename LFSU::Traits::FiniteElementType::
258  Traits::LocalBasisType::Traits::RangeType RangeType;
259 
260  typedef typename LFSU::Traits::SizeType size_type;
261 
262  // dimensions
263  const int dim = EG::Geometry::dimension;
264 
265  // select quadrature rule
266  Dune::GeometryType gt = eg.geometry().type();
267  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
268 
269  // evaluate diffusion tensor at cell center, assume it is constant over elements
270  typename T::Traits::PermTensorType tensor;
271  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
272  tensor = param.D(eg.entity(),localcenter);
273 
274  // evaluate nonlinearity w(x_i); we assume here it is a Lagrange basis!
275  std::vector<typename T::Traits::RangeFieldType> w(lfsu.size());
276  for (size_type i=0; i<lfsu.size(); i++)
277  w[i] = param.w(eg.entity(),localcenter,x(lfsu,i));
278 
279  // loop over quadrature points
280  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
281  {
282  // evaluate basis functions
283  std::vector<RangeType> phi(lfsu.size());
284  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
285 
286  // evaluate u
287  RF u=0.0;
288  for (size_type i=0; i<lfsu.size(); i++)
289  u += w[i]*phi[i];
290 
291  // evaluate source term
292  typename T::Traits::RangeFieldType f = param.f(eg.entity(),it->position(),u);
293 
294  // evaluate flux term
295  typename T::Traits::RangeType q = param.q(eg.entity(),it->position(),u);
296 
297  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
298  std::vector<JacobianType> js(lfsu.size());
299  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
300 
301  // transform gradients of shape functions to real element
302  const typename EG::Geometry::JacobianInverseTransposed jac
303  = eg.geometry().jacobianInverseTransposed(it->position());
304  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
305  for (size_type i=0; i<lfsu.size(); i++)
306  {
307  gradphi[i] = 0.0;
308  jac.umv(js[i][0],gradphi[i]);
309  }
310 
311  // v(u) compute gradient of u
312  Dune::FieldVector<RF,dim> vgradu(0.0);
313  for (size_type i=0; i<lfsu.size(); i++)
314  vgradu.axpy(w[i],gradphi[i]);
315  vgradu *= param.v(eg.entity(),it->position(),u);
316 
317  // compute D * v(u) * gradient of u
318  Dune::FieldVector<RF,dim> Dvgradu(0.0);
319  tensor.umv(vgradu,Dvgradu);
320 
321  // integrate (K grad u)*grad phi_i + a_0*u*phi_i
322  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
323  for (size_type i=0; i<lfsu.size(); i++)
324  r.accumulate(lfsu,i,( Dvgradu*gradphi[i] - q*gradphi[i] - f*phi[i] )*factor);
325  }
326  }
327 
328  // boundary integral
329  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
330  void alpha_boundary (const IG& ig,
331  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
332  R& r_s) const
333  {
334  // domain and range field type
335  typedef typename LFSV::Traits::FiniteElementType::
336  Traits::LocalBasisType::Traits::DomainFieldType DF;
337  typedef typename LFSV::Traits::FiniteElementType::
338  Traits::LocalBasisType::Traits::RangeFieldType RF;
339  typedef typename LFSV::Traits::FiniteElementType::
340  Traits::LocalBasisType::Traits::RangeType RangeType;
341 
342  typedef typename LFSV::Traits::SizeType size_type;
343 
344  // dimensions
345  const int dim = IG::dimension;
346 
347  // select quadrature rule
348  Dune::GeometryType gtface = ig.geometryInInside().type();
349  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
350 
351  // evaluate nonlinearity w(x_i); we assume here it is a Lagrange basis!
352  Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
353  Dune::FieldVector<DF,dim> facecenterinelement = ig.geometryInInside().global( facecenterlocal );
354  std::vector<typename T::Traits::RangeFieldType> w(lfsu_s.size());
355  for (size_type i=0; i<lfsu_s.size(); i++)
356  w[i] = param.w(*(ig.inside()),facecenterinelement,x_s(lfsu_s,i));
357 
358  // loop over quadrature points and integrate normal flux
359  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
360  {
361  // evaluate boundary condition type
362  // skip rest if we are on Dirichlet boundary
363  if( param.isDirichlet( ig.intersection(), it->position() ) )
364  continue;
365 
366  // position of quadrature point in local coordinates of element
367  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
368 
369  // evaluate test shape functions
370  std::vector<RangeType> phi(lfsv_s.size());
371  lfsv_s.finiteElement().localBasis().evaluateFunction(local,phi);
372 
373  // evaluate u
374  RF u=0.0;
375  for (size_type i=0; i<lfsu_s.size(); i++)
376  u += w[i]*phi[i];
377 
378  // evaluate flux boundary condition
379  typename T::Traits::RangeFieldType j;
380  j = param.j(*(ig.inside()),local,u);
381 
382  // integrate j
383  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
384  for (size_type i=0; i<lfsv_s.size(); i++)
385  r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
386  }
387  }
388 
390  void setTime (double t)
391  {
392  param.setTime(t);
393  }
394 
395  private:
396  T& param;
397  int intorder;
398  };
399 
401  } // namespace PDELab
402 } // namespace Dune
403 
404 #endif
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: convectiondiffusion.hh:65
traits class for two phase parameter class
Definition: convectiondiffusion.hh:41
Traits::PermTensorType D(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
tensor diffusion coefficient
Definition: convectiondiffusion.hh:108
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: convectiondiffusion.hh:71
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Dirichlet boundary condition value.
Definition: convectiondiffusion.hh:132
Definition: convectiondiffusion.hh:236
T Traits
Definition: convectiondiffusion.hh:80
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusion.hh:248
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: convectiondiffusion.hh:56
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:118
GV GridViewType
the grid view
Definition: convectiondiffusion.hh:44
VTKWriter & w
Definition: function.hh:1102
Traits::RangeFieldType w(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
nonlinearity under gradient
Definition: convectiondiffusion.hh:92
Definition: common/constraintsparameters.hh:24
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
dimension of the domain
Definition: convectiondiffusion.hh:49
static const int dim
Definition: adaptivity.hh:82
RF RangeFieldType
Export type for range field.
Definition: convectiondiffusion.hh:62
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:49
Definition: convectiondiffusion.hh:240
Traits::RangeFieldType v(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
scalar nonlinearity in diffusion coefficient
Definition: convectiondiffusion.hh:100
base class for parameter class
Definition: convectiondiffusion.hh:77
Dune::FieldMatrix< RangeFieldType, dimDomain, dimDomain > PermTensorType
permeability tensor type
Definition: convectiondiffusion.hh:68
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Traits::RangeFieldType j(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
Neumann boundary condition.
Definition: convectiondiffusion.hh:141
Dune::PDELab::GridFunctionTraits< typename T::Traits::GridViewType, typename T::Traits::RangeFieldType, 1, Dune::FieldVector< typename T::Traits::RangeFieldType, 1 > > Traits
Definition: convectiondiffusion.hh:195
GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: convectiondiffusion.hh:53
GV::Intersection IntersectionType
Definition: convectiondiffusion.hh:72
const IG & ig
Definition: common/constraints.hh:146
leaf of a function tree
Definition: function.hh:577
void setTime(double t)
set time in parameter class
Definition: convectiondiffusion.hh:390
bool isDirichlet(const I &intersection, const Dune::FieldVector< typename I::ctype, I::dimension-1 > &coord) const
Definition: convectiondiffusion.hh:122
DirichletBoundaryCondition_CD(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: convectiondiffusion.hh:198
Definition: convectiondiffusion.hh:225
BCTypeParam_CD(const typename T::Traits::GridViewType &gv_, const T &t_)
Definition: convectiondiffusion.hh:165
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
const F & f
Definition: common/constraints.hh:145
traits class holding the function signature, same as in local function
Definition: function.hh:176
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: convectiondiffusion.hh:59
sparsity pattern generator
Definition: pattern.hh:14
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Evaluate the GridFunction at given position.
Definition: convectiondiffusion.hh:201
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: convectiondiffusion.hh:158
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:115
const EG & eg
Definition: common/constraints.hh:277
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusion.hh:330
ConvectionDiffusion(T &param_, int intorder_=2)
Definition: convectiondiffusion.hh:242
Definition: convectiondiffusion.hh:239
Traits::RangeFieldType f(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
source/reaction term
Definition: convectiondiffusion.hh:84
bool isDirichlet(const I &intersection, const Dune::FieldVector< typename I::ctype, I::dimension-1 > &coord) const
Definition: convectiondiffusion.hh:171
Definition: convectiondiffusion.hh:186
Traits::RangeType q(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
nonlinear flux vector
Definition: convectiondiffusion.hh:115
const E & e
Definition: interpolate.hh:172
const Traits::GridViewType & getGridView() const
Definition: convectiondiffusion.hh:208