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