A shear deformable, finite element penalty formulation is developed for predicting stress fields, particularly interlaminar shear and normal stresses, in thick composite laminates and near a delaminated free edge in a laminated plate. A fully three-dimensional 20-noded finite element is coupled with the feature of integrating the equilibrium equations to compute interlaminar shear stresses and characterize three-dimensional stress fields. As an alternative to using the constitutive equations, the integration of the equilibrium equations provides an improved variation of interlaminar shear stresses through the thickness. Penalty functions are utilized in the formulation to represent continuity or discontinuity between layers of elements. Proficient use of these penalty functions allows simple definition of critical interfaces in a composite laminate where delamination has occurred. The flexibility of the formulation permits the consideration of different sizes of debonds without having to create a new model. In this paper, thick composite laminates are considered to demonstrate the improved variation of interlaminar shear stresses computed by integrating the equilibrium equations. Finally, a [0 deg/90 deg/0 deg] undelaminated plate subjected to uniform axial extension is compared with the same geometry containing a free-edge debond, one ply thickness in depth.