In this work, the boundary face method (BFM) is applied to implement steady-state heat conduction analysis of solids containing a large number of open-ended tubular shaped holes of small diameters. A new meshing scheme is used to discretize the boundary integral equations (BIE) such that the holes can be modeled by a small number of surface elements while keeps the exact geometry, resulting in substantial savings in both modeling effort and computational cost. In the scheme, each tubular pipe surface is represented with a number of curvilinear tube elements similar to the 'hole element' proposed by P.K. Banerjee. To model the end faces that are intersected by the tubular holes, a special triangular element with negative parts is proposed. These elements are defined in the parametric space of the surface, and the exact geometry data can be directly available from CAD models of the solids. Numerical examples show that current implementation is very efficient in modeling of solids with many holes of arbitrary shape. The temperature and flux on the pipe surfaces or inside solids are obtained with high accuracy, even the local thermal concentration on and near the holes can be captured. (C) 2012 Elsevier Ltd. All rights reserved.