A semi-analytical method is developed from existing elements - innovative analytical and numerical solutions - in order to deal with problems of heterogeneous elasto-plastic contact. The Eshelby's equivalent inclusion method allows to describe the effect of heterogeneities (voids, inclusions, fibers or strands). One of the contact bodies contains several heterogeneities with rectangular or ellipsoidal forms, and related degenerated forms (oblate and prolate ellipsoid, sphere, cylinder, flat disk, ... ). This method is modified and improved to take into account the mutual influence between neighboring heterogeneous inclusions, as well as the effect of debonding interface between the heterogeneity and matrix. Ellipsoidal fibers of woven composite could be represented through the use of unit cuboidal heterogeneities. Voxelization algorithms consist in discretizing a solid geometry (three-dimensional) in many prismatic elements. This method is equivalent to the two-dimensional pixelization. Similarly, the fibers could be also discretized by many ellipsoids of representative size with equivalent properties. Through the use of elementary analytical solutions and numerical methods such as Fast Fourier Transforms (FFT-3D and 2D-FFT), only a few minutes are necessary to solve the contact problems on a grid of several thousands points. Fretting simulations (in gross slip or partial slip) are finally performed, taking into account the effects of heterogeneities and/or coating.