a posteriori error estimator - face residuals?

Question

Hi - I'm trying to compute a posteriori error estimators for Poisson demo of the type

$$ \mathcal{E}_T^2 = h_T^2\|-\Delta u - f\|_{L^2(T)} + \sum_{e\in \partial T} h_e\| [[n\nabla u]]\|_{L^2(e)}$$

where $h_T$ is a representative size for an element $T$, and $h_e$ is a representative size for a given element edge/face.

The volume residual is easy to compute using

h = CellSize(mesh) 
DG0 = FunctionSpace(mesh, "DG", 0) # element indicator function
w = TestFunction(DG0)
residual = h**2*w*(-div(grad(u))-f)**2*dx 

and returns an array indexed by cell index, which can be accessed via

for c in cells(mesh):
    cell_energy[c.index()] # pseudocode, just shows accessing

Is there a similar way to compute the sum of jumps over each edge for a given element, and store them in an element-centered fashion?

Thanks!

Answer 1

Something like the following should work (untested):

h = CellSize(mesh)
n = FacetNormal(mesh)
DG0 = FunctionSpace(mesh, "DG", 0)
w = TestFunction(DG0)
residual = h**2*w*(div(grad(u))-f)**2*dx + avg(w)*avg(h)*jump(grad(u),n)**2*dS
cell_residual = Function(DG0)
assemble(residual, tensor=cell_residual.vector())

Comments

Works like magic! Thanks Christian.