SUMMARYThis paper is concerned with the numerical calculation of stress intensity factors in two‐dimensional linear fracture problems by Petrov–Galerkin natural element method, for which the Voronoi polygon‐based Laplace interpolation functions and constant strain finite element basis functions are taken for the trial and test functions, respectively. A local element patch recovery technique is employed to obtain more accurate smoothed strain and stress fields, and the stress intensity factors of edge and angled center cracks are calculated by both the J‐integral and interaction integral methods. The influence of the grid density, integral domain size, and weighting function type on the numerical accuracy is investigated through the numerical experiments. In addition, the present method is compared with the hybrid crack elements in combination with extended FEM and the scaled boundary FEM. It is confirmed that Petrov–Galerkin natural element method provides the numerical accuracy similar to these recent direct methods. Also, it is found that the accuracy of the J‐integral is affected by the grid density but not by the integral domain size. In the case of the interaction integral, the numerical accuracy is influenced by both the type of weighting function and the integral domain size. Copyright © 2014 John Wiley & Sons, Ltd.