An enriched meshless method, using meshless interpolations and a global Galerkin approach, is developed for the analysis of three-dimensional fracture problems. The displacement field which accounts for the stress singularity nearby the crack front and the boundary layer effect at the intersection between the crack front and the free surface of the structure is adopted to enrich the trial functions. The three-dimensional stress intensity factors can thus be treated as independent unknown parameters and calculated with the nodal displacements directly. To estimate the accuracy of the method developed, several representative three-dimensional cracks are analyzed. These include single-edge crack, embedded elliptical and semi-elliptical surface cracks, and quarter-elliptical corner crack, etc. The variation of the three-dimensional stress intensity factors along the crack front is drawn in details. The influence of crack sizes and boundary layer effect on the three-dimensional stress intensity factors is also studied. Excellent agreements between the calculated results and those available in the literature demonstrate the high accuracy and applicability of the method developed.