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