2 #ifndef DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_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> 67 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
68 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 71 auto geo = eg.geometry();
73 auto local_inside = ref_el.position(0,0);
76 auto c = param.c(eg.entity(),local_inside);
79 r.accumulate(lfsu,0,(c*x(lfsu,0))*geo.volume());
83 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
84 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
88 auto geo = eg.geometry();
90 auto local_inside = ref_el.position(0,0);
93 auto c = param.c(eg.entity(),local_inside);
96 mat.accumulate(lfsu,0,lfsu,0,c*geo.volume());
101 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
103 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
104 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
105 R& r_s, R& r_n)
const 108 using RF =
typename LFSU::Traits::FiniteElementType::
109 Traits::LocalBasisType::Traits::RangeFieldType;
112 const auto dim = IG::dimension;
115 auto cell_inside = ig.inside();
116 auto cell_outside = ig.outside();
119 auto geo = ig.geometry();
120 auto geo_inside = cell_inside.geometry();
121 auto geo_outside = cell_outside.geometry();
124 auto geo_in_inside = ig.geometryInInside();
128 auto face_local = ref_el.position(0,0);
131 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
136 auto local_inside = ref_el_inside.position(0,0);
137 auto local_outside = ref_el_outside.position(0,0);
140 auto tensor_inside = param.A(cell_inside,local_inside);
141 auto tensor_outside = param.A(cell_outside,local_outside);
142 auto n_F = ig.centerUnitOuterNormal();
143 Dune::FieldVector<RF,dim> An_F;
144 tensor_inside.mv(n_F,An_F);
145 auto k_inside = n_F*An_F;
146 tensor_outside.mv(n_F,An_F);
147 auto k_outside = n_F*An_F;
148 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
151 auto iplocal_s = geo_in_inside.global(face_local);
152 auto b = param.b(cell_inside,iplocal_s);
155 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
158 auto global_inside = geo_inside.global(local_inside);
159 auto global_outside = geo_outside.global(local_outside);
162 global_inside -= global_outside;
163 auto distance = global_inside.two_norm();
166 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);
167 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);
170 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
172 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
173 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
174 M& mat_ss, M& mat_sn,
175 M& mat_ns, M& mat_nn)
const 178 using RF =
typename LFSU::Traits::FiniteElementType::
179 Traits::LocalBasisType::Traits::RangeFieldType;
182 const auto dim = IG::dimension;
185 auto cell_inside = ig.inside();
186 auto cell_outside = ig.outside();
189 auto geo = ig.geometry();
190 auto geo_inside = cell_inside.geometry();
191 auto geo_outside = cell_outside.geometry();
194 auto geo_in_inside = ig.geometryInInside();
198 auto face_local = ref_el.position(0,0);
201 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
206 auto local_inside = ref_el_inside.position(0,0);
207 auto local_outside = ref_el_outside.position(0,0);
210 auto tensor_inside = param.A(cell_inside,local_inside);
211 auto tensor_outside = param.A(cell_outside,local_outside);
212 auto n_F = ig.centerUnitOuterNormal();
213 Dune::FieldVector<RF,dim> An_F;
214 tensor_inside.mv(n_F,An_F);
215 auto k_inside = n_F*An_F;
216 tensor_outside.mv(n_F,An_F);
217 auto k_outside = n_F*An_F;
218 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
221 auto iplocal_s = geo_in_inside.global(face_local);
222 auto b = param.b(cell_inside,iplocal_s);
226 auto global_inside = geo_inside.global(local_inside);
227 auto global_outside = geo_outside.global(local_outside);
230 global_inside -= global_outside;
231 auto distance = global_inside.two_norm();
234 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, k_avg*face_volume/distance );
235 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -k_avg*face_volume/distance );
236 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, -k_avg*face_volume/distance );
237 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, k_avg*face_volume/distance );
240 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, vn*face_volume );
241 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -vn*face_volume );
245 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, vn*face_volume );
246 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, -vn*face_volume );
252 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
254 const LFSV& lfsv, R& r)
const 256 if (!first_stage)
return;
259 auto geo = eg.geometry();
261 auto local_inside = ref_el.position(0,0);
264 auto cellcapacity = param.d(eg.entity(),local_inside)*geo.volume();
265 auto celldt = cellcapacity/(cellinflux+1E-30);
266 dtmin = std::min(dtmin,celldt);
272 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
274 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
278 using RF =
typename LFSU::Traits::FiniteElementType::
279 Traits::LocalBasisType::Traits::RangeFieldType;
282 const auto dim = IG::dimension;
285 auto cell_inside = ig.inside();
288 auto geo = ig.geometry();
289 auto geo_inside = cell_inside.geometry();
292 auto geo_in_inside = ig.geometryInInside();
296 auto face_local = ref_el.position(0,0);
299 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
303 auto local_inside = ref_el_inside.position(0,0);
306 auto bctype = param.bctype(ig.intersection(),face_local);
312 auto global_inside = geo_inside.global(local_inside);
313 auto global_outside = geo.global(face_local);
314 global_inside -= global_outside;
315 auto distance = global_inside.two_norm();
318 auto tensor_inside = param.A(cell_inside,local_inside);
319 auto n_F = ig.centerUnitOuterNormal();
320 Dune::FieldVector<RF,dim> An_F;
321 tensor_inside.mv(n_F,An_F);
322 auto k_inside = n_F*An_F;
325 auto iplocal_s = geo_in_inside.global(face_local);
326 auto g = param.g(cell_inside,iplocal_s);
329 auto b = param.b(cell_inside,iplocal_s);
330 auto n = ig.centerUnitOuterNormal();
333 r_s.accumulate(lfsu_s,0,(b*n)*g*face_volume + k_inside*(x_s(lfsu_s,0)-g)*face_volume/distance);
344 auto j = param.j(ig.intersection(),face_local);
347 r_s.accumulate(lfsu_s,0,j*face_volume);
355 auto iplocal_s = geo_in_inside.global(face_local);
356 auto b = param.b(cell_inside,iplocal_s);
357 auto n = ig.centerUnitOuterNormal();
360 auto o = param.o(ig.intersection(),face_local);
363 r_s.accumulate(lfsu_s,0,((b*n)*x_s(lfsu_s,0) + o)*face_volume);
369 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
371 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
375 using RF =
typename LFSU::Traits::FiniteElementType::
376 Traits::LocalBasisType::Traits::RangeFieldType;
379 const auto dim = IG::dimension;
382 auto cell_inside = ig.inside();
385 auto geo = ig.geometry();
386 auto geo_inside = cell_inside.geometry();
389 auto geo_in_inside = ig.geometryInInside();
393 auto face_local = ref_el.position(0,0);
396 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
400 auto local_inside = ref_el_inside.position(0,0);
403 auto bctype = param.bctype(ig.intersection(),face_local);
409 auto global_inside = geo_inside.global(local_inside);
410 auto global_outside = geo.global(face_local);
411 global_inside -= global_outside;
412 auto distance = global_inside.two_norm();
415 auto tensor_inside = param.A(cell_inside,local_inside);
416 auto n_F = ig.centerUnitOuterNormal();
417 Dune::FieldVector<RF,dim> An_F;
418 tensor_inside.mv(n_F,An_F);
419 auto k_inside = n_F*An_F;
422 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, k_inside*face_volume/distance );
430 auto iplocal_s = geo_in_inside.global(face_local);
431 auto b = param.b(cell_inside,iplocal_s);
432 auto n = ig.centerUnitOuterNormal();
435 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, (b*n)*face_volume );
442 template<
typename EG,
typename LFSV,
typename R>
446 auto geo = eg.geometry();
448 auto local_inside = ref_el.position(0,0);
451 auto f = param.f(eg.entity(),local_inside);
453 r.accumulate(lfsv,0,-f*geo.volume());
457 void setTime (
typename TP::Traits::RangeFieldType t)
463 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
469 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
476 else first_stage =
false;
485 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const 493 mutable typename TP::Traits::RangeFieldType dtmin;
494 mutable typename TP::Traits::RangeFieldType cellinflux;
527 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
528 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 531 auto geo = eg.geometry();
533 auto local_inside = ref_el.position(0,0);
536 auto capacity = param.d(eg.entity(),local_inside);
539 r.accumulate(lfsu,0,capacity*x(lfsu,0)*geo.volume());
543 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
548 auto geo = eg.geometry();
550 auto local_inside = ref_el.position(0,0);
553 auto capacity = param.d(eg.entity(),local_inside);
556 mat.accumulate(lfsu,0,lfsu,0,capacity*geo.volume());
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:68
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:528
Definition: convectiondiffusionccfv.hh:52
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:273
Definition: convectiondiffusionccfv.hh:59
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusionccfv.hh:58
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:253
const IG & ig
Definition: constraints.hh:148
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: convectiondiffusionccfv.hh:469
static const int dim
Definition: adaptivity.hh:83
Definition: convectiondiffusionccfv.hh:60
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Definition: adaptivity.hh:27
Definition: convectiondiffusionparameter.hh:65
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: convectiondiffusionccfv.hh:457
Type
Definition: convectiondiffusionparameter.hh:65
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:84
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:102
sparsity pattern generator
Definition: pattern.hh:29
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:463
Default flags for all local operators.
Definition: flags.hh:18
ConvectionDiffusionCCFVTemporalOperator(TP ¶m_)
Definition: convectiondiffusionccfv.hh:521
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:443
Definition: convectiondiffusionparameter.hh:65
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:171
Definition: convectiondiffusionccfv.hh:61
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: convectiondiffusionccfv.hh:507
Definition: convectiondiffusionparameter.hh:65
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: convectiondiffusionccfv.hh:485
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:544
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:370
ConvectionDiffusionCCFV(TP ¶m_)
Definition: convectiondiffusionccfv.hh:63
Definition: convectiondiffusionccfv.hh:38
Definition: convectiondiffusionccfv.hh:56
Definition: convectiondiffusionccfv.hh:53
Definition: convectiondiffusionccfv.hh:57
void postStage()
to be called once at the end of each stage
Definition: convectiondiffusionccfv.hh:480