dune-pdelab  2.0.0
errorindicatordg.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_ERRORINDICATORDG_HH
3 #define DUNE_PDELAB_ERRORINDICATORDG_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 #include "convectiondiffusiondg.hh"
22 #include "eval.hh"
23 
24 // Note:
25 // The residual-based error estimator implemented here (for h-refinement only!)
26 // is taken from the PhD thesis
27 // "Robust A Posteriori Error Estimation for Discontinuous Galerkin Methods for Convection Diffusion Problems"
28 // by Liang Zhu (2010)
29 //
30 
31 namespace Dune {
32  namespace PDELab {
33 
48  template<typename T>
51  {
52  enum { dim = T::Traits::GridViewType::dimension };
53 
54  typedef typename T::Traits::RangeFieldType Real;
56 
57  public:
58  // pattern assembly flags
59  enum { doPatternVolume = false };
60  enum { doPatternSkeleton = false };
61 
62  // residual assembly flags
63  enum { doAlphaVolume = true };
64  enum { doAlphaSkeleton = true };
65  enum { doAlphaBoundary = true };
66 
71  Real gamma_=0.0
72  )
73  : param(param_),
74  method(method_),
75  weights(weights_),
76  gamma(gamma_) // interior penalty parameter, same as alpha in ConvectionDiffusionDG
77  {}
78 
79  // volume integral depending on test and ansatz functions
80  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
81  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
82  {
83  // domain and range field type
84  typedef typename LFSU::Traits::FiniteElementType::
85  Traits::LocalBasisType::Traits::DomainFieldType DF;
86  typedef typename LFSU::Traits::FiniteElementType::
87  Traits::LocalBasisType::Traits::RangeFieldType RF;
88  typedef typename LFSU::Traits::FiniteElementType::
89  Traits::LocalBasisType::Traits::RangeType RangeType;
90  typedef typename LFSU::Traits::SizeType size_type;
91 
92  // dimensions
93  const int dim = EG::Geometry::dimension;
94 
95  // pOrder is constant on all grid elements (h-adaptive scheme).
96  const int pOrder = lfsu.finiteElement().localBasis().order();
97  const int intorder = 2 * pOrder;
98 
99  Dune::GeometryType gt = eg.geometry().type();
100  // Diffusion tensor at cell center
101  typename T::Traits::PermTensorType A;
102  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
103  A = param.A(eg.entity(),localcenter);
104  RF epsilon = std::min( A[0][0], A[1][1]);
105  if( dim>2 ) epsilon = std::min( A[2][2], epsilon );
106 
107  // select quadrature rule
108  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
109 
110  // loop over quadrature points
111  RF sum(0.0);
112  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
113  {
114  // evaluate basis functions
115  std::vector<RangeType> phi(lfsu.size());
116  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
117 
118  // evaluate u
119  RF u=0.0;
120  for (size_type i=0; i<lfsu.size(); i++)
121  u += x(lfsu,i)*phi[i];
122 
123  // evaluate reaction term
124  typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
125 
126  // evaluate right hand side parameter function
127  typename T::Traits::RangeFieldType f = param.f(eg.entity(),it->position());
128 
129  // integrate f^2
130  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
131 
132  // evaluate convection term
133  typename T::Traits::RangeType beta
134  = param.b(eg.entity(),it->position());
135 
136  // compute gradient of u_h
137  Dune::FieldVector<RF,dim> gradu(0.0);
138  evalGradient( it->position(), eg.entity(), lfsu, x, gradu );
139 
140  RF square = f - (beta*gradu) - c*u; // + eps * Laplacian_u (TODO for pMax=2)
141  square *= square;
142  sum += square * factor;
143  }
144 
145  // accumulate cell indicator
146  DF h_T = diameter(eg.geometry());
147 
148  // L.Zhu: First term, interior residual squared
149  RF eta_RK = h_T * h_T / pOrder / pOrder / epsilon * sum;
150  r.accumulate( lfsv, 0, eta_RK );
151  }
152 
153 
154  // skeleton integral depending on test and ansatz functions
155  // each face is only visited ONCE!
156  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
157  void alpha_skeleton (const IG& ig,
158  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
159  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
160  R& r_s, R& r_n) const
161  {
162  // domain and range field type
163  typedef typename LFSU::Traits::FiniteElementType::
164  Traits::LocalBasisType::Traits::DomainFieldType DF;
165  typedef typename LFSU::Traits::FiniteElementType::
166  Traits::LocalBasisType::Traits::RangeFieldType RF;
167 
168  // dimensions
169  const int dim = IG::dimension;
170 
171  // evaluate permeability tensors
172  const Dune::FieldVector<DF,dim>&
173  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
174  const Dune::FieldVector<DF,dim>&
175  outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
176  typename T::Traits::PermTensorType A_s, A_n;
177  A_s = param.A(*(ig.inside()),inside_local);
178  A_n = param.A(*(ig.outside()),outside_local);
179 
180  RF epsilon_s = std::min( A_s[0][0], A_s[1][1]);
181  if( dim>2 ) epsilon_s = std::min( A_s[2][2], epsilon_s );
182 
183  RF epsilon_n = std::min( A_n[0][0], A_n[1][1]);
184  if( dim>2 ) epsilon_n = std::min( A_n[2][2], epsilon_n );
185 
186  // select quadrature rule
187  const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
188  const int intorder = 2*pOrder_s;
189  Dune::GeometryType gtface = ig.geometryInInside().type();
190  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
191 
192  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
193 
194  RF flux_jump_L2normSquare(0.0);
195  RF uh_jump_L2normSquare(0.0);
196 
197  // loop over quadrature points and integrate normal flux
198  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
199  {
200  // position of quadrature point in local coordinates of elements
201  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
202  Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
203 
204  // Diffusion tensor at quadrature point
205  typename T::Traits::PermTensorType A_s = param.A( *(ig.inside()), iplocal_s );
206  typename T::Traits::PermTensorType A_n = param.A( *(ig.outside()), iplocal_n);
207 
208  Dune::FieldVector<RF,dim> An_F_s;
209  A_s.mv(n_F,An_F_s);
210  Dune::FieldVector<RF,dim> An_F_n;
211  A_n.mv(n_F,An_F_n);
212 
213  /**********************/
214  /* Evaluate Functions */
215  /**********************/
216 
217  // evaluate uDG_s, uDG_n
218  RF uDG_s=0.0;
219  evalFunction( iplocal_s, lfsu_s, x_s, uDG_s );
220 
221  RF uDG_n=0.0;
222  evalFunction( iplocal_n, lfsu_n, x_n, uDG_n );
223 
224 
225  /**********************/
226  /* Evaluate Gradients */
227  /**********************/
228 
229  Dune::FieldVector<RF,dim> gradu_s(0.0);
230  evalGradient( iplocal_s, *(ig.inside()), lfsu_s, x_s, gradu_s );
231 
232  Dune::FieldVector<RF,dim> gradu_n(0.0);
233  evalGradient( iplocal_n, *(ig.outside()), lfsu_n, x_n, gradu_n );
234 
235 
236  // integrate
237  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
238 
239  // evaluate flux jump term
240  RF flux_jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
241  flux_jump_L2normSquare += flux_jump * flux_jump * factor;
242 
243  // evaluate jump term
244  RF jump_uDG = (uDG_s - uDG_n);
245  uh_jump_L2normSquare += jump_uDG * jump_uDG * factor ;
246  }
247 
248  // accumulate indicator
249  DF h_face = diameter(ig.geometry());
250 
251  // L.Zhu: second term, edge residual
252  RF eta_Ek_s = 0.5 * h_face * flux_jump_L2normSquare;
253  RF eta_Ek_n = eta_Ek_s; // equal on both sides of the intersection!
254 
255  // L.Zhu: third term, edge jumps
256  RF eta_Jk_s = 0.5 * ( gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
257  RF eta_Jk_n = 0.5 * ( gamma / h_face + h_face / epsilon_n) * uh_jump_L2normSquare;
258 
259  // add contributions from both sides of the intersection
260  r_s.accumulate( lfsv_s, 0, eta_Ek_s + eta_Jk_s );
261  r_n.accumulate( lfsv_n, 0, eta_Ek_n + eta_Jk_n );
262 
263  }
264 
265 
266  // boundary integral depending on test and ansatz functions
267  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
268  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
269  void alpha_boundary (const IG& ig,
270  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
271  R& r_s) const
272  {
273  // domain and range field type
274  typedef typename LFSU::Traits::FiniteElementType::
275  Traits::LocalBasisType::Traits::DomainFieldType DF;
276  typedef typename LFSU::Traits::FiniteElementType::
277  Traits::LocalBasisType::Traits::RangeFieldType RF;
278 
279  // dimensions
280  const int dim = IG::dimension;
281 
282  // evaluate permeability tensors
283  const Dune::FieldVector<DF,dim>&
284  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
285  typename T::Traits::PermTensorType A_s;
286  A_s = param.A(*(ig.inside()),inside_local);
287 
288  RF epsilon_s = std::min( A_s[0][0], A_s[1][1]);
289  if( dim>2 ) epsilon_s = std::min( A_s[2][2], epsilon_s );
290 
291  // select quadrature rule
292  const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
293  const int intorder = 2 * pOrder_s;
294 
295  Dune::GeometryType gtface = ig.geometryInInside().type();
296  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
297 
298  // evaluate boundary condition
299  const Dune::FieldVector<DF,dim-1>
300  face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
301  BCType bctype = param.bctype(ig.intersection(),face_local);
303  return;
304 
305  // loop over quadrature points and integrate normal flux
306  RF uh_jump_L2normSquare(0.0);
307  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
308  {
309  // position of quadrature point in local coordinates of elements
310  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
311 
312  // evaluate Dirichlet boundary condition
313  RF gDirichlet = param.g( *(ig.inside()), iplocal_s );
314 
315  /**********************/
316  /* Evaluate Functions */
317  /**********************/
318 
319  // evaluate uDG_s
320  RF uDG_s=0.0;
321  evalFunction( iplocal_s, lfsu_s, x_s, uDG_s );
322  RF jump_uDG = uDG_s - gDirichlet;
323 
324  // integrate
325  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
326 
327  uh_jump_L2normSquare += jump_uDG * jump_uDG * factor;
328 
329  }
330 
331  // accumulate indicator
332  DF h_face = diameter(ig.geometry());
333 
334  // L.Zhu: third term, edge jumps on the Dirichlet boundary
335  RF eta_Jk_s = (gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
336 
337  r_s.accumulate( lfsv_s, 0, eta_Jk_s ); // boundary edge
338  }
339 
340  private:
341  T& param; // two phase parameter class
344  Real gamma; // interior penalty parameter, same as alpha in class TransportOperatorDG
345 
346  template<class GEO>
347  typename GEO::ctype diameter (const GEO& geo) const
348  {
349  typedef typename GEO::ctype DF;
350  DF hmax = -1.0E00;
351  const int dim = GEO::coorddimension;
352  for (int i=0; i<geo.corners(); i++)
353  {
354  Dune::FieldVector<DF,dim> xi = geo.corner(i);
355  for (int j=i+1; j<geo.corners(); j++)
356  {
357  Dune::FieldVector<DF,dim> xj = geo.corner(j);
358  xj -= xi;
359  hmax = std::max(hmax,xj.two_norm());
360  }
361  }
362  return hmax;
363  }
364 
365  };
366 
367 
368  }
369 }
370 #endif
Definition: convectiondiffusiondg.hh:35
Definition: errorindicatordg.hh:49
Default flags for all local operators.
Definition: flags.hh:18
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: errorindicatordg.hh:269
void evalGradient(const DomainType &location, const EG &eg, const LFS &lfs, const LV &xlocal, RangeType &gradu)
Definition: eval.hh:47
Type
Definition: convectiondiffusiondg.hh:40
void evalFunction(const DomainType &location, const LFS &lfs, const LV &xlocal, RangeFieldType &valU)
Definition: eval.hh:20
const IG & ig
Definition: common/constraints.hh:146
Definition: convectiondiffusiondg.hh:40
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: errorindicatordg.hh:81
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: errorindicatordg.hh:157
const F & f
Definition: common/constraints.hh:145
ConvectionDiffusionDG_ErrorIndicator(T &param_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff, Real gamma_=0.0)
constructor: pass parameter object
Definition: errorindicatordg.hh:68
const EG & eg
Definition: common/constraints.hh:277
Type
Definition: convectiondiffusiondg.hh:35
Definition: convectiondiffusionparameter.hh:68
Type
Definition: convectiondiffusionparameter.hh:68