Which class member(s) of a PyFR integrator can I use to access the coordinates of the solution points? (i.e. the coordinates corresponding to the solution values accessible via the ‘soln’ property of a PyFR integrator)

I’ve been digging around the integrator object created when running the 2D incompressible cylinder example, and it seems that perhaps this is a candidate (assume we have an integrator called integ):
integ.pseudointegrator.system.ele_ploc_upts[0]

For the 2D cylinder example, this is a 3D array with the shape (16,2,196) which could imply that integ.pseudointegrator.system.ele_ploc_upts[0][i,j,k] is the jth coordinate of the ith solution point in the kth element, since integ.soln[0] is also a (16,3,196) tensor which is apparently structured such that the first axis is the solution point id, second axis is the variable (u,v,p) and the third axis is the element id.

If integ.pseudointegrator.system.ele_ploc_upts[0] does indeed contain the solution point coordinates, is this also applicable to other solvers in PyFR? If it does not contain the coordinates, is there some other way to obtain these?

My primary motivation for obtaining this information is interpolating the solution to arbitrary coordinates within the domain using Scipy.

The variable system.ele_ploc_upts is a list of NumPy arrays containing the solution point coordiantes for each element type in the simulation. The element type associated with system.ele_ploc_upts[0] is system.ele_types[0].

Interpolating the solution to arbitrary coordinates, however, is something which must be done with care. Ideally, for each point one wishes to interpolate to, you would determine what element that point is associated with (non-trivial for domains with curved elements) and then use iteration to determine a point in reference space at which to evaluate the polynomial basis at in order to obtain the desired solution in physical space (iteration being required as for curved elements it is not possible to simply invert the mapping).

Thanks for the detailed reply. I certainly hadn’t considered the usage of basis polynomials, I’ll look into that direction now. Especially considering the caveats you mentioned in your post, is there any functionality in PyFR to evaluate the basis polynomials directly from within Python?

There is functionality to evaluate the basis polynomials in reference space. Examples can be found in several of the plugins with the underlying code being in polys.py (see nodal_basis_at which will generate an interpolation operator which can then be applied to the solution points of an element). However, we do not currently have code to locate which element a particular point in physical space corresponds to, or given such an element decide where to evaluate it in reference space.