In this paper, the geometrically nonlinear behaviors of functionally graded composite shells under large displacements and large rotations are analyzed. Instead of constructing conventional layer-wised models, equivalent single-layer composite shell elements based on general displacement fields are developed to model the inhomogeneity of composite materials. The mixed interpolation of tensorial components technique is used to eliminate membrane locking and shear locking phenomena. The assumed covariant strain fields are initially constructed before being used to evaluate the second Piola–Kirchhoff stress and construct stiffness matrices in local Cartesian coordinates. The standard full Newton–Raphson method is utilized to solve a system of nonlinear equations. Geometrically nonlinear analyses are performed on different thin-walled composite structures including thin cantilever beams, slit annular plates, pinched cylindrical shells, and hemispherical shells. Obtained results demonstrate good convergence characteristics and modeling capability of the developed quadrilateral composite shell elements in analyzing thin-walled composite structures.