In this article we study the numerical solution of elliptic shape optimization problems with additional constraints, given by domain or boundary integral functionals. A special boundary variational approach combined with a boundary integral formulation of the state equation yields shape gradients and functionals which are expressed only in terms of boundary integrals . Hence, the efficiency of (standard) descent optimization algorithms is considerably increased, especially for the line search. We demonstrate our method for a class of problems from planar elasticity, where the stationary domains are given analytically by Banichuk and Karihaloo in [N.V. Banichuk and B.L. Karihaloo (1976). Minimum-weight design of multi-purpose cylindrical bars. International Journal of Solids and Structures , 12 , 267-273.]. In particular, the boundary integral equation is solved by a wavelet Galerkin scheme which offers a powerful tool. For optimization we apply gradient and Quasi-Newton type methods for the penalty as well as for the augmented Lagrangian functional.