Recover a flux in RT2 space

mathematics
1

Hi, I have a finite element solution u in FunctionSpace(mesh,‘CG’,2) which could be simply be solved as follows. (The code is in Fenics 1.5.0 version.)

"""
FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.
Test problem is chosen to give an exact solution at all nodes of the mesh.
  -Laplace(u) = f    in the unit square
            u = u_D  on the boundary
  u_D = 1 + x^2 + 2y^2
    f = -6
"""

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'CG', 2)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Now, I would like to construct a flux, sigma, in RT2 space that satisfies the following conditions (partially in sudo code)

RT2 = FunctionSpace(mesh,'RT',2)
sigma = Function(sigma)
DG0 = VectorFunctionSpace(mesh,'DG',0)
P0 = TestFunction(DG0)
n = FacetNormal(mesh)
inner(sigma, P0)*dx = inner(grad(u_h), P0)*dx
inner(sigma, n) = inner(avg(grad(u_h)),n) on every interiori facet F (sudo code here.)

I could not set up a global variational formulation for this in Fenics. Therefore, I tried to manually set up the DOFs. However, without success yet. My code is too lengthy thus I simply post the question here.
My sincere thanks in advance for your time!