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>
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 const auto& cell = eg.entity();
72 auto geo = eg.geometry();
73 auto ref_el = referenceElement(geo);
74 auto local_inside = ref_el.position(0,0);
77 auto c = param.c(cell,local_inside);
80 r.accumulate(lfsu,0,(c*x(lfsu,0))*geo.volume());
84 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
91 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
92 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
96 const auto& cell = eg.entity();
99 auto geo = eg.geometry();
100 auto ref_el = referenceElement(geo);
101 auto local_inside = ref_el.position(0,0);
104 auto c = param.c(cell,local_inside);
107 mat.accumulate(lfsu,0,lfsu,0,c*geo.volume());
112 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
114 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
115 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
116 R& r_s, R& r_n)
const
119 using RF =
typename LFSU::Traits::FiniteElementType::
120 Traits::LocalBasisType::Traits::RangeFieldType;
123 const auto dim = IG::Entity::dimension;
126 auto cell_inside =
ig.inside();
127 auto cell_outside =
ig.outside();
130 auto geo =
ig.geometry();
131 auto geo_inside = cell_inside.geometry();
132 auto geo_outside = cell_outside.geometry();
135 auto geo_in_inside =
ig.geometryInInside();
138 auto ref_el = referenceElement(geo);
139 auto face_local = ref_el.position(0,0);
142 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
145 auto ref_el_inside = referenceElement(geo_inside);
146 auto ref_el_outside = referenceElement(geo_outside);
147 auto local_inside = ref_el_inside.position(0,0);
148 auto local_outside = ref_el_outside.position(0,0);
151 auto tensor_inside = param.A(cell_inside,local_inside);
152 auto tensor_outside = param.A(cell_outside,local_outside);
153 auto n_F =
ig.centerUnitOuterNormal();
154 Dune::FieldVector<RF,dim> An_F;
155 tensor_inside.mv(n_F,An_F);
156 auto k_inside = n_F*An_F;
157 tensor_outside.mv(n_F,An_F);
158 auto k_outside = n_F*An_F;
159 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
162 auto iplocal_s = geo_in_inside.global(face_local);
163 auto b = param.b(cell_inside,iplocal_s);
166 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
169 auto global_inside = geo_inside.global(local_inside);
170 auto global_outside = geo_outside.global(local_outside);
173 global_inside -= global_outside;
174 auto distance = global_inside.two_norm();
177 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);
178 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);
182 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
184 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
185 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
186 Y& y_s, Y& y_n)
const
191 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
193 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
194 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
195 M& mat_ss, M& mat_sn,
196 M& mat_ns, M& mat_nn)
const
199 using RF =
typename LFSU::Traits::FiniteElementType::
200 Traits::LocalBasisType::Traits::RangeFieldType;
203 const auto dim = IG::Entity::dimension;
206 auto cell_inside =
ig.inside();
207 auto cell_outside =
ig.outside();
210 auto geo =
ig.geometry();
211 auto geo_inside = cell_inside.geometry();
212 auto geo_outside = cell_outside.geometry();
215 auto geo_in_inside =
ig.geometryInInside();
218 auto ref_el = referenceElement(geo);
219 auto face_local = ref_el.position(0,0);
222 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
225 auto ref_el_inside = referenceElement(geo_inside);
226 auto ref_el_outside = referenceElement(geo_outside);
227 auto local_inside = ref_el_inside.position(0,0);
228 auto local_outside = ref_el_outside.position(0,0);
231 auto tensor_inside = param.A(cell_inside,local_inside);
232 auto tensor_outside = param.A(cell_outside,local_outside);
233 auto n_F =
ig.centerUnitOuterNormal();
234 Dune::FieldVector<RF,dim> An_F;
235 tensor_inside.mv(n_F,An_F);
236 auto k_inside = n_F*An_F;
237 tensor_outside.mv(n_F,An_F);
238 auto k_outside = n_F*An_F;
239 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
242 auto iplocal_s = geo_in_inside.global(face_local);
243 auto b = param.b(cell_inside,iplocal_s);
247 auto global_inside = geo_inside.global(local_inside);
248 auto global_outside = geo_outside.global(local_outside);
251 global_inside -= global_outside;
252 auto distance = global_inside.two_norm();
255 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, k_avg*face_volume/distance );
256 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -k_avg*face_volume/distance );
257 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, -k_avg*face_volume/distance );
258 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, k_avg*face_volume/distance );
261 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, vn*face_volume );
262 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -vn*face_volume );
266 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, vn*face_volume );
267 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, -vn*face_volume );
273 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
275 const LFSV& lfsv, R& r)
const
277 if (not first_stage)
return;
280 const auto& cell = eg.entity();
283 auto geo = eg.geometry();
284 auto ref_el = referenceElement(geo);
285 auto local_inside = ref_el.position(0,0);
288 auto cellcapacity = param.d(cell,local_inside)*geo.volume();
289 auto celldt = cellcapacity/(cellinflux+1E-30);
290 dtmin = std::min(dtmin,celldt);
296 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
298 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
302 using RF =
typename LFSU::Traits::FiniteElementType::
303 Traits::LocalBasisType::Traits::RangeFieldType;
306 const auto dim = IG::Entity::dimension;
309 auto cell_inside =
ig.inside();
312 auto geo =
ig.geometry();
313 auto geo_inside = cell_inside.geometry();
316 auto geo_in_inside =
ig.geometryInInside();
319 auto ref_el = referenceElement(geo);
320 auto face_local = ref_el.position(0,0);
323 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
326 auto ref_el_inside = referenceElement(geo_inside);
327 auto local_inside = ref_el_inside.position(0,0);
330 auto bctype = param.bctype(
ig.intersection(),face_local);
336 auto global_inside = geo_inside.global(local_inside);
337 auto global_outside = geo.global(face_local);
338 global_inside -= global_outside;
339 auto distance = global_inside.two_norm();
342 auto tensor_inside = param.A(cell_inside,local_inside);
343 auto n_F =
ig.centerUnitOuterNormal();
344 Dune::FieldVector<RF,dim> An_F;
345 tensor_inside.mv(n_F,An_F);
346 auto k_inside = n_F*An_F;
349 auto iplocal_s = geo_in_inside.global(face_local);
350 auto g = param.g(cell_inside,iplocal_s);
353 auto b = param.b(cell_inside,iplocal_s);
354 auto n =
ig.centerUnitOuterNormal();
357 r_s.accumulate(lfsu_s,0,(b*n)*g*face_volume + k_inside*(x_s(lfsu_s,0)-g)*face_volume/distance);
368 auto j = param.j(
ig.intersection(),face_local);
371 r_s.accumulate(lfsu_s,0,j*face_volume);
379 auto iplocal_s = geo_in_inside.global(face_local);
380 auto b = param.b(cell_inside,iplocal_s);
381 auto n =
ig.centerUnitOuterNormal();
384 auto o = param.o(
ig.intersection(),face_local);
387 r_s.accumulate(lfsu_s,0,((b*n)*x_s(lfsu_s,0) + o)*face_volume);
393 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
395 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
399 using RF =
typename LFSU::Traits::FiniteElementType::
400 Traits::LocalBasisType::Traits::RangeFieldType;
403 const auto dim = IG::Entity::dimension;
406 auto cell_inside =
ig.inside();
409 auto geo =
ig.geometry();
410 auto geo_inside = cell_inside.geometry();
413 auto geo_in_inside =
ig.geometryInInside();
416 auto ref_el = referenceElement(geo);
417 auto face_local = ref_el.position(0,0);
420 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
423 auto ref_el_inside = referenceElement(geo_inside);
424 auto local_inside = ref_el_inside.position(0,0);
427 auto bctype = param.bctype(
ig.intersection(),face_local);
433 auto global_inside = geo_inside.global(local_inside);
434 auto global_outside = geo.global(face_local);
435 global_inside -= global_outside;
436 auto distance = global_inside.two_norm();
439 auto tensor_inside = param.A(cell_inside,local_inside);
440 auto n_F =
ig.centerUnitOuterNormal();
441 Dune::FieldVector<RF,dim> An_F;
442 tensor_inside.mv(n_F,An_F);
443 auto k_inside = n_F*An_F;
446 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, k_inside*face_volume/distance );
454 auto iplocal_s = geo_in_inside.global(face_local);
455 auto b = param.b(cell_inside,iplocal_s);
456 auto n =
ig.centerUnitOuterNormal();
459 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, (b*n)*face_volume );
466 template<
typename EG,
typename LFSV,
typename R>
470 const auto& cell = eg.entity();
473 auto geo = eg.geometry();
474 auto ref_el = referenceElement(geo);
475 auto local_inside = ref_el.position(0,0);
478 auto f = param.f(cell,local_inside);
480 r.accumulate(lfsv,0,-f*geo.volume());
484 void setTime (
typename TP::Traits::RangeFieldType t)
490 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
496 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
503 else first_stage =
false;
512 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const
520 mutable typename TP::Traits::RangeFieldType dtmin;
521 mutable typename TP::Traits::RangeFieldType cellinflux;
552 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
553 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
556 const auto& cell = eg.entity();
559 auto geo = eg.geometry();
560 auto ref_el = referenceElement(geo);
561 auto local_inside = ref_el.position(0,0);
564 auto capacity = param.d(cell,local_inside);
567 r.accumulate(lfsu,0,capacity*x(lfsu,0)*geo.volume());
571 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
578 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
583 const auto& cell = eg.entity();
586 auto geo = eg.geometry();
587 auto ref_el = referenceElement(geo);
588 auto local_inside = ref_el.position(0,0);
591 auto capacity = param.d(cell,local_inside);
594 mat.accumulate(lfsu,0,lfsu,0,capacity*geo.volume());
606 #endif // DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH