We have succeeded in developing a numerical code to obtain structures of non-axisymmetric polytropes the shapes of which are sustained by internal flows. In order to compute stationary non-axisymmetric polytropes, two important points must be properly treated: (i) how to specify the boundary condition for the flow velocity, and (ii) how to handle the boundary condition on the stellar surface which is not known beforehand. As for the boundary condition for the flow velocity, we specify the form of the distribution of the difference between the theta component of the flow velocity and its mean value at one meridional cross-section. Concerning the surface condition, we introduce a `surface-fitted' coordinate system in numerical computations. By applying our new numerical code we have obtained many stationary non-axisymmetric equilibrium models the shapes of which are non-rotating in the inertial frame. Two main results of our computations are as follows: (i) the theta component of the flow velocity must depend on the height from the equator for polytropes to be deformed to some extent and (ii) the flow is compressible. Unless we include this z-dependence for the flow, it seems that stationary states can exist only for nearly incompressible polytropes and only for very slightly deformed non-axisymmetric configurations as far as `bar'-type non-axisymmetric configurations are concerned.