diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
index 00720d2382..57514a2053 100644
--- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict
+++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict
@@ -167,8 +167,10 @@ FoamFile
 //    source regionToCell;
 //    sourceInfo
 //    {
-//        set         c0;         // name of cellSet giving mesh subset
-//        insidePoint (1 2 3);    // point inside region to select
+//        set         c0;         // optional name of cellSet giving mesh subset
+//        insidePoints ((1 2 3)); // points inside region to select
+//        nErode      0;          // optional number of layers to erode
+//                                // selection
 //    }
 //
 //    // Cells underneath plane such that volume is reached. E.g. for use
diff --git a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
index 77948ca5a5..57c7b79af5 100644
--- a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
+++ b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,12 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "regionToCell.H"
-#include "polyMesh.H"
 #include "regionSplit.H"
-#include "globalMeshData.H"
+#include "emptyPolyPatch.H"
 #include "cellSet.H"
 #include "syncTools.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -49,104 +47,339 @@ addToRunTimeSelectionTable(topoSetSource, regionToCell, istream);
 Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
 (
     regionToCell::typeName,
-    "\n    Usage: regionToCell subCellSet (x y z)\n\n"
-    "    Select all cells in the connected region containing point.\n"
-    "    If started inside the subCellSet keeps to it;\n"
-    "    if started outside stays outside.\n"
+    "\n    Usage: regionToCell subCellSet (pt0 .. ptn)\n\n"
+    "    Select all cells in the connected region containing"
+    " points (pt0..ptn).\n"
 );
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::regionToCell::markRegionFaces
+(
+    const boolList& selectedCell,
+    boolList& regionFace
+) const
+{
+    // Internal faces
+    const labelList& faceOwner = mesh_.faceOwner();
+    const labelList& faceNeighbour = mesh_.faceNeighbour();
+    forAll(faceNeighbour, faceI)
+    {
+        if
+        (
+            selectedCell[faceOwner[faceI]]
+         != selectedCell[faceNeighbour[faceI]]
+        )
+        {
+            regionFace[faceI] = true;
+        }
+    }
+
+    // Swap neighbour selectedCell state
+    boolList nbrSelected;
+    syncTools::swapBoundaryCellList(mesh_, selectedCell, nbrSelected);
+
+    // Boundary faces
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+        const labelUList& faceCells = pp.faceCells();
+        forAll(faceCells, i)
+        {
+            label faceI = pp.start()+i;
+            label bFaceI = faceI-mesh_.nInternalFaces();
+            if
+            (
+                selectedCell[faceCells[i]]
+             != selectedCell[nbrSelected[bFaceI]]
+            )
+            {
+                regionFace[faceI] = true;
+            }
+        }
+    }
+}
+
+
+Foam::boolList Foam::regionToCell::findRegions
+(
+    const bool verbose,
+    const regionSplit& cellRegion
+) const
+{
+    boolList keepRegion(cellRegion.nRegions(), false);
+
+    forAll(insidePoints_, i)
+    {
+        // Find the region containing the insidePoint
+
+        label cellI = mesh_.findCell(insidePoints_[i]);
+
+        label keepRegionI = -1;
+        label keepProcI = -1;
+        if (cellI != -1)
+        {
+            keepRegionI = cellRegion[cellI];
+            keepProcI = Pstream::myProcNo();
+        }
+        reduce(keepRegionI, maxOp<label>());
+        keepRegion[keepRegionI] = true;
+
+        reduce(keepProcI, maxOp<label>());
+
+        if (keepProcI == -1)
+        {
+            FatalErrorIn
+            (
+                "outsideCellSelection::findRegions"
+                "(const bool, const regionSplit&)"
+            )   << "Did not find " << insidePoints_[i]
+                << " in mesh." << " Mesh bounds are " << mesh_.bounds()
+                << exit(FatalError);
+        }
+
+        if (verbose)
+        {
+            Info<< "Found location " << insidePoints_[i]
+                << " in cell " << cellI << " on processor " << keepProcI
+                << " in global region " << keepRegionI
+                << " out of " << cellRegion.nRegions() << " regions." << endl;
+        }
+    }
+
+    return keepRegion;
+}
+
+
+void Foam::regionToCell::unselectOutsideRegions
+(
+    boolList& selectedCell
+) const
+{
+    // Determine faces on the edge of selectedCell
+    boolList blockedFace(mesh_.nFaces(), false);
+    markRegionFaces(selectedCell, blockedFace);
+
+    // Determine regions
+    regionSplit cellRegion(mesh_, blockedFace);
+
+    // Determine regions containing insidePoints_
+    boolList keepRegion(findRegions(true, cellRegion));
+
+    // Go back to bool per cell
+    forAll(cellRegion, cellI)
+    {
+        if (!keepRegion[cellRegion[cellI]])
+        {
+            selectedCell[cellI] = false;
+        }
+    }
+}
+
+
+void Foam::regionToCell::shrinkRegions
+(
+    boolList& selectedCell
+) const
+{
+    // Select points on unselected cells and boundary
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    boolList boundaryPoint(mesh_.nPoints(), false);
+
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    forAll(pbm, patchI)
+    {
+        const polyPatch& pp = pbm[patchI];
+
+        if (!pp.coupled() && !isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                const face& f = pp[i];
+                forAll(f, fp)
+                {
+                    boundaryPoint[f[fp]] = true;
+                }
+            }
+        }
+    }
+
+    forAll(selectedCell, cellI)
+    {
+        if (!selectedCell[cellI])
+        {
+            const labelList& cPoints = mesh_.cellPoints(cellI);
+            forAll(cPoints, i)
+            {
+                boundaryPoint[cPoints[i]] = true;
+            }
+        }
+    }
+
+    syncTools::syncPointList(mesh_, boundaryPoint, orEqOp<bool>(), false);
+
+
+    // Select all cells using these points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    label nChanged = 0;
+    forAll(boundaryPoint, pointI)
+    {
+        if (boundaryPoint[pointI])
+        {
+            const labelList& pCells = mesh_.pointCells(pointI);
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+                if (selectedCell[cellI])
+                {
+                    selectedCell[cellI] = false;
+                    nChanged++;
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::regionToCell::erode
+(
+    boolList& selectedCell
+) const
+{
+    //Info<< "Entering shrinkRegions:" << count(selectedCell) << endl;
+    //generateField("selectedCell_before", selectedCell)().write();
+
+    // Now erode and see which regions get disconnected
+    boolList shrunkSelectedCell(selectedCell);
+
+    for (label iter = 0; iter < nErode_; iter++)
+    {
+        shrinkRegions(shrunkSelectedCell);
+    }
+
+    //Info<< "After shrinking:" << count(shrunkSelectedCell) << endl;
+    //generateField("shrunkSelectedCell", shrunkSelectedCell)().write();
+
+
+
+    // Determine faces on the edge of shrunkSelectedCell
+    boolList blockedFace(mesh_.nFaces(), false);
+    markRegionFaces(shrunkSelectedCell, blockedFace);
+
+    // Find disconnected regions
+    regionSplit cellRegion(mesh_, blockedFace);
+
+    // Determine regions containing insidePoints
+    boolList keepRegion(findRegions(true, cellRegion));
+
+
+    // Extract cells in regions that are not to be kept.
+    boolList removeCell(mesh_.nCells(), false);
+    forAll(cellRegion, cellI)
+    {
+        if (shrunkSelectedCell[cellI] && !keepRegion[cellRegion[cellI]])
+        {
+            removeCell[cellI] = true;
+        }
+    }
+
+    //Info<< "removeCell before:" << count(removeCell) << endl;
+    //generateField("removeCell_before", removeCell)().write();
+
+
+
+    // Grow removeCell
+    for (label iter = 0; iter < nErode_; iter++)
+    {
+        // Grow selected cell in regions that are not for keeping
+        boolList boundaryPoint(mesh_.nPoints(), false);
+        forAll(removeCell, cellI)
+        {
+            if (removeCell[cellI])
+            {
+                const labelList& cPoints = mesh_.cellPoints(cellI);
+                forAll(cPoints, i)
+                {
+                    boundaryPoint[cPoints[i]] = true;
+                }
+            }
+        }
+        syncTools::syncPointList(mesh_, boundaryPoint, orEqOp<bool>(), false);
+
+        // Select all cells using these points
+
+        label nChanged = 0;
+        forAll(boundaryPoint, pointI)
+        {
+            if (boundaryPoint[pointI])
+            {
+                const labelList& pCells = mesh_.pointCells(pointI);
+                forAll(pCells, i)
+                {
+                    label cellI = pCells[i];
+                    if (!removeCell[cellI])
+                    {
+                        removeCell[cellI] = true;
+                        nChanged++;
+                    }
+                }
+            }
+        }
+    }
+
+    //Info<< "removeCell after:" << count(removeCell) << endl;
+    //generateField("removeCell_after", removeCell)().write();
+
+
+    // Unmark removeCell
+    forAll(removeCell, cellI)
+    {
+        if (removeCell[cellI])
+        {
+            selectedCell[cellI] = false;
+        }
+    }
+}
+
+
 void Foam::regionToCell::combine(topoSet& set, const bool add) const
 {
-    label cellI = mesh_.findCell(insidePoint_);
+    // Note: wip. Select cells first
+    boolList selectedCell(mesh_.nCells(), true);
 
-    // Load the subset of cells
-    boolList blockedFace(mesh_.nFaces(), false);
+    if (setName_.size() && setName_ != "none")
     {
         Info<< "Loading subset " << setName_ << " to delimit search region."
             << endl;
         cellSet subSet(mesh_, setName_);
 
-        boolList inSubset(mesh_.nCells(), false);
+        selectedCell = false;
         forAllConstIter(cellSet, subSet, iter)
         {
-            inSubset[iter.key()] = true;
-        }
-
-        if (cellI != -1 && inSubset[cellI])
-        {
-            Pout<< "Point " << insidePoint_ << " is inside cellSet "
-                << setName_ << endl
-                << "Collecting all cells connected to " << cellI
-                << " and inside cellSet " << setName_ << endl;
-        }
-        else
-        {
-            Pout<< "Point " << insidePoint_ << " is outside cellSet "
-                << setName_ << endl
-                << "Collecting all cells connected to " << cellI
-                << " and outside cellSet " << setName_ << endl;
-        }
-
-        // Get coupled cell status
-        label nInt = mesh_.nInternalFaces();
-        boolList neiSet(mesh_.nFaces()-nInt, false);
-        for (label faceI = nInt; faceI < mesh_.nFaces(); faceI++)
-        {
-             neiSet[faceI-nInt] = inSubset[mesh_.faceOwner()[faceI]];
-        }
-        syncTools::swapBoundaryFaceList(mesh_, neiSet);
-
-        // Find faces inbetween subSet and non-subset.
-        for (label faceI = 0; faceI < nInt; faceI++)
-        {
-            bool ownInSet = inSubset[mesh_.faceOwner()[faceI]];
-            bool neiInSet = inSubset[mesh_.faceNeighbour()[faceI]];
-            blockedFace[faceI] = (ownInSet != neiInSet);
-        }
-        for (label faceI = nInt; faceI < mesh_.nFaces(); faceI++)
-        {
-            bool ownInSet = inSubset[mesh_.faceOwner()[faceI]];
-            bool neiInSet = neiSet[faceI-nInt];
-            blockedFace[faceI] = (ownInSet != neiInSet);
+            selectedCell[iter.key()] = true;
         }
     }
 
-    // Find connected regions without crossing boundary of the cellset.
-    regionSplit regions(mesh_, blockedFace);
 
-    // Get the region containing the insidePoint
-    label regionI = -1;
+    unselectOutsideRegions(selectedCell);
 
-    if (cellI != -1)
+    if (nErode_ > 0)
     {
-        // On processor that has found cell.
-        regionI = regions[cellI];
-    }
-
-    reduce(regionI, maxOp<label>());
-
-    if (regionI == -1)
-    {
-        WarningIn
-        (
-            "regionToCell::combine(topoSet&, const bool) const"
-        )   << "Point " << insidePoint_
-            << " is not inside the mesh." << nl
-            << "Bounding box of the mesh:" << mesh_.bounds()
-            << endl;
-        return;
+        erode(selectedCell);
     }
 
 
-    // Pick up the cells of the region
-    const labelList regionCells(findIndices(regions, regionI));
-
-    forAll(regionCells, i)
+    forAll(selectedCell, cellI)
     {
-        addOrDelete(set, regionCells[i], add);
+        if (selectedCell[cellI])
+        {
+            addOrDelete(set, cellI, add);
+        }
     }
 }
 
@@ -158,12 +391,14 @@ Foam::regionToCell::regionToCell
 (
     const polyMesh& mesh,
     const word& setName,
-    const point& insidePoint
+    const pointField& insidePoints,
+    const label nErode
 )
 :
     topoSetSource(mesh),
     setName_(setName),
-    insidePoint_(insidePoint)
+    insidePoints_(insidePoints),
+    nErode_(nErode)
 {}
 
 
@@ -175,8 +410,14 @@ Foam::regionToCell::regionToCell
 )
 :
     topoSetSource(mesh),
-    setName_(dict.lookup("set")),
-    insidePoint_(dict.lookup("insidePoint"))
+    setName_(dict.lookupOrDefault<word>("set", "none")),
+    insidePoints_
+    (
+        dict.found("insidePoints")
+      ? dict.lookup("insidePoints")
+      : dict.lookup("insidePoint")
+    ),
+    nErode_(dict.lookupOrDefault("nErode", 0))
 {}
 
 
@@ -189,7 +430,8 @@ Foam::regionToCell::regionToCell
 :
     topoSetSource(mesh),
     setName_(checkIs(is)),
-    insidePoint_(checkIs(is))
+    insidePoints_(checkIs(is)),
+    nErode_(0)
 {}
 
 
@@ -209,15 +451,15 @@ void Foam::regionToCell::applyToSet
 {
     if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
     {
-        Info<< "    Adding all cells of connected region containing point "
-            << insidePoint_ << " ..." << endl;
+        Info<< "    Adding all cells of connected region containing points "
+            << insidePoints_ << " ..." << endl;
 
         combine(set, true);
     }
     else if (action == topoSetSource::DELETE)
     {
-        Info<< "    Removing all cells of connected region containing point "
-            << insidePoint_ << " ..." << endl;
+        Info<< "    Removing all cells of connected region containing points "
+            << insidePoints_ << " ..." << endl;
 
         combine(set, false);
     }
diff --git a/src/meshTools/sets/cellSources/regionToCell/regionToCell.H b/src/meshTools/sets/cellSources/regionToCell/regionToCell.H
index 94a49bbf64..4d740c0ca3 100644
--- a/src/meshTools/sets/cellSources/regionToCell/regionToCell.H
+++ b/src/meshTools/sets/cellSources/regionToCell/regionToCell.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,20 @@ Class
 
 Description
     TopoSetSource. Select cells belonging to topological connected region
-    (that contains given point)
+    (that contains given points)
+
+    In dictionary input:
+
+        // optional name of cellSet delimiting search
+        set         c0;
+
+        //- Number of cell layers to erode mesh to detect holes in the mesh
+        //  Set to 0 if not used.
+        nErode 3;
+
+        // points inside region to select
+        insidePoints ((1 2 3));
+
 
 SourceFiles
     regionToCell.C
@@ -37,12 +50,15 @@ SourceFiles
 #define regionToCell_H
 
 #include "topoSetSource.H"
+#include "boolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+class regionSplit;
+
 /*---------------------------------------------------------------------------*\
                            Class regionToCell Declaration
 \*---------------------------------------------------------------------------*/
@@ -58,14 +74,38 @@ class regionToCell
         static addToUsageTable usage_;
 
         //- Name of cellSet to keep to
-        word setName_;
+        const word setName_;
 
-        //- Coordinate that is inside connected region
-        point insidePoint_;
+        //- Coordinate(s) that is inside connected region
+        const pointField insidePoints_;
+
+        //- Number of layers to erode
+        const label nErode_;
 
 
     // Private Member Functions
 
+        //- Mark faces inbetween selected and unselected elements
+        void markRegionFaces
+        (
+            const boolList& selectedCell,
+            boolList& regionFace
+        ) const;
+
+        //- Determine for every disconnected region in the mesh whether
+        //  it contains a locationInMesh
+        boolList findRegions(const bool verbose, const regionSplit&) const;
+
+        //- Unselect regions not containing a locationInMesh
+        void unselectOutsideRegions(boolList& selectedCell) const;
+
+        //- Unselect one layer of cells from selectedCell
+        void shrinkRegions(boolList& selectedCell) const;
+
+        //- Erode a given number of layers from selectedCell. Remove any
+        //  region that gets disconnected that way.
+        void erode(boolList& selectedCell) const;
+
         void combine(topoSet& set, const bool add) const;
 
 
@@ -81,7 +121,8 @@ public:
         (
             const polyMesh& mesh,
             const word& setName,
-            const point& insidePoint
+            const pointField& insidePoints,
+            const label nErode
         );
 
         //- Construct from dictionary