Shape optimization problems constrained by partial differential equations (PDEs) are ever-present in the area of scientific computing, for example structural mechanics, computational fluid dynamics and acoustics. The PDEs considered in such applications spans a broad spectrum, ranging from steady state, linear equations, to transient non-linear equations.

To efficiently solve the shape optimization problem, one has to consider shape sensitivities and shape Hessians of the PDEs. A popular approach is the adjoint method. Deriving the adjoint equations by hand is known to be cumbersome and error-prone. We propose an approach where this derivation is computed by a high level algorithmic differentiation tool. The shape sensitivities are computed by discrete derivatives of the mesh node sensitivities.

An advantage of our approach is that only the forward problem has to be postulated as input. Then the algorithm creates a computational tape that tracks the propagation of variables and computes corresponding gradients and Hessians.

By using the operator overloading approach we inherit the parallelism and performance of the software used to solve the PDE. Our software is overloading the framework FEniCS, which has a high level Python interface and generates parallel, optimized C++ code. We illustrate the efficiency and robustness of the code by presenting results for several PDEs with a wide range of solution methods.