Merge branch 'master' into cvm

ReleaseNotes-dev

@@ -66,6 +66,11 @@
       -doc              display application documentation in browser
       -help             print the usage
 
+*** *New* basicSolidThermo solids thermophysical library
+    + Used in all conjugate heat transfer solvers
+      + constant properties
+      + temperature dependent properties
+      + temperature and direction (in local coordinate system) dependent properties
 *** *New* Surface film library
     + Creation of films by particle addition, or initial film distribution
     + Coupled with the lagrangian/intermediate cloud hierarchy library
@@ -86,7 +91,14 @@
 *** *New* ptscotch decomposition method
 *** *Updated* particle tracking algorithm
 *** *Updated* split cyclics into two separate patches. See doc/changed/splitCyclics.txt
+  * *Updated* interpolation (volPointInterpolation) now works without the
+    globalPointPatch. Moving mesh cases can now be run non-parallel and
+    continued in parallel and reconstructed without any limitation.
 *** *New* compact binary I/O for faces and cells. This speeds up reading/writing meshes in binary.
+*** *Updated* runTimeModifiable
+    + on linux uses inotify instead of time stamps - more efficient for large
+      numbers of monitored files. No more fileModificationSkew needed.
+    + single integer reduction instead of one reduction per monitored file.
 * Solvers
   A number of new solvers have been developed for a range of engineering
   applications.  There has been a set of improvements to certain classes of
@@ -125,6 +137,8 @@
         (nonuniformTransform)cyclic <zoneA>_<zoneB>
       + extrudes into master direction (i.e. away from the owner cell
         if flipMap is false)
+    + =topoSet=: replacement of cellSet,faceSet,pointSet utilities.
+      Comparable to a dictionary driven =setSet= utility.
 *** Updated utilities
     + =setFields=: optionally use faceSets to set patch values (see e.g. hotRoom tutorial).
     + =blockMesh=: specify patches via dictionary instead of type only. This
@@ -133,6 +147,8 @@
 * Post-processing
   + =foamToEnsight=: parallel continuous data. new =-nodeValues= option to generate and output nodal
     field data.
+  + =singleCellMesh=: new utility to convert mesh and fields to a single cell
+    mesh. Great for postprocessing.
   + Function objects:
     + =residualControl=: new function object to allow users to terminate steady
     state calculations when the defined residual levels are achieved
@@ -144,7 +160,9 @@
     + =streamLine=: generate streamlines; ouputs both trajectory and field data
 
 * New tutorials
-  There is a large number of new tutorials to support the new solvers in the
+  There is a large number of new tutorials for existing and new solvers in the
   release.
   + =reactingParcelFilmFoam= tutorials:
     + multipleBoxes, hotBoxes, panel, evaporationTest
+  + =interDyMFoam= tutorials:
+    + testTubeMixer: showcases =solidBodyMotionFunction=

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -34,7 +34,6 @@ Description
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
-#include "compressibleCourantNo.H"
 #include "basicSolidThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C

@@ -37,20 +37,14 @@ Foam::scalar Foam::compressibleCourantNo
     scalar CoNum = 0.0;
     scalar meanCoNum = 0.0;
 
-    //- Can have fluid domains with 0 cells so do not test.
-    //if (mesh.nInternalFaces())
-    {
-        surfaceScalarField SfUfbyDelta =
-            mesh.surfaceInterpolation::deltaCoeffs()
-          * mag(phi)
-          / fvc::interpolate(rho);
+    scalarField sumPhi =
+        fvc::surfaceSum(mag(phi))().internalField()
+       /rho.internalField();
 
-        CoNum = max(SfUfbyDelta/mesh.magSf())
-            .value()*runTime.deltaT().value();
+    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 
-        meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-            .value()*runTime.deltaT().value();
-    }
+    meanCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
 
     Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
         << " max: " << CoNum << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C

@@ -37,18 +37,16 @@ Foam::scalar Foam::solidRegionDiffNo
     scalar DiNum = 0.0;
     scalar meanDiNum = 0.0;
 
-    //- Can have fluid domains with 0 cells so do not test.
-    if (mesh.nInternalFaces())
-    {
-        surfaceScalarField KrhoCpbyDelta =
-            mesh.surfaceInterpolation::deltaCoeffs()
-          * fvc::interpolate(K)
-          / fvc::interpolate(Cprho);
+    //- Take care: can have fluid domains with 0 cells so do not test for
+    //  zero internal faces.
+    surfaceScalarField KrhoCpbyDelta =
+        mesh.surfaceInterpolation::deltaCoeffs()
+      * fvc::interpolate(K)
+      / fvc::interpolate(Cprho);
 
-        DiNum = max(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
+    DiNum = gMax(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
 
-        meanDiNum = (average(KrhoCpbyDelta)).value()*runTime.deltaT().value();
-    }
+    meanDiNum = (average(KrhoCpbyDelta)).value()*runTime.deltaT().value();
 
     Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
         << " max: " << DiNum << endl;

applications/solvers/multiphase/bubbleFoam/createRASTurbulence.H

@@ -55,7 +55,7 @@
     (
         dimensionedScalar::lookupOrAddToDict
         (
-            "alphaEps",
+            "alphak",
             kEpsilonDict,
             1.0
         )

applications/solvers/multiphase/interFoam/Allwclean

@@ -6,5 +6,6 @@ wclean
 wclean interDyMFoam
 wclean MRFInterFoam
 wclean porousInterFoam
+wclean LTSInterFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/interFoam/Allwmake

@@ -6,5 +6,6 @@ wmake
 wmake interDyMFoam
 wmake MRFInterFoam
 wmake porousInterFoam
+wmake LTSInterFoam
 
 # ----------------------------------------------------------------- end-of-file

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

applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C

@@ -23,39 +23,60 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "compressibleCourantNo.H"
-#include "fvc.H"
+#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"
 
-Foam::scalar Foam::compressibleCourantNo
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::MULES::explicitLTSSolve
 (
-    const fvMesh& mesh,
-    const Time& runTime,
-    const volScalarField& rho,
-    const surfaceScalarField& phi
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
 )
 {
-    scalar CoNum = 0.0;
-    scalar meanCoNum = 0.0;
+    explicitLTSSolve
+    (
+        geometricOneField(),
+        psi,
+        phi,
+        phiPsi,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
 
-    //- Can have fluid domains with 0 cells so do not test.
-    //if (mesh.nInternalFaces())
-    {
-        surfaceScalarField SfUfbyDelta =
-            mesh.surfaceInterpolation::deltaCoeffs()
-          * mag(phi)
-          / fvc::interpolate(rho);
 
-        CoNum = max(SfUfbyDelta/mesh.magSf())
-            .value()*runTime.deltaT().value();
-
-        meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-            .value()*runTime.deltaT().value();
-    }
-
-    Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
-        << " max: " << CoNum << endl;
-
-    return CoNum;
+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
+    );
 }
 
 

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

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

applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files

@@ -0,0 +1,4 @@
+LTSInterFoam.C
+MULES.C
+
+EXE = $(FOAM_APPBIN)/LTSInterFoam

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

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;
+}

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;

applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H

@@ -0,0 +1,31 @@
+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),
+    zeroGradientFvPatchScalarField::typeName
+);
+
+volScalarField rSubDeltaT
+(
+    IOobject
+    (
+        "rSubDeltaT",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+);

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -0,0 +1,133 @@
+{
+    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)
+    );
+
+    scalar alphaSpreadDiff
+    (
+        piso.lookupOrDefault<label>("alphaSpreadDiff", 0.2)
+    );
+
+    scalar alphaSpreadMax
+    (
+        piso.lookupOrDefault<label>("alphaSpreadMax", 0.99)
+    );
+
+    scalar alphaSpreadMin
+    (
+        piso.lookupOrDefault<label>("alphaSpreadMin", 0.01)
+    );
+
+    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 from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(mag(rhoPhi))().dimensionedInternalField()
+       /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+        rDeltaT.dimensionedInternalField() = max
+        (
+            rDeltaT.dimensionedInternalField(),
+            pos(alpha1.dimensionedInternalField() - 0.01)
+           *pos(0.99 - alpha1.dimensionedInternalField())
+           *fvc::surfaceSum(mag(phi))().dimensionedInternalField()
+           /((2*maxAlphaCo)*mesh.V())
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    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,
+            alphaSpreadDiff,
+            alphaSpreadMax,
+            alphaSpreadMin
+        );
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
+    }
+
+    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;
+}

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/Make/options

@@ -5,7 +5,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
-EXE_LIBS = \
+LIB_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
     -lfiniteVolume

applications/test/DynamicField/DynamicFieldTest.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,9 +25,9 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+#include "point.H"
 #include "DynamicField.H"
 #include "IOstreams.H"
-#include "labelField.H"
 
 using namespace Foam;
 
@@ -36,44 +36,62 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-    {
-        DynamicField<label> dl(10);
-        Pout<< "null construct dl:" << dl << endl;
-        dl.append(3);
-        dl.append(2);
-        dl.append(1);
-        Pout<< "appending : dl:" << dl << endl;
+    DynamicField<point, 0, 10, 11> testField;
+    DynamicField<point, 0, 10, 11> testField2;
 
-        dl[2] *= 10;
-        Pout<< "assigning : dl:" << dl << endl;
-    }
+    testField.setSize(5);
+    testField2.setSize(5);
 
-    {
-        DynamicField<label> dl(IStringStream("(1 2 3)")());
-        Pout<< "reading : dl:" << dl << endl;
-    }
+    testField[0] = testField2[0] = vector(1.0, 4.5, 6.3);
+    testField[1] = testField2[1] = vector(5.2, 2.3, 3.5);
+    testField[2] = testField2[2] = vector(7.5, 4.7, 7.7);
+    testField[3] = testField2[3] = vector(2.8, 8.2, 2.3);
+    testField[4] = testField2[4] = vector(6.1, 1.7, 8.8);
 
-    {
-        labelField lf(3);
-        lf[0] = 1;
-        lf[1] = 2;
-        lf[2] = 3;
-        DynamicField<label> dl;
-        dl = lf;
-        Pout<< "assigning from labelField : dl:" << dl << endl;
-    }
+    Info << "testField:" << testField << endl;
 
-    {
-        labelField lf(3);
-        lf[0] = 1;
-        lf[1] = 2;
-        lf[2] = 3;
-        DynamicField<label> dl(lf);
-        Pout<< "constructing from labelField dl:" << dl << endl;
-    }
+    testField.append(vector(0.5, 4.8, 6.2));
 
+    Info << "testField after appending:" << testField << endl;
 
