Isogeometric analysis (IGA) provides a feasible technique to seamlessly integrate computer-aided design (CAD) into the existing finite element analysis. In this article, a Bézier elements-based IGA method is established to address the stress-related topology optimization of structures, which incorporates geometrical representation, structural analysis and topology optimization into a unified process as well as can naturally handle the curvilinear edge elements for the stress problems. To enhance the smoothness of the structural boundaries, a multilevel isogeometric density field (IDF) defined at control points (CPs) is designed by utilizing the values of the non-uniform rational B-splines (NURBS) basis functions from the previous design iteration. A global stress measure is approximated by a smooth maximum function with an additional error term. The k-means clustering originally from signal processing is used to partition the stress observations with dynamic centroids and indices, which improves the accuracy of the sensitivity results. The sensitivity analyses of the stress-related problems are derived from the multilevel IDF model. Three kinds of numerical results are provided to illustrate the effectiveness of the proposed method to solve stress-related problems via the multilevel IDF model for achieving the optimal solution with smooth boundaries.