Disclaimer: I was searching for 'higher order FE interpolation' when I came across this question.
Univariate Interpolation (i.e., 1-D)
Higher order interpolation could be interpreted in multiple ways. But one of the most common ways of doing higher order interpolation is piecewise polynomial interpolation where the polynomials are of degree 2 or higher (degree is same as order in this case). Degree 3 is most popular, and the method is known as cubic spline interpolation.
Note: As far as I can tell, piecewise polynomial interpolation might as well be called FE interpolation (in fact, a specific version of FE, due to choice of basis/trial and weight/test functions) and I fail to understand why I haven't come across a textbook that points that out.
Note: 'FE interpolation' is different from 'FE for differential equation solving'.
Multivariate Interpolation (i.e., 2-D and higher)
For 2-D, a degree 3 interpolation on structured rectilinear mesh would be called bicubic, for 3-D tricubic.
A structured rectilinear mesh is one in which the grid points are ordered in X and Y dimension, and in case of 3-D, Z dimension as well. Structured or regular. Mesh or grid. Basis or trial. Weight or test. (If you think proliferation of terminology is a pathological problem of this discipline, welcome to the club.)
MATLAB (as well as the python scientific stack) provide ways of doing such interpolation, e.g., in
scipy the API is
RectBivariateSpline for 2-D. However finite element codes deal with simplices, e.g., triangles in 2-D, and tetrahedrons in 3-D, and this poses a problem (the mesh is unstructured). You can use the API
griddata to "griddify" your solution, but my guess is that it would be less accurate than directly interpolating $u_h$, plus my personal experience is that
griddata based interpolation is very finicky in handling boundaries, and generally fails when a bit of extrapolation is involved.
How to interpolate directly from a solution on unstructured mesh? I don't have a very good answer. Based on what I'm learning these days, FE interpolation (or can it also be called unstructured multivariate spline interpolation?) should work but the book I have (Heath, Scientific Computing) doesn't have any mention of FE in the chapter on interpolation (but then it doesn't discuss multivariate interpolation at all). Langtangen's FE texts do start with FE interpolation (that's where I learned about it) before moving on to FE DE solving, but even his discussion is for pedagogical purposes (to help your understand later chapters), not as an application. Plus Langtangen makes no mention of whether or not hat/bubble FE interpolation is just spline interpolation in disguise.
Remark: In computer graphics, (e.g., in Foley, Computer Graphics) there is a similar concept of 'parametric cubic curves' (1-D), and 'parametric bicubic surfaces' (2-D). How does that compare or contrast with interpolation? I haven't come across a text that answers that.
Should we use FE approximation with higher order elements?
As far as I can tell, the authors you mentioned do not mean using higher-order FE for differential equation solving. And as I said, in my non-expert opinion, FE can be used for interpolation (regardless of what you used for your differential equation solving), but again, I can't find textbook level literature on FE interpolation, and I fail to understand why it's not routinely used for interpolation purposes (for multivariate solutions on unstructured meshes).
P.S.: The book you mentioned doesn't have the word "Solving" in its title.
Edit 1: This stackoverflow answer mentions radial basis functions as an option for unstructured multivariate interpolation. Firstly, I would disagree with the statement "splines are radial basis functions". Splines are based on polynomials, while RBF could be any function of $r$ that you want. Linear and cubic spline are definitely not radial in 1-D, since both $x$ and $x^3$ are odd. Secondly, last I checked, the RBF implementation in the python scientific stack was even slower than
griddata followed by
RectBivariateSpline (but I'm using an older version of numpy/scipy so it's possible they've sped it up in the recent versions, see this for more info).
In general, Irregular grid section in wikipedia page on multivariate interpolation is probably a good point of reference. 'Triangular irregular network' looks similar to what I'm talking about (but in 2-D only). However it only mentions two options, natural neighbor, and linear interpolation, which are 'probably', and 'definitely', of degree 1, respectively (nearest neighbor would be degree 0). Hence they are not higher order. For higher order, you would need to look for a triangular irregular network of degree 2 or higher. And in 3-D it would probably be called tetrahedral irregular network.