Merge branch 'master' into clouds

2025-11-28 03:28:01 +00:00 · 2010-10-25 17:53:29 +01:00
parent 1a61449dce cfbc1fcbef
commit 17c15af32b
57 changed files with 4307 additions and 223 deletions
--- a/applications/solvers/multiphase/interFoam/Allwclean
+++ b/applications/solvers/multiphase/interFoam/Allwclean
@ -6,5 +6,6 @@ wclean
 wclean interDyMFoam
 wclean MRFInterFoam
 wclean porousInterFoam
+wclean LTSInterFoam

 # ----------------------------------------------------------------- end-of-file
--- a/applications/solvers/multiphase/interFoam/Allwmake
+++ b/applications/solvers/multiphase/interFoam/Allwmake
@ -6,5 +6,6 @@ wmake
 wmake interDyMFoam
 wmake MRFInterFoam
 wmake porousInterFoam
+wmake LTSInterFoam

 # ----------------------------------------------------------------- end-of-file
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     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/>.
+
+Application
+    interFoam
+
+Description
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+
+    For a two-fluid approach see twoPhaseEulerFoam.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "interfaceProperties.H"
+#include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialrDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "setrDeltaT.H"
+
+        twoPhaseProperties.correct();
+
+        #include "alphaEqnSubCycle.H"
+        turbulence->correct();
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     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 "MULES.H"
+#include "upwind.H"
+#include "uncorrectedSnGrad.H"
+#include "gaussConvectionScheme.H"
+#include "gaussLaplacianScheme.H"
+#include "uncorrectedSnGrad.H"
+#include "surfaceInterpolate.H"
+#include "fvcSurfaceIntegrate.H"
+#include "slicedSurfaceFields.H"
+#include "syncTools.H"
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::MULES::explicitLTSSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    explicitLTSSolve
+    (
+        geometricOneField(),
+        psi,
+        phi,
+        phiPsi,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
+void Foam::MULES::implicitSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    implicitSolve
+    (
+        geometricOneField(),
+        psi,
+        phi,
+        phiPsi,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     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/>.
+
+Global
+    MULES
+
+Description
+    MULES: Multidimensional universal limiter with explicit solution.
+
+    Solve a convective-only transport equation using an explicit universal
+    multi-dimensional limiter.
+
+    Parameters are the variable to solve, the normal convective flux and the
+    actual explicit flux of the variable which is also used to return limited
+    flux used in the bounded-solution.
+
+SourceFiles
+    MULES.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MULES_H
+#define MULES_H
+
+#include "volFields.H"
+#include "surfaceFieldsFwd.H"
+#include "primitiveFieldsFwd.H"
+#include "zeroField.H"
+#include "geometricOneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace MULES
+{
+
+template<class RhoType, class SpType, class SuType>
+void explicitLTSSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void explicitLTSSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+template<class RhoType, class SpType, class SuType>
+void implicitSolve
+(
+    const RhoType& rho,
+    volScalarField& gamma,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void implicitSolve
+(
+    volScalarField& gamma,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+template<class RhoType, class SpType, class SuType>
+void limiter
+(
+    scalarField& allLambda,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+);
+
+} // End namespace MULES
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MULESTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
@ -0,0 +1,602 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     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 "MULES.H"
+#include "upwind.H"
+#include "uncorrectedSnGrad.H"
+#include "gaussConvectionScheme.H"
+#include "gaussLaplacianScheme.H"
+#include "uncorrectedSnGrad.H"
+#include "surfaceInterpolate.H"
+#include "fvcSurfaceIntegrate.H"
+#include "slicedSurfaceFields.H"
+#include "syncTools.H"
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::explicitLTSSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    Info<< "MULES: Solving for " << psi.name() << endl;
+
+    const fvMesh& mesh = psi.mesh();
+    psi.correctBoundaryConditions();
+
+    surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi);
+
+    surfaceScalarField& phiCorr = phiPsi;
+    phiCorr -= phiBD;
+
+    scalarField allLambda(mesh.nFaces(), 1.0);
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    limiter
+    (
+        allLambda,
+        rho,
+        psi,
+        phiBD,
+        phiCorr,
+        Sp.field(),
+        Su.field(),
+        psiMax,
+        psiMin,
+        3
+    );
+
+    phiPsi = phiBD + lambda*phiCorr;
+
+    scalarField& psiIf = psi;
+    const scalarField& psi0 = psi.oldTime();
+
+    const volScalarField& rDeltaT =
+        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+
+    psiIf = 0.0;
+    fvc::surfaceIntegrate(psiIf, phiPsi);
+
+    if (mesh.moving())
+    {
+        psiIf =
+        (
+            mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
+          + Su.field()
+          - psiIf
+        )/(rho*rDeltaT - Sp.field());
+    }
+    else
+    {
+        psiIf =
+        (
+            rho.oldTime()*psi0*rDeltaT
+          + Su.field()
+          - psiIf
+        )/(rho*rDeltaT - Sp.field());
+    }
+
+    psi.correctBoundaryConditions();
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::implicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    label maxIter
+    (
+        readLabel(MULEScontrols.lookup("maxIter"))
+    );
+
+    label nLimiterIter
+    (
+        readLabel(MULEScontrols.lookup("nLimiterIter"))
+    );
+
+    scalar maxUnboundedness
+    (
+        readScalar(MULEScontrols.lookup("maxUnboundedness"))
+    );
+
+    scalar CoCoeff
+    (
+        readScalar(MULEScontrols.lookup("CoCoeff"))
+    );
+
+    scalarField allCoLambda(mesh.nFaces());
+
+    {
+        surfaceScalarField Cof =
+            mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
+           *mag(phi)/mesh.magSf();
+
+        slicedSurfaceScalarField CoLambda
+        (
+            IOobject
+            (
+                "CoLambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimless,
+            allCoLambda,
+            false   // Use slices for the couples
+        );
+
+        CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
+    }
+
+    scalarField allLambda(allCoLambda);
+    //scalarField allLambda(mesh.nFaces(), 1.0);
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    linear<scalar> CDs(mesh);
+    upwind<scalar> UDs(mesh, phi);
+    //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
+
+    fvScalarMatrix psiConvectionDiffusion
+    (
+        fvm::ddt(rho, psi)
+      + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
+        //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
+        //.fvmLaplacian(Dpsif, psi)
+      - fvm::Sp(Sp, psi)
+      - Su
+    );
+
+    surfaceScalarField phiBD = psiConvectionDiffusion.flux();
+
+    surfaceScalarField& phiCorr = phiPsi;
+    phiCorr -= phiBD;
+
+    for (label i=0; i<maxIter; i++)
+    {
+        if (i != 0 && i < 4)
+        {
+            allLambda = allCoLambda;
+        }
+
+        limiter
+        (
+            allLambda,
+            rho,
+            psi,
+            phiBD,
+            phiCorr,
+            Sp.field(),
+            Su.field(),
+            psiMax,
+            psiMin,
+            nLimiterIter
+        );
+
+        solve
+        (
+            psiConvectionDiffusion + fvc::div(lambda*phiCorr),
+            MULEScontrols
+        );
+
+        scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
+        scalar minPsi = gMin(psi.internalField());
+
+        scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
+
+        if (unboundedness < maxUnboundedness)
+        {
+            break;
+        }
+        else
+        {
+            Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
+                << " min(" << psi.name() << ") = " << minPsi << endl;
+
+            phiBD = psiConvectionDiffusion.flux();
+
+            /*
+            word gammaScheme("div(phi,gamma)");
+            word gammarScheme("div(phirb,gamma)");
+
+            const surfaceScalarField& phir =
+                mesh.lookupObject<surfaceScalarField>("phir");
+
+            phiCorr =
+                fvc::flux
+                (
+                    phi,
+                    psi,
+                    gammaScheme
+                )
+              + fvc::flux
+                (
+                    -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
+                    psi,
+                    gammarScheme
+                )
+                - phiBD;
+            */
+        }
+    }
+
+    phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limiter
+(
+    scalarField& allLambda,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+)
+{
+    const scalarField& psiIf = psi;
+    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
+
+    const scalarField& psi0 = psi.oldTime();
+
+    const fvMesh& mesh = psi.mesh();
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighb = mesh.neighbour();
+    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+    const scalarField& V = tVsc();
+
+    const volScalarField& rDeltaT =
+        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+
+    const scalarField& phiBDIf = phiBD;
+    const surfaceScalarField::GeometricBoundaryField& phiBDBf =
+        phiBD.boundaryField();
+
+    const scalarField& phiCorrIf = phiCorr;
+    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
+        phiCorr.boundaryField();
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    scalarField& lambdaIf = lambda;
+    surfaceScalarField::GeometricBoundaryField& lambdaBf =
+        lambda.boundaryField();
+
+    scalarField psiMaxn(psiIf.size(), psiMin);
+    scalarField psiMinn(psiIf.size(), psiMax);
+
+    scalarField sumPhiBD(psiIf.size(), 0.0);
+
+    scalarField sumPhip(psiIf.size(), VSMALL);
+    scalarField mSumPhim(psiIf.size(), VSMALL);
+
+    forAll(phiCorrIf, facei)
+    {
+        label own = owner[facei];
+        label nei = neighb[facei];
+
+        psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
+        psiMinn[own] = min(psiMinn[own], psiIf[nei]);
+
+        psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
+        psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
+
+        sumPhiBD[own] += phiBDIf[facei];
+        sumPhiBD[nei] -= phiBDIf[facei];
+
+        scalar phiCorrf = phiCorrIf[facei];
+
+        if (phiCorrf > 0.0)
+        {
+            sumPhip[own] += phiCorrf;
+            mSumPhim[nei] += phiCorrf;
+        }
+        else
+        {
+            mSumPhim[own] -= phiCorrf;
+            sumPhip[nei] -= phiCorrf;
+        }
+    }
+
+    forAll(phiCorrBf, patchi)
+    {
+        const fvPatchScalarField& psiPf = psiBf[patchi];
+        const scalarField& phiBDPf = phiBDBf[patchi];
+        const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+        const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+        if (psiPf.coupled())
+        {
+            scalarField psiPNf = psiPf.patchNeighbourField();
+
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
+            }
+        }
+        else
+        {
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
+            }
+        }
+
+        forAll(phiCorrPf, pFacei)
+        {
+            label pfCelli = pFaceCells[pFacei];
+
+            sumPhiBD[pfCelli] += phiBDPf[pFacei];
+
+            scalar phiCorrf = phiCorrPf[pFacei];
+
+            if (phiCorrf > 0.0)
+            {
+                sumPhip[pfCelli] += phiCorrf;
+            }
+            else
+            {
+                mSumPhim[pfCelli] -= phiCorrf;
+            }
+        }
+    }
+
+    psiMaxn = min(psiMaxn, psiMax);
+    psiMinn = max(psiMinn, psiMin);
+
+    //scalar smooth = 0.5;
+    //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+    //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+
+    if (mesh.moving())
+    {
+        tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
+
+        psiMaxn =
+            V*((rho*rDeltaT - Sp)*psiMaxn - Su)
+          - (V0()*rDeltaT)*rho.oldTime()*psi0
+          + sumPhiBD;
+
+        psiMinn =
+            V*(Su - (rho*rDeltaT - Sp)*psiMinn)
+          + (V0*rDeltaT)*rho.oldTime()*psi0
+          - sumPhiBD;
+    }
+    else
+    {
+        psiMaxn =
+            V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su)
+          + sumPhiBD;
+
+        psiMinn =
+            V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su)
+          - sumPhiBD;
+    }
+
+    scalarField sumlPhip(psiIf.size());
+    scalarField mSumlPhim(psiIf.size());
+
+    for(int j=0; j<nLimiterIter; j++)
+    {
+        sumlPhip = 0.0;
+        mSumlPhim = 0.0;
+
+        forAll(lambdaIf, facei)
+        {
+            label own = owner[facei];
+            label nei = neighb[facei];
+
+            scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+
+            if (lambdaPhiCorrf > 0.0)
+            {
+                sumlPhip[own] += lambdaPhiCorrf;
+                mSumlPhim[nei] += lambdaPhiCorrf;
+            }
+            else
+            {
+                mSumlPhim[own] -= lambdaPhiCorrf;
+                sumlPhip[nei] -= lambdaPhiCorrf;
+            }
+        }
+
+        forAll(lambdaBf, patchi)
+        {
+            scalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+            forAll(lambdaPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
+
+                if (lambdaPhiCorrf > 0.0)
+                {
+                    sumlPhip[pfCelli] += lambdaPhiCorrf;
+                }
+                else
+                {
+                    mSumlPhim[pfCelli] -= lambdaPhiCorrf;
+                }
+            }
+        }
+
+        forAll (sumlPhip, celli)
+        {
+            sumlPhip[celli] =
+                max(min
+                (
+                    (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
+                    1.0), 0.0
+                );
+
+            mSumlPhim[celli] =
+                max(min
+                (
+                    (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
+                    1.0), 0.0
+                );
+        }
+
+        const scalarField& lambdam = sumlPhip;
+        const scalarField& lambdap = mSumlPhim;
+
+        forAll(lambdaIf, facei)
+        {
+            if (phiCorrIf[facei] > 0.0)
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdap[owner[facei]], lambdam[neighb[facei]])
+                );
+            }
+            else
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdam[owner[facei]], lambdap[neighb[facei]])
+                );
+            }
+        }
+
+
+        forAll(lambdaBf, patchi)
+        {
+            fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+            forAll(lambdaPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                if (phiCorrfPf[pFacei] > 0.0)
+                {
+                    lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
+                }
+                else
+                {
+                    lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
+                }
+            }
+        }
+
+        syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
+    }
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
@ -0,0 +1,4 @@
+LTSInterFoam.C
+MULES.C
+
+EXE = $(FOAM_APPBIN)/LTSInterFoam
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
@ -0,0 +1,15 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -ltwoPhaseInterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
@ -0,0 +1,36 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic = mag(phi/mesh.magSf());
+    phic = min(interface.cAlpha()*phic, max(phic));
+    surfaceScalarField phir = phic*interface.nHatf();
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        surfaceScalarField phiAlpha =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
+        //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+
+        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+    }
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+}
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@ -0,0 +1,35 @@
+label nAlphaCorr
+(
+    readLabel(piso.lookup("nAlphaCorr"))
+);
+
+label nAlphaSubCycles
+(
+    readLabel(piso.lookup("nAlphaSubCycles"))
+);
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+#       include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+#       include "alphaEqn.H"
+}
+
+interface.correct();
+
+rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
@ -0,0 +1,30 @@
+scalar maxDeltaT
+(
+    piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+);
+
+volScalarField rDeltaT
+(
+    IOobject
+    (
+        "rDeltaT",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+);
+
+volScalarField rSubDeltaT
+(
+    IOobject
+    (
+        "rSubDeltaT",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+);
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@ -0,0 +1,109 @@
+{
+    scalar maxCo
+    (
+        piso.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    scalar maxAlphaCo
+    (
+        piso.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    label nAlphaSpreadIter
+    (
+        piso.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    label nAlphaSweepIter
+    (
+        piso.lookupOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        piso.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalarField rDeltaT0 = rDeltaT;
+
+    // Set the reciprocal time-step using an effective maximum Courant number
+    rDeltaT = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum
+        (
+            mag(rhoPhi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf())
+        )/rho
+    );
+
+    // Limit the time-step further in the region of the interface
+    {
+        surfaceScalarField alphaf = fvc::interpolate(alpha1);
+
+        surfaceScalarField SfUfbyDelta =
+            pos(alphaf - 0.01)*pos(0.99 - alphaf)
+           *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
+
+        rDeltaT = max
+        (
+            rDeltaT,
+            fvc::surfaceSum(mag(SfUfbyDelta/(maxAlphaCo*mesh.magSf())))
+        );
+    }
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread(rDeltaT, alpha1, nAlphaSpreadIter);
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter);
+    }
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        Info<< "Damping rDeltaT" << endl;
+        rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff);
+    }
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    label nAlphaSubCycles
+    (
+        readLabel(piso.lookup("nAlphaSubCycles"))
+    );
+
+    rSubDeltaT = rDeltaT*nAlphaSubCycles;
+}
--- a/applications/utilities/mesh/advanced/selectCells/selectCells.C
+++ b/applications/utilities/mesh/advanced/selectCells/selectCells.C
@ -55,6 +55,7 @@ Description
 #include "edgeStats.H"
 #include "treeDataTriSurface.H"
 #include "indexedOctree.H"
+#include "globalMeshData.H"

 using namespace Foam;

@ -298,7 +299,7 @@ label selectOutsideCells
        mesh,
        outsideFaces.shrink(),
        outsideFacesInfo.shrink(),
-        mesh.nCells()  // max iterations
+        mesh.globalData().nTotalCells()+1  // max iterations
    );

    // Now regionCalc should hold info on cells that are reachable from
--- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
+++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
@ -43,7 +43,7 @@ Description
 #include "polyMeshZipUpCells.H"
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "Swap.H"
 #include "IFstream.H"
 #include "readHexLabel.H"
@ -900,7 +900,7 @@ int main(int argc, char *argv[])
    fluentToFoamType.insert("interface", polyPatch::typeName);
    fluentToFoamType.insert("internal", polyPatch::typeName);
    fluentToFoamType.insert("solid", polyPatch::typeName);
-    fluentToFoamType.insert("fan", cyclicPolyPatch::typeName);
+    fluentToFoamType.insert("fan", oldCyclicPolyPatch::typeName);
    fluentToFoamType.insert("radiator", polyPatch::typeName);
    fluentToFoamType.insert("porous-jump", polyPatch::typeName);

--- a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C
+++ b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C
@ -41,7 +41,7 @@ Description
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "wedgePolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "unitConversion.H"

 using namespace Foam;
--- a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H
+++ b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H
@ -196,7 +196,7 @@ const word* kivaPatchTypes[nBCs] =
    &polyPatch::typeName,
    &polyPatch::typeName,
    &symmetryPolyPatch::typeName,
-    &cyclicPolyPatch::typeName
+    &oldCyclicPolyPatch::typeName
 };

 enum patchTypeNames
