We consider shape optimization problems, where the state is governed by elliptic partial differential equations. Using a regularization technique, unknown shapes are encoded via shape functions, turning the shape optimization into optimal control problems for the unknown functions. The method is studied for elliptic PDEs to be solved in an unknown region (to be optimized), where the regularization technique together with a penalty method extends the PDE to a larger fixed domain. Additionally, the method is studied for the optimal layout problem, where the unknown regions determine the coefficients of the state equation. In both cases, the existence of optimal shapes is established for the regularized and for the original problem, with convergence of optimal shapes if the regularization parameter tends to zero. Error estimates are proved for the layout problem. In the context of finite element approximations, convergence and differentiability properties are shown. The method is designed to allow topological changes in a natural way, which is illustrated in a series of numerical experiments, applying the method to an elliptic PDE arising from an oil industry application with two unknown shapes, one giving the region where the PDE is solved, and the other determining the PDE's coefficients.