A double-layer interpolation method (DLIM) is proposed by combining the conventional polynomial element interpolation and the radial basis function (RBF) interpolation. The boundary face method (BFM) is coupled with the DLIM to solve potential problems in two dimensions. The nodes of a DLIM element can be classified into two categories: source nodes (which are inside the element) and virtual nodes (which are on the vertices or edges of the element). If we ignore the virtual nodes, this element exactly becomes a standard discontinuous element. However, source nodes together with virtual nodes comprise the DLIM element which is equivalent to a conventional continuous element. This continuous element (first layer interpolation) will be used to interpolate boundary variables, but the boundary integral equations (BIEs) are collocated at the source nodes only. By building relationships between source nodes and virtual nodes using the RBF (second layer interpolation), additional constraint equations are provided. With these equations, a square matrix for the overall system of linear algebraic equations can be finally obtained. For the same number of source nodes, our method will significantly increase the computational accuracy compared with using discontinuous element. Moreover, both continuous fields and discontinuous fields can be accurately and naturally interpolated by manipulating the supporting domains of some specific virtual nodes in the RBF. Based on this feature, the DLIM can be easily extended to solve problems involving discontinuity, such as elastoplasticity problems, crack propagation, etc. Numerical results on curve fitting indicate that the computational accuracy and the convergence rate of the DLIM are much higher than those of the conventional continuous element used in the BFM.