-    Info<< "\nEnd\n";
+    testField.append(vector(2.7, 2.3, 6.1));
+
+    Info << "testField after appending:" << testField << endl;
+
+    vector elem = testField.remove();
+
+    Info << "removed element:" << elem << endl;
+    Info << "testField:" << testField << endl;
+
+    testField.append(vector(3.0, 1.3, 9.2));
+
+    Info << "testField:" << testField << endl;
+
+    testField.setSize(10, vector(1.5, 0.6, -1.0));
+
+    Info << "testField after setSize:" << testField << endl;
+
+    testField.append(testField2);
+
+    Info << "testField after appending testField2:" << testField << endl;
+
+    testField = testField2;
+
+    Info << "testField after assignment:" << testField << endl;
+
+    testField += testField2;
+
+    Info << "testField after field algebra:" << testField << endl;
+
+    testField.clear();
+
+    testField.append(vector(3.0, 1.3, 9.2));
+
+    Info << "testField after clear and append:" << testField << endl;
+
+    testField.clearStorage();
+
+    Info << "testField after clearStorage:" << testField << endl;
 
     return 0;
 }

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -934,18 +934,31 @@ int main(int argc, char *argv[])
                 {
                     char* linePtr = readline("readline>");
 
-                    rawLine = string(linePtr);
-
-                    if (*linePtr)
+                    if (linePtr)
                     {
-                        add_history(linePtr);
-                        write_history(historyFile);
-                    }
+                        rawLine = string(linePtr);
 
-                    free(linePtr);   // readline uses malloc, not new.
+                        if (*linePtr)
+                        {
+                            add_history(linePtr);
+                            write_history(historyFile);
+                        }
+
+                        free(linePtr);   // readline uses malloc, not new.
+                    }
+                    else
+                    {
+                        break;
+                    }
                 }
 #               else
                 {
+                    if (!std::cin.good())
+                    {
+                        Info<< "End of cin" << endl;
+                        // No error.
+                        break;
+                    }
                     Info<< "Command>" << flush;
                     std::getline(std::cin, rawLine);
                 }

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -34,8 +34,12 @@ Description
     - any face inbetween differing cellZones (-cellZones)
 
     Output is:
-    - volScalarField with regions as different scalars (-detectOnly) or
-    - mesh with multiple regions or
+    - volScalarField with regions as different scalars (-detectOnly)
+            or
+    - mesh with multiple regions and directMapped patches. These patches
+      either cover the whole interface between two region (default) or
+      only part according to faceZones (-useFaceZones)
+            or
     - mesh with cells put into cellZones (-makeCellZones)
 
     Note:
@@ -495,18 +499,70 @@ labelList getNonRegionCells(const labelList& cellRegion, const label regionI)
 }
 
 
-// Get per region-region interface the sizes. If sumParallel sums sizes.
+void addToInterface
+(
+    const polyMesh& mesh,
+    const label zoneID,
+    const label ownRegion,
+    const label neiRegion,
+    EdgeMap<Map<label> >& regionsToSize
+)
+{
+    edge interface
+    (
+        min(ownRegion, neiRegion),
+        max(ownRegion, neiRegion)
+    );
+
+    EdgeMap<Map<label> >::iterator iter = regionsToSize.find
+    (
+        interface
+    );
+
+    if (iter != regionsToSize.end())
+    {
+        // Check if zone present
+        Map<label>::iterator zoneFnd = iter().find(zoneID);
+        if (zoneFnd != iter().end())
+        {
+            // Found zone. Increment count.
+            zoneFnd()++;
+        }
+        else
+        {
+            // New or no zone. Insert with count 1.
+            iter().insert(zoneID, 1);
+        }
+    }
+    else
+    {
+        // Create new interface of size 1.
+        Map<label> zoneToSize;
+        zoneToSize.insert(zoneID, 1);
+        regionsToSize.insert(interface, zoneToSize);
+    }
+}
+
+
+// Get region-region interface name and sizes.
 // Returns interfaces as straight list for looping in identical order.
 void getInterfaceSizes
 (
     const polyMesh& mesh,
+    const bool useFaceZones,
     const labelList& cellRegion,
-    const bool sumParallel,
+    const wordList& regionNames,
 
     edgeList& interfaces,
-    EdgeMap<label>& interfaceSizes
+    List<Pair<word> >& interfaceNames,
+    labelList& interfaceSizes,
+    labelList& faceToInterface
 )
 {
+    // From region-region to faceZone (or -1) to number of faces.
+
+    EdgeMap<Map<label> > regionsToSize;
+
 
     // Internal faces
     // ~~~~~~~~~~~~~~
@@ -518,22 +574,14 @@ void getInterfaceSizes
 
         if (ownRegion != neiRegion)
         {
-            edge interface
+            addToInterface
             (
-                min(ownRegion, neiRegion),
-                max(ownRegion, neiRegion)
+                mesh,
+                (useFaceZones ? mesh.faceZones().whichZone(faceI) : -1),
+                ownRegion,
+                neiRegion,
+                regionsToSize
             );
-
-            EdgeMap<label>::iterator iter = interfaceSizes.find(interface);
-
-            if (iter != interfaceSizes.end())
-            {
-                iter()++;
-            }
-            else
-            {
-                interfaceSizes.insert(interface, 1);
-            }
         }
     }
 
@@ -558,27 +606,19 @@ void getInterfaceSizes
 
         if (ownRegion != neiRegion)
         {
-            edge interface
+            addToInterface
             (
-                min(ownRegion, neiRegion),
-                max(ownRegion, neiRegion)
+                mesh,
+                (useFaceZones ? mesh.faceZones().whichZone(faceI) : -1),
+                ownRegion,
+                neiRegion,
+                regionsToSize
             );
-
-            EdgeMap<label>::iterator iter = interfaceSizes.find(interface);
-
-            if (iter != interfaceSizes.end())
-            {
-                iter()++;
-            }
-            else
-            {
-                interfaceSizes.insert(interface, 1);
-            }
         }
     }
 
 
-    if (sumParallel && Pstream::parRun())
+    if (Pstream::parRun())
     {
         if (Pstream::master())
         {
@@ -592,57 +632,169 @@ void getInterfaceSizes
             {
                 IPstream fromSlave(Pstream::blocking, slave);
 
-                EdgeMap<label> slaveSizes(fromSlave);
+                EdgeMap<Map<label> > slaveSizes(fromSlave);
 
-                forAllConstIter(EdgeMap<label>, slaveSizes, slaveIter)
+                forAllConstIter(EdgeMap<Map<label> >, slaveSizes, slaveIter)
                 {
-                    EdgeMap<label>::iterator masterIter =
-                        interfaceSizes.find(slaveIter.key());
+                    EdgeMap<Map<label> >::iterator masterIter =
+                        regionsToSize.find(slaveIter.key());
 
-                    if (masterIter != interfaceSizes.end())
+                    if (masterIter != regionsToSize.end())
                     {
-                        masterIter() += slaveIter();
+                        // Same inter-region
+                        const Map<label>& slaveInfo = slaveIter();
+                        Map<label>& masterInfo = masterIter();
+
+                        forAllConstIter(Map<label>, slaveInfo, iter)
+                        {
+                            label zoneID = iter.key();
+                            label slaveSize = iter();
+
+                            Map<label>::iterator zoneFnd = masterInfo.find
+                            (
+                                zoneID
+                            );
+                            if (zoneFnd != masterInfo.end())
+                            {
+                                zoneFnd() += slaveSize;
+                            }
+                            else
+                            {
+                                masterInfo.insert(zoneID, slaveSize);
+                            }
+                        }
                     }
                     else
                     {
-                        interfaceSizes.insert(slaveIter.key(), slaveIter());
+                        regionsToSize.insert(slaveIter.key(), slaveIter());
                     }
                 }
             }
-
-            // Distribute
-            for
-            (
-                int slave=Pstream::firstSlave();
-                slave<=Pstream::lastSlave();
-                slave++
-            )
-            {
-                // Receive the edges using shared points from the slave.
-                OPstream toSlave(Pstream::blocking, slave);
-                toSlave << interfaceSizes;
-            }
         }
         else
         {
             // Send to master
             {
                 OPstream toMaster(Pstream::blocking, Pstream::masterNo());
-                toMaster << interfaceSizes;
-            }
-            // Receive from master
-            {
-                IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
-                fromMaster >> interfaceSizes;
+                toMaster << regionsToSize;
             }
         }
     }
 
-    // Make sure all processors have interfaces in same order
-    interfaces = interfaceSizes.toc();
-    if (sumParallel)
+    // Rework 
+
+    Pstream::scatter(regionsToSize);
+
+
+
+    // Now we have the global sizes of all inter-regions.
+    // Invert this on master and distribute.
+    label nInterfaces = 0;
+    forAllConstIter(EdgeMap<Map<label> >, regionsToSize, iter)
     {
-        Pstream::scatter(interfaces);
+        const Map<label>& info = iter();
+        nInterfaces += info.size();
+    }
+
+    interfaces.setSize(nInterfaces);
+    interfaceNames.setSize(nInterfaces);
+    interfaceSizes.setSize(nInterfaces);
+    EdgeMap<Map<label> > regionsToInterface(nInterfaces);
+
+    nInterfaces = 0;
+    forAllConstIter(EdgeMap<Map<label> >, regionsToSize, iter)
+    {
+        const edge& e = iter.key();
+        const word& name0 = regionNames[e[0]];
+        const word& name1 = regionNames[e[1]];
+
+        const Map<label>& info = iter();
+        forAllConstIter(Map<label>, info, infoIter)
+        {
+            interfaces[nInterfaces] = iter.key();
+            label zoneID = infoIter.key();
+            if (zoneID == -1)
+            {
+                interfaceNames[nInterfaces] = Pair<word>
+                (
+                    name0 + "_to_" + name1,
+                    name1 + "_to_" + name0
+                );
+            }
+            else
+            {
+                const word& zoneName = mesh.faceZones()[zoneID].name();
+                interfaceNames[nInterfaces] = Pair<word>
+                (
+                    zoneName + "_" + name0 + "_to_" + name1,
+                    zoneName + "_" + name1 + "_to_" + name0
+                );
+            }
+            interfaceSizes[nInterfaces] = infoIter();
+
+            Map<label> zoneAndInterface;
+            zoneAndInterface.insert(zoneID, nInterfaces);
+            regionsToInterface.insert(e, zoneAndInterface);
+
+            nInterfaces++;
+        }
+    }
+
+
+    // Now all processor have consistent interface information
+
+    Pstream::scatter(interfaces);
+    Pstream::scatter(interfaceNames);
+    Pstream::scatter(interfaceSizes);
+    Pstream::scatter(regionsToInterface);
+
+    // Mark all inter-region faces.
+    faceToInterface.setSize(mesh.nFaces(), -1);
+
+    forAll(mesh.faceNeighbour(), faceI)
+    {
+        label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
+        label neiRegion = cellRegion[mesh.faceNeighbour()[faceI]];
+
+        if (ownRegion != neiRegion)
+        {
+            label zoneID = -1;
+            if (useFaceZones)
+            {
+                zoneID = mesh.faceZones().whichZone(faceI);
+            }
+
+            edge interface
+            (
+                min(ownRegion, neiRegion),
+                max(ownRegion, neiRegion)
+            );
+
+            faceToInterface[faceI] = regionsToInterface[interface][zoneID];
+        }
+    }
+    forAll(coupledRegion, i)
+    {
+        label faceI = i+mesh.nInternalFaces();
+        label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
+        label neiRegion = coupledRegion[i];
+
+        if (ownRegion != neiRegion)
+        {
+            label zoneID = -1;
+            if (useFaceZones)
+            {
+                zoneID = mesh.faceZones().whichZone(faceI);
+            }
+
+            edge interface
+            (
+                min(ownRegion, neiRegion),
+                max(ownRegion, neiRegion)
+            );
+
+            faceToInterface[faceI] = regionsToInterface[interface][zoneID];
+        }
     }
 }
 
