This presentation introduces a numerical method for calculating the energy release rates and their higher order derivatives for a multiply cracked body under general mixed-node conditions in two and three dimensions. This work generalizes the analytical virtual crack extension method for linear elastic fracture mechanics presented by Lin and Abel, who introduced the direct integral forms of the energy release rate and its derivatives for a structure containing a two dimensional single crack. Here Lin and Abel’s method is generalized and derivations are provided for verification of the following: extension to the general case of a system of interacting cracks in two dimensions, extension to the axisymmetric case, extension to three-dimensional crack with an arbitrarily curved front under general mixed-mode loading conditions, inclusion of non-uniform crack-face pressure and thermal loading, and an evaluation of the second order derivative of the energy release rate. The method provides the direct integral forms of stiffness derivatives, and thus there is no need for the analyst to specify a finite length of virtual crack extension. The salient feature of this method is that the energy release rates and their higher derivatives for multiple cracks in two and three dimensions can be computed in a single analysis. It is shown that the number of rings of elements surrounding the crack tip that are involved in the mesh perturbation due to the virtual crack extension has an effect on the solution accuracy.