When Hertzian theory is not available and it can be considered contact area is existing on the plane, numerical calculation is performed with superposing Boussinesq or Cerruti solution. In that calculation, generally, uniform pressure is assumed in element. On the other hand, MLMI method has spread, and contact pressure can be able to calculate with high accuracy because of using high-order interpolation of influence coefficient. But concerning sub-surface stress calculation, there are little announcement of solutions. Author developed both quadrangle linear pressure element and subsurface shear stress calculation method which uses linear pressure element. In this paper, author explains the derivation process of those. In addition, one case are taken up, and sub-surface shear stress by linear pressure element is compared with that by conventional uniform pressure element. As a result, in the case of uniform pressure element, the waving phenomenon of shear stress was seen near the contact surface. On the other hand, in the case of linear pressure element, correct solution was obtained in whole area.