Thank you for the suggestion.
I have thought about it, but it does not solve my problem, because u and uh are on different meshes and the compiler complains.
The meshes are “equivalent” in the sense that one is movemesh of the other. In this case mapu and mapt work, but I prefer to build the matrix in a direct way.
It seems that the modification of the matrix M from varf is the way. But this requires knowing the ordering of dofs, ie to which dof corresponds the entry M(i,j), which in the case of vector FE spaces, such as [P1b,P1b,P1] seems not available in the docs - I have posted a message with this regard at Order of variables in matrix generated with varf with a [P1, P1, P0] fespace - #6 by pi.3.141.