The computational accuracy and efficiency of finite element method and spectral element method (SEM) are investigated thoroughly in time-domain elastic wavefield modeling. The diagonal mass matrices of the FEM and SEM free from matrix inversion are compared comprehensively by making full use of the mass-lumped technique and quadrature rules. We investigate the FEM and SEM based, respectively, on quadrilateral with the polynomials of degrees one and two, and on triangular grids with the polynomials of degrees one and three. Generally, the numerical solutions based on quadrilateral grids have a higher precision than those computed on triangular grids when the same order of polynomials is used. The FEM has a comparable accuracy to the SEM with the same number of interpolant points. In view of the triangular and quadrilateral SEMs, the former suffers from larger computational costs and relatively lower accuracy compared with the latter. Furthermore, the convergence study proves that the triangular SEM produces consistently larger errors than the quadrilateral SEM for any order and element sizes. However, the triangular SEM can adapt to arbitrary complex geometries effectively. In terms of efficiency, the FEM has an efficiency comparable with the SEM on condition that the order of interpolation polynomials is identical. In addition, a perfectly matched layer (PML) boundary condition in variational form is deduced. By introducing four intermediate variables in frequency domain, the PML avoids convolution calculation and obtains an exact solution through inverse Fourier transform in time domain. The numerical examples verify the validity and effectiveness of the PML.