Dear Lawrence, 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? Thank you, Sarraf