In engineering applications, it is pretty often to have domain heat source involved inside. This article proposes an approach using the boundary element method to study thermal stresses in 3D anisotropic solids when internal domain heat source is involved. As has been well noticed, thermal effect will give rise to a volume integral, where its direct evaluation will need domain discretization. This shall definitely destroy the most distinctive notion of the boundary element method that only boundary discretization is required. The present work presents an analytical transformation of the volume integral in the boundary integral equation due to the presence of internal volume heat source. For simplicity, distribution of the heat source is modeled by a quadratic function. When needed, the formulations can be further extended to treat higher-ordered volume heat sources. Indeed, the present work has completely restored the boundary discretization feature of the boundary element method for treating 3D anisotropic thermoelasticity involving volume heat source.