--- a/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C
+++ b/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C
@ -29,7 +29,7 @@ Description
 #include "sammMesh.H"
 #include "Time.H"
 #include "wallPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "preservePatchTypes.H"
 #include "IFstream.H"
@ -208,7 +208,7 @@ void sammMesh::readBoundary()
            {
                // incorrect. should be cyclicPatch but this
                // requires info on connected faces.
-                patchTypes_[patchLabel] = cyclicPolyPatch::typeName;
+                patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName;
            }
            else
            {
--- a/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C
+++ b/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C
@ -29,7 +29,7 @@ Description
 #include "starMesh.H"
 #include "Time.H"
 #include "wallPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "preservePatchTypes.H"
 #include "IFstream.H"
@ -206,7 +206,7 @@ void starMesh::readBoundary()
            {
                // incorrect. should be cyclicPatch but this
                // requires info on connected faces.
-                patchTypes_[patchLabel] = cyclicPolyPatch::typeName;
+                patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName;
            }
            else
            {
--- a/applications/utilities/mesh/manipulation/subsetMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/subsetMesh/Make/options
@ -4,4 +4,5 @@ EXE_INC = \

 EXE_LIBS = \
    -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lgenericPatchFields
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
@ -575,99 +575,99 @@ void ensightPointField
    }


-     label ensightPatchI = eMesh.patchPartOffset();
+    label ensightPatchI = eMesh.patchPartOffset();
 
-     forAll(allPatchNames, patchi)
-     {
-         const word& patchName = allPatchNames[patchi];
+    forAll(allPatchNames, patchi)
+    {
+        const word& patchName = allPatchNames[patchi];
 
-         eMesh.barrier();
+        eMesh.barrier();
 
-         if (patchNames.empty() || patchNames.found(patchName))
-         {
-             const fvPatch& p = mesh.boundary()[patchi];
-             if
-             (
-                 returnReduce(p.size(), sumOp<label>())
-               > 0
-             )
-             {
-                 // Renumber the patch points/faces into unique points
-                 labelList pointToGlobal;
-                 labelList uniqueMeshPointLabels;
-                 autoPtr<globalIndex> globalPointsPtr =
-                 mesh.globalData().mergePoints
-                 (
-                     p.patch().meshPoints(),
-                     p.patch().meshPointMap(),
-                     pointToGlobal,
-                     uniqueMeshPointLabels
-                 );
+        if (patchNames.empty() || patchNames.found(patchName))
+        {
+            const fvPatch& p = mesh.boundary()[patchi];
+            if
+            (
+                returnReduce(p.size(), sumOp<label>())
+              > 0
+            )
+            {
+                // Renumber the patch points/faces into unique points
+                labelList pointToGlobal;
+                labelList uniqueMeshPointLabels;
+                autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    p.patch().meshPoints(),
+                    p.patch().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
 
-                 if (Pstream::master())
-                 {
-                     ensightFile.writePartHeader(ensightPatchI);
-                 }
+                if (Pstream::master())
+                {
+                    ensightFile.writePartHeader(ensightPatchI);
+                }
 
-                 writeField
-                 (
-                     "coordinates",
-                     Field<Type>(pf.internalField(), uniqueMeshPointLabels),
-                     ensightFile
-                 );
+                writeField
+                (
+                    "coordinates",
+                    Field<Type>(pf.internalField(), uniqueMeshPointLabels),
+                    ensightFile
+                );
 
-                 ensightPatchI++;
-             }
-         }
-     }
+                ensightPatchI++;
+            }
+        }
+    }
 
-     // write faceZones, if requested
-     if (faceZoneNames.size())
-     {
-         forAllConstIter(wordHashSet, faceZoneNames, iter)
-         {
-             const word& faceZoneName = iter.key();
- 
-             eMesh.barrier();
- 
-             label zoneID = mesh.faceZones().findZoneID(faceZoneName);
- 
-             const faceZone& fz = mesh.faceZones()[zoneID];
- 
-             if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
-             {
-                 // Renumber the faceZone points/faces into unique points
-                 labelList pointToGlobal;
-                 labelList uniqueMeshPointLabels;
-                 autoPtr<globalIndex> globalPointsPtr =
-                 mesh.globalData().mergePoints
-                 (
-                     fz().meshPoints(),
-                     fz().meshPointMap(),
-                     pointToGlobal,
-                     uniqueMeshPointLabels
-                 );
- 
-                 if (Pstream::master())
-                 {
-                     ensightFile.writePartHeader(ensightPatchI);
-                 }
- 
-                 writeField
-                 (
-                     "coordinates",
-                     Field<Type>
-                     (
-                         pf.internalField(),
-                         uniqueMeshPointLabels
-                     ),
-                     ensightFile
-                 );
- 
-                 ensightPatchI++;
-             }
-         }
-     }
+    // write faceZones, if requested
+    if (faceZoneNames.size())
+    {
+        forAllConstIter(wordHashSet, faceZoneNames, iter)
+        {
+            const word& faceZoneName = iter.key();
+
+            eMesh.barrier();
+
+            label zoneID = mesh.faceZones().findZoneID(faceZoneName);
+
+            const faceZone& fz = mesh.faceZones()[zoneID];
+
+            if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
+            {
+                // Renumber the faceZone points/faces into unique points
+                labelList pointToGlobal;
+                labelList uniqueMeshPointLabels;
+                autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    fz().meshPoints(),
+                    fz().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
+
+                if (Pstream::master())
+                {
+                    ensightFile.writePartHeader(ensightPatchI);
+                }
+
+                writeField
+                (
+                    "coordinates",
+                    Field<Type>
+                    (
+                        pf.internalField(),
+                        uniqueMeshPointLabels
+                    ),
+                    ensightFile
+                );
+
+                ensightPatchI++;
+            }
+        }
+    }

    if (Pstream::master())
    {
--- a/bin/tools/foamConfigurePaths
+++ b/bin/tools/foamConfigurePaths
@ -34,10 +34,10 @@ usage() {
    cat<<USAGE

 usage: ${0##*/}
-  --foamInstall dir       specify installation directory (e.g. /opt)
-  --projectName name      specify project name (e.g. openfoam170)
-  --archOption  arch      specify architecture option (only 32 or 64 applicable)
-  --paraviewInstall dir   specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
+  --foamInstall dir         specify installation directory (e.g. /opt)
+  --projectName name        specify project name (e.g. openfoam170)
+  --archOption  arch        specify architecture option (only 32 or 64 applicable)
+  --paraviewInstall dir     specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)

 * hardcode paths to installation

@ -46,6 +46,17 @@ USAGE
 }


+# Function to do replacement on file. Checks if any replacement has been done.
+# inlineSed <file> <sedCommand> <description>
+_inlineSed()
+{
+    backup=`tempfile`
+    cp $1 $backup
+    sed -i -e "$2" $1
+    cmp $1 $backup || usage "Failed : $3"
+}
+
+
 [ -f etc/bashrc -a -f etc/settings.sh ] || usage "Please run from top-level directory of installation"

 unset foamInstall projectName archOption paraviewInstall
@ -59,30 +70,40 @@ do
        ;;
    --foamInstall)
        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        foamInstall="$2"
-        echo "Replacing foamInstall setting by $foamInstall"
-        sed -i -e '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" etc/bashrc
+    foamInstall="$2"
+        echo "** foamInstall:$foamInstall"
+
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" \
+            "Replacing foamInstall setting by $foamInstall"
        shift 2
        ;;
   --projectName)
        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        projectName="$2"
-        echo "Replacing WM_PROJECT_DIR setting by $projectName"
-        sed -i -e '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" etc/bashrc
+    projectName="$2"
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" \
+            "Replacing WM_PROJECT_DIR setting by $projectName"
        shift 2
        ;;
    --archOption)
        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        archOption="$2"
-        echo "Replacing WM_ARCH_OPTION setting by $archOption"
-        sed -i -e '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" etc/bashrc
+    archOption="$2"
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" \
+             "Replacing WM_ARCH_OPTION setting by $archOption"
        shift 2
        ;;
    --paraviewInstall)
        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        paraviewInstall="$2"
-        echo "Replacing ParaView_DIR setting by $paraviewInstall"
-        sed -i -e '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" etc/apps/paraview3/bashrc
+    paraviewInstall="$2"
+        _inlineSed \
+            etc/apps/paraview3/bashrc \
+            '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" \
+             "Replacing ParaView_DIR setting by $paraviewInstall"
        shift 2
        ;;
    *)
@ -97,11 +118,14 @@ done
 #sed -i -e 's@WM_PROJECT=.*@WM_PROJECT='"$projectName@" etc/bashrc

 # Replace the WM_MPLIB always
