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> 66 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
67 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 70 const auto& cell = eg.entity();
73 auto geo = eg.geometry();
75 auto local_inside = ref_el.position(0,0);
78 auto c = param.c(cell,local_inside);
81 r.accumulate(lfsu,0,(c*x(lfsu,0))*geo.volume());
85 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
92 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
93 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
97 const auto& cell = eg.entity();
100 auto geo = eg.geometry();
102 auto local_inside = ref_el.position(0,0);
105 auto c = param.c(cell,local_inside);
108 mat.accumulate(lfsu,0,lfsu,0,c*geo.volume());
113 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
115 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
116 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
117 R& r_s, R& r_n)
const 120 using RF =
typename LFSU::Traits::FiniteElementType::
121 Traits::LocalBasisType::Traits::RangeFieldType;
124 const auto dim = IG::dimension;
127 auto cell_inside = ig.inside();
128 auto cell_outside = ig.outside();
131 auto geo = ig.geometry();
132 auto geo_inside = cell_inside.geometry();
133 auto geo_outside = cell_outside.geometry();
136 auto geo_in_inside = ig.geometryInInside();
140 auto face_local = ref_el.position(0,0);
143 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
148 auto local_inside = ref_el_inside.position(0,0);
149 auto local_outside = ref_el_outside.position(0,0);
152 auto tensor_inside = param.A(cell_inside,local_inside);
153 auto tensor_outside = param.A(cell_outside,local_outside);
154 auto n_F = ig.centerUnitOuterNormal();
155 Dune::FieldVector<RF,dim> An_F;
156 tensor_inside.mv(n_F,An_F);
157 auto k_inside = n_F*An_F;
158 tensor_outside.mv(n_F,An_F);
159 auto k_outside = n_F*An_F;
160 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
163 auto iplocal_s = geo_in_inside.global(face_local);
164 auto b = param.b(cell_inside,iplocal_s);
167 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
170 auto global_inside = geo_inside.global(local_inside);
171 auto global_outside = geo_outside.global(local_outside);
174 global_inside -= global_outside;
175 auto distance = global_inside.two_norm();
178 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);
179 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);
183 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
185 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
186 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
187 Y& y_s, Y& y_n)
const 192 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
194 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
195 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
196 M& mat_ss, M& mat_sn,
197 M& mat_ns, M& mat_nn)
const 200 using RF =
typename LFSU::Traits::FiniteElementType::
201 Traits::LocalBasisType::Traits::RangeFieldType;
204 const auto dim = IG::dimension;
207 auto cell_inside = ig.inside();
208 auto cell_outside = ig.outside();
211 auto geo = ig.geometry();
212 auto geo_inside = cell_inside.geometry();
213 auto geo_outside = cell_outside.geometry();
216 auto geo_in_inside = ig.geometryInInside();
220 auto face_local = ref_el.position(0,0);
223 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
228 auto local_inside = ref_el_inside.position(0,0);
229 auto local_outside = ref_el_outside.position(0,0);
232 auto tensor_inside = param.A(cell_inside,local_inside);
233 auto tensor_outside = param.A(cell_outside,local_outside);
234 auto n_F = ig.centerUnitOuterNormal();
235 Dune::FieldVector<RF,dim> An_F;
236 tensor_inside.mv(n_F,An_F);
237 auto k_inside = n_F*An_F;
238 tensor_outside.mv(n_F,An_F);
239 auto k_outside = n_F*An_F;
240 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
243 auto iplocal_s = geo_in_inside.global(face_local);
244 auto b = param.b(cell_inside,iplocal_s);
248 auto global_inside = geo_inside.global(local_inside);
249 auto global_outside = geo_outside.global(local_outside);
252 global_inside -= global_outside;
253 auto distance = global_inside.two_norm();
256 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, k_avg*face_volume/distance );
257 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -k_avg*face_volume/distance );
258 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, -k_avg*face_volume/distance );
259 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, k_avg*face_volume/distance );
262 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, vn*face_volume );
263 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -vn*face_volume );
267 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, vn*face_volume );
268 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, -vn*face_volume );
274 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
276 const LFSV& lfsv, R& r)
const 278 if (not first_stage)
return;
281 const auto& cell = eg.entity();
284 auto geo = eg.geometry();
286 auto local_inside = ref_el.position(0,0);
289 auto cellcapacity = param.d(cell,local_inside)*geo.volume();
290 auto celldt = cellcapacity/(cellinflux+1E-30);
291 dtmin = std::min(dtmin,celldt);
297 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
299 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
303 using RF =
typename LFSU::Traits::FiniteElementType::
304 Traits::LocalBasisType::Traits::RangeFieldType;
307 const auto dim = IG::dimension;
310 auto cell_inside = ig.inside();
313 auto geo = ig.geometry();
314 auto geo_inside = cell_inside.geometry();
317 auto geo_in_inside = ig.geometryInInside();
321 auto face_local = ref_el.position(0,0);
324 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
328 auto local_inside = ref_el_inside.position(0,0);
331 auto bctype = param.bctype(ig.intersection(),face_local);
337 auto global_inside = geo_inside.global(local_inside);
338 auto global_outside = geo.global(face_local);
339 global_inside -= global_outside;
340 auto distance = global_inside.two_norm();
343 auto tensor_inside = param.A(cell_inside,local_inside);
344 auto n_F = ig.centerUnitOuterNormal();
345 Dune::FieldVector<RF,dim> An_F;
346 tensor_inside.mv(n_F,An_F);
347 auto k_inside = n_F*An_F;
350 auto iplocal_s = geo_in_inside.global(face_local);
351 auto g = param.g(cell_inside,iplocal_s);
354 auto b = param.b(cell_inside,iplocal_s);
355 auto n = ig.centerUnitOuterNormal();
358 r_s.accumulate(lfsu_s,0,(b*n)*g*face_volume + k_inside*(x_s(lfsu_s,0)-g)*face_volume/distance);
369 auto j = param.j(ig.intersection(),face_local);
372 r_s.accumulate(lfsu_s,0,j*face_volume);
380 auto iplocal_s = geo_in_inside.global(face_local);
381 auto b = param.b(cell_inside,iplocal_s);
382 auto n = ig.centerUnitOuterNormal();
385 auto o = param.o(ig.intersection(),face_local);
388 r_s.accumulate(lfsu_s,0,((b*n)*x_s(lfsu_s,0) + o)*face_volume);
394 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
396 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
400 using RF =
typename LFSU::Traits::FiniteElementType::
401 Traits::LocalBasisType::Traits::RangeFieldType;
404 const auto dim = IG::dimension;
407 auto cell_inside = ig.inside();
410 auto geo = ig.geometry();
411 auto geo_inside = cell_inside.geometry();
414 auto geo_in_inside = ig.geometryInInside();
418 auto face_local = ref_el.position(0,0);
421 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
425 auto local_inside = ref_el_inside.position(0,0);
428 auto bctype = param.bctype(ig.intersection(),face_local);
434 auto global_inside = geo_inside.global(local_inside);
435 auto global_outside = geo.global(face_local);
436 global_inside -= global_outside;
437 auto distance = global_inside.two_norm();
440 auto tensor_inside = param.A(cell_inside,local_inside);
441 auto n_F = ig.centerUnitOuterNormal();
442 Dune::FieldVector<RF,dim> An_F;
443 tensor_inside.mv(n_F,An_F);
444 auto k_inside = n_F*An_F;
447 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, k_inside*face_volume/distance );
455 auto iplocal_s = geo_in_inside.global(face_local);
456 auto b = param.b(cell_inside,iplocal_s);
457 auto n = ig.centerUnitOuterNormal();
460 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, (b*n)*face_volume );
467 template<
typename EG,
typename LFSV,
typename R>
471 const auto& cell = eg.entity();
474 auto geo = eg.geometry();
476 auto local_inside = ref_el.position(0,0);
479 auto f = param.f(cell,local_inside);
481 r.accumulate(lfsv,0,-f*geo.volume());
485 void setTime (
typename TP::Traits::RangeFieldType t)
491 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
497 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
504 else first_stage =
false;
513 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const 521 mutable typename TP::Traits::RangeFieldType dtmin;
522 mutable typename TP::Traits::RangeFieldType cellinflux;
553 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
554 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 557 const auto& cell = eg.entity();
560 auto geo = eg.geometry();
562 auto local_inside = ref_el.position(0,0);
565 auto capacity = param.d(cell,local_inside);
568 r.accumulate(lfsu,0,capacity*x(lfsu,0)*geo.volume());
572 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
579 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
584 const auto& cell = eg.entity();
587 auto geo = eg.geometry();
589 auto local_inside = ref_el.position(0,0);
592 auto capacity = param.d(cell,local_inside);
595 mat.accumulate(lfsu,0,lfsu,0,capacity*geo.volume());
607 #endif // DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH Definition: convectiondiffusionccfv.hh:56
const IG & ig
Definition: constraints.hh:148
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:468
static const int dim
Definition: adaptivity.hh:83
const Entity & e
Definition: localfunctionspace.hh:111
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:114
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:491
Definition: convectiondiffusionccfv.hh:535
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: convectiondiffusionccfv.hh:485
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
Definition: convectiondiffusionccfv.hh:50
ConvectionDiffusionCCFV(TP ¶m_)
Definition: convectiondiffusionccfv.hh:60
Definition: convectiondiffusionparameter.hh:65
Definition: convectiondiffusionccfv.hh:49
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: convectiondiffusionccfv.hh:86
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: convectiondiffusionccfv.hh:513
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Definition: convectiondiffusionccfv.hh:54
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:580
Definition: convectiondiffusionparameter.hh:65
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:298
Definition: convectiondiffusionccfv.hh:58
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:554
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:67
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: convectiondiffusionccfv.hh:497
sparsity pattern generator
Definition: pattern.hh:29
Definition: convectiondiffusionccfv.hh:55
ConvectionDiffusionCCFVTemporalOperator(TP ¶m_)
Definition: convectiondiffusionccfv.hh:547
void jacobian_apply_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, Y &y_s, Y &y_n) const
Definition: convectiondiffusionccfv.hh:184
void postStage()
to be called once at the end of each stage
Definition: convectiondiffusionccfv.hh:508
Definition: convectiondiffusionparameter.hh:65
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: convectiondiffusionccfv.hh:573
sparsity pattern generator
Definition: pattern.hh:13
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:275
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:395
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Definition: convectiondiffusionccfv.hh:57
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:93
Definition: convectiondiffusionccfv.hh:38
Definition: convectiondiffusionccfv.hh:53
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:193