Elastic full-waveform inversion (EFWI) is a process used to estimate subsurface properties by fitting seismic data while satisfying wave propagation physics. The problem is formulated as a least-squares data fitting minimization problem with two sets of constraints: partial-differential equation (PDE) constraints governing elastic wave propagation and physical model constraints implementing prior information. The alternating direction method of multipliers is used to solve the problem, resulting in an iterative algorithm with well-conditioned subproblems. Although wavefield reconstruction is the most challenging part of the iteration, sparse linear algebra techniques can be used for moderate-sized problems and frequency domain formulations. The Hessian matrix is blocky with diagonal blocks, making model updates fast. Gradient ascent is used to update Lagrange multipliers by summing PDE violations. Various numerical examples are used to investigate algorithmic components, including model parameterizations, physical model constraints, the role of the Hessian matrix in suppressing interparameter cross-talk, computational efficiency with the source sketching method, and the effect of noise and near-surface effects.