A six-node co-rotational curved triangular shell finite element with a novel rotation treatment for folded and multi-shell structures is presented. Different from other co-rotational triangular element formulations, rotations are not represented by axial (pseudo) vectors, but by components of polar (proper) vectors, of which additivity and commutativity lead to symmetry of the tangent stiffness matrices in both local and global coordinate systems. In the co-rotational local coordinate system, the two smallest components of the shell director are defined as the nodal rotational variables. Similarly, the two smallest components of each director in the global coordinate system are adopted as the global rotational variables for nodes located either on smooth shells or away from non-smooth shell intersections. At intersections of folded and multi-shells, global rotational variables are defined as three selected components of an orthogonal triad initially oriented along the global coordinate system axes. As such, the vectorial rotational variables enable simple additive update of all nodal variables in an incremental-iterative procedure, resulting in significant enhancement in computational efficiency for large deformation analysis. To alleviate membrane and shear locking phenomena, an assumed strain method is employed in obtaining the element tangent stiffness matrices and the internal force vector. The effectiveness of the presented co-rotational triangular shell element formulation is verified by analyzing several benchmark problems of smooth, folded and multi-shell structures undergoing large displacements and large rotations.