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