The paper presents a general solver based on the smoothed finite element method (S-FEM) to solve three dimensional steady and transient heat transfer problems. The present model uses surface integral method in computing the smoothed strains, strictly according to the theory of S-FEM, instead of the volume-weighted average method preferred by other researchers. The present S-FEM solver can thus use quadrilateral elements or even any higher order elements or interpolation techniques. Different kinds of smoothing domains are constructed and all the necessary connectivity is obtained during the construction, which includes node-based smoothing domains, edge-based smoothing domains and face-based smoothing domains. When dealing with the time, both explicit and implicit schemes are implemented. In validation of our new solver, we solve 3D steady and transient heat transfer problems with different boundary condition and complex geometry. From these examples, we can conclude that the present solver can efficiently solve heat transfer problems using any S-FEM models for stable and accurate solutions, as well as solutions of desired properties such as the important upper bound property.