We have developed a 3D elastic full-waveform inversion (FWI) algorithm with forward modeling and inversion performed in the frequency domain. The Helmholtz equation is solved with a second-order finite-difference method using an iterative solver equipped with an efficient complex-shifted incomplete LU-based preconditioner. The inversion is based on the minimization of the data misfit functional and a total variation regularization for the unknown model parameters. We implement the Gauss-Newton method as the optimization engine for the inversions. The codes are parallelized with a message passing interface based on the number of shots and receivers. We examine the performance of this elastic FWI algorithm and workflow on synthetic examples including surface seismic and vertical seismic profile configurations. With various initial models, we manage to obtain high-quality velocity images for 3D earth models.