An effective boundary element method (BEM) for three-dimensional elastostatic contact problems is developed by using the penalty method. Contact conditions of the slip, stick and opening states on the contact surface of the two bodies are assembled into the total virtual work functional and a system integral equation of boundary values is formulated without any increase of the unknown parameters by introducing a penalty coefficient. A BEM based on this formulation is applied to typical three-dimensional contact problems such as the contacts between a cylinder and a block, and between cylinders. The shapes and the displacements are also measured by the experiment of the silicon rubber. The numerical results agree well with the analytical solutions and the experimental results.