@@ -650,11 +802,15 @@ void getInterfaceSizes
 // Create mesh for region.
 autoPtr<mapPolyMesh> createRegionMesh
 (
-    const labelList& cellRegion,
-    const EdgeMap<label>& interfaceToPatch,
     const fvMesh& mesh,
+    // Region info
+    const labelList& cellRegion,
     const label regionI,
     const word& regionName,
+    // Interface info
+    const labelList& interfacePatches,
+    const labelList& faceToInterface,
+
     autoPtr<fvMesh>& newMesh
 )
 {
@@ -739,6 +895,7 @@ autoPtr<mapPolyMesh> createRegionMesh
     forAll(exposedFaces, i)
     {
         label faceI = exposedFaces[i];
+        label interfaceI = faceToInterface[faceI];
 
         label ownRegion = cellRegion[mesh.faceOwner()[faceI]];
         label neiRegion = -1;
@@ -752,6 +909,10 @@ autoPtr<mapPolyMesh> createRegionMesh
             neiRegion = coupledRegion[faceI-mesh.nInternalFaces()];
         }
 
+
+        // Check which side is being kept - determines which of the two
+        // patches will be used.
+
         label otherRegion = -1;
 
         if (ownRegion == regionI && neiRegion != regionI)
@@ -773,19 +934,14 @@ autoPtr<mapPolyMesh> createRegionMesh
                 << exit(FatalError);
         }
 
-        if (otherRegion != -1)
+        // Find the patch.
+        if (regionI < otherRegion)
         {
-            edge interface(regionI, otherRegion);
-
-            // Find the patch.
-            if (regionI < otherRegion)
-            {
-                exposedPatchIDs[i] = interfaceToPatch[interface];
-            }
-            else
-            {
-                exposedPatchIDs[i] = interfaceToPatch[interface]+1;
-            }
+            exposedPatchIDs[i] = interfacePatches[interfaceI];
+        }
+        else
+        {
+            exposedPatchIDs[i] = interfacePatches[interfaceI]+1;
         }
     }
 
@@ -821,7 +977,8 @@ void createAndWriteRegion
     const fvMesh& mesh,
     const labelList& cellRegion,
     const wordList& regionNames,
-    const EdgeMap<label>& interfaceToPatch,
+    const labelList& faceToInterface,
+    const labelList& interfacePatches,
     const label regionI,
     const word& newMeshInstance
 )
@@ -832,21 +989,22 @@ void createAndWriteRegion
     autoPtr<fvMesh> newMesh;
     autoPtr<mapPolyMesh> map = createRegionMesh
     (
-        cellRegion,
-        interfaceToPatch,
         mesh,
+        cellRegion,
         regionI,
         regionNames[regionI],
+        interfacePatches,
+        faceToInterface,
         newMesh
     );
 
 
     // Make map of all added patches
-    labelHashSet addedPatches(2*interfaceToPatch.size());
-    forAllConstIter(EdgeMap<label>, interfaceToPatch, iter)
+    labelHashSet addedPatches(2*interfacePatches.size());
+    forAll(interfacePatches, interfaceI)
     {
-        addedPatches.insert(iter());
-        addedPatches.insert(iter()+1);
+        addedPatches.insert(interfacePatches[interfaceI]);
+        addedPatches.insert(interfacePatches[interfaceI]+1);
     }
 
     Info<< "Mapping fields" << endl;
@@ -1074,70 +1232,67 @@ void createAndWriteRegion
 // First one is for minimumregion to maximumregion.
 // Note that patches get created in same order on all processors (if parallel)
 // since looping over synchronised 'interfaces'.
-EdgeMap<label> addRegionPatches
+labelList addRegionPatches
 (
     fvMesh& mesh,
-    const labelList& cellRegion,
-    const label nCellRegions,
+    const wordList& regionNames,
     const edgeList& interfaces,
-    const EdgeMap<label>& interfaceSizes,
-    const wordList& regionNames
+    const List<Pair<word> >& interfaceNames
 )
 {
-    // Check that all patches are present in same order.
-    mesh.boundaryMesh().checkParallelSync(true);
-
     Info<< nl << "Adding patches" << nl << endl;
 
-    EdgeMap<label> interfaceToPatch(nCellRegions);
+    labelList interfacePatches(interfaces.size());
 
     forAll(interfaces, interI)
     {
         const edge& e = interfaces[interI];
+        const Pair<word>& names = interfaceNames[interI];
 
-        if (interfaceSizes[e] > 0)
-        {
-            const word inter1 = regionNames[e[0]] + "_to_" + regionNames[e[1]];
-            const word inter2 = regionNames[e[1]] + "_to_" + regionNames[e[0]];
+        //Info<< "For interface " << interI
+        //    << " between regions " << e
+        //    << " trying to add patches " << names << endl;
 
-            directMappedWallPolyPatch patch1
-            (
-                inter1,
-                0,                  // overridden
-                0,                  // overridden
-                0,                  // overridden
-                regionNames[e[1]],  // sampleRegion
-                directMappedPatchBase::NEARESTPATCHFACE,
-                inter2,             // samplePatch
-                point::zero,        // offset
-                mesh.boundaryMesh()
-            );
 
-            label patchI = addPatch(mesh, patch1);
+        directMappedWallPolyPatch patch1
+        (
+            names[0],
+            0,                  // overridden
+            0,                  // overridden
+            0,                  // overridden
+            regionNames[e[1]],  // sampleRegion
+            directMappedPatchBase::NEARESTPATCHFACE,
+            names[1],           // samplePatch
+            point::zero,        // offset
+            mesh.boundaryMesh()
+        );
 
-            directMappedWallPolyPatch patch2
-            (
-                inter2,
-                0,
-                0,
-                0,
-                regionNames[e[0]],  // sampleRegion
-                directMappedPatchBase::NEARESTPATCHFACE,
-                inter1,
-                point::zero,        // offset
-                mesh.boundaryMesh()
-            );
-            addPatch(mesh, patch2);
+        interfacePatches[interI] = addPatch(mesh, patch1);
 
-            Info<< "For interface between region " << e[0]
-                << " and " << e[1] << " added patch " << patchI
-                << " " << mesh.boundaryMesh()[patchI].name()
-                << endl;
+        directMappedWallPolyPatch patch2
+        (
+            names[1],
+            0,
+            0,
+            0,
+            regionNames[e[0]],  // sampleRegion
+            directMappedPatchBase::NEARESTPATCHFACE,
+            names[0],
+            point::zero,        // offset
+            mesh.boundaryMesh()
+        );
+        addPatch(mesh, patch2);
 
-            interfaceToPatch.insert(e, patchI);
-        }
+        Info<< "For interface between region " << regionNames[e[0]]
+            << " and " << regionNames[e[1]] << " added patches" << endl
+            << "    " << interfacePatches[interI]
+            << "\t" << mesh.boundaryMesh()[interfacePatches[interI]].name()
+            << endl
+            << "    " << interfacePatches[interI]+1
+            << "\t" << mesh.boundaryMesh()[interfacePatches[interI]+1].name()
+            << endl;
     }
-    return interfaceToPatch;
+    return interfacePatches;
 }
 
 
@@ -1195,76 +1350,6 @@ label findCorrespondingRegion
 }
 
 
-//// Checks if cellZone has corresponding cellRegion.
-//label findCorrespondingRegion
-//(
-//    const cellZoneMesh& cellZones,
-//    const labelList& existingZoneID,    // per cell the (unique) zoneID
-//    const labelList& cellRegion,
-//    const label nCellRegions,
-//    const label zoneI
-//)
-//{
-//    // Region corresponding to zone. Start off with special value: no
-//    // corresponding region.
-//    label regionI = labelMax;
-//
-//    const cellZone& cz = cellZones[zoneI];
-//
-//    if (cz.empty())
-//    {
-//        // My local portion is empty. Maps to any empty cellZone. Mark with
-//        // special value which can get overwritten by other processors.
-//        regionI = -1;
-//    }
-//    else
-//    {
-//        regionI = cellRegion[cz[0]];
-//
-//        forAll(cz, i)
-//        {
-//            if (cellRegion[cz[i]] != regionI)
-//            {
-//                regionI = labelMax;
-//                break;
-//            }
-//        }
-//    }
-//
-//    // Determine same zone over all processors.
-//    reduce(regionI, maxOp<label>());
-//
-//
-//    // 2. All of region present?
-//
-//    if (regionI == labelMax)
-//    {
-//        regionI = -1;
-//    }
-//    else if (regionI != -1)
-//    {
-//        forAll(cellRegion, cellI)
-//        {
-//            if
-//            (
-//                cellRegion[cellI] == regionI
-//             && existingZoneID[cellI] != zoneI
-//            )
-//            {
-//                // cellI in regionI but not in zoneI
-//                regionI = -1;
-//                break;
-//            }
-//        }
-//        // If one in error, all should be in error. Note that branch
-//        // gets taken on all procs.
-//        reduce(regionI, minOp<label>());
-//    }
-//
-//    return regionI;
-//}
-
-
 // Get zone per cell
 // - non-unique zoning
 // - coupled zones
@@ -1484,6 +1569,7 @@ void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
 }
 
 
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -1541,6 +1627,11 @@ int main(int argc, char *argv[])
         "sloppyCellZones",
         "try to match heuristically regions to existing cell zones"
     );
+    argList::addBoolOption
+    (
+        "useFaceZones",
+        "use faceZones to patch inter-region faces instead of single patch"
+    );
 
     #include "setRootCase.H"
     #include "createTime.H"
@@ -1564,6 +1655,8 @@ int main(int argc, char *argv[])
     const bool overwrite        = args.optionFound("overwrite");
     const bool detectOnly       = args.optionFound("detectOnly");
     const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
+    const bool useFaceZones     = args.optionFound("useFaceZones");
+
     if
     (
         (useCellZonesOnly || useCellZonesFile)
@@ -1579,6 +1672,20 @@ int main(int argc, char *argv[])
     }
 
 
