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 .
+
+\*---------------------------------------------------------------------------*/
+
+#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
+ template
+ class ListPlusEqOp
+ {
+ public:
+ void operator()
+ (
+ FixedList& x,
+ const FixedList& 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(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::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(obr_);
+
+ autoPtr alphaPtr;
+ if (obr_.foundObject(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(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 townFld(fvp.patchInternalField());
+ const scalarField& ownFld = townFld();
+ tmp 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 > regionVolume(regions.nRegions()/Pstream::nProcs());
+ forAll(alpha, cellI)
+ {
+ scalar cellVol = mesh.V()[cellI];
+ scalar alphaVol = alpha[cellI]*cellVol;
+
+ label regionI = regions[cellI];
+
+ Map >::iterator fnd = regionVolume.find(regionI);
+ if (fnd == regionVolume.end())
+ {
+ regionVolume.insert
+ (
+ regionI,
+ Pair(cellVol, alphaVol)
+ );
+ }
+ else
+ {
+ fnd().first() += cellVol;
+ fnd().second() += alphaVol;
+ }
+ }
+ Pstream::mapCombineGather(regionVolume, ListPlusEqOp());
+ 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 >, 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