Hi Tuomas,
Are you sure you want to use RT1xP1, as opposed to RT1xP0[DG]? This would be the usual element for 'mimetic' discretisations.
The code for RT1xP0 would be [possibly with typos...]
---------------------------
from firedrake import *
mesh2d = UnitSquareMesh(10, 10)
layers = 6
mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers)
Uh_elt = FiniteElement('RT', triangle, 1)
Uv_elt = FiniteElement('DG', interval, 0)
U_elt = HDiv(OuterProductElement(Uh_elt, Uv_elt))
U = FunctionSpace(mesh, U_elt)
P1 = FunctionSpace(mesh, 'CG', 1) # automatically gives the appropriate P1 x P1 element
f = Function(P1)
f.interpolate(Expression('-x[0]'))
sol = Function(U)
test = TestFunction(U)
# sol + grad(f) = 0
F = inner(sol, test)*dx + dot(grad(f), test)*dx
solve(F == 0, sol)
print sol.dat.data.min(), sol.dat.data.max()
---------------------------
FIAT will throw an error if you try to do this with RT1xP1, however (though it will let you build a curl-conforming 'edge' element by using HCurl). If you genuinely need RT1xP1 then I can probably turn off some FIAT validation so that it can be used. I think it's more likely that you want RT1xP0, however, and a corresponding vertical velocity element of P0xP1, HDiv'd.