Cardiac activation mapping requires computing wavefront arrival times on patient-specific geometries with fiber anisotropy, often from multiple stimulation sites firing at different times. The heat method computes geodesic distances via two symmetric positive-definite (SPD) linear solves, but is designed for isotropic media and only supports a single activation time. We extend the heat method to handle both anisotropic conduction tensors and heterogeneous Dirichlet conditions by connecting the eikonal-diffusion equation to a screened Poisson problem. We call the resulting algorithm the multi-front heat method. Both linear systems remain SPD, enabling fast Cholesky or algebraic multigrid solvers. On a 2D anisotropic ring with five activation sites, the method achieves 5.5% relative L2 error against a monodomain reference, compared to 8.1% for the fast iterative method. On a 3D left-ventricular mesh (111K nodes) with fiber-based 3:1 anisotropy and three activation sites, it achieves 7.2% error versus 13.1\% for FIM, with an solve time of 2.1 s. The proposed method offers a competitive alternative for the fast simulation of cardiac activation times.