Failure in geotechnical engineering is often related to tension-induced cracking in geomaterials. In this paper, a coupled meshless method and FEM is developed to analyze the problem of three-dimensional cracking. The radial point interpolation method (RPIM) is used to model cracks in the smeared crack framework with an isotropic damage model. The identification of the meshless region is based on the stress state computed by FEM, and the adaptive coupling of RPIM and FEM is achieved by a direct algorithm. Mesh-bias dependency, which poses difficulties in FEM-based cracking simulations, is circumvented by a crack tracking algorithm. The performance of our scheme is demonstrated by two numerical examples, that is, the four-point bending test on concrete beam and the surface cracks caused by tunnel excavation. progress in the field of smeared crack modeling, some major issues still remain unsolved. Two of these important issues are the local remeshing and the modeling of complex three-dimensional crack topology.As crack is represented by strain localization in the crack band, the crack propagation is assumed to proceed elementwise in the smeared crack model. A conundrum in problems involving strain localization is the mesh dependence of the numerical solutions. Although the mesh size dependence can be mitigated by making use of the fracture energy conservation, the element size still has some influence on the numerical result. If a relatively large element size is used, which is more likely in three-dimensional problems, the crack initiation and propagation path cannot be properly captured. Therefore, the mesh used in crack simulation should be fine enough to obtain a satisfactory result. Because of the requirement of node connectivity, however, the three-dimensional local remeshing in FEM can be a cumbersome work. Some recently developed meshless methods need only node information for the establishment of the system of equations, thus are advantageous in adaptive analysis and crack simulation [13,14]. However, the computational efficiency of most meshless methods is rather low compared with FEM. Obviously, a coupling of these two methods may offer the optimal choice by making full use of their advantages [15,16]. Among the meshless methods, the radial point interpolation method (RPIM) possesses the property of the Kronecker delta function, which makes the enforcement of essential boundary condition and the coupling with FEM feasible [17,18]. Wang et al. [19] proposed a direct coupling approach for the simulation of crack propagation in rock and compacted clay, in which the crack region is treated by RPIM, while FEM is used in the rest of the area. In the coupling method, a crucial issue is the identification of the crack region in order to decide where RPIM should be used.Another problem is the modeling of cracks with complex topology. Within the smeared crack framework, the mesh-bias dependency makes it difficult to handle curved cracks properly because the crack propagation path is found to be strongly affected ...