A subdomain-interface finite element method is suggested to solve a class of fully- coupled thermomechanical problems with contact boundaries. The penalty method is used for connecting subdomains that satisfy interface compatibility conditions. As a result, effective stiffness matrices are always positive definite, and computational efficiency can be improved to a considerable degree. Moreover, any complex-shaped domain can be divided into independently modeled subdomains without considering the conformity of meshes on interfaces. Using a computer code based on the present method, these advantageous features are shown through a set of numerical studies.