This paper presents an efficient numerical technique for stress analysis of three-dimensional infinite media containing cracks and localized complex regions. To enhance the computational efficiency of the boundary element methods generally found inefficient to treat nonlinearities and non-homogeneous data present within a domain and the finite element method (FEM) potentially demanding substantial computational cost in the modeling of an unbounded medium containing cracks, a coupling procedure exploiting positive features of both the FEM and a symmetric Galerkin boundary element method (SGBEM) is proposed. The former is utilized to model a finite, small part of the domain containing a complex region whereas the latter is employed to treat the remaining unbounded part possibly containing cracks. Use of boundary integral equations to form the key governing equation for the unbounded region offers essential benefits including the reduction of the spatial dimension and the corresponding discretization effort without the domain truncation. In addition, all involved boundary integral equations contain only weakly singular kernels thus allowing continuous interpolation functions to be utilized in the approximation and also easing the numerical integration. Nonlinearities and other complex behaviors within the localized regions are efficiently modeled by utilizing vast features of the FEM. A selected set of results is then reported to demonstrate the accuracy and capability of the technique.