I could solve the problem. It seems that project does not work well as indicated here so I used interpolate instead. In addition I defined the identity as an expression-although I am not sure if that was also affecting the code-as:
Identity = Expression(("x[0]","x[1]"), degree=1)
Idt = interpolate(Identity, V)
Finally, here Rogens indicates that the approximation only works with linear elements, so I changed the degree of the vector space to 1 and it worked. I am going to run some extra simulations to see how it goes with more complex geometries.
I also want to add that the fact that the implementation only works with linear elements is a big issue in this particular case. Heine showed that the order of convergence of this isoparametric approach is degree-1 (degree of Lagrange elements). So in the case of linear elements, it is not convergent.
If anyone finds a way to use higher-order elements I would really appreciate if they share it with me. I think this issue can be overcome if we could import a quadratic mesh including the mid nodes of each element. I guess it might be already possible but I have not found how to do it.