As Tormod says there is a fast version in fenicstools.
However, I would also like to point out that your above code is actually
correct. But you should be aware that you will obtain the gradient in
the "right-hand side representation". What you probably wanted was the
representation used on the left-hand side. To get to that you would have
to multiply your gradient with the inverse of a mass matrix. That is,
you need to solve as follows.
(u, v are vectors spaces )
m = uvdx
L = inner(grad(f), v)*dx
solve (a == L)
This will often be the case for finite element methods or in general
methods that employ weak formulations. That is, the scaling of vectors
on the left- and right-hand side vectors are different.