dune-pdelab  2.7-git
twophaseccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
4 
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/geometry/referenceelements.hh>
8 #include<dune/localfunctions/raviartthomas/raviartthomascube.hh>
9 
12 
13 #include"../common/function.hh"
14 #include"pattern.hh"
15 #include"flags.hh"
16 #include"idefault.hh"
17 
18 namespace Dune {
19  namespace PDELab {
20 
22  template<typename GV, typename RF>
24  {
26  using GridViewType = GV;
27 
29  enum {
31  dimDomain = GV::dimension
32  };
33 
35  using DomainFieldType = typename GV::Grid::ctype;
36 
38  using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
39 
41  using IntersectionDomainType = Dune::FieldVector<DomainFieldType,dimDomain-1>;
42 
44  using RangeFieldType = RF;
45 
47  using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
48 
51 
53  using ElementType = typename GV::Traits::template Codim<0>::Entity;
54  using IntersectionType = typename GV::Intersection;
55  };
56 
57  template<typename GV, typename RF>
59  {
62 
64  using PermTensorType = Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>;
65  };
66 
68  template<class T, class Imp>
70  {
71  public:
72  using Traits = T;
73 
75  typename Traits::RangeFieldType
76  phi (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
77  {
78  return asImp().phi(e,x);
79  }
80 
82  typename Traits::RangeFieldType
83  pc (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
84  typename Traits::RangeFieldType s_l) const
85  {
86  return asImp().pc(e,x,s_l);
87  }
88 
90  typename Traits::RangeFieldType
91  s_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
92  typename Traits::RangeFieldType pc) const
93  {
94  return asImp().s_l(e,x,pc);
95  }
96 
98  typename Traits::RangeFieldType
99  kr_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
100  typename Traits::RangeFieldType s_l) const
101  {
102  return asImp().kr_l(e,x,s_l);
103  }
104 
106  typename Traits::RangeFieldType
107  kr_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
108  typename Traits::RangeFieldType s_g) const
109  {
110  return asImp().kr_g(e,x,s_g);
111  }
112 
114  typename Traits::RangeFieldType
115  mu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
116  typename Traits::RangeFieldType p_l) const
117  {
118  return asImp().mu_l(e,x,p_l);
119  }
120 
122  typename Traits::RangeFieldType
123  mu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
124  typename Traits::RangeFieldType p_g) const
125  {
126  return asImp().mu_l(e,x,p_g);
127  }
128 
130  typename Traits::PermTensorType
131  k_abs (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
132  {
133  return asImp().k_abs(e,x);
134  }
135 
137  const typename Traits::RangeType& gravity () const
138  {
139  return asImp().gravity();
140  }
141 
143  template<typename E>
144  typename Traits::RangeFieldType
145  nu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
146  typename Traits::RangeFieldType p_l) const
147  {
148  return asImp().nu_l(e,x,p_l);
149  }
150 
152  typename Traits::RangeFieldType
153  nu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
154  typename Traits::RangeFieldType p_g) const
155  {
156  return asImp().nu_g(e,x,p_g);
157  }
158 
160  typename Traits::RangeFieldType
161  rho_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
162  typename Traits::RangeFieldType p_l) const
163  {
164  return asImp().rho_l(e,x,p_l);
165  }
166 
168  typename Traits::RangeFieldType
169  rho_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
170  typename Traits::RangeFieldType p_g) const
171  {
172  return asImp().rho_g(e,x,p_g);
173  }
174 
176  int
177  bc_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
178  {
179  return asImp().bc_l(is,x,time);
180  }
181 
183  int
184  bc_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
185  {
186  return asImp().bc_g(is,x,time);
187  }
188 
190  typename Traits::RangeFieldType
191  g_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
192  {
193  return asImp().g_l(is,x,time);
194  }
195 
197  typename Traits::RangeFieldType
198  g_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
199  {
200  return asImp().g_g(is,x,time);
201  }
202 
204  typename Traits::RangeFieldType
205  j_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
206  {
207  return asImp().j_l(is,x,time);
208  }
209 
211  typename Traits::RangeFieldType
212  j_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
213  {
214  return asImp().j_g(is,x,time);
215  }
216 
218  typename Traits::RangeFieldType
219  q_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
220  typename Traits::RangeFieldType time) const
221  {
222  return asImp().q_l(e,x,time);
223  }
224 
226  typename Traits::RangeFieldType
227  q_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
228  typename Traits::RangeFieldType time) const
229  {
230  return asImp().q_g(e,x,time);
231  }
232 
233  private:
234  Imp& asImp () {return static_cast<Imp &> (*this);}
235  const Imp& asImp () const {return static_cast<const Imp &>(*this);}
236  };
237 
238 
239  // a local operator for solving the two-phase flow in pressure-pressure formulation
240  // with two-point flux approximation
241  // TP : parameter class, see above
242  // V : Vector holding last time step
243  template<typename TP>
245  : public NumericalJacobianSkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
246  public NumericalJacobianApplySkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
247 
248  public NumericalJacobianBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
249  public NumericalJacobianApplyBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
250 
251  public FullSkeletonPattern,
252  public FullVolumePattern,
254 
255  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
256  {
257  enum { dim = TP::Traits::GridViewType::dimension };
258  enum { liquid = 0 };
259  enum { gas = 1 };
260 
261  using Real = typename TP::Traits::RangeFieldType;
262  public:
263  // pattern assembly flags
264  enum { doPatternVolume = true };
265  enum { doPatternSkeleton = true };
266 
267  // residual assembly flags
268  enum { doAlphaSkeleton = true };
269  enum { doAlphaBoundary = true };
270  enum { doLambdaVolume = true };
271  enum { doLambdaBoundary = true };
272 
274  TwoPhaseTwoPointFluxOperator (const TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
275  : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
276  {}
277 
278  // volume integral depending only on test functions
279  template<typename EG, typename LFSV, typename R>
280  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
281  {
282  // Reference to cell
283  const auto& cell = eg.entity();
284 
285  // get geometry
286  auto geo = eg.geometry();
287 
288  // cell geometry
289  auto ref_el = referenceElement(geo);
290  auto cell_center_local = ref_el.position(0,0);
291  auto cell_volume = geo.volume();
292 
293  // contribution from source term
294  r.accumulate(lfsv, liquid, -scale_l * tp.q_l(cell,cell_center_local,time) * cell_volume);
295  r.accumulate(lfsv, gas, -scale_g * tp.q_g(cell,cell_center_local,time) * cell_volume);
296  }
297 
298  // skeleton integral depending on test and ansatz functions
299  // each face is only visited ONCE!
300  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
301  void alpha_skeleton (const IG& ig,
302  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
303  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
304  R& r_s, R& r_n) const
305  {
306  // Define types
307  using namespace Indices;
308  using PLSpace = TypeTree::Child<LFSV,_0>;
309  using RF = typename PLSpace::Traits::FiniteElementType::
310  Traits::LocalBasisType::Traits::RangeFieldType;
311 
312  // References to inside and outside cells
313  const auto& cell_inside = ig.inside();
314  const auto& cell_outside = ig.outside();
315 
316  // get geometries
317  auto geo = ig.geometry();
318  auto geo_inside = cell_inside.geometry();
319  auto geo_outside = cell_outside.geometry();
320 
321  // cell geometries
322  auto ref_el_inside = referenceElement(geo_inside);
323  auto ref_el_outside = referenceElement(geo_outside);
324  auto inside_cell_center_local = ref_el_inside.position(0,0);
325  auto outside_cell_center_local = ref_el_outside.position(0,0);
326  auto inside_cell_center_global = geo_inside.center();
327  auto outside_cell_center_global = geo_outside.center();
328 
329  // distance of cell centers
330  auto d = outside_cell_center_global;
331  d -= inside_cell_center_global;
332  auto distance = d.two_norm();
333 
334  // face geometry
335  auto ref_el = referenceElement(geo);
336  auto face_local = ref_el.position(0,0);
337  auto face_volume = geo.volume();
338 
339  // absolute permeability
340  auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
341  auto k_abs_outside = tp.k_abs(cell_outside,outside_cell_center_local);
342 
343  // gravity times normal
344  auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
345 
346  // liquid phase calculation
347  auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
348  auto rho_l_outside = tp.rho_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
349  auto w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction
350  RF pc_upwind, s_l_upwind, s_g_upwind;
351  auto nu_l = aavg(tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)),
352  tp.nu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid)));
353  if (w_l>=0) // upwind capillary pressure on face
354  {
355  pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
356  s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
357  }
358  else
359  {
360  pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
361  s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
362  }
363  s_g_upwind = 1-s_l_upwind;
364  auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l_upwind)/
365  tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
366  auto lambda_l_outside = tp.kr_l(cell_outside,outside_cell_center_local,s_l_upwind)/
367  tp.mu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
368  auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
369 
370  r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
371  r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume);
372 
373  // gas phase calculation
374  auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
375  auto rho_g_outside = tp.rho_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
376  auto w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction
377  auto nu_g = aavg(tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)),
378  tp.nu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas)));
379  if (w_l*w_g<=0) // new evaluation necessary only if signs differ
380  {
381  if (w_g>=0) // upwind capillary pressure on face
382  {
383  pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
384  s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
385  }
386  else
387  {
388  pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
389  s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
390  }
391  s_g_upwind = 1-s_l_upwind;
392  }
393  auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g_upwind)/
394  tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
395  auto lambda_g_outside = tp.kr_g(cell_outside,outside_cell_center_local,s_g_upwind)/
396  tp.mu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
397  auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
398 
399  r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume);
400  r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume);
401  }
402 
403  // skeleton integral depending on test and ansatz functions
404  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
405  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
406  void alpha_boundary (const IG& ig,
407  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
408  R& r_s) const
409  {
410  // References to inside cell
411  const auto& cell_inside = ig.inside();
412 
413  // get geometries
414  auto geo = ig.geometry();
415  auto geo_inside = cell_inside.geometry();
416 
417  // face geometry
418  auto ref_el = referenceElement(geo);
419  auto face_local = ref_el.position(0,0);
420  auto face_volume = geo.volume();
421 
422  // evaluate boundary condition type
423  auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
424  auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
425  if (bc_l!=1 && bc_g!=1) return; // no Dirichlet boundary conditions
426 
427  // cell geometry
428  auto ref_el_inside = referenceElement(geo_inside);
429  auto inside_cell_center_local = ref_el_inside.position(0,0);
430  auto inside_cell_center_global = geo_inside.center();
431 
432  // distance of cell center to boundary
433  auto d = geo.global(face_local);
434  d -= inside_cell_center_global;
435  auto distance = d.two_norm();
436 
437  // absolute permeability
438  auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
439 
440  // gravity times normal
441  auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
442 
443  // liquid phase Dirichlet boundary
444  if (bc_l==1)
445  {
446  auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
447  auto g_l = tp.g_l(ig.intersection(),face_local,time);
448  auto w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn;
449  auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
450  auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l)/
451  tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
452  auto sigma_l = lambda_l_inside*k_abs_inside;
453  auto nu_l = tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
454  r_s.accumulate(lfsv_s, liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
455  }
456 
457  // gas phase Dirichlet boundary
458  if (bc_g==1)
459  {
460  auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
461  auto g_g = tp.g_g(ig.intersection(),face_local,time);
462  auto w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn;
463  auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
464  auto s_g = 1-s_l;
465  auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g)/
466  tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
467  auto sigma_g = lambda_g_inside*k_abs_inside;
468  auto nu_g = tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
469  r_s.accumulate(lfsv_s, gas, scale_l * nu_g * sigma_g * w_g * face_volume);
470  }
471  }
472 
473  // boundary integral independent of ansatz functions
474  template<typename IG, typename LFSV, typename R>
475  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r_s) const
476  {
477  // get geometries
478  auto geo = ig.geometry();
479 
480  // face geometry
481  auto ref_el = referenceElement(geo);
482  auto face_local = ref_el.position(0,0);
483  auto face_volume = geo.volume();
484 
485  // evaluate boundary condition type
486  auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
487  auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
488  if (bc_l!=0 && bc_g!=0) return; // no Neumann boundary conditions
489 
490  // liquid phase Neumann boundary
491  if (bc_l==0)
492  {
493  auto j_l = tp.j_l(ig.intersection(),face_local,time);
494  r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
495  }
496 
497  // gas phase Neumann boundary
498  if (bc_g==0)
499  {
500  auto j_g = tp.j_g(ig.intersection(),face_local,time);
501  r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
502  }
503  }
504 
506  void setTime (typename TP::Traits::RangeFieldType t)
507  {
508  time = t;
509  }
510 
511  private:
512  const TP& tp; // two phase parameter class
513  typename TP::Traits::RangeFieldType time;
514  Real scale_l, scale_g;
515 
516  template<typename T>
517  T aavg (T a, T b) const
518  {
519  return 0.5*(a+b);
520  }
521 
522  template<typename T>
523  T havg (T a, T b) const
524  {
525  T eps = 1E-30;
526  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
527  }
528  };
529 
530 
537  template<class TP>
539  : public NumericalJacobianVolume<TwoPhaseOnePointTemporalOperator<TP> >,
540  public NumericalJacobianApplyVolume<TwoPhaseOnePointTemporalOperator<TP> >,
541  public FullVolumePattern,
543  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
544  {
545  enum { dim = TP::Traits::GridViewType::dimension };
546  enum { liquid = 0 };
547  enum { gas = 1 };
548 
549  using Real = typename TP::Traits::RangeFieldType;
550 
551  public:
552  // pattern assembly flags
553  enum { doPatternVolume = true };
554 
555  // residual assembly flags
556  enum { doAlphaVolume = true };
557 
558  TwoPhaseOnePointTemporalOperator (TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
559  : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
560  {
561  }
562 
564  void setTime (typename TP::Traits::RangeFieldType t)
565  {
566  time = t;
567  }
568 
569  // volume integral depending on test and ansatz functions
570  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
571  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
572  {
573  // Reference to cell
574  const auto& cell = eg.entity();
575 
576  // get geometry
577  auto geo = eg.geometry();
578 
579  // cell geometry
580  auto ref_el = referenceElement(geo);
581  auto cell_center_local = ref_el.position(0,0);
582  auto cell_volume = geo.volume();
583 
584  auto phi = tp.phi(cell,cell_center_local);
585  auto s_l = tp.s_l(cell,cell_center_local,x(lfsu,gas)-x(lfsu,liquid));
586 
587  r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(cell,cell_center_local,x(lfsu,liquid)) * cell_volume);
588  r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(cell,cell_center_local,x(lfsu,gas)) * cell_volume);
589  }
590 
591  private:
592  TP& tp;
593  typename TP::Traits::RangeFieldType time;
594  Real scale_l, scale_g;
595  };
596 
597 
606  template<typename TP, typename PL, typename PG>
607  class V_l
608  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
609  typename PL::Traits::RangeFieldType,
610  PL::Traits::GridViewType::dimension,
611  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
612  V_l<TP,PL,PG> >
613  {
614  // extract useful types
615  using GV = typename PL::Traits::GridViewType;
616  using DF = typename GV::Grid::ctype;
617  using RF = typename PL::Traits::RangeFieldType;
618  using RangeType = typename PL::Traits::RangeType;
619  enum { dim = PL::Traits::GridViewType::dimension };
620  using Element = typename GV::Traits::template Codim<0>::Entity;
621  using IntersectionIterator = typename GV::IntersectionIterator;
622  using Intersection = typename IntersectionIterator::Intersection;
623 
624  const TP& tp;
625  const PL& pl;
626  const PG& pg;
627  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
628  typename TP::Traits::RangeFieldType time;
629 
630 
631  using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
632 
633  public:
635 
637 
638  V_l (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
639 
640  // set time where operator is to be evaluated (i.e. end of the time intervall)
641  void set_time (typename TP::Traits::RangeFieldType time_)
642  {
643  time = time_;
644  }
645 
646  inline void evaluate (const typename Traits::ElementType& e,
647  const typename Traits::DomainType& x,
648  typename Traits::RangeType& y) const
649  {
650  // get geometries
651  auto geo = e.geometry();
652 
653  // face geometry
654  auto ref_el = referenceElement(geo);
655  auto face_local = ref_el.position(0,0);
656  auto face_volume = geo.volume();
657 
658  // cell geometry
659  auto inside_cell_center_local = ref_el.position(0,0);
660  auto inside_cell_center_global = geo.global(inside_cell_center_local);
661 
662  // absolute permeability
663  auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
664 
665  // pressure evaluation
666  typename PL::Traits::RangeType pl_inside, pg_inside;
667  pl.evaluate(e,inside_cell_center_local,pl_inside);
668  pg.evaluate(e,inside_cell_center_local,pg_inside);
669 
670  // density
671  auto nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
672 
673  // for coefficient computation
674  RF vn[2*dim]; // normal velocities
675  RF coeff[2*dim]; // RT0 coefficient
676  auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear
677  auto determinant = B.determinant();
678 
679  // loop over cell neighbors
680  for (const auto& intersection : intersections(pl.getGridView(),e))
681  {
682  // set to zero for processor boundary
683  vn[intersection.indexInInside()] = 0.0;
684 
685  // get geometry
686  auto geo_intersection = intersection.geometry();
687 
688  // face geometry
689  const auto& face_local = referenceElement(geo_intersection).position(0,0);
690 
691  // interior face
692  if (intersection.neighbor())
693  {
694  // get geometry
695  auto geo_outside = intersection.outside().geometry();
696  auto ref_el_outside = referenceElement(geo_outside);
697  auto outside_cell_center_local = ref_el_outside.position(0,0);
698  auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
699 
700  // distance of cell centers
701  auto d = outside_cell_center_global;
702  d -= inside_cell_center_global;
703  auto distance = d.two_norm();
704 
705  // absolute permeability
706  auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
707 
708  // gravity times normal
709  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
710 
711  // pressure evaluation
712  typename PL::Traits::RangeType pl_outside, pg_outside;
713  pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
714  pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
715 
716  // density
717  auto nu_l_outside = tp.nu_l(*(intersection.outside()),outside_cell_center_local,pg_outside);
718 
719  // liquid phase calculation
720  auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
721  auto rho_l_outside = tp.rho_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
722  auto w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction
723  RF pc_upwind, s_l_upwind, s_g_upwind;
724  if (w_l>=0) // upwind capillary pressure on face
725  {
726  pc_upwind = pg_inside-pl_inside;
727  s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
728  }
729  else
730  {
731  pc_upwind = pg_outside-pl_outside;
732  s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
733  }
734  s_g_upwind = 1-s_l_upwind;
735  auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/
736  tp.mu_l(e,inside_cell_center_local,pl_inside);
737  auto lambda_l_outside = tp.kr_l(*(intersection.outside()),outside_cell_center_local,s_l_upwind)/
738  tp.mu_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
739  auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
740  auto nu_l = aavg(nu_l_inside,nu_l_outside);
741 
742  // set coefficient
743  vn[intersection.indexInInside()] = nu_l * sigma_l * w_l;
744  }
745 
746  // boundary face
747  if (intersection.boundary())
748  {
749  // distance of cell center to boundary
750  auto d = geo_intersection.global(face_local);
751  d -= inside_cell_center_global;
752  auto distance = d.two_norm();
753 
754  // evaluate boundary condition type
755  auto bc_l = tp.bc_l(intersection,face_local,time);
756 
757  // liquid phase Dirichlet boundary
758  if (bc_l==1)
759  {
760  auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
761  auto g_l = tp.g_l(intersection,face_local,time);
762  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
763  auto w_l = (pl_inside-g_l)/distance + rho_l_inside*gn;
764  auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
765  auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/
766  tp.mu_l(e,inside_cell_center_local,pl_inside);
767  auto sigma_l = lambda_l_inside*k_abs_inside;
768  vn[intersection.indexInInside()] = nu_l_inside * sigma_l * w_l;
769  }
770 
771  // liquid phase Neumann boundary
772  if (bc_l==0)
773  {
774  auto j_l = tp.j_l(intersection,face_local,time);
775  vn[intersection.indexInInside()] = j_l;
776  }
777  }
778 
779  // compute coefficient
780  auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element
781  vstar *= vn[intersection.indexInInside()];
782  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
783  if (intersection.indexInInside()%2==0)
784  normalhat[intersection.indexInInside()/2] = -1.0;
785  else
786  normalhat[intersection.indexInInside()/2] = 1.0;
787  Dune::FieldVector<DF,dim> vstarhat(0);
788  B.umtv(vstar,vstarhat); // Piola backward transformation
789  vstarhat *= determinant;
790  coeff[intersection.indexInInside()] = vstarhat*normalhat;
791  }
792 
793  // compute velocity on reference element
794  std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
795  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
796  typename Traits::RangeType yhat(0);
797  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
798  yhat.axpy(coeff[i],rt0vectors[i]);
799 
800  // apply Piola transformation
801  B.invert();
802  y = 0;
803  B.umtv(yhat,y);
804  y /= determinant;
805 
806  // std::cout << "vn= " ;
807  // for (int i=0; i<2*dim; i++) std::cout << vn[i] << " ";
808  // std::cout << std::endl;
809  // std::cout << "V_l: x=" << x << " y=" << y << std::endl;
810  }
811 
812  inline const typename Traits::GridViewType& getGridView ()
813  {
814  return pl.getGridView();
815  }
816 
817  private:
818 
819  template<typename T>
820  T aavg (T a, T b) const
821  {
822  return 0.5*(a+b);
823  }
824 
825  template<typename T>
826  T havg (T a, T b) const
827  {
828  T eps = 1E-30;
829  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
830  }
831  };
832 
841  template<typename TP, typename PL, typename PG>
842  class V_g
843  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
844  typename PL::Traits::RangeFieldType,
845  PL::Traits::GridViewType::dimension,
846  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
847  V_g<TP,PL,PG> >
848  {
849  // extract useful types
850  using GV = typename PL::Traits::GridViewType;
851  using DF = typename GV::Grid::ctype;
852  using RF = typename PL::Traits::RangeFieldType;
853  using RangeType = typename PL::Traits::RangeType;
854  enum { dim = PL::Traits::GridViewType::dimension };
855  using Element = typename GV::Traits::template Codim<0>::Entity;
856  using IntersectionIterator = typename GV::IntersectionIterator;
857  using Intersection = typename IntersectionIterator::Intersection;
858 
859  const TP& tp;
860  const PL& pl;
861  const PG& pg;
862  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
863  typename TP::Traits::RangeFieldType time;
864 
865 
866  using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
867 
868  public:
870 
872 
873  V_g (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
874 
875  // set time where operator is to be evaluated (i.e. end of the time intervall)
876  void set_time (typename TP::Traits::RangeFieldType time_)
877  {
878  time = time_;
879  }
880 
881  inline void evaluate (const typename Traits::ElementType& e,
882  const typename Traits::DomainType& x,
883  typename Traits::RangeType& y) const
884  {
885  // get geometries
886  auto geo = e.geometry();
887 
888  // face geometry
889  auto ref_el = referenceElement(geo);
890  auto face_local = ref_el.position(0,0);
891  auto face_volume = geo.volume();
892 
893  // cell geometry
894  auto inside_cell_center_local = ref_el.position(0,0);
895  auto inside_cell_center_global = geo.global(inside_cell_center_local);
896 
897  // absolute permeability
898  auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
899 
900  // pressure evaluation
901  typename PL::Traits::RangeType pl_inside, pg_inside;
902  pl.evaluate(e,inside_cell_center_local,pl_inside);
903  pg.evaluate(e,inside_cell_center_local,pg_inside);
904 
905  // density evaluation
906  auto nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
907 
908  // for coefficient computation
909  RF vn[2*dim]; // normal velocities
910  RF coeff[2*dim]; // RT0 coefficient
911  auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear
912  auto determinant = B.determinant();
913 
914  // loop over cell neighbors
915  for (const auto& intersection : intersections(pl.getGridView(),e))
916  {
917  // set to zero for processor boundary
918  vn[intersection.indexInInside()] = 0.0;
919 
920  // get geometry
921  auto geo_intersection = intersection.geometry();
922 
923  // face geometry
924  const auto& face_local = referenceElement(geo_intersection).position(0,0);
925 
926  // interior face
927  if (intersection.neighbor())
928  {
929  // get geometry
930  auto geo_outside = intersection.outside().geometry();
931  auto ref_el_outside = referenceElement(geo_outside);
932  auto outside_cell_center_local = ref_el_outside.position(0,0);
933  auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
934 
935  // distance of cell centers
936  auto d = outside_cell_center_global;
937  d -= inside_cell_center_global;
938  auto distance = d.two_norm();
939 
940  // absolute permeability
941  auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
942 
943  // gravity times normal
944  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
945 
946  // pressure evaluation
947  typename PL::Traits::RangeType pl_outside, pg_outside;
948  pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
949  pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
950 
951  // needed for both phases
952  RF pc_upwind, s_l_upwind, s_g_upwind;
953 
954  // gas phase calculation
955  auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
956  auto rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside);
957  auto w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction
958  if (w_g>=0) // upwind capillary pressure on face
959  {
960  pc_upwind = pg_inside-pl_inside;
961  s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
962  }
963  else
964  {
965  pc_upwind = pg_outside-pl_outside;
966  s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
967  }
968  s_g_upwind = 1-s_l_upwind;
969  auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/
970  tp.mu_g(e,inside_cell_center_local,pg_inside);
971  auto lambda_g_outside = tp.kr_g(*(intersection.outside()),outside_cell_center_local,s_g_upwind)/
972  tp.mu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
973  auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
974 
975  auto nu_g_outside = tp.nu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
976 
977  // set coefficient
978  vn[intersection.indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
979  }
980 
981  // boundary face
982  if (intersection.boundary())
983  {
984  // distance of cell center to boundary
985  auto d = geo_intersection.global(face_local);
986  d -= inside_cell_center_global;
987  auto distance = d.two_norm();
988 
989  // evaluate boundary condition type
990  auto bc_g = tp.bc_g(intersection,face_local,time);
991 
992  // gas phase Dirichlet boundary
993  if (bc_g==1)
994  {
995  auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
996  auto g_g = tp.g_g(intersection,face_local,time);
997  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
998  auto w_g = (pg_inside-g_g)/distance + rho_g_inside*gn;
999  auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
1000  auto s_g = 1-s_l;
1001  auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/
1002  tp.mu_g(e,inside_cell_center_local,pg_inside);
1003  auto sigma_g = lambda_g_inside*k_abs_inside;
1004 
1005  vn[intersection.indexInInside()] = nu_g_inside * sigma_g * w_g;
1006  }
1007 
1008  // gas phase Neumann boundary
1009  if (bc_g==0)
1010  {
1011  auto j_g = tp.j_g(intersection,face_local,time);
1012  vn[intersection.indexInInside()] = j_g; /* /nu_g_inside*/;
1013  }
1014  }
1015 
1016  // compute coefficient
1017  auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element
1018  vstar *= vn[intersection.indexInInside()];
1019  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
1020  if (intersection.indexInInside()%2==0)
1021  normalhat[intersection.indexInInside()/2] = -1.0;
1022  else
1023  normalhat[intersection.indexInInside()/2] = 1.0;
1024  Dune::FieldVector<DF,dim> vstarhat(0);
1025  B.umtv(vstar,vstarhat); // Piola backward transformation
1026  vstarhat *= determinant;
1027  coeff[intersection.indexInInside()] = vstarhat*normalhat;
1028  }
1029 
1030  // compute velocity on reference element
1031  std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1032  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1033  typename Traits::RangeType yhat(0);
1034  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1035  yhat.axpy(coeff[i],rt0vectors[i]);
1036 
1037  // apply Piola transformation
1038  B.invert();
1039  y = 0;
1040  B.umtv(yhat,y);
1041  y /= determinant;
1042  }
1043 
1044  inline const typename Traits::GridViewType& getGridView ()
1045  {
1046  return pl.getGridView();
1047  }
1048 
1049  private:
1050 
1051  template<typename T>
1052  T aavg (T a, T b) const
1053  {
1054  return 0.5*(a+b);
1055  }
1056 
1057  template<typename T>
1058  T havg (T a, T b) const
1059  {
1060  T eps = 1E-30;
1061  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
1062  }
1063  };
1064 
1065 
1066  }
1067 }
1068 
1069 #endif // DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
const IG & ig
Definition: constraints.hh:149
const Entity & e
Definition: localfunctionspace.hh:123
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
IntersectionType
Enum describing the type of an intersection.
Definition: intersectiontype.hh:15
traits class holding the function signature, same as in local function
Definition: function.hh:183
T Traits
Export type traits.
Definition: function.hh:193
leaf of a function tree
Definition: function.hh:302
Default flags for all local operators.
Definition: flags.hh:19
@ doAlphaVolume
Definition: flags.hh:59
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:157
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton()
Definition: numericaljacobianapply.hh:182
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
sparsity pattern generator
Definition: pattern.hh:14
sparsity pattern generator
Definition: pattern.hh:30
traits class for two phase parameter class
Definition: twophaseccfv.hh:24
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: twophaseccfv.hh:47
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:44
GV GridViewType
the grid view
Definition: twophaseccfv.hh:26
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:35
@ dimDomain
dimension of the domain
Definition: twophaseccfv.hh:31
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:50
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: twophaseccfv.hh:41
typename GV::Intersection IntersectionType
Definition: twophaseccfv.hh:54
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: twophaseccfv.hh:38
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:53
base class for parameter class
Definition: twophaseccfv.hh:70
int bc_g(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
gas phase boundary condition type
Definition: twophaseccfv.hh:184
Traits::RangeFieldType nu_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase molar density
Definition: twophaseccfv.hh:153
Traits::RangeFieldType pc(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_l) const
capillary pressure function
Definition: twophaseccfv.hh:83
Traits::RangeFieldType j_l(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
liquid phase Neumann boundary condition
Definition: twophaseccfv.hh:205
Traits::RangeFieldType g_l(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
liquid phase Dirichlet boundary condition
Definition: twophaseccfv.hh:191
Traits::RangeFieldType kr_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_g) const
gas phase relative permeability
Definition: twophaseccfv.hh:107
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:131
Traits::RangeFieldType nu_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_l) const
liquid phase molar density
Definition: twophaseccfv.hh:145
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:76
int bc_l(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
liquid phase boundary condition type
Definition: twophaseccfv.hh:177
Traits::RangeFieldType mu_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase viscosity
Definition: twophaseccfv.hh:123
Traits::RangeFieldType rho_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase mass density
Definition: twophaseccfv.hh:169
Traits::RangeFieldType kr_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_l) const
liquid phase relative permeability
Definition: twophaseccfv.hh:99
Traits::RangeFieldType q_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType time) const
gas phase source term
Definition: twophaseccfv.hh:227
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:137
Traits::RangeFieldType mu_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_l) const
liquid phase viscosity
Definition: twophaseccfv.hh:115
Traits::RangeFieldType q_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType time) const
liquid phase source term
Definition: twophaseccfv.hh:219
T Traits
Definition: twophaseccfv.hh:72
Traits::RangeFieldType rho_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_l) const
liquid phase mass density
Definition: twophaseccfv.hh:161
Traits::RangeFieldType s_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType pc) const
inverse capillary pressure function
Definition: twophaseccfv.hh:91
Traits::RangeFieldType j_g(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
gas phase Neumann boundary condition
Definition: twophaseccfv.hh:212
Traits::RangeFieldType g_g(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
gas phase Dirichlet boundary condition
Definition: twophaseccfv.hh:198
Definition: twophaseccfv.hh:256
@ doLambdaVolume
Definition: twophaseccfv.hh:270
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: twophaseccfv.hh:301
@ doPatternSkeleton
Definition: twophaseccfv.hh:265
@ doLambdaBoundary
Definition: twophaseccfv.hh:271
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:280
@ doPatternVolume
Definition: twophaseccfv.hh:264
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r_s) const
Definition: twophaseccfv.hh:475
TwoPhaseTwoPointFluxOperator(const TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
constructor: pass parameter object
Definition: twophaseccfv.hh:274
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:506
@ doAlphaSkeleton
Definition: twophaseccfv.hh:268
@ doAlphaBoundary
Definition: twophaseccfv.hh:269
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: twophaseccfv.hh:406
Definition: twophaseccfv.hh:544
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:571
TwoPhaseOnePointTemporalOperator(TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
Definition: twophaseccfv.hh:558
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:564
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:613
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:812
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:646
V_l(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:638
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:641
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:848
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:1044
V_g(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:873
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:881
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:876