Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/rhoReactingFoam/Make/files

@@ -0,0 +1,3 @@
+rhoReactingFoam.C
+
+EXE = $(FOAM_APPBIN)/rhoReactingFoam

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -0,0 +1,19 @@
+EXE_INC = \
+    -I../XiFoam \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lreactionThermophysicalModels \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lchemistryModel \
+    -lODE \
+    -lfiniteVolume

applications/solvers/combustion/rhoReactingFoam/YEqn.H

@@ -0,0 +1,43 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    for (label i=0; i<Y.size(); i++)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                kappa*chemistry.RR(i),
+                mesh.solver("Yi")
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/combustion/rhoReactingFoam/chemistry.H

@@ -0,0 +1,24 @@
+{
+    Info << "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaT().value(),
+        runTime.deltaT().value()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        volScalarField tk =
+                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        volScalarField tc = chemistry.tc();
+
+        // Chalmers PaSR model
+        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+}

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -0,0 +1,85 @@
+Info<< nl << "Reading thermophysicalProperties" << endl;
+autoPtr<rhoChemistryModel> pChemistry
+(
+    rhoChemistryModel::New(mesh)
+);
+rhoChemistryModel& chemistry = pChemistry();
+
+hReactionThermo& thermo = chemistry.thermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+word inertSpecie(thermo.lookup("inertSpecie"));
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+volScalarField& h = thermo.h();
+
+
+#include "compressibleCreatePhi.H"
+
+volScalarField kappa
+(
+    IOobject
+    (
+        "kappa",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimless, 0.0)
+);
+
+Info << "Creating turbulence model.\n" << nl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+Info<< "Creating field DpDt\n" << endl;
+volScalarField DpDt =
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll (Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(h);
+

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -0,0 +1,93 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rUA = 1.0/UEqn.A();
+    U = rUA*UEqn.H();
+
+    if (transonic)
+    {
+        surfaceScalarField phiv =
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi);
+
+        phi = fvc::interpolate(rho)*phiv;
+
+        surfaceScalarField phid
+        (
+            "phid",
+            fvc::interpolate(thermo.psi())*phiv
+        );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::ddt(rho) + fvc::div(phi)
+              + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
+            {
+                pEqn.solve(mesh.solver(p.name() + "Final"));
+            }
+            else
+            {
+                pEqn.solve();
+            }
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+    else
+    {
+        phi =
+            fvc::interpolate(rho)
+           *(
+                (fvc::interpolate(U) & mesh.Sf())
+              + fvc::ddtPhiCorr(rUA, rho, U, phi)
+            );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+              + fvc::div(phi)
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
+            {
+                pEqn.solve(mesh.solver(p.name() + "Final"));
+            }
+            else
+            {
+                pEqn.solve();
+            }
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+
+    U -= rUA*fvc::grad(p);
+    U.correctBoundaryConditions();
+
+    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+}

applications/solvers/combustion/rhoReactingFoam/readChemistryProperties.H

@@ -0,0 +1,23 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rhoReactingFoam
+
+Description
+    Chemical reaction code using density based thermodynamics package.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hReactionThermo.H"
+#include "turbulenceModel.H"
+#include "rhoChemistryModel.H"
+#include "chemistrySolver.H"
+#include "multivariateScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
+
+        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        {
+            #include "YEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        turbulence->correct();
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -9,35 +9,10 @@
       + parcels.SU()
     );
 
-    UEqn.relax();
-
-    tmp<volScalarField> trAU;
-    tmp<volTensorField> trTU;
-
-    if (pressureImplicitPorosity)
-    {
-        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
-        pZones.addResistance(UEqn, tTU());
-        trTU = inv(tTU());
-        trTU().rename("rAU");
-
-        volVectorField gradp = fvc::grad(p);
-
-        for (int UCorr=0; UCorr<nUCorr; UCorr++)
-        {
-            U = trTU() & (UEqn.H() - gradp);
-        }
-        U.correctBoundaryConditions();
-    }
-    else
-    {
     pZones.addResistance(UEqn);
 
-        trAU = 1.0/UEqn.A();
-        trAU().rename("rAU");
-
-        solve
-        (
-            UEqn == -fvc::grad(p)
-        );
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
     }
+

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -43,5 +43,4 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
-
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -73,13 +73,7 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt
-    (
-        "DpDt",
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
-    );
-
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll (Y, i)

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H

@@ -1,21 +1,3 @@
     Info<< "Creating porous zones" << nl << endl;
 
     porousZones pZones(mesh);
-    Switch pressureImplicitPorosity(false);
-
-    label nUCorr = 0;
-    if (pZones.size())
-    {
-        // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
-        {
-            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
-                >> nUCorr;
-        }
-
-        if (nUCorr > 0)
-        {
-            pressureImplicitPorosity = true;
-        }
-    }
-

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H

@@ -1,20 +1,51 @@
 {
-    fvScalarMatrix hEqn
+    tmp<volScalarField> pWork
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "pWork",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
+        )
+    );
+
+    if (dpdt)
+    {
+        pWork() += fvc::ddt(p);
+    }
+    if (eWork)
+    {
+        pWork() = -p*fvc::div(phi/fvc::interpolate(rho));
+    }
+    if (hWork)
+    {
+        pWork() += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
+    }
+
+    {
+        solve
         (
             fvm::ddt(rho, h)
           + mvConvection->fvmDiv(phi, h)
           - fvm::laplacian(turbulence->alphaEff(), h)
          ==
-        DpDt
+            pWork()
           + parcels.Sh()
           + radiation->Sh(thermo)
         );
 
-    hEqn.relax();
-
-    hEqn.solve();
-
         thermo.correct();
 
         radiation->correct();
+
+        Info<< "T gas min/max   = " << min(T).value() << ", "
+            << max(T).value() << endl;
+    }
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -5,65 +5,13 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    if (pressureImplicitPorosity)
-    {
-        U = trTU()&UEqn.H();
-    }
-    else
-    {
-        U = trAU()*UEqn.H();
-    }
+    volScalarField rAU = 1.0/UEqn.A();
+    U = rAU*UEqn.H();
 
-    if (transonic)
+    if (pZones.size() > 0)
     {
-        surfaceScalarField phiv = fvc::interpolate(U) & mesh.Sf();
-
-        phi = fvc::interpolate(rho)*phiv;
-
-        surfaceScalarField phid("phid", fvc::interpolate(psi)*phiv);
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            tmp<fvScalarMatrix> lapTerm;
-
-            if (pressureImplicitPorosity)
-            {
-                lapTerm = fvm::laplacian(rho*trTU(), p);
-            }
-            else
-            {
-                lapTerm = fvm::laplacian(rho*trAU(), p);
-            }
-
-            fvScalarMatrix pEqn
-            (
-                fvc::ddt(rho) + fvc::div(phi)
-              + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-              - lapTerm()
-            ==
-                parcels.Srho()
-              + pointMassSources.Su()
-            );
-
-            if
-            (
-                oCorr == nOuterCorr-1
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-            {
-                pEqn.solve(mesh.solver("pFinal"));
-            }
-            else
-            {
-                pEqn.solve();
-            }
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi == pEqn.flux();
-            }
-        }
+        // ddtPhiCorr not well defined for cases with porosity
+        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     }
     else
     {
@@ -71,38 +19,23 @@
             fvc::interpolate(rho)
            *(
                 (fvc::interpolate(U) & mesh.Sf())
-//            + fvc::ddtPhiCorr(trAU(), rho, U, phi)
+              + fvc::ddtPhiCorr(rAU, rho, U, phi)
             );
+    }
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-            tmp<fvScalarMatrix> lapTerm;
-
-            if (pressureImplicitPorosity)
-            {
-                lapTerm = fvm::laplacian(rho*trTU(), p);
-            }
-            else
-            {
-                lapTerm = fvm::laplacian(rho*trAU(), p);
-            }
-
         fvScalarMatrix pEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phi)
-              - lapTerm()
+          - fvm::laplacian(rho*rAU, p)
         ==
             parcels.Srho()
           + pointMassSources.Su()
         );
 
-            if
-            (
-                oCorr == nOuterCorr-1
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
             pEqn.solve(mesh.solver("pFinal"));
         }
@@ -116,7 +49,6 @@
             phi += pEqn.flux();
         }
     }
-    }
 
     // Second part of thermodynamic density update
     thermo.rho() += psi*p;
@@ -124,15 +56,6 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    if (pressureImplicitPorosity)
-    {
-        U -= trTU()&fvc::grad(p);
-    }
-    else
-    {
-        U -= trAU()*fvc::grad(p);
-    }
+    U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
-
-    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -23,15 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    trackedReactingParcelFoam
+    porousExplicitSourceReactingParcelFoam
 
 Description
-    - reacting parcel cloud tracking
+    - reacting parcel cloud
     - porous media
     - point mass sources
     - polynomial based, incompressible thermodynamics (f(T))
 
-    Note: ddtPhiCorr not used here - not well defined for porous calcs
+    Note: ddtPhiCorr not used here when porous zones are active
+    - not well defined for porous calcs
 
 \*---------------------------------------------------------------------------*/
 
@@ -74,6 +75,7 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
         #include "readPISOControls.H"
