Get drainage basin of a polygon (with QGIS, GRASS or other FOSS)

Using a polygon shapefile and a raster DEM, I'd like to find a way to get the drainage basin of a polygon area — that is, the area of all cells which drain into or through the cells contained in the polygon. This is essentially the function provided by r.water.outlet in GRASS except that it only looks at the drainage basin of a single cell.

I figure there must be some way to do this, if only by summing together the drainage basins of all cells within the polygon. However this would require a means of identifying and iterating through all cells contained in a polygon and running r.water.outlet with each cell as a unique input. Ideally I would be able to automate this process with python scripting. Ultimately, I'd like to run it for a large number of polygons in a shapefile layer.

Is this possible?

EDIT: This work relates to semi-automated assessment of sinkholes, karst features and groundwater susceptibility with GRASS GIS & other FOSS. To illustrate what I am hoping to accomplish here, I've taken some screenshots.

A slopeshade raster overlaid on a map of low-lying areas and depressions (shown in blue, derived with r.neighbors). You can see a number of sinkholes here.

Given a vector layer containing polygon features representing the outlines of sinkholes, I'd like to calculate the geometry of the basin that drains into the sinkhole, which will generally be a very small fragment of the larger watershed basin.

Here is a map showing drainage basins for each sinkhole, calculated manually with r.water.outlet

Here you can see a split basin. The sinkhole contains two cells which are separate drainage endpoints in the DEM. I want to get the 'complete' drainage basin for a sinkhole polygon feature, that is the composite of the drainage basin for all points within the sinkhole polygon.