A numerical technique for wave-propagation simulation in 2-D heterogeneous anisotropic structures is presented. The scheme is flexible in incorporating arbitrary surface topography, inner openings, liquid/solid boundaries, and irregular interfaces, and it naturally satisfies the free-surface conditions of complex geometrical boundaries. The algorithm, based on a discretization mesh of triangles and quadrilaterals, solves the problem using integral equilibrium around each node instead of satisfying elastodynamic differential equations at each node as in the finite-difference method. This study is an extension of previous work for the elastic-isotropic case. Besides accounting for anisotropy, a simplified quadrilateral grid cell with low computational cost is introduced. The transversely isotropic medium with a symmetry axis on the horizontal or vertical plane, as typically caused by a system of parallel cracks or fine layers, is discussed in detail. A 2-D algorithm is presented that can handle the situation where the symmetry axis of the anisotropy does not lie in the 2-D plane. The proposed scheme is successfully tested against an analytical solution for Lamb's problem with a symmetry axis normal to the surface and agrees well with a numerical solution of the reflectivity method for a plane-layered model in the isotropic case. Computed radiation patterns show characteristics such as shear-wave splitting and triplications of quasi-SV wavefronts, as predicted by the theory. Examples of surface-wave propagation in an anisotropic half-space with a semicylindrical pit on the surface and mixed liquid/(anisotropic) solid model with an inclined liquid/solid interface are presented. Moreover, seismograms are modeled for dome-layered and plane-layered anisotropic structures.