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:\[\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}\]
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):
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]]>