A method for generating constrained Voronoi grids in a plane with internal features and boundaries is disclosed. The disclosed method generally includes approximation of internal features and boundaries with polylines based on plane geometry. Protected polygons or points are generated around the polylines, and Delaunay triangulation of protected points or protected polygon vertices is constructed. Delaunay triangulation that honors protected polygons or points is generated in the rest of the gridding domain. The constrained Voronoi grid is then generated from the Delaunay triangulation, which resolves all of the approximated features and boundaries with the edges of Voronoi cells. Constrained Voronoi grids may be generated with adaptive cell sizes based on specified density criterion.