Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

etc/bashrc

@@ -137,7 +137,7 @@ _foamSource()
 {
     while [ $# -ge 1 ]
     do
-        [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "Sourcing: $1"
+        [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "Sourcing: $1" 1>&2
         . $1
        shift
     done
@@ -155,12 +155,12 @@ _foamEval()
             ;;
         *=)
             # name=       -> unset name
-            [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "unset ${1%=}"
+            [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "unset ${1%=}" 1>&2
             eval "unset ${1%=}"
             ;;
         *=*)
             # name=value  -> export name=value
-            [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "export $1"
+            [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "export $1" 1>&2
             eval "export $1"
             ;;
         *)

etc/config/aliases.sh

@@ -53,7 +53,7 @@ unset foamPV
 foamPV()
 {
     . $WM_PROJECT_DIR/etc/config/paraview.sh ParaView_VERSION=$1
-    echo "paraview-$ParaView_VERSION  (major: $ParaView_MAJOR)"
+    echo "paraview-$ParaView_VERSION  (major: $ParaView_MAJOR)" 1>&2
 }
 
 

etc/config/settings.sh

@@ -93,7 +93,7 @@ Linux)
             export WM_LDFLAGS='-m64'
             ;;
         *)
-            echo "Unknown WM_ARCH_OPTION '$WM_ARCH_OPTION', should be 32 or 64"
+            echo "Unknown WM_ARCH_OPTION '$WM_ARCH_OPTION', should be 32 or 64" 1>&2
             ;;
         esac
         ;;
@@ -135,7 +135,7 @@ Linux)
         ;;
 
     *)
-        echo Unknown processor type `uname -m` for Linux
+        echo Unknown processor type `uname -m` for Linux 1>&2
         ;;
     esac
     ;;
@@ -152,7 +152,7 @@ SunOS)
     ;;
 
 *)    # an unsupported operating system
-    cat <<USAGE
+    /bin/cat <<USAGE 1>&2
 
     Your "$WM_ARCH" operating system is not supported by this release
     of OpenFOAM. For further assistance, please contact www.OpenFOAM.org
@@ -268,11 +268,11 @@ OpenFOAM | ThirdParty)
         #clang_version=llvm-svn
         ;;
     *)
-        echo
+        echo 1>&2
-        echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:"
+        echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:" 1>&2
-        echo "    Unknown OpenFOAM compiler type '$WM_COMPILER'"
+        echo "    Unknown OpenFOAM compiler type '$WM_COMPILER'" 1>&2
-        echo "    Please check your settings"
+        echo "    Please check your settings" 1>&2
-        echo
+        echo 1>&2
         ;;
     esac
 
@@ -288,11 +288,11 @@ OpenFOAM | ThirdParty)
 
         # Check that the compiler directory can be found
         [ -d "$gccDir" ] || {
-            echo
+            echo 1>&2
-            echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:"
+            echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:" 1>&2
-            echo "    Cannot find $gccDir installation."
+            echo "    Cannot find $gccDir installation." 1>&2
-            echo "    Please install this compiler version or if you wish to use the system compiler,"
+            echo "    Please install this compiler version or if you wish to use the system compiler," 1>&2
-            echo "    change the 'foamCompiler' setting to 'system'"
+            echo "    change the 'foamCompiler' setting to 'system'" 1>&2
             echo
         }
 
@@ -325,12 +325,12 @@ OpenFOAM | ThirdParty)
 
         # Check that the compiler directory can be found
         [ -d "$clangDir" ] || {
-            echo
+            echo 1>&2
-            echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:"
+            echo "Warning in $WM_PROJECT_DIR/etc/config/settings.sh:" 1>&2
-            echo "    Cannot find $clangDir installation."
+            echo "    Cannot find $clangDir installation." 1>&2
-            echo "    Please install this compiler version or if you wish to use the system compiler,"
+            echo "    Please install this compiler version or if you wish to use the system compiler," 1>&2
-            echo "    change the 'foamCompiler' setting to 'system'"
+            echo "    change the 'foamCompiler' setting to 'system'" 1>&2
-            echo
+            echo 1>&2
         }
 
         _foamAddMan     $clangDir/share/man
@@ -415,10 +415,10 @@ SYSTEMOPENMPI)
 
     if [ "$FOAM_VERBOSE" -a "$PS1" ]
     then
