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):
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?
Thanks in advance!
Best regards,
Rob
Attachment: amr-diffu-reac-2D.py