diff --git a/applications/utilities/miscellaneous/patchSummary/patchSummary.C b/applications/utilities/miscellaneous/patchSummary/patchSummary.C
index d344845fae..d23acc6a69 100644
--- a/applications/utilities/miscellaneous/patchSummary/patchSummary.C
+++ b/applications/utilities/miscellaneous/patchSummary/patchSummary.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
@@ -28,6 +28,11 @@ Description
     Writes fields and boundary condition info for each patch at each requested
     time instance.
 
+    Default action is to write a single entry for patches/patchGroups with the
+    same boundary conditions. Use the -expand option to print every patch
+    separately. In case of multiple groups matching it will print only the
+    first one.
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
@@ -42,12 +47,22 @@ int main(int argc, char *argv[])
     timeSelector::addOptions();
 
 #   include "addRegionOption.H"
+    argList::addBoolOption
+    (
+        "expand",
+        "Do not combine patches"
+    );
 #   include "setRootCase.H"
 #   include "createTime.H"
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
+    const bool expand = args.optionFound("expand");
+
+
 #   include "createNamedMesh.H"
+    const polyBoundaryMesh& bm = mesh.boundaryMesh();
+
 
     forAll(timeDirs, timeI)
     {
@@ -55,6 +70,14 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        // Update the mesh if changed
+        if (mesh.readUpdate() == polyMesh::TOPO_PATCH_CHANGE)
+        {
+            Info<< "Detected changed patches. Recreating patch group table."
+                << endl;
+        }
+
+
         const IOobjectList fieldObjs(mesh, runTime.timeName());
         const wordList objNames = fieldObjs.names();
 
@@ -88,16 +111,87 @@ int main(int argc, char *argv[])
 
         Info<< endl;
 
-        const polyBoundaryMesh& bm = mesh.boundaryMesh();
-        forAll(bm, patchI)
+
+        if (expand)
         {
-            Info<< bm[patchI].type() << ": " << bm[patchI].name() << nl;
-            outputFieldList<scalar>(vsf, patchI);
-            outputFieldList<vector>(vvf, patchI);
-            outputFieldList<sphericalTensor>(vsptf, patchI);
-            outputFieldList<symmTensor>(vsytf, patchI);
-            outputFieldList<tensor>(vtf, patchI);
-            Info<< endl;
+            // Print each patch separately
+
+            forAll(bm, patchI)
+            {
+                Info<< bm[patchI].type() << "\t: " << bm[patchI].name() << nl;
+                outputFieldList<scalar>(vsf, patchI);
+                outputFieldList<vector>(vvf, patchI);
+                outputFieldList<sphericalTensor>(vsptf, patchI);
+                outputFieldList<symmTensor>(vsytf, patchI);
+                outputFieldList<tensor>(vtf, patchI);
+                Info<< endl;
+            }
+        }
+        else
+        {
+            // Collect for each patch the bc type per field. Merge similar
+            // patches.
+
+            // Per 'group', the map from fieldname to patchfield type
+            DynamicList<HashTable<word> > fieldToTypes(bm.size());
+            // Per 'group' the patches
+            DynamicList<DynamicList<label> > groupToPatches(bm.size());
+            forAll(bm, patchI)
+            {
+                HashTable<word> fieldToType;
+                collectFieldList<scalar>(vsf, patchI, fieldToType);
+                collectFieldList<vector>(vvf, patchI, fieldToType);
+                collectFieldList<sphericalTensor>(vsptf, patchI, fieldToType);
+                collectFieldList<symmTensor>(vsytf, patchI, fieldToType);
+                collectFieldList<tensor>(vtf, patchI, fieldToType);
+
+                label groupI = findIndex(fieldToTypes, fieldToType);
+                if (groupI == -1)
+                {
+                    DynamicList<label> group(1);
+                    group.append(patchI);
+                    groupToPatches.append(group);
+                    fieldToTypes.append(fieldToType);
+                }
+                else
+                {
+                    groupToPatches[groupI].append(patchI);
+                }
+            }
+
+
+            forAll(groupToPatches, groupI)
+            {
+                const DynamicList<label>& patchIDs = groupToPatches[groupI];
+
+                if (patchIDs.size() > 1)
+                {
+                    // Check if part of a group
+                    wordList groups;
+                    labelHashSet nonGroupPatches;
+                    bm.matchGroups(patchIDs, groups, nonGroupPatches);
+
+                    const labelList sortedPatches(nonGroupPatches.sortedToc());
+                    forAll(sortedPatches, i)
+                    {
+                        Info<< bm[sortedPatches[i]].type()
+                            << "\t: " << bm[sortedPatches[i]].name() << nl;
+                    }
+                    if (groups.size())
+                    {
+                        forAll(groups, i)
+                        {
+                            Info<< "group\t: " << groups[i] << nl;
+                        }
+                    }
+                    outputFieldList<scalar>(vsf, patchIDs[0]);
+                    outputFieldList<vector>(vvf, patchIDs[0]);
+                    outputFieldList<sphericalTensor>(vsptf, patchIDs[0]);
+                    outputFieldList<symmTensor>(vsytf, patchIDs[0]);
+                    outputFieldList<tensor>(vtf, patchIDs[0]);
+                    Info<< endl;
+                }
+            }
         }
     }
 
