This article focuses on the stochastic modeling of nonlinear systems featuring discontinuities in their response surface. More specifically, original developments are presented to improve the efficiency of multi-element polynomial chaos expansion for such systems. First, an automated detection procedure of the discontinuities is proposed. It relies on an iterative algorithm with a polynomial annihilation edge detection method and support vector machine classification algorithms leading to the representation of the discontinuity as a B-spline curve. Based on this curve, an ad-hoc decomposition of the variable space is considered and an original approach for the application of polynomial chaos expansion on each subdomain yields a robust and versatile way to compute the response surface of the system of interest. The proposed methodology is detailed and applied to several nonlinear academic systems such as a circular discontinuity and the Duffing oscillator including one or two discontinuities in its response surface. Through these applications, it is evidenced that, compared to the classical polynomial chaos and multi-element polynomial chaos expansions, the proposed methodology yields an approximation of the discontinuous responses that is both more accurate and less computationally expensive. The influence of the main parameters of the proposed methodology is also analyzed in details. This parametric analysis underlines the robustness of the methodology and highlights the key parameters in terms of computational cost and accuracy of the discontinuities. The proposed methodology is finally applied to an industrial application, it allows to efficiently compute the surface response of an industrial compressor blade undergoing structural contacts.