-echo "Replacing WM_MPLIB setting by SYSTEMOPENMPI"
-sed -i -e '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' etc/bashrc
-
+_inlineSed \
+    etc/bashrc \
+    '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' \
+    "Replacing WM_MPLIB setting by SYSTEMOPENMPI"
 # Replace the compilerInstall always
-echo "Replacing compilerInstall setting by system"
-sed -i -e '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' etc/settings.sh
+_inlineSed \
+    etc/settings.sh \
+    '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' \
+    "Replacing compilerInstall setting by system"

 #------------------------------------------------------------------------------
--- a/src/OSspecific/POSIX/fileMonitor.C
+++ b/src/OSspecific/POSIX/fileMonitor.C
@ -39,7 +39,7 @@ Class
 #else
 #   include <sys/inotify.h>
 #   include <sys/ioctl.h>
-
+#   include <errno.h>
 #   define EVENT_SIZE  ( sizeof (struct inotify_event) )
 #   define EVENT_LEN   (EVENT_SIZE + 16)
 #   define EVENT_BUF_LEN     ( 1024 * EVENT_LEN )
@ -144,7 +144,7 @@ namespace Foam

 #else
        //- File descriptor for the inotify instance
-        int fd;
+        int inotifyFd_;

        //- Current watchIDs and corresponding directory id
        DynamicList<label> dirWatches_;
@ -153,23 +153,45 @@ namespace Foam
        //- initialise inotify
        inline fileMonitorWatcher(const label sz = 20)
        :
-            fd(inotify_init()),
+            inotifyFd_(inotify_init()),
            dirWatches_(sz),
            dirFiles_(sz)
