We consider the one dimensional boundary driven harmonic model and its continuous version, both introduced in \cite{FGK}. By combining duality and integrability the authors of \cite{FG} obtained the invariant measures in a combinatorial representation. Here we give an integral representation of the invariant measures which turns out to be a convex combination of inhomogeneous product of geometric distributions for the discrete model and a convex combination of inhomogeneous product of exponential distributions for the continuous one. The mean values of the geometric and of the exponential variables are distributed according to the order statistics of i.i.d. uniform random variables on a suitable interval fixed by the boundary sources. The result is obtained solving exactly the stationary condition written in terms of the joint generating function. The method has an interest in itself and can be generalized to study other models. We briefly discuss some applications.