Direct simulations of the dynamics of a large number of deformable droplets are necessary for a more accurate prediction of rheological properties and the microstructure of liquid-liquid systems that arise in a wide range of industrial applications, such as enhanced oil recovery, advanced material processing, and biotechnology. The present study is dedicated to the development of efficient computational methods and tools for understanding the behavior of three dimensional emulsion flows in Stokes regime. The numerical approach is based on the accelerated boundary element method (BEM) both via an efficient scalable algorithm, the fast multipole method (FMM), and via the utilization of advanced hardware, particularly, heterogeneous computing architecture (multicore CPUs and graphics processors). Example computations are conducted for 3D dynamics of systems of tens of thousands of deformable drops and several droplets with very high discretization of the interface in shear flow. The results of simulations and details of the method and accuracy/performance of the algorithm are discussed. The developed approach can be used for the solution of a wide range of problems related to emulsion flows in micro- and nanoscales.