An efficient domain/boundary decomposition method is applied for heat transfer problems with non-linear thermal radiation boundaries. The whole domain of solids or structures is considered as set of subdomains, an interface, and radiation interfaces. In a variational formulation, simple penalty functions are introduced to connect an interface or radiation interfaces with neighboring subdomains that satisfy continuity conditions. As a result, non-linear finite element computations due to the thermal radiation boundaries can be localized within a few subdomains or radiation interfaces. Therefore, by setting up suitable solution algorithms for the governing finite element equations, the computational efficiency can be improved considerably. Through a set of numerical examples, these distinguishing characteristics of the present method are investigated in detail.