Merge remote branch 'OpenCFD/master' into olesenm

ReleaseNotes-dev

@@ -29,13 +29,29 @@
   =metis-5.0pre2=, =ParMetis-3.1= and =ParMGridGen-1.0=.
 
 * Library developments
+
   There have been a number of developments to the libraries to support the
   extension of functionality in solver and utility applications.
 *** Core library
     + Large number of code refinements and consistency improvements to support
       other developments.
 *** Turbulence modelling
-*** Thermo-physical Models
+*** *Updated* thermophysical libraries
+    + Now use the more descriptive dictionary input format, e.g. for a single
+      component mixture:
+
+    thermoType      hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+
+    mixture
+    {
+        nMoles      1;
+        molWeight   28.9;
+        Cp          1007;
+        Hf          0;
+        As          1.4792e-06;
+        Ts          116;
+    }
+
 *** DSMC
 *** Dynamic Mesh
 *** Numerics
@@ -50,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
@@ -67,8 +88,17 @@
     + Additional wall functions for primary region momentum and temperature
       taking film into account
     + Parallel aware
+*** *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
@@ -107,12 +137,18 @@
         (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
+      makes rereading the boundary superfluous. see
+      e.g. pitzDailyDirectMapped tutorial.
 * Post-processing
-  + =foamToEnsight=: new =-nodeValues= option to generate and output nodal
+  + =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
@@ -124,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/advanced/selectCells/selectCells.C

@@ -55,6 +55,7 @@ Description
 #include "edgeStats.H"
 #include "treeDataTriSurface.H"
 #include "indexedOctree.H"
+#include "globalMeshData.H"
 
 using namespace Foam;
 
@@ -298,7 +299,7 @@ label selectOutsideCells
         mesh,
         outsideFaces.shrink(),
         outsideFacesInfo.shrink(),
-        mesh.nCells()  // max iterations
+        mesh.globalData().nTotalCells()+1  // max iterations
     );
 
     // Now regionCalc should hold info on cells that are reachable from

applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L

@@ -43,7 +43,7 @@ Description
 #include "polyMeshZipUpCells.H"
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "Swap.H"
 #include "IFstream.H"
 #include "readHexLabel.H"
@@ -900,7 +900,7 @@ int main(int argc, char *argv[])
     fluentToFoamType.insert("interface", polyPatch::typeName);
     fluentToFoamType.insert("internal", polyPatch::typeName);
     fluentToFoamType.insert("solid", polyPatch::typeName);
-    fluentToFoamType.insert("fan", cyclicPolyPatch::typeName);
+    fluentToFoamType.insert("fan", oldCyclicPolyPatch::typeName);
     fluentToFoamType.insert("radiator", polyPatch::typeName);
     fluentToFoamType.insert("porous-jump", polyPatch::typeName);
 

applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C

@@ -41,7 +41,7 @@ Description
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "wedgePolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "unitConversion.H"
 
 using namespace Foam;

applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H

@@ -196,7 +196,7 @@ const word* kivaPatchTypes[nBCs] =
     &polyPatch::typeName,
     &polyPatch::typeName,
     &symmetryPolyPatch::typeName,
-    &cyclicPolyPatch::typeName
+    &oldCyclicPolyPatch::typeName
 };
 
 enum patchTypeNames

applications/utilities/mesh/conversion/sammToFoam/readBoundary.C

@@ -29,7 +29,7 @@ Description
 #include "sammMesh.H"
 #include "Time.H"
 #include "wallPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "preservePatchTypes.H"
 #include "IFstream.H"