-        {}
+        {
+            if (inotifyFd_ < 0)
+            {
+                static bool hasWarned = false;
+                if (!hasWarned)
+                {
+                    hasWarned = true;
+                    WarningIn("fileMonitorWatcher(const label)")
+                        << "Failed allocating an inotify descriptor : "
+                        << string(strerror(errno)) << endl
+                        << "    Please increase the number of allowable "
+                        << "inotify instances" << endl
+                        << "    (/proc/sys/fs/inotify/max_user_instances"
+                        << " on Linux)" << endl
+                        << "    or switch off runTimeModifiable." << endl
+                        << "    Continuing without additional file monitoring."
+                        << endl;
+                }
+            }
+        }

        //- remove all watches
        inline ~fileMonitorWatcher()
        {
-            forAll(dirWatches_, i)
+            if (inotifyFd_ >= 0)
            {
-                if (dirWatches_[i] >= 0)
+                forAll(dirWatches_, i)
                {
-                    if (inotify_rm_watch(fd, int(dirWatches_[i])))
+                    if (dirWatches_[i] >= 0)
                    {
-                        WarningIn("fileMonitor::~fileMonitor()")
-                            << "Failed deleting directory watch "
-                            << dirWatches_[i] << endl;
+                        if (inotify_rm_watch(inotifyFd_, int(dirWatches_[i])))
+                        {
+                            WarningIn("fileMonitor::~fileMonitor()")
+                                << "Failed deleting directory watch "
+                                << dirWatches_[i] << endl;
+                        }
                    }
                }
            }
@ -177,10 +199,15 @@ namespace Foam

        inline bool addWatch(const label watchFd, const fileName& fName)
        {
+            if (inotifyFd_ < 0)
+            {
+                return false;
+            }
+
            // Add/retrieve watch on directory containing file
            label dirWatchID = inotify_add_watch
            (
-                fd,
+                inotifyFd_,
                fName.path().c_str(),
                IN_CLOSE_WRITE
            );
@ -189,7 +216,8 @@ namespace Foam
            {
                FatalErrorIn("addWatch(const label, const fileName&)")
                    << "Failed adding watch " << watchFd
-                    << " to directory " << fName
+                    << " to directory " << fName << " due to "
+                    << string(strerror(errno))
                    << exit(FatalError);
            }

@ -209,6 +237,11 @@ namespace Foam

        inline bool removeWatch(const label watchFd)
        {
+            if (inotifyFd_ < 0)
+            {
+                return false;
+            }
+
            dirWatches_[watchFd] = -1;
            return true;
        }
@ -263,11 +296,11 @@ void Foam::fileMonitor::checkFiles() const
        fd_set fdSet;
        // Add notify descriptor to select fd_set
        FD_ZERO(&fdSet);
-        FD_SET(watcher_->fd, &fdSet);
+        FD_SET(watcher_->inotifyFd_, &fdSet);

        int ready = select
        (
-            watcher_->fd+1,     // num filedescriptors in fdSet
+            watcher_->inotifyFd_+1,     // num filedescriptors in fdSet
            &fdSet,             // fdSet with only inotifyFd
            NULL,               // No writefds
            NULL,               // No errorfds
@ -280,15 +313,15 @@ void Foam::fileMonitor::checkFiles() const
                << "Problem in issuing select."
                << abort(FatalError);
        }
-        else if (FD_ISSET(watcher_->fd, &fdSet))
+        else if (FD_ISSET(watcher_->inotifyFd_, &fdSet))
        {
            // Read events
-            ssize_t nBytes = read(watcher_->fd, buffer, EVENT_BUF_LEN);
+            ssize_t nBytes = read(watcher_->inotifyFd_, buffer, EVENT_BUF_LEN);

            if (nBytes < 0)
            {
                FatalErrorIn("fileMonitor::updateStates(const fileName&)")
-                    << "read of " << watcher_->fd
+                    << "read of " << watcher_->inotifyFd_
                    << " failed with " << label(nBytes)
                    << abort(FatalError);
            }
@ -364,6 +397,7 @@ Foam::label Foam::fileMonitor::addWatch(const fileName& fName)
    label watchFd;

    label sz = freeWatchFds_.size();
+
    if (sz)
    {
        watchFd = freeWatchFds_[sz-1];
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@ -347,6 +347,7 @@ $(basicPolyPatches)/generic/genericPolyPatch.C
 constraintPolyPatches = $(polyPatches)/constraint
 $(constraintPolyPatches)/cyclic/cyclicPolyPatch.C
 $(constraintPolyPatches)/cyclicSlip/cyclicSlipPolyPatch.C
+$(constraintPolyPatches)/oldCyclic/oldCyclicPolyPatch.C
 $(constraintPolyPatches)/empty/emptyPolyPatch.C
 $(constraintPolyPatches)/nonuniformTransformCyclic/nonuniformTransformCyclicPolyPatch.C
 $(constraintPolyPatches)/processorCyclic/processorCyclicPolyPatch.C
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@ -245,13 +245,14 @@ Foam::Time::Time
    readLibs_(controlDict_, "libs"),
    functionObjects_(*this)
 {
+    setControls();
+
    // Time objects not registered so do like objectRegistry::checkIn ourselves.
    if (runTimeModifiable_)
    {
+        monitorPtr_.reset(new fileMonitor());
        controlDict_.watchIndex() = addWatch(controlDict_.filePath());
    }
-
-    setControls();
 }


@ -307,14 +308,20 @@ Foam::Time::Time
    readLibs_(controlDict_, "libs"),
    functionObjects_(*this)
 {
+    setControls();

    // Time objects not registered so do like objectRegistry::checkIn ourselves.
    if (runTimeModifiable_)
    {
-        controlDict_.watchIndex() = addWatch(controlDict_.filePath());
-    }
+        monitorPtr_.reset(new fileMonitor());

-    setControls();
+        // File might not exist yet.
+        fileName f(controlDict_.filePath());
+        if (f != fileName::null)
+        {
+            controlDict_.watchIndex() = addWatch(f);
+        }
+    }
 }


@ -388,18 +395,18 @@ Foam::Time::~Time()

 Foam::label Foam::Time::addWatch(const fileName& fName) const
 {
-    return monitor_.addWatch(fName);
+    return monitorPtr_().addWatch(fName);
 }


 bool Foam::Time::removeWatch(const label watchIndex) const
 {
-    return monitor_.removeWatch(watchIndex);
+    return monitorPtr_().removeWatch(watchIndex);
 }

 const Foam::fileName& Foam::Time::getFile(const label watchIndex) const
 {
-    return monitor_.getFile(watchIndex);
+    return monitorPtr_().getFile(watchIndex);
 }


@ -408,13 +415,13 @@ Foam::fileMonitor::fileState Foam::Time::getState
    const label watchFd
 ) const
 {
-    return monitor_.getState(watchFd);
+    return monitorPtr_().getState(watchFd);
 }


 void Foam::Time::setUnmodified(const label watchFd) const
 {
-    monitor_.setUnmodified(watchFd);
+    monitorPtr_().setUnmodified(watchFd);
 }


--- a/src/OpenFOAM/db/Time/Time.H
+++ b/src/OpenFOAM/db/Time/Time.H
@ -73,7 +73,7 @@ class Time
    // Private data

        //- file-change monitor for all registered files
-        mutable fileMonitor monitor_;
+        mutable autoPtr<fileMonitor> monitorPtr_;

        //- The controlDict
        IOdictionary controlDict_;
--- a/src/OpenFOAM/db/Time/TimeIO.C
+++ b/src/OpenFOAM/db/Time/TimeIO.C
@ -211,10 +211,7 @@ void Foam::Time::readModifiedObjects()
        // valid filePath).
        // Note: requires same ordering in objectRegistries on different
        // processors!
-        monitor_.updateStates(Pstream::parRun());
-
-//Pout<< "Time : runTimeModifiable_ and watchIndex:"
-//    << controlDict_.watchIndex() << endl;
+        monitorPtr_().updateStates(Pstream::parRun());

        // Time handling is special since controlDict_ is the one dictionary
        // that is not registered to any database.
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H
@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     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::oldCyclicPolyPatch
+
+Description
+    'old' style cyclic polyPatch with all faces in single patch. Does ordering
+    but cannot be used to run. Writes 'type cyclic' so foamUpgradeCyclics
+    can be run afterwards.
+    Used to get cyclics from mesh converters that assume cyclics in single
+    patch (e.g. fluent3DMeshToFoam)
+
+SourceFiles
+    oldCyclicPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oldCyclicPolyPatch_H
+#define oldCyclicPolyPatch_H
+
+#include "coupledPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class oldCyclicPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class oldCyclicPolyPatch
+:
+    public coupledPolyPatch
+{
+public:
+
+    enum transformType
+    {
+        UNKNOWN,
+        ROTATIONAL,
+        TRANSLATIONAL
+    };
+    static const NamedEnum<transformType, 3> transformTypeNames;
+
+
+private:
+
+    // Private data
+
+        //- Morph:angle between normals of neighbouring faces.
+        //  Used to split cyclic into halves.
+        scalar featureCos_;
+
+        //- Type of transformation - rotational or translational
+        transformType transform_;
+
+        // For rotation
+
+            //- Axis of rotation for rotational cyclics
+            vector rotationAxis_;
+
+            //- point on axis of rotation for rotational cyclics
+            point rotationCentre_;
+
+        // For translation
+
+            //- Translation vector
+            vector separationVector_;
+
+
+    // Private member functions
+
+        //- Find amongst selected faces the one with the largest area
+        static label findMaxArea(const pointField&, const faceList&);
+
+        void calcTransforms();
+
+        //- Calculate face centres
+        static pointField calcFaceCentres
+        (
+            const UList<face>&,
+            const pointField&
+        );
+
+        //- Get f[0] for all faces
+        static pointField getAnchorPoints
+        (
+            const UList<face>&,
+            const pointField&
+        );
+
+        // Face ordering
+
+            //- Find the two parts of the faces of pp using feature edges.
+            //  Returns true if successfull.
+            bool getGeometricHalves
+            (
+                const primitivePatch&,
+                labelList&,
+                labelList&
+            ) const;
+
+            //- Calculate geometric factors of the two halves.
+            void getCentresAndAnchors
+            (
+                const primitivePatch&,
+                const faceList& half0Faces,
+                const faceList& half1Faces,
+
+                pointField& ppPoints,
+                pointField& half0Ctrs,
+                pointField& half1Ctrs,
+                pointField& anchors0,
+                scalarField& tols
+            ) const;
+
+            //- Given matched faces matches the anchor point. Sets faceMap,
+            //  rotation. Returns true if all matched.
+            bool matchAnchors
+            (
+                const bool report,
+                const primitivePatch&,
+                const labelList&,
+                const pointField&,
+                const labelList&,
+                const faceList&,
+                const labelList&,
+                const scalarField&,
+
+                labelList& faceMap,
+                labelList& rotation
+            ) const;
+
+            //- For rotational cases, try to find a unique face on each side
+            //  of the cyclic.
+            label getConsistentRotationFace
+            (
+                const pointField& faceCentres
+            ) const;
+
+
+protected:
+
+    // Protected Member functions
+
+        //- Initialise the calculation of the patch geometry
+        virtual void initGeometry(PstreamBuffers&);
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry(PstreamBuffers&);
+
+        //- Initialise the patches for moving points
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
+
+        //- Correct patches after moving points
+        virtual void movePoints(PstreamBuffers&, const pointField&);
+
+        //- Initialise the update of the patch topology
+        virtual void initUpdateMesh(PstreamBuffers&);
+
+        //- Update of the patch topology
+        virtual void updateMesh(PstreamBuffers&);
+
+public:
+
+    //- Runtime type information
+    TypeName("oldCyclic");
+
+
+    // Constructors
+
+        //- Construct from components
+        oldCyclicPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from dictionary
+        oldCyclicPolyPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct as copy, resetting the boundary mesh
+        oldCyclicPolyPatch(const oldCyclicPolyPatch&, const polyBoundaryMesh&);
+
+        //- Construct given the original patch and resetting the
+        //  face list and boundary mesh information
+        oldCyclicPolyPatch
+        (
+            const oldCyclicPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        );
+
+        //- Construct and return a clone, resetting the boundary mesh
+        virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
+        {
+            return autoPtr<polyPatch>(new oldCyclicPolyPatch(*this, bm));
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new oldCyclicPolyPatch(*this, bm, index, newSize, newStart)
+            );
+        }
+
+
+    // Destructor
+
+        virtual ~oldCyclicPolyPatch();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Does this side own the patch ?
+            virtual bool owner() const
+            {
+                notImplemented("oldCyclicPolyPatch::owner()");
+                return true;
+            }
+
+            //- Transform a patch-based position from other side to this side
+            virtual void transformPosition(pointField& l) const
+            {
+                notImplemented("transformPosition(pointField&)");
+            }
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry
+        (
+            const primitivePatch& referPatch,
+            const UList<point>& thisCtrs,
+            const UList<point>& thisAreas,
+            const UList<point>& thisCc,
+            const UList<point>& nbrCtrs,
+            const UList<point>& nbrAreas,
+            const UList<point>& nbrCc
+        );
+
+        //- Initialize ordering for primitivePatch. Does not
+        //  refer to *this (except for name() and type() etc.)
+        virtual void initOrder
+        (
+            PstreamBuffers&,
+            const primitivePatch&
+        ) const;
+
+        //- Return new ordering for primitivePatch.
+        //  Ordering is -faceMap: for every face
+        //  index of the new face -rotation:for every new face the clockwise
+        //  shift of the original face. Return false if nothing changes
+        //  (faceMap is identity, rotation is 0), true otherwise.
+        virtual bool order
+        (
+            PstreamBuffers&,
+            const primitivePatch&,
+            labelList& faceMap,
+            labelList& rotation
+        ) const;
+
+        //- Write the polyPatch data as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/dynamicMesh/meshCut/directions/directions.C
+++ b/src/dynamicMesh/meshCut/directions/directions.C
@ -32,6 +32,7 @@ License
 #include "meshTools.H"
 #include "hexMatcher.H"
 #include "Switch.H"
+#include "globalMeshData.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

@ -206,7 +207,7 @@ Foam::vectorField Foam::directions::propagateDirection
        mesh,
        changedFaces,
        changedFacesInfo,
-        mesh.nCells()
+        mesh.globalData().nTotalCells()+1
    );

    const List<directionInfo>& cellInfo = directionCalc.allCellInfo();
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@ -2332,7 +2332,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
        seedFacesInfo.clear();

        // Iterate until no change. Now 2:1 face difference should be satisfied
-        levelCalc.iterate(mesh_.globalData().nTotalFaces());
+        levelCalc.iterate(mesh_.globalData().nTotalFaces()+1);


        // Now check point-connected cells (face-connected cells already ok):
@ -2836,7 +2836,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
        seedFacesInfo,
        allFaceInfo,
        allCellInfo,
-        mesh_.globalData().nTotalCells()
+        mesh_.globalData().nTotalCells()+1
    );


--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@ -337,6 +337,7 @@ $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C
 $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C

 finiteVolume/fvc/fvcMeshPhi.C
+finiteVolume/fvc/fvcSmooth/fvcSmooth.C

 general = cfdTools/general
 $(general)/findRefCell/findRefCell.C
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
@ -171,6 +171,23 @@ tmp<Field<Type> > cyclicFvPatchField<Type>::patchNeighbourField() const
 }


