dune-pdelab  2.0.0
convectiondiffusionfem.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONFEM_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSIONFEM_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 
19 
21 
22 namespace Dune {
23  namespace PDELab {
24 
39  template<typename T, typename FiniteElementMap>
41  public Dune::PDELab::NumericalJacobianApplyVolume<ConvectionDiffusionFEM<T,FiniteElementMap> >,
42  public Dune::PDELab::NumericalJacobianApplyBoundary<ConvectionDiffusionFEM<T,FiniteElementMap> >,
43  public Dune::PDELab::NumericalJacobianVolume<ConvectionDiffusionFEM<T,FiniteElementMap> >,
44  public Dune::PDELab::NumericalJacobianBoundary<ConvectionDiffusionFEM<T,FiniteElementMap> >,
47  public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
48  {
49  public:
50  // pattern assembly flags
51  enum { doPatternVolume = true };
52 
53  // residual assembly flags
54  enum { doAlphaVolume = true };
55  enum { doAlphaBoundary = true };
56 
57  ConvectionDiffusionFEM (T& param_, int intorderadd_=0)
58  : param(param_), intorderadd(intorderadd_)
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  typedef typename LFSU::Traits::FiniteElementType::
70  Traits::LocalBasisType::Traits::RangeFieldType RF;
71  typedef typename LFSU::Traits::FiniteElementType::
72  Traits::LocalBasisType::Traits::JacobianType JacobianType;
73  typedef typename LFSU::Traits::FiniteElementType::
74  Traits::LocalBasisType::Traits::RangeType RangeType;
75 
76  typedef typename LFSU::Traits::SizeType size_type;
77 
78  // dimensions
79  const int dim = EG::Geometry::dimension;
80 
81  // select quadrature rule
82  Dune::GeometryType gt = eg.geometry().type();
83  const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
84  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
85 
86  // evaluate diffusion tensor at cell center, assume it is constant over elements
87  typename T::Traits::PermTensorType tensor;
88  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
89  tensor = param.A(eg.entity(),localcenter);
90 
91  // loop over quadrature points
92  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
93  {
94  // evaluate basis functions
95  // std::vector<RangeType> phi(lfsu.size());
96  // lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
97  const std::vector<RangeType>& phi = cache.evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
98 
99  // evaluate u
100  RF u=0.0;
101  for (size_type i=0; i<lfsu.size(); i++)
102  u += x(lfsu,i)*phi[i];
103 
104  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
105  // std::vector<JacobianType> js(lfsu.size());
106  // lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
107  const std::vector<JacobianType>& js = cache.evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
108 
109  // transform gradients of shape functions to real element
110  const typename EG::Geometry::JacobianInverseTransposed jac =
111  eg.geometry().jacobianInverseTransposed(it->position());
112  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
113  for (size_type i=0; i<lfsu.size(); i++)
114  jac.mv(js[i][0],gradphi[i]);
115 
116  // compute gradient of u
117  Dune::FieldVector<RF,dim> gradu(0.0);
118  for (size_type i=0; i<lfsu.size(); i++)
119  gradu.axpy(x(lfsu,i),gradphi[i]);
120 
121  // compute A * gradient of u
122  Dune::FieldVector<RF,dim> Agradu(0.0);
123  tensor.umv(gradu,Agradu);
124 
125  // evaluate velocity field, sink term and source term
126  typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
127  typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
128  typename T::Traits::RangeFieldType f = param.f(eg.entity(),it->position());
129 
130  // integrate (A grad u)*grad phi_i - u b*grad phi_i + c*u*phi_i
131  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
132  for (size_type i=0; i<lfsu.size(); i++)
133  r.accumulate(lfsu,i,( Agradu*gradphi[i] - u*(b*gradphi[i]) + (c*u-f)*phi[i] )*factor);
134  }
135  }
136 
137  // jacobian of volume term
138  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
139  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
140  M& mat) const
141  {
142  // domain and range field type
143  typedef typename LFSU::Traits::FiniteElementType::
144  Traits::LocalBasisType::Traits::DomainFieldType DF;
145  typedef typename LFSU::Traits::FiniteElementType::
146  Traits::LocalBasisType::Traits::RangeFieldType RF;
147  typedef typename LFSU::Traits::FiniteElementType::
148  Traits::LocalBasisType::Traits::JacobianType JacobianType;
149  typedef typename LFSU::Traits::FiniteElementType::
150  Traits::LocalBasisType::Traits::RangeType RangeType;
151  typedef typename LFSU::Traits::SizeType size_type;
152 
153  // dimensions
154  const int dim = EG::Geometry::dimension;
155 
156  // select quadrature rule
157  Dune::GeometryType gt = eg.geometry().type();
158  const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
159  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
160 
161  // evaluate diffusion tensor at cell center, assume it is constant over elements
162  typename T::Traits::PermTensorType tensor;
163  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
164  tensor = param.A(eg.entity(),localcenter);
165 
166  // loop over quadrature points
167  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
168  {
169  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
170  // std::vector<JacobianType> js(lfsu.size());
171  // lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
172  const std::vector<JacobianType>& js = cache.evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
173 
174  // transform gradient to real element
175  const typename EG::Geometry::JacobianInverseTransposed jac
176  = eg.geometry().jacobianInverseTransposed(it->position());
177  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
178  std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
179  for (size_type i=0; i<lfsu.size(); i++)
180  {
181  jac.mv(js[i][0],gradphi[i]);
182  tensor.mv(gradphi[i],Agradphi[i]);
183  }
184 
185  // evaluate basis functions
186  // std::vector<RangeType> phi(lfsu.size());
187  // lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
188  const std::vector<RangeType>& phi = cache.evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
189 
190  // evaluate velocity field, sink term and source te
191  typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
192  typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
193 
194  // integrate (A grad phi_j)*grad phi_i - phi_j b*grad phi_i + c*phi_j*phi_i
195  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
196  for (size_type j=0; j<lfsu.size(); j++)
197  for (size_type i=0; i<lfsu.size(); i++)
198  mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i]-phi[j]*(b*gradphi[i])+c*phi[j]*phi[i] )*factor);
199  }
200  }
201 
202  // boundary integral
203  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
204  void alpha_boundary (const IG& ig,
205  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
206  R& r_s) const
207  {
208  // domain and range field type
209  typedef typename LFSV::Traits::FiniteElementType::
210  Traits::LocalBasisType::Traits::DomainFieldType DF;
211  typedef typename LFSV::Traits::FiniteElementType::
212  Traits::LocalBasisType::Traits::RangeFieldType RF;
213  typedef typename LFSV::Traits::FiniteElementType::
214  Traits::LocalBasisType::Traits::RangeType RangeType;
215 
216  typedef typename LFSV::Traits::SizeType size_type;
217 
218  // dimensions
219  const int dim = IG::dimension;
220 
221  // evaluate boundary condition type
222  Dune::GeometryType gtface = ig.geometryInInside().type();
223  Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
225  bctype = param.bctype(ig.intersection(),facecenterlocal);
226 
227  // skip rest if we are on Dirichlet boundary
229 
230  // select quadrature rule
231  const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
232  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
233 
234  // loop over quadrature points and integrate normal flux
235  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
236  {
237  // position of quadrature point in local coordinates of element
238  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
239 
240  // evaluate shape functions (assume Galerkin method)
241  // std::vector<RangeType> phi(lfsu_s.size());
242  // lfsu_s.finiteElement().localBasis().evaluateFunction(local,phi);
243  const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
244 
246  {
247  // evaluate flux boundary condition
248  typename T::Traits::RangeFieldType j = param.j(ig.intersection(),it->position());
249 
250  // integrate j
251  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
252  for (size_type i=0; i<lfsu_s.size(); i++)
253  r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
254  }
255 
257  {
258  // evaluate u
259  RF u=0.0;
260  for (size_type i=0; i<lfsu_s.size(); i++)
261  u += x_s(lfsu_s,i)*phi[i];
262 
263  // evaluate velocity field and outer unit normal
264  typename T::Traits::RangeType b = param.b(*(ig.inside()),local);
265  const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(it->position());
266 
267  // evaluate outflow boundary condition
268  typename T::Traits::RangeFieldType o = param.o(ig.intersection(),it->position());
269 
270  // integrate o
271  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
272  for (size_type i=0; i<lfsu_s.size(); i++)
273  r_s.accumulate(lfsu_s,i,( (b*n)*u + o)*phi[i]*factor);
274  }
275  }
276  }
277 
278  // jacobian contribution from boundary
279  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
280  void jacobian_boundary (const IG& ig,
281  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
282  M& mat_s) const
283  {
284  // domain and range field type
285  typedef typename LFSV::Traits::FiniteElementType::
286  Traits::LocalBasisType::Traits::DomainFieldType DF;
287  typedef typename LFSV::Traits::FiniteElementType::
288  Traits::LocalBasisType::Traits::RangeFieldType RF;
289  typedef typename LFSV::Traits::FiniteElementType::
290  Traits::LocalBasisType::Traits::RangeType RangeType;
291 
292  typedef typename LFSV::Traits::SizeType size_type;
293 
294  // dimensions
295  const int dim = IG::dimension;
296 
297  // evaluate boundary condition type
298  Dune::GeometryType gtface = ig.geometryInInside().type();
299  Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
301  bctype = param.bctype(ig.intersection(),facecenterlocal);
302 
303  // skip rest if we are on Dirichlet boundary
306 
307  // select quadrature rule
308  const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
309  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
310 
311  // loop over quadrature points and integrate normal flux
312  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
313  {
314  // position of quadrature point in local coordinates of element
315  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
316 
317  // evaluate shape functions (assume Galerkin method)
318  // std::vector<RangeType> phi(lfsu_s.size());
319  // lfsu_s.finiteElement().localBasis().evaluateFunction(local,phi);
320  const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
321 
322  // evaluate velocity field and outer unit normal
323  typename T::Traits::RangeType b = param.b(*(ig.inside()),local);
324  const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(it->position());
325 
326  // integrate
327  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
328  for (size_type j=0; j<lfsu_s.size(); j++)
329  for (size_type i=0; i<lfsu_s.size(); i++)
330  mat_s.accumulate(lfsu_s,i,lfsu_s,j,(b*n)*phi[j]*phi[i]*factor);
331  }
332  }
333 
334 
336  void setTime (double t)
337  {
338  param.setTime(t);
339  }
340 
341  private:
342  T& param;
343  int intorderadd;
344  typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
346  };
347 
348 
349 
365  template<typename T>
368  {
369  enum { dim = T::Traits::GridViewType::dimension };
370 
371  typedef typename T::Traits::RangeFieldType Real;
373 
374  public:
375  // pattern assembly flags
376  enum { doPatternVolume = false };
377  enum { doPatternSkeleton = false };
378 
379  // residual assembly flags
380  enum { doAlphaVolume = true };
381  enum { doAlphaSkeleton = true };
382  enum { doAlphaBoundary = true };
383 
386  : param(param_)
387  {}
388 
389  // volume integral depending on test and ansatz functions
390  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
391  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
392  {
393  // domain and range field type
394  typedef typename LFSU::Traits::FiniteElementType::
395  Traits::LocalBasisType::Traits::DomainFieldType DF;
396  typedef typename LFSU::Traits::FiniteElementType::
397  Traits::LocalBasisType::Traits::RangeFieldType RF;
398  typedef typename LFSU::Traits::FiniteElementType::
399  Traits::LocalBasisType::Traits::RangeType RangeType;
400  typedef typename LFSU::Traits::SizeType size_type;
401 
402  // dimensions
403  const int dim = EG::Geometry::dimension;
404  const int intorder = 2*lfsu.finiteElement().localBasis().order();
405 
406  // select quadrature rule
407  Dune::GeometryType gt = eg.geometry().type();
408  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
409 
410  // loop over quadrature points
411  RF sum(0.0);
412  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
413  {
414  // evaluate basis functions
415  std::vector<RangeType> phi(lfsu.size());
416  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
417 
418  // evaluate u
419  RF u=0.0;
420  for (size_type i=0; i<lfsu.size(); i++)
421  u += x(lfsu,i)*phi[i];
422 
423  // evaluate reaction term
424  typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
425 
426  // evaluate right hand side parameter function
427  typename T::Traits::RangeFieldType f = param.f(eg.entity(),it->position());
428 
429  // integrate f^2
430  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
431  sum += (f*f-c*c*u*u)*factor;
432  }
433 
434  // accumulate cell indicator
435  DF h_T = diameter(eg.geometry());
436  r.accumulate(lfsv,0,h_T*h_T*sum);
437  }
438 
439 
440  // skeleton integral depending on test and ansatz functions
441  // each face is only visited ONCE!
442  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
443  void alpha_skeleton (const IG& ig,
444  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
445  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
446  R& r_s, R& r_n) const
447  {
448  // domain and range field type
449  typedef typename LFSU::Traits::FiniteElementType::
450  Traits::LocalBasisType::Traits::DomainFieldType DF;
451  typedef typename LFSU::Traits::FiniteElementType::
452  Traits::LocalBasisType::Traits::RangeFieldType RF;
453  typedef typename LFSU::Traits::FiniteElementType::
454  Traits::LocalBasisType::Traits::JacobianType JacobianType;
455  typedef typename LFSU::Traits::SizeType size_type;
456 
457  // dimensions
458  const int dim = IG::dimension;
459 
460  // evaluate permeability tensors
461  const Dune::FieldVector<DF,dim>&
462  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
463  const Dune::FieldVector<DF,dim>&
464  outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
465  typename T::Traits::PermTensorType A_s, A_n;
466  A_s = param.A(*(ig.inside()),inside_local);
467  A_n = param.A(*(ig.outside()),outside_local);
468 
469  // select quadrature rule
470  const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
471  Dune::GeometryType gtface = ig.geometryInInside().type();
472  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
473 
474  // transformation
475  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
476 
477  // tensor times normal
478  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
479  Dune::FieldVector<RF,dim> An_F_s;
480  A_s.mv(n_F,An_F_s);
481  Dune::FieldVector<RF,dim> An_F_n;
482  A_n.mv(n_F,An_F_n);
483 
484  // loop over quadrature points and integrate normal flux
485  RF sum(0.0);
486  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
487  {
488  // position of quadrature point in local coordinates of elements
489  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
490  Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
491 
492  // evaluate gradient of basis functions
493  std::vector<JacobianType> gradphi_s(lfsu_s.size());
494  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
495  std::vector<JacobianType> gradphi_n(lfsu_n.size());
496  lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
497 
498  // transform gradients of shape functions to real element
499  jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
500  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
501  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
502  jac = ig.outside()->geometry().jacobianInverseTransposed(iplocal_n);
503  std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
504  for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
505 
506  // compute gradient of u
507  Dune::FieldVector<RF,dim> gradu_s(0.0);
508  for (size_type i=0; i<lfsu_s.size(); i++)
509  gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
510  Dune::FieldVector<RF,dim> gradu_n(0.0);
511  for (size_type i=0; i<lfsu_n.size(); i++)
512  gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
513 
514  // integrate
515  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
516  RF jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
517  sum += 0.25*jump*jump*factor;
518  }
519 
520  // accumulate indicator
521  // DF h_T = diameter(ig.geometry());
522  DF h_T = std::max(diameter(ig.inside()->geometry()),diameter(ig.outside()->geometry()));
523  r_s.accumulate(lfsv_s,0,h_T*sum);
524  r_n.accumulate(lfsv_n,0,h_T*sum);
525  }
526 
527 
528  // boundary integral depending on test and ansatz functions
529  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
530  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
531  void alpha_boundary (const IG& ig,
532  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
533  R& r_s) const
534  {
535  // domain and range field type
536  typedef typename LFSU::Traits::FiniteElementType::
537  Traits::LocalBasisType::Traits::DomainFieldType DF;
538  typedef typename LFSU::Traits::FiniteElementType::
539  Traits::LocalBasisType::Traits::RangeFieldType RF;
540  typedef typename LFSU::Traits::FiniteElementType::
541  Traits::LocalBasisType::Traits::JacobianType JacobianType;
542  typedef typename LFSU::Traits::SizeType size_type;
543 
544  // dimensions
545  const int dim = IG::dimension;
546 
547  // evaluate permeability tensors
548  const Dune::FieldVector<DF,dim>&
549  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
550  typename T::Traits::PermTensorType A_s;
551  A_s = param.A(*(ig.inside()),inside_local);
552  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
553  Dune::FieldVector<RF,dim> An_F_s;
554  A_s.mv(n_F,An_F_s);
555 
556  // select quadrature rule
557  const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
558  Dune::GeometryType gtface = ig.geometryInInside().type();
559  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
560 
561  // transformation
562  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
563 
564  // evaluate boundary condition
565  const Dune::FieldVector<DF,dim-1>
566  face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
567  BCType bctype = param.bctype(ig.intersection(),face_local);
569  return;
570 
571  // loop over quadrature points and integrate normal flux
572  RF sum(0.0);
573  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
574  {
575  // position of quadrature point in local coordinates of elements
576  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
577 
578  // evaluate gradient of basis functions
579  std::vector<JacobianType> gradphi_s(lfsu_s.size());
580  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
581 
582  // transform gradients of shape functions to real element
583  jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
584  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
585  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
586 
587  // compute gradient of u
588  Dune::FieldVector<RF,dim> gradu_s(0.0);
589  for (size_type i=0; i<lfsu_s.size(); i++)
590  gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
591 
592  // evaluate flux boundary condition
593  RF j = param.j(ig.intersection(),it->position());
594 
595  // integrate
596  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
597  RF jump = j+(An_F_s*gradu_s);
598  sum += jump*jump*factor;
599  }
600 
601  // accumulate indicator
602  //DF h_T = diameter(ig.geometry());
603  DF h_T = diameter(ig.inside()->geometry());
604  r_s.accumulate(lfsv_s,0,h_T*sum);
605  }
606 
607  private:
608  T& param; // two phase parameter class
609 
610  template<class GEO>
611  typename GEO::ctype diameter (const GEO& geo) const
612  {
613  typedef typename GEO::ctype DF;
614  DF hmax = -1.0E00;
615  const int dim = GEO::coorddimension;
616  for (int i=0; i<geo.corners(); i++)
617  {
618  Dune::FieldVector<DF,dim> xi = geo.corner(i);
619  for (int j=i+1; j<geo.corners(); j++)
620  {
621  Dune::FieldVector<DF,dim> xj = geo.corner(j);
622  xj -= xi;
623  hmax = std::max(hmax,xj.two_norm());
624  }
625  }
626  return hmax;
627  }
628 
629  };
630 
631 
647  template<typename T>
651  {
652  enum { dim = T::Traits::GridViewType::dimension };
653 
654  typedef typename T::Traits::RangeFieldType Real;
656 
657  public:
658  // pattern assembly flags
659  enum { doPatternVolume = false };
660  enum { doPatternSkeleton = false };
661 
662  // residual assembly flags
663  enum { doAlphaVolume = true };
664 
666  // supply time step from implicit Euler scheme
667  ConvectionDiffusionTemporalResidualEstimator1 (T& param_, double time_, double dt_)
668  : param(param_), time(time_), dt(dt_), cmax(0)
669  {}
670 
671  // volume integral depending on test and ansatz functions
672  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
673  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
674  {
675  // domain and range field type
676  typedef typename LFSU::Traits::FiniteElementType::
677  Traits::LocalBasisType::Traits::DomainFieldType DF;
678  typedef typename LFSU::Traits::FiniteElementType::
679  Traits::LocalBasisType::Traits::RangeFieldType RF;
680  typedef typename LFSU::Traits::FiniteElementType::
681  Traits::LocalBasisType::Traits::RangeType RangeType;
682  typedef typename LFSU::Traits::SizeType size_type;
683 
684  // dimensions
685  const int dim = EG::Geometry::dimension;
686  const int intorder = 2*lfsu.finiteElement().localBasis().order();
687 
688  // select quadrature rule
689  Dune::GeometryType gt = eg.geometry().type();
690  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
691 
692  // loop over quadrature points
693  RF sum(0.0);
694  RF fsum_up(0.0);
695  RF fsum_mid(0.0);
696  RF fsum_down(0.0);
697  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
698  {
699  // evaluate basis functions
700  std::vector<RangeType> phi(lfsu.size());
701  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
702 
703  // evaluate u
704  RF u=0.0;
705  for (size_type i=0; i<lfsu.size(); i++)
706  u += x(lfsu,i)*phi[i];
707 
708  // integrate f^2
709  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
710  sum += u*u*factor;
711 
712  // evaluate right hand side parameter function
713  param.setTime(time);
714  typename T::Traits::RangeFieldType f_down = param.f(eg.entity(),it->position());
715  param.setTime(time+0.5*dt);
716  typename T::Traits::RangeFieldType f_mid = param.f(eg.entity(),it->position());
717  param.setTime(time+dt);
718  typename T::Traits::RangeFieldType f_up = param.f(eg.entity(),it->position());
719  RF f_average = (1.0/6.0)*f_down + (2.0/3.0)*f_mid + (1.0/6.0)*f_up;
720 
721  // integrate f-f_average
722  fsum_down += (f_down-f_average)*(f_down-f_average)*factor;
723  fsum_mid += (f_mid-f_average)*(f_mid-f_average)*factor;
724  fsum_up += (f_up-f_average)*(f_up-f_average)*factor;
725  }
726 
727  // accumulate cell indicator
728  DF h_T = diameter(eg.geometry());
729  r.accumulate(lfsv,0,(h_T*h_T/dt)*sum); // h^2*k_n||jump/k_n||^2
730  r.accumulate(lfsv,0,h_T*h_T * dt*((1.0/6.0)*fsum_down+(2.0/3.0)*fsum_mid+(1.0/6.0)*fsum_up) ); // h^2*||f-time_average(f)||^2_0_s_t
731  }
732 
733  void clearCmax ()
734  {
735  cmax = 0;
736  }
737 
738  double getCmax () const
739  {
740  return cmax;
741  }
742 
743  private:
744  T& param; // two phase parameter class
745  double time;
746  double dt;
747  mutable double cmax;
748 
749  template<class GEO>
750  typename GEO::ctype diameter (const GEO& geo) const
751  {
752  typedef typename GEO::ctype DF;
753  DF hmax = -1.0E00;
754  const int dim = GEO::coorddimension;
755  for (int i=0; i<geo.corners(); i++)
756  {
757  Dune::FieldVector<DF,dim> xi = geo.corner(i);
758  for (int j=i+1; j<geo.corners(); j++)
759  {
760  Dune::FieldVector<DF,dim> xj = geo.corner(j);
761  xj -= xi;
762  hmax = std::max(hmax,xj.two_norm());
763  }
764  }
765  return hmax;
766  }
767 
768  };
769 
770  // a functor that can be used to evaluate rhs parameter function in interpolate
771  template<typename T, typename EG>
773  {
774  public:
775  CD_RHS_LocalAdapter (const T& t_, const EG& eg_) : t(t_), eg(eg_)
776  {}
777 
778  template<typename X, typename Y>
779  inline void evaluate (const X& x, Y& y) const
780  {
781  y[0] = t.f(eg.entity(),x);
782  }
783 
784  private:
785  const T& t;
786  const EG& eg;
787  };
788 
804  template<typename T>
808  {
809  enum { dim = T::Traits::GridViewType::dimension };
810 
811  typedef typename T::Traits::RangeFieldType Real;
813 
814  public:
815  // pattern assembly flags
816  enum { doPatternVolume = false };
817  enum { doPatternSkeleton = false };
818 
819  // residual assembly flags
820  enum { doAlphaVolume = true };
821  enum { doAlphaBoundary = true };
822 
824  // supply time step from implicit Euler scheme
825  ConvectionDiffusionTemporalResidualEstimator (T& param_, double time_, double dt_)
826  : param(param_), time(time_), dt(dt_)
827  {}
828 
829  // volume integral depending on test and ansatz functions
830  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
831  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
832  {
833  // domain and range field type
834  typedef typename LFSU::Traits::FiniteElementType::
835  Traits::LocalBasisType::Traits::DomainFieldType DF;
836  typedef typename LFSU::Traits::FiniteElementType::
837  Traits::LocalBasisType::Traits::RangeFieldType RF;
838  typedef typename LFSU::Traits::FiniteElementType::
839  Traits::LocalBasisType::Traits::RangeType RangeType;
840  typedef typename LFSU::Traits::SizeType size_type;
841  typedef typename LFSU::Traits::FiniteElementType::
842  Traits::LocalBasisType::Traits::JacobianType JacobianType;
843 
844  // dimensions
845  const int dim = EG::Geometry::dimension;
846  const int intorder = 2*lfsu.finiteElement().localBasis().order();
847 
848  // select quadrature rule
849  Dune::GeometryType gt = eg.geometry().type();
850  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
851 
852  // interpolate f as finite element function to compute the gradient
853  CD_RHS_LocalAdapter<T,EG> f_adapter(param,eg);
854  std::vector<RF> f_up, f_down, f_mid;
855  param.setTime(time);
856  lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_down);
857  param.setTime(time+0.5*dt);
858  lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_mid);
859  param.setTime(time+dt);
860  lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_up);
861 
862  // loop over quadrature points
863  RF sum(0.0);
864  RF sum_grad(0.0);
865  RF fsum_grad_up(0.0);
866  RF fsum_grad_mid(0.0);
867  RF fsum_grad_down(0.0);
868  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
869  {
870  // evaluate basis functions
871  std::vector<RangeType> phi(lfsu.size());
872  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
873 
874  // evaluate u
875  RF u=0.0;
876  for (size_type i=0; i<lfsu.size(); i++)
877  u += x(lfsu,i)*phi[i];
878 
879  // integrate jump
880  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
881  sum += u*u*factor;
882 
883  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
884  std::vector<JacobianType> js(lfsu.size());
885  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
886 
887  // transform gradients of shape functions to real element
888  const typename EG::Geometry::JacobianInverseTransposed jac =
889  eg.geometry().jacobianInverseTransposed(it->position());
890  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
891  for (size_type i=0; i<lfsu.size(); i++)
892  jac.mv(js[i][0],gradphi[i]);
893 
894  // compute gradient of u
895  Dune::FieldVector<RF,dim> gradu(0.0);
896  for (size_type i=0; i<lfsu.size(); i++)
897  gradu.axpy(x(lfsu,i),gradphi[i]);
898 
899  // integrate jump of gradient
900  sum_grad += (gradu*gradu)*factor;
901 
902  // compute gradients of f
903  Dune::FieldVector<RF,dim> gradf_down(0.0);
904  for (size_type i=0; i<lfsu.size(); i++) gradf_down.axpy(f_down[i],gradphi[i]);
905  Dune::FieldVector<RF,dim> gradf_mid(0.0);
906  for (size_type i=0; i<lfsu.size(); i++) gradf_mid.axpy(f_mid[i],gradphi[i]);
907  Dune::FieldVector<RF,dim> gradf_up(0.0);
908  for (size_type i=0; i<lfsu.size(); i++) gradf_up.axpy(f_up[i],gradphi[i]);
909  Dune::FieldVector<RF,dim> gradf_average(0.0);
910  for (size_type i=0; i<lfsu.size(); i++)
911  gradf_average.axpy((1.0/6.0)*f_down[i]+(2.0/3.0)*f_mid[i]+(1.0/6.0)*f_up[i],gradphi[i]);
912 
913  // integrate grad(f-f_average)
914  gradf_down -= gradf_average;
915  fsum_grad_down += (gradf_down*gradf_down)*factor;
916  gradf_mid -= gradf_average;
917  fsum_grad_mid += (gradf_mid*gradf_mid)*factor;
918  gradf_up -= gradf_average;
919  fsum_grad_up += (gradf_up*gradf_up)*factor;
920  }
921 
922  // accumulate cell indicator
923  DF h_T = diameter(eg.geometry());
924  r.accumulate(lfsv,0,dt * sum_grad); // k_n*||grad(jump)||^2
925  r.accumulate(lfsv,0,dt*dt * dt*((1.0/6.0)*fsum_grad_down+(2.0/3.0)*fsum_grad_mid+(1.0/6.0)*fsum_grad_up)); // k_n^2*||grad(f-time_average(f))||^2_s_t
926  }
927 
928  // boundary integral depending on test and ansatz functions
929  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
930  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
931  void alpha_boundary (const IG& ig,
932  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
933  R& r_s) const
934  {
935  // domain and range field type
936  typedef typename LFSU::Traits::FiniteElementType::
937  Traits::LocalBasisType::Traits::DomainFieldType DF;
938  typedef typename LFSU::Traits::FiniteElementType::
939  Traits::LocalBasisType::Traits::RangeFieldType RF;
940 
941  // dimensions
942  const int dim = IG::dimension;
943 
944  // select quadrature rule
945  const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
946  Dune::GeometryType gtface = ig.geometryInInside().type();
947  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
948 
949  // evaluate boundary condition
950  const Dune::FieldVector<DF,dim-1>
951  face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
952  BCType bctype = param.bctype(ig.intersection(),face_local);
954  return;
955 
956  // loop over quadrature points and integrate
957  RF sum_up(0.0);
958  RF sum_down(0.0);
959  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
960  {
961  // evaluate flux boundary condition
962  param.setTime(time);
963  RF j_down = param.j(ig.intersection(),it->position());
964  param.setTime(time+0.5*dt);
965  RF j_mid = param.j(ig.intersection(),it->position());
966  param.setTime(time+dt);
967  RF j_up = param.j(ig.intersection(),it->position());
968 
969  // integrate
970  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
971  sum_down += (j_down-j_mid)*(j_down-j_mid)*factor;
972  sum_up += (j_up-j_mid)*(j_up-j_mid)*factor;
973  }
974 
975  // accumulate indicator
976  //DF h_T = diameter(ig.geometry());
977  DF h_T = diameter(ig.inside()->geometry());
978  r_s.accumulate(lfsv_s,0,(h_T+dt*dt/h_T)*dt*0.5*(sum_down+sum_up));
979  }
980 
981  private:
982  T& param; // two phase parameter class
983  double time;
984  double dt;
985 
986  template<class GEO>
987  typename GEO::ctype diameter (const GEO& geo) const
988  {
989  typedef typename GEO::ctype DF;
990  DF hmax = -1.0E00;
991  const int dim = GEO::coorddimension;
992  for (int i=0; i<geo.corners(); i++)
993  {
994  Dune::FieldVector<DF,dim> xi = geo.corner(i);
995  for (int j=i+1; j<geo.corners(); j++)
996  {
997  Dune::FieldVector<DF,dim> xj = geo.corner(j);
998  xj -= xi;
999  hmax = std::max(hmax,xj.two_norm());
1000  }
1001  }
1002  return hmax;
1003  }
1004 
1005  };
1006 
1007 
1008 
1009  }
1010 }
1011 #endif
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionfem.hh:51
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: convectiondiffusionfem.hh:443
ConvectionDiffusionTemporalResidualEstimator1(T &param_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:667
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionfem.hh:139
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:391
Definition: convectiondiffusionfem.hh:772
Definition: convectiondiffusionfem.hh:40
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_s) const
Definition: convectiondiffusionfem.hh:280
void evaluate(const X &x, Y &y) const
Definition: convectiondiffusionfem.hh:779
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:64
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
static const int dim
Definition: adaptivity.hh:82
Definition: convectiondiffusionfem.hh:648
void clearCmax()
Definition: convectiondiffusionfem.hh:733
ConvectionDiffusionFEMResidualEstimator(T &param_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:385
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
const IG & ig
Definition: common/constraints.hh:146
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:204
void setTime(double t)
set time in parameter class
Definition: convectiondiffusionfem.hh:336
ConvectionDiffusionFEM(T &param_, int intorderadd_=0)
Definition: convectiondiffusionfem.hh:57
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:931
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:531
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:16
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
const F & f
Definition: common/constraints.hh:145
Definition: convectiondiffusionfem.hh:805
Definition: convectiondiffusionfem.hh:54
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
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:673
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:831
Definition: convectiondiffusionfem.hh:55
Definition: convectiondiffusionparameter.hh:68
ConvectionDiffusionTemporalResidualEstimator(T &param_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:825
CD_RHS_LocalAdapter(const T &t_, const EG &eg_)
Definition: convectiondiffusionfem.hh:775
Type
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionparameter.hh:68
double getCmax() const
Definition: convectiondiffusionfem.hh:738
Definition: convectiondiffusionfem.hh:366