A three-dimensional boundary element method (BEM) implementation of the interaction integral methodology for the numerical analysis of mixed-mode three-dimensional cracks is presented in this paper. The interaction integral is evaluated from a domain representation naturally compatible with the BEM, since stresses, strains and derivatives of displacements at internal points can be evaluated using their appropriate boundary integral equations. Special emphasis is put in the selection of the auxiliary function that represents the virtual crack advance in the domain integral. This is found to be a key feature to obtain reliable results at the intersection of crack fronts with free surfaces. Several examples are analysed to demonstrate the efficiency and accuracy of the implementation. © 2004 Elsevier B.V. All rights reserved.