The paper presents an isogeometric shape optimization method that is based on Bézier triangles. Bézier triangles are used to represent both the geometry and physical fields. For a given physical domain defined by B-spline boundary, triangular Bézier parameterization can be automatically generated. This shape optimization method is thus applicable to structures of complex topology. Due to the use of B-spline parameterization of the boundary, the optimized shape can be compactly represented with a relatively small number of optimization variables. In order to ensure mesh validity during shape optimization, we adopt a bi-level mesh representation, where the coarse mesh is used to maintain mesh quality through positivity of Jacobian ordinates of the Bézier triangles. The fine mesh is used in isogeometric analysis for numerical accuracy. Numerical examples are presented to demonstrate the efficacy of the proposed method.