+        #include "readAdditionalSolutionControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -88,22 +90,15 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
-
-        // --- PIMPLE loop
-        for (int oCorr=1; oCorr<=nOuterCorr; oCorr++)
-        {
         #include "YEqn.H"
         #include "hEqn.H"
 
         // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+        for (int corr=0; corr<nCorr; corr++)
         {
             #include "pEqn.H"
         }
 
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
-        }
 
         turbulence->correct();
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H

@@ -0,0 +1,20 @@
+dictionary additional = mesh.solutionDict().subDict("additional");
+
+bool dpdt = true;
+if (additional.found("dpdt"))
+{
+    additional.lookup("dpdt") >> dpdt;
+}
+
+bool eWork = true;
+if (additional.found("eWork"))
+{
+    additional.lookup("eWork") >> eWork;
+}
+
+bool hWork = true;
+if (additional.found("hWork"))
+{
+    additional.lookup("hWork") >> hWork;
+}
+

applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C

@@ -78,6 +78,7 @@ int main(int argc, char *argv[])
     if (sourceType == "patch")
     {
         fileName sourceCasePath(dict.lookup("sourceCase"));
+        sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
         word patchName(dict.lookup("sourcePatch"));

applications/utilities/mesh/generation/extrudeMesh/setRoots.H

@@ -1,16 +0,0 @@
-    argList::validArgs.clear();
-    argList::noParallel();
-
-    argList::validOptions.insert("sourceCase", "source case");
-    argList::validOptions.insert("sourcePatch", "source patch");
-
-    argList::validOptions.insert("surface", "surface file");
-
-    argList::validOptions.insert("mergeFaces", "");
-
-    argList args(argc, argv);
-
-    if (!args.check())
-    {
-         FatalError.exit();
-    }

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -209,6 +209,8 @@ snapControls
 // Settings for the layer addition.
 addLayersControls
 {
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
     // Per final patch (so not geometry!) the layer information
@@ -230,14 +232,13 @@ addLayersControls
 
 
     //- Wanted thickness of final added cell layer. If multiple layers
-    //  is the
-    //  thickness of the layer furthest away from the wall.
-    //  Relative to undistorted size of cell outside layer.
+    //  is the thickness of the layer furthest away from the wall.
+    //  See relativeSizes parameter.
     finalLayerThickness 0.3;
 
     //- Minimum thickness of cell layer. If for any reason layer
     //  cannot be above minThickness do not add layer.
-    //  Relative to undistorted size of cell outside layer.
+    //  See relativeSizes parameter.
     minThickness 0.25;
 
     //- If points get not extruded do nGrow layers of connected faces that are

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C

@@ -38,7 +38,7 @@ Description
 #include "Time.H"
 #include "timeSelector.H"
 #include "OFstream.H"
-#include "passiveParticle.H"
+#include "passiveParticleCloud.H"
 
 using namespace Foam;
 
@@ -61,7 +61,8 @@ int main(int argc, char *argv[])
     fileName vtkPath(runTime.path()/"VTK");
     mkDir(vtkPath);
 
-    Info<< "Scanning times to determine track data" << nl << endl;
+    Info<< "Scanning times to determine track data for cloud " << cloudName
+        << nl << endl;
 
     labelList maxIds(Pstream::nProcs(), -1);
     forAll(timeDirs, timeI)
@@ -69,23 +70,12 @@ int main(int argc, char *argv[])
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
 
-        IOobject positionsHeader
-        (
-            "positions",
-            runTime.timeName(),
-            cloud::prefix/cloudName,
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        );
-
-        if (positionsHeader.headerOk())
-        {
         Info<< "    Reading particle positions" << endl;
-            Cloud<passiveParticle> myCloud(mesh, cloudName, false);
+        passiveParticleCloud myCloud(mesh, cloudName);
+        Info<< "    Read " << returnReduce(myCloud.size(), sumOp<label>())
+            << " particles" << endl;
 
-            forAllConstIter(Cloud<passiveParticle>, myCloud, iter)
+        forAllConstIter(passiveParticleCloud, myCloud, iter)
         {
             label origId = iter().origId();
             label origProc = iter().origProc();
@@ -93,11 +83,20 @@ int main(int argc, char *argv[])
             maxIds[origProc] = max(maxIds[origProc], origId);
         }
     }
-    }
-    Pstream::listCombineGather(maxIds, maxOp<label>());
+    Pstream::listCombineGather(maxIds, maxEqOp<label>());
     Pstream::listCombineScatter(maxIds);
+
     labelList numIds = maxIds + 1;
 
+    Info<< nl << "Particle statistics:" << endl;
+    forAll(maxIds, procI)
+    {
+        Info<< "    Found " << numIds[procI] << " particles originating"
+            << " from processor " << procI << endl;
+    }
+    Info<< nl << endl;
+
+
     // calc starting ids for particles on each processor
     List<label> startIds(numIds.size(), 0);
     for (label i = 0; i < numIds.size()-1; i++)
@@ -106,56 +105,49 @@ int main(int argc, char *argv[])
     }
     label nParticle = startIds[startIds.size()-1] + numIds[startIds.size()-1];
 
+
+
     // number of tracks to generate
     label nTracks = nParticle/sampleFrequency;
 
     // storage for all particle tracks
     List<DynamicList<vector> > allTracks(nTracks);
 
-    Info<< "\nGenerating " << nTracks << " particle tracks" << nl << endl;
+    Info<< "\nGenerating " << nTracks << " particle tracks for cloud "
+        << cloudName << nl << endl;
+
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
 
-        IOobject positionsHeader
-        (
-            "positions",
-            runTime.timeName(),
-            cloud::prefix/cloudName,
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        );
-
-        if (positionsHeader.headerOk())
-        {
-            Info<< "    Reading particle positions" << endl;
-            Cloud<passiveParticle> myCloud(mesh, cloudName, false);
-
-            // collect the track data on the master processor
         List<pointField> allPositions(Pstream::nProcs());
+        List<labelField> allOrigIds(Pstream::nProcs());
+        List<labelField> allOrigProcs(Pstream::nProcs());
+
+        // Read particles. Will be size 0 if no particles.
+        Info<< "    Reading particle positions" << endl;
+        passiveParticleCloud myCloud(mesh, cloudName);
+
+        // collect the track data on all processors that have positions
         allPositions[Pstream::myProcNo()].setSize
         (
             myCloud.size(),
             point::zero
         );
-
-            List<labelField> allOrigIds(Pstream::nProcs());
         allOrigIds[Pstream::myProcNo()].setSize(myCloud.size(), 0);
-
-            List<labelField> allOrigProcs(Pstream::nProcs());
         allOrigProcs[Pstream::myProcNo()].setSize(myCloud.size(), 0);
 
         label i = 0;
-            forAllConstIter(Cloud<passiveParticle>, myCloud, iter)
+        forAllConstIter(passiveParticleCloud, myCloud, iter)
         {
             allPositions[Pstream::myProcNo()][i] = iter().position();
             allOrigIds[Pstream::myProcNo()][i] = iter().origId();
             allOrigProcs[Pstream::myProcNo()][i] = iter().origProc();
             i++;
         }
+
+        // collect the track data on the master processor
         Pstream::gatherList(allPositions);
         Pstream::gatherList(allOrigIds);
         Pstream::gatherList(allOrigProcs);
@@ -186,18 +178,14 @@ int main(int argc, char *argv[])
             }
         }
     }
-        else
-        {
-            Info<< "    No particles read" << nl << endl;
-        }
-    }
 
     if (Pstream::master())
     {
-        Info<< "\nWriting particle tracks" << nl << endl;
-
         OFstream vtkTracks(vtkPath/"particleTracks.vtk");
 
+        Info<< "\nWriting particle tracks to " << vtkTracks.name()
+            << nl << endl;
+
         // Total number of points in tracks + 1 per track
         label nPoints = 0;
         forAll(allTracks, trackI)

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

@@ -63,17 +63,28 @@ class OutputFilterFunctionObject
 {
     // Private data
 
+        //- Output filter name
         word name_;
+
+        //- Reference to the time database
         const Time& time_;
+
+        //- Input dictionary
         dictionary dict_;
+
+        //- Name of region
         word regionName_;
+
+        //- Optional dictionary name to supply required inputs
         word dictName_;
 
         //- Switch for the execution of the functionObject
         bool enabled_;
 
+        //- Output controls
         outputFilterOutputControl outputControl_;
 
+        //- Pointer to the output filter
         autoPtr<OutputFilter> ptr_;
 
 
@@ -108,12 +119,59 @@ public:
 
     // Member Functions
 
+        // Access
+
             //- Return name
             virtual const word& name() const
             {
                 return name_;
             }
 
+            //- Return time database
+            virtual const Time& time() const
+            {
+                return time_;
+            }
+
+            //- Return the input dictionary
+            virtual const dictionary& dict() const
+            {
+                return dict_;
+            }
+
+            //- Return the region name
+            virtual const word& regionName() const
+            {
+                return regionName_;
+            }
+
+            //- Return the optional dictionary name
+            virtual const word& dictName() const
+            {
+                return dictName_;
+            }
+
+            //- Return the enabled flag
+            virtual bool enabled() const
+            {
+                return enabled_;
+            }
+
+            //- Return the output control object
+            virtual const outputFilterOutputControl& outputControl() const
+            {
+                return outputControl_;
+            }
+
+            //- Return the output filter
+            virtual const OutputFilter& outputFilter() const
+            {
+                return ptr_();
+            }
+
+
+        // Function object control
+
             //- Switch the function object on
             virtual void on();
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -113,7 +113,7 @@ Foam::label Foam::meshRefinement::createBaffle
                 true,                       // face flip
                 neiPatch,                   // patch for face
                 zoneID,                     // zone for face
-                zoneFlip                    // face flip in zone
+                !zoneFlip                   // face flip in zone
             )
         );
     }

