Precision of the Interpolation and Gradient operators of Flux Reconstruction methods

Dear PyFR developers

I have a question regarding regarding the precision of the interpolation and gradient operators of the FR/SD methods in hexahedral meshes.

Let’s say that I initialize a variable in the solution points using a general polynomial expression f:

f(x, y, z) = (x + y + z)^p

Where p is the order of the polynomial and x, y, z is the position in the physical space of a given solution point. Then, I use the interpolation operators (Legendre/Lagrange interpolation basis based on the elements’ solution points) to compute the numerical interpolated values of f in the flux/Gauss points f_g. I have observed that, if the flux reconstruction scheme is of order m then the interpolation operators will interpolate exactly a function f of order p = m - 1 in any type of hexahedral mesh, am I right? Moreover, regarding the approximation of the gradients \nabla f at the flux/gauss points \nabla f_g I have observed that this extrapolated gradient will be numerically exact for p = m - 1 in cartesian hexahedral meshes. However, for non-cartesian hexahedral meshes I observed that I can only compute numerically exact the gradient of a function f of order p = m - 2. This implies 0th truncation errors of the gradient computation for the FR2/SD2 method in general non-cartesian hexahedral meshes…

Do you know of any references which discuss about the topics that I have described beforehand?

Thank you very much for your help and insight,

Gonzalo