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