2#ifndef DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
3#define DUNE_PDELAB_LOCALOPERATOR_ERRORINDICATORDG_HH
10#include <dune/geometry/referenceelements.hh>
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>
48 enum { dim = T::Traits::GridViewType::dimension };
50 using Real =
typename T::Traits::RangeFieldType;
51 using BCType =
typename ConvectionDiffusionBoundaryConditions::Type;
55 enum { doPatternVolume =
false };
56 enum { doPatternSkeleton =
false };
59 enum { doAlphaVolume =
true };
60 enum { doAlphaSkeleton =
true };
61 enum { doAlphaBoundary =
true };
65 ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG,
66 ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff,
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();
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());
118 const int intorder = 2 * pOrder;
119 for (
const auto &qp : quadratureRule(geo,intorder))
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());
143 evalGradient( qp.position(), cell, lfsu, x, gradu );
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>
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
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();
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);
228 const int intorder = 2*pOrder_s;
229 for (
const auto &qp : quadratureRule(geo,intorder))
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 );
248 evalFunction( iplocal_s, lfsu_s, x_s, uDG_s );
251 evalFunction( iplocal_n, lfsu_n, x_n, uDG_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>
295 void alpha_boundary (
const IG& ig,
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();
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);
340 if (bctype != ConvectionDiffusionBoundaryConditions::Dirichlet)
343 RF uh_jump_L2normSquare(0.0);
346 for (
const auto &qp : quadratureRule(geo,intorder))
349 const auto &iplocal_s = geo_in_inside.global(qp.position());
352 auto gDirichlet = param.g( cell_inside, iplocal_s );
360 evalFunction( iplocal_s, lfsu_s, x_s, uDG_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 );
382 ConvectionDiffusionDGMethod::Type method;
383 ConvectionDiffusionDGWeights::Type weights;
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++)
395 for (
int j=i+1; j<geo.corners(); j++)
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 ¶m_, 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