Thanks for your help.

We can't directly integrate like this. You can do a whole-field operation to get the fluxes into a DG0 trace space:

Vt = FunctionSpace(mesh, "HDiv Trace", 0)

fluxes = Function(Vt)

# you have to code this trace projection by hand:

v = TestFunction(Vt)

F = (fluxes*v)('+')*dS + fluxes*v*ds - avg(inner(grad(sol), n))*v('+')*dS - inner(grad(sol), n))*v*ds

The inner(grad(sol),n)*v in last two terms give me an error:

UFL:ERROR Shapes do not match: <Grad id=4986788992> and <Zero id=4998074176>.

But the projection was okay when changed them to sol*v. Do you know where is this error coming from?

(FYI: my sol Function is defined on DG1 functionspace)

solve(F == 0, fluxes)

# Now fluxes contains your integrals.

Lawrence

Also, how this globally defined fluxes are connected to cells? In other words for a triangular mesh, is it possible to know which 3 components of fluxes.vector().array() are related to cell E?