2 #ifndef DUNE_PDELAB_ERRORINDICATORDG_HH
3 #define DUNE_PDELAB_ERRORINDICATORDG_HH
8 #include <dune/common/fvector.hh>
10 #include <dune/geometry/quadraturerules.hh>
11 #include <dune/geometry/referenceelements.hh>
48 enum { dim = T::Traits::GridViewType::dimension };
50 typedef typename T::Traits::RangeFieldType Real;
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
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;
89 const int dim = EG::Geometry::dimension;
92 auto cell = eg.entity();
93 auto geo = eg.geometry();
95 const auto > = geo.type();
97 const auto &ref = Dune::ReferenceElements<DF,dim>::general(gt);
99 const auto &localcenter = ref.position(0,0);
102 auto A = param.A(cell,localcenter);
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 );
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);
119 std::vector<RangeType> phi(lfsu.size());
122 for (
const auto &qp : rule)
125 lfsu.finiteElement().localBasis().evaluateFunction(qp.position(),phi);
129 for (size_type i=0; i<lfsu.size(); i++)
130 u += x(lfsu,i)*phi[i];
133 auto c = param.c(cell,qp.position());
136 auto f = param.f(cell,qp.position());
139 auto beta = param.b(cell,qp.position());
146 Dune::FieldVector<RF,dim> gradu(0.0);
151 RF factor = qp.weight() * geo.integrationElement(qp.position());
153 RF square = f - (beta*gradu) - c*u;
155 sum += square * factor;
159 DF h_T = diameter(geo);
162 RF eta_RK = h_T * h_T / pOrder / pOrder / epsilon * sum;
165 r.accumulate( lfsv, 0, eta_RK );
171 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
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
178 typedef typename LFSU::Traits::FiniteElementType::
179 Traits::LocalBasisType::Traits::DomainFieldType DF;
180 typedef typename LFSU::Traits::FiniteElementType::
181 Traits::LocalBasisType::Traits::RangeFieldType RF;
184 const int dim = IG::dimension;
187 auto cell_inside = ig.inside();
188 auto cell_outside = ig.outside();
190 auto geo = ig.geometry();
191 auto geo_in_inside = ig.geometryInInside();
192 auto geo_in_outside = ig.geometryInOutside();
194 const auto >face = geo.type();
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());
201 const auto &inside_local = insideRef.position(0,0);
202 const auto &outside_local = outsideRef.position(0,0);
205 auto A_s = param.A(cell_inside,inside_local);
206 auto A_n = param.A(cell_outside,outside_local);
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 );
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 );
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);
224 const auto &n_F = ig.centerUnitOuterNormal();
226 RF flux_jump_L2normSquare(0.0);
227 RF uh_jump_L2normSquare(0.0);
230 for (
const auto &qp : rule)
233 const auto &iplocal_s = geo_in_inside .global(qp.position());
234 const auto &iplocal_n = geo_in_outside.global(qp.position());
237 A_s = param.A( cell_inside, iplocal_s );
238 A_n = param.A( cell_outside, iplocal_n );
240 Dune::FieldVector<RF,dim> An_F_s;
242 Dune::FieldVector<RF,dim> An_F_n;
261 Dune::FieldVector<RF,dim> gradu_s(0.0);
262 evalGradient( iplocal_s, cell_inside, lfsu_s, x_s, gradu_s );
264 Dune::FieldVector<RF,dim> gradu_n(0.0);
265 evalGradient( iplocal_n, cell_outside, lfsu_n, x_n, gradu_n );
269 RF factor = qp.weight() * geo.integrationElement(qp.position());
272 RF flux_jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
273 flux_jump_L2normSquare += flux_jump * flux_jump * factor;
276 RF jump_uDG = (uDG_s - uDG_n);
277 uh_jump_L2normSquare += jump_uDG * jump_uDG * factor ;
281 DF h_face = diameter(ig.geometry());
284 RF eta_Ek_s = 0.5 * h_face * flux_jump_L2normSquare;
285 RF eta_Ek_n = eta_Ek_s;
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;
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 );
299 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
301 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
305 typedef typename LFSU::Traits::FiniteElementType::
306 Traits::LocalBasisType::Traits::DomainFieldType DF;
307 typedef typename LFSU::Traits::FiniteElementType::
308 Traits::LocalBasisType::Traits::RangeFieldType RF;
311 const int dim = IG::dimension;
314 auto cell_inside = ig.inside();
316 auto geo = ig.geometry();
317 auto geo_in_inside = ig.geometryInInside();
319 const auto >face = geo.type();
322 Dune::ReferenceElements<DF,dim-1>::general(gtface);
323 const auto &insideRef =
324 Dune::ReferenceElements<DF,dim>::general(cell_inside.type());
326 const auto &inside_local = insideRef.position(0,0);
329 auto A_s = param.A(cell_inside,inside_local);
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 );
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);
345 const auto &face_local = ref.position(0,0);
346 BCType bctype = param.bctype(ig.intersection(),face_local);
350 RF uh_jump_L2normSquare(0.0);
353 for (
const auto &qp : rule)
356 const auto &iplocal_s = geo_in_inside .global(qp.position());
359 RF gDirichlet = param.g( cell_inside, iplocal_s );
370 RF factor = qp.weight() * geo.integrationElement(qp.position());
373 RF jump_uDG = uDG_s - gDirichlet;
374 uh_jump_L2normSquare += jump_uDG * jump_uDG * factor;
378 DF h_face = diameter(ig.geometry());
381 RF eta_Jk_s = (gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
384 r_s.accumulate( lfsv_s, 0, eta_Jk_s );
394 typename GEO::ctype diameter (
const GEO& geo)
const
396 typedef typename GEO::ctype DF;
398 const int dim = GEO::coorddimension;
399 for (
int i=0; i<geo.corners(); i++)
401 Dune::FieldVector<DF,dim> xi = geo.corner(i);
402 for (
int j=i+1; j<geo.corners(); j++)
404 Dune::FieldVector<DF,dim> xj = geo.corner(j);
406 hmax = std::max(hmax,xj.two_norm());
Definition: errorindicatordg.hh:55
Definition: convectiondiffusiondg.hh:34
a local operator for residual-based error estimation
Definition: errorindicatordg.hh:45
Definition: errorindicatordg.hh:60
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 ¶m_, 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
Definition: errorindicatordg.hh:61
Definition: errorindicatordg.hh:56
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
Definition: errorindicatordg.hh:59
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