2#ifndef DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH
3#define DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH
8#include<dune/geometry/referenceelements.hh>
10#include<dune/pdelab/common/quadraturerules.hh>
11#include<dune/pdelab/common/geometrywrapper.hh>
12#include<dune/pdelab/localoperator/pattern.hh>
13#include<dune/pdelab/localoperator/flags.hh>
14#include<dune/pdelab/localoperator/idefault.hh>
15#include<dune/pdelab/localoperator/defaultimp.hh>
16#include<dune/pdelab/localoperator/convectiondiffusionparameter.hh>
45 using BCType =
typename ConvectionDiffusionBoundaryConditions::Type;
49 enum { doPatternVolume =
true };
50 enum { doPatternSkeleton =
true };
53 enum { doAlphaVolume =
true };
54 enum { doAlphaSkeleton =
true };
55 enum { doAlphaBoundary =
true };
56 enum { doLambdaVolume =
true };
57 enum { doLambdaSkeleton =
false };
58 enum { doLambdaBoundary =
false };
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>
86 void jacobian_apply_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, Y& y)
const
88 alpha_volume(eg,lfsu,x,lfsv,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>
114 void alpha_skeleton (
const IG& ig,
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::Entity::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();
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>
184 void jacobian_apply_skeleton (
const IG& ig,
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
189 alpha_skeleton(ig,lfsu_s,x_s,lfsv_s,lfsu_n,x_n,lfsv_n,y_s,y_n);
192 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
193 void jacobian_skeleton (
const IG& ig,
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::Entity::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();
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>
275 void alpha_volume_post_skeleton(
const EG& eg,
const LFSU& lfsu,
const X& x,
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);
297 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
298 void alpha_boundary (
const IG& ig,
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::Entity::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);
333 if (bctype==ConvectionDiffusionBoundaryConditions::Dirichlet)
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();
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);
363 if (bctype==ConvectionDiffusionBoundaryConditions::Neumann)
369 auto j = param.j(ig.intersection(),face_local);
372 r_s.accumulate(lfsu_s,0,j*face_volume);
377 if (bctype==ConvectionDiffusionBoundaryConditions::Outflow)
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>
395 void jacobian_boundary (
const IG& ig,
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::Entity::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);
430 if (bctype==ConvectionDiffusionBoundaryConditions::Dirichlet)
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();
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 );
452 if (bctype==ConvectionDiffusionBoundaryConditions::Outflow)
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>
468 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
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;
542 enum { doPatternVolume =
true };
545 enum { doAlphaVolume =
true };
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>
573 void jacobian_apply_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, Y& y)
const
575 alpha_volume(eg,lfsu,x,lfsv,y);
579 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
580 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
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());
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Definition: convectiondiffusionccfv.hh:539
Definition: convectiondiffusionccfv.hh:44
void postStage()
to be called once at the end of each stage
Definition: convectiondiffusionccfv.hh:508
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: convectiondiffusionccfv.hh:485
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: convectiondiffusionccfv.hh:513
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
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:30
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Default flags for all local operators.
Definition: flags.hh:19
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
A few common exception classes.
Traits for type conversions and type information.
Implements a vector constructed from a given type representing a field and a compile-time given size.
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Dune namespace.
Definition: alignedallocator.hh:13