+template<class Type>
+const cyclicFvPatchField<Type>& cyclicFvPatchField<Type>::neighbourPatchField()
+const
+{
+    const GeometricField<Type, fvPatchField, volMesh>& fld =
+    static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
+    (
+        this->internalField()
+    );
+
+    return refCast<const cyclicFvPatchField<Type> >
+    (
+        fld.boundaryField()[this->cyclicPatch().neighbPatchID()]
+    );
+}
+
+
 template<class Type>
 void cyclicFvPatchField<Type>::updateInterfaceMatrix
 (
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
@ -150,9 +150,12 @@ public:

        // Evaluation functions

-            //- Return neighbour coupled given internal cell data
+            //- Return neighbour coupled internal cell data
            tmp<Field<Type> > patchNeighbourField() const;

+            //- Return reference to neighbour patchField
+            const cyclicFvPatchField<Type>& neighbourPatchField() const;
+
            //- Update result field based on interface functionality
            virtual void updateInterfaceMatrix
            (
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@ -140,7 +140,17 @@ public:
            //- Return the "jump" across the patch.
            virtual tmp<Field<Type> > jump() const
            {
-                return jump_;
+                if (this->cyclicPatch().owner())
+                {
+                    return jump_;
+                }
+                else
+                {
+                    return refCast<const fanFvPatchField<Type> >
+                    (
+                        this->neighbourPatchField()
+                    ).jump();
+                }
            }


--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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 "fvcSmooth.H"
+#include "volFields.H"
+#include "FaceCellWave.H"
+#include "smoothData.H"
+#include "sweepData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::smoothData::maxRatio = 1.0;
+
+void Foam::fvc::smooth
+(
+    volScalarField& field,
+    const scalar coeff
+)
+{
+    const fvMesh& mesh = field.mesh();
+    smoothData::maxRatio = 1 + coeff;
+
+    DynamicList<label> changedFaces(mesh.nFaces()/100 + 100);
+    DynamicList<smoothData> changedFacesInfo(changedFaces.size());
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighbour = mesh.neighbour();
+
+    forAll(owner, facei)
+    {
+        const label own = owner[facei];
+        const label nbr = neighbour[facei];
+
+        // Check if owner value much larger than neighbour value or vice versa
+        if (field[own] > smoothData::maxRatio*field[nbr])
+        {
+            changedFaces.append(facei);
+            changedFacesInfo.append(smoothData(field[own]));
+        }
+        else if (field[nbr] > smoothData::maxRatio*field[own])
+        {
+            changedFaces.append(facei);
+            changedFacesInfo.append(smoothData(field[nbr]));
+        }
+    }
+
+    // Insert all faces of coupled patches - FaceCellWave will correct them
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFacei)
+            {
+                label facei = patch.start() + patchFacei;
+                label own = mesh.faceOwner()[facei];
+
+                changedFaces.append(facei);
+                changedFacesInfo.append(smoothData(field[own]));
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+
+    // Set initial field on cells
+    List<smoothData> cellData(mesh.nCells());
+
+    forAll(field, celli)
+    {
+        cellData[celli] = field[celli];
+    }
+
+    // Set initial field on faces
+    List<smoothData> faceData(mesh.nFaces());
+
+
+    // Propagate information over whole domain
+    FaceCellWave<smoothData > smoothData
+    (
+        mesh,
+        changedFaces,
+        changedFacesInfo,
+        faceData,
+        cellData,
+        mesh.globalData().nTotalCells()    // max iterations
+    );
+
+    forAll(field, celli)
+    {
+        field[celli] = cellData[celli].value();
+    }
+
+    field.correctBoundaryConditions();
+}
+
+
+void Foam::fvc::spread
+(
+    volScalarField& field,
+    const volScalarField& alpha,
+    const label nLayers
+)
+{
+    const fvMesh& mesh = field.mesh();
+    smoothData::maxRatio = 1;
+
+    DynamicList<label> changedFaces(mesh.nFaces()/100 + 100);
+    DynamicList<smoothData> changedFacesInfo(changedFaces.size());
+
+    // Set initial field on cells
+    List<smoothData> cellData(mesh.nCells());
+
+    forAll(field, celli)
+    {
+        cellData[celli] = field[celli];
+    }
+
+    // Set initial field on faces
+    List<smoothData> faceData(mesh.nFaces());
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighbour = mesh.neighbour();
+
+    forAll(owner, facei)
+    {
+        const label own = owner[facei];
+        const label nbr = neighbour[facei];
+
+        if
+        (
+            (alpha[own] > 0.01 && alpha[own] < 0.99)
+         || (alpha[nbr] > 0.01 && alpha[nbr] < 0.99)
+        )
+        {
+            if (mag(alpha[own] - alpha[nbr]) > 0.2)
+            {
+                changedFaces.append(facei);
+                changedFacesInfo.append
+                (
+                    smoothData(max(field[own], field[nbr]))
+                );
+            }
+        }
+    }
+
+    // Insert all faces of coupled patches - FaceCellWave will correct them
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFacei)
+            {
+                label facei = patch.start() + patchFacei;
+                label own = mesh.faceOwner()[facei];
+
+                scalarField alphapn =
+                    alpha.boundaryField()[patchi].patchNeighbourField();
+
+                if
+                (
+                    (alpha[own] > 0.01 && alpha[own] < 0.99)
+                 || (alphapn[patchFacei] > 0.01 && alphapn[patchFacei] < 0.99)
+                )
+                {
+                    if (mag(alpha[own] - alphapn[patchFacei]) > 0.2)
+                    {
+                        changedFaces.append(facei);
+                        changedFacesInfo.append(smoothData(field[own]));
+                    }
+                }
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+
+    // Propagate information over whole domain
+    FaceCellWave<smoothData> smoothData
+    (
+        mesh,
+        faceData,
+        cellData
+    );
+
+    smoothData.setFaceInfo(changedFaces, changedFacesInfo);
+
+    smoothData.iterate(nLayers);
+
+    forAll(field, celli)
+    {
+        field[celli] = cellData[celli].value();
+    }
+
+    field.correctBoundaryConditions();
+}
+
+
+void Foam::fvc::sweep
+(
+    volScalarField& field,
+    const volScalarField& alpha,
+    const label nLayers
+)
+{
+    const fvMesh& mesh = field.mesh();
+
+    DynamicList<label> changedFaces(mesh.nFaces()/100 + 100);
+    DynamicList<sweepData> changedFacesInfo(changedFaces.size());
+
+    // Set initial field on cells
+    List<sweepData> cellData(mesh.nCells());
+
+    // Set initial field on faces
+    List<sweepData> faceData(mesh.nFaces());
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighbour = mesh.neighbour();
+    const vectorField& Cf = mesh.faceCentres();
+
+    forAll(owner, facei)
+    {
+        const label own = owner[facei];
+        const label nbr = neighbour[facei];
+
+        if (mag(alpha[own] - alpha[nbr]) > 0.2)
+        {
+            changedFaces.append(facei);
+            changedFacesInfo.append
+            (
+                sweepData(max(field[own], field[nbr]),
+                Cf[facei])
+            );
+        }
+    }
+
+    // Insert all faces of coupled patches - FaceCellWave will correct them
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFacei)
+            {
+                label facei = patch.start() + patchFacei;
+                label own = mesh.faceOwner()[facei];
+
+                scalarField alphapn =
+                    alpha.boundaryField()[patchi].patchNeighbourField();
+
+                if (mag(alpha[own] - alphapn[patchFacei]) > 0.2)
+                {
+                    changedFaces.append(facei);
+                    changedFacesInfo.append
+                    (
+                        sweepData(field[own], Cf[facei])
+                    );
+                }
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+
+    // Propagate information over whole domain
+    FaceCellWave<sweepData> sweepData
+    (
+        mesh,
+        faceData,
+        cellData
+    );
+
+    sweepData.setFaceInfo(changedFaces, changedFacesInfo);
+
+    sweepData.iterate(nLayers);
+
+    forAll(field, celli)
+    {
+        if (cellData[celli].valid())
+        {
+            field[celli] = max(field[celli], cellData[celli].value());
+        }
+    }
+
+    field.correctBoundaryConditions();
+}
+
+
+// ************************************************************************* //
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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/>.
+
+function
+    Foam::fvc::smooth
+
+Description
+    Function that uses smoothData to apply spatial smoothing of a
+    volume field using the FaceCellWave algorithm
+
+SourceFiles
+    fvcSmooth.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fvcSmooth_H
+#define fvcSmooth_H
+
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fvc
+{
+    void smooth
+    (
+        volScalarField& field,
+        const scalar coeff
+    );
+
+    void spread
+    (
+        volScalarField& field,
+        const volScalarField& alpha,
+        const label nLayers
+    );
+
+    void sweep
+    (
+        volScalarField& field,
+        const volScalarField& alpha,
+        const label nLayers
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H
@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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::smoothData
+
+Description
+    Helper class used by the smoothVolField class
+
+Files
+    smoothData.H
+    smoothDataI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef smoothData_H
+#define smoothData_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class smoothData Declaration
+\*---------------------------------------------------------------------------*/
+
+class smoothData
+{
+    scalar value_;
+
+    // Private Member Functions
+
+        //- Update - gets information from neighbouring face/cell and
+        //  uses this to update itself (if necessary) and return true
+        inline bool update
+        (
+            const smoothData& svf,
+            const scalar scale,
+            const scalar tol
+        );
+
+
+public:
+
+    // Static data members
+
+        //- Field fraction
+        static scalar maxRatio;
+
+
+    // Constructors
+
+        //- Construct null
+        inline smoothData();
+
+        //- Construct from inverse field value
+        inline smoothData(const scalar value);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return value
+            scalar value() const
+            {
+                return value_;
+            }
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value
+            inline bool valid() const;
+
+            //- Check for identical geometrical data
+            //  Used for cyclics checking
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const smoothData&,
+                const scalar
+            ) const;
+
+            //- Convert any absolute coordinates into relative to
+            //  (patch)face centre
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Reverse of leaveDomain
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Apply rotation matrix to any coordinates
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&
+            );
+
+            //- Influence of neighbouring face
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const smoothData& svf,
+                const scalar tol
+            );
+
+            //- Influence of neighbouring cell
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const smoothData& svf,
+                const scalar tol
+            );
+
+            //- Influence of different value on same face
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const smoothData& svf,
+                const scalar tol
+            );
+
+
+        // Member Operators
+
+            inline void operator=(const scalar value);
+
+            // Needed for List IO
+            inline bool operator==(const smoothData&) const;
+
+            inline bool operator!=(const smoothData&) const;
+
+
+        // IOstream Operators
+
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const smoothData& svf
+            )
+            {
+                return os  << svf.value_;
+            }
+
+            friend Istream& operator>>(Istream& is, smoothData& svf)
+            {
+                return is  >> svf.value_;
+            }
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "smoothDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H
@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline bool Foam::smoothData::update
+(
+    const smoothData::smoothData& svf,
+    const scalar scale,
+    const scalar tol
+)
+{
+    if (!valid() || (value_ < VSMALL))
+    {
+        // My value not set - take over neighbour
+        value_ = svf.value()/scale;
+
+        // Something changed - let caller know
+        return true;
+    }
+    else if (svf.value() > (1 + tol)*scale*value_)
+    {
+        // Neighbour is too big for me - Up my value
+        value_ = svf.value()/scale;
+
+        // Something changed - let caller know
+        return true;
+    }
+    else
+    {
+        // Neighbour is not too big for me or change is too small
+        // Nothing changed
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::smoothData::smoothData()
+:
+    value_(-GREAT)
+{}
+
+
+inline Foam::smoothData::smoothData(const scalar value)
+:
+    value_(value)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::smoothData::valid() const
+{
+    return value_ > -SMALL;
+}
+
+
+inline bool Foam::smoothData::sameGeometry
+(
+    const polyMesh&,
+    const smoothData&,
+    const scalar
+) const
+{
+    return true;
+}
+
+
+inline void Foam::smoothData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point&
+)
+{}
+
+
+inline void Foam::smoothData::transform
+(
+    const polyMesh&,
+    const tensor&
+)
+{}
+
+
+inline void Foam::smoothData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point&
+)
+{}
+
+
+inline bool Foam::smoothData::updateCell
+(
+    const polyMesh&,
+    const label,
+    const label,
+    const smoothData& svf,
+    const scalar tol
+)
+{
+    // Take over info from face if more than deltaRatio larger
+    return update(svf, maxRatio, tol);
+}
+
+
+inline bool Foam::smoothData::updateFace
+(
+    const polyMesh&,
+    const label,
+    const label,
+    const smoothData& svf,
+    const scalar tol
+)
+{
+    // Take over information from cell without any scaling (scale = 1.0)
+    return update(svf, 1.0, tol);
+}
+
+
+// Update this (face) with coupled face information.
+inline bool Foam::smoothData::updateFace
+(
+    const polyMesh&,
+    const label,
+    const smoothData& svf,
+    const scalar tol
+)
+{
+    // Take over information from coupled face without any scaling (scale = 1.0)
+    return update(svf, 1.0, tol);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::smoothData::operator=
+(
+    const scalar value
+)
+{
+    value_ = value;
+}
+
+
+inline bool Foam::smoothData::operator==
+(
+    const smoothData& rhs
+) const
+{
+    return value_ == rhs.value();
+}
+
+
+inline bool Foam::smoothData::operator!=
+(
+    const smoothData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H
@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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::sweepData
+
+Description
+    Helper class used by fvcSmooth
+
+Files
+    sweepData.H
+    sweepDataI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sweepData_H
+#define sweepData_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sweepData Declaration
+\*---------------------------------------------------------------------------*/
+
+class sweepData
+{
+    scalar value_;
+    point origin_;
+
+    // Private Member Functions
+
+        //- Update - gets information from neighbouring face/cell and
+        //  uses this to update itself (if necessary) and return true
+        inline bool update
+        (
+            const sweepData& svf,
+            const point& position,
+            const scalar tol
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline sweepData();
+
+        //- Construct from component
+        inline sweepData(const scalar value, const point& origin);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return value
+            scalar value() const
+            {
+                return value_;
+            }
+
+            //- Return origin
+            const point& origin() const
+            {
+                return origin_;
+            }
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value
+            inline bool valid() const;
+
+            //- Check for identical geometrical data
+            //  Used for cyclics checking
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const sweepData&,
+                const scalar
+            ) const;
+
+            //- Convert any absolute coordinates into relative to
+            //  (patch)face centre
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Reverse of leaveDomain
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Apply rotation matrix to any coordinates
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&
+            );
+
+            //- Influence of neighbouring face
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const sweepData& svf,
+                const scalar tol
+            );
+
+            //- Influence of neighbouring cell
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const sweepData& svf,
+                const scalar tol
+            );
+
+            //- Influence of different value on same face
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const sweepData& svf,
+                const scalar tol
+            );
+
+
+        // Member Operators
+
+            inline void operator=(const scalar value);
+
+            // Needed for List IO
+            inline bool operator==(const sweepData&) const;
+
+            inline bool operator!=(const sweepData&) const;
+
+
+        // IOstream Operators
+
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const sweepData& svf
+            )
+            {
+                return os << svf.value_ << svf.origin_;
+            }
+
+            friend Istream& operator>>(Istream& is, sweepData& svf)
+            {
+                return is >> svf.value_ >> svf.origin_;
+            }
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sweepDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H
@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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 "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline bool Foam::sweepData::update
+(
+    const sweepData::sweepData& svf,
+    const point& position,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        operator=(svf);
+        return true;
+    }
+
+    scalar myDist2 = magSqr(position - origin());
+
+    if (myDist2 < SMALL)
+    {
+        if (svf.value() > value())
+        {
+            operator=(svf);
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    scalar dist2 = magSqr(position - svf.origin());
+
+    if (dist2 < myDist2)
+    {
+        operator=(svf);
+        return true;
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::sweepData::sweepData()
+:
+    value_(-GREAT),
+    origin_(vector::max)
+{}
+
+
+inline Foam::sweepData::sweepData(const scalar value, const point& origin)
+:
+    value_(value),
+    origin_(origin)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::sweepData::valid() const
+{
+    return value_ > -SMALL;
+}
+
+
+inline bool Foam::sweepData::sameGeometry
+(
+    const polyMesh&,
+    const sweepData&,
+    const scalar
+) const
+{
+    return true;
+}
+
+
+inline void Foam::sweepData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point& faceCentre
+)
+{
+    origin_ -= faceCentre;
+}
+
+
+inline void Foam::sweepData::transform
+(
+    const polyMesh&,
+    const tensor& rotTensor
+)
+{
+    origin_ = Foam::transform(rotTensor, origin_);
+}
+
+
+inline void Foam::sweepData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point& faceCentre
+)
+{
+    // back to absolute form
+    origin_ += faceCentre;
+}
+
+
+inline bool Foam::sweepData::updateCell
+(
+    const polyMesh& mesh,
+    const label thisCellI,
+    const label,
+    const sweepData& svf,
+    const scalar tol
+)
+{
+    return update(svf, mesh.cellCentres()[thisCellI], tol);
+}
+
+
+inline bool Foam::sweepData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const label,
+    const sweepData& svf,
+    const scalar tol
+)
+{
+    return update(svf, mesh.faceCentres()[thisFaceI], tol);
+}
+
+
+// Update this (face) with coupled face information.
+inline bool Foam::sweepData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const sweepData& svf,
+    const scalar tol
+)
+{
+    return update(svf, mesh.faceCentres()[thisFaceI], tol);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::sweepData::operator=
+(
+    const scalar value
+)
+{
+    value_ = value;
+}
+
+
+inline bool Foam::sweepData::operator==
+(
+    const sweepData& rhs
+) const
+{
+    return origin() == rhs.origin();
+}
+
+
+inline bool Foam::sweepData::operator!=
+(
+    const sweepData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
--- a/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C
+++ b/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C
@ -115,7 +115,7 @@ void Foam::inverseFaceDistanceDiffusivity::correct()
        mesh,
        changedFaces,
        faceDist,
-        mesh.globalData().nTotalCells() // max iterations
+        mesh.globalData().nTotalCells()+1 // max iterations
    );

    const List<wallPoint>& faceInfo = waveInfo.allFaceInfo();
--- a/src/meshTools/cellClassification/cellClassification.C
+++ b/src/meshTools/cellClassification/cellClassification.C
@ -322,7 +322,7 @@ void Foam::cellClassification::markCells
        changedFaces,                       // Labels of changed faces
        changedFacesInfo,                   // Information on changed faces
        cellInfoList,                       // Information on all cells
-        mesh_.globalData().nTotalCells()    // max iterations
+        mesh_.globalData().nTotalCells()+1  // max iterations
    );

    // Get information out of cellInfoList
