2 #ifndef DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
8 #include <dune/common/fvector.hh>
10 #include <dune/geometry/referenceelements.hh>
48 enum { dim = T::Traits::GridViewType::dimension };
50 using Real =
typename T::Traits::RangeFieldType;
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 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;
87 const int dim = EG::Geometry::mydimension;
90 const auto& cell = eg.entity();
93 auto geo = eg.geometry();
96 auto ref_el = referenceElement(geo);
97 auto localcenter = ref_el.position(0,0);
98 auto A = param.A(cell,localcenter);
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 );
110 const int pOrder = lfsu.finiteElement().localBasis().order();
113 std::vector<RangeType> phi(lfsu.size());
114 Dune::FieldVector<RF,dim> gradu(0.0);
118 const int intorder = 2 * pOrder;
122 lfsu.finiteElement().localBasis().evaluateFunction(qp.position(),phi);
126 for (size_type i=0; i<lfsu.size(); i++)
127 u += x(lfsu,i)*phi[i];
130 auto c = param.c(cell,qp.position());
133 auto f = param.f(cell,qp.position());
136 auto beta = param.b(cell,qp.position());
147 auto factor = qp.weight() * geo.integrationElement(qp.position());
149 auto square = f - (beta*gradu) - c*u;
151 sum += square * factor;
155 auto h_T = diameter(geo);
158 auto eta_RK = h_T * h_T / pOrder / pOrder / epsilon * sum;
161 r.accumulate( lfsv, 0, eta_RK );
167 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
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
174 using RF =
typename LFSU::Traits::FiniteElementType::
175 Traits::LocalBasisType::Traits::RangeFieldType;
178 const int dim = IG::Entity::dimension;
181 const auto& cell_inside =
ig.inside();
182 const auto& cell_outside =
ig.outside();
185 auto geo =
ig.geometry();
186 auto geo_inside = cell_inside.geometry();
187 auto geo_outside = cell_outside.geometry();
190 auto geo_in_inside =
ig.geometryInInside();
191 auto geo_in_outside =
ig.geometryInOutside();
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);
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ö");
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 );
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 );
214 const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
216 auto n_F =
ig.centerUnitOuterNormal();
218 RF flux_jump_L2normSquare(0.0);
219 RF uh_jump_L2normSquare(0.0);
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;
228 const int intorder = 2*pOrder_s;
232 const auto &iplocal_s = geo_in_inside .global(qp.position());
233 const auto &iplocal_n = geo_in_outside.global(qp.position());
236 A_s = param.A( cell_inside, iplocal_s );
237 A_n = param.A( cell_outside, iplocal_n );
258 evalGradient( iplocal_s, cell_inside, lfsu_s, x_s, gradu_s );
260 evalGradient( iplocal_n, cell_outside, lfsu_n, x_n, gradu_n );
264 auto factor = qp.weight() * geo.integrationElement(qp.position());
267 auto flux_jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
268 flux_jump_L2normSquare += flux_jump * flux_jump * factor;
271 auto jump_uDG = (uDG_s - uDG_n);
272 uh_jump_L2normSquare += jump_uDG * jump_uDG * factor ;
276 auto h_face = diameter(geo);
279 auto eta_Ek_s = 0.5 * h_face * flux_jump_L2normSquare;
280 auto eta_Ek_n = eta_Ek_s;
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;
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 );
294 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
296 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
300 using RF =
typename LFSU::Traits::FiniteElementType::
301 Traits::LocalBasisType::Traits::RangeFieldType;
304 const int dim = IG::Entity::dimension;
307 const auto& cell_inside =
ig.inside();
310 auto geo =
ig.geometry();
311 auto geo_inside = cell_inside.geometry();
314 auto geo_in_inside =
ig.geometryInInside();
317 auto ref_el = referenceElement(geo);
318 auto ref_el_inside = referenceElement(geo_inside);
321 auto local_inside = ref_el_inside.position(0,0);
322 auto A_s = param.A(cell_inside,local_inside);
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ö");
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 );
334 const int pOrder_s = lfsu_s.finiteElement().localBasis().order();
335 const int intorder = 2*pOrder_s;
338 auto face_local = ref_el.position(0,0);
339 auto bctype = param.bctype(
ig.intersection(),face_local);
343 RF uh_jump_L2normSquare(0.0);
349 const auto &iplocal_s = geo_in_inside.global(qp.position());
352 auto gDirichlet = param.g( cell_inside, iplocal_s );
363 auto factor = qp.weight() * geo.integrationElement(qp.position());
366 auto jump_uDG = uDG_s - gDirichlet;
367 uh_jump_L2normSquare += jump_uDG * jump_uDG * factor;
371 auto h_face = diameter(geo);
374 auto eta_Jk_s = (gamma / h_face + h_face / epsilon_s) * uh_jump_L2normSquare;
377 r_s.accumulate( lfsv_s, 0, eta_Jk_s );
387 typename GEO::ctype diameter (
const GEO& geo)
const
389 using DF =
typename GEO::ctype;
391 const int dim = GEO::coorddimension;
392 for (
int i=0; i<geo.corners(); i++)
394 Dune::FieldVector<DF,dim> xi = geo.corner(i);
395 for (
int j=i+1; j<geo.corners(); j++)
397 Dune::FieldVector<DF,dim> xj = geo.corner(j);
399 hmax = std::max(hmax,xj.two_norm());
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 ¶m_, 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