The stability, accuracy, and computational efficiency of higher order Galerkin (hp-type) finite elements for steady flow of viscoelastic fluids past square arrays of cylinders and through a corrugated tube with two different polynomial approximating spaces, namely truncated and product, have been investigated. It has been shown that both spaces produce a stable discretization with an exponential convergence rate toward the exact solution without an upper Weissenberg number limitation. Based on global deviation from mass and momentum conservation, it is shown that the truncated space provides a better quality solution at a reduced computational cost for a given number of degrees of freedom. In addition, the restrictions on the relative order of approximating polynomials for stresses and velocities have been examined. It has been shown that without splitting of the stress into purely viscous and elastic components the most cost efficient solution is obtained when the stresses and velocities are approximated by the same order polynomial, while with this splitting the best computational performance is obtained when the stresses are discretized with a polynomial order of one less than the velocities.