One question:
when I get the fluid forces px,py,pz,vx,vy,vz to compute actual lift I’ve done vy-py, but to compute the lift coefficient for a 3D problem, namely Cl = lift/(0.5 rho v^2 * Area), is there a simple way to find the area of the surface (wing) via pyFR or python scripts where the integration of this forces has been done such that one can normalize correctly the lift force?

This is not something which PyFR can easily calculate and is often up for debate. As an example, for a wing, there is the choice between the wetted area and the planform area. Both are widely used.

Regarding this question, I want to know if my assumption is correct:

to compute lift coefficient with pyFR I have to do:

(vy + py) / ( 0.5 rho v*v * Area)

(previously I said I was doing vy - py)
I have to sum both viscous stresses and pressure force because in the solver you already compute the scalar product between the surface outward normal to airfoil and the unit versor which defines the lift direction (y in this case).
Is this true?
Because aerodynamics formula generally refers to

but the dot product has to be computed to account for both pressure and suction side