2 #ifndef DUNE_PDELAB_TRANSPORTCCFV_HH
3 #define DUNE_PDELAB_TRANSPORTCCFV_HH
5 #include<dune/common/fvector.hh>
6 #include<dune/geometry/referenceelements.hh>
10 #include"../common/function.hh"
19 template<
typename GV,
typename RF>
35 typedef Dune::FieldVector<DomainFieldType,dimDomain>
DomainType;
44 typedef Dune::FieldVector<RF,GV::dimensionworld>
RangeType;
47 typedef typename GV::Traits::template Codim<0>::Entity
ElementType;
52 template<
class T,
class Imp>
59 typename Traits::RangeType
60 v (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
62 return asImp().v(e,x);
66 typename Traits::RangeFieldType
67 D (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
69 return asImp().D(e,x);
73 typename Traits::RangeFieldType
74 q (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
76 return asImp().q(e,x);
86 bc (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x)
const
88 return asImp().bc(is,x);
92 typename Traits::RangeFieldType
93 g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
95 return asImp().g(e,x);
99 typename Traits::RangeFieldType
100 j (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
102 return asImp().j(e,x);
106 Imp& asImp () {
return static_cast<Imp &
> (*this);}
107 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
118 BoundaryGridFunctionTraits<typename T::Traits::GridViewType,int,1,
119 Dune::FieldVector<int,1> >,
120 BoundaryConditionType_Transport<T> >
122 typename T::Traits::GridViewType gv;
154 :
public GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
155 typename T::Traits::RangeFieldType,
156 1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >
157 ,DirichletBoundaryCondition_Transport<T> >
161 typename T::Traits::RangeFieldType,
162 1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >
Traits;
198 template<
typename TP>
211 enum { dim = TP::Traits::GridViewType::dimension };
233 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
234 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
240 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
249 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
251 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
252 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
253 R& r_s, R& r_n)
const
256 typedef typename LFSU::Traits::FiniteElementType::
257 Traits::LocalBasisType::Traits::DomainFieldType DF;
258 typedef typename LFSU::Traits::FiniteElementType::
259 Traits::LocalBasisType::Traits::RangeFieldType RF;
262 const Dune::FieldVector<DF,IG::dimension-1>&
263 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
264 RF face_volume = ig.geometry().volume();
267 Dune::FieldVector<DF,IG::dimension> face_center_in_element = ig.geometryInInside().global(face_local);
270 typename TP::Traits::RangeType v(tp.v(*(ig.inside()),face_center_in_element));
273 RF vn = v*ig.centerUnitOuterNormal();
277 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
278 r_s.accumulate(lfsu_s,0,(u_upwind*vn)*face_volume);
280 cellinflux += vn*face_volume;
283 const Dune::FieldVector<DF,IG::dimension>&
284 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
285 const Dune::FieldVector<DF,IG::dimension>&
286 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
287 typename TP::Traits::RangeFieldType D_inside = tp.D(*(ig.inside()),inside_local);
288 typename TP::Traits::RangeFieldType D_outside = tp.D(*(ig.outside()),outside_local);
289 typename TP::Traits::RangeFieldType D_avg = 2.0/(1.0/(D_inside+1E-30) + 1.0/(D_outside+1E-30));
292 Dune::FieldVector<DF,IG::dimension>
293 inside_global = ig.inside()->geometry().center();
294 Dune::FieldVector<DF,IG::dimension>
295 outside_global = ig.outside()->geometry().center();
296 inside_global -= outside_global;
297 RF distance = inside_global.two_norm();
301 r_s.accumulate(lfsu_s,0,-(D_avg*(x_n(lfsu_n,0)-x_s(lfsu_s,0))/distance)*face_volume);
305 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
307 const LFSV& lfsv, R& r)
const
310 typedef typename LFSU::Traits::FiniteElementType::
311 Traits::LocalBasisType::Traits::DomainFieldType DF;
312 const int dim = EG::Geometry::dimension;
314 if (!first_stage)
return;
317 const Dune::FieldVector<DF,dim>&
318 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
321 typename TP::Traits::RangeFieldType cellcapacity = tp.c(eg.entity(),inside_local)*eg.geometry().volume();
322 typename TP::Traits::RangeFieldType celldt = cellcapacity/(cellinflux+1E-30);
323 dtmin = std::min(dtmin,celldt);
328 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
330 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
334 typedef typename LFSU::Traits::FiniteElementType::
335 Traits::LocalBasisType::Traits::DomainFieldType DF;
336 typedef typename LFSU::Traits::FiniteElementType::
337 Traits::LocalBasisType::Traits::RangeFieldType RF;
340 const Dune::FieldVector<DF,IG::dimension-1>&
341 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
342 RF face_volume = ig.geometry().volume();
343 Dune::FieldVector<DF,dim> face_center_in_element = ig.geometryInInside().global(face_local);
346 int bc = tp.bc(ig.intersection(),face_local);
351 typename TP::Traits::RangeFieldType j = tp.j(*(ig.inside()),face_center_in_element);
352 r_s.accumulate(lfsu_s,0,j*face_volume);
357 typename TP::Traits::RangeType v(tp.v(*(ig.inside()),face_center_in_element));
360 RF vn = v*ig.centerUnitOuterNormal();
362 cellinflux += vn*face_volume;
366 r_s.accumulate(lfsu_s,0,vn*x_s(lfsu_s,0)*face_volume);
372 typename TP::Traits::RangeFieldType g;
373 if (vn>=0) g=x_s(lfsu_s,0);
else g=tp.g(*(ig.inside()),face_center_in_element);
374 const Dune::FieldVector<DF,IG::dimension>&
375 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
376 typename TP::Traits::RangeFieldType D_inside = tp.D(*(ig.inside()),inside_local);
377 Dune::FieldVector<DF,IG::dimension>
378 inside_global = ig.inside()->geometry().center();
379 Dune::FieldVector<DF,IG::dimension>
380 outside_global = ig.geometry().center();
381 inside_global -= outside_global;
382 RF distance = inside_global.two_norm();
383 r_s.accumulate(lfsu_s,0,(g*vn - D_inside*(g-x_s(lfsu_s,0))/distance)*face_volume);
389 template<
typename EG,
typename LFSV,
typename R>
393 typedef typename LFSV::Traits::FiniteElementType::
394 Traits::LocalBasisType::Traits::DomainFieldType DF;
395 const int dim = EG::Geometry::dimension;
398 const Dune::FieldVector<DF,dim>&
399 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
402 typename TP::Traits::RangeFieldType q = tp.q(eg.entity(),inside_local);
404 r.accumulate(lfsv,0,-q*eg.geometry().volume());
408 void setTime (
typename TP::Traits::RangeFieldType t)
413 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
419 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
426 else first_stage =
false;
435 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const
443 mutable typename TP::Traits::RangeFieldType dtmin;
444 mutable typename TP::Traits::RangeFieldType cellinflux;
449 template<
class T,
class Imp>
456 typename Traits::RangeFieldType
457 c (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
459 return asImp().c(e,x);
463 Imp& asImp () {
return static_cast<Imp &
> (*this);}
464 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
492 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
493 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
496 typedef typename LFSU::Traits::FiniteElementType::
497 Traits::LocalBasisType::Traits::DomainFieldType DF;
500 const int dim = EG::Geometry::dimension;
503 const Dune::FieldVector<DF,dim>&
504 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
507 typename TP::Traits::RangeFieldType c = tp.c(eg.entity(),inside_local);
510 r.accumulate(lfsu,0,c*x(lfsu,0)*eg.geometry().volume());
514 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
519 typedef typename LFSU::Traits::FiniteElementType::
520 Traits::LocalBasisType::Traits::DomainFieldType DF;
523 const int dim = EG::Geometry::dimension;
526 const Dune::FieldVector<DF,dim>&
527 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
530 typename TP::Traits::RangeFieldType c = tp.c(eg.entity(),inside_local);
533 mat.accumulate(lfsu,0,lfsu,0,c*eg.geometry().volume());
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: transportccfv.hh:35
Definition: transportccfv.hh:223
GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: transportccfv.hh:32
Definition: transportccfv.hh:474
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:390
Definition: transportccfv.hh:484
Dune::PDELab::BoundaryGridFunctionTraits< typename T::Traits::GridViewType, int, 1, Dune::FieldVector< int, 1 > > Traits
Definition: transportccfv.hh:127
Traits::RangeFieldType c(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
source term
Definition: transportccfv.hh:457
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const Traits::GridViewType & getGridView() const
Definition: transportccfv.hh:175
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: transportccfv.hh:329
Traits::RangeFieldType q(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
source term
Definition: transportccfv.hh:74
traits class holding function signature, same as in local function
Definition: function.hh:230
BoundaryConditionType_Transport(const typename T::Traits::GridViewType &gv_, const T &t_)
Definition: transportccfv.hh:130
traits class for two phase parameter class
Definition: transportccfv.hh:20
void preStep(typename TP::Traits::RangeFieldType time, typename TP::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: transportccfv.hh:413
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:114
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:493
Definition: defaultimp.hh:96
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: transportccfv.hh:435
const T::Traits::GridViewType & getGridView()
get a reference to the GridView
Definition: transportccfv.hh:141
Dune::PDELab::GridFunctionTraits< typename T::Traits::GridViewType, typename T::Traits::RangeFieldType, 1, Dune::FieldVector< typename T::Traits::RangeFieldType, 1 > > Traits
Definition: transportccfv.hh:162
const IG & ig
Definition: constraints.hh:147
Definition: transportccfv.hh:481
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: transportccfv.hh:47
Definition: transportccfv.hh:216
leaf of a function tree
Definition: function.hh:576
Traits::RangeType v(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
velocity field
Definition: transportccfv.hh:60
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: transportccfv.hh:38
Definition: transportccfv.hh:222
void evaluate(const Dune::PDELab::IntersectionGeometry< I > &ig, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: transportccfv.hh:133
Default flags for all local operators.
Definition: flags.hh:18
int bc(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x) const
boundary condition type function
Definition: transportccfv.hh:86
RF RangeFieldType
Export type for range field.
Definition: transportccfv.hh:41
sparsity pattern generator
Definition: pattern.hh:13
Definition: transportccfv.hh:219
Definition: transportccfv.hh:220
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: transportccfv.hh:515
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: transportccfv.hh:44
T Traits
Definition: transportccfv.hh:56
T Traits
Definition: transportccfv.hh:453
GV GridViewType
the grid view
Definition: transportccfv.hh:23
A local operator for a cell-centered finite volume scheme for the transport equation.
Definition: transportccfv.hh:199
Dune::FieldVector< GV::Grid::ctype, GV::dimension-1 > DomainType
domain type in dim-size coordinates
Definition: function.hh:48
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: transportccfv.hh:419
const I & intersection() const
Definition: geometrywrapper.hh:268
Definition: transportccfv.hh:221
CCFVTemporalOperator(TP &tp_)
Definition: transportccfv.hh:486
traits class holding the function signature, same as in local function
Definition: function.hh:175
Dune::PDELab::BoundaryGridFunctionBase< Traits, BoundaryConditionType_Transport< T > > BaseT
Definition: transportccfv.hh:128
dimension of the domain
Definition: transportccfv.hh:28
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: adaptivity.hh:27
Wrap intersection.
Definition: geometrywrapper.hh:56
Traits::RangeFieldType D(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
scalar diffusion coefficient
Definition: transportccfv.hh:67
leaf of a function tree
Definition: function.hh:596
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:306
Definition: transportccfv.hh:215
static const int dim
Definition: adaptivity.hh:83
Definition: transportccfv.hh:153
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: transportccfv.hh:408
DirichletBoundaryCondition_Transport(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: transportccfv.hh:165
base class for parameter class
Definition: transportccfv.hh:450
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:117
GV::Intersection IntersectionType
Definition: transportccfv.hh:48
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Evaluate the GridFunction at given position.
Definition: transportccfv.hh:168
sparsity pattern generator
Definition: pattern.hh:29
Traits::RangeFieldType j(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Neumann boundary condition.
Definition: transportccfv.hh:100
CCFVSpatialTransportOperator(TP &tp_)
Definition: transportccfv.hh:227
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:234
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: transportccfv.hh:250
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Dirichlet boundary condition on inflow.
Definition: transportccfv.hh:93
R RangeType
range type
Definition: function.hh:60
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Definition: transportccfv.hh:116
Definition: transportccfv.hh:225
Definition: defaultimp.hh:384
void postStage()
to be called once at the end of each stage
Definition: transportccfv.hh:430
base class for parameter class
Definition: transportccfv.hh:53
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: transportccfv.hh:241