src/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C

@@ -56,7 +56,7 @@ Foam::trackedParticle::trackedParticle
     bool readFields
 )
 :
-    ExactParticle<trackedParticle>(c, is)
+    ExactParticle<trackedParticle>(c, is, readFields)
 {
     if (readFields)
     {

src/lagrangian/basic/Cloud/Cloud.C

@@ -42,15 +42,12 @@ Foam::Cloud<ParticleType>::Cloud
 )
 :
     cloud(pMesh),
-    IDLList<ParticleType>(particles),
+    IDLList<ParticleType>(),
     polyMesh_(pMesh),
-    allFaces_(pMesh.faces()),
-    points_(pMesh.points()),
-    cellFaces_(pMesh.cells()),
-    allFaceCentres_(pMesh.faceCentres()),
-    owner_(pMesh.faceOwner()),
-    neighbour_(pMesh.faceNeighbour())
-{}
+    particleCount_(0)
+{
+    IDLList<ParticleType>::operator=(particles);
+}
 
 
 template<class ParticleType>
@@ -62,19 +59,31 @@ Foam::Cloud<ParticleType>::Cloud
 )
 :
     cloud(pMesh, cloudName),
-    IDLList<ParticleType>(particles),
+    IDLList<ParticleType>(),
     polyMesh_(pMesh),
-    allFaces_(pMesh.faces()),
-    points_(pMesh.points()),
-    cellFaces_(pMesh.cells()),
-    allFaceCentres_(pMesh.faceCentres()),
-    owner_(pMesh.faceOwner()),
-    neighbour_(pMesh.faceNeighbour())
-{}
+    particleCount_(0)
+{
+    IDLList<ParticleType>::operator=(particles);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ParticleType>
+Foam::label Foam::Cloud<ParticleType>::getNewParticleID() const
+{
+    label id = particleCount_++;
+
+    if (id == labelMax)
+    {
+        WarningIn("Cloud<ParticleType>::getNewParticleID() const")
+            << "Particle counter has overflowed. This might cause problems"
+            << " when reconstructing particle tracks." << endl;
+    }
+    return id;
+}
+
+
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::addParticle(ParticleType* pPtr)
 {

src/lagrangian/basic/Cloud/Cloud.H

@@ -50,6 +50,9 @@ namespace Foam
 template<class ParticleType>
 class Cloud;
 
+template<class ParticleType>
+class IOPosition;
+
 template<class ParticleType>
 Ostream& operator<<
 (
@@ -71,12 +74,9 @@ class Cloud
     // Private data
 
         const polyMesh& polyMesh_;
-        const faceList& allFaces_;
-        const vectorField& points_;
-        const cellList& cellFaces_;
-        const vectorField& allFaceCentres_;
-        const unallocLabelList& owner_;
-        const unallocLabelList& neighbour_;
+
+        //- Overall count of particles ever created. Never decreases.
+        mutable label particleCount_;
 
         //- Temporary storage for addressing. Used in findFaces.
         mutable DynamicList<label> labels_;
@@ -92,6 +92,8 @@ public:
 
     template<class ParticleT>
     friend class Particle;
+    template<class ParticleT>
+    friend class IOPosition;
 
     typedef ParticleType particleType;
 
@@ -218,6 +220,9 @@ public:
                 return IDLList<ParticleType>::clear();
             };
 
+            //- Get unique particle creation id
+            label getNewParticleID() const;
+
             //- Transfer particle to cloud
             void addParticle(ParticleType* pPtr);
 
@@ -251,13 +256,15 @@ public:
                 const IOField<DataType>& data
             ) const;
 
-            //- Read the field data for the cloud of particles
-            void readFields();
+            //- Read the field data for the cloud of particles. Dummy at
+            //  this level.
+            virtual void readFields();
 
 
         // Write
 
-            //- Write the field data for the cloud of particles
+            //- Write the field data for the cloud of particles Dummy at
+            //  this level.
             virtual void writeFields() const;
 
             //- Write using given format, version and compression.

src/lagrangian/basic/Cloud/CloudIO.C

@@ -67,12 +67,7 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh),
     polyMesh_(pMesh),
-    allFaces_(pMesh.faces()),
-    points_(pMesh.points()),
-    cellFaces_(pMesh.cells()),
-    allFaceCentres_(pMesh.faceCentres()),
-    owner_(pMesh.faceOwner()),
-    neighbour_(pMesh.faceNeighbour())
+    particleCount_(0)
 {
     initCloud(checkClass);
 }
@@ -88,12 +83,7 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh, cloudName),
     polyMesh_(pMesh),
-    allFaces_(pMesh.faces()),
-    points_(pMesh.points()),
-    cellFaces_(pMesh.cells()),
-    allFaceCentres_(pMesh.faceCentres()),
-    owner_(pMesh.faceOwner()),
-    neighbour_(pMesh.faceNeighbour())
+    particleCount_(0)
 {
     initCloud(checkClass);
 }

src/lagrangian/basic/IOPosition/IOPosition.C

@@ -34,6 +34,7 @@ Foam::word Foam::IOPosition<ParticleType>::particlePropertiesName
     "particleProperties"
 );
 
+
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
 template<class ParticleType>
@@ -58,7 +59,7 @@ void Foam::IOPosition<ParticleType>::readParticleProperties()
         if (propsDict.found(procName))
         {
             propsDict.subDict(procName).lookup("particleCount")
-                >> Particle<ParticleType>::particleCount;
+                >> cloud_.particleCount_;
         }
     }
 }
@@ -83,11 +84,7 @@ void Foam::IOPosition<ParticleType>::writeParticleProperties() const
 
     word procName("processor" + Foam::name(Pstream::myProcNo()));
     propsDict.add(procName, dictionary());
-    propsDict.subDict(procName).add
-    (
-        "particleCount",
-        Particle<ParticleType>::particleCount
-    );
+    propsDict.subDict(procName).add("particleCount", cloud_.particleCount_);
 
     propsDict.regIOobject::write();
 }
@@ -135,13 +132,20 @@ bool Foam::IOPosition<ParticleType>::write() const
 template<class ParticleType>
 bool Foam::IOPosition<ParticleType>::writeData(Ostream& os) const
 {
+    // Write global cloud data
     writeParticleProperties();
 
     os<< cloud_.size() << nl << token::BEGIN_LIST << nl;
 
     forAllConstIter(typename Cloud<ParticleType>, cloud_, iter)
     {
-        os<< static_cast<const Particle<ParticleType>&>(iter()) << nl;
+        // Prevent writing additional fields
+        static_cast<const Particle<ParticleType>&>(iter()).write
+        (
+            os,
+            false
+        );
+        os  << nl;
     }
 
     os<< token::END_LIST << endl;
@@ -157,6 +161,7 @@ void Foam::IOPosition<ParticleType>::readData
     bool checkClass
 )
 {
+    // Read global cloud data. Resets count on cloud.
     readParticleProperties();
 
     Istream& is = readStream(checkClass ? typeName : "");
@@ -172,6 +177,7 @@ void Foam::IOPosition<ParticleType>::readData
 
         for (label i=0; i<s; i++)
         {
+            // Do not read any fields, position only
             c.append(new ParticleType(c, is, false));
         }
 
@@ -202,6 +208,7 @@ void Foam::IOPosition<ParticleType>::readData
         )
         {
             is.putBack(lastToken);
+            // Do not read any fields, position only
             c.append(new ParticleType(c, is, false));
             is >> lastToken;
         }

src/lagrangian/basic/IOPosition/IOPosition.H

@@ -26,7 +26,7 @@ Class
     Foam::IOPosition
 
 Description
-    Helper IO class to write particle positions
+    Helper IO class to read and write particle positions
 
 SourceFiles
     IOPosition.C
@@ -72,7 +72,7 @@ public:
 
     // Static data
 
-        //- Runtime type name information
+        //- Runtime type name information. Use cloud type.
         virtual const word& type() const
         {
             return cloud_.type();
@@ -90,11 +90,11 @@ public:
 
     // Member functions
 
-        void readData(Cloud<ParticleType>& c, bool checkClass);
+        virtual void readData(Cloud<ParticleType>& c, bool checkClass);
 
-        bool write() const;
+        virtual bool write() const;
 
-        bool writeData(Ostream& os) const;
+        virtual bool writeData(Ostream& os) const;
 };
 
 

src/lagrangian/basic/Particle/Particle.C

@@ -33,12 +33,6 @@ License
 #include "wallPolyPatch.H"
 #include "transform.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-template<class ParticleType>
