This paper presents the development of a highly efficient user-defined finite element for modelling the bolt-load distribution in large-scale composite structures. The method is a combined analytical/numerical approach and is capable of representing the full non-linear load-displacement behaviour of bolted composite joints both up to, and including, joint failure. In the elastic range, the method is generic and is a numerical extension of a closedform method capable of modelling the load distribution in single-column joints. A semiempirical approach is used to model failure initiation and energy absorption in the joint and this has been successfully applied in models of single-bolt, single-lap joints. In terms of large-scale applications, the method is validated against an experimental study of complex load distributions in multi-row, multi-column joints. The method is robust, accurate and highly efficient, thus demonstrating its potential as a time/cost saving design tool for the aerospace industry and indeed other industries utilising bolted composite structures.