We consider discretized two-dimensional PDE-constrained shape optimization problems, in which shapes are represented by triangular meshes. Our contribution is two-fold. First, we study the space of admissible vertex positions of a given connectivity that a triangular mesh can attain. It turns out it is indeed a smooth manifold, and we endow it with a complete Riemannian metric, which allows large mesh deformations without jeopardizing mesh quality. Second, we consider the discrete shape optimization problem of finding optimal vertex positions, which generally does not possess a globally optimal solution. To overcome this ill-posedness, we propose to add a mesh quality penalization term to the objective function. This penalization allows us to simultaneously render the shape optimization problem solvable and keep track of the mesh quality. We also conduct a numerical study of the impact of the choice of Riemannian metric on the steepest descent method.