1. A computer-implemented method for generating a constrained Voronoi grid for a planar domain with boundaries and internal features, the planar domain representing at least part of a hydrocarbon reservoir, the method comprising the following steps, at least one of which is implemented using a computer: approximating the boundaries and internal features of the domain with polylines; constructing constrained Delaunay triangulation constrained to the polylines; constructing protected points around the polylines from the vertices of the constrained Delaunay triangulation lying on the polylines; constructing new unconstrained Delaunay triangulation of the set of protected points and interior vertices of the constrained Delaunay triangulation; and constructing the constrained Voronoi grid from the new unconstrained Delaunay triangulation; wherein the construction of protected points comprises: replacing each point lying on the feature by a mirror image pair of protected points; and replacing each point lying on a feature intersection by a set of points placed on a circle which is centered at the intersection.