--- a/src/meshTools/cellDist/patchWave/patchDataWave.C
+++ b/src/meshTools/cellDist/patchWave/patchDataWave.C
@ -225,7 +225,7 @@ void Foam::patchDataWave<TransferType>::correct()
        mesh(),
        changedFaces,
        faceDist,
-        mesh().globalData().nTotalCells() // max iterations
+        mesh().globalData().nTotalCells()+1 // max iterations
    );


--- a/src/meshTools/cellDist/patchWave/patchWave.C
+++ b/src/meshTools/cellDist/patchWave/patchWave.C
@ -186,7 +186,7 @@ void Foam::patchWave::correct()
        mesh(),
        changedFaces,
        faceDist,
-        mesh().globalData().nTotalCells()   // max iterations
+        mesh().globalData().nTotalCells()+1 // max iterations
    );


--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
@ -148,7 +148,7 @@ Foam::labelList Foam::structuredDecomp::decompose
        patchData,
        faceData,
        cellData,
-        mesh.globalData().nTotalCells()
+        mesh.globalData().nTotalCells()+1
    );

    // And extract
--- a/src/parallel/decompose/ptscotchDecomp/Make/options
+++ b/src/parallel/decompose/ptscotchDecomp/Make/options
@ -3,7 +3,7 @@ sinclude $(RULES)/mplib$(WM_MPLIB)

 EXE_INC = \
    $(PFLAGS) $(PINC) \