+    if (useFaceZones)
+    {
+        Info<< "Using current faceZones to divide inter-region interfaces"
+            << " into multiple patches."
+            << nl << endl;
+    }
+    else
+    {
+        Info<< "Creating single patch per inter-region interface."
+            << nl << endl;
+    }
+
+
+
     if (insidePoint && largestOnly)
     {
         FatalErrorIn(args.executable())
@@ -1768,6 +1875,7 @@ int main(int argc, char *argv[])
     writeCellToRegion(mesh, cellRegion);
 
 
+
     // Sizes per region
     // ~~~~~~~~~~~~~~~~
 
@@ -1805,34 +1913,48 @@ int main(int argc, char *argv[])
 
 
 
-    // Since we're going to mess with patches make sure all non-processor ones
-    // are on all processors.
+    // Since we're going to mess with patches and zones make sure all
+    // is synchronised
     mesh.boundaryMesh().checkParallelSync(true);
+    mesh.faceZones().checkParallelSync(true);
 
 
-    // Sizes of interface between regions. From pair of regions to number of
-    // faces.
+    // Interfaces
+    // ----------
+    // per interface:
+    // - the two regions on either side
+    // - the name
+    // - the (global) size
     edgeList interfaces;
-    EdgeMap<label> interfaceSizes;
+    List<Pair<word> > interfaceNames;
+    labelList interfaceSizes;
+    // per face the interface
+    labelList faceToInterface;
+
     getInterfaceSizes
     (
         mesh,
+        useFaceZones,
         cellRegion,
-        true,      // sum in parallel?
+        regionNames,
 
         interfaces,
-        interfaceSizes
+        interfaceNames,
+        interfaceSizes,
+        faceToInterface
     );
 
-    Info<< "Sizes inbetween regions:" << nl << nl
-        << "Region\tRegion\tFaces" << nl
-        << "------\t------\t-----" << endl;
+    Info<< "Sizes of interfaces between regions:" << nl << nl
+        << "Interface\tRegion\tRegion\tFaces" << nl
+        << "---------\t------\t------\t-----" << endl;
 
     forAll(interfaces, interI)
     {
         const edge& e = interfaces[interI];
 
-        Info<< e[0] << '\t' << e[1] << '\t' << interfaceSizes[e] << nl;
+        Info<< interI
+            << "\t\t" << e[0] << '\t' << e[1]
+            << '\t' << interfaceSizes[interI] << nl;
     }
     Info<< endl;
 
@@ -1982,16 +2104,14 @@ int main(int argc, char *argv[])
         // Add all possible patches. Empty ones get filtered later on.
         Info<< nl << "Adding patches" << nl << endl;
 
-        EdgeMap<label> interfaceToPatch
+        labelList interfacePatches
         (
             addRegionPatches
             (
                 mesh,
-                cellRegion,
-                nCellRegions,
+                regionNames,
                 interfaces,
-                interfaceSizes,
-                regionNames
+                interfaceNames
             )
         );
 
@@ -2041,7 +2161,8 @@ int main(int argc, char *argv[])
                 mesh,
                 cellRegion,
                 regionNames,
-                interfaceToPatch,
+                faceToInterface,
+                interfacePatches,
                 regionI,
                 (overwrite ? oldInstance : runTime.timeName())
             );
@@ -2059,7 +2180,8 @@ int main(int argc, char *argv[])
                 mesh,
                 cellRegion,
                 regionNames,
-                interfaceToPatch,
+                faceToInterface,
+                interfacePatches,
                 regionI,
                 (overwrite ? oldInstance : runTime.timeName())
             );
@@ -2078,7 +2200,8 @@ int main(int argc, char *argv[])
                     mesh,
                     cellRegion,
                     regionNames,
-                    interfaceToPatch,
+                    faceToInterface,
+                    interfacePatches,
                     regionI,
                     (overwrite ? oldInstance : runTime.timeName())
                 );

applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C

@@ -68,13 +68,13 @@ namespace Foam
 
             functionObjectList fol(runTime, dict);
             fol.start();
-            fol.execute();
+            fol.execute(true);  // override outputControl - force writing
         }
         else
         {
             functionObjectList fol(runTime);
             fol.start();
-            fol.execute();
+            fol.execute(true);  // override outputControl - force writing
         }
     }
 }

bin/mpirunDebug

@@ -162,12 +162,14 @@ do
     procCmdFile="$PWD/processor${proc}.sh"
     procLog="processor${proc}.log"
     geom="-geometry 120x20+$xpos+$ypos"
-    node=""
 
-    if [ "$WM_MPLIB" = OPENMPI ]
-    then
+    case "$WM_MPLIB" in
+    *OPENMPI)
         node="-np 1 "
-    fi
+        ;;
+    *)
+        node=""
+    esac
 
     echo "#!/bin/sh" > $procCmdFile
     case "$method" in

bin/tools/CleanFunctions

@@ -65,6 +65,7 @@ cleanCase()
     rm -rf processor* > /dev/null 2>&1
     rm -rf probes* > /dev/null 2>&1
     rm -rf forces* > /dev/null 2>&1
+    rm -rf sets > /dev/null 2>&1
     rm -rf system/machines > /dev/null 2>&1
     (cd constant/polyMesh && \
         rm -rf \

doc/changes/splitCyclic.txt

@@ -59,6 +59,14 @@ It will check if anything needs to be converted, backup the current file to .old
 and split any cyclic patchFields into two entries.
 
 
+Mesh converters
+---------------
+Most mesh formats use cyclics in a single patch (i.e. the old way).
+The converters have been adapted to use the patch 'oldCyclic' instead of
+'cyclic'. oldCyclic uses the 17x automatic ordering but writes 'type cyclic'
+so afterwards foamUpgradeCyclics can be run to upgrade.
+
+
 decomposePar
 ------------
 Decomposes cyclics into processorCyclic:

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -37,8 +37,11 @@ namespace Foam
     const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
 
     const scalar
-        SIBS::safe1 = 0.25, SIBS::safe2 = 0.7,
-        SIBS::redMax = 1.0e-5, SIBS::redMin = 0.0, SIBS::scaleMX = 0.1;
+        SIBS::safe1 = 0.25,
+        SIBS::safe2 = 0.7,
+        SIBS::redMax = 1.0e-5,
+        SIBS::redMin = 0.7,
+        SIBS::scaleMX = 0.1;
 };
 
 

src/OSspecific/POSIX/fileMonitor.C

@@ -170,7 +170,9 @@ namespace Foam
                         << "inotify instances" << endl
                         << "    (/proc/sys/fs/inotify/max_user_instances"
                         << " on Linux)" << endl
-                        << "    or switch off runTimeModifiable." << endl
+                        << "    , switch off runTimeModifiable." << endl
+                        << "    or compile this file with FOAM_USE_STAT to use"
+                        << " time stamps instead of inotify." << endl
                         << "    Continuing without additional file monitoring."
                         << endl;
                 }

src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H

@@ -28,8 +28,8 @@ License
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList()
 :
-    List<T>(SizeInc),
-    capacity_(SizeInc)
+    List<T>(0),
+    capacity_(0)
 {
     List<T>::size(0);
 }

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C

@@ -128,7 +128,10 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::start()
 
 
 template<class OutputFilter>
