Understanding the mechanisms of rock fracture is of key importance to the mining and petroleum industries. Rock masses feature complicated geometry and structure including joints, heterogeneities of different scales (e.g. grains, pores, macroscopic inhomogeneities, etc.) and may be subject to various effects of injected fluid pressure, temperature gradient, etc. Therefore, comprehensive three-dimensional computational models that would allow to adequately treat complex behavior of a rock mass are required. Among the industrial applications of such models are the quantification of safety of underground workings and simulation of hydraulic fracturing. The dissertation presents a new Boundary Element Method-based technique for analysis of a three-dimensional elastic medium containing multiple cracks and/or openings of arbitrary shapes. The technique employs planar triangular elements to discretize the boundaries and quadratic polynomials to approximate the boundary unknowns, with two options of the arrangement of the nodal points on the elements. The novel features of the technique include the following: • the use of complex variable formalism involving various combinations of the geometrical parameters and the elastic fields, e.g. components of tractions and displacement discontinuities in the plane of the considered element; • analytical integration with the use of Cauchy-Pompeiu formula to reduce the surface integrals to the contour ones; • “limit after integration" approach, i.e. enforcing the boundary conditions after the discretization and analytical handling of the internal fields, by allowing the field point to reach the boundaries. The method can still capture the behavior of stress field near the crack fronts (tips) although no special approximating functions (tip asypmtotics) are used. The solutions of some benchmark problems are provided to demonstrate the capabilities of the proposed method.