-        echo "Using system installed MPI:"
+        echo "Using system installed MPI:" 1>&2
-        echo "    compile flags : $PINC"
+        echo "    compile flags : $PINC" 1>&2
-        echo "    link flags    : $PLIBS"
+        echo "    link flags    : $PLIBS" 1>&2
-        echo "    libmpi dir    : $libDir"
+        echo "    libmpi dir    : $libDir" 1>&2
     fi
 
     _foamAddLib     $libDir
@@ -492,7 +492,7 @@ HPMPI)
         _foamAddLib $MPI_ARCH_PATH/lib/linux_ia64
         ;;
     *)
-        echo Unknown processor type `uname -m` for Linux
+        echo Unknown processor type `uname -m` 1>&2
         ;;
     esac
     ;;
@@ -526,13 +526,10 @@ QSMPI)
     ;;
 
 SGIMPI)
-    lastCharID=$(( ${#MPI_ROOT} - 1 ))
+    # no trailing slash
-    if [ "${MPI_ROOT:$lastCharID:1}" == '/' ]
+    [ "${MPI_ROOT%/}" = "${MPI_ROOT}" ] || MPI_ROOT="${MPI_ROOT%/}"
-    then
-        MPI_ROOT=${MPI_ROOT:0:$lastCharID}
-    fi
 
-    export FOAM_MPI=${MPI_ROOT##*/}
+    export FOAM_MPI="${MPI_ROOT##*/}"
     export MPI_ARCH_PATH=$MPI_ROOT
 
     if [ ! -d "$MPI_ROOT" -o -z "$MPI_ARCH_PATH" ]
@@ -545,9 +542,9 @@ SGIMPI)
 
     if [ "$FOAM_VERBOSE" -a "$PS1" ]
     then
-        echo "Using SGI MPT:"
+        echo "Using SGI MPT:" 1>&2
-        echo "    MPI_ROOT : $MPI_ROOT"
+        echo "    MPI_ROOT : $MPI_ROOT" 1>&2
-        echo "    FOAM_MPI : $FOAM_MPI"
+        echo "    FOAM_MPI : $FOAM_MPI" 1>&2
     fi
 
     _foamAddPath    $MPI_ARCH_PATH/bin
@@ -555,13 +552,10 @@ SGIMPI)
     ;;
 
 INTELMPI)
-    lastCharID=$(( ${#MPI_ROOT} - 1 ))
+    # no trailing slash
-    if [ "${MPI_ROOT:$lastCharID:1}" == '/' ]
+    [ "${MPI_ROOT%/}" = "${MPI_ROOT}" ] || MPI_ROOT="${MPI_ROOT%/}"
-    then
-        MPI_ROOT=${MPI_ROOT:0:$lastCharID}
-    fi
 
-    export FOAM_MPI=${MPI_ROOT##*/}
+    export FOAM_MPI="${MPI_ROOT##*/}"
     export MPI_ARCH_PATH=$MPI_ROOT
 
     if [ ! -d "$MPI_ROOT" -o -z "$MPI_ARCH_PATH" ]
@@ -574,9 +568,9 @@ INTELMPI)
 
     if [ "$FOAM_VERBOSE" -a "$PS1" ]
     then
-        echo "Using INTEL MPI:"
+        echo "Using INTEL MPI:" 1>&2
-        echo "    MPI_ROOT : $MPI_ROOT"
+        echo "    MPI_ROOT : $MPI_ROOT" 1>&2
-        echo "    FOAM_MPI : $FOAM_MPI"
+        echo "    FOAM_MPI : $FOAM_MPI" 1>&2
     fi
 
     _foamAddPath    $MPI_ARCH_PATH/bin64

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C

@@ -369,7 +369,7 @@ void Foam::rotorDiskSource::constructGeometry()
 
             // blade flap angle [radians]
             scalar beta =
-                flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi);
+                flap_.beta0 - flap_.beta1c*cos(psi) - flap_.beta2s*sin(psi);
 
             // determine rotation tensor to convert from planar system into the
             // rotor cone system
@@ -659,11 +659,11 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
 
         const dictionary& flapCoeffs(coeffs_.subDict("flapCoeffs"));
         flapCoeffs.lookup("beta0") >> flap_.beta0;
-        flapCoeffs.lookup("beta1") >> flap_.beta1;
+        flapCoeffs.lookup("beta1c") >> flap_.beta1c;
-        flapCoeffs.lookup("beta2") >> flap_.beta2;
+        flapCoeffs.lookup("beta2s") >> flap_.beta2s;
         flap_.beta0 = degToRad(flap_.beta0);
-        flap_.beta1 = degToRad(flap_.beta1);
+        flap_.beta1c = degToRad(flap_.beta1c);
-        flap_.beta2 = degToRad(flap_.beta2);
+        flap_.beta2s = degToRad(flap_.beta2s);
 
 
         // create co-ordinate system

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H

@@ -51,8 +51,8 @@ Description
             flapCoeffs
             {
                 beta0           0;  // coning angle [deg]
-                beta1           0;  // lateral flapping coeff
+                beta1c          0;  // lateral flapping coeff (cos coeff)
-                beta2           0;  // longitudinal flapping coeff
+                beta2s          0;  // longitudinal flapping coeff (sin coeff)
             }
 
             blade
@@ -78,7 +78,6 @@ Description
             local         : use local flow conditions
 
 
-
 SourceFiles
     rotorDiskSource.C
     rotorDiskSourceTemplates.C
@@ -138,8 +137,8 @@ protected:
         struct flapData
         {
             scalar beta0;   // coning angle
-            scalar beta1;   // lateral flapping coeff
+            scalar beta1c;  // lateral flapping coeff (cos coeff)
-            scalar beta2;   // longitudinal flapping coeff
+            scalar beta2s;  // longitudinal flapping coeff (sin coeff)
         };
 
 
@@ -148,7 +147,7 @@ protected:
         //- Name of density field
         word rhoName_;
 
-        //- Reference density if rhoName = 'none'
+        //- Reference density for rhoName = 'none'
         scalar rhoRef_;
 
         //- Rotational speed [rad/s]
@@ -254,6 +253,9 @@ public:
 
         // Access
 
+            //- Return the reference density for rhoName = 'none'
+            inline scalar rhoRef() const;
+
             //- Return the rotational speed [rad/s]
             //  Positive anti-clockwise when looking along -ve lift direction
             inline scalar omega() const;

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H

@@ -27,11 +27,18 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+Foam::scalar Foam::rotorDiskSource::rhoRef() const
+{
+    return rhoRef_;
+}
+
+
 Foam::scalar Foam::rotorDiskSource::omega() const
 {
     return omega_;
 }
 
+
 const Foam::List<Foam::point>& Foam::rotorDiskSource::x() const
 {
     return x_;

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C

@@ -70,25 +70,33 @@ Foam::vector Foam::targetCoeffTrim::calcCoeffs
     {
         label cellI = cells[i];
 
-        // create normalisation coefficient
-        scalar radius = x[i].x();
-        scalar coeff2 = coeff1*pow4(radius);
-        if (compressible)
-        {
-            coeff2 *= trho()[cellI];
-        }
-
         vector fc = force[cellI];
         vector mc = fc^(C[cellI] - origin);
 
-        // add to coefficient vector
+        if (useCoeffs_)
-        cf[0] += (fc & yawAxis)/(coeff2 + ROOTVSMALL);
+        {
-        cf[1] += (mc & pitchAxis)/(coeff2*radius + ROOTVSMALL);
+            scalar radius = x[i].x();
-        cf[2] += (mc & rollAxis)/(coeff2*radius + ROOTVSMALL);
+            scalar coeff2 = coeff1*pow4(radius);
+            if (compressible)
+            {
+                coeff2 *= trho()[cellI];
+            }
+
+            // add to coefficient vector
+            cf[0] += (fc & yawAxis)/(coeff2 + ROOTVSMALL);
+            cf[1] += (mc & pitchAxis)/(coeff2*radius + ROOTVSMALL);
+            cf[2] += (mc & rollAxis)/(coeff2*radius + ROOTVSMALL);
+        }
+        else
+        {
+            cf[0] += fc & yawAxis;
+            cf[1] += mc & pitchAxis;
+            cf[2] += mc & rollAxis;
+        }
     }
 
     reduce(cf, sumOp<vector>());
-    
+
     return cf;
 }
 
@@ -103,6 +111,7 @@ Foam::targetCoeffTrim::targetCoeffTrim
 :
     trimModel(rotor, dict, typeName),
     calcFrequency_(-1),
+    useCoeffs_(true),
     target_(vector::zero),
     theta_(vector::zero),
     nIter_(50),
@@ -128,9 +137,16 @@ void Foam::targetCoeffTrim::read(const dictionary& dict)
     trimModel::read(dict);
 
     const dictionary& targetDict(coeffs_.subDict("target"));
-    target_[0] = readScalar(targetDict.lookup("thrustCoeff"));
+    useCoeffs_ = targetDict.lookupOrDefault<bool>("useCoeffs", true);
-    target_[1] = readScalar(targetDict.lookup("pitchCoeff"));
+    word ext = "";
-    target_[2] = readScalar(targetDict.lookup("rollCoeff"));
+    if (useCoeffs_)
+    {
+        ext = "Coeff";
+    }
+
+    target_[0] = readScalar(targetDict.lookup("thrust" + ext));
+    target_[1] = readScalar(targetDict.lookup("pitch" + ext));
+    target_[2] = readScalar(targetDict.lookup("roll" + ext));
 
     const dictionary& pitchAngleDict(coeffs_.subDict("pitchAngles"));
     theta_[0] = degToRad(readScalar(pitchAngleDict.lookup("theta0Ini")));
@@ -173,8 +189,16 @@ void Foam::targetCoeffTrim::correct(const vectorField& U, vectorField& force)
 {
     if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
     {
+        word calcType = "forces";
+        if (useCoeffs_)
+        {
+            calcType = "coefficients";
+        }
+
         Info<< type() << ":" << nl
-            << "    solving for target trim coefficients" << nl;
+            << "    solving for target trim " << calcType << nl;
+
+        const scalar rhoRef = rotor_.rhoRef();
 
         // iterate to find new pitch angles to achieve target force
         scalar err = GREAT;
@@ -210,7 +234,7 @@ void Foam::targetCoeffTrim::correct(const vectorField& U, vectorField& force)
             }
 
             // calculate the change in pitch angle vector
-            vector dt = inv(J) & (target_ - old);
+            vector dt = inv(J) & (target_/rhoRef - old);
 
             // update pitch angles
             vector thetaNew = theta_ + relax_*dt;
@@ -235,10 +259,10 @@ void Foam::targetCoeffTrim::correct(const vectorField& U, vectorField& force)
                 << "), iterations = " << iter << endl;
         }
 
-        Info<< "    current and target coefficients:" << nl
+        Info<< "    current and target " << calcType << nl
-            << "        thrust  = " << old[0] << ", " << target_[0] << nl
+            << "        thrust  = " << old[0]*rhoRef << ", " << target_[0] << nl
-            << "        pitch   = " << old[1] << ", " << target_[1] << nl
+            << "        pitch   = " << old[1]*rhoRef << ", " << target_[1] << nl
-            << "        roll    = " << old[2] << ", " << target_[2] << nl
+            << "        roll    = " << old[2]*rhoRef << ", " << target_[2] << nl
             << "    new pitch angles [deg]:" << nl
             << "        theta0  = " << radToDeg(theta_[0]) << nl
             << "        theta1c = " << radToDeg(theta_[1]) << nl

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H

@@ -25,7 +25,7 @@ Class
     Foam::targetCoeffTrim
 
 Description
-    Target trim coefficients
+    Target trim forces/coefficients
 
     Solves:
 
@@ -43,8 +43,8 @@ Description
     Newton-Raphson iterative method.  The solver tolerance can be user-input,
     using the 'tol' entry.
 
-
+    If coefficients are requested (useCoeffs = true), the force and moments
-    The coefficients are determined using:
+    are normalised using:
 
                          force
         c = ---------------------------------
@@ -97,6 +97,9 @@ protected:
         //- Number of iterations between calls to 'correct'
         label calcFrequency_;
 
+        //- Flag to indicate whether to solve coeffs (true) or forces (false)
+        bool useCoeffs_;
+
         //- Target coefficient vector (thrust force, roll moment, pitch moment)
         vector target_;
 

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.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
@@ -246,7 +246,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
     const scalar mass = p.mass();
     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
-    scalar tMom = 1.0/(Fcp.Sp() + Fncp.Sp());
+    scalar tMom = mass/(Fcp.Sp() + Fncp.Sp());
 
     const vector g = td.cloud().g().value();
 
@@ -301,7 +301,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
         child->tc() = -GREAT;
         child->ms() = 0.0;
         child->injector() = this->injector();
-        child->tMom() = 1.0/(Fcp.Sp() + Fncp.Sp());
+        child->tMom() = massChild/(Fcp.Sp() + Fncp.Sp());
         child->user() = 0.0;
         child->setCellValues(td, dt, cellI);
 

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

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);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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
+    );
+}
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/LES/Make/files

@@ -34,5 +34,4 @@ wallFunctions=derivedFvPatchFields/wallFunctions
 nuSgsWallFunctions=$(wallFunctions)/nuSgsWallFunctions
 $(nuSgsWallFunctions)/nuSgsUSpaldingWallFunction/nuSgsUSpaldingWallFunctionFvPatchScalarField.C
 
-
 LIB = $(FOAM_LIBBIN)/libincompressibleLESModels