@@ -208,7 +208,7 @@ void sammMesh::readBoundary()
             {
                 // incorrect. should be cyclicPatch but this
                 // requires info on connected faces.
-                patchTypes_[patchLabel] = cyclicPolyPatch::typeName;
+                patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName;
             }
             else
             {

applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C

@@ -29,7 +29,7 @@ Description
 #include "starMesh.H"
 #include "Time.H"
 #include "wallPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "preservePatchTypes.H"
 #include "IFstream.H"
@@ -206,7 +206,7 @@ void starMesh::readBoundary()
             {
                 // incorrect. should be cyclicPatch but this
                 // requires info on connected faces.
-                patchTypes_[patchLabel] = cyclicPolyPatch::typeName;
+                patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName;
             }
             else
             {

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/mesh/manipulation/subsetMesh/Make/options

@@ -4,4 +4,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lgenericPatchFields

applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H

@@ -57,7 +57,6 @@ public:
         labelList prisms;
         labelList wedges;
         labelList hexes;
-        labelList hexesWedges;
         labelList polys;
 
 
@@ -77,7 +76,6 @@ public:
             prisms(nCells),
             wedges(nCells),
             hexes(nCells),
-            hexesWedges(nCells),
             polys(nCells)
         {}
 };

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C

@@ -38,6 +38,32 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+Field<Type> map
+(
+    const Field<Type>& vf,
+    const labelList& map1,
+    const labelList& map2
+)
+{
+    Field<Type> mf(map1.size() + map2.size());
+
+    forAll(map1, i)
+    {
+        mf[i] = vf[map1[i]];
+    }
+
+    label offset = map1.size();
+
+    forAll(map2, i)
+    {
+        mf[i + offset] = vf[map2[i]];
+    }
+
+    return mf;
+}
+
+
 template<class Type>
 void writeField
 (
@@ -274,7 +300,8 @@ void ensightField
     const labelList& tets = meshCellSets.tets;
     const labelList& pyrs = meshCellSets.pyrs;
     const labelList& prisms = meshCellSets.prisms;
-    const labelList& hexesWedges = meshCellSets.hexesWedges;
+    const labelList& wedges = meshCellSets.wedges;
+    const labelList& hexes = meshCellSets.hexes;
     const labelList& polys = meshCellSets.polys;
 
     ensightStream* ensightFilePtr = NULL;
@@ -328,7 +355,7 @@ void ensightField
         writeField
         (
             "hexa8",
-            Field<Type>(vf, hexesWedges),
+            map(vf, hexes, wedges),
             ensightFile
         );
 
@@ -548,99 +575,99 @@ void ensightPointField
     }
 
 
-     label ensightPatchI = eMesh.patchPartOffset();
+    label ensightPatchI = eMesh.patchPartOffset();
  
-     forAll(allPatchNames, patchi)
-     {
-         const word& patchName = allPatchNames[patchi];
+    forAll(allPatchNames, patchi)
+    {
+        const word& patchName = allPatchNames[patchi];
  
-         eMesh.barrier();
+        eMesh.barrier();
  
-         if (patchNames.empty() || patchNames.found(patchName))
-         {
-             const fvPatch& p = mesh.boundary()[patchi];
-             if
-             (
-                 returnReduce(p.size(), sumOp<label>())
-               > 0
-             )
-             {
-                 // Renumber the patch points/faces into unique points
-                 labelList pointToGlobal;
-                 labelList uniqueMeshPointLabels;
-                 autoPtr<globalIndex> globalPointsPtr =
-                 mesh.globalData().mergePoints
-                 (
-                     p.patch().meshPoints(),
-                     p.patch().meshPointMap(),
-                     pointToGlobal,
-                     uniqueMeshPointLabels
-                 );
+        if (patchNames.empty() || patchNames.found(patchName))
+        {
+            const fvPatch& p = mesh.boundary()[patchi];
+            if
+            (
+                returnReduce(p.size(), sumOp<label>())
+              > 0
+            )
+            {
+                // Renumber the patch points/faces into unique points
+                labelList pointToGlobal;
+                labelList uniqueMeshPointLabels;
+                autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    p.patch().meshPoints(),
+                    p.patch().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
  
-                 if (Pstream::master())
-                 {
-                     ensightFile.writePartHeader(ensightPatchI);
-                 }
+                if (Pstream::master())
+                {
+                    ensightFile.writePartHeader(ensightPatchI);
+                }
  
-                 writeField
-                 (
-                     "coordinates",
-                     Field<Type>(pf.internalField(), uniqueMeshPointLabels),
-                     ensightFile
-                 );
+                writeField
+                (
+                    "coordinates",
+                    Field<Type>(pf.internalField(), uniqueMeshPointLabels),
+                    ensightFile
+                );
  
-                 ensightPatchI++;
-             }
-         }
-     }
+                ensightPatchI++;
+            }
+        }
+    }
  
-     // write faceZones, if requested
-     if (faceZoneNames.size())
-     {
-         forAllConstIter(wordHashSet, faceZoneNames, iter)
-         {
-             const word& faceZoneName = iter.key();
- 
-             eMesh.barrier();
- 
-             label zoneID = mesh.faceZones().findZoneID(faceZoneName);
- 
-             const faceZone& fz = mesh.faceZones()[zoneID];
- 
-             if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
-             {
-                 // Renumber the faceZone points/faces into unique points
-                 labelList pointToGlobal;
-                 labelList uniqueMeshPointLabels;
-                 autoPtr<globalIndex> globalPointsPtr =
-                 mesh.globalData().mergePoints
-                 (
-                     fz().meshPoints(),
-                     fz().meshPointMap(),
-                     pointToGlobal,
-                     uniqueMeshPointLabels
-                 );
- 
-                 if (Pstream::master())
-                 {
-                     ensightFile.writePartHeader(ensightPatchI);
-                 }
- 
-                 writeField
-                 (
-                     "coordinates",
-                     Field<Type>
-                     (
-                         pf.internalField(),
-                         uniqueMeshPointLabels
-                     ),
-                     ensightFile
-                 );
- 
-                 ensightPatchI++;
-             }
-         }
-     }
+    // write faceZones, if requested
+    if (faceZoneNames.size())
+    {
+        forAllConstIter(wordHashSet, faceZoneNames, iter)
+        {
+            const word& faceZoneName = iter.key();
+
+            eMesh.barrier();
+
+            label zoneID = mesh.faceZones().findZoneID(faceZoneName);
+
+            const faceZone& fz = mesh.faceZones()[zoneID];
+
+            if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
+            {
+                // Renumber the faceZone points/faces into unique points
+                labelList pointToGlobal;
+                labelList uniqueMeshPointLabels;
+                autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    fz().meshPoints(),
+                    fz().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
+
+                if (Pstream::master())
+                {
+                    ensightFile.writePartHeader(ensightPatchI);
+                }
+
+                writeField
+                (
+                    "coordinates",
+                    Field<Type>
+                    (
+                        pf.internalField(),
+                        uniqueMeshPointLabels
+                    ),
+                    ensightFile
+                );
+
+                ensightPatchI++;
+            }
+        }
+    }
 
     if (Pstream::master())
     {

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C

@@ -70,12 +70,15 @@ void Foam::ensightMesh::correct()
         // Patches are output. Check that they're synced.
         mesh_.boundaryMesh().checkParallelSync(true);
 
-        allPatchNames_ = wordList::subList
-        (
-            mesh_.boundaryMesh().names(),
-            mesh_.boundary().size()
-          - mesh_.globalData().processorPatches().size()
-        );
+        allPatchNames_ = mesh_.boundaryMesh().names();
+        if (Pstream::parRun())
+        {
+            allPatchNames_.setSize
+            (
+                mesh_.boundary().size()
+              - mesh_.globalData().processorPatches().size()
+            );
+        }
 
         if (patches_)
         {
@@ -114,7 +117,6 @@ void Foam::ensightMesh::correct()
         labelList& prisms = meshCellSets_.prisms;
         labelList& wedges = meshCellSets_.wedges;
         labelList& hexes = meshCellSets_.hexes;
-        labelList& hexesWedges = meshCellSets_.hexesWedges;
         labelList& polys = meshCellSets_.polys;
 
         label nTets = 0;
@@ -122,7 +124,6 @@ void Foam::ensightMesh::correct()
         label nPrisms = 0;
         label nWedges = 0;
         label nHexes = 0;
-        label nHexesWedges = 0;
         label nPolys = 0;
 
         forAll(cellShapes, cellI)
@@ -145,12 +146,10 @@ void Foam::ensightMesh::correct()
             else if (cellModel == wedge)
             {
                 wedges[nWedges++] = cellI;
-                hexesWedges[nHexesWedges++] = cellI;
             }
             else if (cellModel == hex)
             {
                 hexes[nHexes++] = cellI;
-                hexesWedges[nHexesWedges++] = cellI;
             }
             else
             {
@@ -163,7 +162,6 @@ void Foam::ensightMesh::correct()
         prisms.setSize(nPrisms);
         wedges.setSize(nWedges);
         hexes.setSize(nHexes);
-        hexesWedges.setSize(nHexesWedges);
         polys.setSize(nPolys);
 
         meshCellSets_.nTets = nTets;
@@ -175,7 +173,7 @@ void Foam::ensightMesh::correct()
         meshCellSets_.nPrisms = nPrisms;
         reduce(meshCellSets_.nPrisms, sumOp<label>());
 
-        meshCellSets_.nHexesWedges = nHexesWedges;
+        meshCellSets_.nHexesWedges = nWedges+nHexes;
         reduce(meshCellSets_.nHexesWedges, sumOp<label>());
 
         meshCellSets_.nPolys = nPolys;

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

applications/utilities/preProcessing/setFields/setFields.C

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Description
-    Selects a cell set through a dictionary.
+    Set values on a selected set of cells/patchfaces through a dictionary.
 
 \*---------------------------------------------------------------------------*/
 
@@ -32,12 +32,13 @@ Description
 #include "fvMesh.H"
 #include "topoSetSource.H"
 #include "cellSet.H"
+#include "faceSet.H"
 #include "volFields.H"
 
 using namespace Foam;
 
 template<class Type>
-bool setFieldType
+bool setCellFieldType
 (
     const word& fieldTypeDesc,
     const fvMesh& mesh,
@@ -65,7 +66,8 @@ bool setFieldType
     // Check field exists
     if (fieldHeader.headerOk())
     {
-        Info<< "    Setting " << fieldHeader.headerClassName()
+        Info<< "    Setting internal values of "
+            << fieldHeader.headerClassName()
             << " " << fieldName << endl;
 
         fieldType field(fieldHeader, mesh);
@@ -96,7 +98,7 @@ bool setFieldType
     {
         WarningIn
         (
-            "void setFieldType"
+            "void setCellFieldType"
             "(const fvMesh& mesh, const labelList& selectedCells,"
             "Istream& fieldValueStream)"
         ) << "Field " << fieldName << " not found" << endl;
@@ -106,17 +108,17 @@ bool setFieldType
 }
 
 
-class setField
+class setCellField
 {
 
 public:
 
-    setField()
+    setCellField()
     {}
 
-    autoPtr<setField> clone() const
+    autoPtr<setCellField> clone() const
     {
-        return autoPtr<setField>(new setField());
+        return autoPtr<setCellField>(new setCellField());
     }
 
     class iNew
@@ -132,37 +134,200 @@ public:
             selectedCells_(selectedCells)
         {}
 
-        autoPtr<setField> operator()(Istream& fieldValues) const
+        autoPtr<setCellField> operator()(Istream& fieldValues) const
         {
             word fieldType(fieldValues);
 
             if
             (
                !(
-                    setFieldType<scalar>
+                    setCellFieldType<scalar>
                         (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setFieldType<vector>
+                 || setCellFieldType<vector>
                         (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setFieldType<sphericalTensor>
+                 || setCellFieldType<sphericalTensor>
                         (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setFieldType<symmTensor>
+                 || setCellFieldType<symmTensor>
                         (fieldType, mesh_, selectedCells_, fieldValues)
-                 || setFieldType<tensor>
+                 || setCellFieldType<tensor>
                         (fieldType, mesh_, selectedCells_, fieldValues)
                 )
             )
             {
-                WarningIn("setField::iNew::operator()(Istream& is)")
+                WarningIn("setCellField::iNew::operator()(Istream& is)")
                     << "field type " << fieldType << " not currently supported"
                     << endl;
             }
 
-            return autoPtr<setField>(new setField());
+            return autoPtr<setCellField>(new setCellField());
         }
     };
 };
 
 
+template<class Type>
+bool setFaceFieldType
+(
+    const word& fieldTypeDesc,
+    const fvMesh& mesh,
+    const labelList& selectedFaces,
+    Istream& fieldValueStream
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (fieldTypeDesc != fieldType::typeName + "Value")
+    {
+        return false;
+    }
+
+    word fieldName(fieldValueStream);
+
+    IOobject fieldHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    // Check field exists
+    if (fieldHeader.headerOk())
+    {
+        Info<< "    Setting patchField values of "
+            << fieldHeader.headerClassName()
+            << " " << fieldName << endl;
+
+        fieldType field(fieldHeader, mesh);
+
+        const Type& value = pTraits<Type>(fieldValueStream);
+
+        // Create flat list of selected faces and their value.
+        Field<Type> allBoundaryValues(mesh.nFaces()-mesh.nInternalFaces());
+        forAll(field.boundaryField(), patchi)
+        {
+            SubField<Type>
+            (
+                allBoundaryValues,
+                field.boundaryField()[patchi].size(),
+                field.boundaryField()[patchi].patch().start()
+              - mesh.nInternalFaces()
+            ).assign(field.boundaryField()[patchi]);
+        }
+
+        // Override
+        labelList nChanged(field.boundaryField().size(), 0);
+        forAll(selectedFaces, i)
+        {
+            label facei = selectedFaces[i];
+            if (mesh.isInternalFace(facei))
+            {
+                WarningIn("setFaceFieldType(..)")
+                    << "Ignoring internal face " << facei << endl;
+            }
+            else
+            {
+                label bFaceI = facei-mesh.nInternalFaces();
+                allBoundaryValues[bFaceI] = value;
+                nChanged[mesh.boundaryMesh().patchID()[bFaceI]]++;
+            }
+        }
+
+        Pstream::listCombineGather(nChanged, plusEqOp<label>());
+        Pstream::listCombineScatter(nChanged);
+
+        // Reassign.
+        forAll(field.boundaryField(), patchi)
+        {
+            if (nChanged[patchi] > 0)
+            {
+                Info<< "    On patch "
+                    << field.boundaryField()[patchi].patch().name()
+                    << " set " << nChanged[patchi] << " values" << endl;
+                field.boundaryField()[patchi] == SubField<Type>
+                (
+                    allBoundaryValues,
+                    field.boundaryField()[patchi].size(),
+                    field.boundaryField()[patchi].patch().start()
+                  - mesh.nInternalFaces()
+                );
+            }
+        }
+
+        field.write();
+    }
+    else
+    {
+        WarningIn
+        (
+            "void setFaceFieldType"
+            "(const fvMesh& mesh, const labelList& selectedFaces,"
+            "Istream& fieldValueStream)"
+        ) << "Field " << fieldName << " not found" << endl;
+    }
+
+    return true;
+}
+
+
+class setFaceField
+{
+
+public:
+
+    setFaceField()
+    {}
+
+    autoPtr<setFaceField> clone() const
+    {
+        return autoPtr<setFaceField>(new setFaceField());
+    }
+
+    class iNew
+    {
+        const fvMesh& mesh_;
+        const labelList& selectedFaces_;
+
+    public:
+
+        iNew(const fvMesh& mesh, const labelList& selectedFaces)
+        :
+            mesh_(mesh),
+            selectedFaces_(selectedFaces)
+        {}
+
+        autoPtr<setFaceField> operator()(Istream& fieldValues) const
+        {
+            word fieldType(fieldValues);
+
+            if
+            (
+               !(
+                    setFaceFieldType<scalar>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<vector>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<sphericalTensor>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<symmTensor>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<tensor>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                )
+            )
+            {
+                WarningIn("setFaceField::iNew::operator()(Istream& is)")
+                    << "field type " << fieldType << " not currently supported"
+                    << endl;
+            }
+
+            return autoPtr<setFaceField>(new setFaceField());
+        }
+    };
+};
+
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -194,10 +359,10 @@ int main(int argc, char *argv[])
     if (setFieldsDict.found("defaultFieldValues"))
     {
         Info<< "Setting field default values" << endl;
-        PtrList<setField> defaultFieldValues
+        PtrList<setCellField> defaultFieldValues
         (
             setFieldsDict.lookup("defaultFieldValues"),
-            setField::iNew(mesh, labelList(mesh.nCells()))
+            setCellField::iNew(mesh, labelList(mesh.nCells()))
         );
         Info<< endl;
     }
@@ -211,29 +376,54 @@ int main(int argc, char *argv[])
     {
         const entry& region = regions[regionI];
 
-        autoPtr<topoSetSource> cellSelector =
+        autoPtr<topoSetSource> source =
             topoSetSource::New(region.keyword(), mesh, region.dict());
 
-        cellSet selectedCellSet
-        (
-            mesh,
-            "cellSet",
-            mesh.nCells()/10+1  // Reasonable size estimate.
-        );
+        if (source().setType() == topoSetSource::CELLSETSOURCE)
+        {
+            cellSet selectedCellSet
+            (
+                mesh,
+                "cellSet",
+                mesh.nCells()/10+1  // Reasonable size estimate.
+            );
 
-        cellSelector->applyToSet
-        (
-            topoSetSource::NEW,
-            selectedCellSet
-        );
+            source->applyToSet
+            (
+                topoSetSource::NEW,
+                selectedCellSet
+            );
 
-        PtrList<setField> fieldValues
-        (
-            region.dict().lookup("fieldValues"),
-            setField::iNew(mesh, selectedCellSet.toc())
-        );
+            PtrList<setCellField> fieldValues
+            (
+                region.dict().lookup("fieldValues"),
+                setCellField::iNew(mesh, selectedCellSet.toc())
+            );
+        }
+        else if (source().setType() == topoSetSource::FACESETSOURCE)
+        {
+            faceSet selectedFaceSet
+            (
+                mesh,
+                "faceSet",
+                (mesh.nFaces()-mesh.nInternalFaces())/10+1
+            );
+
+            source->applyToSet
+            (
+                topoSetSource::NEW,
+                selectedFaceSet
+            );
+
+            PtrList<setFaceField> fieldValues
+            (
+                region.dict().lookup("fieldValues"),
+                setFaceField::iNew(mesh, selectedFaceSet.toc())
+            );
+        }
     }
 
+
     Info<< "\nEnd\n" << endl;
 
     return 0;

applications/utilities/preProcessing/setFields/setFieldsDict

@@ -23,6 +23,8 @@ defaultFieldValues
 
 regions
 (
+    // Set cell values
+    // (does zerogradient on boundaries)
     boxToCell
     {
         box (0 0 -1) (0.1461 0.292 1);
@@ -32,6 +34,17 @@ regions
             volScalarFieldValue gamma 1
         );
     }
+
+    // Set patch values (using ==)
+    boxToFace
+    {
+        box (0 0 -1) (0.1461 0.292 1);
+
+        fieldValues
+        (
+            volScalarFieldValue gamma 1
+        );
+    }
 );
 
 // ************************************************************************* //

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

@@ -56,7 +56,7 @@ cleanTimeDirectories()
         zeros=`printf %0${nZeros}d 0`
         nZeros=$(($nZeros + 1))
     done
-    rm -rf ./{[1-9]*,-[1-9]*,log,log.*,log-*,logSummary.*,.fxLock,*.xml,ParaView*,paraFoam*,*.OpenFOAM} > /dev/null 2>&1
+    rm -rf ./[1-9]* ./-[1-9]* ./log ./log.* ./log-* ./logSummary.* ./.fxLock ./*.xml ./ParaView* ./paraFoam* ./*.OpenFOAM > /dev/null 2>&1
 }
 
 cleanCase()
@@ -65,16 +65,21 @@ cleanCase()
     rm -rf processor* > /dev/null 2>&1
     rm -rf probes* > /dev/null 2>&1
     rm -rf forces* > /dev/null 2>&1
-
-    rm -rf system/machines \
-           constant/polyMesh/{allOwner*,cell*,face*,meshModifiers*} \
-           constant/polyMesh/{owner*,neighbour*,point*,edge*} \
-           constant/polyMesh/{cellLevel*,pointLevel*,refinementHistory*,surfaceIndex*} \
-           constant/{cellToRegion,cellLevel*,pointLevel*} \
-           constant/polyMesh/sets/ \
-           VTK \
-           sets/streamLines \
-           > /dev/null 2>&1
+    rm -rf sets > /dev/null 2>&1
+    rm -rf system/machines > /dev/null 2>&1
+    (cd constant/polyMesh && \
+        rm -rf \
+        allOwner* cell* face* meshModifiers* \
+        owner* neighbour* point* edge* \
+        cellLevel* pointLevel* refinementHistory* surfaceIndex* sets \
+        > /dev/null 2>&1 \
+    )
+    (cd constant && \
+        rm -rf \
+        cellToRegion cellLevel* pointLevel* \
+        > /dev/null 2>&1 \
+    )
+    rm -rf VTK > /dev/null 2>&1
 
     if [ -e constant/polyMesh/blockMeshDict.m4 ]
     then
@@ -96,12 +101,12 @@ removeCase()
 
 cleanSamples()
 {
-    rm -rf {sets,samples,sampleSurfaces} > /dev/null 2>&1
+    rm -rf ./sets ./samples ./sampleSurfaces > /dev/null 2>&1
 }
 
 cleanUcomponents()
 {
-    rm -rf 0/{Ux,Uy,Uz} > /dev/null 2>&1
+    rm -rf 0/Ux 0/Uy 0/Uz > /dev/null 2>&1
 }
 
 cleanApplication()

bin/tools/foamConfigurePaths

@@ -34,10 +34,10 @@ usage() {
     cat<<USAGE
 
 usage: ${0##*/}
-  --foamInstall dir       specify installation directory (e.g. /opt)
-  --projectName name      specify project name (e.g. openfoam170)
-  --archOption  arch      specify architecture option (only 32 or 64 applicable)
-  --paraviewInstall dir   specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
+  --foamInstall dir         specify installation directory (e.g. /opt)
+  --projectName name        specify project name (e.g. openfoam170)
+  --archOption  arch        specify architecture option (only 32 or 64 applicable)
+  --paraviewInstall dir     specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
 
 * hardcode paths to installation
 
@@ -46,6 +46,17 @@ USAGE
 }
 
 
+# Function to do replacement on file. Checks if any replacement has been done.
+# inlineSed <file> <sedCommand> <description>
+_inlineSed()
+{
+    backup=`tempfile`
+    cp $1 $backup
+    sed -i -e "$2" $1
+    cmp $1 $backup || usage "Failed : $3"
+}
+
+
 [ -f etc/bashrc -a -f etc/settings.sh ] || usage "Please run from top-level directory of installation"
 
 unset foamInstall projectName archOption paraviewInstall
@@ -59,30 +70,40 @@ do
         ;;
     --foamInstall)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        foamInstall="$2"
-        echo "Replacing foamInstall setting by $foamInstall"
-        sed -i -e '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" etc/bashrc
+    foamInstall="$2"
+        echo "** foamInstall:$foamInstall"
+
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" \
+            "Replacing foamInstall setting by $foamInstall"
         shift 2
         ;;
    --projectName)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        projectName="$2"
-        echo "Replacing WM_PROJECT_DIR setting by $projectName"
-        sed -i -e '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" etc/bashrc
+    projectName="$2"
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" \
+            "Replacing WM_PROJECT_DIR setting by $projectName"
         shift 2
         ;;
     --archOption)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        archOption="$2"
-        echo "Replacing WM_ARCH_OPTION setting by $archOption"
-        sed -i -e '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" etc/bashrc
+    archOption="$2"
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" \
+             "Replacing WM_ARCH_OPTION setting by $archOption"
         shift 2
         ;;
     --paraviewInstall)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        paraviewInstall="$2"
-        echo "Replacing ParaView_DIR setting by $paraviewInstall"
-        sed -i -e '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" etc/apps/paraview3/bashrc
+    paraviewInstall="$2"
+        _inlineSed \
+            etc/apps/paraview3/bashrc \
+            '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" \
+             "Replacing ParaView_DIR setting by $paraviewInstall"
         shift 2
         ;;
     *)
