The present paper suggests a modification of boundary elements numerical method, which uses displacement discontinuity method in 3D space. The current method can effectively decrease the computational elements number in 3D space because final elements are placed only on the surface of a fracture simulating the discontinuity of elastic medium. Thus the dimension of the problem is decreased. The new fundamental solutions were developed based on double layer potentials. The method proved to be effective for hydraulic fractures modeling and their interaction with natural faults in host rock formations.