  This work develops an algorithm for PDE-constrained shape optimization based
on Lipschitz transformations. Building on previous work in this field, the
$p$-Laplace operator is utilized to approximate a descent method for Lipschitz
shapes. In particular, it is shown how geometric constraints are
algorithmically incorporated avoiding penalty terms by assigning them to the
subproblem of finding a suitable descent direction. A special focus is placed
on the scalability of the proposed methods for large scale parallel computers
via the application of multigrid solvers. The preservation of mesh quality
under large deformations, where shape singularities have to be smoothed or
generated within the optimization process, is also discussed. It is shown that
the interaction of hierarchically refined grids and shape optimization can be
realized by the choice of appropriate descent directions. The performance of
the proposed methods is demonstrated for energy dissipation minimization in
fluid dynamics applications.