@@ -97,11 +118,14 @@ done
 #sed -i -e 's@WM_PROJECT=.*@WM_PROJECT='"$projectName@" etc/bashrc
 
 # Replace the WM_MPLIB always
-echo "Replacing WM_MPLIB setting by SYSTEMOPENMPI"
-sed -i -e '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' etc/bashrc
-
+_inlineSed \
+    etc/bashrc \
+    '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' \
+    "Replacing WM_MPLIB setting by SYSTEMOPENMPI"
 # Replace the compilerInstall always
-echo "Replacing compilerInstall setting by system"
-sed -i -e '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' etc/settings.sh
+_inlineSed \
+    etc/settings.sh \
+    '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' \
+    "Replacing compilerInstall setting by system"
 
 #------------------------------------------------------------------------------

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

@@ -39,7 +39,7 @@ Class
 #else
 #   include <sys/inotify.h>
 #   include <sys/ioctl.h>
-
+#   include <errno.h>
 #   define EVENT_SIZE  ( sizeof (struct inotify_event) )
 #   define EVENT_LEN   (EVENT_SIZE + 16)
 #   define EVENT_BUF_LEN     ( 1024 * EVENT_LEN )
@@ -144,7 +144,7 @@ namespace Foam
 
 #else
         //- File descriptor for the inotify instance
