Computation of the boundary terms on a residual type estimator

Hello,

I am trying to create an adaptive algorithm by means of a simple residual type estimator
on H1 norm, for an elliptic PDE with Robin boundary conditions,

\[
\begin{aligned}
-\Delta u&=f \ \ \textit{in}\ \ \Omega,
\\ \frac{\partial u}{\partial n}+\gamma u&=g \ \ \textit{in}\ \ \Gamma,
\end{aligned}
\]

with $\Gamma = \partial \Omega$ and and $S_h$ the set of internal interelement sides of $\mathcal{T}_h$.
The mark and refinement will be done with the following local estimator

\[
\eta_{1, K}=\left\{h_K^2\left\|f+\Delta u_h\right\|_K^2+\frac{1}{2} \sum_{S \in \partial K \cap S_h} h_S\left\|\left[ \nabla u_h\cdot \vec{n}\right]\right\|_S^2+\sum_{S\in \partial K\cap \Gamma} h_S\left\|g- \nabla u_h\cdot \vec{n} -\gamma u_h\right\|_S^2\right\}^{\frac{1}{2}}.
\]

Therefore, I need the contribution of the above terms for every single element of the triangulation. 
To the best of my knowledge, the implementation of the first and second term of the above formula will be 

Code:
#rhd is the force term f, gfu is the finite element solution and lapgfu is its Laplacian
ee=Integrate(h**2*(rhd + lapgfu)**2, mesh, VOL, element_wise=True) # first term 
ej=Integrate (1/2*h*(ni*(grad(gfu)-grad(gfu).Other()))**2*dx(element_boundary=True), mesh, element_wise=True) #second term

My question is how I can compute the contribution of the third term. For example (I attach rtest.py file) if we start with
the unit square with 34 Elements and 16 Boundary Edges, then the above 2 lines will give us 2 vectors (ee and ej)
with 34 values each, that correspond to the "interior and jump contribution" of the 34 elements. Nevertheless, for the
elements that belong to the boundary, we also need to include the final term of the estimator.

So my first question is how I can compute this final term for every boundary element. I tried some potential ways, like
the one I describe below, but Ι am not sure about them, because none of them seems to work perfectly.

The second one, is that if we use a similar technique for this term, e.g.,

Code: 
#rhb is g, that is, the right hand side of the boundary and gamma is the robin coefficient 
eb=Integrate(h*(rhb-grad(gfu)*ni-gamma*gfu)**2, mesh, BND, element_wise=True) #potential third term

then we will have the vector eb with 16 values (error estimates) that correspond to the 16 boundary edges. So my
question is how  can someone match these 16 estimates to the proper elements of the mesh in order to have the
right contribution of each element according to the local estimator.

My first thought was to make a function (see BEOEV in rtest.py) which has the mesh as an input and as the output
a list with the number of boundary edges of every element through the loop "for el in mesh.Elements():". Then
add the values of eb to ee+ej accordingly to this list. However, after some experiments, I think that the
enumeration that occurs from the "element_wise=True" of the Integrate function and the enumeration of the
loop "for el in mesh.Elements():", are not the same, hence, this method will not work.

The rtest.py file
Code:
import netgen.gui
from ngsolve import *
from netgen.geom2d import SplineGeometry
import numpy as np

hh=1/4
mesh=Mesh(unit_square.GenerateMesh(maxh=hh))
Draw(mesh)

ni = specialcf.normal(mesh.dim)
gamma=1
rd=1
fes=H1(mesh,order=rd)
u,v=fes.TnT()

#Exact Solution
uex=CoefficientFunction((1-x)*(1-y)*sin(x)*sin(y))
uexfes=GridFunction(fes)
uexfes.Set(uex)
Draw(uexfes)

#Compute the right hand side 
gradu=CoefficientFunction((uex.Diff(x),uex.Diff(y)))
grx=gradu[0]
gry=gradu[1]
lap=CoefficientFunction(grx.Diff(x)+gry.Diff(y))
rhd=-lap
rhb=gradu*ni+gamma*uex

#Variational Formulation
a=BilinearForm(fes)
a+=grad(u)*grad(v)*dx+gamma*u*v*ds
a.Assemble()
f=LinearForm(fes)
f+=rhd*v*dx+rhb*v*ds
f.Assemble()

#Solve
gfu=GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec
Draw(gfu)

#L2 space
L2S = L2(mesh, order = rd-1)
gradgfu_x = GridFunction(L2S)
gradgfu_y = GridFunction(L2S)
gradgfu_x.Set(gfu.Deriv()[0])
gradgfu_y.Set(gfu.Deriv()[1])
h = specialcf.mesh_size
lapgfu = gradgfu_x.Deriv()[0] + gradgfu_y.Deriv()[1]

#Interior contribution
ee=Integrate(h**2*(rhd + lapgfu)**2, mesh, VOL, element_wise=True)
#print(ee)
print(max(ee))
#Jump contribution
ej=Integrate (1/2*h*(ni*(grad(gfu)-grad(gfu).Other()))**2*dx(element_boundary=True), mesh, element_wise=True)
#print(ej)
print(max(ej))
#Boundary Contribution
eb=Integrate(h*(rhb-grad(gfu)*ni-gamma*gfu)**2, mesh, BND, element_wise=True)
print(max(eb))
print(len(eb))
#print(eb)


#Function BEOEV
def BEOEV(mesh):
    #[be] is a list with all boundary edges of the elements of the mesh
    be=[]
    for el in mesh.Elements(BND):
        elm=mesh[el]
        be.append(elm.edges[0])
    
    #VO1 is a list that for every element, checks if there is a boundary edge with dim num_el*num_be*num_el_edge
    VO1=[]
    for el in mesh.Elements():
        elm2=mesh[el]
        pr=1
        for ed in be:
            for i in range(0,len(elm2.edges)):
                if ed == elm2.edges:
                    VO1.append(pr)
                    pr+=1
                else:
                    VO1.append(0)
    #VO2 is list that contains the number of boundary edges of an an element
    VO2=[]
    for i in range(0,mesh.ne):
        nbe=len(be)
        rdb=mesh.dim+1
        pr=max(VO1[i*nbe*rdb:(i+1)*nbe*rdb])
        VO2.append(pr)
    return VO2
    
    
bei=BEOEV(mesh)
print(bei)
#Matching the boundary terms to the estimator
eta=np.zeros((len(ee),1))
ki=0
for i in range(0,len(eta)):
    eta=ee+ej
    if bei!=0:
        for j in range(0,bei):
            eta+=eb[ki]
            ki+=1
        
#The estimator
print(eta)


Best regards,
Vassilis