-    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1.10a/include \
    -I/usr/include/scotch \
    -I../decompositionMethods/lnInclude

--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@ -180,17 +180,31 @@ Foam::label Foam::ptscotchDecomp::decompose
 //            Info<< "Dumping Scotch graph file to " << str.name() << endl
 //                << "Use this in combination with gpart." << endl;
 //
-//            label version = 0;
+//            // Distributed graph file (.grf)
+//            label version = 2;
 //            str << version << nl;
-//            // Numer of vertices
-//            str << xadj.size()-1 << ' ' << adjncy.size() << nl;
+//            // Number of files
+
+//            // Number of files (procglbnbr)
+//            str << Pstream::nProcs();
+//            // My file number (procloc)
+//            str << ' ' << Pstream::myProcNo() << nl;
+//
+//            // Total number of vertices (vertglbnbr)
+//            str << returnReduce(mesh.nCells(), sumOp<label>());
+//            // Total number of connections (edgeglbnbr)
+//            str << ' ' << returnReduce(xadj[mesh.nCells()], sumOp<label>())
+//                << nl;
+//            // Local number of vertices (vertlocnbr)
+//            str << mesh.nCells();
+//            // Local number of connections (edgelocnbr)
+//            str << ' ' << xadj[mesh.nCells()] << nl;
 //            // Numbering starts from 0
 //            label baseval = 0;
-//            // Has weights?
-//            label hasEdgeWeights = 0;
-//            label hasVertexWeights = 0;
-//            label numericflag = 10*hasEdgeWeights+hasVertexWeights;
-//            str << baseval << ' ' << numericflag << nl;
+//            // Start of my global vertices (procdsptab[proclocnum])
+//            str << ' ' << globalCells.toGlobal(0);
+//            100*hasVertlabels+10*hasEdgeWeights+1*hasVertWeighs
+//            str << ' ' << "0" << nl;
 //            for (label cellI = 0; cellI < xadj.size()-1; cellI++)
 //            {
 //                label start = xadj[cellI];
--- a/src/parallel/decompose/scotchDecomp/Make/options
+++ b/src/parallel/decompose/scotchDecomp/Make/options
@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1.10a/include \
    -I/usr/include/scotch \
    -I../decompositionMethods/lnInclude

--- a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C
+++ b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C
@ -137,7 +137,7 @@ void smoothDelta::calcDelta()
        changedFacesInfo,
        faceDeltaData,
        cellDeltaData,
-        mesh_.globalData().nTotalCells()    // max iterations
+        mesh_.globalData().nTotalCells()+1   // max iterations
    );

    forAll(delta_, cellI)
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@ -172,7 +172,7 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
    (
        IOobject
        (
-            "epsilon",
+            "nut",
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U
@ -15,46 +15,46 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-dimensions      [ 0 1 -1 0 0 0 0 ];
+dimensions      [0 1 -1 0 0 0 0];

-internalField   uniform ( 0 0 0 );
+internalField   uniform (0 0 0);

 boundaryField
 {
    inlet
    {
        type            pressureInletOutletVelocity;
-        value           uniform ( 0 0 0 );
+        value           uniform (0 0 0);
    }
    outlet1
    {
        type            inletOutlet;
-        inletValue      uniform ( 0 0 0 );
-        value           uniform ( 0 0 0 );
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
    }
    outlet2
    {
        type            inletOutlet;
-        inletValue      uniform ( 0 0 0 );
-        value           uniform ( 0 0 0 );
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
    }
    baffles
    {
        type            fixedValue;
-        value           uniform ( 0 0 0 );
+        value           uniform (0 0 0);
    }
    fan_half0
    {
        type            cyclic;
    }
+    fan_half1
+    {
+        type            cyclic;
+    }
    defaultFaces
    {
        type            fixedValue;
-        value           uniform ( 0 0 0 );
-    }
-    fan_half1
-    {
-        type            cyclic;
+        value           uniform (0 0 0);
    }
 }

--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon
@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-dimensions      [ 0 2 -3 0 0 0 0 ];
+dimensions      [0 2 -3 0 0 0 0];

 internalField   uniform 1;

@ -42,21 +42,27 @@ boundaryField
    baffles
    {
        type            epsilonWallFunction;
-        value           uniform 1;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
    }
    fan_half0
    {
        type            cyclic;
    }
+    fan_half1
+    {
+        type            cyclic;
+    }
    defaultFaces
    {
        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
        value           uniform 1;
    }
-    fan_half1
-    {
-        type            cyclic;
-    }
 }


--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k
@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-dimensions      [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];

 internalField   uniform 1;

@ -42,21 +42,21 @@ boundaryField
    baffles
    {
        type            kqRWallFunction;
-        value           uniform 1;
+        value           uniform 0;
    }
    fan_half0
    {
        type            cyclic;
    }
+    fan_half1
+    {
+        type            cyclic;
+    }
    defaultFaces
    {
        type            kqRWallFunction;
        value           uniform 1;
    }
-    fan_half1
-    {
-        type            cyclic;
-    }
 }


--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda
@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-dimensions      [ 0 2 -1 0 0 0 0 ];
+dimensions      [0 2 -1 0 0 0 0];

 internalField   uniform 0;

@ -41,14 +41,14 @@ boundaryField
    {
        type            cyclic;
    }
-    defaultFaces
-    {
-        type            zeroGradient;
-    }
    fan_half1
    {
        type            cyclic;
    }
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
 }


--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut
@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-dimensions      [ 0 2 -1 0 0 0 0 ];
+dimensions      [0 2 -1 0 0 0 0];

 internalField   uniform 0;

@ -39,21 +39,27 @@ boundaryField
    baffles
    {
        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
        value           uniform 0;
    }
    fan_half0
    {
        type            cyclic;
    }
-    defaultFaces
-    {
-        type            nutkWallFunction;
-        value           uniform 0;
-    }
    fan_half1
    {
        type            cyclic;
    }
+    defaultFaces
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
 }


--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p
@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

-dimensions      [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];

 internalField   uniform 100000;

@ -24,14 +24,12 @@ boundaryField
    inlet
    {
        type            timeVaryingTotalPressure;
-        p0              100040;
-        outOfBounds     clamp;
-        fileName        "$FOAM_CASE/constant/p0vsTime";
-        U               U;
-        phi             phi;
        rho             none;
        psi             none;
        gamma           1;
+        p0              100040;
+        fileName        "$FOAM_CASE/constant/p0vsTime";
+        outOfBounds     clamp;
        value           uniform 100040;
    }
    outlet1
@ -52,19 +50,19 @@ boundaryField
    {
        type            fan;
        patchType       cyclic;
-        f               2 ( 50 -0.1 );
-        value           $internalField;
-    }
-    defaultFaces
-    {
-        type            zeroGradient;
+        f               2(100 -0.1);
+        value           uniform 0;
    }
    fan_half1
    {
        type            fan;
        patchType       cyclic;
-        f               2 ( 50 -0.1 );
-        value           $internalField;
+        f               2(100 -0.1);
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            zeroGradient;
    }
 }

--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict
@ -28,7 +28,6 @@ endTime         1;
 deltaT          0.001;

 writeControl    adjustableRunTime;
-
 writeInterval   0.1;

 purgeWrite      0;
@ -47,7 +46,7 @@ runTimeModifiable true;

 adjustTimeStep  yes;

-maxCo           5;
+maxCo           3;

 libs
 (