-        int fd;
+        int inotifyFd_;
 
         //- Current watchIDs and corresponding directory id
         DynamicList<label> dirWatches_;
@@ -153,23 +153,47 @@ namespace Foam
         //- initialise inotify
         inline fileMonitorWatcher(const label sz = 20)
         :
-            fd(inotify_init()),
+            inotifyFd_(inotify_init()),
             dirWatches_(sz),
             dirFiles_(sz)
-        {}
+        {
+            if (inotifyFd_ < 0)
+            {
+                static bool hasWarned = false;
+                if (!hasWarned)
+                {
+                    hasWarned = true;
+                    WarningIn("fileMonitorWatcher(const label)")
+                        << "Failed allocating an inotify descriptor : "
+                        << string(strerror(errno)) << endl
+                        << "    Please increase the number of allowable "
+                        << "inotify instances" << endl
+                        << "    (/proc/sys/fs/inotify/max_user_instances"
+                        << " on Linux)" << endl
+                        << "    , 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;
+                }
+            }
+        }
 
         //- remove all watches
         inline ~fileMonitorWatcher()
         {
-            forAll(dirWatches_, i)
+            if (inotifyFd_ >= 0)
             {
-                if (dirWatches_[i] >= 0)
+                forAll(dirWatches_, i)
                 {
-                    if (inotify_rm_watch(fd, int(dirWatches_[i])))
+                    if (dirWatches_[i] >= 0)
                     {
-                        WarningIn("fileMonitor::~fileMonitor()")
-                            << "Failed deleting directory watch "
-                            << dirWatches_[i] << endl;
+                        if (inotify_rm_watch(inotifyFd_, int(dirWatches_[i])))
+                        {
+                            WarningIn("fileMonitor::~fileMonitor()")
+                                << "Failed deleting directory watch "
+                                << dirWatches_[i] << endl;
+                        }
                     }
                 }
             }
@@ -177,10 +201,15 @@ namespace Foam
 
         inline bool addWatch(const label watchFd, const fileName& fName)
         {
+            if (inotifyFd_ < 0)
+            {
+                return false;
+            }
+
             // Add/retrieve watch on directory containing file
             label dirWatchID = inotify_add_watch
             (
-                fd,
+                inotifyFd_,
                 fName.path().c_str(),
                 IN_CLOSE_WRITE
             );
@@ -189,7 +218,8 @@ namespace Foam
             {
                 FatalErrorIn("addWatch(const label, const fileName&)")
                     << "Failed adding watch " << watchFd
-                    << " to directory " << fName
+                    << " to directory " << fName << " due to "
+                    << string(strerror(errno))
                     << exit(FatalError);
             }
 
@@ -209,6 +239,11 @@ namespace Foam
 
         inline bool removeWatch(const label watchFd)
         {
+            if (inotifyFd_ < 0)
+            {
+                return false;
+            }
+
             dirWatches_[watchFd] = -1;
             return true;
         }
@@ -263,11 +298,11 @@ void Foam::fileMonitor::checkFiles() const
         fd_set fdSet;
         // Add notify descriptor to select fd_set
         FD_ZERO(&fdSet);
-        FD_SET(watcher_->fd, &fdSet);
+        FD_SET(watcher_->inotifyFd_, &fdSet);
 
         int ready = select
         (
-            watcher_->fd+1,     // num filedescriptors in fdSet
+            watcher_->inotifyFd_+1,     // num filedescriptors in fdSet
             &fdSet,             // fdSet with only inotifyFd
             NULL,               // No writefds
             NULL,               // No errorfds
@@ -280,15 +315,15 @@ void Foam::fileMonitor::checkFiles() const
                 << "Problem in issuing select."
                 << abort(FatalError);
         }
-        else if (FD_ISSET(watcher_->fd, &fdSet))
+        else if (FD_ISSET(watcher_->inotifyFd_, &fdSet))
         {
             // Read events
-            ssize_t nBytes = read(watcher_->fd, buffer, EVENT_BUF_LEN);
+            ssize_t nBytes = read(watcher_->inotifyFd_, buffer, EVENT_BUF_LEN);
 
             if (nBytes < 0)
             {
                 FatalErrorIn("fileMonitor::updateStates(const fileName&)")
-                    << "read of " << watcher_->fd
+                    << "read of " << watcher_->inotifyFd_
                     << " failed with " << label(nBytes)
                     << abort(FatalError);
             }
@@ -364,6 +399,7 @@ Foam::label Foam::fileMonitor::addWatch(const fileName& fName)
     label watchFd;
 
     label sz = freeWatchFds_.size();
+
     if (sz)
     {
         watchFd = freeWatchFds_[sz-1];

src/OpenFOAM/Make/files

@@ -346,6 +346,7 @@ $(basicPolyPatches)/generic/genericPolyPatch.C
 constraintPolyPatches = $(polyPatches)/constraint
 $(constraintPolyPatches)/cyclic/cyclicPolyPatch.C
 $(constraintPolyPatches)/cyclicSlip/cyclicSlipPolyPatch.C
+$(constraintPolyPatches)/oldCyclic/oldCyclicPolyPatch.C
 $(constraintPolyPatches)/empty/emptyPolyPatch.C
 $(constraintPolyPatches)/nonuniformTransformCyclic/nonuniformTransformCyclicPolyPatch.C
 $(constraintPolyPatches)/processorCyclic/processorCyclicPolyPatch.C

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/Time/Time.C

@@ -245,13 +245,14 @@ Foam::Time::Time
     readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    setControls();
+
     // Time objects not registered so do like objectRegistry::checkIn ourselves.
     if (runTimeModifiable_)
     {
+        monitorPtr_.reset(new fileMonitor());
         controlDict_.watchIndex() = addWatch(controlDict_.filePath());
     }
-
-    setControls();
 }
 
 
@@ -307,14 +308,20 @@ Foam::Time::Time
     readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    setControls();
 
     // Time objects not registered so do like objectRegistry::checkIn ourselves.
     if (runTimeModifiable_)
     {
-        controlDict_.watchIndex() = addWatch(controlDict_.filePath());
-    }
+        monitorPtr_.reset(new fileMonitor());
 
-    setControls();
+        // File might not exist yet.
+        fileName f(controlDict_.filePath());
+        if (f != fileName::null)
+        {
+            controlDict_.watchIndex() = addWatch(f);
+        }
+    }
 }
 
 
@@ -388,18 +395,18 @@ Foam::Time::~Time()
 
 Foam::label Foam::Time::addWatch(const fileName& fName) const
 {
-    return monitor_.addWatch(fName);
+    return monitorPtr_().addWatch(fName);
 }
 
 
 bool Foam::Time::removeWatch(const label watchIndex) const
 {
-    return monitor_.removeWatch(watchIndex);
+    return monitorPtr_().removeWatch(watchIndex);
 }
 
 const Foam::fileName& Foam::Time::getFile(const label watchIndex) const
 {
-    return monitor_.getFile(watchIndex);
+    return monitorPtr_().getFile(watchIndex);
 }
 
 
@@ -408,13 +415,13 @@ Foam::fileMonitor::fileState Foam::Time::getState
     const label watchFd
 ) const
 {
-    return monitor_.getState(watchFd);
+    return monitorPtr_().getState(watchFd);
 }
 
 
 void Foam::Time::setUnmodified(const label watchFd) const
 {
-    monitor_.setUnmodified(watchFd);
+    monitorPtr_().setUnmodified(watchFd);
 }
 
 

