Indeed this would be nice, but on the other hand it is trivial to hack in FreeFEM itself and then save the solution. Whenever better graphics are needed (or when using more exotic elements than Lagrangian), I do an h-refinement and the interpolation in FreeFEM itself, with trunk and u=u, like in the snippet below. Maybe not a real solution to your feature request, but it actually works, and you can make a macro out of this to get a certain degree of automation.
load "iovtk"
string tag = "u";
int[int] ord = [1];
border a001(t = 0.0, 2.0*pi){ x = cos(t); y = sin(t);}
mesh T = buildmesh(a001(10));
fespace S(T, P2);
S u, v;
solve Laplace(u, v) = int2d(T)( dx(u)*dx(v) + dy(u)*dy(v) ) + on(a001, u=x*y);
plot(u, fill=true, nbiso=64, wait=true);
savevtk("u-wo.vtu", T, u, dataname=tag, order=ord);
T = trunc(T, true, split=3);
u = u;
plot(u, fill=true, nbiso=64, wait=true);
savevtk("u-w.vtu", T, u, dataname=tag, order=ord);
Actually I use the same approach when I need to see how vector fields look.
mesh T = square(1, 1);
fespace S(T, P0);
fespace V(T, RT0Ortho);
V [ux, uy] = [0, 0];
ux[][2] = 1;
T = trunc(T, true, split = 50);
[ux, uy] = [ux, uy];
S unorm = sqrt(ux*ux + uy*uy);
plot(unorm, [ux, uy], fill=true, nbiso=128, wait=true);