-Foam::label Foam::Particle<ParticleType>::particleCount = 0;
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class ParticleType>
@@ -75,7 +69,7 @@ void Foam::Particle<ParticleType>::findFaces
     DynamicList<label>& faceList
 ) const
 {
-    const polyMesh& mesh = cloud_.mesh();
+    const polyMesh& mesh = cloud_.pMesh();
     const labelList& faces = mesh.cells()[celli];
     const vector& C = mesh.cellCentres()[celli];
 
@@ -184,7 +178,7 @@ Foam::Particle<ParticleType>::Particle
     facei_(-1),
     stepFraction_(0.0),
     origProc_(Pstream::myProcNo()),
-    origId_(particleCount++)
+    origId_(cloud_.getNewParticleID())
 {}
 
 
@@ -313,13 +307,13 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
         // change cell
         if (internalFace) // Internal face
         {
-            if (celli_ == cloud_.owner_[facei_])
+            if (celli_ == mesh.faceOwner()[facei_])
             {
-                celli_ = cloud_.neighbour_[facei_];
+                celli_ = mesh.faceNeighbour()[facei_];
             }
-            else if (celli_ == cloud_.neighbour_[facei_])
+            else if (celli_ == mesh.faceNeighbour()[facei_])
             {
-                celli_ = cloud_.owner_[facei_];
+                celli_ = mesh.faceOwner()[facei_];
             }
             else
             {

src/lagrangian/basic/Particle/Particle.H

@@ -279,9 +279,6 @@ public:
         //- String representation of properties
         static string propHeader;
 
-        //- Cumulative particle count used for particle id
-        static label particleCount;
-
 
     // Constructors
 
@@ -330,7 +327,15 @@ public:
 
             autoPtr<ParticleType> operator()(Istream& is) const
             {
-                return autoPtr<ParticleType>(new ParticleType(cloud_, is));
+                return autoPtr<ParticleType>
+                (
+                    new ParticleType
+                    (
+                        cloud_,
+                        is,
+                        true
+                    )
+                );
             }
         };
 
@@ -454,9 +459,14 @@ public:
 
     // I-O
 
+        //- Read the fields associated with the owner cloud
+        static void readFields(Cloud<ParticleType>& c);
+
         //- Write the fields associated with the owner cloud
         static void writeFields(const Cloud<ParticleType>& c);
 
+        //- Write the particle data
+        void write(Ostream& os, bool writeFields) const;
 
     // Ostream Operator
 

src/lagrangian/basic/Particle/ParticleIO.C

@@ -50,13 +50,26 @@ Foam::Particle<ParticleType>::Particle
     origProc_(Pstream::myProcNo()),
     origId_(-1)
 {
+
+    // readFields : read additional data. Should be consistent with writeFields.
+
     if (is.format() == IOstream::ASCII)
     {
-        is >> position_ >> celli_ >> origProc_ >> origId_;
+        is >> position_ >> celli_;
+        if (readFields)
+        {
+            is >> origProc_ >> origId_;
+        }
+        else
+        {
+            origId_ = cloud_.getNewParticleID();
+        }
     }
     else
     {
         // In binary read all particle data - needed for parallel transfer
+        if (readFields)
+        {
             is.read
             (
                 reinterpret_cast<char*>(&position_),
@@ -68,6 +81,19 @@ Foam::Particle<ParticleType>::Particle
               + sizeof(origId_)
             );
         }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&position_),
+                sizeof(position_)
+              + sizeof(celli_)
+              + sizeof(facei_)
+              + sizeof(stepFraction_)
+            );
+            origId_ = cloud_.getNewParticleID();
+        }
+    }
 
     if (celli_ == -1)
     {
@@ -79,6 +105,39 @@ Foam::Particle<ParticleType>::Particle
 }
 
 
+template<class ParticleType>
+void Foam::Particle<ParticleType>::readFields
+(
+    Cloud<ParticleType>& c
+)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    IOobject procIO(c.fieldIOobject("origProcId", IOobject::MUST_READ));
+
+    if (procIO.headerOk())
+    {
+        IOField<label> origProcId(procIO);
+        c.checkFieldIOobject(c, origProcId);
+        IOField<label> origId(c.fieldIOobject("origId", IOobject::MUST_READ));
+        c.checkFieldIOobject(c, origId);
+
+        label i = 0;
+        forAllIter(typename Cloud<ParticleType>, c, iter)
+        {
+            ParticleType& p = iter();
+
+            p.origProc_ = origProcId[i];
+            p.origId_ = origId[i];
+            i++;
+        }
+    }
+}
+
+
 template<class ParticleType>
 void Foam::Particle<ParticleType>::writeFields
 (
@@ -88,36 +147,91 @@ void Foam::Particle<ParticleType>::writeFields
     // Write the cloud position file
     IOPosition<ParticleType> ioP(c);
     ioP.write();
+
+    label np =  c.size();
+
+    IOField<label> origProc
+    (
+        c.fieldIOobject
+        (
+            "origProcId",
+            IOobject::NO_READ
+        ),
+        np
+    );
+    IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np);
+
+    label i = 0;
+    forAllConstIter(typename Cloud<ParticleType>, c, iter)
+    {
+        origProc[i] = iter().origProc_;
+        origId[i] = iter().origId_;
+        i++;
+    }
+
+    origProc.write();
+    origId.write();
+}
+
+
+template<class ParticleType>
+void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        if (writeFields)
+        {
+            // Write the additional entries
+            os << position_
+               << token::SPACE << celli_
+               << token::SPACE << origProc_
+               << token::SPACE << origId_;
+        }
+        else
+        {
+            os << position_
+               << token::SPACE << celli_;
+        }
+    }
+    else
+    {
+        // In binary write both celli_ and facei_, needed for parallel transfer
+        if (writeFields)
+        {
+            os.write
+            (
+                reinterpret_cast<const char*>(&position_),
+                sizeof(position_)
+              + sizeof(celli_)
+              + sizeof(facei_)
+              + sizeof(stepFraction_)
+              + sizeof(origProc_)
+              + sizeof(origId_)
+            );
+        }
+        else
+        {
+            os.write
+            (
+                reinterpret_cast<const char*>(&position_),
+                sizeof(position_)
+              + sizeof(celli_)
+              + sizeof(facei_)
+              + sizeof(stepFraction_)
+            );
+        }
+    }
+
+    // Check state of Ostream
+    os.check("Particle<ParticleType>::write(Ostream& os, bool) const");
 }
 
 
 template<class ParticleType>
 Foam::Ostream& Foam::operator<<(Ostream& os, const Particle<ParticleType>& p)
 {
-    if (os.format() == IOstream::ASCII)
-    {
-        os << p.position_
-           << token::SPACE << p.celli_
-           << token::SPACE << p.origProc_
-           << token::SPACE << p.origId_;
-    }
-    else
-    {
-        // In binary write both celli_ and facei_, needed for parallel transfer
-        os.write
-        (
-            reinterpret_cast<const char*>(&p.position_),
-            sizeof(p.position_)
-          + sizeof(p.celli_)
-          + sizeof(p.facei_)
-          + sizeof(p.stepFraction_)
-          + sizeof(p.origProc_)
-          + sizeof(p.origId_)
-        );
-    }
-
-    // Check state of Ostream
-    os.check("Ostream& operator<<(Ostream&, const Particle<ParticleType>&)");
+    // Write all data
+    p.write(os, true);
 
     return os;
 }

src/lagrangian/basic/indexedParticle/indexedParticle.H

@@ -84,7 +84,7 @@ public:
             bool readFields = true
         )
         :
-            Particle<indexedParticle>(c, is)
+            Particle<indexedParticle>(c, is, readFields)
         {}
 
         //- Construct as a copy

src/lagrangian/basic/indexedParticle/indexedParticleCloud.C

