Residual-type a posteriori error estimator

Hello,

Assuming that I have the 2D Poisson problem with pure Neumann b.c., is there a way to obtain a residual type estimator for H1
norm, in order to use it in an adaptive algorithm? Namely, I would like to compute the following quantity for every element K
of the mesh,

$${\eta }_{1,K}={\{h_K^2{\Vert f+\mathrm{\Delta }{u}_h\Vert }_K^2+\frac{1}{2}\sum _{S\in \partial K\cap S_{int}}{h}_S{\Vert [\mathrm{\nabla }{u}_h\cdot \overrightarrow{n}]\Vert }_S^2+\sum _{S\in \partial K\cap \partial \mathrm{\Omega }}{h}_S{\Vert g-\mathrm{\nabla }{u}_h\cdot \overrightarrow{n}\Vert }_S^2\}}^{\frac{1}{2}}.$$

Here $f$ and $g$ are the right hand side on the domain and the boundary, respectively, while $S$ is an edge and $S_{int}$ is the set of all  interior edges.

Best regards,
Vex

Hi Vex,
a residual error estimator for a Dirichlet Poisson problem is given by

h = specialcf.mesh_size
n = specialcf.normal(2)

resT = h*h*(Trace(gfu.Operator("hesse"))+f)**2*dx
resE = h*( (grad(gfu)-grad(gfu).Other())*n )**2 *dx(element_vb=BND)
eta1 = Integrate(resT+resE, mesh, element_wise=True)


To adapt to a Neumann problem I think you have to use some kind of indicator function to identify the boundary edges

indBnd = GridFunction(FacetFESpace(mesh,order=0, dirichlet=".*", autoupdate=True), autoupdate=True)
indBnd.Set(1, BND)
resEbnd = indBnd* h*(grad(gfu)*n-g )**2 *dx(element_vb=BND)


indBnd is 1 at the boundary edges and 0 at the inner edges.

Best,
Michael