We present a novel framework for topological shape optimization of curved thick-shells subjected to external loads. Our method integrates the Level Set Method (LSM) with a diffuse interface, a Hadamard shape derivative, and multi-patch isogeometric analysis (IGA) into a gradient descent algorithm to systematically capture the evolution of the shape. This integration enables us to directly manipulate CAD-compatible geometries and analysis techniques and to obtain the results as a CAD surface. The novelty lies in the utilization of multi-patch IGA models based on NURBS functions, which allows us to maximize the stiffness or minimize the von Mises Lp norm while simultaneously minimizing the volume of the shell by searching for an optimal material distribution within its mid-surface. The material was modeled under a small strain assumption in linear elasticity using a Reissner-Mindlin kinematics shell model in plane stress, as in [1]. The effectiveness of the present approach is demonstrated on several curved multi-patch geometries in 3D. Let us recall that IGA defines the physical domain by means of parameterized NURBS [2]. This natural parametrization allows solving a series of optimization problems using the direct differentiation of the problem. However, the recursively deformed shape during the optimization can conduct to ill-defined geometries during the iterative optimization process. To overcome this difficulty, the LSM defined over a fixed domain and the Hadamard derivative have been coupled with IGA to track the evolution of the shape in the descent direction of the cost function. Topological shape design for shells has already been addressed in the literature. For example, the compliance and maximum stress minimization with a density-based approach and a finite element method has been addressed for 2D plates in [3], minimizing the maximum stress, while considering the curvature of three-dimensional shells was developed in [4]. Optimal microstructure for accounting simultaneously for in-plane stiffness, out-of-plane stiffness, and the extension–bending coupling effects in panels obtained by inverse homogenization, the Hadamard shape derivative and a LSM in [5], etc. Furthermore, while previous research [6] showcased the effectiveness of topology optimization coupled with IGA, our approach specifically focuses on bridging the existing gap between topology optimization and the manipulation of CAD geometries for 3D multi-patch structures within a shell