This talk considers the numerical solution of elastodynamic wave propagation problems in a bounded polygonal domain with mixed boundary conditions. The dynamics is described by the Navier equation whose unknown represents the displacement field of the elastic body. The solution to this problem exhibits singularities at both corners and points where the boundary condition change [1,2]. The asymptotic behavior is theoretically derived from a detailed study of the Dirichlet trace of the solution and the traction at the boundary. The analysis leads to quasi-optimal error estimates for piecewise polynomial approximations. The numerical solution is performed by a Boundary Element Method, which is based on the representation of the relevant physical fields in terms of layer potentials. The resulting Boundary Integral Equations (BIEs) involve the single layer, double layer and hypersingular boundary integrals operators and provide a convenient reduction of the dimensional complexity. Energy arguments are used to obtain a weak formulation of the boundary integral equations, which are then numerically solved using a space-time Galerkin discretization (E-BEM). The linear system has a block Toeplitz structure with symmetric blocks, arising from a time discretization with Lagrange basis functions. The matrix entries are computed using numerical quadrature strategies for the computation of the well-known space integrals for the BIE system. Results for 2D interior problems [3] and 2D contact problems [4] show the accuracy and long-time stability obtained by the energetic Galerkin approach. Boundary meshes graded towards the singularity for both displacement and traction lead to quasi-optimal numerical approximations, improving the slow convergence and demanding computational costs on uniform meshes. Numerical examples, obtained with the E-BEM implementation, illustrate the theoretical results for 2D polygonal geometries, confirming the expected quasi-optimal convergence rates.