2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONCCFV_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSIONCCFV_HH
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/typetraits.hh>
8 #include<dune/geometry/referenceelements.hh>
65 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
66 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
69 typedef typename LFSU::Traits::FiniteElementType::
70 Traits::LocalBasisType::Traits::DomainFieldType DF;
71 typedef typename LFSU::Traits::FiniteElementType::
72 Traits::LocalBasisType::Traits::RangeFieldType RF;
73 typedef typename LFSU::Traits::FiniteElementType::
74 Traits::LocalBasisType::Traits::JacobianType JacobianType;
75 typedef typename LFSU::Traits::FiniteElementType::
76 Traits::LocalBasisType::Traits::RangeType RangeType;
79 const int dim = EG::Geometry::mydimension;
82 const Dune::FieldVector<DF,dim>
83 inside_local(Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0));
86 typename TP::Traits::RangeFieldType c = param.c(eg.entity(),inside_local);
89 r.accumulate(lfsu,0,(c*x(lfsu,0))*eg.geometry().volume());
93 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
98 typedef typename LFSU::Traits::FiniteElementType::
99 Traits::LocalBasisType::Traits::DomainFieldType DF;
100 typedef typename LFSU::Traits::FiniteElementType::
101 Traits::LocalBasisType::Traits::RangeFieldType RF;
102 typedef typename LFSU::Traits::FiniteElementType::
103 Traits::LocalBasisType::Traits::JacobianType JacobianType;
104 typedef typename LFSU::Traits::FiniteElementType::
105 Traits::LocalBasisType::Traits::RangeType RangeType;
106 typedef typename LFSU::Traits::SizeType size_type;
109 const int dim = EG::Geometry::mydimension;
112 const Dune::FieldVector<DF,dim>
113 inside_local(Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0));
116 typename TP::Traits::RangeFieldType c = param.c(eg.entity(),inside_local);
119 mat.accumulate(lfsu,0,lfsu,0,c*eg.geometry().volume());
124 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
126 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
127 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
128 R& r_s, R& r_n)
const
131 typedef typename LFSU::Traits::FiniteElementType::
132 Traits::LocalBasisType::Traits::DomainFieldType DF;
133 typedef typename LFSU::Traits::FiniteElementType::
134 Traits::LocalBasisType::Traits::RangeFieldType RF;
135 const int dim = IG::dimension;
138 const Dune::FieldVector<DF,dim-1>
139 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
142 RF face_volume = ig.geometry().integrationElement(face_local)
143 *Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).volume();
145 auto cell_inside = ig.inside();
146 auto cell_outside = ig.outside();
149 const Dune::FieldVector<DF,dim>
150 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(cell_inside.type()).position(0,0);
151 const Dune::FieldVector<DF,dim>
152 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(cell_outside.type()).position(0,0);
155 typename TP::Traits::PermTensorType tensor_inside;
156 tensor_inside = param.A(cell_inside,inside_local);
157 typename TP::Traits::PermTensorType tensor_outside;
158 tensor_outside = param.A(cell_outside,outside_local);
159 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
160 Dune::FieldVector<RF,dim> An_F;
161 tensor_inside.mv(n_F,An_F);
162 RF k_inside = n_F*An_F;
163 tensor_outside.mv(n_F,An_F);
164 RF k_outside = n_F*An_F;
165 RF k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
168 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(face_local);
169 typename TP::Traits::RangeType b = param.b(cell_inside,iplocal_s);
172 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
175 Dune::FieldVector<DF,IG::dimension>
176 inside_global = cell_inside.geometry().global(inside_local);
177 Dune::FieldVector<DF,IG::dimension>
178 outside_global = cell_outside.geometry().global(outside_local);
181 inside_global -= outside_global;
182 RF distance = inside_global.two_norm();
185 r_s.accumulate(lfsu_s,0,(u_upwind*vn)*face_volume+k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
186 r_n.accumulate(lfsu_n,0,-(u_upwind*vn)*face_volume-k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
189 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
191 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
192 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
193 M& mat_ss, M& mat_sn,
194 M& mat_ns, M& mat_nn)
const
197 typedef typename LFSU::Traits::FiniteElementType::
198 Traits::LocalBasisType::Traits::DomainFieldType DF;
199 typedef typename LFSU::Traits::FiniteElementType::
200 Traits::LocalBasisType::Traits::RangeFieldType RF;
201 const int dim = IG::dimension;
204 const Dune::FieldVector<DF,dim-1>
205 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
208 RF face_volume = ig.geometry().integrationElement(face_local)
209 *Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).volume();
211 auto cell_inside = ig.inside();
212 auto cell_outside = ig.outside();
215 const Dune::FieldVector<DF,dim>
216 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(cell_inside.type()).position(0,0);
217 const Dune::FieldVector<DF,dim>
218 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(cell_outside.type()).position(0,0);
221 typename TP::Traits::PermTensorType tensor_inside;
222 tensor_inside = param.A(cell_inside,inside_local);
223 typename TP::Traits::PermTensorType tensor_outside;
224 tensor_outside = param.A(cell_outside,outside_local);
225 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
226 Dune::FieldVector<RF,dim> An_F;
227 tensor_inside.mv(n_F,An_F);
228 RF k_inside = n_F*An_F;
229 tensor_outside.mv(n_F,An_F);
230 RF k_outside = n_F*An_F;
231 RF k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
234 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(face_local);
235 typename TP::Traits::RangeType b = param.b(cell_inside,iplocal_s);
239 Dune::FieldVector<DF,IG::dimension>
240 inside_global = cell_inside.geometry().global(inside_local);
241 Dune::FieldVector<DF,IG::dimension>
242 outside_global = cell_outside.geometry().global(outside_local);
245 inside_global -= outside_global;
246 RF distance = inside_global.two_norm();
249 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, k_avg*face_volume/distance );
250 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -k_avg*face_volume/distance );
251 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, -k_avg*face_volume/distance );
252 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, k_avg*face_volume/distance );
255 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, vn*face_volume );
256 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -vn*face_volume );
260 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, vn*face_volume );
261 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, -vn*face_volume );
269 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
271 const LFSV& lfsv, R& r)
const
274 typedef typename LFSU::Traits::FiniteElementType::
275 Traits::LocalBasisType::Traits::DomainFieldType DF;
276 const int dim = EG::Geometry::mydimension;
278 if (!first_stage)
return;
281 const Dune::FieldVector<DF,dim>&
282 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
285 typename TP::Traits::RangeFieldType cellcapacity = param.d(eg.entity(),inside_local)*eg.geometry().volume();
286 typename TP::Traits::RangeFieldType celldt = cellcapacity/(cellinflux+1E-30);
287 dtmin = std::min(dtmin,celldt);
294 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
296 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
300 typedef typename LFSU::Traits::FiniteElementType::
301 Traits::LocalBasisType::Traits::DomainFieldType DF;
302 typedef typename LFSU::Traits::FiniteElementType::
303 Traits::LocalBasisType::Traits::RangeFieldType RF;
304 const int dim = IG::dimension;
307 const Dune::FieldVector<DF,dim-1>
308 face_local = Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).position(0,0);
311 RF face_volume = ig.geometry().integrationElement(face_local)
312 *Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).volume();
314 auto cell_inside = ig.inside();
317 const Dune::FieldVector<DF,dim>
318 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(cell_inside.type()).position(0,0);
322 bctype = param.bctype(ig.intersection(),face_local);
328 Dune::FieldVector<DF,dim> inside_global = cell_inside.geometry().global(inside_local);
329 Dune::FieldVector<DF,dim> outside_global = ig.geometry().global(face_local);
330 inside_global -= outside_global;
331 RF distance = inside_global.two_norm();
334 typename TP::Traits::PermTensorType tensor_inside;
335 tensor_inside = param.A(cell_inside,inside_local);
336 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
337 Dune::FieldVector<RF,dim> An_F;
338 tensor_inside.mv(n_F,An_F);
339 RF k_inside = n_F*An_F;
342 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(face_local);
343 RF g = param.g(cell_inside,iplocal_s);
346 typename TP::Traits::RangeType b = param.b(cell_inside,iplocal_s);
347 const Dune::FieldVector<DF,dim> n = ig.centerUnitOuterNormal();
350 r_s.accumulate(lfsu_s,0,(b*n)*g*face_volume + k_inside*(x_s(lfsu_s,0)-g)*face_volume/distance);
361 typename TP::Traits::RangeFieldType j = param.j(ig.intersection(),face_local);
364 r_s.accumulate(lfsu_s,0,j*face_volume);
372 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(face_local);
373 typename TP::Traits::RangeType b = param.b(cell_inside,iplocal_s);
374 const Dune::FieldVector<DF,dim> n = ig.centerUnitOuterNormal();
377 typename TP::Traits::RangeFieldType o = param.o(ig.intersection(),face_local);
380 r_s.accumulate(lfsu_s,0,((b*n)*x_s(lfsu_s,0) + o)*face_volume);
386 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
388 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
392 typedef typename LFSU::Traits::FiniteElementType::
393 Traits::LocalBasisType::Traits::DomainFieldType DF;
394 typedef typename LFSU::Traits::FiniteElementType::
395 Traits::LocalBasisType::Traits::RangeFieldType RF;
396 const int dim = IG::dimension;
399 const Dune::FieldVector<DF,dim-1>
400 face_local = Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).position(0,0);
403 RF face_volume = ig.geometry().integrationElement(face_local)
404 *Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).volume();
406 auto cell_inside = ig.inside();
409 const Dune::FieldVector<DF,dim>
410 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(cell_inside.type()).position(0,0);
414 bctype = param.bctype(ig.intersection(),face_local);
421 Dune::FieldVector<DF,dim> inside_global = cell_inside.geometry().global(inside_local);
422 Dune::FieldVector<DF,dim> outside_global = ig.geometry().global(face_local);
423 inside_global -= outside_global;
424 RF distance = inside_global.two_norm();
427 typename TP::Traits::PermTensorType tensor_inside;
428 tensor_inside = param.A(cell_inside,inside_local);
429 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
430 Dune::FieldVector<RF,dim> An_F;
431 tensor_inside.mv(n_F,An_F);
432 RF k_inside = n_F*An_F;
435 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, k_inside*face_volume/distance );
443 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(face_local);
444 typename TP::Traits::RangeType b = param.b(cell_inside,iplocal_s);
445 const Dune::FieldVector<DF,dim> n = ig.centerUnitOuterNormal();
448 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, (b*n)*face_volume );
455 template<
typename EG,
typename LFSV,
typename R>
459 typedef typename LFSV::Traits::FiniteElementType::
460 Traits::LocalBasisType::Traits::DomainFieldType DF;
461 const int dim = EG::Geometry::mydimension;
464 const Dune::FieldVector<DF,dim>&
465 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
468 typename TP::Traits::RangeFieldType f = param.f(eg.entity(),inside_local);
470 r.accumulate(lfsv,0,-f*eg.geometry().volume());
474 void setTime (
typename TP::Traits::RangeFieldType t)
480 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
486 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
493 else first_stage =
false;
502 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const
510 mutable typename TP::Traits::RangeFieldType dtmin;
511 mutable typename TP::Traits::RangeFieldType cellinflux;
544 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
545 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
548 typedef typename LFSU::Traits::FiniteElementType::
549 Traits::LocalBasisType::Traits::DomainFieldType DF;
552 const int dim = EG::Geometry::mydimension;
555 const Dune::FieldVector<DF,dim>&
556 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
558 typename TP::Traits::RangeFieldType capacity = param.d(eg.entity(),inside_local);
561 r.accumulate(lfsu,0,capacity*x(lfsu,0)*eg.geometry().volume());
565 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
570 typedef typename LFSU::Traits::FiniteElementType::
571 Traits::LocalBasisType::Traits::DomainFieldType DF;
574 const int dim = EG::Geometry::mydimension;
577 const Dune::FieldVector<DF,dim>&
578 inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
581 typename TP::Traits::RangeFieldType capacity = param.d(eg.entity(),inside_local);
584 mat.accumulate(lfsu,0,lfsu,0,capacity*eg.geometry().volume());
Definition: convectiondiffusionccfv.hh:36
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:456
void jacobian_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, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: convectiondiffusionccfv.hh:190
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: convectiondiffusionccfv.hh:387
Definition: convectiondiffusionccfv.hh:55
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: convectiondiffusionccfv.hh:502
Definition: convectiondiffusionccfv.hh:56
ConvectionDiffusionCCFV(TP ¶m_)
Definition: convectiondiffusionccfv.hh:61
const IG & ig
Definition: constraints.hh:147
Type
Definition: convectiondiffusionparameter.hh:67
Default flags for all local operators.
Definition: flags.hh:18
Definition: convectiondiffusionparameter.hh:67
Definition: convectiondiffusionccfv.hh:50
Definition: convectiondiffusionccfv.hh:58
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusionparameter.hh:67
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:566
Definition: convectiondiffusionparameter.hh:67
ConvectionDiffusionCCFVTemporalOperator(TP ¶m_)
Definition: convectiondiffusionccfv.hh:538
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionccfv.hh:295
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:545
void postStage()
to be called once at the end of each stage
Definition: convectiondiffusionccfv.hh:497
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:66
Definition: convectiondiffusionccfv.hh:533
Definition: adaptivity.hh:27
Definition: convectiondiffusionccfv.hh:59
Definition: convectiondiffusionccfv.hh:524
static const int dim
Definition: adaptivity.hh:83
Definition: convectiondiffusionccfv.hh:54
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: convectiondiffusionccfv.hh:474
Definition: convectiondiffusionccfv.hh:536
Definition: convectiondiffusionccfv.hh:57
sparsity pattern generator
Definition: pattern.hh:29
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:94
Definition: convectiondiffusionccfv.hh:51
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: convectiondiffusionccfv.hh:125
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:270
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: convectiondiffusionccfv.hh:486
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void preStep(typename TP::Traits::RangeFieldType time, typename TP::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: convectiondiffusionccfv.hh:480