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:
[tex]\begin{equation}
\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}},
\end{equation}[/tex]
where [tex]\mathcal{E}_{K,\Omega}[/tex] represents all the interior edges and [tex]\mathcal{E}_{K,\Gamma_N}[/tex] represents all the Neumann boundary edges.

The code I have implemented so far (attached below) computes the interior jump term using (line 105): 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 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?

Thanks in advance!

Best regards,
Rob