## 1 Introduction  
The propagation of waves in elastic solids is governed by a second-order, hyperbolic differential equation. Waves therefore travel inside bodies with a finite speed, and cause an abrupt change of the local velocity and stress fields across their fronts.  
Numerical methods have to accurately track such moving inner discontinuities. Since standard dis placement-based finite element (FE) schemes adopt continuous interpolation fields to mimic the discontinuous ones inside the modeled domain, non-physical high-frequency oscillations show up around wave fronts; these oscillations are a numerical artifact and need to be filtered out of the solution. Algorithmic treatments of this issue have been proposed in the literature, and they typically consist in artificial viscosity or mesh adaption. Focusing on time integration, algorithms like the \(\alpha-\)method [1] or the generalized \(\alpha-\)method [2] were also devised to damp oscillations. All these methods can reduce in size such spurious effects, but sometimes entail energy dissipation which is, again, non-physical.  
An alternative treatment, whose roots can be traced back to the seminal work of [PERSON] [3; 4], is the flux-corrected transport (FCT) method. The FCT algorithm consists of post-processing a standard FE solution: diffusive and antidiffusive fluxes, the latter being appropriately limited in size, are handled to improve the discrete velocity field around the discontinuities, and to filter out spurious oscillations. This method has been extensively used to simulate the propagation of shock waves in fluids [5; 6]; recently, it has been adopted to simulate traveling stress waves in solids [7]. Noteworthy results have been obtained in [8] through a coupling of a displacement-based FE solution and the FCT algorithm; owing to the adopted structured meshes, results for bodies of arbitrary shapes were based on a partition-of-unity enrichment of the nodal shape functions.  
In this work, to study the dynamics of heterogeneous bodies we propose two enhancements to the frame developed in [8; 9]. First, to simulate the propagation of stress waves inside domains of arbitrary shape, an algorithmic length-scale \(\ell_{\text{FCT}}\) is introduced: this length-scale allows to define local supports of finite size, which are independent of the space discretization, wherein diffusive/antidiffusive fluxes are computed. Second, the rationale behind the computation of diffusive/antidiffusive fluxes close to body or phase boundaries is revisited, so as to permit the treatment of compliant interphases confined along loci of zero measure (i.e. surfaces in three-dimensional frames and lines in two-dimensional frames).  
These two enhancements are of paramount importance when dynamic failure of quasi-brittle polycrystals (like, e.g., polysilicon) needs to be modeled, since damaging phenomena at the micro-scale are incepted as soon as the tensile strength is locally attained. Because of the polycrystal micro-structure, traveling waves are partially reflected by each grain boundary, and eventually lead to complex stress patterns. If the numerical solution does contain the aforementioned fictitious oscillations, the amplitude of the local stress field may be artificially increased and, therefore, damage wrongly started. Through simulations of stress waves traveling along homogeneous as well as heterogeneous (bimaterial) rods, we show that the proposed FE-FCT method can be used to accurately study the evolution of the stress field in quasi-brittle polycrystals.  
As far as notation is concerned, a matrix one will be adopted throughout the whole paper with upper-case and lowercase bold symbols respectively denoting matrices and vectors, a superscript \({}^{\text{T}}\) standing for transpose, and a superposed dot representing time rates.