When using boundary element methods for problems with domain effects such as initial/residual stresses, body forces and heating, domain integrals are usually involved. Different numerical algorithms have been developed in the past to deal with these domain integrals, such as interior cell mapping, DRM, MRM, etc., which require considerable code development efforts. To take advantage of the boundary element method (BEM) for simulating fracturing processes, but without additional effort for domain integral calculations or transformations, a hybrid approach is suggested in this paper to simulate the excavation and thermo-mechanical process of rock structures. Equivalent boundary tractions reflecting the combined effects of the initial and redistributed thermo-mechanical stresses in the domain of interest at multiple excavation and heating steps are produced by an algorithm of inverse stress reconstruction, in a sense of least-square optimization. The algorithm utilizes resultant stress fields at each excavation and heating step, either measured or produced by other numerical methods of domain-types such as finite element method/finite difference method (FEM/FDM), and calculates the equivalent boundary tractions for further fracturing analysis. In this paper, the inverse stress reconstruction algorithm developed for this purpose, using a truncated singular value decomposition (SVD) technique, is presented. Three examples of verification, two closed-form solutions of problems of elasticity with simple geometry and one example from a predictive modeling of a planned in situ heating experiment in a rock pillar at the Äspö Hard Rock Laboratory (HRL), Southern Sweden, are presented.