Contact phenomena between deformable bodies are a common problem in engineering. The surface stress distribution and subsurface stresses depend on many parameters. In view of increasingly strict tolerances the effects of local variations in material properties need to be accurately predicted. In this paper a multigrid solution method is presented for contact problems between three dimensional elastic heterogeneous materials. The contact problem is incorporated as boundary condition in the multigrid solution of the displacement equations for the volume. First, validation results are presented. Subsequently a study is presented for soft and hard clusters of inclusions. Finally, results are presented for a contact problem involving a realistic case of a polycrystalline material representative for applications with ceramic materials.