Interpolation of RT1 Function on manifold. Does it work?

Question

I am trying to interpolate a vector Expression on a manifold surface.
Maybe I am doing something stupid, but I thought this should work.
What does happen, looks like some kind of indexing error, as the values are OK, but probably in the wrong places.

M = UnitSquareMesh(10,10)

# Convert 2D mesh to 2D mesh in 3D space
mesh = Mesh()
e = MeshEditor()
e.open(mesh, 2, 3)
e.init_vertices(M.num_vertices())
for (i, x) in enumerate(M.coordinates()):
    e.add_vertex(i, x[0], x[1], 0.0)
e.init_cells(M.num_cells())
for (i, cell) in enumerate(M.cells()):
    e.add_cell(i, cell[0], cell[1], cell[2])
e.close()

mesh.init_cell_orientations(Constant((0, 0, 1)))

E = Expression(("x[0]*x[1]","0.0","0.0"))
Q = FunctionSpace(mesh, "RT", 1)
F = Function(Q)

F.interpolate(E)

output = File("out.pvd")
output << F

Comments

Looking at the output again, I am seeing a z-component, which shouldn't be there.

What I was trying to create was a 2D mesh in the X-Y plane, with a vector field confined to that plane. I realise some of the problem may be due to the interpolation of the RT elements with compute_vertex_values(), but really, there should not be a z-component...

I think this is a bug, so have raised on bitbucket as Issue #75