We present a method to generate polyhedral meshes in arbitrary domains that are suitable for use in numerical analysis involving methods such as finite elements or virtual elements. Previous work on polyhedral mesh generation consisted of computing Voronoi tessellations and using Lloyd's algorithm to obtain mesh regularity, i.e., removing excessive element distortion. However, this procedure may produce short edges in 2D (or small facets in 3D) in the Voronoi diagram resulting in numerical instabilities. To address this issue, we optimize the mesh by minimizing a carefully designed energy functional that effectively removes the short edges (and the small facets) in the mesh. We verify the effectiveness of the proposed technique by reducing the condition number of the stiffness matrix associated to a linear elasticity problem. Moreover, we demonstrate that these meshes are effective in the solution of topology optimization problems.