dune-pdelab  2.7-git
errorindicatordg.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
4 
5 #include <algorithm>
6 #include <vector>
7 
8 #include <dune/common/fvector.hh>
9 
10 #include <dune/geometry/referenceelements.hh>
11 
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  using Real = typename T::Traits::RangeFieldType;
51  using BCType = typename ConvectionDiffusionBoundaryConditions::Type;
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  // define types
80  using RF = typename LFSU::Traits::FiniteElementType::
81  Traits::LocalBasisType::Traits::RangeFieldType;
82  using RangeType = typename LFSU::Traits::FiniteElementType::
83  Traits::LocalBasisType::Traits::RangeType;
84  using size_type = typename LFSU::Traits::SizeType;
85 
86  // dimensions
87  const int dim = EG::Geometry::mydimension;
88 
89  // Reference to cell
90  const auto& cell = eg.entity();
91 
92  // Get geometry
93  auto geo = eg.geometry();
94 
95  // evaluate diffusion tensor at cell center, assume it is constant over elements
96  auto ref_el = referenceElement(geo);
97  auto localcenter = ref_el.position(0,0);
98  auto A = param.A(cell,localcenter);
99 
100  static_assert(dim == 2 || dim == 3,
101  "The computation of epsilon looks very "
102  "much like it will only work in 2D or 3D. If you think "
103  "otherwise, replace this static assert with a comment "
104  "that explains why. --Jö");
105  auto epsilon = std::min( A[0][0], A[1][1]);
106  if( dim>2 ) epsilon = std::min( A[2][2], epsilon );
107 
108  // select quadrature rule
109  // pOrder is constant on all grid elements (h-adaptive scheme).
110  const int pOrder = lfsu.finiteElement().localBasis().order();
111 
112  // Initialize vectors outside for loop
113  std::vector<RangeType> phi(lfsu.size());
114  Dune::FieldVector<RF,dim> gradu(0.0);
115 
116  // loop over quadrature points
117  RF sum(0.0);
118  const int intorder = 2 * pOrder;
119  for (const auto &qp : quadratureRule(geo,intorder))
120  {
121  // evaluate basis functions
122  lfsu.finiteElement().localBasis().evaluateFunction(qp.position(),phi);
123 
124  // evaluate u
125  RF u=0.0;
126  for (size_type i=0; i<lfsu.size(); i++)
127  u += x(lfsu,i)*phi[i];
128 
129  // evaluate reaction term
130  auto c = param.c(cell,qp.position());
131 
132  // evaluate right hand side parameter function
133  auto f = param.f(cell,qp.position());
134 
135  // evaluate convection term
136  auto beta = param.b(cell,qp.position());
137 
138 
139  /**********************/
140  /* Evaluate Gradients */
141  /**********************/
142  gradu = 0.0;
143  evalGradient( qp.position(), cell, lfsu, x, gradu );
144 
145 
146  // integrate f^2
147  auto factor = qp.weight() * geo.integrationElement(qp.position());
148 
149  auto square = f - (beta*gradu) - c*u; // + eps * Laplacian_u (TODO for pMax=2)
150  square *= square;
151  sum += square * factor;
152  }
153 
154  // accumulate cell indicator
155  auto h_T = diameter(geo);
156 
157  // L.Zhu: First term, interior residual squared
158  auto eta_RK = h_T * h_T / pOrder / pOrder / epsilon * sum;
159 
160  // add contributions
161  r.accumulate( lfsv, 0, eta_RK );
162  }
163 
164 
165  // skeleton integral depending on test and ansatz functions
166  // each face is only visited ONCE!
167  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
168  void alpha_skeleton (const IG& ig,
169  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
170  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
171  R& r_s, R& r_n) const
172  {
173  // Define types
174  using RF = typename LFSU::Traits::FiniteElementType::
175  Traits::LocalBasisType::Traits::RangeFieldType;
176 
177  // Dimension
178  const int dim = IG::Entity::dimension;
179 
180  // Refererences to inside and outside cells
181  const auto& cell_inside = ig.inside();
182  const auto& cell_outside = ig.outside();
183 
184  // Get geometries
185  auto geo = ig.geometry();
186  auto geo_inside = cell_inside.geometry();
187  auto geo_outside = cell_outside.geometry();
188 
189  // Get geometry of intersection in local coordinates of inside_cell and outside_cell
190  auto geo_in_inside = ig.geometryInInside();
191  auto geo_in_outside = ig.geometryInOutside();
192 
193 
194  // Evaluate permeability tensors
195  auto ref_el_inside = referenceElement(geo_inside);
196  auto ref_el_outside = referenceElement(geo_outside);
197  auto local_inside = ref_el_inside.position(0,0);
198  auto local_outside = ref_el_outside.position(0,0);
199  auto A_s = param.A(cell_inside,local_inside);
200  auto A_n = param.A(cell_outside,local_outside);
201 
202  static_assert(dim == 2 || dim == 3,
203  "The computation of epsilon_s and epsilon_n looks very "
204  "much like it will only work in 2D or 3D. If you think "
205  "otherwise, replace this static assert with a comment "
206  "that explains why. --Jö");
207 
208  auto epsilon_s = std::min( A_s[0][0], A_s[1][1]);
209  if( dim>2 ) epsilon_s = std::min( A_s[2][2], epsilon_s );
210 
211  auto epsilon_n = std::min( A_n[0][0], A_n[1][1]);
212  if( dim>2 ) epsilon_n = std::min( A_n[2][2], epsilon_n );
213 
214  const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
215 
216  auto n_F = ig.centerUnitOuterNormal();
217 
218  RF flux_jump_L2normSquare(0.0);
219  RF uh_jump_L2normSquare(0.0);
220 
221  // Declare vectors outside for loop
222  Dune::FieldVector<RF,dim> An_F_s;
223  Dune::FieldVector<RF,dim> An_F_n;
224  Dune::FieldVector<RF,dim> gradu_s;
225  Dune::FieldVector<RF,dim> gradu_n;
226 
227  // loop over quadrature points and integrate normal flux
228  const int intorder = 2*pOrder_s;
229  for (const auto &qp : quadratureRule(geo,intorder))
230  {
231  // position of quadrature point in local coordinates of elements
232  const auto &iplocal_s = geo_in_inside .global(qp.position());
233  const auto &iplocal_n = geo_in_outside.global(qp.position());
234 
235  // Diffusion tensor at quadrature point
236  A_s = param.A( cell_inside, iplocal_s );
237  A_n = param.A( cell_outside, iplocal_n );
238 
239  A_s.mv(n_F,An_F_s);
240  A_n.mv(n_F,An_F_n);
241 
242  /**********************/
243  /* Evaluate Functions */
244  /**********************/
245 
246  // evaluate uDG_s, uDG_n
247  RF uDG_s=0.0;
248  evalFunction( iplocal_s, lfsu_s, x_s, uDG_s );
249 
250  RF uDG_n=0.0;
251  evalFunction( iplocal_n, lfsu_n, x_n, uDG_n );
252 
253 
254  /**********************/
255  /* Evaluate Gradients */
256  /**********************/
257  gradu_s = 0.0;
258  evalGradient( iplocal_s, cell_inside, lfsu_s, x_s, gradu_s );
259  gradu_n = 0.0;
260  evalGradient( iplocal_n, cell_outside, lfsu_n, x_n, gradu_n );
261 
262 
263  // integrate
264  auto factor = qp.weight() * geo.integrationElement(qp.position());
265 
266  // evaluate flux jump term
267  auto flux_jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
268  flux_jump_L2normSquare += flux_jump * flux_jump * factor;
269 
270  // evaluate jump term
271  auto jump_uDG = (uDG_s - uDG_n);
272  uh_jump_L2normSquare += jump_uDG * jump_uDG * factor ;
273  }
274 
275  // accumulate indicator
276  auto h_face = diameter(geo);
277 
278  // L.Zhu: second term, edge residual
279  auto eta_Ek_s = 0.5 * h_face * flux_jump_L2normSquare;
280  auto eta_Ek_n = eta_Ek_s; // equal on both sides of the intersection!
281 
282  // L.Zhu: third term, edge jumps
283  auto eta_Jk_s = 0.5 * ( gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
284  auto eta_Jk_n = 0.5 * ( gamma / h_face + h_face / epsilon_n) * uh_jump_L2normSquare;
285 
286  // add contributions from both sides of the intersection
287  r_s.accumulate( lfsv_s, 0, eta_Ek_s + eta_Jk_s );
288  r_n.accumulate( lfsv_n, 0, eta_Ek_n + eta_Jk_n );
289  }
290 
291 
292  // boundary integral depending on test and ansatz functions
293  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
294  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
295  void alpha_boundary (const IG& ig,
296  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
297  R& r_s) const
298  {
299  // define types
300  using RF = typename LFSU::Traits::FiniteElementType::
301  Traits::LocalBasisType::Traits::RangeFieldType;
302 
303  // dimensions
304  const int dim = IG::Entity::dimension;
305 
306  // Reference to inside cell
307  const auto& cell_inside = ig.inside();
308 
309  // Get geometries
310  auto geo = ig.geometry();
311  auto geo_inside = cell_inside.geometry();
312 
313  // Get geometry of intersection in local coordinates of inside_cell
314  auto geo_in_inside = ig.geometryInInside();
315 
316  // reference elements
317  auto ref_el = referenceElement(geo);
318  auto ref_el_inside = referenceElement(geo_inside);
319 
320  // evaluate permeability tensors
321  auto local_inside = ref_el_inside.position(0,0);
322  auto A_s = param.A(cell_inside,local_inside);
323 
324  static_assert(dim == 2 || dim == 3,
325  "The computation of epsilon_s looks very "
326  "much like it will only work in 2D or 3D. If you think "
327  "otherwise, replace this static assert with a comment "
328  "that explains why. --Jö");
329 
330  auto epsilon_s = std::min( A_s[0][0], A_s[1][1]);
331  if( dim>2 ) epsilon_s = std::min( A_s[2][2], epsilon_s );
332 
333  // select quadrature rule
334  const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
335  const int intorder = 2*pOrder_s;
336 
337  // evaluate boundary condition
338  auto face_local = ref_el.position(0,0);
339  auto bctype = param.bctype(ig.intersection(),face_local);
341  return;
342 
343  RF uh_jump_L2normSquare(0.0);
344 
345  // loop over quadrature points and integrate normal flux
346  for (const auto &qp : quadratureRule(geo,intorder))
347  {
348  // position of quadrature point in local coordinates of elements
349  const auto &iplocal_s = geo_in_inside.global(qp.position());
350 
351  // evaluate Dirichlet boundary condition
352  auto gDirichlet = param.g( cell_inside, iplocal_s );
353 
354  /**********************/
355  /* Evaluate Functions */
356  /**********************/
357 
358  // evaluate uDG_s
359  RF uDG_s=0.0;
360  evalFunction( iplocal_s, lfsu_s, x_s, uDG_s );
361 
362  // integrate
363  auto factor = qp.weight() * geo.integrationElement(qp.position());
364 
365  // evaluate jump term
366  auto jump_uDG = uDG_s - gDirichlet;
367  uh_jump_L2normSquare += jump_uDG * jump_uDG * factor;
368  }
369 
370  // accumulate indicator
371  auto h_face = diameter(geo);
372 
373  // L.Zhu: third term, edge jumps on the Dirichlet boundary
374  auto eta_Jk_s = (gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
375 
376  // add contributions
377  r_s.accumulate( lfsv_s, 0, eta_Jk_s ); // boundary edge
378  }
379 
380  private:
381  T& param; // two phase parameter class
384  Real gamma; // interior penalty parameter, same as alpha in class TransportOperatorDG
385 
386  template<class GEO>
387  typename GEO::ctype diameter (const GEO& geo) const
388  {
389  using DF = typename GEO::ctype;
390  DF hmax = -1.0E00;
391  const int dim = GEO::coorddimension;
392  for (int i=0; i<geo.corners(); i++)
393  {
394  Dune::FieldVector<DF,dim> xi = geo.corner(i);
395  for (int j=i+1; j<geo.corners(); j++)
396  {
397  Dune::FieldVector<DF,dim> xj = geo.corner(j);
398  xj -= xi;
399  hmax = std::max(hmax,xj.two_norm());
400  }
401  }
402  return hmax;
403  }
404 
405  };
406 
407 
408  }
409 }
410 #endif // DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void evalFunction(const DomainType &location, const LFS &lfs, const LV &xlocal, RangeFieldType &valU)
Definition: eval.hh:20
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
void evalGradient(const DomainType &location, const EG &eg, const LFS &lfs, const LV &xlocal, RangeType &gradu)
Definition: eval.hh:47
Type
Definition: convectiondiffusiondg.hh:31
@ NIPG
Definition: convectiondiffusiondg.hh:31
Type
Definition: convectiondiffusiondg.hh:36
@ weightsOff
Definition: convectiondiffusiondg.hh:36
Type
Definition: convectiondiffusionparameter.hh:113
@ Dirichlet
Definition: convectiondiffusionparameter.hh:113
a local operator for residual-based error estimation
Definition: errorindicatordg.hh:47
@ doAlphaBoundary
Definition: errorindicatordg.hh:61
@ doAlphaVolume
Definition: errorindicatordg.hh:59
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:168
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:295
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: errorindicatordg.hh:77
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
@ doAlphaSkeleton
Definition: errorindicatordg.hh:60
@ doPatternSkeleton
Definition: errorindicatordg.hh:56
@ doPatternVolume
Definition: errorindicatordg.hh:55
Default flags for all local operators.
Definition: flags.hh:19