Two different formulations for shape design sensitivity analysis are provided for heat conducting solids. The material derivative concept and the adjoint variable method of analysis are the common techniques utilized in the two formulations. It is found that complementary boundary integral equations arise in the case of the integral sensitivity analysis. Numerical solutions of the primary and adjoint equations are achieved by the BEM, whilst shape configurations are discretized by means of ray and height functions. Several 2-D constrained shape optimization problems are solved by the proposed numerical solution techniques.