Seismic exploration is the primary tool for finding and mapping out hydrocarbon deposits. We have considered the 2D forward modeling problem. The subsurface structures were assumed to be known, and the task was to simulate elastic-wave propagation throughout the medium. Many traditional approaches to solving this problem account for material interfaces by allowing model parameter values (such as wave speed and density) to vary either sharply or smoothly between adjacent data points. Although simple to implement, these strategies typically produce solutions that feature a suboptimal first-order convergence to the true solution. Spatial discretization based on radial-basis-function-generated finite differences (RBF-FDs) has previously been shown to offer high accuracy and algebraic simplicity when using scattered layouts of computational nodes. We have now developed this method further, to provide third-order accuracy not only throughout smooth regions, but also for wave reflections and transmissions at arbitrarily curved material interfaces. The key step is supplementing the RBFs that underlie the RBF-FD approximation with a specific space of piecewise polynomials. The nonsmoothness of these polynomials across an interface is designed to enforce continuity of traction and motion at that interface. The high-order accuracy of the method is illustrated on a couple of test problems for the 2D isotropic elastic-wave equation.