We consider shape optimization problems with elliptic partial differential state equations.Using regularization and penalization, unknown shapes are encoded via shape functions, turning the shape optimization into optimal control problems for the unknown functions. The method is designed to allow topological changes in a natural way. Based on convergence and differentiability results, numerical algorithms are formulated, using different descent directions and projections. The algorithms are assessed in a series of numerical experiments, applied 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.