A finite element method based on the penalized subdomain-interface framework is proposed for fully-coupled, nonlinear thermomechanical analyses with thermal contact anuor radiation boundaries. In the variational formulation, a well-known penalty functional scheme is adopted for connecting subdomains and interfaces that satisfy various continuity requirements. As a logical consequence, the whole domain can be arbitrarily divided into independently-modeled subdomains without considering the conformity of meshes along their interfaces. Since the nonlinearities due to the contact and radiation boundaries can be localized within a few subdomains, the computational efficiency of the present method is greatly increased with appropriate solution algorithms. By solving some numerical problems, these advantageous features are confirmed carefully.