@@ -22,12 +22,9 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#include "indexedParticle.H"
-#include "Cloud.H"
+#include "indexedParticleCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,4 +40,25 @@ defineTemplateTypeNameAndDebug(Cloud<indexedParticle>, 0);
 
 } // End namespace Foam
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::indexedParticleCloud::indexedParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName
+)
+:
+    Cloud<indexedParticle>(mesh, cloudName, false)
+{
+    indexedParticle::readFields(*this);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::indexedParticleCloud::writeFields() const
+{
+    indexedParticle::writeFields(*this);
+}
+
 // ************************************************************************* //

src/lagrangian/basic/indexedParticle/indexedParticleCloud.H

@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::indexedParticleCloud
+
+Description
+    A Cloud of particles carrying an additional index.
+
+SourceFiles
+    indexedParticleCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedParticleCloud_H
+#define indexedParticleCloud_H
+
+#include "Cloud.H"
+#include "indexedParticle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class indexedParticleCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class indexedParticleCloud
+:
+    public Cloud<indexedParticle>
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        indexedParticleCloud(const indexedParticleCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const indexedParticleCloud&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given mesh
+        indexedParticleCloud
+        (
+            const polyMesh&,
+            const word& cloudName = "defaultCloud"
+        );
+
+
+    // Member Functions
+
+        //- Write fields
+        virtual void writeFields() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/basic/passiveParticle/passiveParticle.H

@@ -28,9 +28,6 @@ Class
 Description
 
 SourceFiles
-    passiveParticleI.H
-    passiveParticle.C
-    passiveParticleIO.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -78,7 +75,7 @@ public:
             bool readFields = true
         )
         :
-            Particle<passiveParticle>(c, is)
+            Particle<passiveParticle>(c, is, readFields)
         {}
 
         //- Construct as copy

src/lagrangian/basic/passiveParticle/passiveParticleCloud.C

@@ -22,12 +22,9 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#include "passiveParticle.H"
-#include "Cloud.H"
+#include "passiveParticleCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,8 +36,34 @@ namespace Foam
 defineParticleTypeNameAndDebug(passiveParticle, 0);
 defineTemplateTypeNameAndDebug(Cloud<passiveParticle>, 0);
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+};
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::passiveParticleCloud::passiveParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName
+)
+:
+    Cloud<passiveParticle>(mesh, cloudName, false)
+{
+    readFields();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::passiveParticleCloud::readFields()
+{
+    passiveParticle::readFields(*this);
+}
+
+
+void Foam::passiveParticleCloud::writeFields() const
+{
+    passiveParticle::writeFields(*this);
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //

src/lagrangian/basic/passiveParticle/passiveParticleCloud.H

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::passiveParticleCloud
+
+Description
+    A Cloud of passive particles
+
+SourceFiles
+    passiveParticleCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef passiveParticleCloud_H
+#define passiveParticleCloud_H
+
+#include "Cloud.H"
+#include "passiveParticle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class passiveParticleCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class passiveParticleCloud
+:
+    public Cloud<passiveParticle>
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        passiveParticleCloud(const passiveParticleCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const passiveParticleCloud&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given mesh
+        passiveParticleCloud
+        (
+            const polyMesh&,
+            const word& cloudName = "defaultCloud"
+        );
+
+
+    // Member Functions
+
+        //- Read fields
+        virtual void readFields();
+
+        //- Write fields
+        virtual void writeFields() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/CoalCloud/CoalCloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/coalCombustion/CoalCloud/CoalCloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,7 +87,7 @@ public:
     // Member Functions
 
         //- Write fields
-        void writeFields() const;
+        virtual void writeFields() const;
 };
 
 

src/lagrangian/coalCombustion/CoalParcel/CoalParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/coalCombustion/CoalParcel/CoalParcel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/coalCombustion/CoalParcel/defineCoalParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/dieselSpray/parcel/parcelIO.C

@@ -36,7 +36,7 @@ Foam::parcel::parcel
     bool readFields
 )
 :
-    Particle<parcel>(cloud, is),
+    Particle<parcel>(cloud, is, readFields),
 
     liquidComponents_
     (

src/lagrangian/dieselSpray/spray/spray.H

@@ -345,7 +345,7 @@ public:
     // I/O
 
             //- Write fields
-            void writeFields() const;
+            virtual void writeFields() const;
 };
 
 

src/lagrangian/dieselSpray/spraySubModels/wallModel/reflectParcel/reflectParcel.C2

@@ -1,211 +0,0 @@
-// The FOAM Project // File: reflectParcel.C
-/*
--------------------------------------------------------------------------------
- =========         | Class Implementation
- \\      /         |
-  \\    /          | Name:   reflectParcel
-   \\  /           | Family: wallModel
-    \\/            |
-    F ield         | FOAM version: 2.3
-    O peration     |
-    A and          | Copyright (C) 1991-2009 OpenCFD Ltd.
-    M anipulation  |          All Rights Reserved.
--------------------------------------------------------------------------------
-DESCRIPTION
-
-AUTHOR
-    Henry G Weller.
-
--------------------------------------------------------------------------------
-*/
-
-#include "reflectParcel.H"
-#include "addToRunTimeSelectionTable.H"
-#include "wallPolyPatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(reflectParcel, 0);
-
-addToRunTimeSelectionTable
-(
-    wallModel,
-    reflectParcel,
-    dictionary
-);
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-reflectParcel::reflectParcel
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    spray& sm
-)
-:
-    wallModel(dict, mesh, sm),
-    coeffsDict_(dict.subDict(typeName + "Coeffs")),
-    elasticity_(readScalar(coeffsDict_.lookup("elasticity")))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-reflectParcel::~reflectParcel()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// Return 'keepParcel'
-bool reflectParcel::wallTreatment
-(
-    parcel& p
-) const
-{
-    label patchi = p.patch();
-    label facei = p.patchFace(patchi);
-
-    const polyMesh& mesh = spray_.mesh();
-
-    if (typeid(mesh_.boundaryMesh()[patchi]) == typeid(wallPolyPatch))
-    {
-        // wallNormal defined to point outwards of domain
-        vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
-        Sf /= mag(Sf);
-
-        // adjust velocity when wall is moving
-        if (!mesh.moving())
-        {
-            scalar Un = p.U() & Sf;
-            
-            if (Un > 0)
-            {
-                p.U() -= (1.0 + elasticity_)*Un*Sf;
-            }
-
-        }
-        else
-        {
-            label facep = p.face();
-
-            scalar dt = spray_.runTime().deltaT().value();
-            scalar fraction = p.t0()/dt;
-            const vectorField& oldPoints = mesh.oldPoints();
-            const vector& Cf1 = mesh.faceCentres()[facep];
-
-            vector Cf0 = mesh.faces()[facep].centre(oldPoints);
-            vector Cf = Cf0 + fraction*(Cf1 - Cf0);
-            vector Sf0 = mesh.faces()[facep].normal(oldPoints);
-            Sf0 /= mag(Sf0);
-
-            scalar magSfDiff = mag(Sf - Sf0);
-
-            if (magSfDiff > SMALL)
-            {
-                bool pureRotation = false;
-
-                if (pureRotation)
-                {
-                    // rotation
-                    
-                    // find center of rotation
-                    vector omega = Sf0 ^ Sf;
-                    scalar magOmega = mag(omega);
-                    omega /= magOmega+SMALL;
-                    vector n0 = omega ^ Sf0;
-                    scalar lam = ((Cf1 - Cf0) & Sf)/(n0 & Sf);
-                    vector r0 = Cf0 + lam*n0;
-                    
-                    scalar phiVel = ::asin(magOmega)/dt;
-                    
-                    vector pos = p.position() - r0;
-                    vector v = phiVel*(omega ^ pos);
-
-                    vector Sfp = Sf0 + fraction*(Sf - Sf0);
-                    //vector Sfp = Sf;
-
-                    vector Ur = p.U() - v;
-                    scalar Urn = Ur & Sfp;
-                    /*
-                    scalar dd = (p.position() - r0) & Sfp;
-                    Info << "Urn = " << Urn
-                        << ", dd = " << dd
-                        << ", pos = " << p.position()
-                        << ", Sfp = " << Sfp
-                        << ", omega = " << omega
-                        << ", r0 = " << r0
-                        << endl;
-                        */
-                    if (Urn > 0.0)
-                    {
-                        p.U() -= (1.0 + elasticity_)*Urn*Sfp;
-                    }
-                }
-                else
-                {
-                    vector Sfp = Sf0 + fraction*(Sf - Sf0);
-                    //vector Sfp = Sf;
-                    vector omega = Sf0 ^ Sf;
-                    scalar magOmega = mag(omega);
-                    omega /= magOmega+SMALL;
-                    
-                    scalar phiVel = ::asin(magOmega)/dt;
-                    
-                    scalar dist = (p.position() - Cf) & Sfp;
-                    vector pos = p.position() - dist*Sfp;
-                    vector vrot = phiVel*(omega ^ (pos - Cf));
-                    vector vtrans = (Cf1 - Cf0)/dt;
-                    vector v = vtrans + vrot;
-
-                    scalar Un = ((p.U() - v) & Sfp);
-                    
-                    if (Un > 0.0)
-                    {
-                        p.U() -= (1.0 + elasticity_)*Un*Sfp;
-                    }
-                }
-            }
-            else
-            {
-                // translation
-
-                vector v = (Cf1 - Cf0)/dt;
-                vector Ur = p.U() - v;
-                scalar Urn = Ur & Sf;
-
-                if (Urn > 0.0)
-                {
-                    Ur -= (1.0 + elasticity_)*Urn*Sf;
-                    p.U() = Ur + v;
-                }
-
-
-            }
-        }
-
-    }
-    else
-    {
-        FatalError
-            << "bool reflectParcel::wallTreatment(parcel& parcel) const "
-                << " parcel has hit a boundary " 
-                << mesh_.boundary()[patchi].type()
-                << " which not yet has been implemented."
-            << abort(FatalError);
-    }
-    return true;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C

@@ -299,7 +299,7 @@ void Foam::DsmcCloud<ParcelType>::collisions()
     // Temporary storage for subCells
     List<DynamicList<label> > subCells(8);
 
-    scalar deltaT = mesh_.time().deltaT().value();
+    scalar deltaT = cachedDeltaT();
 
     label collisionCandidates = 0;
 
@@ -778,6 +778,9 @@ Foam::DsmcCloud<ParcelType>::~DsmcCloud()
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::evolve()
 {
+    // cache the value of deltaT for this timestep
+    storeDeltaT();
+
     typename ParcelType::trackData td(*this);
 
     // Reset the surface data collection fields

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H

@@ -116,6 +116,10 @@ class DsmcCloud
         //- Random number generator
         Random rndGen_;
 
+        //- In-cloud cache of deltaT, lookup in submodels and parcel is
+        //  expensive
+        scalar cachedDeltaT_;
+
 
         // References to the macroscopic fields
 
@@ -243,6 +247,12 @@ public:
                 //- Return refernce to the random object
                 inline Random& rndGen();
 
+                //- Store (cache) the current value of deltaT
+                inline void storeDeltaT();
+
+                //- Return the cached value of deltaT
+                inline scalar cachedDeltaT() const;
+
 
             // References to the surface data collection fields
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H

@@ -120,6 +120,20 @@ inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
 }
 
 
+template<class ParcelType>
+inline void Foam::DsmcCloud<ParcelType>::storeDeltaT()
+{
+    cachedDeltaT_ = mesh().time().deltaT().value();
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::cachedDeltaT() const
+{
+    return cachedDeltaT_;
+}
+
+
 template<class ParcelType>
 inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
 {

src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H

@@ -91,7 +91,7 @@ public:
     // Member functions
 
         //- Write fields
-        void writeFields() const;
+        virtual void writeFields() const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C

@@ -44,7 +44,7 @@ bool Foam::DsmcParcel<ParcelType>::move
     const polyMesh& mesh = td.cloud().pMesh();
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
 
-    const scalar deltaT = mesh.time().deltaT().value();
+    const scalar deltaT = td.cloud().cachedDeltaT();
     scalar tEnd = (1.0 - p.stepFraction())*deltaT;
     const scalar dtMax = tEnd;
 
@@ -145,7 +145,7 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
 
     const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
 
-    const scalar deltaT = td.cloud().mesh().time().deltaT().value();
+    const scalar deltaT = td.cloud().cachedDeltaT();
 
     scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);
 

src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

@@ -126,7 +126,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
     const polyMesh& mesh(cloud.mesh());
 
-    const scalar deltaT = mesh.time().deltaT().value();
+    const scalar deltaT = cloud.cachedDeltaT();
 
     Random& rndGen(cloud.rndGen());
 

src/lagrangian/intermediate/Make/files

@@ -59,6 +59,7 @@ submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
 /* data entries */
 submodels/IO/DataEntry/makeDataEntries.C
 submodels/IO/DataEntry/polynomial/polynomial.C
+submodels/IO/DataEntry/polynomial/polynomialIO.C
 
 
 /* integration schemes */

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -173,9 +173,34 @@ void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
 }
 
 
+template<class ParcelType>
+void Foam::KinematicCloud<ParcelType>::preEvolve()
+{
+    this->dispersion().cacheFields(true);
+    forces_.cacheFields(true);
+}
+
+
+template<class ParcelType>
+void Foam::KinematicCloud<ParcelType>::postEvolve()
+{
+    if (debug)
+    {
+        this->writePositions();
+    }
+
+    this->dispersion().cacheFields(false);
+    forces_.cacheFields(false);
+
+    this->postProcessing().post();
+}
+
+
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::evolve()
 {
+    preEvolve();
+
     autoPtr<interpolation<scalar> > rhoInterpolator =
         interpolation<scalar>::New
         (
@@ -209,11 +234,6 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 
     this->injection().inject(td);
 
-    if (debug)
-    {
-        this->writePositions();
-    }
-
     if (coupled_)
     {
         resetSourceTerms();
@@ -221,7 +241,7 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 
     Cloud<ParcelType>::move(td);
 
-    this->postProcessing().post();
+    postEvolve();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -184,6 +184,15 @@ protected:
             DimensionedField<vector, volMesh> UTrans_;
 
 
+        // Cloud evolution functions
+
+            //- Pre-evolve
+            void preEvolve();
+
+            //- Post-evolve
+            void postEvolve();
+
+
 public:
 
     // Constructors
@@ -265,18 +274,6 @@ public:
                 inline const dictionary& interpolationSchemes() const;
 
 
-            // Forces to include in particle motion evaluation
-
-                //- Return reference to the gravity force flag
-                inline Switch forceGravity() const;
-
-                //- Return reference to the virtual mass force flag
-                inline Switch forceVirtualMass() const;
-
-                //- Return reference to the pressure gradient force flag
-                inline Switch forcePressureGradient() const;
-
-
             // Sub-models
 
                 //- Return const-access to the dispersion model

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -175,9 +175,25 @@ void Foam::ReactingCloud<ParcelType>::resetSourceTerms()
 }
 
 
+template<class ParcelType>
+void Foam::ReactingCloud<ParcelType>::preEvolve()
+{
+    ThermoCloud<ParcelType>::preEvolve();
+}
+
+
+template<class ParcelType>
+void Foam::ReactingCloud<ParcelType>::postEvolve()
+{
+    ThermoCloud<ParcelType>::postEvolve();
+}
+
+
 template<class ParcelType>
 void Foam::ReactingCloud<ParcelType>::evolve()
 {
+    preEvolve();
+
     const volScalarField& T = this->carrierThermo().T();
     const volScalarField cp = this->carrierThermo().Cp();
     const volScalarField& p = this->carrierThermo().p();
@@ -233,11 +249,6 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 
     this->injection().inject(td);
 
-    if (debug)
-    {
-        this->writePositions();
-    }
-
     if (this->coupled())
     {
         resetSourceTerms();
@@ -245,7 +256,7 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 
     Cloud<ParcelType>::move(td);
 
-    this->postProcessing().post();
+    postEvolve();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -132,6 +132,15 @@ protected:
             );
 
 
+        // Cloud evolution functions
+
+            //- Pre-evolve
+            void preEvolve();
+
+            //- Post-evolve
+            void postEvolve();
+
+
 public:
 
     // Constructors

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C

@@ -129,9 +129,25 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::resetSourceTerms()
 }
 
 
+template<class ParcelType>
+void Foam::ReactingMultiphaseCloud<ParcelType>::preEvolve()
+{
+    ReactingCloud<ParcelType>::preEvolve();
+}
+
+
+template<class ParcelType>
+void Foam::ReactingMultiphaseCloud<ParcelType>::postEvolve()
+{
+    ReactingCloud<ParcelType>::postEvolve();
+}
+
+
 template<class ParcelType>
 void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
 {
+    preEvolve();
+
     const volScalarField& T = this->carrierThermo().T();
     const volScalarField cp = this->carrierThermo().Cp();
     const volScalarField& p = this->carrierThermo().p();
@@ -187,11 +203,6 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
 
     this->injection().inject(td);
 
-    if (debug)
-    {
-        this->writePositions();
-    }
-
     if (this->coupled())
     {
         resetSourceTerms();
@@ -199,7 +210,7 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
 
     Cloud<ParcelType>::move(td);
 
-    this->postProcessing().post();
+    postEvolve();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H

@@ -68,7 +68,7 @@ class ReactingMultiphaseCloud
     public ReactingCloud<ParcelType>,
     public reactingMultiphaseCloud
 {
-    // Private Member Functions
+    // Private member functions
 
         //- Disallow default bitwise copy construct
         ReactingMultiphaseCloud(const ReactingMultiphaseCloud&);
@@ -112,6 +112,17 @@ protected:
             scalar dMassSurfaceReaction_;
 
 
+    // Protected member functions
+
+        // Cloud evolution functions
+
+            //- Pre-evolve
+            void preEvolve();
+
+            //- Post-evolve
+            void postEvolve();
+
+
 public:
 
     // Constructors

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -142,9 +142,25 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
 }
 
 
+template<class ParcelType>
+void Foam::ThermoCloud<ParcelType>::preEvolve()
+{
+    KinematicCloud<ParcelType>::preEvolve();
+}
+
+
+template<class ParcelType>
+void Foam::ThermoCloud<ParcelType>::postEvolve()
+{
+    KinematicCloud<ParcelType>::postEvolve();
+}
+
+
 template<class ParcelType>
 void Foam::ThermoCloud<ParcelType>::evolve()
 {
+    preEvolve();
+
     const volScalarField& T = carrierThermo_.T();
     const volScalarField cp = carrierThermo_.Cp();
 
@@ -192,11 +208,6 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 
     this->injection().inject(td);
 
-    if (debug)
-    {
-        this->writePositions();
-    }
-
     if (this->coupled())
     {
         resetSourceTerms();
@@ -204,7 +215,7 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 
     Cloud<ParcelType>::move(td);
 
-    this->postProcessing().post();
+    postEvolve();
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -118,6 +118,17 @@ protected:
             DimensionedField<scalar, volMesh> hcTrans_;
 
 
+    // Protected member functions
+
+        // Cloud evolution functions
+
+            //- Pre-evolve
+            void preEvolve();
+
+            //- Post-evolve
+            void postEvolve();
+
+
 public:
 
     // Constructors

src/lagrangian/intermediate/clouds/derived/BasicReactingCloud/BasicReactingCloud.H

@@ -92,7 +92,7 @@ public:
     // Member Functions
 
         //- Write fields
-        void writeFields() const;
+        virtual void writeFields() const;
 };
 
 

src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H

@@ -92,7 +92,7 @@ public:
     // Member Functions
 
         //- Write fields
-        void writeFields() const;
+        virtual void writeFields() const;
 };
 
 

src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.H

@@ -87,7 +87,7 @@ public:
     // Member functions
 
         //- Write fields
-        void writeFields() const;
+        virtual void writeFields() const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H

@@ -90,7 +90,7 @@ public:
     // Member Functions
 
         //- Write fields
-        void writeFields() const;
+        virtual void writeFields() const;
 };
 
 

src/lagrangian/intermediate/parcels/include/makeParcelInjectionModels.H

@@ -37,6 +37,7 @@ License
 #include "KinematicLookupTableInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"
+#include "PatchInjection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -79,6 +80,12 @@ License
         NoInjection,                                                          \
         KinematicCloud,                                                       \
         ParcelType                                                            \
+    );                                                                        \
+    makeInjectionModelType                                                    \
+    (                                                                         \
+        PatchInjection,                                                       \
+        KinematicCloud,                                                       \
+        ParcelType                                                            \
     );
 
 

src/lagrangian/intermediate/parcels/include/makeReactingMultiphaseParcelInjectionModels.H

@@ -37,6 +37,7 @@ License
 #include "FieldActivatedInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"
+#include "PatchInjection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -100,6 +101,13 @@ License
         ParcelType,                                                           \
         ThermoType                                                            \
     );                                                                        \
+    makeInjectionModelThermoType                                              \
+    (                                                                         \
+        PatchInjection,                                                       \
+        KinematicCloud,                                                       \
+        ParcelType,                                                           \
+        ThermoType                                                            \
+    );                                                                        \
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/parcels/include/makeReactingParcelInjectionModels.H

@@ -37,6 +37,7 @@ License
 #include "FieldActivatedInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"
+#include "PatchInjection.H"
 #include "ReactingLookupTableInjection.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -102,6 +103,13 @@ License
         ThermoType                                                            \
     );                                                                        \
     makeInjectionModelThermoType                                              \
+    (                                                                         \
+        PatchInjection,                                                       \
+        KinematicCloud,                                                       \
+        ParcelType,                                                           \
+        ThermoType                                                            \
+    );                                                                        \
+    makeInjectionModelThermoType                                              \
     (                                                                         \
         ReactingLookupTableInjection,                                         \
         KinematicCloud,                                                       \

src/lagrangian/intermediate/particleForces/particleForces.C

@@ -27,6 +27,7 @@ License
 #include "particleForces.H"
 #include "fvMesh.H"
 #include "volFields.H"
+#include "fvcGrad.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -40,21 +41,17 @@ Foam::particleForces::particleForces
     mesh_(mesh),
     dict_(dict.subDict("particleForces")),
     g_(g),
+    gradUPtr_(NULL),
     gravity_(dict_.lookup("gravity")),
     virtualMass_(dict_.lookup("virtualMass")),
     Cvm_(0.0),
     pressureGradient_(dict_.lookup("pressureGradient")),
-    gradUName_("unknown_gradUName")
+    UName_(dict_.lookupOrDefault<word>("U", "U"))
 {
-    if (gravity_)
+    if (virtualMass_)
     {
         dict_.lookup("Cvm") >> Cvm_;
     }
-
-    if (pressureGradient_)
-    {
-        dict_.lookup("gradU") >> gradUName_;
-    }
 }
 
 
@@ -63,18 +60,21 @@ Foam::particleForces::particleForces(const particleForces& f)
     mesh_(f.mesh_),
     dict_(f.dict_),
     g_(f.g_),
+    gradUPtr_(f.gradUPtr_),
     gravity_(f.gravity_),
     virtualMass_(f.virtualMass_),
     Cvm_(f.Cvm_),
     pressureGradient_(f.pressureGradient_),
-    gradUName_(f.gradUName_)
+    UName_(f.UName_)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::particleForces::~particleForces()
-{}
+{
+    cacheFields(false);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -109,9 +109,26 @@ Foam::Switch Foam::particleForces::pressureGradient() const
 }
 
 
-const Foam::word& Foam::particleForces::gradUName() const
+const Foam::word& Foam::particleForces::UName() const
 {
-    return gradUName_;
+    return UName_;
+}
+
+
+void Foam::particleForces::cacheFields(const bool store)
+{
+    if (store && pressureGradient_)
+    {
+        const volVectorField U = mesh_.lookupObject<volVectorField>(UName_);
+        gradUPtr_ = fvc::grad(U).ptr();
+    }
+    else
+    {
+        if (gradUPtr_)
+        {
+            delete gradUPtr_;
+        }
+    }
 }
 
 
@@ -130,15 +147,17 @@ Foam::vector Foam::particleForces::calcCoupled
     // Virtual mass force
     if (virtualMass_)
     {
-        notImplemented("Foam::particleForces::calc(...) - virtualMass force");
+        notImplemented
+        (
+            "Foam::particleForces::calcCoupled(...) - virtual mass force"
+        );
 //        Ftot += Cvm_*rhoc/rho*d(Uc - U)/dt;
     }
 
     // Pressure gradient force
     if (pressureGradient_)
     {
-        const volSymmTensorField& gradU =
-            mesh_.lookupObject<volSymmTensorField>(gradUName_);
+        const volTensorField& gradU = *gradUPtr_;
         Ftot += rhoc/rho*(U & gradU[cellI]);
     }
 

src/lagrangian/intermediate/particleForces/particleForces.H

@@ -39,6 +39,7 @@ SourceFiles
 #include "dictionary.H"
 #include "Switch.H"
 #include "vector.H"
+#include "volFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,6 +66,9 @@ class particleForces
         //- Gravity
         const vector g_;
 
+        //- Velocity gradient field
+        const volTensorField* gradUPtr_;
+
 
         // Forces to include in particle motion evaluation
 
@@ -80,8 +84,11 @@ class particleForces
             //- Pressure gradient
             Switch pressureGradient_;
 
-            //- Name of velocity gradient field for pressure gradient force
-            word gradUName_;
+
+        // Additional info
+
+            //- Name of velucity field - default = "U"
+            const word UName_;
 
 
 public:
@@ -126,12 +133,15 @@ public:
             //- Return pressure gradient force activate switch
             Switch pressureGradient() const;
 
-            //- Return the name of the velocity gradient field
-            const word& gradUName() const;
+            //- Return name of velocity field
+            const word& UName() const;
 
 
        // Evaluation
 
+            //- Cache carrier fields
+            void cacheFields(const bool store);
+
             //- Calculate action/reaction forces between carrier and particles
             vector calcCoupled
             (

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionModel/DispersionModel.H

@@ -121,6 +121,9 @@ public:
         //- Flag to indicate whether model activates injection model
         virtual bool active() const = 0;
 
+        //- Cache carrier fields
+        virtual void cacheFields(const bool store) = 0;
+
         //- Update (disperse particles)
         virtual vector update
         (

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionRASModel/DispersionRASModel.C

@@ -42,7 +42,11 @@ Foam::DispersionRASModel<CloudType>::DispersionRASModel
         (
             "RASProperties"
         )
-    )
+    ),
+    kPtr_(NULL),
+    ownK_(false),
+    epsilonPtr_(NULL),
+    ownEpsilon_(false)
 {}
 
 
@@ -50,7 +54,56 @@ Foam::DispersionRASModel<CloudType>::DispersionRASModel
 
 template<class CloudType>
 Foam::DispersionRASModel<CloudType>::~DispersionRASModel()
-{}
+{
+    cacheFields(false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::DispersionRASModel<CloudType>::cacheFields(const bool store)
+{
+    if (store)
+    {
+        tmp<volScalarField> tk = this->turbulence().k();
+        if (tk.isTmp())
+        {
+            kPtr_ = tk.ptr();
+            ownK_ = true;
+        }
+        else
+        {
+            kPtr_ = tk.operator->();
+            ownK_ = false;
+        }
+
+        tmp<volScalarField> tepsilon = this->turbulence().epsilon();
+        if (tepsilon.isTmp())
+        {
+            epsilonPtr_ = tepsilon.ptr();
+            ownEpsilon_ = true;
+        }
+        else
+        {
+            epsilonPtr_ = tepsilon.operator->();
+            ownEpsilon_ = false;
+        }
+    }
+    else
+    {
+        if (ownK_ && kPtr_)
+        {
+            delete kPtr_;
+            ownK_ = false;
+        }
+        if (ownEpsilon_ && epsilonPtr_)
+        {
+            delete epsilonPtr_;
+            ownEpsilon_ = false;
+        }
+    }
+}
 
 
 // ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionRASModel/DispersionRASModel.H

@@ -50,11 +50,27 @@ class DispersionRASModel
 :
     public DispersionModel<CloudType>
 {
-    // Private data
+protected:
+
+    // Protected data
 
         //- Reference to the compressible turbulence model
         const compressible::RASModel& turbulence_;
 
+        // Locally cached turbulence fields
+
+            //- Turbulence k
+            const volScalarField* kPtr_;
+
+            //- Take ownership of the k field
+            bool ownK_;
+
+            //- Turbulence epsilon
+            const volScalarField* epsilonPtr_;
+
+            //- Take ownership of the epsilon field
+            bool ownEpsilon_;
+
 
 public:
 
@@ -78,6 +94,9 @@ public:
 
     // Member Functions
 
+        //- Cache carrier fields
+        virtual void cacheFields(const bool store);
+
         //- Return const access to the turbulence model
         const compressible::RASModel& turbulence() const
         {

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C

@@ -35,7 +35,8 @@ Foam::GradientDispersionRAS<CloudType>::GradientDispersionRAS
     CloudType& owner
 )
 :
-    DispersionRASModel<CloudType>(dict, owner)
+    DispersionRASModel<CloudType>(dict, owner),
+    gradkPtr_(NULL)
 {}
 
 
@@ -43,7 +44,9 @@ Foam::GradientDispersionRAS<CloudType>::GradientDispersionRAS
 
 template<class CloudType>
 Foam::GradientDispersionRAS<CloudType>::~GradientDispersionRAS()
-{}
+{
+    cacheFields(false);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -55,6 +58,25 @@ bool Foam::GradientDispersionRAS<CloudType>::active() const
 }
 
 
+template<class CloudType>
+void Foam::GradientDispersionRAS<CloudType>::cacheFields(const bool store)
+{
+    DispersionRASModel<CloudType>::cacheFields(store);
+
+    if (store)
+    {
+        gradkPtr_ = fvc::grad(*this->kPtr_).ptr();
+    }
+    else
+    {
+        if (gradkPtr_)
+        {
+            delete gradkPtr_;
+        }
+    }
+}
+
+
 template<class CloudType>
 Foam::vector Foam::GradientDispersionRAS<CloudType>::update
 (
@@ -68,9 +90,9 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
 {
     const scalar cps = 0.16432;
 
-    const volScalarField& k = this->turbulence().k();
-    const volScalarField& epsilon = this->turbulence().epsilon();
-    const volVectorField gradk = fvc::grad(k);
+    const volScalarField& k = *this->kPtr_;
+    const volScalarField& epsilon = *this->epsilonPtr_;
+    const volVectorField& gradk = *this->gradkPtr_;
 
     const scalar UrelMag = mag(U - Uc - UTurb);
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.H

@@ -51,6 +51,14 @@ class GradientDispersionRAS
 :
     public DispersionRASModel<CloudType>
 {
+protected:
+
+    // Locally cached turbulence fields
+
+        //- Gradient of k
+        const volVectorField* gradkPtr_;
+
+
 public:
 
     //- Runtime type information
@@ -76,8 +84,11 @@ public:
         //- Flag to indicate whether model activates injection model
         bool active() const;
 
+        //- Cache carrier fields
+        virtual void cacheFields(const bool store);
+
         //- Update (disperse particles)
-        vector update
+        virtual vector update
         (
             const scalar dt,
             const label celli,

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/NoDispersion/NoDispersion.C

@@ -55,6 +55,13 @@ bool Foam::NoDispersion<CloudType>::active() const
 }
 
 
+template<class CloudType>
+void Foam::NoDispersion<CloudType>::cacheFields(const bool)
+{
+// do nothing
+}
+
+
 template<class CloudType>
 Foam::vector Foam::NoDispersion<CloudType>::update
 (

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/NoDispersion/NoDispersion.H

@@ -72,10 +72,13 @@ public:
     // Member Functions
 
         //- Flag to indicate whether model activates injection model
-        bool active() const;
+        virtual bool active() const;
+
+        //- Cache carrier fields
+        virtual void cacheFields(const bool store);
 
         //- Update (disperse particles)
-        vector update
+        virtual vector update
         (
             const scalar dt,
             const label celli,

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C

@@ -68,8 +68,8 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
 {
     const scalar cps = 0.16432;
 
-    const volScalarField& k = this->turbulence().k();
-    const volScalarField& epsilon = this->turbulence().epsilon();
+    const volScalarField& k = *this->kPtr_;
+    const volScalarField& epsilon = *this->epsilonPtr_;
 
     const scalar UrelMag = mag(U - Uc - UTurb);
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.H

@@ -74,10 +74,10 @@ public:
     // Member Functions
 
         //- Flag to indicate whether model activates injection model
-        bool active() const;
+        virtual bool active() const;
 
         //- Update (disperse particles)
-        vector update
+        virtual vector update
         (
             const scalar dt,
             const label celli,

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -233,7 +233,8 @@ void Foam::InjectionModel<CloudType>::postInjectCheck(const label parcelsAdded)
 {
     if (parcelsAdded > 0)
     {
-        Pout<< "\n--> Cloud: " << owner_.name() << nl
+        Pout<< nl
+            << "--> Cloud: " << owner_.name() << nl
             << "    Added " << parcelsAdded
             << " new parcels" << nl << endl;
     }

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C

@@ -138,7 +138,7 @@ void Foam::ManualInjection<CloudType>::setPositionAndCell
 (
     const label parcelI,
     const label,
-    const scalar time,
+    const scalar,
     vector& position,
     label& cellOwner
 )

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchInjection.H"
+#include "DataEntry.H"
+#include "pdf.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+) const
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return round(fraction_*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::PatchInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+) const
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return fraction_*volumeFlowRate_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PatchInjection<CloudType>::PatchInjection
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    InjectionModel<CloudType>(dict, owner, typeName),
+    patchName_(this->coeffDict().lookup("patchName")),
+    duration_(readScalar(this->coeffDict().lookup("duration"))),
+    parcelsPerSecond_
+    (
+        readScalar(this->coeffDict().lookup("parcelsPerSecond"))
+    ),
+    U0_(this->coeffDict().lookup("U0")),
+    volumeFlowRate_
+    (
+        DataEntry<scalar>::New
+        (
+            "volumeFlowRate",
+            this->coeffDict()
+        )
+    ),
+    parcelPDF_
+    (
+        pdf::New
+        (
+            this->coeffDict().subDict("parcelPDF"),
+            owner.rndGen()
+        )
+    ),
+    cellOwners_(),
+    fraction_(1.0)
+{
+    label patchId = owner.mesh().boundaryMesh().findPatchID(patchName_);
+
+    if (patchId < 0)
+    {
+        FatalErrorIn
+        (
+            "PatchInjection<CloudType>::PatchInjection"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Requested patch " << patchName_ << " not found" << nl
+            << "Available patches are: " << owner.mesh().boundaryMesh().names()
+            << nl << exit(FatalError);
+    }
+
+    const polyPatch& patch = owner.mesh().boundaryMesh()[patchId];
+
+    cellOwners_ = patch.faceCells();
+
+    label patchSize = cellOwners_.size();
+    label totalPatchSize = patchSize;
+    reduce(totalPatchSize, sumOp<scalar>());
+    fraction_ = patchSize/totalPatchSize;
+
+    // Set total volume to inject
+    this->volumeTotal_ = fraction_*volumeFlowRate_().integrate(0.0, duration_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PatchInjection<CloudType>::~PatchInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::PatchInjection<CloudType>::active() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::PatchInjection<CloudType>::timeEnd() const
+{
+    return this->SOI_ + duration_;
+}
+
+
+template<class CloudType>
+void Foam::PatchInjection<CloudType>::setPositionAndCell
+(
+    const label,
+    const label,
+    const scalar,
+    vector& position,
+    label& cellOwner
+)
+{
+    label cellI = this->owner().rndGen().integer(0, cellOwners_.size() - 1);
+
+    cellOwner = cellOwners_[cellI];
+    position = this->owner().mesh().C()[cellOwner];
+}
+
+
+template<class CloudType>
+void Foam::PatchInjection<CloudType>::setProperties
+(
+    const label,
+    const label,
+    const scalar,
+    typename CloudType::parcelType& parcel
+)
+{
+    // set particle velocity
+    parcel.U() = U0_;
+
+    // set particle diameter
+    parcel.d() = parcelPDF_->sample();
+}
+
+
+template<class CloudType>
+bool Foam::PatchInjection<CloudType>::fullyDescribed() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::PatchInjection<CloudType>::validInjection(const label)
+{
+    return true;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::PatchInjection
+
+Description
+    Patch injection
+
+    - User specifies
+      - Total mass to inject
+      - Name of patch
+      - Injection duration
+      - Initial parcel velocity
+      - Injection volume flow rate
+    - Parcel diameters obtained by PDF model
+    - Parcels injected at cell centres adjacent to patch
+
+SourceFiles
+    PatchInjection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PatchInjection_H
+#define PatchInjection_H
+
+#include "InjectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type>
+class DataEntry;
+
+class pdf;
+
+/*---------------------------------------------------------------------------*\
+                       Class PatchInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class PatchInjection
+:
+    public InjectionModel<CloudType>
+{
+    // Private data
+
+        //- Name of patch
+        const word patchName_;
+
+        //- Injection duration [s]
+        const scalar duration_;
+
+        //- Number of parcels to introduce per second []
+        const label parcelsPerSecond_;
+
+        //- Initial parcel velocity [m/s]
+        const vector U0_;
+
+        //- Volume flow rate of parcels to introduce relative to SOI [m^3/s]
+        const autoPtr<DataEntry<scalar> > volumeFlowRate_;
+
+        //- Parcel size PDF model
+        const autoPtr<pdf> parcelPDF_;
+
+        //- Cell owners
+        labelList cellOwners_;
+
+        //- Fraction of injection controlled by this processor
+        scalar fraction_;
+
+
+protected:
+
+    // Protected member functions
+
+        //- Number of parcels to introduce over the time step relative to SOI
+        label parcelsToInject
+        (
+            const scalar time0,
+            const scalar time1
+        ) const;
+
+        //- Volume of parcels to introduce over the time step relative to SOI
+        scalar volumeToInject
+        (
+            const scalar time0,
+            const scalar time1
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PatchInjection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        PatchInjection
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+
+    //- Destructor
+    virtual ~PatchInjection();
+
+
+    // Member Functions
+
+        //- Flag to indicate whether model activates injection model
+        bool active() const;
+
+        //- Return the end-of-injection time
+        scalar timeEnd() const;
+
+
+        // Injection geometry
+
+            //- Set the injection position and owner cell
+            virtual void setPositionAndCell
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                vector& position,
+                label& cellOwner
+            );
+
+            virtual void setProperties
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                typename CloudType::parcelType& parcel
+            );
+
+            //- Flag to identify whether model fully describes the parcel
+            virtual bool fullyDescribed() const;
+
+            //- Return flag to identify whether or not injection of parcelI is
+            //  permitted
+            virtual bool validInjection(const label parcelI);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PatchInjection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //