In this paper, a new type of cohesive element that employs multiple subdomain integration (MSDI) for improved cohesive stress integration accuracy of bonded plate/shell elements has been formulated. Within each subdomain, stress integration is compatible with existing schemes such as Gaussian integration (GI), Newton–Cotes integration, or the mixed Gaussian and subdomain integration (mixed GI&SDI). The numerical accuracy, efficiency, and robustness of this element when employing three integration methods for MSD cohesive stress integration have been evaluated and compared through a benchmark mode-I fracture problem of bonded double-cantilever plates. The MSDI offers at least 50% improvement of numerical accuracy as compared to the best integration method in literature and has the best numerical robustness. This significant improvement pushes the structural mesh size restriction from limiting size of 1/3–1/5 cohesive zone length to 1.5–2 times the cohesive zone length. The formulation is very easy to be implemented into any finite element programs including commercial packages. Furthermore, this formulation enables the use of dual-mesh for delamination analyses of bonded structural shells/plates, which is of practical importance because it greatly reduces the burden of mesh generation for complicated composite structures. It has also been demonstrated that using high-order shell/plate elements can improve the numerical accuracy in general because the nonlinear deformation profile permitted by this type of elements can better describe the nonlinear deformation in the crack-tip element (partially bonded elements).