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