A novel computational framework is established for modeling normal modes, surface waves and time-harmonic body waves with a combination of several highly parallel algorithms. To deal with complex geological features such as topography and interior discontinuities, both forward modeling and inverse modeling are performed on unstructured, deformable tetrahedral meshes. To study the inverse boundary value problem for time-harmonic elastic waves, for the recovery of P- and S-wave speeds from vibroseis data or the Neumann-to-Dirichlet map, a procedure for full waveform inversion with iterative regularization is designed. The multi-level iterative regularization is implemented by projecting gradients, after scaling, onto subspaces to avoid over-parametrization yielding conditional Lipschitz stability. The procedure is illustrated in computational experiments recovering the rough shapes and wave speeds of geological bodies from simple starting models, near and far from the boundary, that is the free surface. To study the seismic spectra, a Rayleigh-Ritz with mixed continuous Galerkin finite-element method based approach is developed to compute the normal modes of a planet in the presence of an essential spectrum. The relevant generalized eigenvalue problem is solved by a Lanczos approach combined with polynomial filtering and separation of the essential spectrum. Self-gravitation is treated as an N-body problem and the relevant gravitational potential is evaluated directly and efficiently utilizing the fast multipole method. In contrast with the standard shift-and-invert and the full-mode coupling algorithms, the polynomial filtering technique is ideally suited for solving large-scale three-dimensional interior eigenvalue problems since it significantly enhances the memory and computational efficiency without loss of accuracy. The parallel efficiency and scalability of the proposed approach are demonstrated on several world-class supercomputers. To include the effects of the rotation and solve the resulting quadratic eigenvalue problem, the extended Lanczos vectors computed from a non-rotating planet are utilized as a subspace to reduce the dimension of the original problem significantly. The reduced system can further be solved via a standard eigensolver. Several computational experiments are performed to study the effects of the normal modes due to three-dimensional fine variations, rotation, and Coriolis effects.