dune-pdelab  2.0.0
diffusionccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_DIFFUSIONCCFV_HH
3 #define DUNE_PDELAB_DIFFUSIONCCFV_HH
4 
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/static_assert.hh>
8 #include<dune/geometry/referenceelements.hh>
9 
10 #include"defaultimp.hh"
11 #include"pattern.hh"
12 #include"flags.hh"
13 #include"diffusionparam.hh"
14 
15 
16 namespace Dune {
17  namespace PDELab {
18 
19  // a local operator for solving the Laplace equation with Dirichlet boundary conditions
20  // - div (k(x) grad u) + a0*u = f in \Omega,
21  // u = g on \partial\Omega_D
22  // - k(x) grad u * \nu = j on \partial\Omega_N
23  // with cell centered finite volumes on axiparallel cube grids
24  // K : scalar permeability field
25  // A0: Helmholtz Term
26  // F : source term
27  // B : boundary condition function
28  // J : flux boundary condition
29  // G : grid function for Dirichlet boundary conditions
30  template<typename K, typename A0, typename F, typename B, typename J, typename G>
31  class DiffusionCCFV : public NumericalJacobianApplySkeleton<DiffusionCCFV<K,A0,F,B,J,G> >,
32  public NumericalJacobianApplyBoundary<DiffusionCCFV<K,A0,F,B,J,G> >,
33  public NumericalJacobianApplyVolume<DiffusionCCFV<K,A0,F,B,J,G> >,
34  public NumericalJacobianSkeleton<DiffusionCCFV<K,A0,F,B,J,G> >,
35  public NumericalJacobianBoundary<DiffusionCCFV<K,A0,F,B,J,G> >,
36  public NumericalJacobianVolume<DiffusionCCFV<K,A0,F,B,J,G> >,
37  public FullSkeletonPattern,
38  public FullVolumePattern,
40 
41  {
42  public:
43  // pattern assembly flags
44  enum { doPatternVolume = true };
45  enum { doPatternSkeleton = true };
46 
47  // residual assembly flags
48  enum { doAlphaVolume = true };
49  enum { doAlphaSkeleton = true };
50  enum { doAlphaBoundary = true };
51  enum { doLambdaVolume = false };
52  enum { doLambdaSkeleton = false };
53  enum { doLambdaBoundary = false };
54 
55  DiffusionCCFV (const K& k_, const A0& a0_, const F& f_, const B& b_, const J& j_, const G& g_)
56  : k(k_), a0(a0_), f(f_), b(b_), j(j_), g(g_)
57  {}
58 
59 
60  // volume integral depending on test and ansatz functions
61  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
62  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
63  {
64  // domain and range field type
65  typedef typename LFSU::Traits::FiniteElementType::
66  Traits::LocalBasisType::Traits::DomainFieldType DF;
67 
68  // dimensions
69  const int dim = EG::Geometry::dimension;
70 
71  // cell center
72  const Dune::FieldVector<DF,dim>&
73  inside_local = Dune::ReferenceElements<DF,dim>::general(eg.entity().type()).position(0,0);
74 
75  // evaluate Helmholtz term
76  typename A0::Traits::RangeType a0value;
77  a0.evaluate(eg.entity(),inside_local,a0value);
78 
79  r.accumulate(lfsu,0,a0value*x(lfsu,0)*eg.geometry().volume());
80 
81  // evaluate source term
82  typename F::Traits::RangeType fvalue;
83  f.evaluate(eg.entity(),inside_local,fvalue);
84 
85  r.accumulate(lfsu,0,-fvalue*eg.geometry().volume());
86  }
87 
88  // skeleton integral depending on test and ansatz functions
89  // each face is only visited ONCE!
90  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
91  void alpha_skeleton (const IG& ig,
92  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
93  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
94  R& r_s, R& r_n) const
95  {
96  // domain and range field type
97  typedef typename LFSU::Traits::FiniteElementType::
98  Traits::LocalBasisType::Traits::DomainFieldType DF;
99  typedef typename LFSU::Traits::FiniteElementType::
100  Traits::LocalBasisType::Traits::RangeFieldType RF;
101 
102  // center in face's reference element
103  const Dune::FieldVector<DF,IG::dimension-1>&
104  face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
105 
106  // face volume for integration
107  RF face_volume = ig.geometry().integrationElement(face_local)
108  *Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).volume();
109 
110  // cell centers in references elements
111  const Dune::FieldVector<DF,IG::dimension>&
112  inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
113  const Dune::FieldVector<DF,IG::dimension>&
114  outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
115 
116  // evaluate diffusion coefficient
117  typename K::Traits::RangeType k_inside, k_outside;
118  k.evaluate(*(ig.inside()),inside_local,k_inside);
119  k.evaluate(*(ig.outside()),outside_local,k_outside);
120  typename K::Traits::RangeType k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
121 
122  // cell centers in global coordinates
123  Dune::FieldVector<DF,IG::dimension>
124  inside_global = ig.inside()->geometry().global(inside_local);
125  Dune::FieldVector<DF,IG::dimension>
126  outside_global = ig.outside()->geometry().global(outside_local);
127 
128  // distance between the two cell centers
129  inside_global -= outside_global;
130  RF distance = inside_global.two_norm();
131 
132  // contribution to residual on inside element, other residual is computed by symmetric call
133  r_s.accumulate(lfsu_s,0,k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
134  r_n.accumulate(lfsu_n,0,-k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
135  }
136 
137  // skeleton integral depending on test and ansatz functions
138  // We put the Dirchlet evaluation also in the alpha term to save som e geometry evaluations
139  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
140  void alpha_boundary (const IG& ig,
141  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
142  R& r_s) const
143  {
144  // domain and range field type
145  typedef typename LFSU::Traits::FiniteElementType::
146  Traits::LocalBasisType::Traits::DomainFieldType DF;
147  typedef typename LFSU::Traits::FiniteElementType::
148  Traits::LocalBasisType::Traits::RangeFieldType RF;
149 
150  // center in face's reference element
151  const Dune::FieldVector<DF,IG::dimension-1>&
152  face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
153 
154  // face volume for integration
155  RF face_volume = ig.geometry().integrationElement(face_local)
156  *Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).volume();
157 
158  // cell center in reference element
159  const Dune::FieldVector<DF,IG::dimension>&
160  inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
161 
162  // evaluate boundary condition type
163  typename B::Traits::RangeType bctype;
164  b.evaluate(ig,face_local,bctype);
165 
167  {
168  // Dirichlet boundary
169  // distance between cell center and face center
170  Dune::FieldVector<DF,IG::dimension>
171  inside_global = ig.inside()->geometry().global(inside_local);
172  Dune::FieldVector<DF,IG::dimension>
173  outside_global = ig.geometry().global(face_local);
174  inside_global -= outside_global;
175  RF distance = inside_global.two_norm();
176 
177  // evaluate diffusion coefficient
178  typename K::Traits::RangeType k_inside;
179  k.evaluate(*(ig.inside()),inside_local,k_inside);
180 
181  // evaluate boundary condition function
182  typename G::Traits::DomainType x = ig.geometryInInside().global(face_local);
183  typename G::Traits::RangeType y;
184  g.evaluate(*(ig.inside()),x,y);
185 
186  // contribution to residual on inside element
187  r_s.accumulate(lfsu_s,0,k_inside*(x_s(lfsu_s,0)-y[0])*face_volume/distance);
188  }
189  else // if (DiffusionBoundaryCondition::isNeumann(bctype))
190  {
191  // Neumann boundary
192  // evaluate flux boundary condition
193 
194  //evaluate boundary function
195  typename J::Traits::DomainType x = ig.geometryInInside().global(face_local);
196  typename J::Traits::RangeType jvalue;
197  j.evaluate(*(ig.inside()),x,jvalue);
198 
199  // contribution to residual on inside element
200  r_s.accumulate(lfsu_s,0,jvalue*face_volume);
201  }
202  }
203 
204  private:
205  const K& k;
206  const A0& a0;
207  const F& f;
208  const B& b;
209  const J& j;
210  const G& g;
211  };
212 
214  } // namespace PDELab
215 } // namespace Dune
216 
217 #endif
DiffusionCCFV(const K &k_, const A0 &a0_, const F &f_, const B &b_, const J &j_, const G &g_)
Definition: diffusionccfv.hh:55
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: diffusionccfv.hh:91
sparsity pattern generator
Definition: pattern.hh:30
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: diffusionccfv.hh:140
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusionccfv.hh:31
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
Definition: diffusionccfv.hh:48
static const int dim
Definition: adaptivity.hh:82
Definition: diffusionccfv.hh:50
static bool isDirichlet(Type i)
Test for Dirichlet boundary condition.
Definition: diffusionparam.hh:25
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
const IG & ig
Definition: common/constraints.hh:146
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionccfv.hh:62
Definition: diffusionccfv.hh:44
Definition: diffusionccfv.hh:53
Definition: diffusionccfv.hh:45
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
Definition: diffusionccfv.hh:49
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Definition: diffusionccfv.hh:51
Definition: diffusionccfv.hh:52