2#ifndef DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH
3#define DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH
11#include<dune/geometry/referenceelements.hh>
13#include<dune/pdelab/common/quadraturerules.hh>
14#include <dune/pdelab/common/function.hh>
15#include<dune/pdelab/constraints/common/constraintsparameters.hh>
17#include"defaultimp.hh"
39 template<
typename GV,
typename RF>
71 using IntersectionType =
typename GV::Intersection;
75 template<
class T,
class Imp>
82 typename Traits::RangeFieldType
83 f (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
84 typename Traits::RangeFieldType u)
const
86 return asImp().f(e,x,u);
90 typename Traits::RangeFieldType
91 w (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
92 typename Traits::RangeFieldType u)
const
94 return asImp().w(e,x,u);
98 typename Traits::RangeFieldType
99 v (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
100 typename Traits::RangeFieldType u)
const
102 return asImp().v(e,x,u);
106 typename Traits::PermTensorType
107 D (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x)
const
109 return asImp().D(e,x);
113 typename Traits::RangeType
114 q (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
115 typename Traits::RangeFieldType u)
const
117 return asImp().q(e,x,u);
122 const I & intersection,
126 return asImp().isDirichlet( intersection, coord );
130 typename Traits::RangeFieldType
131 g (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x)
const
133 return asImp().g(e,x);
139 typename Traits::RangeFieldType
140 j (
const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
141 typename Traits::RangeFieldType u)
const
143 return asImp().j(e,x);
147 Imp& asImp () {
return static_cast<Imp &
> (*this);}
148 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
160 const typename T::Traits::GridViewType gv;
164 BCTypeParam_CD(
const typename T::Traits::GridViewType& gv_,
const T& t_ )
171 const I & intersection,
175 return t.isDirichlet( intersection, coord );
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> >
193 typename T::Traits::RangeFieldType,
235 enum { doPatternVolume =
true };
238 enum { doAlphaVolume =
true };
239 enum { doAlphaBoundary =
true };
242 : param(param_), intorder(intorder_)
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
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;
259 const int dim = EG::Geometry::mydimension;
262 const auto& cell = eg.entity();
265 auto geo = eg.geometry();
269 auto localcenter = ref_el.position(0,0);
270 auto tensor = param.D(cell,localcenter);
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));
278 typename EG::Geometry::JacobianInverseTransposed jac;
281 std::vector<RangeType> phi(lfsu.size());
282 std::vector<JacobianType> js(lfsu.size());
283 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
288 for (
const auto& ip : quadratureRule(geo,intorder))
291 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
295 for (size_type i=0; i<lfsu.size(); i++)
299 auto f = param.f(cell,ip.position(),u);
302 auto q = param.q(cell,ip.position(),u);
305 lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
308 jac = geo.jacobianInverseTransposed(ip.position());
309 for (size_type i=0; i<lfsu.size(); i++)
312 jac.umv(js[i][0],gradphi[i]);
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);
323 tensor.umv(vgradu,Dvgradu);
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);
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,
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;
346 const auto& cell_inside = ig.inside();
349 auto geo = ig.geometry();
352 auto geo_in_inside = ig.geometryInInside();
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));
363 std::vector<RangeType> phi(lfsv_s.size());
366 for (
const auto& ip : quadratureRule(geo,intorder))
370 if( param.isDirichlet( ig.intersection(), ip.position() ) )
374 auto local = geo_in_inside.global(ip.position());
377 lfsv_s.finiteElement().localBasis().evaluateFunction(local,phi);
381 for (size_type i=0; i<lfsu_s.size(); i++)
385 auto j = param.j(cell_inside,local,u);
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);
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Definition: nonlinearconvectiondiffusionfem.hh:159
Definition: nonlinearconvectiondiffusionfem.hh:190
DirichletBoundaryCondition_CD(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: nonlinearconvectiondiffusionfem.hh:197
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
sparsity pattern generator
Definition: pattern.hh:14
leaf of a function tree
Definition: function.hh:302
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Default flags for all local operators.
Definition: flags.hh:19
Definition: nonlinearconvectiondiffusionfem.hh:232
void setTime(double t)
set time in parameter class
Definition: nonlinearconvectiondiffusionfem.hh:395
base class for parameter class
Definition: nonlinearconvectiondiffusionfem.hh:77
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Dirichlet boundary condition value.
Definition: nonlinearconvectiondiffusionfem.hh:131
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
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
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
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
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
Traits::PermTensorType D(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
tensor diffusion coefficient
Definition: nonlinearconvectiondiffusionfem.hh:107
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune namespace.
Definition: alignedallocator.hh:13
Static tag representing a codimension.
Definition: dimension.hh:24
Definition: constraintsparameters.hh:26
R RangeType
range type
Definition: function.hh:62
traits class holding the function signature, same as in local function
Definition: function.hh:183
traits class for two phase parameter class
Definition: nonlinearconvectiondiffusionfem.hh:41
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: nonlinearconvectiondiffusionfem.hh:52
@ dimDomain
dimension of the domain
Definition: nonlinearconvectiondiffusionfem.hh:48
GV GridViewType
the grid view
Definition: nonlinearconvectiondiffusionfem.hh:43
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: nonlinearconvectiondiffusionfem.hh:70
RF RangeFieldType
Export type for range field.
Definition: nonlinearconvectiondiffusionfem.hh:61
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:119
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:116
A unique label for each type of element that can occur in a grid.