Many calculating methods have been already developed for solving contact problems of parts such as gears, cams, and followers under fluid film lubrication conditions considering the temperature and pressure dependence. Similarly, the determination of the elasto-hydrodynamic pressure distribution the processes taking place in the lubricant and the contacting bodies, as well as in their environment, have to be dealt with simultaneously for the determination of the temperature field. A system of equation for the modelling of thermo-elastohydrodynamic lubrication between two contacting bodies containing hydrodynamic, thermodynamic, and strength problems is a highly non-linear system which becomes even more so if the temperature and pressure dependence of the material properties are considered. To solve this system, scientists started to use the finite element formulation in the 1960s and it was found to be a promising and reliable method. Earlier, the lubrication analysts used only the h-version finite element method (h-FEM) till 1991, when the first usage of the p-version finite element method (p-FEM) was published in the literature. In order to reduce the problem, in case of point or line contact, the contact bodies can be handled as semi-infinite ones. Following this simplification that had been successfully applied for the gap size determination, a substructure model was defined using analytical solution of the moving heat source. Instead of an iterative way between the solid and fluid problem, in this paper we present an efficient solution when thermal model for lubricant and surfaces were coupled and solved by a direct numerical method in one step.