-bool Foam::OutputFilterFunctionObject<OutputFilter>::execute()
+bool Foam::OutputFilterFunctionObject<OutputFilter>::execute
+(
+    const bool forceWrite
+)
 {
     if (enabled_)
     {
@@ -139,7 +142,7 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::execute()
 
         ptr_->execute();
 
-        if (outputControl_.output())
+        if (forceWrite || outputControl_.output())
         {
             ptr_->write();
         }

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H

@@ -183,7 +183,7 @@ public:
             virtual bool start();
 
             //- Called at each ++ or += of the time-loop
-            virtual bool execute();
+            virtual bool execute(const bool forceWrite);
 
             //- Called when Time::run() determines that the time-loop exits
             virtual bool end();

src/OpenFOAM/db/functionObjects/functionObject/functionObject.C

@@ -112,7 +112,7 @@ const Foam::word& Foam::functionObject::name() const
 
 bool Foam::functionObject::end()
 {
-    return execute();
+    return execute(false);
 }
 
 

src/OpenFOAM/db/functionObjects/functionObject/functionObject.H

@@ -147,8 +147,9 @@ public:
         //- Called at the start of the time-loop
         virtual bool start() = 0;
 
-        //- Called at each ++ or += of the time-loop
-        virtual bool execute() = 0;
+        //- Called at each ++ or += of the time-loop. forceWrite overrides the
+        //  outputControl behaviour.
+        virtual bool execute(const bool forceWrite) = 0;
 
         //- Called when Time::run() determines that the time-loop exits.
         //  By default it simply calls execute().

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.C

@@ -144,7 +144,7 @@ bool Foam::functionObjectList::start()
 }
 
 
-bool Foam::functionObjectList::execute()
+bool Foam::functionObjectList::execute(const bool forceWrite)
 {
     bool ok = true;
 
@@ -157,7 +157,7 @@ bool Foam::functionObjectList::execute()
 
         forAll(*this, objectI)
         {
-            ok = operator[](objectI).execute() && ok;
+            ok = operator[](objectI).execute(forceWrite) && ok;
         }
     }
 

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.H

@@ -153,8 +153,10 @@ public:
         //- Called at the start of the time-loop
         virtual bool start();
 
-        //- Called at each ++ or += of the time-loop
-        virtual bool execute();
+        //- Called at each ++ or += of the time-loop. forceWrite overrides
+        //  the usual outputControl behaviour and forces writing always
+        //  (used in postprocessing mode)
+        virtual bool execute(const bool forceWrite = false);
 
         //- Called when Time::run() determines that the time-loop exits
         virtual bool end();

src/OpenFOAM/fields/Fields/DynamicField/DynamicField.C

@@ -25,87 +25,53 @@ License
 
 #include "DynamicField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Static Members  * * * * * * * * * * * * * * //
-
-template<class Type>
-const char* const DynamicField<Type>::typeName("DynamicField");
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class Type>
-DynamicField<Type>::DynamicField(Istream& is)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField(Istream& is)
 :
-    Field<Type>(is),
-    capacity_(Field<Type>::size())
+    Field<T>(is),
+    capacity_(Field<T>::size())
 {}
 
 
-template<class Type>
-tmp<DynamicField<Type> > DynamicField<Type>::clone() const
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+Foam::tmp<Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv> >
+Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::clone() const
 {
-    return tmp<DynamicField<Type> >(new DynamicField<Type>(*this));
+    return tmp<DynamicField<T, SizeInc, SizeMult, SizeDiv> >
+    (
+        new DynamicField<T, SizeInc, SizeMult, SizeDiv>(*this)
+    );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-void DynamicField<Type>::setSize(const label nElem)
-{
-    // allocate more capacity?
-    if (nElem > capacity_)
-    {
-        capacity_ = max(nElem, label(1 + capacity_*2));
-
-        Field<Type>::setSize(capacity_);
-    }
-
-    // adjust addressed size
-    Field<Type>::size(nElem);
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * IOstream Operator * * * * * * * * * * * * * //
 
-template<class Type>
-Ostream& operator<<(Ostream& os, const DynamicField<Type>& f)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const DynamicField<T, SizeInc, SizeMult, SizeDiv>& lst
+)
 {
-    os << static_cast<const Field<Type>&>(f);
+    os << static_cast<const Field<T>&>(lst);
     return os;
 }
 
 
-template<class Type>
-Ostream& operator<<(Ostream& os, const tmp<DynamicField<Type> >& tf)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    DynamicField<T, SizeInc, SizeMult, SizeDiv>& lst
+)
 {
-    os << tf();
-    tf.clear();
-    return os;
-}
-
-
-template<class Type>
-Istream& operator>>(Istream& is, DynamicField<Type>& lst)
-{
-    is >> static_cast<Field<Type>&>(lst);
-    lst.capacity_ = lst.Field<Type>::size();
+    is >> static_cast<Field<T>&>(lst);
+    lst.capacity_ = lst.Field<T>::size();
 
     return is;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-
 // ************************************************************************* //

src/OpenFOAM/fields/Fields/DynamicField/DynamicField.H

@@ -25,9 +25,10 @@ Class
     Foam::DynamicField
 
 Description
-    Dynamically sized Field. WIP.
+    Dynamically sized Field.
 
 SourceFiles
+    DynamicFieldI.H
     DynamicField.C
 
 \*---------------------------------------------------------------------------*/
@@ -44,89 +45,78 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class Type>
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 class DynamicField;
 
-template<class Type>
-Ostream& operator<<(Ostream&, const DynamicField<Type>&);
-
-template<class Type>
-Ostream& operator<<(Ostream&, const tmp<DynamicField<Type> >&);
-
-template<class Type>
-Istream& operator>>(Istream&, DynamicField<Type>&);
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+Ostream& operator<<
+(
+    Ostream&,
+    const DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+);
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+Istream& operator>>
+(
+    Istream&,
+    DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+);
 
 
 /*---------------------------------------------------------------------------*\
                            Class DynamicField Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class Type>
+template<class T, unsigned SizeInc=0, unsigned SizeMult=2, unsigned SizeDiv=1>
 class DynamicField
 :
-    public Field<Type>
+    public Field<T>
 {
-
     // Private data
 
         //- The capacity (allocated size) of the underlying field.
         label capacity_;
 
 
-        //- Construct given size and initial value
-        DynamicField(const label, const Type&);
-
-        //- Construct as copy of tmp<DynamicField>
-#       ifdef ConstructFromTmp
-        DynamicField(const tmp<DynamicField<Type> >&);
-#       endif
-
-        //- Construct from a dictionary entry
-        DynamicField(const word&, const dictionary&, const label);
-
 public:
 
-    // Static data members
-
-        static const char* const typeName;
-
-
     // Static Member Functions
 
         //- Return a null field
-        inline static const DynamicField<Type>& null()
+        inline static const DynamicField<T, SizeInc, SizeMult, SizeDiv>& null()
         {
-            return *reinterpret_cast< DynamicField<Type>* >(0);
+            return *reinterpret_cast
+            <
+                DynamicField<T, SizeInc, SizeMult, SizeDiv>*
+            >(0);
         }
 
 
     // Constructors
 
         //- Construct null
-        //  Used for temporary fields which are initialised after construction
-        DynamicField();
+        inline DynamicField();
 
-        //- Construct given size
-        //  Used for temporary fields which are initialised after construction
+        //- Construct given size.
         explicit inline DynamicField(const label);
 
-        //- Construct as copy of a UList\<Type\>
-        explicit inline DynamicField(const UList<Type>&);
+        //- Construct from UList. Size set to UList size.
+        //  Also constructs from DynamicField with different sizing parameters.
+        explicit inline DynamicField(const UList<T>&);
 
-        //- Construct by transferring the List contents
-        explicit inline DynamicField(const Xfer<List<Type> >&);
+        //- Construct by transferring the parameter contents
+        explicit inline DynamicField(const Xfer<List<T> >&);
 
         //- Construct by 1 to 1 mapping from the given field
         inline DynamicField
         (
-            const UList<Type>& mapF,
+            const UList<T>& mapF,
             const labelList& mapAddressing
         );
 
         //- Construct by interpolative mapping from the given field
         inline DynamicField
         (
-            const UList<Type>& mapF,
+            const UList<T>& mapF,
             const labelListList& mapAddressing,
             const scalarListList& weights
         );
@@ -134,59 +124,129 @@ public:
         //- Construct by mapping from the given field
         inline DynamicField
         (
-            const UList<Type>& mapF,
+            const UList<T>& mapF,
             const FieldMapper& map
         );
 
-        //- Construct as copy
-        inline DynamicField(const DynamicField<Type>&);
-
-        //- Construct as copy or re-use as specified.
-        inline DynamicField(DynamicField<Type>&, bool reUse);
+        //- Construct copy
+        inline DynamicField(const DynamicField<T, SizeInc, SizeMult, SizeDiv>&);
 
         //- Construct by transferring the Field contents
-        inline DynamicField(const Xfer<DynamicField<Type> >&);
+        inline DynamicField
+        (
+            const Xfer<DynamicField<T, SizeInc, SizeMult, SizeDiv> >&
+        );
 
-        //- Construct from Istream
-        inline DynamicField(Istream&);
+        //- Construct from Istream. Size set to size of list read.
+        explicit DynamicField(Istream&);
 
         //- Clone
-        tmp<DynamicField<Type> > clone() const;
+        tmp<DynamicField<T, SizeInc, SizeMult, SizeDiv> > clone() const;
 
 
     // Member Functions
 
-        //- Size of the underlying storage.
-        inline label capacity() const;
+        // Access
 
-        //- Append an element at the end of the list
-        inline void append(const Type&);
+            //- Size of the underlying storage.
+            inline label capacity() const;
 
-        //- Alter the addressed list size.
-        //  New space will be allocated if required.
-        //  Use this to resize the list prior to using the operator[] for
-        //  setting values (as per List usage).
-        void setSize(const label nElem);
+        // Edit
 
-    // Member operators
+            //- Alter the size of the underlying storage.
+            //  The addressed size will be truncated if needed to fit, but will
+            //  remain otherwise untouched.
+            //  Use this or reserve() in combination with append().
+            inline void setCapacity(const label);
 
-        inline void operator=(const DynamicField<Type>&);
-        inline void operator=(const UList<Type>&);
-        inline void operator=(const tmp<DynamicField<Type> >&);
+            //- Alter the addressed list size.
+            //  New space will be allocated if required.
+            //  Use this to resize the list prior to using the operator[] for
+            //  setting values (as per List usage).
+            inline void setSize(const label);
 
-        //- Return element of Field.
-        using Field<Type>::operator[];
+            //- Alter the addressed list size and fill new space with a
+            //  constant.
+            inline void setSize(const label, const T&);
 
-    // IOstream operators
+            //- Alter the addressed list size.
+            //  New space will be allocated if required.
+            //  Use this to resize the list prior to using the operator[] for
+            //  setting values (as per List usage).
+            inline void resize(const label);
 
-        friend Ostream& operator<< <Type>
-        (Ostream&, const DynamicField<Type>&);
+            //- Alter the addressed list size and fill new space with a
+            //  constant.
+            inline void resize(const label, const T&);
 
-        friend Ostream& operator<< <Type>
-        (Ostream&, const tmp<DynamicField<Type> >&);
+            //- Reserve allocation space for at least this size.
+            //  Never shrinks the allocated size, use setCapacity() for that.
+            inline void reserve(const label);
 
-        friend Istream& operator>> <Type>
-        (Istream&, DynamicField<Type>&);
+            //- Clear the addressed list, i.e. set the size to zero.
+            //  Allocated size does not change
+            inline void clear();
+
+            //- Clear the list and delete storage.
+            inline void clearStorage();
+
+            //- Shrink the allocated space to the number of elements used.
+            //  Returns a reference to the DynamicField.
+            inline DynamicField<T, SizeInc, SizeMult, SizeDiv>& shrink();
+
+            //- Transfer contents to the Xfer container as a plain List
+            inline Xfer<List<T> > xfer();
+
+
+        // Member Operators
+
+            //- Append an element at the end of the list
+            inline DynamicField<T, SizeInc, SizeMult, SizeDiv>& append
+            (
+                const T&
+            );
+
+            //- Append a List at the end of this list
+            inline DynamicField<T, SizeInc, SizeMult, SizeDiv>& append
+            (
+                const UList<T>&
+            );
+
+            //- Remove and return the top element
+            inline T remove();
+
+            //- Return non-const access to an element, resizing list if
+            //  necessary
+            inline T& operator()(const label);
+
+            //- Assignment of all addressed entries to the given value
+            inline void operator=(const T&);
+
+            //- Assignment from DynamicField
+            inline void operator=
+            (
+                const DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+            );
+
+            //- Assignment from UList
+            inline void operator=(const UList<T>&);
+
+
+        // IOstream operators
+
+            // Write DynamicField to Ostream.
+            friend Ostream& operator<< <T, SizeInc, SizeMult, SizeDiv>
+            (
+                Ostream&,
+                const DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+            );
+
+            //- Read from Istream, discarding contents of existing DynamicField.
+            friend Istream& operator>> <T, SizeInc, SizeMult, SizeDiv>
+            (
+                Istream&,
+                DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+            );
 };
 
 

src/OpenFOAM/fields/Fields/DynamicField/DynamicFieldI.H

@@ -23,175 +23,441 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "DynamicField.H"
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::DynamicField<Type>::DynamicField()
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField()
 :
-    Field<Type>(),
-    capacity_(0)
+    Field<T>(0),
+    capacity_(Field<T>::size())
 {}
 
 
-template<class Type>
-Foam::DynamicField<Type>::DynamicField(const label size)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
+(
+    const label nElem
+)
 :
-    Field<Type>(size),
-    capacity_(Field<Type>::size())
+    Field<T>(nElem),
+    capacity_(Field<T>::size())
 {
-    Field<Type>::size(0);
+    // we could also enforce SizeInc granularity when (!SizeMult || !SizeDiv)
+    Field<T>::size(0);
 }
 
 
-template<class Type>
-inline Foam::DynamicField<Type>::DynamicField
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
 (
-    const UList<Type>& lst
+    const UList<T>& lst
 )
 :
-    Field<Type>(lst),
-    capacity_(Field<Type>::size())
+    Field<T>(lst),
+    capacity_(Field<T>::size())
 {}
 
 
-template<class Type>
-inline Foam::DynamicField<Type>::DynamicField
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
 (
-    const Xfer<List<Type> >& lst
+    const Xfer<List<T> >& lst
 )
 :
-    Field<Type>(lst),
-    capacity_(Field<Type>::size())
+    Field<T>(lst),
+    capacity_(Field<T>::size())
 {}
 
 
-template<class Type>
-Foam::DynamicField<Type>::DynamicField
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
 (
-    const UList<Type>& mapF,
+    const UList<T>& mapF,
     const labelList& mapAddressing
 )
 :
-    Field<Type>(mapF, mapAddressing),
-    capacity_(Field<Type>::size())
+    Field<T>(mapF, mapAddressing),
+    capacity_(Field<T>::size())
 {}
 
 
-template<class Type>
-Foam::DynamicField<Type>::DynamicField
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
 (
-    const UList<Type>& mapF,
+    const UList<T>& mapF,
     const labelListList& mapAddressing,
     const scalarListList& weights
 )
 :
-    Field<Type>(mapF, mapAddressing, weights),
-    capacity_(Field<Type>::size())
+    Field<T>(mapF, mapAddressing, weights),
+    capacity_(Field<T>::size())
 {}
 
 
 //- Construct by mapping from the given field
-template<class Type>
-Foam::DynamicField<Type>::DynamicField
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
 (
-    const UList<Type>& mapF,
+    const UList<T>& mapF,
     const FieldMapper& map
 )
 :
-    DynamicField<Type>(mapF, map),
-    capacity_(Field<Type>::size())
+    Field<T>(mapF, map),
+    capacity_(Field<T>::size())
 {}
 
 
-template<class Type>
-Foam::DynamicField<Type>::DynamicField(const DynamicField<Type>& f)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
+(
+    const DynamicField<T, SizeInc, SizeMult, SizeDiv>& lst
+)
 :
-    Field<Type>(f),
-    capacity_(Field<Type>::size())
+    Field<T>(lst),
+    capacity_(lst.capacity())
 {}
 
 
-template<class Type>
-Foam::DynamicField<Type>::DynamicField(DynamicField<Type>& f, bool reUse)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::DynamicField
+(
+    const Xfer<DynamicField<T, SizeInc, SizeMult, SizeDiv> >& lst
+)
 :
-    Field<Type>(f, reUse),
-    capacity_(Field<Type>::size())
-{}
-
-
-template<class Type>
-Foam::DynamicField<Type>::DynamicField(const Xfer<DynamicField<Type> >& f)
-:
-    Field<Type>(f),
-    capacity_(Field<Type>::size())
+    Field<T>(lst),
+    capacity_(Field<T>::size())
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::label Foam::DynamicField<Type>::capacity() const
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::label Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::capacity()
+const
 {
     return capacity_;
 }
 
 
-template<class Type>
-void Foam::DynamicField<Type>::append(const Type& t)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::setCapacity
+(
+    const label nElem
+)
 {
-    label elemI = Field<Type>::size();
+    label nextFree = Field<T>::size();
+    capacity_ = nElem;
+
+    if (nextFree > capacity_)
+    {
+        // truncate addressed sizes too
+        nextFree = capacity_;
+    }
+    // we could also enforce SizeInc granularity when (!SizeMult || !SizeDiv)
+
+    Field<T>::setSize(capacity_);
+    Field<T>::size(nextFree);
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::reserve
+(
+    const label nElem
+)
+{
+    // allocate more capacity?
+    if (nElem > capacity_)
+    {
+// TODO: convince the compiler that division by zero does not occur
+//        if (SizeInc && (!SizeMult || !SizeDiv))
+//        {
+//            // resize with SizeInc as the granularity
+//            capacity_ = nElem;
+//            unsigned pad = SizeInc - (capacity_ % SizeInc);
+//            if (pad != SizeInc)
+//            {
+//                capacity_ += pad;
+//            }
+//        }
+//        else
+        {
+            capacity_ = max
+            (
+                nElem,
+                label(SizeInc + capacity_ * SizeMult / SizeDiv)
+            );
+        }
+
+        // adjust allocated size, leave addressed size untouched
+        label nextFree = Field<T>::size();
+        Field<T>::setSize(capacity_);
+        Field<T>::size(nextFree);
+    }
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::setSize
+(
+    const label nElem
+)
+{
+    // allocate more capacity?
+    if (nElem > capacity_)
+    {
+// TODO: convince the compiler that division by zero does not occur
+//        if (SizeInc && (!SizeMult || !SizeDiv))
+//        {
+//            // resize with SizeInc as the granularity
+//            capacity_ = nElem;
+//            unsigned pad = SizeInc - (capacity_ % SizeInc);
+//            if (pad != SizeInc)
+//            {
+//                capacity_ += pad;
+//            }
+//        }
+//        else
+        {
+            capacity_ = max
+            (
+                nElem,
+                label(SizeInc + capacity_ * SizeMult / SizeDiv)
+            );
+        }
+
+        Field<T>::setSize(capacity_);
+    }
+
+    // adjust addressed size
+    Field<T>::size(nElem);
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::setSize
+(
+    const label nElem,
+    const T& t
+)
+{
+    label nextFree = Field<T>::size();
+    setSize(nElem);
+
+    // set new elements to constant value
+    while (nextFree < nElem)
+    {
+        this->operator[](nextFree++) = t;
+    }
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::resize
+(
+    const label nElem
+)
+{
+    this->setSize(nElem);
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::resize
+(
+    const label nElem,
+    const T& t
+)
+{
+    this->setSize(nElem, t);
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::clear()
+{
+    Field<T>::size(0);
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::clearStorage()
+{
+    Field<T>::clear();
+    capacity_ = 0;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::shrink()
+{
+    label nextFree = Field<T>::size();
+    if (capacity_ > nextFree)
+    {
+        // use the full list when resizing
+        Field<T>::size(capacity_);
+
+        // the new size
+        capacity_ = nextFree;
+        Field<T>::setSize(capacity_);
+        Field<T>::size(nextFree);
+    }
+    return *this;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::Xfer<Foam::List<T> >
+Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::xfer()
+{
+    return xferMoveTo< List<T> >(*this);
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::append
+(
+    const T& t
+)
+{
+    const label elemI = List<T>::size();
     setSize(elemI + 1);
 
     this->operator[](elemI) = t;
+    return *this;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>&
+Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::append
+(
+    const UList<T>& lst
+)
+{
+    if (this == &lst)
+    {
+        FatalErrorIn
+        (
+            "DynamicField<T, SizeInc, SizeMult, SizeDiv>::append"
+            "(const UList<T>&)"
+        )   << "attempted appending to self" << abort(FatalError);
+    }
+
+    label nextFree = List<T>::size();
+    setSize(nextFree + lst.size());
+
+    forAll(lst, elemI)
+    {
+        this->operator[](nextFree++) = lst[elemI];
+    }
+    return *this;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline T Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::remove()
+{
+    const label elemI = List<T>::size() - 1;
+
+    if (elemI < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::remove()"
+        )   << "List is empty" << abort(FatalError);
+    }
+
+    const T& val = List<T>::operator[](elemI);
+
+    List<T>::size(elemI);
+
+    return val;
 }
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class Type>
-void Foam::DynamicField<Type>::operator=(const DynamicField<Type>& rhs)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline T& Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::operator()
+(
+    const label elemI
+)
 {
-    if (this == &rhs)
+    if (elemI >= Field<T>::size())
     {
-        FatalErrorIn("DynamicField<Type>::operator=(const DynamicField<Type>&)")
-            << "attempted assignment to self"
-            << abort(FatalError);
+        setSize(elemI + 1);
     }
 
-    Field<Type>::operator=(rhs);
-    capacity_ = Field<Type>::size();
+    return this->operator[](elemI);
 }
 
 
-template<class Type>
-void Foam::DynamicField<Type>::operator=(const UList<Type>& rhs)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::operator=
+(
+    const T& t
+)
 {
-    Field<Type>::operator=(rhs);
-    capacity_ = Field<Type>::size();
+    UList<T>::operator=(t);
 }
 
 
-template<class Type>
-void Foam::DynamicField<Type>::operator=(const tmp<DynamicField>& rhs)
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::operator=
+(
+    const DynamicField<T, SizeInc, SizeMult, SizeDiv>& lst
+)
 {
-    if (this == &(rhs()))
+    if (this == &lst)
     {
-        FatalErrorIn("DynamicField<Type>::operator=(const tmp<DynamicField>&)")
-            << "attempted assignment to self"
-            << abort(FatalError);
+        FatalErrorIn
+        (
+            "DynamicField<T, SizeInc, SizeMult, SizeDiv>::operator="
+            "(const DynamicField<T, SizeInc, SizeMult, SizeDiv>&)"
+        )   << "attempted assignment to self" << abort(FatalError);
     }
 
-    // This is dodgy stuff, don't try it at home.
-    DynamicField* fieldPtr = rhs.ptr();
-    List<Type>::transfer(*fieldPtr);
-    delete fieldPtr;
-    capacity_ = Field<Type>::size();
+    if (capacity_ >= lst.size())
+    {
+        // can copy w/o reallocating, match initial size to avoid reallocation
+        Field<T>::size(lst.size());
+        Field<T>::operator=(lst);
+    }
+    else
+    {
+        // make everything available for the copy operation
+        Field<T>::size(capacity_);
+
+        Field<T>::operator=(lst);
+        capacity_ = Field<T>::size();
+    }
 }
 
 
-// * * * * * * * * * * * * * * * IOstream Operator * * * * * * * * * * * * * //
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicField<T, SizeInc, SizeMult, SizeDiv>::operator=
+(
+    const UList<T>& lst
+)
+{
+    if (capacity_ >= lst.size())
+    {
+        // can copy w/o reallocating, match initial size to avoid reallocation
+        Field<T>::size(lst.size());
+        Field<T>::operator=(lst);
+    }
+    else
+    {
+        // make everything available for the copy operation
+        Field<T>::size(capacity_);
+
+        Field<T>::operator=(lst);
+        capacity_ = Field<T>::size();
+    }
+}
 
 
 // ************************************************************************* //

src/finiteVolume/Make/files

@@ -287,6 +287,7 @@ $(ddtSchemes)/steadyStateDdtScheme/steadyStateDdtSchemes.C
 $(ddtSchemes)/EulerDdtScheme/EulerDdtSchemes.C
 $(ddtSchemes)/CoEulerDdtScheme/CoEulerDdtSchemes.C
 $(ddtSchemes)/SLTSDdtScheme/SLTSDdtSchemes.C
+$(ddtSchemes)/localEulerDdtScheme/localEulerDdtSchemes.C
 $(ddtSchemes)/backwardDdtScheme/backwardDdtSchemes.C
 $(ddtSchemes)/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
 $(ddtSchemes)/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
@@ -337,6 +338,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

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C

@@ -0,0 +1,584 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "localEulerDdtScheme.H"
+#include "surfaceInterpolate.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+const volScalarField& localEulerDdtScheme<Type>::localRDeltaT() const
+{
+    return mesh().objectRegistry::lookupObject<volScalarField>(rDeltaTName_);
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+localEulerDdtScheme<Type>::fvcDdt
+(
+    const dimensioned<Type>& dt
+)
+{
+    const volScalarField& rDeltaT = localRDeltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddt(" + dt.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if (mesh().moving())
+    {
+        tmp<GeometricField<Type, fvPatchField, volMesh> > tdtdt
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                mesh(),
+                dimensioned<Type>
+                (
+                    "0",
+                    dt.dimensions()/dimTime,
+                    pTraits<Type>::zero
+                )
+            )
+        );
+
+        tdtdt().internalField() =
+            rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
+
+        return tdtdt;
+    }
+    else
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                mesh(),
+                dimensioned<Type>
+                (
+                    "0",
+                    dt.dimensions()/dimTime,
+                    pTraits<Type>::zero
+                ),
+                calculatedFvPatchField<Type>::typeName
+            )
+        );
+    }
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+localEulerDdtScheme<Type>::fvcDdt
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const volScalarField& rDeltaT = localRDeltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddt(" + vf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if (mesh().moving())
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                mesh(),
+                rDeltaT.dimensions()*vf.dimensions(),
+                rDeltaT.internalField()*
+                (
+                    vf.internalField()
+                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                ),
+                rDeltaT.boundaryField()*
+                (
+                    vf.boundaryField() - vf.oldTime().boundaryField()
+                )
+            )
+        );
+    }
+    else
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                rDeltaT*(vf - vf.oldTime())
+            )
+        );
+    }
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+localEulerDdtScheme<Type>::fvcDdt
+(
+    const dimensionedScalar& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const volScalarField& rDeltaT = localRDeltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddt(" + rho.name() + ',' + vf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if (mesh().moving())
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                mesh(),
+                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+                rDeltaT.internalField()*rho.value()*
+                (
+                    vf.internalField()
+                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                ),
+                rDeltaT.boundaryField()*rho.value()*
+                (
+                    vf.boundaryField() - vf.oldTime().boundaryField()
+                )
+            )
+        );
+    }
+    else
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                rDeltaT*rho*(vf - vf.oldTime())
+            )
+        );
+    }
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+localEulerDdtScheme<Type>::fvcDdt
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const volScalarField& rDeltaT = localRDeltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddt(" + rho.name() + ',' + vf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if (mesh().moving())
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                mesh(),
+                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
+                rDeltaT.internalField()*
+                (
+                    rho.internalField()*vf.internalField()
+                  - rho.oldTime().internalField()
+                   *vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                ),
+                rDeltaT.boundaryField()*
+                (
+                    rho.boundaryField()*vf.boundaryField()
+                  - rho.oldTime().boundaryField()
+                   *vf.oldTime().boundaryField()
+                )
+            )
+        );
+    }
+    else
+    {
+        return tmp<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            new GeometricField<Type, fvPatchField, volMesh>
+            (
+                ddtIOobject,
+                rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
+            )
+        );
+    }
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type> >
+localEulerDdtScheme<Type>::fvmDdt
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            vf.dimensions()*dimVol/dimTime
+        )
+    );
+
+    fvMatrix<Type>& fvm = tfvm();
+
+    const scalarField& rDeltaT = localRDeltaT().internalField();
+
+    fvm.diag() = rDeltaT*mesh().V();
+
+    if (mesh().moving())
+    {
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
+    }
+    else
+    {
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
+    }
+
+    return tfvm;
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type> >
+localEulerDdtScheme<Type>::fvmDdt
+(
+    const dimensionedScalar& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            rho.dimensions()*vf.dimensions()*dimVol/dimTime
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm();
+
+    const scalarField& rDeltaT = localRDeltaT().internalField();
+
+    fvm.diag() = rDeltaT*rho.value()*mesh().V();
+
+    if (mesh().moving())
+    {
+        fvm.source() = rDeltaT
+            *rho.value()*vf.oldTime().internalField()*mesh().V0();
+    }
+    else
+    {
+        fvm.source() = rDeltaT
+            *rho.value()*vf.oldTime().internalField()*mesh().V();
+    }
+
+    return tfvm;
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type> >
+localEulerDdtScheme<Type>::fvmDdt
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type> > tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            rho.dimensions()*vf.dimensions()*dimVol/dimTime
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm();
+
+    const scalarField& rDeltaT = localRDeltaT().internalField();
+
+    fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
+
+    if (mesh().moving())
+    {
+        fvm.source() = rDeltaT
+            *rho.oldTime().internalField()
+            *vf.oldTime().internalField()*mesh().V0();
+    }
+    else
+    {
+        fvm.source() = rDeltaT
+            *rho.oldTime().internalField()
+            *vf.oldTime().internalField()*mesh().V();
+    }
+
+    return tfvm;
+}
+
+
+template<class Type>
+tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
+localEulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rA,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if (mesh().moving())
+    {
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                mesh(),
+                dimensioned<typename flux<Type>::type>
+                (
+                    "0",
+                    rA.dimensions()*phi.dimensions()/dimTime,
+                    pTraits<typename flux<Type>::type>::zero
+                )
+            )
+        );
+    }
+    else
+    {
+        const volScalarField& rDeltaT = localRDeltaT();
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
+                (
+                    fvc::interpolate(rDeltaT*rA)*phi.oldTime()
+                  - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
+                )
+            )
+        );
+    }
+}
+
+
+template<class Type>
+tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
+localEulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rA,
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtPhiCorr("
+      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if (mesh().moving())
+    {
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                mesh(),
+                dimensioned<typename flux<Type>::type>
+                (
+                    "0",
+                    rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
+                    pTraits<typename flux<Type>::type>::zero
+                )
+            )
+        );
+    }
+    else
+    {
+        const volScalarField& rDeltaT = localRDeltaT();
+
+        if
+        (
+            U.dimensions() == dimVelocity
+         && phi.dimensions() == dimVelocity*dimArea
+        )
+        {
+            return tmp<fluxFieldType>
+            (
+                new fluxFieldType
+                (
+                    ddtIOobject,
+                    fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+                   *(
+                        fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
+                      - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
+                      & mesh().Sf())
+                    )
+                )
+            );
+        }
+        else if
+        (
+            U.dimensions() == dimVelocity
+         && phi.dimensions() == dimDensity*dimVelocity*dimArea
+        )
+        {
+            return tmp<fluxFieldType>
+            (
+                new fluxFieldType
+                (
+                    ddtIOobject,
+                    fvcDdtPhiCoeff
+                    (
+                        U.oldTime(),
+                        phi.oldTime()/fvc::interpolate(rho.oldTime())
+                    )
+                   *(
+                        fvc::interpolate(rDeltaT*rA*rho.oldTime())
+                       *phi.oldTime()/fvc::interpolate(rho.oldTime())
+                      - (
+                            fvc::interpolate
+                            (
+                                rDeltaT*rA*rho.oldTime()*U.oldTime()
+                            ) & mesh().Sf()
+                        )
+                    )
+                )
+            );
+        }
+        else if
+        (
+            U.dimensions() == dimDensity*dimVelocity
+         && phi.dimensions() == dimDensity*dimVelocity*dimArea
+        )
+        {
+            return tmp<fluxFieldType>
+            (
+                new fluxFieldType
+                (
+                    ddtIOobject,
+                    fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())
+                   *(
+                        fvc::interpolate(rDeltaT*rA)*phi.oldTime()
+                      - (
+                            fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
+                        )
+                    )
+                )
+            );
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "localEulerDdtScheme<Type>::fvcDdtPhiCorr"
+            )   << "dimensions of phi are not correct"
+                << abort(FatalError);
+
+            return fluxFieldType::null();
+        }
+    }
+}
+
+
+template<class Type>
+tmp<surfaceScalarField> localEulerDdtScheme<Type>::meshPhi
+(
+    const GeometricField<Type, fvPatchField, volMesh>&
+)
+{
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "meshPhi",
+                mesh().time().timeName(),
+                mesh()
+            ),
+            mesh(),
+            dimensionedScalar("0", dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::fv::localEulerDdtScheme
+
+Description
+    Local time-step first-order Euler implicit/explicit ddt.
+    The reciprocal of the local time-step field is looked-up from the
+    database with the name provided.
+
+    This scheme should only be used for steady-state computations
+    using transient codes where local time-stepping is preferably to
+    under-relaxation for transport consistency reasons.
+
+    See also CoEulerDdtScheme.
+
+SourceFiles
+    localEulerDdtScheme.C
+    localEulerDdtSchemes.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef localEulerDdtScheme_H
+#define localEulerDdtScheme_H
+
+#include "ddtScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class localEulerDdtScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class localEulerDdtScheme
+:
+    public fv::ddtScheme<Type>
+{
+    // Private Data
+
+        //- Name of the reciprocal local time-step field
+        word rDeltaTName_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        localEulerDdtScheme(const localEulerDdtScheme&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const localEulerDdtScheme&);
+
+        //- Return the reciprocal of the local time-step
+        const volScalarField& localRDeltaT() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("localEuler");
+
+
+    // Constructors
+
+        //- Construct from mesh and Istream
+        localEulerDdtScheme(const fvMesh& mesh, Istream& is)
+        :
+            ddtScheme<Type>(mesh, is),
+            rDeltaTName_(is)
+        {}
+
+
+    // Member Functions
+
+        //- Return mesh reference
+        const fvMesh& mesh() const
+        {
+            return fv::ddtScheme<Type>::mesh();
+        }
+
+        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDdt
+        (
+            const dimensioned<Type>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDdt
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDdt
+        (
+            const dimensionedScalar&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDdt
+        (
+            const volScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<fvMatrix<Type> > fvmDdt
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<fvMatrix<Type> > fvmDdt
+        (
+            const dimensionedScalar&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<fvMatrix<Type> > fvmDdt
+        (
+            const volScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const volScalarField& rA,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const volScalarField& rA,
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<surfaceScalarField> meshPhi
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+};
+
+
+template<>
+tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
+(
+    const volScalarField& rA,
+    const volScalarField& U,
+    const surfaceScalarField& phi
+);
+
+
+template<>
+tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
+(
+    const volScalarField& rA,
+    const volScalarField& rho,
+    const volScalarField& U,
+    const surfaceScalarField& phi
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "localEulerDdtScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtSchemes.C

@@ -21,28 +21,19 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Calculates and outputs the mean and maximum Courant Numbers for the fluid
-    regions
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef compressibleCourantNo_H
-#define compressibleCourantNo_H
-
+#include "localEulerDdtScheme.H"
 #include "fvMesh.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 namespace Foam
 {
-    scalar compressibleCourantNo
-    (
-        const fvMesh& mesh,
-        const Time& runTime,
-        const volScalarField& rho,
-        const surfaceScalarField& phi
-    );
+namespace fv
+{
+    makeFvDdtScheme(localEulerDdtScheme)
+}
 }
 
-#endif
-
 // ************************************************************************* //

src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C

@@ -0,0 +1,324 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 scalar alphaDiff,
+    const scalar alphaMax,
+    const scalar alphaMin
+)
+{
+    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] > alphaMin && alpha[own] < alphaMax)
+         || (alpha[nbr] > alphaMin && alpha[nbr] < alphaMax)
+        )
+        {
+            if (mag(alpha[own] - alpha[nbr]) > alphaDiff)
+            {
+                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] > alphaMin && alpha[own] < alphaMax)
+                 || (
+                        alphapn[patchFacei] > alphaMin
+                     && alphapn[patchFacei] < alphaMax
+                    )
+                )
+                {
+                    if (mag(alpha[own] - alphapn[patchFacei]) > alphaDiff)
+                    {
+                        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 scalar alphaDiff
+)
+{
+    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]) > alphaDiff)
+        {
+            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]) > alphaDiff)
+                {
+                    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();
+}
+
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+InNamespace
+    Foam::fvc
+
+Description
+    Provides functions smooth spread and sweep which use the FaceCellWave
+    algorithm to smooth and redistribute the first field argument.
+
+    smooth: smooths the field by ensuring the values in neighbouring cells are
+            at least coeff* the cell value.
+
+    spread: redistributes the field by spreading the maximum value within the
+            region defined by the value (being between alphaMax and alphaMin)
+            and gradient of alpha (where the difference between the values in
+            neighbouring cells is larger than alphaDiff).
+
+    sweep: redistributes the field by sweeping the maximum value where the
+           gradient of alpha is large (where the difference between the values
+           in neighbouring cells is larger than alphaDiff) away from that
+           starting point of the sweep.
+
+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,
+        const scalar alphaDiff = 0.2,
+        const scalar alphaMax = 0.99,
+        const scalar alphaMin = 0.01
+    );
+
+    void sweep
+    (
+        volScalarField& field,
+        const volScalarField& alpha,
+        const label nLayers,
+        const scalar alphaDiff = 0.2
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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 fvc::smooth and fvc::spread functions.
+
+SourceFiles
+    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
+
+// ************************************************************************* //

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

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 fvc::sweep function.
+
+SourceFiles
+    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
+
+// ************************************************************************* //

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

src/meshTools/meshSearch/meshSearch.C

@@ -617,30 +617,25 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
         {
             label faceI = cFaces[i];
 
-            const face& f = mesh_.faces()[faceI];
+            pointHit inter = mesh_.faces()[faceI].ray
+            (
+                ctr,
+                dir,
+                mesh_.points(),
+                intersection::HALF_RAY,
+                intersection::VECTOR
+            );
 
-            forAll(f, fp)
+            if (inter.hit())
             {
-                pointHit inter = f.ray
-                (
-                    ctr,
-                    dir,
-                    mesh_.points(),
-                    intersection::HALF_RAY,
-                    intersection::VECTOR
-                );
+                scalar dist = inter.distance();
 
-                if (inter.hit())
+                if (dist < magDir)
                 {
-                    scalar dist = inter.distance();
+                    // Valid hit. Hit face so point is not in cell.
+                    intersection::setPlanarTol(oldTol);
 
-                    if (dist < magDir)
-                    {
-                        // Valid hit. Hit face so point is not in cell.
-                        intersection::setPlanarTol(oldTol);
-
-                        return false;
-                    }
+                    return false;
                 }
             }
         }

src/postProcessing/functionObjects/field/streamLine/streamLine.C

@@ -53,16 +53,10 @@ void Foam::streamLine::track()
         initialParticles
     );
 
-    //Pout<< "Seeding particles." << endl;
-
     const sampledSet& seedPoints = sampledSetPtr_();
 
     forAll(seedPoints, i)
     {
-        //Pout<< "Seeded particle at " << seedPoints[i]
-        //    << " at cell:" << seedPoints.cells()[i]
-        //    << endl;
-
         particles.addParticle
         (
             new streamLineParticle
@@ -228,15 +222,14 @@ void Foam::streamLine::track()
         vvInterp,
         UIndex,         // index of U in vvInterp
         trackForward_,  // track in +u direction
+        nSubCycle_,     // step through cells in steps?
         allTracks_,
         allScalars_,
         allVectors_
     );
 
     // Track
-    //Pout<< "Tracking particles." << endl;
     particles.move(td);
-    //Pout<< "Finished tracking particles." << endl;
 }
 
 
@@ -296,6 +289,19 @@ void Foam::streamLine::read(const dictionary& dict)
         UName_ = dict.lookupOrDefault<word>("U", "U");
         dict.lookup("trackForward") >> trackForward_;
         dict.lookup("lifeTime") >> lifeTime_;
+        if (lifeTime_ < 1)
+        {
+            FatalErrorIn(":streamLine::read(const dictionary&)")
+                << "Illegal value " << lifeTime_ << " for lifeTime"
+                << exit(FatalError);
+        }
+
+        dict.lookup("nSubCycle") >> nSubCycle_;
+        if (nSubCycle_ < 1)
+        {
+            nSubCycle_ = 1;
+        }
+
         cloudName_ = dict.lookupOrDefault<word>("cloudName", "streamLine");
         dict.lookup("seedSampleSet") >> seedSet_;
 
@@ -323,7 +329,6 @@ void Foam::streamLine::execute()
 {
 //    const Time& runTime = const_cast<Time&>(obr_.time());
 //    Pout<< "**streamLine::execute : time:" << runTime.timeName() << endl;
-//    Pout<< "**streamLine::execute : time:" << runTime.timeIndex() << endl;
 //
 //    bool isOutputTime = false;
 //

src/postProcessing/functionObjects/field/streamLine/streamLine.H

@@ -92,6 +92,9 @@ class streamLine
         //- Maximum lifetime (= number of cells) of particle
         label lifeTime_;
 
+        //- Number of subcycling steps
+        label nSubCycle_;
+
         //- Optional specified name of particles
         word cloudName_;
 

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C

@@ -36,6 +36,25 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+// Estimate dt to cross cell in a few steps
+Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
+(
+    streamLineParticle::trackData& td,
+    const scalar dt,
+    const vector& U
+) const
+{
+    streamLineParticle testParticle(*this);
+    bool oldKeepParticle = td.keepParticle;
+    bool oldSwitchProcessor = td.switchProcessor;
+    scalar fraction = testParticle.trackToFace(position()+dt*U, td);
+    td.keepParticle = oldKeepParticle;
+    td.switchProcessor = oldSwitchProcessor;
+    // Adapt the dt to subdivide the trajectory into 4 substeps.
+    return dt*fraction/td.nSubCycle_;
+}
+
+
 Foam::vector Foam::streamLineParticle::interpolateFields
 (
     const streamLineParticle::trackData& td,
@@ -170,34 +189,60 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
     && lifeTime_ > 0
     )
     {
-        // TBD: implement subcycling so step through cells in more than
-        //      one step.
-
         // set the lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        // Store current position and sampled velocity.
-        --lifeTime_;
-        sampledPositions_.append(position());
-        vector U = interpolateFields(td, position(), cell());
-
-        if (!td.trackForward_)
+        // Cross cell in steps:
+        // - at subiter 0 calculate dt to cross cell in nSubCycle steps
+        // - at the last subiter do all of the remaining track
+        // - do a few more subiters than nSubCycle since velocity might
+        //   be decreasing
+        for (label subIter = 0; subIter < 2*td.nSubCycle_; subIter++)
         {
-            U = -U;
-        }
+            --lifeTime_;
 
-        dt *= trackToFace(position()+dt*U, td);
+            // Store current position and sampled velocity.
+            sampledPositions_.append(position());
+            vector U = interpolateFields(td, position(), cell());
 
-        tEnd -= dt;
-        stepFraction() = 1.0 - tEnd/deltaT;
+            if (!td.trackForward_)
+            {
+                U = -U;
+            }
 
-        if (tEnd <= ROOTVSMALL)
-        {
-            // Force removal
-            lifeTime_ = 0;
+
+            if (subIter == 0 && td.nSubCycle_ > 1)
+            {
+                // Adapt dt to cross cell in a few steps
+                dt = calcSubCycleDeltaT(td, dt, U);
+            }
+            else if (subIter == td.nSubCycle_-1)
+            {
+                // Do full step on last subcycle
+                dt = min(dtMax, tEnd);
+            }
+
+
+            scalar fraction = trackToFace(position()+dt*U, td);
+            dt *= fraction;
+
+            tEnd -= dt;
+            stepFraction() = 1.0 - tEnd/deltaT;
+
+            if (tEnd <= ROOTVSMALL)
+            {
+                // Force removal
+                lifeTime_ = 0;
+            }
+
+            if (!td.keepParticle || td.switchProcessor || lifeTime_ == 0)
+            {
+                break;
+            }
         }
     }
 
+
     if (!td.keepParticle || lifeTime_ == 0)
     {
         if (lifeTime_ == 0)

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H

@@ -72,6 +72,7 @@ public:
         const PtrList<interpolationCellPoint<vector> >& vvInterp_;
         const label UIndex_;
         const bool trackForward_;
+        const label nSubCycle_;
 
         DynamicList<vectorList>& allPositions_;
         List<DynamicList<scalarList> >& allScalars_;
@@ -87,6 +88,7 @@ public:
                 const PtrList<interpolationCellPoint<vector> >& vvInterp,
                 const label UIndex,
                 const bool trackForward,
+                const label nSubCycle,
                 DynamicList<List<point> >& allPositions,
                 List<DynamicList<scalarList> >& allScalars,
                 List<DynamicList<vectorList> >& allVectors
@@ -97,6 +99,7 @@ public:
                 vvInterp_(vvInterp),
                 UIndex_(UIndex),
                 trackForward_(trackForward),
+                nSubCycle_(nSubCycle),
                 allPositions_(allPositions),
                 allScalars_(allScalars),
                 allVectors_(allVectors)
@@ -123,6 +126,15 @@ private:
 
     // Private Member Functions
 
+        //- Estimate dt to cross from current face to next one in nSubCycle
+        //  steps.
+        scalar calcSubCycleDeltaT
+        (
+            streamLineParticle::trackData& td,
+            const scalar dt,
+            const vector& U
+        ) const;
+
         //- Interpolate all quantities; return interpolated velocity.
         vector interpolateFields
         (

tutorials/incompressible/simpleFoam/motorBike/system/controlDict

@@ -82,8 +82,11 @@ functions
         // Names of fields to sample. Should contain above velocity field!
         fields (p U k);
 
-        // Cells particles can travel before being removed
-        lifeTime        1000;
+        // Steps particles can travel before being removed
+        lifeTime        10000;
+
+        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+        nSubCycle 5;
 
         // Cloud name to use
         cloudName       particleTracks;

tutorials/incompressible/simpleFoam/pitzDaily/system/controlDict

@@ -84,8 +84,11 @@ functions
         // Names of fields to sample. Should contain above velocity field!
         fields (p k U);
 
-        // Cells particles can travel before being removed
-        lifeTime        1000;
+        // Steps particles can travel before being removed
+        lifeTime        10000;
+
+        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+        nSubCycle 5;
 
         // Cloud name to use
         cloudName       particleTracks;