Inclusions are common in bearing materials and are a primary site for subsurface fatigue crack initiation in rolling element bearings. This paper presents a new approach for computing the pressure, film thickness, and subsurface stresses in an elastohydrodynamic lubrication (EHL) contact when inclusions are present in the elastic half-space. The approach is based on using the discrete element method to determine the surface elastic deformation in the EHL film thickness equation. The model is validated through comparison with the smooth EHL line contact results generated using linear elasticity. Studies are then carried out to investigate the effects of size, location, orientation, and elastic properties of inclusions on the EHL pressure and film thickness profiles. Both inclusions that are stiffer than and/or softer than the base material are seen to have effects on the pressure distribution within the lubricant film and to give rise to stress concentrations. For inclusions that are stiffer than the base material (hard inclusions), the pressure distribution within the lubricant film behaves as though there is a bump on the surface, whereas for inclusions that are less stiff than the base material (soft inclusions), the pressure distribution behaves in a manner similar to that of a dented surface. Inclusions close to the surface cause significant changes in the contact stresses that are very significant considering the stress life relationship. For inclusions that are located deep within the surface, there is little change in the EHL pressure and film thickness.