2 #ifndef DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH 7 #include<dune/common/exceptions.hh> 8 #include<dune/common/fvector.hh> 9 #include<dune/geometry/type.hh> 11 #include<dune/geometry/referenceelements.hh> 12 #include<dune/geometry/quadraturerules.hh> 40 template<
typename GV,
typename RF>
42 struct ConvectionDiffusionParameterTraits
57 using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
66 using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
69 using PermTensorType = Dune::FieldMatrix<RangeFieldType,dimDomain,dimDomain>;
72 using ElementType =
typename GV::Traits::template Codim<0>::Entity;
77 template<
class T,
class Imp>
84 typename Traits::RangeFieldType
85 f (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
86 typename Traits::RangeFieldType u)
const 88 return asImp().f(e,x,u);
92 typename Traits::RangeFieldType
93 w (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
94 typename Traits::RangeFieldType u)
const 96 return asImp().w(e,x,u);
100 typename Traits::RangeFieldType
101 v (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
102 typename Traits::RangeFieldType u)
const 104 return asImp().v(e,x,u);
108 typename Traits::PermTensorType
109 D (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const 111 return asImp().D(e,x);
115 typename Traits::RangeType
116 q (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
117 typename Traits::RangeFieldType u)
const 119 return asImp().q(e,x,u);
124 const I & intersection,
125 const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
128 return asImp().isDirichlet( intersection, coord );
132 typename Traits::RangeFieldType
133 g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const 135 return asImp().g(e,x);
141 typename Traits::RangeFieldType
142 j (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
143 typename Traits::RangeFieldType u)
const 145 return asImp().j(e,x);
149 Imp& asImp () {
return static_cast<Imp &
> (*this);}
150 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
162 const typename T::Traits::GridViewType gv;
173 const I & intersection,
174 const Dune::FieldVector<typename I::ctype, I::dimension-1> & coord
177 return t.isDirichlet( intersection, coord );
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> >
195 typename T::Traits::RangeFieldType,
196 1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >;
237 enum { doPatternVolume =
true };
240 enum { doAlphaVolume =
true };
241 enum { doAlphaBoundary =
true };
244 : param(param_), intorder(intorder_)
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 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;
261 const int dim = EG::Geometry::mydimension;
264 const auto& cell = eg.entity();
267 auto geo = eg.geometry();
271 auto localcenter = ref_el.position(0,0);
272 auto tensor = param.D(cell,localcenter);
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));
280 typename EG::Geometry::JacobianInverseTransposed jac;
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);
293 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
297 for (size_type i=0; i<lfsu.size(); i++)
301 auto f = param.f(cell,ip.position(),u);
304 auto q = param.q(cell,ip.position(),u);
307 lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
310 jac = geo.jacobianInverseTransposed(ip.position());
311 for (size_type i=0; i<lfsu.size(); i++)
314 jac.umv(js[i][0],gradphi[i]);
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);
325 tensor.umv(vgradu,Dvgradu);
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);
335 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
337 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
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;
348 const int dim = IG::dimension;
351 const auto& cell_inside = ig.inside();
354 auto geo = ig.geometry();
357 auto geo_in_inside = ig.geometryInInside();
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));
368 std::vector<RangeType> phi(lfsv_s.size());
375 if( param.isDirichlet( ig.intersection(), ip.position() ) )
379 auto local = geo_in_inside.global(ip.position());
382 lfsv_s.finiteElement().localBasis().evaluateFunction(local,phi);
386 for (size_type i=0; i<lfsu_s.size(); i++)
390 auto j = param.j(cell_inside,local,u);
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);
414 #endif // DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH 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
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
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
leaf of a function tree
Definition: function.hh:575
Definition: nonlinearconvectiondiffusionfem.hh:159
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
BCTypeParam_CD(const typename T::Traits::GridViewType &gv_, const T &t_)
Definition: nonlinearconvectiondiffusionfem.hh:166
GV::Intersection IntersectionType
Definition: convectiondiffusionparameter.hh:58
NonLinearConvectionDiffusionFEM(T ¶m_, int intorder_=2)
Definition: nonlinearconvectiondiffusionfem.hh:243
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:117
Definition: nonlinearconvectiondiffusionfem.hh:226
base class for parameter class
Definition: nonlinearconvectiondiffusionfem.hh:78
Definition: constraintsparameters.hh:24
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: convectiondiffusionparameter.hh:51
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
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
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Dirichlet boundary condition value.
Definition: nonlinearconvectiondiffusionfem.hh:133
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:114
RF RangeFieldType
Export type for range field.
Definition: convectiondiffusionparameter.hh:48
bool isDirichlet(const I &intersection, const Dune::FieldVector< typename I::ctype, I::dimension-1 > &coord) const
Definition: nonlinearconvectiondiffusionfem.hh:123
bool isDirichlet(const I &intersection, const Dune::FieldVector< typename I::ctype, I::dimension-1 > &coord) const
Definition: nonlinearconvectiondiffusionfem.hh:172
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
T Traits
Definition: nonlinearconvectiondiffusionfem.hh:81
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
GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: convectiondiffusionparameter.hh:39
GV GridViewType
the grid view
Definition: convectiondiffusionparameter.hh:30
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:48
dimension of the domain
Definition: convectiondiffusionparameter.hh:35
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: convectiondiffusionparameter.hh:42
sparsity pattern generator
Definition: pattern.hh:13
const Entity & e
Definition: localfunctionspace.hh:111
void setTime(double t)
set time in parameter class
Definition: nonlinearconvectiondiffusionfem.hh:400
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: convectiondiffusionparameter.hh:57
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
Definition: nonlinearconvectiondiffusionfem.hh:187
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: nonlinearconvectiondiffusionfem.hh:249
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
Dune::FieldMatrix< RangeFieldType, dimDomain, dimDomain > PermTensorType
permeability tensor type
Definition: convectiondiffusionparameter.hh:54
DirichletBoundaryCondition_CD(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: nonlinearconvectiondiffusionfem.hh:199
const Traits::GridViewType & getGridView() const
Definition: nonlinearconvectiondiffusionfem.hh:209
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:1100
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: convectiondiffusionparameter.hh:45
Traits::PermTensorType D(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
tensor diffusion coefficient
Definition: nonlinearconvectiondiffusionfem.hh:109
R RangeType
range type
Definition: function.hh:60
traits class holding the function signature, same as in local function
Definition: function.hh:175
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Default flags for all local operators.
Definition: flags.hh:18