Fault slip estimation from crustal deformation is essential for understanding earthquake mechanisms, and three-dimensional subsurface models have become incorporated in the estimation to account for the geometric and material heterogeneity in the crust. In the estimation process, a precomputed Green’s function library (GFL), which represents the displacement field due to unit point loads at each receiver on ground surface, can be employed for computing arbitrary source responses through convolution of the source terms and the GFL. However, challenges remain in generating meshes adapted to the singularity of the point loads and in assessing the accuracy of the GFL. We developed a GFL computation method, Adaptive Finite Element Method for Green’s Function Library (AFEM-GFL). By combining initial meshes with good element quality and a mesh refinement algorithm that is resistant to element quality degradation, our method can generate meshes well adapted to computing GFLs. The accuracy of the results is assessed through convergence of the solution and comparisons between the convergent solutions from different initial meshes. In numerical experiments on a two-layered half-space and a realistic crustal structure of Japan, accurate and convergent GFLs were obtained with moderate amount of computational resources in the settings where it was difficult to achieve with uniform meshes even using a massively parallel supercomputer. These quality GFLs will serve as a robust foundation for reliable fault slip estimation.