New numerical methods for PDE resolution (such as Finite Volumes (FV) or Virtual Element Method (VEM)) open new needs in terms of meshing of domains of interest, and in particular polyhedral meshes have many advantages. One way to build such meshes consists in constructing Restricted Voronoi Diagrams (RVDs) whose boundaries respect the domain of interest. By minimizing a function defined for RVDs, the shapes of cells can be controlled, i.e. elongated according to user-defined directions or adjusted to comply with given aspect ratios (anisotropy) and density variations. In this paper, our contribution is threefold: first, we present a gradient formula for the Voronoi tessellation energy under a continuous anisotropy field. Second, we describe a meshing algorithm based on the optimisation of this function that we validate against state-of-the-art approaches. Finally, we propose a hierarchical approach to speed up our meshing algorithm.