For some time now Paraview has supported arbitrary order Lagrange finite elements, see
Currently the iovtk module only supports P0 and P1 elements. On the other hand deal.II and mfem already support this:
It would be nice to see this on the FreeFem roadmap too.
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);