diff --git a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C
index 7d03b999a4..9c9842bd32 100644
--- a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.C
+++ b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.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
@@ -55,7 +55,7 @@ void Foam::addToFieldList
 template<class Type>
 void Foam::outputFieldList
 (
-    PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+    const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
     const label patchI
 )
 {
@@ -71,4 +71,26 @@ void Foam::outputFieldList
 }
 
 
+template<class Type>
+void Foam::collectFieldList
+(
+    const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+    const label patchI,
+    HashTable<word>& fieldToType
+)
+{
+    forAll(fieldList, fieldI)
+    {
+        if (fieldList.set(fieldI))
+        {
+            fieldToType.insert
+            (
+                fieldList[fieldI].name(),
+                fieldList[fieldI].boundaryField()[patchI].type()
+            );
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.H b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.H
index 707c7d7517..d61a416a52 100644
--- a/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.H
+++ b/applications/utilities/miscellaneous/patchSummary/patchSummaryTemplates.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
@@ -45,9 +45,17 @@ namespace Foam
     template<class Type>
     void outputFieldList
     (
-        PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+        const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
         const label patchI
     );
+
+    template<class Type>
+    void collectFieldList
+    (
+        const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fieldList,
+        const label patchI,
+        HashTable<word>& fieldToType
+    );
 } // End namespace Foam
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index d5670fa331..a54438f105 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.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
@@ -792,6 +792,52 @@ Foam::labelHashSet Foam::polyBoundaryMesh::patchSet
 }
 
 
+void Foam::polyBoundaryMesh::matchGroups
+(
+    const labelUList& patchIDs,
+    wordList& groups,
+    labelHashSet& nonGroupPatches
+) const
+{
+    // Current matched groups
+    DynamicList<word> matchedGroups(1);
+
+    // Current set of unmatched patches
+    nonGroupPatches = labelHashSet(patchIDs);
+
+    const HashTable<labelList, word>& groupPatchIDs = this->groupPatchIDs();
+    for
+    (
+        HashTable<labelList,word>::const_iterator iter =
+            groupPatchIDs.begin();
+        iter != groupPatchIDs.end();
+        ++iter
+    )
+    {
+        // Store currently unmatched patches so we can restore
+        labelHashSet oldNonGroupPatches(nonGroupPatches);
+
+        // Match by deleting patches in group from the current set and seeing
+        // if all have been deleted.
+        labelHashSet groupPatchSet(iter());
+
+        label nMatch = nonGroupPatches.erase(groupPatchSet);
+
+        if (nMatch == groupPatchSet.size())
+        {
+            matchedGroups.append(iter.key());
+        }
+        else if (nMatch != 0)
+        {
+            // No full match. Undo.
+            nonGroupPatches.transfer(oldNonGroupPatches);
+        }
+    }
+
+    groups.transfer(matchedGroups);
+}
+
+
 bool Foam::polyBoundaryMesh::checkParallelSync(const bool report) const
 {
     if (!Pstream::parRun())
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 999b3815bf..34f0a6d409 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.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
@@ -195,6 +195,15 @@ public:
             const bool usePatchGroups = true
         ) const;
 
+        //- Match the patches to groups. Returns all the (fully matched) groups
+        //  and any remaining unmatched patches.
+        void matchGroups
+        (
+            const labelUList& patchIDs,
+            wordList& groups,
+            labelHashSet& nonGroupPatches
+        ) const;
+
         //- Check whether all procs have all patches and in same order. Return
         //  true if in error.
         bool checkParallelSync(const bool report = false) const;
diff --git a/src/postProcessing/functionObjects/field/Make/files b/src/postProcessing/functionObjects/field/Make/files
index e351784bdc..0abf314f96 100644
--- a/src/postProcessing/functionObjects/field/Make/files
+++ b/src/postProcessing/functionObjects/field/Make/files
@@ -41,4 +41,7 @@ wallBoundedStreamLine/wallBoundedParticle.C
 surfaceInterpolateFields/surfaceInterpolateFields.C
 surfaceInterpolateFields/surfaceInterpolateFieldsFunctionObject.C
 
+regionSizeDistribution/regionSizeDistribution.C
+regionSizeDistribution/regionSizeDistributionFunctionObject.C
+
 LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
new file mode 100644
index 0000000000..b3aaac0dbf
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -0,0 +1,588 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSizeDistribution.H"
+#include "volFields.H"
+#include "regionSplit.H"
+#include "fvcVolumeIntegrate.H"
+#include "Histogram.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(regionSizeDistribution, 0);
+
+    //- plus op for FixedList<scalar>
+    template<class T, unsigned Size>
+    class ListPlusEqOp
+    {
+        public:
+        void operator()
+        (
+            FixedList<T, Size>& x,
+            const FixedList<T, Size>& y
+        ) const
+        {
+            forAll(x, i)
+            {
+                x[i] += y[i];
+            }
+        }
+    };
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionSizeDistribution::regionSizeDistribution
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    name_(name),
+    obr_(obr),
+    active_(true),
+    alphaName_(dict.lookup("field")),
+    patchNames_(dict.lookup("patches"))
+{
+    // Check if the available mesh is an fvMesh, otherwise deactivate
+    if (!isA<fvMesh>(obr_))
+    {
+        active_ = false;
+        WarningIn
+        (
+            "regionSizeDistribution::regionSizeDistribution"
+            "(const objectRegistry&, const dictionary&)"
+        )   << "No fvMesh available, deactivating." << nl
+            << endl;
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionSizeDistribution::~regionSizeDistribution()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::regionSizeDistribution::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        dict.lookup("field") >> alphaName_;
+        dict.lookup("patches") >> patchNames_;
+        dict.lookup("threshold") >> threshold_;
+        dict.lookup("volFraction") >> volFraction_;
+        dict.lookup("nBins") >> nBins_;
+
+        word format(dict.lookup("setFormat"));
+        formatterPtr_ = writer<scalar>::New(format);
+    }
+}
+
+
+void Foam::regionSizeDistribution::execute()
+{
+    // Do nothing - only valid on write
+}
+
+
+void Foam::regionSizeDistribution::end()
+{
+    // Do nothing - only valid on write
+}
+
+
+void Foam::regionSizeDistribution::write()
+{
+    if (active_)
+    {
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+        autoPtr<volScalarField> alphaPtr;
+        if (obr_.foundObject<volScalarField>(alphaName_))
+        {
+            Info<< "Looking up field " << alphaName_ << endl;
+        }
+        else
+        {
+            Info<< "Reading field " << alphaName_ << endl;
+            alphaPtr.reset
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        alphaName_,
+                        mesh.time().timeName(),
+                        mesh,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    mesh
+                )
+            );
+        }
+
+
+        const volScalarField& alpha =
+        (
+             alphaPtr.valid()
+           ? alphaPtr()
+           : obr_.lookupObject<volScalarField>(alphaName_)
+        );
+
+        Info<< "Volume of alpha = "
+            << fvc::domainIntegrate(alpha).value()
+            << endl;
+
+        const scalar meshVol = gSum(mesh.V());
+        Info<< "Mesh volume = " << meshVol << endl;
+        Info<< "Background region volume limit = " << volFraction_*meshVol
+            << endl;
+
+
+        // Determine blocked faces
+        boolList blockedFace(mesh.nFaces(), false);
+        label nBlocked = 0;
+
+        {
+            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+            {
+                scalar ownVal = alpha[mesh.faceOwner()[faceI]];
+                scalar neiVal = alpha[mesh.faceNeighbour()[faceI]];
+
+                if
+                (
+                    (ownVal < threshold_ && neiVal > threshold_)
+                 || (ownVal > threshold_ && neiVal < threshold_)
+                )
+                {
+                    blockedFace[faceI] = true;
+                    nBlocked++;
+                }
+            }
+
+            // Block coupled faces
+            forAll(alpha.boundaryField(), patchI)
+            {
+                const fvPatchScalarField& fvp = alpha.boundaryField()[patchI];
+                if (fvp.coupled())
+                {
+                    tmp<scalarField> townFld(fvp.patchInternalField());
+                    const scalarField& ownFld = townFld();
+                    tmp<scalarField> tnbrFld(fvp.patchNeighbourField());
+                    const scalarField& nbrFld = tnbrFld();
+
+                    label start = fvp.patch().patch().start();
+
+                    forAll(ownFld, i)
+                    {
+                        scalar ownVal = ownFld[i];
+                        scalar neiVal = nbrFld[i];
+
+                        if
+                        (
+                            (ownVal < threshold_ && neiVal > threshold_)
+                         || (ownVal > threshold_ && neiVal < threshold_)
+                        )
+                        {
+                            blockedFace[start+i] = true;
+                            nBlocked++;
+                        }
+                    }
+                }
+            }
+        }
+
+
+        regionSplit regions(mesh, blockedFace);
+
+        Info<< "Determined " << regions.nRegions() << " disconnected regions"
+            << endl;
+
+
+        if (debug)
+        {
+            volScalarField region
+            (
+                IOobject
+                (
+                    "region",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimless, 0)
+            );
+            Info<< "Dumping region as volScalarField to " << region.name()
+                << endl;
+
+            forAll(regions, cellI)
+            {
+                region[cellI] = regions[cellI];
+            }
+            region.correctBoundaryConditions();
+            region.write();
+        }
+
+
+        // Sum all regions
+        Map<Pair<scalar> > regionVolume(regions.nRegions()/Pstream::nProcs());
+        forAll(alpha, cellI)
+        {
+            scalar cellVol = mesh.V()[cellI];
+            scalar alphaVol = alpha[cellI]*cellVol;
+
+            label regionI = regions[cellI];
+
+            Map<Pair<scalar> >::iterator fnd = regionVolume.find(regionI);
+            if (fnd == regionVolume.end())
+            {
+                regionVolume.insert
+                (
+                    regionI,
+                    Pair<scalar>(cellVol, alphaVol)
+                );
+            }
+            else
+            {
+                fnd().first() += cellVol;
+                fnd().second() += alphaVol;
+            }
+        }
+        Pstream::mapCombineGather(regionVolume, ListPlusEqOp<scalar, 2>());
+        Pstream::mapCombineScatter(regionVolume);
+
+
+        if (debug)
+        {
+            Info<< token::TAB << "Region"
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            scalar meshSumVol = 0.0;
+            scalar alphaSumVol = 0.0;
+
+            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            {
+                Info<< token::TAB << iter.key()
+                    << token::TAB << iter().first()
+                    << token::TAB << iter().second() << endl;
+
+                meshSumVol += iter().first();
+                alphaSumVol += iter().second();
+            }
+            Info<< token::TAB << "Total:"
+                << token::TAB << meshSumVol
+                << token::TAB << alphaSumVol << endl;
+            Info<< endl;
+        }
+
+
+
+        // Mark all regions starting at patches
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Count number of patch faces (just for initial sizing)
+        label nPatchFaces = 0;
+        forAll(patchNames_, i)
+        {
+            const word& pName = patchNames_[i];
+            label patchI = mesh.boundaryMesh().findPatchID(pName);
+            if (patchI == -1)
+            {
+                WarningIn("regionSizeDistribution::write()")
+                    << "Cannot find patch " << pName << ". Valid patches are "
+                    << mesh.boundaryMesh().names()
+                    << endl;
+            }
+            else
+            {
+                nPatchFaces += mesh.boundaryMesh()[patchI].size();
+            }
+        }
+
+        Map<label> keepRegions(nPatchFaces);
+        forAll(patchNames_, i)
+        {
+            const word& pName = patchNames_[i];
+
+            label patchI = mesh.boundaryMesh().findPatchID(pName);
+            if (patchI != -1)
+            {
+                const polyPatch& pp = mesh.boundaryMesh()[patchI];
+
+                // Collect all regions on the patch
+                const labelList& faceCells = pp.faceCells();
+
+                forAll(faceCells, i)
+                {
+                    keepRegions.insert
+                    (
+                        regions[faceCells[i]],
+                        Pstream::myProcNo()
+                    );
+                }
+            }
+        }
+
+
+        // Make sure all the processors have the same set of regions
+        Pstream::mapCombineGather(keepRegions, minEqOp<label>());
+        Pstream::mapCombineScatter(keepRegions);
+
+        Info<< "Patch connected regions (liquid core):" << endl;
+        forAllConstIter(Map<label>, keepRegions, iter)
+        {
+            label regionI = iter.key();
+            Pair<scalar>& vols = regionVolume[regionI];
+            Info<< token::TAB << iter.key()
+                << token::TAB << vols.first()
+                << token::TAB << vols.second() << endl;
+
+        }
+        Info<< endl;
+
+        Info<< "Background regions:" << endl;
+        forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+        {
+            if
+            (
+               !keepRegions.found(iter.key())
+             && iter().first() >= volFraction_*meshVol
+            )
+            {
+                Info<< token::TAB << iter.key()
+                    << token::TAB << iter().first()
+                    << token::TAB << iter().second() << endl;
+            }
+        }
+        Info<< endl;
+
+
+        // Split alpha field
+        // ~~~~~~~~~~~~~~~~~
+        // Split into
+        //  - liquidCore            : region connected to inlet patches
+        //  - per region a volume   : for all other regions
+        //  - backgroundAlpha       : remaining alpha
+
+
+        // Construct field
+        volScalarField liquidCore
+        (
+            IOobject
+            (
+                alphaName_ + "_liquidCore",
+                obr_.time().timeName(),
+                obr_,
+                IOobject::NO_READ
+            ),
+            alpha,
+            fvPatchField<scalar>::calculatedType()
+        );
+
+        volScalarField backgroundAlpha
+        (
+            IOobject
+            (
+                alphaName_ + "_background",
+                obr_.time().timeName(),
+                obr_,
+                IOobject::NO_READ
+            ),
+            alpha,
+            fvPatchField<scalar>::calculatedType()
+        );
+
+
+        // Knock out any cell not in keepRegions
+        forAll(liquidCore, cellI)
+        {
+            label regionI = regions[cellI];
+            if (keepRegions.found(regionI))
+            {
+                backgroundAlpha[cellI] = 0;
+            }
+            else
+            {
+                liquidCore[cellI] = 0;
+
+                scalar regionVol = regionVolume[regionI].first();
+                if (regionVol < volFraction_*meshVol)
+                {
+                    backgroundAlpha[cellI] = 0;
+                }
+            }
+        }
+        liquidCore.correctBoundaryConditions();
+        backgroundAlpha.correctBoundaryConditions();
+
+        Info<< "Volume of liquid-core = "
+            << fvc::domainIntegrate(liquidCore).value()
+            << endl;
+
+        Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
+        liquidCore.write();
+
+        Info<< "Volume of background = "
+            << fvc::domainIntegrate(backgroundAlpha).value()
+            << endl;
+
+        Info<< "Writing background field to " << backgroundAlpha.name() << endl;
+        backgroundAlpha.write();
+
+
+
+        // Collect histogram
+        if (Pstream::master())
+        {
+            DynamicList<scalar> diameters(regionVolume.size());
+            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            {
+                if (!keepRegions.found(iter.key()))
+                {
+                    if (iter().first() < volFraction_*meshVol)
+                    {
+                        scalar v = iter().second();
+                      //scalar diam = Foam::cbrt(v*6/mathematicalConstant::pi);
+                        scalar diam =
+                            Foam::cbrt(v*6/constant::mathematical::pi);
+                        diameters.append(diam);
+                    }
+                }
+            }
+
+            if (diameters.size())
+            {
+                scalar maxDiam = max(diameters);
+                scalar minDiam = 0.0;
+
+                Info<< "Maximum diameter:" << maxDiam << endl;
+
+                Histogram<List<scalar> > bins
+                (
+                    minDiam,
+                    maxDiam,
+                    nBins_,
+                    diameters
+                );
+
+                /* 1.7.x
+                scalarField xBin(nBins_);
+
+                scalar dx = (maxDiam-minDiam)/nBins_;
+                scalar x = 0.5*dx;
+                forAll(bins.counts(), i)
+                {
+                    xBin[i] = x;
+                    x += dx;
+                }
+
+                scalarField normalisedCount(bins.counts().size());
+                forAll(bins.counts(), i)
+                {
+                    normalisedCount[i] = 1.0*bins.counts()[i];
+                }
+
+                const coordSet coords
+                (
+                    "diameter",
+                    "x",
+                    xBin
+                );
+                */
+
+                pointField xBin(nBins_);
+                scalar dx = (maxDiam - minDiam)/nBins_;
+                scalar x = 0.5*dx;
+                forAll(bins.counts(), i)
+                {
+                    xBin[i] = point(x, 0, 0);
+                    x += dx;
+                }
+
+                scalarField normalisedCount(bins.counts().size());
+                forAll(bins.counts(), i)
+                {
+                    normalisedCount[i] = 1.0*bins.counts()[i];
+                }
+
+                const coordSet coords
+                (
+                    "diameter",
+                    "x",
+                    xBin,
+                    mag(xBin)
+                );
+                const wordList valNames(1, "count");
+
+
+                fileName outputPath;
+                if (Pstream::parRun())
+                {
+                    outputPath = mesh.time().path()/".."/name_;
+                }
+                else
+                {
+                    outputPath = mesh.time().path()/name_;
+                }
+
+                if (mesh.name() != fvMesh::defaultRegion)
+                {
+                    outputPath = outputPath/mesh.name();
+                }
+
+                mkDir(outputPath/mesh.time().timeName());
+                OFstream str
+                (
+                    outputPath
+                  / mesh.time().timeName()
+                  / formatterPtr_().getFileName(coords, valNames)
+                );
+                Info<< "Writing distribution to " << str.name() << endl;
+
+                List<const scalarField*> valPtrs(1);
+                valPtrs[0] = &normalisedCount;
+                formatterPtr_().write(coords, valNames, valPtrs, str);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
new file mode 100644
index 0000000000..033b5f9aec
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionSizeDistribution
+
+Description
+    Looks up a field, interpolates it to the faces and determines a connected
+    region from a patch where the field is above a certain value.
+    - Writes a field containing all regions starting at given patch
+      ('liquid core')
+    - All other regions are summed for volume and a histogram is calculated.
+
+SourceFiles
+    regionSizeDistribution.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionSizeDistribution_H
+#define regionSizeDistribution_H
+
+#include "pointFieldFwd.H"
+#include "writer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class regionSizeDistribution Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionSizeDistribution
+{
+    // Private data
+
+        //- Name of this set of regionSizeDistribution objects
+        word name_;
+
+        const objectRegistry& obr_;
+
+        //- on/off switch
+        bool active_;
+
+        //- Name of field
+        word alphaName_;
+
+        //- Patches to walk from
+        wordList patchNames_;
+
+        //- Clip value
+        scalar threshold_;
+
+        //- Background region volFraction
+        scalar volFraction_;
+
+        //- Mumber of bins
+        label nBins_;
+
+        //- Output formatter to write
+        autoPtr<writer<scalar> > formatterPtr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        regionSizeDistribution(const regionSizeDistribution&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const regionSizeDistribution&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("regionSizeDistribution");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        regionSizeDistribution
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    // Destructor
+
+        virtual ~regionSizeDistribution();
+
+
+    // Member Functions
+
+        //- Return name of the set of regionSizeDistribution
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the regionSizeDistribution data
+        virtual void read(const dictionary&);
+
+        //- Execute, currently does nothing
+        virtual void execute();
+
+        //- Execute at the final time-loop, currently does nothing
+        virtual void end();
+
+        //- Calculate the regionSizeDistribution and write
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&)
+        {}
+
+        //- Update for changes of mesh
+        virtual void movePoints(const pointField&)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.C
new file mode 100644
index 0000000000..21e25fd6b2
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.C
@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSizeDistributionFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug
+    (
+        regionSizeDistributionFunctionObject,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        regionSizeDistributionFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.H b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.H
new file mode 100644
index 0000000000..97526d0c3a
--- /dev/null
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionFunctionObject.H
@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::regionSizeDistributionFunctionObject
+
+Description
+    FunctionObject wrapper around regionSizeDistribution to allow it to be
+    created via the functions list within controlDict.
+
+SourceFiles
+    regionSizeDistributionFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionSizeDistributionFunctionObject_H
+#define regionSizeDistributionFunctionObject_H
+
+#include "regionSizeDistribution.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<regionSizeDistribution>
+        regionSizeDistributionFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/incompressible/LES/Make/files b/src/turbulenceModels/incompressible/LES/Make/files
index ff2bde02f2..eadcf1a99f 100644
--- a/src/turbulenceModels/incompressible/LES/Make/files
+++ b/src/turbulenceModels/incompressible/LES/Make/files
@@ -34,5 +34,4 @@ wallFunctions=derivedFvPatchFields/wallFunctions
 nuSgsWallFunctions=$(wallFunctions)/nuSgsWallFunctions
 $(nuSgsWallFunctions)/nuSgsUSpaldingWallFunction/nuSgsUSpaldingWallFunctionFvPatchScalarField.C
 
-
 LIB = $(FOAM_LIBBIN)/libincompressibleLESModels