This paper presented a multi-domain boundary face method (MD-BFM) to solve steady-state heat conduction problems in large-scale engineering structures. Based on BFM, and using a similar approach as that in single-domain problems, this paper derived the matrix assembly process in multi-domain problems. Then, temperature and other unknown quantities can be drawn from the discrete boundary integral equation. The multi-domain BFM was developed and applied to the steady state heat conduction analysis of an actual gravity dam, which contained 62 concrete layers (totally 125 domains). In order to verify the result, the temperature contours of the dam were compared with the result obtained from FEM (using Abaqus). The comparison chart has shown that the maximum temperature value of both results is 20.7 ℃. Moreover, the contours are highly consistent. Numerical results have demonstrated that our method can achieve comparable accuracy than other methods (e.g. the FEM) at much lower cost in terms of both computer resources and human labor.