dune-pdelab  2.0.0
diffusion.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_DIFFUSION_HH
3 #define DUNE_PDELAB_DIFFUSION_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/common/static_assert.hh>
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/referenceelements.hh>
12 #include<dune/geometry/quadraturerules.hh>
13 
14 #include"defaultimp.hh"
15 #include"pattern.hh"
16 #include"flags.hh"
17 #include"idefault.hh"
18 #include "diffusionparam.hh"
19 
20 namespace Dune {
21  namespace PDELab {
25 
38  template<typename K, typename A0, typename F, typename B, typename J>
39  class Diffusion : public NumericalJacobianApplyVolume<Diffusion<K,A0,F,B,J> >,
40  public FullVolumePattern,
43  //,public NumericalJacobianVolume<Diffusion<K,A0,F,B,J> >
44  {
45  public:
46  // pattern assembly flags
47  enum { doPatternVolume = true };
48 
49  // residual assembly flags
50  enum { doAlphaVolume = true };
51  enum { doLambdaVolume = true };
52  enum { doLambdaBoundary = true };
53 
54  Diffusion (const K& k_, const A0& a0_, const F& f_, const B& bctype_, const J& j_, int intorder_=2)
55  : k(k_), a0(a0_), f(f_), bctype(bctype_), j(j_), intorder(intorder_)
56  {}
57 
58  // volume integral depending on test and ansatz functions
59  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
60  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
61  {
62  // domain and range field type
63  typedef typename LFSU::Traits::FiniteElementType::
64  Traits::LocalBasisType::Traits::DomainFieldType DF;
65  typedef typename LFSU::Traits::FiniteElementType::
66  Traits::LocalBasisType::Traits::RangeFieldType RF;
67  typedef typename LFSU::Traits::FiniteElementType::
68  Traits::LocalBasisType::Traits::JacobianType JacobianType;
69  typedef typename LFSU::Traits::FiniteElementType::
70  Traits::LocalBasisType::Traits::RangeType RangeType;
71 
72  typedef typename LFSU::Traits::SizeType size_type;
73 
74  // dimensions
75  const int dim = EG::Geometry::dimension;
76 
77  // select quadrature rule
78  Dune::GeometryType gt = eg.geometry().type();
79  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
80 
81  // evaluate diffusion tensor at cell center, assume it is constant over elements
82  typename K::Traits::RangeType tensor(0.0);
83  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
84  k.evaluate(eg.entity(),localcenter,tensor);
85 
86  // loop over quadrature points
87  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
88  {
89  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
90  std::vector<JacobianType> js(lfsu.size());
91  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
92 
93  // transform gradient to real element
94  const typename EG::Geometry::JacobianInverseTransposed jac =
95  eg.geometry().jacobianInverseTransposed(it->position());
96  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
97  for (size_type i=0; i<lfsu.size(); i++)
98  {
99  gradphi[i] = 0.0;
100  jac.umv(js[i][0],gradphi[i]);
101  }
102 
103  // compute gradient of u
104  Dune::FieldVector<RF,dim> gradu(0.0);
105  for (size_type i=0; i<lfsu.size(); i++)
106  gradu.axpy(x[i],gradphi[i]);
107 
108  // compute K * gradient of u
109  Dune::FieldVector<RF,dim> Kgradu(0.0);
110  tensor.umv(gradu,Kgradu);
111 
112  // evaluate basis functions
113  std::vector<RangeType> phi(lfsu.size());
114  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
115 
116  // evaluate u
117  RF u=0.0;
118  for (size_type i=0; i<lfsu.size(); i++)
119  u += x[i]*phi[i];
120 
121  // evaluate Helmholtz term
122  typename A0::Traits::RangeType y;
123  a0.evaluate(eg.entity(),it->position(),y);
124 
125  // integrate (K grad u)*grad phi_i + a_0*u*phi_i
126  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
127  for (size_type i=0; i<lfsu.size(); i++)
128  r[i] += ( Kgradu*gradphi[i] + y*u*phi[i] )*factor;
129  }
130  }
131 
132  // jacobian of volume term
133  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
134  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
135  LocalMatrix<R>& mat) const
136  {
137  // domain and range field type
138  typedef typename LFSU::Traits::FiniteElementType::
139  Traits::LocalBasisType::Traits::DomainFieldType DF;
140  typedef typename LFSU::Traits::FiniteElementType::
141  Traits::LocalBasisType::Traits::RangeFieldType RF;
142  typedef typename LFSU::Traits::FiniteElementType::
143  Traits::LocalBasisType::Traits::JacobianType JacobianType;
144  typedef typename LFSU::Traits::FiniteElementType::
145  Traits::LocalBasisType::Traits::RangeType RangeType;
146  typedef typename LFSU::Traits::SizeType size_type;
147 
148  // dimensions
149  const int dim = EG::Geometry::dimension;
150 
151  // select quadrature rule
152  Dune::GeometryType gt = eg.geometry().type();
153  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
154 
155  // evaluate diffusion tensor at cell center, assume it is constant over elements
156  typename K::Traits::RangeType tensor(0.0);
157  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
158  k.evaluate(eg.entity(),localcenter,tensor);
159 
160  // loop over quadrature points
161  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
162  {
163  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
164  std::vector<JacobianType> js(lfsu.size());
165  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
166 
167  // transform gradient to real element
168  const typename EG::Geometry::JacobianInverseTransposed jac =
169  eg.geometry().jacobianInverseTransposed(it->position());
170  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
171  for (size_type i=0; i<lfsu.size(); i++)
172  {
173  gradphi[i] = 0.0;
174  jac.umv(js[i][0],gradphi[i]);
175  }
176 
177  // compute K * gradient of shape functions
178  std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
179  for (size_type i=0; i<lfsu.size(); i++)
180  tensor.mv(gradphi[i],Kgradphi[i]);
181 
182  // evaluate basis functions
183  std::vector<RangeType> phi(lfsu.size());
184  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
185 
186  // evaluate Helmholtz term
187  typename A0::Traits::RangeType y;
188  a0.evaluate(eg.entity(),it->position(),y);
189 
190  // integrate (K grad phi_j)*grad phi_i + a_0*phi_j*phi_i
191  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
192  for (size_type j=0; j<lfsu.size(); j++)
193  for (size_type i=0; i<lfsu.size(); i++)
194  mat(i,j) += ( Kgradphi[j]*gradphi[i] + y*phi[j]*phi[i] )*factor;
195  }
196  }
197 
198  // volume integral depending only on test functions
199  template<typename EG, typename LFSV, typename R>
200  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
201  {
202  // domain and range field type
203  typedef typename LFSV::Traits::FiniteElementType::
204  Traits::LocalBasisType::Traits::DomainFieldType DF;
205  typedef typename LFSV::Traits::FiniteElementType::
206  Traits::LocalBasisType::Traits::RangeFieldType RF;
207  typedef typename LFSV::Traits::FiniteElementType::
208  Traits::LocalBasisType::Traits::RangeType RangeType;
209 
210  typedef typename LFSV::Traits::SizeType size_type;
211 
212  // dimensions
213  const int dim = EG::Geometry::dimension;
214 
215  // select quadrature rule
216  Dune::GeometryType gt = eg.geometry().type();
217  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
218 
219  // loop over quadrature points
220  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
221  {
222  // evaluate shape functions
223  std::vector<RangeType> phi(lfsv.size());
224  lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
225 
226  // evaluate right hand side parameter function
227  typename F::Traits::RangeType y;
228  f.evaluate(eg.entity(),it->position(),y);
229 
230  // integrate f
231  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
232  for (size_type i=0; i<lfsv.size(); i++)
233  r[i] -= y*phi[i]*factor;
234  }
235  }
236 
237  // boundary integral independen of ansatz functions
238  template<typename IG, typename LFSV, typename R>
239  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
240  {
241  // domain and range field type
242  typedef typename LFSV::Traits::FiniteElementType::
243  Traits::LocalBasisType::Traits::DomainFieldType DF;
244  typedef typename LFSV::Traits::FiniteElementType::
245  Traits::LocalBasisType::Traits::RangeFieldType RF;
246  typedef typename LFSV::Traits::FiniteElementType::
247  Traits::LocalBasisType::Traits::RangeType RangeType;
248 
249  typedef typename LFSV::Traits::SizeType size_type;
250 
251  // dimensions
252  const int dim = IG::dimension;
253 
254  // select quadrature rule
255  Dune::GeometryType gtface = ig.geometryInInside().type();
256  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
257 
258  // loop over quadrature points and integrate normal flux
259  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
260  {
261  // evaluate boundary condition type
262  // skip rest if we are on Dirichlet boundary
263  if( bctype.isDirichlet( ig,it->position() ) )
264  continue;
265 
266  // position of quadrature point in local coordinates of element
267  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
268 
269  // evaluate test shape functions
270  std::vector<RangeType> phi(lfsv.size());
271  lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
272 
273  // evaluate flux boundary condition
274  typename J::Traits::RangeType y;
275  j.evaluate(*(ig.inside()),local,y);
276 
277  // integrate J
278  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
279  for (size_type i=0; i<lfsv.size(); i++)
280  r[i] += y*phi[i]*factor;
281  }
282  }
283 
284  private:
285  const K& k;
286  const A0& a0;
287  const F& f;
288  const B& bctype;
289  const J& j;
290  int intorder;
291  };
292 
294  } // namespace PDELab
295 } // namespace Dune
296 
297 #endif
Definition: diffusion.hh:51
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusion.hh:47
static const int dim
Definition: adaptivity.hh:82
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:239
const IG & ig
Definition: common/constraints.hh:146
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:200
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix< R > &mat) const
Definition: diffusion.hh:134
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:60
Definition: diffusion.hh:39
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: diffusion.hh:50
Diffusion(const K &k_, const A0 &a0_, const F &f_, const B &bctype_, const J &j_, int intorder_=2)
Definition: diffusion.hh:54