2 #ifndef DUNE_PDELAB_ERRORINDICATORDG_HH
3 #define DUNE_PDELAB_ERRORINDICATORDG_HH
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>
52 enum { dim = T::Traits::GridViewType::dimension };
54 typedef typename T::Traits::RangeFieldType Real;
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
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;
93 const int dim = EG::Geometry::dimension;
96 const int pOrder = lfsu.finiteElement().localBasis().order();
97 const int intorder = 2 * pOrder;
99 Dune::GeometryType gt = eg.geometry().type();
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 );
108 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
112 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
115 std::vector<RangeType> phi(lfsu.size());
116 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
120 for (size_type i=0; i<lfsu.size(); i++)
121 u += x(lfsu,i)*phi[i];
124 typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
127 typename T::Traits::RangeFieldType
f = param.f(eg.entity(),it->position());
130 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
133 typename T::Traits::RangeType beta
134 = param.b(eg.entity(),it->position());
137 Dune::FieldVector<RF,dim> gradu(0.0);
138 evalGradient( it->position(), eg.entity(), lfsu, x, gradu );
140 RF square = f - (beta*gradu) - c*u;
142 sum += square * factor;
146 DF h_T = diameter(eg.geometry());
149 RF eta_RK = h_T * h_T / pOrder / pOrder / epsilon * sum;
150 r.accumulate( lfsv, 0, eta_RK );
156 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
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
163 typedef typename LFSU::Traits::FiniteElementType::
164 Traits::LocalBasisType::Traits::DomainFieldType DF;
165 typedef typename LFSU::Traits::FiniteElementType::
166 Traits::LocalBasisType::Traits::RangeFieldType RF;
169 const int dim = IG::dimension;
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);
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 );
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 );
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);
192 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
194 RF flux_jump_L2normSquare(0.0);
195 RF uh_jump_L2normSquare(0.0);
198 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
201 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
202 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
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);
208 Dune::FieldVector<RF,dim> An_F_s;
210 Dune::FieldVector<RF,dim> An_F_n;
229 Dune::FieldVector<RF,dim> gradu_s(0.0);
230 evalGradient( iplocal_s, *(ig.inside()), lfsu_s, x_s, gradu_s );
232 Dune::FieldVector<RF,dim> gradu_n(0.0);
233 evalGradient( iplocal_n, *(ig.outside()), lfsu_n, x_n, gradu_n );
237 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
240 RF flux_jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
241 flux_jump_L2normSquare += flux_jump * flux_jump * factor;
244 RF jump_uDG = (uDG_s - uDG_n);
245 uh_jump_L2normSquare += jump_uDG * jump_uDG * factor ;
249 DF h_face = diameter(ig.geometry());
252 RF eta_Ek_s = 0.5 * h_face * flux_jump_L2normSquare;
253 RF eta_Ek_n = eta_Ek_s;
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;
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 );
268 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
270 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
274 typedef typename LFSU::Traits::FiniteElementType::
275 Traits::LocalBasisType::Traits::DomainFieldType DF;
276 typedef typename LFSU::Traits::FiniteElementType::
277 Traits::LocalBasisType::Traits::RangeFieldType RF;
280 const int dim = IG::dimension;
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);
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 );
292 const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
293 const int intorder = 2 * pOrder_s;
295 Dune::GeometryType gtface = ig.geometryInInside().type();
296 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
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);
306 RF uh_jump_L2normSquare(0.0);
307 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
310 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
313 RF gDirichlet = param.g( *(ig.inside()), iplocal_s );
322 RF jump_uDG = uDG_s - gDirichlet;
325 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
327 uh_jump_L2normSquare += jump_uDG * jump_uDG * factor;
332 DF h_face = diameter(ig.geometry());
335 RF eta_Jk_s = (gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
337 r_s.accumulate( lfsv_s, 0, eta_Jk_s );
347 typename GEO::ctype diameter (
const GEO& geo)
const
349 typedef typename GEO::ctype DF;
351 const int dim = GEO::coorddimension;
352 for (
int i=0; i<geo.corners(); i++)
354 Dune::FieldVector<DF,dim> xi = geo.corner(i);
355 for (
int j=i+1; j<geo.corners(); j++)
357 Dune::FieldVector<DF,dim> xj = geo.corner(j);
359 hmax = std::max(hmax,xj.two_norm());
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
Definition: errorindicatordg.hh:64
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 ¶m_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff, Real gamma_=0.0)
constructor: pass parameter object
Definition: errorindicatordg.hh:68
Definition: errorindicatordg.hh:63
const EG & eg
Definition: common/constraints.hh:277
Definition: errorindicatordg.hh:65
Type
Definition: convectiondiffusiondg.hh:35
Definition: convectiondiffusionparameter.hh:68
Definition: errorindicatordg.hh:60
Definition: errorindicatordg.hh:59
Type
Definition: convectiondiffusionparameter.hh:68