It is hard to give a direct answer to this, I talked about it a lot in this post: Mesh resolution, order, and curved boundaries - #3 by WillT
Depending on the settings you are using in Pointwise, the mesh will only be curved near smooth boundaries. As it is doing a polynomial fit, the convergence of the polynomial modes will likely be spectral. I.e. if the polynomial approximating the geometry in an element has the form:
where c are the modes, \phi is some polynomial basis, and (\zeta,\eta) are the coordinates in the reference domain. Then we know that as the degree of the individual bases increases then the scale of c will decrease quite rapidly (if the shape of the geometry is in some H^p space). This is really just a restatement of Taylor’s theorem.
What I am trying to get across is that yes the geometry might be better approximated by a higher-order polynomial, but how much difference is that really going to make? This is also probably why Pointwise is struggling to generate the mesh, the high order coefficients in the geometry approximation are probably small and difficult to converge.
On the flux points, Pointwise will be placing its own set of points on the boundary. PyFR reads the location of those, builds a fit of the surface and then places its own points. Therefore, I don’t think there is a guarantee that the flux points will be on the physical surface, but rather that they will be on the projected surface provided by Pointwise. Abe etal. (2018) is a good starting point.
The point I’m trying to make is to be pragmatic, using a lower order geometry representation is prefered as it makes things easier on the meshing front, and reduces aliasing based problems. But if you find that your solution is faceted then maybe you need to up the geometry order. Others will have views on this as well, so feel free to disagree.