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;
197 template<
typename TP>
210 enum { dim = TP::Traits::GridViewType::dimension };
232 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
233 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
239 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
248 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
250 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
251 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
252 R& r_s, R& r_n)
const
255 typedef typename LFSU::Traits::FiniteElementType::
256 Traits::LocalBasisType::Traits::DomainFieldType DF;
257 typedef typename LFSU::Traits::FiniteElementType::
258 Traits::LocalBasisType::Traits::RangeFieldType RF;
261 const Dune::FieldVector<DF,IG::dimension-1>&
262 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
263 RF face_volume = ig.geometry().volume();
266 Dune::FieldVector<DF,IG::dimension> face_center_in_element = ig.geometryInInside().global(face_local);
269 typename TP::Traits::RangeType v(tp.v(*(ig.inside()),face_center_in_element));
272 RF vn = v*ig.centerUnitOuterNormal();
276 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
277 r_s.accumulate(lfsu_s,0,(u_upwind*vn)*face_volume);
279 cellinflux += vn*face_volume;
282 const Dune::FieldVector<DF,IG::dimension>&
283 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
284 const Dune::FieldVector<DF,IG::dimension>&
285 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
286 typename TP::Traits::RangeFieldType D_inside = tp.D(*(ig.inside()),inside_local);
287 typename TP::Traits::RangeFieldType D_outside = tp.D(*(ig.outside()),outside_local);
288 typename TP::Traits::RangeFieldType D_avg = 2.0/(1.0/(D_inside+1E-30) + 1.0/(D_outside+1E-30));
291 Dune::FieldVector<DF,IG::dimension>
292 inside_global = ig.inside()->geometry().center();
293 Dune::FieldVector<DF,IG::dimension>
294 outside_global = ig.outside()->geometry().center();
295 inside_global -= outside_global;
296 RF distance = inside_global.two_norm();
300 r_s.accumulate(lfsu_s,0,-(D_avg*(x_n(lfsu_n,0)-x_s(lfsu_s,0))/distance)*face_volume);
304 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
306 const LFSV& lfsv, R& r)
const
309 typedef typename LFSU::Traits::FiniteElementType::
310 Traits::LocalBasisType::Traits::DomainFieldType DF;
311 const int dim = EG::Geometry::dimension;
313 if (!first_stage)
return;
316 const Dune::FieldVector<DF,dim>&
317 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
320 typename TP::Traits::RangeFieldType cellcapacity = tp.c(eg.entity(),inside_local)*eg.geometry().volume();
321 typename TP::Traits::RangeFieldType celldt = cellcapacity/(cellinflux+1E-30);
322 dtmin = std::min(dtmin,celldt);
327 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
329 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
333 typedef typename LFSU::Traits::FiniteElementType::
334 Traits::LocalBasisType::Traits::DomainFieldType DF;
335 typedef typename LFSU::Traits::FiniteElementType::
336 Traits::LocalBasisType::Traits::RangeFieldType RF;
339 const Dune::FieldVector<DF,IG::dimension-1>&
340 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
341 RF face_volume = ig.geometry().volume();
342 Dune::FieldVector<DF,dim> face_center_in_element = ig.geometryInInside().global(face_local);
345 int bc = tp.bc(ig.intersection(),face_local);
350 typename TP::Traits::RangeFieldType j = tp.j(*(ig.inside()),face_center_in_element);
351 r_s.accumulate(lfsu_s,0,j*face_volume);
356 typename TP::Traits::RangeType v(tp.v(*(ig.inside()),face_center_in_element));
359 RF vn = v*ig.centerUnitOuterNormal();
361 cellinflux += vn*face_volume;
365 r_s.accumulate(lfsu_s,0,vn*x_s(lfsu_s,0)*face_volume);
371 typename TP::Traits::RangeFieldType g;
372 if (vn>=0) g=x_s(lfsu_s,0);
else g=tp.g(*(ig.inside()),face_center_in_element);
373 const Dune::FieldVector<DF,IG::dimension>&
374 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
375 typename TP::Traits::RangeFieldType D_inside = tp.D(*(ig.inside()),inside_local);
376 Dune::FieldVector<DF,IG::dimension>
377 inside_global = ig.inside()->geometry().center();
378 Dune::FieldVector<DF,IG::dimension>
379 outside_global = ig.geometry().center();
380 inside_global -= outside_global;
381 RF distance = inside_global.two_norm();
382 r_s.accumulate(lfsu_s,0,(g*vn - D_inside*(g-x_s(lfsu_s,0))/distance)*face_volume);
388 template<
typename EG,
typename LFSV,
typename R>
392 typedef typename LFSV::Traits::FiniteElementType::
393 Traits::LocalBasisType::Traits::DomainFieldType DF;
394 const int dim = EG::Geometry::dimension;
397 const Dune::FieldVector<DF,dim>&
398 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
401 typename TP::Traits::RangeFieldType q = tp.q(eg.entity(),inside_local);
403 r.accumulate(lfsv,0,-q*eg.geometry().volume());
407 void setTime (
typename TP::Traits::RangeFieldType t)
412 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
418 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
425 else first_stage =
false;
434 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const
442 mutable typename TP::Traits::RangeFieldType dtmin;
443 mutable typename TP::Traits::RangeFieldType cellinflux;
448 template<
class T,
class Imp>
455 typename Traits::RangeFieldType
456 c (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
458 return asImp().c(e,x);
462 Imp& asImp () {
return static_cast<Imp &
> (*this);}
463 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
491 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
492 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
495 typedef typename LFSU::Traits::FiniteElementType::
496 Traits::LocalBasisType::Traits::DomainFieldType DF;
499 const int dim = EG::Geometry::dimension;
502 const Dune::FieldVector<DF,dim>&
503 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
506 typename TP::Traits::RangeFieldType c = tp.c(eg.entity(),inside_local);
509 r.accumulate(lfsu,0,c*x(lfsu,0)*eg.geometry().volume());
513 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
518 typedef typename LFSU::Traits::FiniteElementType::
519 Traits::LocalBasisType::Traits::DomainFieldType DF;
522 const int dim = EG::Geometry::dimension;
525 const Dune::FieldVector<DF,dim>&
526 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
529 typename TP::Traits::RangeFieldType c = tp.c(eg.entity(),inside_local);
532 mat.accumulate(lfsu,0,lfsu,0,c*eg.geometry().volume());
Definition: transportccfv.hh:219
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:30
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
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:249
Definition: transportccfv.hh:221
GV GridViewType
the grid view
Definition: transportccfv.hh:23
Traits::RangeFieldType c(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
source term
Definition: transportccfv.hh:456
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:305
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
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
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:118
T Traits
Definition: transportccfv.hh:56
Definition: transportccfv.hh:198
Definition: transportccfv.hh:220
Definition: defaultimp.hh:96
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: transportccfv.hh:47
Traits::RangeFieldType D(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
scalar diffusion coefficient
Definition: transportccfv.hh:67
Definition: transportccfv.hh:214
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
DirichletBoundaryCondition_Transport(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: transportccfv.hh:165
Definition: transportccfv.hh:218
static const int dim
Definition: adaptivity.hh:82
const I & intersection() const
Definition: geometrywrapper.hh:239
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: transportccfv.hh:514
base class for parameter class
Definition: transportccfv.hh:53
Dune::FieldVector< GV::Grid::ctype, GV::dimension-1 > DomainType
domain type in dim-size coordinates
Definition: function.hh:49
Definition: transportccfv.hh:224
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:412
BoundaryConditionType_Transport(const typename T::Traits::GridViewType &gv_, const T &t_)
Definition: transportccfv.hh:130
Definition: transportccfv.hh:473
Traits::RangeFieldType j(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Neumann boundary condition.
Definition: transportccfv.hh:100
const T::Traits::GridViewType & getGridView()
get a reference to the GridView
Definition: transportccfv.hh:141
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:328
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: transportccfv.hh:407
Definition: transportccfv.hh:480
Definition: transportccfv.hh:215
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: transportccfv.hh:38
Definition: transportccfv.hh:153
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:492
void evaluate(const Dune::PDELab::IntersectionGeometry< I > &ig, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: transportccfv.hh:133
Traits::RangeType v(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
velocity field
Definition: transportccfv.hh:60
base class for parameter class
Definition: transportccfv.hh:449
R RangeType
range type
Definition: function.hh:61
const IG & ig
Definition: common/constraints.hh:146
Definition: defaultimp.hh:384
leaf of a function tree
Definition: function.hh:577
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: transportccfv.hh:44
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: transportccfv.hh:418
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: transportccfv.hh:240
leaf of a function tree
Definition: function.hh:597
Dune::PDELab::BoundaryGridFunctionBase< Traits, BoundaryConditionType_Transport< T > > BaseT
Definition: transportccfv.hh:128
const Traits::GridViewType & getGridView() const
Definition: transportccfv.hh:175
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
CCFVSpatialTransportOperator(TP &tp_)
Definition: transportccfv.hh:226
traits class holding the function signature, same as in local function
Definition: function.hh:176
traits class for two phase parameter class
Definition: transportccfv.hh:20
traits class holding function signature, same as in local function
Definition: function.hh:231
sparsity pattern generator
Definition: pattern.hh:14
Dune::PDELab::BoundaryGridFunctionTraits< typename T::Traits::GridViewType, int, 1, Dune::FieldVector< int, 1 > > Traits
Definition: transportccfv.hh:127
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:115
void postStage()
to be called once at the end of each stage
Definition: transportccfv.hh:429
const EG & eg
Definition: common/constraints.hh:277
GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: transportccfv.hh:32
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
GV::Intersection IntersectionType
Definition: transportccfv.hh:48
Definition: transportccfv.hh:116
Definition: transportccfv.hh:483
Traits::RangeFieldType q(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
source term
Definition: transportccfv.hh:74
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:389
const E & e
Definition: interpolate.hh:172
CCFVTemporalOperator(TP &tp_)
Definition: transportccfv.hh:485
Definition: transportccfv.hh:222
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: transportccfv.hh:35
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: transportccfv.hh:233
Dune::PDELab::GridFunctionTraits< typename T::Traits::GridViewType, typename T::Traits::RangeFieldType, 1, Dune::FieldVector< typename T::Traits::RangeFieldType, 1 > > Traits
Definition: transportccfv.hh:162
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Dirichlet boundary condition on inflow.
Definition: transportccfv.hh:93
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: transportccfv.hh:434
T Traits
Definition: transportccfv.hh:452
Wrap intersection.
Definition: geometrywrapper.hh:56
dimension of the domain
Definition: transportccfv.hh:28