In order to evaluate changes in the geomagnetic field due to various tectonic models, a new method using the surface integrals over magnetized two dimensional bodies is presented. This method is based on combination of the boundary element method (BEM) with the surface integral formula. BEM allows calculations for arbitrarily shaped magnetoelastic bodies. The technique using surface integrals saves CPU time compared with a code using volume integral formula. In particular, BEM method can deal with complicated tectonomagnetic problems including actual surface topography. In order to test the new method, we make comparative calculations for what is called Yukutake model whose solution has been well investigated. The result for Yukutake model obtained using the present method is essentially equivalent to those obtained using numerical volume integrals. The method is also applied to tectonomagnetic problems in which topographic effect is taken into account. Initial calculations of topographic effect indicate model corrections are about 1nT or so.