In geophysical numerical models using the finite-element method or its variant, the spectral-element method, to solve seismic wave equations, a mesh is employed to discretize the domain. Generating or adapting a mesh to complex geological properties is a challenging task. To tackle this challenge, we develop an r-adaptivity method to generate or adapt a two-dimensional mesh to a seismic velocity field. Our scheme relies on the optimal transport theory to perform vertices relocation, which generates good-shaped meshes and prevents tangled elements. The mesh adaptation can delineate different regions of interest, like sharp interfaces, salt bodies, and discontinuities. The algorithm has a few user-defined parameters that control the mesh density. With typical seismic velocity examples (e.g., Camembert, SEAM Phase, Marmousi-2), mesh adaptation capability is illustrated within meshes with triangular and quadrilateral elements, commonly employed in seismic codes. Besides its potential use in mesh generation, the method developed can be embedded in seismic inversion workflows like multiscale full waveform inversion to adapt the mesh to the field being inverted without incurring the I/O cost of re-meshing and load rebalancing in parallel computations. The method can be extended to three-dimensional meshes.