A shape design problem (or inverse geometry problem) in determining the geometry of interfacial surface between two conductive bodies in a three-dimensional multiple region domains, based on the desired system heat flux and domain volume, is examined in this study. The design algorithm utilized the Levenberg-Marquardt method (LMM), B-spline surface generation and the commercial software CFD-ACE+. The validity of this shape design analysis is examined using the numerical experiments. Different desired system heat fluxes are considered in the numerical test cases to justify the validity of the present algorithm in solving the three-dimensional shape design problems. Finally, the results show that for the two different cases considered in this work, the maximum increasing in the system heat flux is obtained as 11.3% and 14.1%, respectively. It is also concluded that when the boundary control points of interfacial surface are free to move, maximum system heat flux can be obtained by the present algorithm since it has more degree of freedom in describing the interfacial surface. (c) 2011 Elsevier Ltd. All rights reserved.