This paper deals with the antiplane wave propagation in a 2D heterogeneous dissipative medium with complex layer interfaces and irregular topography. The initial boundary value problem which represents the viscoelastic dynamics driving 2D antiplane wave propagation is formulated. The discretization scheme is based on the finite-difference technique. Our approach presents some innovative features. First, the introduction of the forcing term into the equation of motion offers the advantage of an easier handling of different inputs such as general functions of spatial coordinates and time. Second, in the case of a straight-line source, the symmetry of the incident plane wave allows us to solve the problem of oblique incidence simply by rotating the 2D model. This artifice reduces the oblique incidence to the vertical one. Third, the conventional rheological model of the generalized Maxwell body has been extended to include the stress-free boundary condition. For this reason we solve explicitly the stress-free boundary condition, not following the most popular technique called vacuum formalism. Finally, our numerical code has been constructed to model the seismic response of complex geological structures: real geological interfaces are automatically digitized and easily introduced in the input model. Three numerical applications are discussed. To validate our numerical model, the first test compares the results of our code with others shown in the literature. The second application rotates the input model to simulate the oblique incidence. The third one deals with a real high-complexity 2D geological structure.