A non-hypersingular traction boundary integral equation method (BIEM) is proposed for the treatment of crack systems in piezoelectric or anisotropic plane domains loaded by time-harmonic waves. The solution is based on the frequency dependent fundamental solution obtained by Radon transform. The proposed method is flexible, numerically efficient and has virtually no limitations regarding the material type, crack geometry and type of wave loading. The accuracy and convergence of the BIEM solution for stress intensity factors is validated by comparison with existing results from the literature. Simulations for different crack configurations such as coplanar, collinear or cracks in arbitrary position to each other are presented and discussed. They demonstrate among others the strong effect of electromechanical coupling, show the frequency dependent shielding and amplification resulting from crack interaction and reveal the sensitivity of the K-factors on the complex influence of both wave-crack and crack-crack interaction.