This paper extends our recent work  to evaluate the elastic fields and effective modulus of a composite containing arbitrarily shaped inhomogeneities for both two-dimensional (2D) and three-dimensional (3D) problems. Based on Eshelby’s equivalent inclusion method (EIM), the material mismatch between inhomogeneities and matrix phases is represented with a continuously distributed eigenstrain on inclusions. Since there exists singularities at the vertices of polyhedral inhomogeneities, domain discretization of inhomogeneities is used to interpolate the eigenstrain distribution, and high accuracy is obtained by using the closed-form integrals of the source field. The influence of the singularity decays in and for 2D and 3D problems respectively. Because Eshelby’s tensors depend on the shape instead of the size, the iBEM is particularly suitable for cross scale modeling of composites with a wide range of particle sizes. Although a number of elements are required to provide high fidelity results of the local field around the vertices, in general, very few elements or a single element are enough for each inhomogeneity to obtain the convergent solutions of the effective material properties, which enables virtual experiments of a composite containing many inhomogeneities. This novel numerical method has been verified with the finite element method (FEM) with much more elements. Virtual experiments of composites with many particles demonstrate its versatile capability and great potentials.