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