src/OpenFOAM/db/Time/Time.H

@@ -73,7 +73,7 @@ class Time
     // Private data
 
         //- file-change monitor for all registered files
-        mutable fileMonitor monitor_;
+        mutable autoPtr<fileMonitor> monitorPtr_;
 
         //- The controlDict
         IOdictionary controlDict_;

src/OpenFOAM/db/Time/TimeIO.C

@@ -211,10 +211,7 @@ void Foam::Time::readModifiedObjects()
         // valid filePath).
         // Note: requires same ordering in objectRegistries on different
         // processors!
-        monitor_.updateStates(Pstream::parRun());
-
-//Pout<< "Time : runTimeModifiable_ and watchIndex:"
-//    << controlDict_.watchIndex() << endl;
+        monitorPtr_().updateStates(Pstream::parRun());
 
         // Time handling is special since controlDict_ is the one dictionary
         // that is not registered to any database.

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/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H

@@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::oldCyclicPolyPatch
+
+Description
+    'old' style cyclic polyPatch with all faces in single patch. Does ordering
+    but cannot be used to run. Writes 'type cyclic' so foamUpgradeCyclics
+    can be run afterwards.
+    Used to get cyclics from mesh converters that assume cyclics in single
+    patch (e.g. fluent3DMeshToFoam)
+
+SourceFiles
+    oldCyclicPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oldCyclicPolyPatch_H
+#define oldCyclicPolyPatch_H
+
+#include "coupledPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class oldCyclicPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class oldCyclicPolyPatch
+:
+    public coupledPolyPatch
+{
+public:
+
+    enum transformType
+    {
+        UNKNOWN,
+        ROTATIONAL,
+        TRANSLATIONAL
+    };
+    static const NamedEnum<transformType, 3> transformTypeNames;
+
+
+private:
+
+    // Private data
+
+        //- Morph:angle between normals of neighbouring faces.
+        //  Used to split cyclic into halves.
+        scalar featureCos_;
+
+        //- Type of transformation - rotational or translational
+        transformType transform_;
+
+        // For rotation
+
+            //- Axis of rotation for rotational cyclics
+            vector rotationAxis_;
+
+            //- point on axis of rotation for rotational cyclics
+            point rotationCentre_;
+
+        // For translation
+
+            //- Translation vector
+            vector separationVector_;
+
+
+    // Private member functions
+
+        //- Find amongst selected faces the one with the largest area
+        static label findMaxArea(const pointField&, const faceList&);
+
+        void calcTransforms();
+
+        //- Calculate face centres
+        static pointField calcFaceCentres
+        (
+            const UList<face>&,
+            const pointField&
+        );
+
+        //- Get f[0] for all faces
+        static pointField getAnchorPoints
+        (
+            const UList<face>&,
+            const pointField&
+        );
+
+        // Face ordering
+
+            //- Find the two parts of the faces of pp using feature edges.
+            //  Returns true if successfull.
+            bool getGeometricHalves
+            (
+                const primitivePatch&,
+                labelList&,
+                labelList&
+            ) const;
+
+            //- Calculate geometric factors of the two halves.
+            void getCentresAndAnchors
+            (
+                const primitivePatch&,
+                const faceList& half0Faces,
+                const faceList& half1Faces,
+
+                pointField& ppPoints,
+                pointField& half0Ctrs,
+                pointField& half1Ctrs,
+                pointField& anchors0,
+                scalarField& tols
+            ) const;
+
+            //- Given matched faces matches the anchor point. Sets faceMap,
+            //  rotation. Returns true if all matched.
+            bool matchAnchors
+            (
+                const bool report,
+                const primitivePatch&,
+                const labelList&,
+                const pointField&,
+                const labelList&,
+                const faceList&,
+                const labelList&,
+                const scalarField&,
+
+                labelList& faceMap,
+                labelList& rotation
+            ) const;
+
+            //- For rotational cases, try to find a unique face on each side
+            //  of the cyclic.
+            label getConsistentRotationFace
+            (
+                const pointField& faceCentres
+            ) const;
+
+
+protected:
+
+    // Protected Member functions
+
+        //- Initialise the calculation of the patch geometry
+        virtual void initGeometry(PstreamBuffers&);
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry(PstreamBuffers&);
+
+        //- Initialise the patches for moving points
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
+
+        //- Correct patches after moving points
+        virtual void movePoints(PstreamBuffers&, const pointField&);
+
+        //- Initialise the update of the patch topology
+        virtual void initUpdateMesh(PstreamBuffers&);
+
+        //- Update of the patch topology
+        virtual void updateMesh(PstreamBuffers&);
+
+public:
+
+    //- Runtime type information
+    TypeName("oldCyclic");
+
+
+    // Constructors
+
+        //- Construct from components
+        oldCyclicPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from dictionary
+        oldCyclicPolyPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct as copy, resetting the boundary mesh
+        oldCyclicPolyPatch(const oldCyclicPolyPatch&, const polyBoundaryMesh&);
+
+        //- Construct given the original patch and resetting the
+        //  face list and boundary mesh information
+        oldCyclicPolyPatch
+        (
+            const oldCyclicPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        );
+
+        //- Construct and return a clone, resetting the boundary mesh
+        virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
+        {
+            return autoPtr<polyPatch>(new oldCyclicPolyPatch(*this, bm));
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new oldCyclicPolyPatch(*this, bm, index, newSize, newStart)
+            );
+        }
+
+
+    // Destructor
+
+        virtual ~oldCyclicPolyPatch();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Does this side own the patch ?
+            virtual bool owner() const
+            {
+                notImplemented("oldCyclicPolyPatch::owner()");
+                return true;
+            }
+
+            //- Transform a patch-based position from other side to this side
+            virtual void transformPosition(pointField& l) const
+            {
+                notImplemented("transformPosition(pointField&)");
+            }
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry
+        (
+            const primitivePatch& referPatch,
+            const UList<point>& thisCtrs,
+            const UList<point>& thisAreas,
+            const UList<point>& thisCc,
+            const UList<point>& nbrCtrs,
+            const UList<point>& nbrAreas,
+            const UList<point>& nbrCc
+        );
+
+        //- Initialize ordering for primitivePatch. Does not
+        //  refer to *this (except for name() and type() etc.)
+        virtual void initOrder
+        (
+            PstreamBuffers&,
+            const primitivePatch&
+        ) const;
+
+        //- Return new ordering for primitivePatch.
+        //  Ordering is -faceMap: for every face
+        //  index of the new face -rotation:for every new face the clockwise
+        //  shift of the original face. Return false if nothing changes
+        //  (faceMap is identity, rotation is 0), true otherwise.
+        virtual bool order
+        (
+            PstreamBuffers&,
+            const primitivePatch&,
+            labelList& faceMap,
+            labelList& rotation
+        ) const;
+
+        //- Write the polyPatch data as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/meshCut/directions/directions.C

