There are no questions here, so I'll just give some comments.
A matrix corresponding to the gradient can be computed as
u = TrialFunction(...)
v = TestFunction(...)
A = assemble(inner(grad(u), v)*dx
Here v is a field of Nedelec type and u needs to be in
a compatible space.
Notice that the matrix (as always in the FEM method) will be
a mapping between the nodal (left hand side) representation
and the dual (right hand side) representation. If you want
the gradient in the nodal representation you should invert
the (lumped) mass matrix or something similar.