DUNE PDELab (git)

errorindicatordg.hh
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
9
10#include <dune/geometry/referenceelements.hh>
11
12#include <dune/pdelab/common/quadraturerules.hh>
13#include <dune/pdelab/localoperator/convectiondiffusiondg.hh>
14#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
15#include <dune/pdelab/localoperator/eval.hh>
16#include <dune/pdelab/localoperator/flags.hh>
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
26namespace 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
65 ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG,
66 ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff,
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
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);
340 if (bctype != ConvectionDiffusionBoundaryConditions::Dirichlet)
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
382 ConvectionDiffusionDGMethod::Type method;
383 ConvectionDiffusionDGWeights::Type weights;
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
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
vector space out of a tensor product of fields.
Definition: fvector.hh:91
a local operator for residual-based error estimation
Definition: errorindicatordg.hh:47
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
Default flags for all local operators.
Definition: flags.hh:19
Implements a vector constructed from a given type representing a field and a compile-time given size.
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)