@@ -32,6 +32,7 @@ License
 #include "meshTools.H"
 #include "hexMatcher.H"
 #include "Switch.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -206,7 +207,7 @@ Foam::vectorField Foam::directions::propagateDirection
         mesh,
         changedFaces,
         changedFacesInfo,
-        mesh.nCells()
+        mesh.globalData().nTotalCells()+1
     );
 
     const List<directionInfo>& cellInfo = directionCalc.allCellInfo();

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -2332,7 +2332,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
         seedFacesInfo.clear();
 
         // Iterate until no change. Now 2:1 face difference should be satisfied
-        levelCalc.iterate(mesh_.globalData().nTotalFaces());
+        levelCalc.iterate(mesh_.globalData().nTotalFaces()+1);
 
 
         // Now check point-connected cells (face-connected cells already ok):
@@ -2836,7 +2836,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
         seedFacesInfo,
         allFaceInfo,
         allCellInfo,
-        mesh_.globalData().nTotalCells()
+        mesh_.globalData().nTotalCells()+1
     );
 
 

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/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C

@@ -171,6 +171,23 @@ tmp<Field<Type> > cyclicFvPatchField<Type>::patchNeighbourField() const
 }
 
 
+template<class Type>
+const cyclicFvPatchField<Type>& cyclicFvPatchField<Type>::neighbourPatchField()
+const
+{
+    const GeometricField<Type, fvPatchField, volMesh>& fld =
+    static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
+    (
+        this->internalField()
+    );
+
+    return refCast<const cyclicFvPatchField<Type> >
+    (
+        fld.boundaryField()[this->cyclicPatch().neighbPatchID()]
+    );
+}
+
+
 template<class Type>
 void cyclicFvPatchField<Type>::updateInterfaceMatrix
 (

src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H

@@ -150,9 +150,12 @@ public:
 
         // Evaluation functions
 
-            //- Return neighbour coupled given internal cell data
+            //- Return neighbour coupled internal cell data
             tmp<Field<Type> > patchNeighbourField() const;
 
+            //- Return reference to neighbour patchField
+            const cyclicFvPatchField<Type>& neighbourPatchField() const;
+
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (

src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.H

@@ -25,7 +25,13 @@ Class
     Foam::directMappedFixedValueFvPatchField
 
 Description
-    Recycles the value at a set of internal faces back to *this.
+    Recycles the value at a set of cells or patch faces back to *this. Can not
+    sample internal faces (since volField not defined on faces).
+
+    mode = NEARESTCELL : sample nearest cell
+    mode = NEARESTPATCHFACE : sample nearest face on selected patch
+    mode = NEARESTFACE : sample nearest face on any patch. Note: does not
+                         warn if nearest actually is on internal face!
 
 SourceFiles
     directMappedFixedValueFvPatchField.C

src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.H

@@ -27,6 +27,11 @@ Class
 Description
     Recycles the velocity and flux at a patch to this patch
 
+    mode = NEARESTCELL : sample nearest cell
+    mode = NEARESTPATCHFACE : sample nearest face on selected patch
+    mode = NEARESTFACE : sample nearest face on any patch. Note: does not
+                         warn if nearest actually is on internal face!
+
 SourceFiles
     directMappedVelocityFluxFixedValueFvPatchField.C
 

src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H

@@ -140,7 +140,17 @@ public:
             //- Return the "jump" across the patch.
             virtual tmp<Field<Type> > jump() const
             {
-                return jump_;
+                if (this->cyclicPatch().owner())
+                {
+                    return jump_;
+                }
+                else
+                {
+                    return refCast<const fanFvPatchField<Type> >
+                    (
+                        this->neighbourPatchField()
+                    ).jump();
+                }
             }
 
 

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/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H

@@ -108,7 +108,7 @@ public:
             }
 
             //- Return processor number
-            virtual const cyclicLduInterface& neighbPatch() const
+            virtual const cyclicFvPatch& neighbPatch() const
             {
                 return refCast<const cyclicFvPatch>
                 (

src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C

@@ -115,7 +115,7 @@ void Foam::inverseFaceDistanceDiffusivity::correct()
         mesh,
         changedFaces,
         faceDist,
-        mesh.globalData().nTotalCells() // max iterations
+        mesh.globalData().nTotalCells()+1 // max iterations
     );
 
     const List<wallPoint>& faceInfo = waveInfo.allFaceInfo();

src/meshTools/cellClassification/cellClassification.C

@@ -322,7 +322,7 @@ void Foam::cellClassification::markCells
         changedFaces,                       // Labels of changed faces
         changedFacesInfo,                   // Information on changed faces
         cellInfoList,                       // Information on all cells
-        mesh_.globalData().nTotalCells()    // max iterations
+        mesh_.globalData().nTotalCells()+1  // max iterations
     );
 
     // Get information out of cellInfoList

src/meshTools/cellDist/patchWave/patchDataWave.C

@@ -225,7 +225,7 @@ void Foam::patchDataWave<TransferType>::correct()
         mesh(),
         changedFaces,
         faceDist,
-        mesh().globalData().nTotalCells() // max iterations
+        mesh().globalData().nTotalCells()+1 // max iterations
     );
 
 

