# local edge estimator separately defined on Neumann boundaries and interfaces

Hi,

I was trying to implement an AMR code for a diffusion-reaction problem using residual-based estimator. The problem I’m dealing with has a non-homogeneous Neumann boundary condition. Hence the edge residual splits into two different cases, i.e., the second and the third term in the following formula of the estimator:
$$\eta_{R,K} = \left{ h_K^2{\left\lVert f_K+\Delta u_K-u_K\right\rVert}^2_K + \frac{1}{2}\sum_{E\in \mathcal{E}{K,\Omega}}h_E {\left\lVert \mathbb{J}(\mathbf{n}\cdot \nabla u_K) \right\rVert}E^2 + \sum {E\in \mathcal{E} {K,\Gamma_N}}h_E {\left\lVert g - \mathbf{n}\cdot \nabla u_K \right\rVert}E^2 \right}^{\frac{1}{2}},$$
where $$\mathcal{E} {K,\Omega}$$ represents all the interior edges and $$\mathcal{E}_{K,\Gamma_N}$$ represents all the Neumann boundary edges.

The code I have implemented so far (attached below) computes the interior jump term using (line 105):

etaEdge = Integrate(h*(n*(grad(gfu)-grad(gfu).Other()))**2*dx(element_boundary=True), mesh, element_wise=True) 
This clearly does not account for the Neumann b.c. part which I believe is the reason why the results show a large discrepancy between true error and the residual estimates.

So I was wondering how to distinguish two of them in the implementation? I tried the following

etaEdge = Integrate(h*(n*grad(gfu)-g2Bottom)**2*dx(element_boundary=True, definedon=mesh.Boundaries("neumBottom")), mesh, element_wise=True) 
in order to compute the contribution from part of the Neumann boundary. Weirdly this seems to traverse over all the elements and outputs a list of the size of the # of all elements. I thought it should generate a list in which only the elements on the specified boundary are computed?