Boundary element method (BEM) is a powerful computational tool for analysing piezoelectric problems. The BEM has been so well developed during the past 40 years that it has been considered as a very popular computational tool. This method consists of formulating the engineering problem in terms of an integral equations relating only boundary values and determining its solutions numerically. Thus, it requires only a surface discretization, rather than a full-domain discretization with “domain-type techniques, such as the finite element method. If interior domain values are needed, however, it can be subsequently calculated from boundary data. Remarkable advances in this area during last decades can be found, for example, in Manolis and Beskos[1] and Qin [2,3]. For piezoelectric problems, BEM, as an important complementary tool for analytical method [4-8], has been the subject of fruitful scientific attention by many a distinguished researcher e.g. Lee and Jiang [9], Denda and Mansukh [10], Sanz, Ariza and Dominguez [11], Lee and Jiang [9] derived the boundary integral equation of piezoelectric media by the method of weighted residuals for plane piezoelectricity. Lu and Mahrenholtz [12] presented a variational boundary integral equation for the same problem. Ding, Wang and Chen [13] developed a boundary integral formulation which is efficient for analysing crack problems in piezoelectric material. Xu and Rajapakse [14,15] presented the formulations for problems of piezoelectric solids with various defects (cavities, inclusions, cracks, etc.). Pan [16] derived a single domain BE formulation for 2D static crack problems. Denda and Lua [17] developed a BEM formulation using Stroh’s formalism to derive the fundamental solution without numerical results. Qin [18] presented a BEM formulation for cracked piezoelectric materials with half-plane boundary. Later, Qin and Lu [19] extended the model presented in [18] to the case of piezoelectric materials containing both crack and inclusion problems. Davi and Molazo [20] used the known sub domain method to formulate a multi-domain BEM, well suited for crack problems, by modelling crack faces as boundaries of the different sub domains. Groh and Kuma [21] developed a direct collocation boundary element code with a sub domain technique for analysing crack problems and calculating stress intensity factors. Qin an Mai [22] investigate crack-hole interaction of piezoelectric materials using BEM. Zhao et al. [23], presented a boundary integral-differential model for interfacial cracks in 3D piezoelectric solids. Qin [24,25] imbedded BEM into micromechanics model for evaluating effective properties of piezoelectric solids. Analysis of Piezoelectric Solids through Boundary Element Method