src/meshTools/cellDist/patchWave/patchWave.C

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

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H

@@ -66,9 +66,9 @@ public:
         //- Mesh items to sample
         enum sampleMode
         {
-            NEARESTCELL,
-            NEARESTPATCHFACE,
-            NEARESTFACE
+            NEARESTCELL,        // nearest cell
+            NEARESTPATCHFACE,   // faces on selected patch
+            NEARESTFACE         // nearest face
         };
 
 private:
@@ -201,7 +201,7 @@ public:
             return sampleRegion_;
         }
 
-        //- Patch (only if NEARESTBOUNDARY)
+        //- Patch (only if NEARESTPATCHFACE)
         const word& samplePatch() const
         {
             return samplePatch_;

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/meshTools/sets/cellSources/boxToCell/boxToCell.H

@@ -102,6 +102,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/cellToCell/cellToCell.H

@@ -93,6 +93,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/cylinderToCell/cylinderToCell.H

@@ -109,6 +109,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/faceToCell/faceToCell.H

@@ -119,6 +119,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.H

@@ -116,6 +116,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/fieldToCell/fieldToCell.H

@@ -115,6 +115,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/labelToCell/labelToCell.H

@@ -100,6 +100,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/nbrToCell/nbrToCell.H

@@ -101,6 +101,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/nearestToCell/nearestToCell.H

@@ -100,6 +100,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/pointToCell/pointToCell.H

@@ -115,6 +115,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/regionToCell/regionToCell.H

@@ -105,6 +105,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet(const topoSetSource::setAction action, topoSet&)
          const;
 

src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.H

@@ -111,6 +111,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/shapeToCell/shapeToCell.H

@@ -109,6 +109,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/sphereToCell/sphereToCell.H

@@ -105,6 +105,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.H

@@ -195,6 +195,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellSources/zoneToCell/zoneToCell.H

@@ -101,6 +101,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLSETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/cellZoneSources/setToCellZone/setToCellZone.H

@@ -93,6 +93,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return CELLZONESOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/faceSources/boundaryToFace/boundaryToFace.H

@@ -92,6 +92,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return FACESETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/faceSources/boxToFace/boxToFace.H

@@ -101,6 +101,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return FACESETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,

src/meshTools/sets/faceSources/cellToFace/cellToFace.H

@@ -120,6 +120,11 @@ public:
 
     // Member Functions
 
+        virtual sourceType setType() const
+        {
+            return FACESETSOURCE;
+        }
+
         virtual void applyToSet
         (
             const topoSetSource::setAction action,