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

applications/solvers/combustion/fireFoam/pEqn.H

@@ -5,8 +5,8 @@ surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 phi.boundaryField() =
-    fvc::interpolate(rho.boundaryField())
+    fvc::interpolate(rho.boundaryField()*U.boundaryField())
-   *(fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField());
+  & mesh.Sf().boundaryField();
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 

applications/solvers/multiphase/bubbleFoam/pEqn.H

@@ -14,12 +14,6 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
-    U.boundaryField() =
-        alpha1.boundaryField()*U1.boundaryField()
-      + alpha2.boundaryField()*U2.boundaryField();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
-
     surfaceScalarField phiDrag1
     (
         fvc::interpolate(alpha2/rho1*dragCoef*rAU1)*phi2
@@ -44,16 +38,18 @@
     (
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + phiDrag1
     );
 
     surfaceScalarField phiHbyA2
     (
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + phiDrag2
     );
 
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+
+    phiHbyA1 += phiDrag1;
+    phiHbyA2 += phiDrag2;
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
     surfaceScalarField Dp

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -29,8 +29,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -38,6 +36,7 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -2,11 +2,6 @@
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
-    mrfZones.absoluteFlux(phi1.oldTime());
-    mrfZones.absoluteFlux(phi1);
-    mrfZones.absoluteFlux(phi2.oldTime());
-    mrfZones.absoluteFlux(phi2);
-
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
 
@@ -22,11 +17,10 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
-    U.boundaryField() =
+    mrfZones.absoluteFlux(phi1.oldTime());
-        alpha1.boundaryField()*U1.boundaryField()
+    mrfZones.absoluteFlux(phi1);
-      + alpha2.boundaryField()*U2.boundaryField();
+    mrfZones.absoluteFlux(phi2.oldTime());
-    phi.boundaryField() =
+    mrfZones.absoluteFlux(phi2);
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
@@ -47,18 +41,29 @@
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
-      + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
-      + ppDrag
-      + rAlphaAU1f*(g & mesh.Sf())
     );
-    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
         "phiHbyA2",
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
-      + fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
+    );
+
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+    mrfZones.relativeFlux(phi);
+
+    phiHbyA1 +=
+    (
+        fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + ppDrag
+      + rAlphaAU1f*(g & mesh.Sf())
+    );
+    mrfZones.relativeFlux(phiHbyA1);
+
+    phiHbyA2 +=
+    (
+        fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
       + rAlphaAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);

applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H

@@ -5,20 +5,15 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    mrfZones.absoluteFlux(phi);
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
-    mrfZones.relativeFlux(phi);
-
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiHbyA);
-
     adjustPhi(phiHbyA, U, p_rgh);
+    mrfZones.relativeFlux(phiHbyA);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -21,6 +19,8 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
+    phiAbs = phiHbyA;
+
     surfaceScalarField phig
     (
         (

applications/solvers/multiphase/interFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -15,6 +13,7 @@
     );
 
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -13,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -58,6 +58,8 @@
         dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
     );
 
+    phi = dimensionedScalar("phi", phi.dimensions(), 0);
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -74,14 +76,17 @@
         (
             (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
-          + rAlphaAUfs[phasei]
+        );
+        mrfZones.relativeFlux(phiHbyAs[phasei]);
+
+        phi += alphafs[phasei]*phiHbyAs[phasei];
+
+        phiHbyAs[phasei] +=
+            rAlphaAUfs[phasei]
            *(
                fluid.surfaceTension(phase)*mesh.magSf()/phase.rho()
              + (g & mesh.Sf())
-            )
+            );
-        );
-
-        mrfZones.relativeFlux(phiHbyAs[phasei]);
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
             fluid.dragModels().begin();

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H

@@ -5,20 +5,15 @@
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    mrfZones.absoluteFlux(phi);
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
-    mrfZones.relativeFlux(phi);
-
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiHbyA);
-
     adjustPhi(phiHbyA, U, p_rgh);
+    mrfZones.relativeFlux(phiHbyA);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -13,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -5,9 +5,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(rho.boundaryField()*U.boundaryField())
-      & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -18,6 +15,7 @@
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -4,8 +4,6 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (
@@ -13,8 +11,8 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     adjustPhi(phiHbyA, U, p_rgh);
+    phi = phiHbyA;
 
     surfaceScalarField phig
     (

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -18,13 +18,6 @@
     mrfZones.absoluteFlux(phi1);
     mrfZones.absoluteFlux(phi2.oldTime());
     mrfZones.absoluteFlux(phi2);
-    mrfZones.absoluteFlux(phi);
-
-    U.boundaryField() =
-        alpha1.boundaryField()*U1.boundaryField()
-      + alpha2.boundaryField()*U2.boundaryField();
-    phi.boundaryField() =
-        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
 
@@ -43,18 +36,29 @@
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, U1, phi1)
-      + fvc::interpolate(alpha2/rho1*K*rAU1)*phi2
-      + ppDrag
-      + rAU1f*(g & mesh.Sf())
     );
-    mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
     (
         "phiHbyA2",
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, U2, phi2)
-      + fvc::interpolate(alpha1/rho2*K*rAU2)*phi1
+    );
+
+    phi = alpha1f*phiHbyA1 + alpha2f*phiHbyA2;
+    mrfZones.relativeFlux(phi);
+
+    phiHbyA1 +=
+    (
+        fvc::interpolate(alpha2/rho1*K*rAU1)*phi2
+      + ppDrag
+      + rAU1f*(g & mesh.Sf())
+    );
+    mrfZones.relativeFlux(phiHbyA1);
+
+    phiHbyA2 +=
+    (
+        fvc::interpolate(alpha1/rho2*K*rAU2)*phi1
       + rAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);
@@ -63,7 +67,6 @@
     mrfZones.relativeFlux(phi1);
     mrfZones.relativeFlux(phi2.oldTime());
     mrfZones.relativeFlux(phi2);
-    mrfZones.relativeFlux(phi);
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 

applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H

@@ -161,7 +161,7 @@
             )
         );
     }
-    else if(nuType == "field")
+    else if (nuType == "field")
     {
         nuIO.readOpt() = IOobject::MUST_READ;
         nuPtr.reset

applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H

@@ -81,7 +81,7 @@ if (thermalStress)
         );
 
     }
-    else if(CType == "field")
+    else if (CType == "field")
     {
         CIO.readOpt() = IOobject::MUST_READ;
 
@@ -214,7 +214,7 @@ if (thermalStress)
     {
         FatalErrorIn
         (
-             "readThermalProperties.H"
+            "readThermalProperties.H"
         )   << "Valid type entries are uniform or field for alpha"
             << abort(FatalError);
     }
@@ -228,5 +228,4 @@ if (thermalStress)
 
     threeKalpha = threeK*alpha;
     DT = k/C;
-
 }

src/OpenFOAM/meshes/polyMesh/polyMeshClear.C

@@ -97,6 +97,16 @@ void Foam::polyMesh::clearAddressing()
     geometricD_ = Vector<label>::zero;
     solutionD_ = Vector<label>::zero;
 
+    // Update zones
+    pointZones_.clearAddressing();
+    faceZones_.clearAddressing();
+    cellZones_.clearAddressing();
+
+    // Remove the stored tet base points
+    tetBasePtIsPtr_.clear();
+    // Remove the cell tree
+    cellTreePtr_.clear();
+
     pointMesh::Delete(*this);
 }
 

src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryFwd.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DataEntryFws_H
+#define DataEntryFws_H
+
+#include "DataEntry.H"
+#include "vector.H"
+#include "symmTensor.H"
+#include "sphericalTensor.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef DataEntry<label> labelDataEntry;
+    typedef DataEntry<scalar> scalarDataEntry;
+    typedef DataEntry<vector> vectorDataEntry;
+    typedef DataEntry<symmTensor> symmTensorDataEntry;
+    typedef DataEntry<sphericalTensor> sphericalTensorDataEntry;
+    typedef DataEntry<tensor> tensorDataEntry;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.C

@@ -25,14 +25,14 @@ License
 
 #include "polynomial.H"
 #include "Time.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(polynomial, 0);
-    DataEntry<scalar>::adddictionaryConstructorToTable<polynomial>
+    addToRunTimeSelectionTable(scalarDataEntry, polynomial, dictionary);
-        addpolynomialConstructorToTable_;
 }
 
 
@@ -40,7 +40,7 @@ namespace Foam
 
 Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
 :
-    DataEntry<scalar>(entryName),
+    scalarDataEntry(entryName),
     coeffs_(),
     canIntegrate_(true),
     dimensions_(dimless)
@@ -52,15 +52,17 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
     is.putBack(firstToken);
     if (firstToken == token::BEGIN_SQR)
     {
-        is >> this->dimensions_;
+        is  >> this->dimensions_;
     }
 
     is  >> coeffs_;
 
     if (!coeffs_.size())
     {
-        FatalErrorIn("Foam::polynomial::polynomial(const word&, dictionary&)")
+        FatalErrorIn
-            << "polynomial coefficients for entry " << this->name_
+        (
+            "Foam::polynomial::polynomial(const word&, const dictionary&)"
+        )   << "polynomial coefficients for entry " << this->name_
             << " are invalid (empty)" << nl << exit(FatalError);
     }
 
@@ -77,8 +79,10 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
     {
         if (!canIntegrate_)
         {
-            WarningIn("Foam::polynomial::polynomial(const word&, dictionary&)")
+            WarningIn
-                << "Polynomial " << this->name_ << " cannot be integrated"
+            (
+                "Foam::polynomial::polynomial(const word&, const dictionary&)"
+            )   << "Polynomial " << this->name_ << " cannot be integrated"
                 << endl;
         }
     }
@@ -91,7 +95,7 @@ Foam::polynomial::polynomial
     const List<Tuple2<scalar, scalar> >& coeffs
 )
 :
-    DataEntry<scalar>(entryName),
+    scalarDataEntry(entryName),
     coeffs_(coeffs),
     canIntegrate_(true),
     dimensions_(dimless)
@@ -101,7 +105,7 @@ Foam::polynomial::polynomial
         FatalErrorIn
         (
             "Foam::polynomial::polynomial"
-            "(const word&, const List<Tuple2<scalar, scalar> >&&)"
+            "(const word&, const List<Tuple2<scalar, scalar> >&)"
         )   << "polynomial coefficients for entry " << this->name_
             << " are invalid (empty)" << nl << exit(FatalError);
     }
@@ -122,7 +126,7 @@ Foam::polynomial::polynomial
             WarningIn
             (
                 "Foam::polynomial::polynomial"
-                "(const word&, const List<Tuple2<scalar, scalar> >&&)"
+                "(const word&, const List<Tuple2<scalar, scalar> >&)"
             )   << "Polynomial " << this->name_ << " cannot be integrated"
                 << endl;
         }
@@ -132,7 +136,7 @@ Foam::polynomial::polynomial
 
 Foam::polynomial::polynomial(const polynomial& poly)
 :
-    DataEntry<scalar>(poly),
+    scalarDataEntry(poly),
     coeffs_(poly.coeffs_),
     canIntegrate_(poly.canIntegrate_),
     dimensions_(poly.dimensions_)
@@ -201,7 +205,8 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimValue
 
 Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
 (
-    const scalar x1, const scalar x2
+    const scalar x1,
+    const scalar x2
 ) const
 {
     return dimensioned<scalar>
@@ -212,4 +217,5 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
     );
 }
 
+
 // ************************************************************************* //

src/OpenFOAM/primitives/functions/DataEntry/polynomial/polynomial.H

@@ -48,6 +48,7 @@ SourceFiles
 #include "DataEntry.H"
 #include "Tuple2.H"
 #include "dimensionSet.H"
+#include "DataEntryFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,7 +71,7 @@ Ostream& operator<<
 
 class polynomial
 :
-    public DataEntry<scalar>
+    public scalarDataEntry
 {
     // Private data
 
@@ -107,9 +108,9 @@ public:
         polynomial(const polynomial& poly);
 
         //- Construct and return a clone
-        virtual tmp<DataEntry<scalar> > clone() const
+        virtual tmp<scalarDataEntry> clone() const
         {
-            return tmp<DataEntry<scalar> >(new polynomial(*this));
+            return tmp<scalarDataEntry>(new polynomial(*this));
         }
 
 

src/combustionModels/FSD/FSD.C

@@ -272,7 +272,7 @@ void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
 
     forAll(ft_, cellI)
     {
-        if(ft_[cellI] < ftStoich)
+        if (ft_[cellI] < ftStoich)
         {
             pc[cellI] = ft_[cellI]*(YprodTotal/ftStoich);
         }

src/combustionModels/Make/files

@@ -10,6 +10,7 @@ rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
 rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
 rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 
+diffusion/diffusions.C
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 
 PaSR/PaSRs.C

src/combustionModels/PaSR/PaSR.C

@@ -138,7 +138,7 @@ void Foam::combustionModels::PaSR<Type>::correct()
 
 template<class Type>
 Foam::tmp<Foam::fvScalarMatrix>
-Foam::combustionModels::PaSR<Type>::R(const volScalarField& Y) const
+Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
 {
     tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
 

src/combustionModels/PaSR/PaSR.H

@@ -103,7 +103,7 @@ public:
             virtual void correct();
 
             //- Fuel consumption rate matrix.
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;

src/combustionModels/combustionModel/combustionModel.H

@@ -135,7 +135,7 @@ public:
         virtual void correct() = 0;
 
         //- Fuel consumption rate matrix, i.e. source term for fuel equation
-        virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const = 0;
+        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const = 0;
 
         //- Heat release rate calculated from fuel consumption rate matrix
         virtual tmp<volScalarField> dQ() const = 0;

src/combustionModels/diffusion/diffusion.C

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "diffusion.H"
+#include "fvcGrad.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::diffusion
+(
+    const word& modelType, const fvMesh& mesh
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
+    C_(readScalar(this->coeffs().lookup("C"))),
+    oxidantName_(this->coeffs().template lookupOrDefault<word>("oxidant", "O2"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+diffusion<CombThermoType, ThermoType>::~diffusion()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void diffusion<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        this->singleMixturePtr_->fresCorrect();
+
+        const label fuelI = this->singleMixturePtr_->fuelIndex();
+
+        const volScalarField& YFuel =
+            this->thermoPtr_->composition().Y()[fuelI];
+
+        if (this->thermoPtr_->composition().contains(oxidantName_))
+        {
+            const volScalarField& YO2 =
+                this->thermoPtr_->composition().Y(oxidantName_);
+
+            this->wFuel_ ==
+                C_*this->turbulence().muEff()
+               *mag(fvc::grad(YFuel) & fvc::grad(YO2))
+               *pos(YFuel)*pos(YO2);
+        }
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool diffusion<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("C") >> C_ ;
+        this->coeffs().readIfPresent("oxidant", oxidantName_);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/diffusion/diffusion.H

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::diffusion
+
+Description
+    Simple diffusion-based combustion model based on the principle mixed is
+    burnt. Additional parameter C is used to distribute the heat release rate
+    in time.
+
+SourceFiles
+    diffusion.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef diffusion_H
+#define diffusion_H
+
+#include "singleStepCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class diffusion Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class diffusion
+:
+    public singleStepCombustion<CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Model constant
+        scalar C_;
+
+        //- Name of oxidant - default is "O2"
+        word oxidantName_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        diffusion(const diffusion&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const diffusion&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diffusion");
+
+
+    // Constructors
+
+        //- Construct from components
+        diffusion(const word& modelType, const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~diffusion();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+
+        // I-O
+
+            //- Update properties
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "diffusion.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/diffusion/diffusions.C

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiThermoCombustion.H"
+#include "rhoThermoCombustion.H"
+#include "diffusion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        gasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        constGasThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        gasThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        constGasThermoPhysics,
+        rhoCombustionModel
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/noCombustion/noCombustion.C

@@ -59,7 +59,7 @@ template<class CombThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::combustionModels::noCombustion<CombThermoType>::R
 (
-    const volScalarField& Y
+    volScalarField& Y
 ) const
 {
     tmp<fvScalarMatrix> tSu

src/combustionModels/noCombustion/noCombustion.H

@@ -84,8 +84,8 @@ public:
             //- Correct combustion rate
             virtual void correct();
 
-            //- Fuel consumption rate matrix.
+            //- Fuel consumption rate matrix
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;

src/combustionModels/singleStepCombustion/singleStepCombustion.C

@@ -53,7 +53,8 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
         ),
         this->mesh(),
         dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-    )
+    ),
+    semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
 {
     if (isA<singleStepReactingMixture<ThermoType> >(this->thermo()))
     {
@@ -79,6 +80,15 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
             << "Please select a thermo package based on "
             << "singleStepReactingMixture" << exit(FatalError);
     }
+
+    if (semiImplicit_)
+    {
+        Info<< "Combustion mode: semi-implicit" << endl;
+    }
+    else
+    {
+        Info<< "Combustion mode: explicit" << endl;
+    }
 }
 
 
@@ -95,17 +105,28 @@ template<class CombThermoType, class ThermoType>
 Foam::tmp<Foam::fvScalarMatrix>
 singleStepCombustion<CombThermoType, ThermoType>::R
 (
-    const volScalarField& Y
+    volScalarField& Y
 ) const
 {
     const label specieI = this->thermoPtr_->composition().species()[Y.name()];
 
-    const volScalarField wSpecie
+    volScalarField wSpecie
     (
         wFuel_*singleMixturePtr_->specieStoichCoeffs()[specieI]
     );
 
-    return wSpecie + fvm::Sp(0.0*wSpecie, Y);
+    if (semiImplicit_)
+    {
+        const label fNorm = singleMixturePtr_->specieProd()[specieI];
+        const volScalarField fres(singleMixturePtr_->fres(specieI));
+        wSpecie /= max(fNorm*(Y - fres), 1e-2);
+
+        return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
+    }
+    else
+    {
+        return wSpecie + fvm::Sp(0.0*wSpecie, Y);
+    }
 }
 
 
@@ -114,7 +135,8 @@ Foam::tmp<Foam::volScalarField>
 singleStepCombustion< CombThermoType, ThermoType>::Sh() const
 {
     const label fuelI = singleMixturePtr_->fuelIndex();
-    const volScalarField& YFuel = this->thermoPtr_->composition().Y(fuelI);
+    volScalarField& YFuel =
+        const_cast<volScalarField&>(this->thermoPtr_->composition().Y(fuelI));
 
     return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
 }

src/combustionModels/singleStepCombustion/singleStepCombustion.H

@@ -71,6 +71,9 @@ protected:
         //- Fuel consumption rate
         volScalarField wFuel_;
 
+        //- Semi-implicit (true) or explicit (false) treatment
+        bool semiImplicit_;
+
 
 public:
 
@@ -89,7 +92,7 @@ public:
         // Evolution
 
             //- Fuel consumption rate matrix
-            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
             //- Heat release rate calculated from fuel consumption rate matrix
             virtual tmp<volScalarField> dQ() const;

src/fieldSources/basicSource/interRegionHeatTransferModel/interRegionHeatTransferModel/interRegionHeatTransferModel.C

@@ -68,17 +68,11 @@ void Foam::interRegionHeatTransferModel::check()
         }
     }
 
-    if(!secSourceFound)
+    if (!secSourceFound)
     {
         FatalErrorIn
         (
-            "constantHeatTransfer::interRegionHeatTransferModel"
+            "constantHeatTransfer::interRegionHeatTransferModel::check()"
-            "("
-            "   const word& name,"
-            "   const word& modelType,"
-            "   const dictionary& dict,"
-            "   const fvMesh& mesh"
-            ")"
         )   << "Secondary source name not found" << secondarySourceName_
             << " in region " << secondaryMesh.name()
             << nl
@@ -125,9 +119,12 @@ Foam::interRegionHeatTransferModel::interRegionHeatTransferModel
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
 Foam::interRegionHeatTransferModel::~interRegionHeatTransferModel()
 {}
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
@@ -222,7 +219,7 @@ void Foam::interRegionHeatTransferModel::addSup
                     << endl;
             }
         }
-        else if(h.dimensions() == dimTemperature)
+        else if (h.dimensions() == dimTemperature)
         {
             eEqn += htc_*Tmapped - fvm::Sp(htc_, h);
 
@@ -269,4 +266,6 @@ bool Foam::interRegionHeatTransferModel::read(const dictionary& dict)
         return false;
     }
 }
+
+
 // ************************************************************************* //

src/lagrangian/coalCombustion/Make/files

@@ -1,4 +1,5 @@
 /* Coal parcel and sub-models */
 coalParcel/makeCoalParcelSubmodels.C
+coalCloudList/coalCloudList.C
 
 LIB = $(FOAM_LIBBIN)/libcoalCombustion

src/lagrangian/coalCombustion/coalCloudList/coalCloudList.C

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "coalCloudList.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::coalCloudList::coalCloudList
+(
+    const volScalarField& rho,
+    const volVectorField& U,
+    const dimensionedVector& g,
+    const SLGThermo& slgThermo
+)
+:
+    PtrList<coalCloud>(),
+    mesh_(rho.mesh())
+{
+    IOdictionary props
+    (
+        IOobject
+        (
+            "coalCloudList",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ
+        )
+    );
+
+    const wordHashSet cloudNames(wordList(props.lookup("clouds")));
+
+    setSize(cloudNames.size());
+
+    label i = 0;
+    forAllConstIter(wordHashSet, cloudNames, iter)
+    {
+        const word& name = iter.key();
+
+        Info<< "creating cloud: " << name << endl;
+
+        set
+        (
+            i++,
+            new coalCloud
+            (
+                name,
+                rho,
+                U,
+                g,
+                slgThermo
+            )
+        );
+
+        Info<< endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::coalCloudList::evolve()
+{
+    forAll(*this, i)
+    {
+        operator[](i).evolve();
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/coalCloudList/coalCloudList.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef coalCloudList_H
+#define coalCloudList_H
+
+#include "coalCloud.H"
+#include "volFieldsFwd.H"
+#include "fvMatricesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class coalCloudList Declaration
+\*---------------------------------------------------------------------------*/
+
+class coalCloudList
+:
+    public PtrList<coalCloud>
+{
+private:
+
+    //- Reference to the mesh
+    const fvMesh& mesh_;
+
+
+public:
+
+    // Constructor
+
+        coalCloudList
+        (
+            const volScalarField& rho,
+            const volVectorField& U,
+            const dimensionedVector& g,
+            const SLGThermo& slgThermo
+        );
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Evolve the cloud collection
+            void evolve();
+
+
+        // Source terms
+
+            //- Return const reference to momentum source
+            inline tmp<DimensionedField<vector, volMesh> > UTrans() const;
+
+            //- Return tmp momentum source term
+            inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
+
+            //- Sensible enthalpy transfer [J/kg]
+            inline tmp<DimensionedField<scalar, volMesh> > hsTrans() const;
+
+            //- Return sensible enthalpy source term [J/kg/m3/s]
+            inline tmp<fvScalarMatrix> Sh(volScalarField& hs) const;
+
+
+            //- Return mass source term for specie i - specie eqn
+            inline tmp<fvScalarMatrix> SYi
+            (
+                const label i,
+                volScalarField& Yi
+            ) const;
+
+            //- Return total mass transfer [kg/m3]
+            inline tmp<DimensionedField<scalar, volMesh> > rhoTrans() const;
+
+            //- Return tmp total mass source for carrier phase
+            //  - fully explicit
+            inline tmp<DimensionedField<scalar, volMesh> > Srho() const;
+
+            //- Return tmp total mass source for carrier phase specie i
+            //  - fully explicit
+            inline tmp<DimensionedField<scalar, volMesh> > Srho
+            (
+                const label i
+            ) const;
+
+            //- Return total mass source term [kg/m3/s]
+            inline tmp<fvScalarMatrix> Srho(volScalarField& rho) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "coalCloudListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/coalCloudList/coalCloudListI.H

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMatrices.H"
+#include "volFields.H"
+#include "DimensionedField.H"
+
+Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
+Foam::coalCloudList::UTrans() const
+{
+    tmp<DimensionedField<vector, volMesh> > tfld
+    (
+        new DimensionedField<vector, volMesh>
+        (
+            IOobject
+            (
+                "UTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
+        )
+    );
+
+    DimensionedField<vector, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).UTrans();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
+(
+    volVectorField& U
+) const
+{
+    tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
+    fvVectorMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).SU(U);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::hsTrans() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "hsTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).hsTrans();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
+(
+    volScalarField& hs
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).Sh(hs);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
+(
+    const label ii,
+    volScalarField& Yi
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).SYi(ii, Yi);
+    }
+
+    return tfvm;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::rhoTrans() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        forAll(operator[](i).rhoTrans(), j)
+        {
+            fld += operator[](i).rhoTrans()[j];
+        }
+    }
+
+    return tfld;
+}
+
+
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::Srho() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, i)
+    {
+        fld += operator[](i).Srho();
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::coalCloudList::Srho
+(
+    const label i
+) const
+{
+    tmp<DimensionedField<scalar, volMesh> > tfld
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "rhoTransEff",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& fld = tfld();
+
+    forAll(*this, j)
+    {
+        fld += operator[](j).Srho(i);
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
+(
+    volScalarField& rho
+) const
+{
+    tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
+    fvScalarMatrix& fvm = tfvm();
+
+    forAll(*this, i)
+    {
+        fvm += operator[](i).Srho(rho);
+    }
+
+    return tfvm;
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ License
 #include "NoSurfaceReaction.H"
 #include "COxidationDiffusionLimitedRate.H"
 #include "COxidationKineticDiffusionLimitedRate.H"
+#include "COxidationHurtMitchell.H"
 #include "COxidationMurphyShaddix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -43,6 +44,7 @@ License
         COxidationKineticDiffusionLimitedRate,                                \
         CloudType                                                             \
     );                                                                        \
+    makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType);          \
     makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "COxidationHurtMitchell.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    SurfaceReactionModel<CloudType>(dict, owner, typeName),
+    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
+    CsLocalId_(-1),
+    ashLocalId_(-1),
+    O2GlobalId_(owner.composition().globalCarrierId("O2")),
+    CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
+    WC_(0.0),
+    WO2_(0.0),
+    HcCO2_(0.0),
+    heatOfReaction_(-1.0)
+{
+    // Determine Cs and ash ids
+    label idSolid = owner.composition().idSolid();
+    CsLocalId_ = owner.composition().localId(idSolid, "C");
+    ashLocalId_ = owner.composition().localId(idSolid, "ash", true);
+
+    // Set local copies of thermo properties
+    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
+    const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
+    WC_ = WCO2 - WO2_;
+
+    HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+
+    if (this->coeffDict().readIfPresent("heatOfReaction", heatOfReaction_))
+    {
+        Info<< "    Using user specified heat of reaction: "
+            << heatOfReaction_ << " [J/kg]" << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
+(
+    const COxidationHurtMitchell<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    Sb_(srm.Sb_),
+    CsLocalId_(srm.CsLocalId_),
+    ashLocalId_(srm.ashLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    CO2GlobalId_(srm.CO2GlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_),
+    HcCO2_(srm.HcCO2_),
+    heatOfReaction_(srm.heatOfReaction_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::COxidationHurtMitchell<CloudType>::~COxidationHurtMitchell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::COxidationHurtMitchell<CloudType>::calculate
+(
+    const scalar dt,
+    const label cellI,
+    const scalar d,
+    const scalar T,
+    const scalar Tc,
+    const scalar pc,
+    const scalar rhoc,
+    const scalar mass,
+    const scalarField& YGas,
+    const scalarField& YLiquid,
+    const scalarField& YSolid,
+    const scalarField& YMixture,
+    const scalar N,
+    scalarField& dMassGas,
+    scalarField& dMassLiquid,
+    scalarField& dMassSolid,
+    scalarField& dMassSRCarrier
+) const
+{
+    const label idGas = CloudType::parcelType::GAS;
+    const label idSolid = CloudType::parcelType::SLD;
+    const scalar Ychar = YMixture[idSolid]*YSolid[CsLocalId_];
+
+    // Surface combustion until combustible fraction is consumed
+    if (Ychar < SMALL)
+    {
+        return 0.0;
+    }
+
+    const SLGThermo& thermo = this->owner().thermo();
+
+    // Local mass fraction of O2 in the carrier phase
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
+
+    // No combustion if no oxygen present
+    if (YO2 < SMALL)
+    {
+        return 0.0;
+    }
+
+    // Conversion from [g/cm^2) to [kg/m^2]
+    const scalar convSI = 1000.0/10000.0;
+
+    // Universal gas constant in [kcal/mol/K]
+    const scalar RRcal = 1985.877534;
+
+    // Dry mass fraction
+    scalar Ydaf = YMixture[idGas] + YMixture[idSolid];
+    if (ashLocalId_ != -1)
+    {
+        Ydaf -= YMixture[idSolid]*YSolid[ashLocalId_];
+    }
+
+    // Char percentage
+    const scalar charPrc = Ychar/Ydaf*100.0;
+
+    // Particle surface area
+    const scalar Ap = constant::mathematical::pi*sqr(d);
+
+    // Far field partial pressure O2 [Pa]
+    // Note: Should really use the surface partial pressure
+    const scalar ppO2 = max(0.0, rhoc*YO2/WO2_*specie::RR*Tc);
+
+    // Activation energy [kcal/mol]
+    const scalar E = -5.94 + 0.355*charPrc;
+
+    // Pre-exponential factor [g/(cm^2.s.atm^0.5)]
+    const scalar lnK1750 = 2.8 - 0.0758*charPrc;
+    const scalar A = exp(lnK1750 + E/RRcal/1750.0);
+
+    // Kinetic rate of char oxidation [g/(cm^2.s.atm^0.5)]
+    const scalar Rk = A*exp(-E/(RRcal*T));
+
+    // Molar reaction rate per unit surface area [kmol/(m^2.s)]
+    const scalar qCsLim = mass*Ychar/(WC_*Ap*dt);
+    const scalar qCs = min(convSI*Rk*Foam::sqrt(ppO2/101325.0), qCsLim);
+
+    // Calculate the number of molar units reacted [kmol]
+    const scalar dOmega = qCs*Ap*dt;
+
+    // Add to carrier phase mass transfer
+    dMassSRCarrier[O2GlobalId_] += -dOmega*Sb_*WO2_;
+    dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + Sb_*WO2_);
+
+    // Add to particle mass transfer
+    dMassSolid[CsLocalId_] += dOmega*WC_;
+
+
+    // Return the heat of reaction [J]
+    // note: carrier sensible enthalpy exchange handled via change in mass
+    if (heatOfReaction_ < 0)
+    {
+        const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
+        return dOmega*(WC_*HsC - (WC_ + Sb_*WO2_)*HcCO2_);
+    }
+    else
+    {
+        return dOmega*WC_*heatOfReaction_;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationHurtMitchell/COxidationHurtMitchell.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    COxidationHurtMitchell
+
+Description
+    Char oxidation model given by Hurt and Mitchell:
+
+    Based on the reference:
+        Hurt R. and Mitchell R., "Unified high-temperature char combustion
+        kinetics for a suite of coals of various rank", 24th Symposium in
+        Combustion, The Combustion Institute, 1992, p 1243-1250
+
+    Model specifies the rate of char combustion.
+
+        C(s) + Sb*O2 -> CO2
+
+    where Sb is the stoichiometry of the reaction
+
+    Model validity:
+        Gas temperature: Tc > 1500 K
+        Particle sizes:  75 um -> 200 um
+        Pox > 0.3 atm        
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef COxidationHurtMitchell_H
+#define COxidationHurtMitchell_H
+
+#include "SurfaceReactionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+template<class CloudType>
+class COxidationHurtMitchell;
+
+/*---------------------------------------------------------------------------*\
+                  Class COxidationHurtMitchell Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class COxidationHurtMitchell
+:
+    public SurfaceReactionModel<CloudType>
+{
+    // Private data
+
+        // Model constants
+
+            //- Stoichiometry of reaction
+            const scalar Sb_;
+
+
+        // Addressing
+
+            //- Cs position in global/local lists
+            label CsLocalId_;
+
+            //- Ash position in global/local lists
+            label ashLocalId_;
+
+            //- O2 position in global list
+            label O2GlobalId_;
+
+            //- CO2 positions in global list
+            label CO2GlobalId_;
+
+
+        // Local copies of thermo properties
+
+            //- Molecular weight of C [kg/kmol]
+            scalar WC_;
+
+            //- Molecular weight of O2 [kg/kmol]
+            scalar WO2_;
+
+            //- Formation enthalpy for CO2 [J/kg]
+            scalar HcCO2_;
+
+            //- Heat of reaction [J/kg] (optional)
+            scalar heatOfReaction_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("COxidationHurtMitchell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        COxidationHurtMitchell
+        (
+            const dictionary& dict,
+            CloudType& owner
+        );
+
+        //- Construct copy
+        COxidationHurtMitchell
+        (
+            const COxidationHurtMitchell<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationHurtMitchell<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~COxidationHurtMitchell();
+
+
+    // Member Functions
+
+        //- Update surface reactions
+        virtual scalar calculate
+        (
+            const scalar dt,
+            const label cellI,
+            const scalar d,
+            const scalar T,
+            const scalar Tc,
+            const scalar pc,
+            const scalar rhoc,
+            const scalar mass,
+            const scalarField& YGas,
+            const scalarField& YLiquid,
+            const scalarField& YSolid,
+            const scalarField& YMixture,
+            const scalar N,
+            scalarField& dMassGas,
+            scalarField& dMassLiquid,
+            scalarField& dMassSolid,
+            scalarField& dMassSRCarrier
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "COxidationHurtMitchell.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/dsmc/submodels/BinaryCollisionModel/NoBinaryCollision/NoBinaryCollision.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,15 +68,13 @@ Foam::scalar Foam::NoBinaryCollision<CloudType>::sigmaTcR
     (
         "Foam::scalar Foam::NoBinaryCollision<CloudType>::sigmaTcR"
         "("
-            "label typeIdP,"
+            "const typename CloudType::parcelType&, "
-            "label typeIdQ,"
+            "const typename CloudType::parcelType"
-            "const vector& UP,"
-            "const vector& UQ"
         ") const"
     )
         << "sigmaTcR called on NoBinaryCollision model, this should "
         << "not happen, this is not an actual model." << nl
-        << "Enclose calls to sigmaTcR within a if(binaryCollision().active()) "
+        << "Enclose calls to sigmaTcR within a if (binaryCollision().active()) "
         << " check."
         << abort(FatalError);
 

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

@@ -661,7 +661,7 @@ void Foam::KinematicCloud<CloudType>::evolve()
 
 template<class CloudType>
 template<class TrackData>
-void  Foam::KinematicCloud<CloudType>::motion(TrackData& td)
+void Foam::KinematicCloud<CloudType>::motion(TrackData& td)
 {
     td.part() = TrackData::tpLinearTrack;
     CloudType::move(td,  solution_.trackTime());
@@ -670,6 +670,159 @@ void  Foam::KinematicCloud<CloudType>::motion(TrackData& td)
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::patchData
+(
+    const parcelType& p,
+    const polyPatch& pp,
+    const scalar trackFraction,
+    const tetIndices& tetIs,
+    vector& nw,
+    vector& Up
+) const
+{
+    label patchI = pp.index();
+    label patchFaceI = pp.whichFace(p.face());
+
+    vector n = tetIs.faceTri(mesh_).normal();
+    n /= mag(n);
+
+    vector U = U_.boundaryField()[patchI][patchFaceI];
+
+    // Unless the face is rotating, the required normal is n;
+    nw = n;
+
+    if (!mesh_.moving())
+    {
+        // Only wall patches may have a non-zero wall velocity from
+        // the velocity field when the mesh is not moving.
+
+        if (isA<wallPolyPatch>(pp))
+        {
+            Up = U;
+        }
+        else
+        {
+            Up = vector::zero;
+        }
+    }
+    else
+    {
+        vector U00 = U_.oldTime().boundaryField()[patchI][patchFaceI];
+
+        vector n00 = tetIs.oldFaceTri(mesh_).normal();
+
+        // Difference in normal over timestep
+        vector dn = vector::zero;
+
+        if (mag(n00) > SMALL)
+        {
+            // If the old normal is zero (for example in layer
+            // addition) then use the current normal, meaning that the
+            // motion can only be translational, and dn remains zero,
+            // otherwise, calculate dn:
+
+            n00 /= mag(n00);
+
+            dn = n - n00;
+        }
+
+        // Total fraction through the timestep of the motion,
+        // including stepFraction before the current tracking step
+        // and the current trackFraction
+        // i.e.
+        // let s = stepFraction, t = trackFraction
+        // Motion of x in time:
+        // |-----------------|---------|---------|
+        // x00               x0        xi        x
+        //
+        // where xi is the correct value of x at the required
+        // tracking instant.
+        //
+        // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
+        //
+        // i.e. the motion covered by previous tracking portions
+        // within this timestep, and
+        //
+        // xi = x0 + t*(x - x0)
+        //    = t*x + (1 - t)*x0
+        //    = t*x + (1 - t)*(s*x + (1 - s)*x00)
+        //    = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
+        //
+        // let m = (s + t - s*t)
+        //
+        // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
+        //
+        // In the same form as before.
+
+        scalar m =
+            p.stepFraction()
+          + trackFraction
+          - (p.stepFraction()*trackFraction);
+
+        // When the mesh is moving, the velocity field on wall patches
+        // will contain the velocity associated with the motion of the
+        // mesh, in which case it is interpolated in time using m.
+        // For other patches the face velocity will need to be
+        // reconstructed from the face centre motion.
+
+        const vector& Cf = mesh_.faceCentres()[p.face()];
+
+        vector Cf00 = mesh_.faces()[p.face()].centre(mesh_.oldPoints());
+
+        if (isA<wallPolyPatch>(pp))
+        {
+            Up = U00 + m*(U - U00);
+        }
+        else
+        {
+            Up = (Cf - Cf00)/mesh_.time().deltaTValue();
+        }
+
+        if (mag(dn) > SMALL)
+        {
+            // Rotational motion, nw requires interpolation and a
+            // rotational velocity around face centre correction to Up
+            // is required.
+
+            nw = n00 + m*dn;
+
+            // Cf at tracking instant
+            vector Cfi = Cf00 + m*(Cf - Cf00);
+
+            // Normal vector cross product
+            vector omega = (n00 ^ n);
+
+            scalar magOmega = mag(omega);
+
+            // magOmega = sin(angle between unit normals)
+            // Normalise omega vector by magOmega, then multiply by
+            // angle/dt to give the correct angular velocity vector.
+            omega *= Foam::asin(magOmega)/(magOmega*mesh_.time().deltaTValue());
+
+            // Project position onto face and calculate this position
+            // relative to the face centre.
+            vector facePos =
+                p.position()
+              - ((p.position() - Cfi) & nw)*nw
+              - Cfi;
+
+            Up += (omega ^ facePos);
+        }
+
+        // No further action is required if the motion is
+        // translational only, nw and Up have already been set.
+    }
+}
+
+
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::updateMesh()
+{
+    injectors_.updateMesh();
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
 {
@@ -678,6 +831,8 @@ void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    updateMesh();
 }
 
 

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

@@ -554,9 +554,24 @@ public:
             template<class TrackData>
             void motion(TrackData& td);
 
+            //- Calculate the patch normal and velocity to interact with,
+            //  accounting for patch motion if required.
+            void patchData
+            (
+                const parcelType& p,
+                const polyPatch& pp,
+                const scalar trackFraction,
+                const tetIndices& tetIs,
+                vector& normal,
+                vector& Up
+            ) const;
+
 
         // Mapping
 
+            //- Update mesh
+            void updateMesh();
+
             //- Remap the cells of particles corresponding to the
             //  mesh topology change with a default tracking data object
             virtual void autoMap(const mapPolyMesh&);

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

@@ -87,8 +87,6 @@ void Foam::ReactingCloud<CloudType>::cloudReset(ReactingCloud<CloudType>& c)
 
     compositionModel_.reset(c.compositionModel_.ptr());
     phaseChangeModel_.reset(c.phaseChangeModel_.ptr());
-
-    dMassPhaseChange_ = c.dMassPhaseChange_;
 }
 
 
@@ -111,8 +109,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     constProps_(this->particleProperties(), this->solution().active()),
     compositionModel_(NULL),
     phaseChangeModel_(NULL),
-    rhoTrans_(thermo.carrier().species().size()),
+    rhoTrans_(thermo.carrier().species().size())
-    dMassPhaseChange_(0.0)
 {
     if (this->solution().active())
     {
@@ -167,8 +164,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     constProps_(c.constProps_),
     compositionModel_(c.compositionModel_->clone()),
     phaseChangeModel_(c.phaseChangeModel_->clone()),
-    rhoTrans_(c.rhoTrans_.size()),
+    rhoTrans_(c.rhoTrans_.size())
-    dMassPhaseChange_(c.dMassPhaseChange_)
 {
     forAll(c.rhoTrans_, i)
     {
@@ -209,8 +205,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
     compositionModel_(c.compositionModel_->clone()),
 //    compositionModel_(NULL),
     phaseChangeModel_(NULL),
-    rhoTrans_(0),
+    rhoTrans_(0)
-    dMassPhaseChange_(0.0)
 {}
 
 
@@ -350,6 +345,8 @@ void Foam::ReactingCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    this->updateMesh();
 }
 
 

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

@@ -121,12 +121,6 @@ protected:
             PtrList<DimensionedField<scalar, volMesh> > rhoTrans_;
 
 
-        // Check
-
-            //- Total mass transferred to continuous phase via phase change
-            scalar dMassPhaseChange_;
-
-
     // Protected Member Functions
 
         // New parcel helper functions

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

@@ -264,6 +264,8 @@ void Foam::ReactingMultiphaseCloud<CloudType>::autoMap
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    this->updateMesh();
 }
 
 

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

@@ -51,7 +51,69 @@ void Foam::ThermoCloud<CloudType>::setModels()
         ).ptr()
     );
 
-    this->subModelProperties().lookup("radiation") >> radiation_;
+    if (this->solution().coupled())
+    {
+        this->subModelProperties().lookup("radiation") >> radiation_;
+    }
+
+    if (radiation_)
+    {
+        radAreaP_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaP",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimArea, 0.0)
+            )
+        );
+
+        radT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", pow4(dimTemperature), 0.0)
+            )
+        );
+
+        radAreaPT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaPT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar
+                (
+                    "zero",
+                    sqr(dimLength)*pow4(dimTemperature),
+                    0.0
+                )
+            )
+        );
+    }
 }
 
 
@@ -98,13 +160,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_
     (
         new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
-                this->name() + "hsTrans",
+                this->name() + "::hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -120,7 +185,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "hsCoeff",
+                this->name() + "::hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -165,13 +230,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(c.heatTransferModel_->clone()),
     TIntegrator_(c.TIntegrator_->clone()),
     radiation_(c.radiation_),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_
     (
         new DimensionedField<scalar, volMesh>
         (
             IOobject
             (
-                this->name() + "hsTrans",
+                this->name() + "::hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -187,7 +255,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "hsCoeff",
+                this->name() + "::hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -197,7 +265,61 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             c.hsCoeff()
         )
     )
-{}
+{
+    if (radiation_)
+    {
+        radAreaP_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaP",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radAreaP()
+            )
+        );
+
+        radT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radT4()
+            )
+        );
+
+        radAreaPT4_.reset
+        (
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    this->name() + "::radAreaPT4",
+                    this->db().time().timeName(),
+                    this->db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                c.radAreaPT4()
+            )
+        );
+    }
+}
 
 
 template<class CloudType>
@@ -218,6 +340,9 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
     heatTransferModel_(NULL),
     TIntegrator_(NULL),
     radiation_(false),
+    radAreaP_(NULL),
+    radT4_(NULL),
+    radAreaPT4_(NULL),
     hsTrans_(NULL),
     hsCoeff_(NULL)
 {}
@@ -285,6 +410,13 @@ void Foam::ThermoCloud<CloudType>::resetSourceTerms()
     CloudType::resetSourceTerms();
     hsTrans_->field() = 0.0;
     hsCoeff_->field() = 0.0;
+
+    if (radiation_)
+    {
+        radAreaP_->field() = 0.0;
+        radT4_->field() = 0.0;
+        radAreaPT4_->field() = 0.0;
+    }
 }
 
 
@@ -298,6 +430,13 @@ void Foam::ThermoCloud<CloudType>::relaxSources
 
     this->relax(hsTrans_(), cloudOldTime.hsTrans(), "h");
     this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "h");
+
+    if (radiation_)
+    {
+        this->relax(radAreaP_(), cloudOldTime.radAreaP(), "radiation");
+        this->relax(radT4_(), cloudOldTime.radT4(), "radiation");
+        this->relax(radAreaPT4_(), cloudOldTime.radAreaPT4(), "radiation");
+    }
 }
 
 
@@ -308,6 +447,13 @@ void Foam::ThermoCloud<CloudType>::scaleSources()
 
     this->scale(hsTrans_(), "h");
     this->scale(hsCoeff_(), "h");
+
+    if (radiation_)
+    {
+        this->scale(radAreaP_(), "radiation");
+        this->scale(radT4_(), "radiation");
+        this->scale(radAreaPT4_(), "radiation");
+    }
 }
 
 
@@ -341,6 +487,8 @@ void Foam::ThermoCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
     tdType td(*this);
 
     Cloud<parcelType>::template autoMap<tdType>(td, mapper);
+
+    this->updateMesh();
 }
 
 

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

@@ -132,6 +132,15 @@ protected:
             //- Include radiation
             Switch radiation_;
 
+            //- Radiation sum of parcel projected areas
+            autoPtr<DimensionedField<scalar, volMesh> > radAreaP_;
+
+            //- Radiation sum of parcel temperature^4
+            autoPtr<DimensionedField<scalar, volMesh> > radT4_;
+
+            //- Radiation sum of parcel projected areas * temperature^4
+            autoPtr<DimensionedField<scalar, volMesh> > radAreaPT4_;
+
 
         // Sources
 
@@ -244,6 +253,26 @@ public:
                  //- Radiation flag
                 inline bool radiation() const;
 
+                //- Radiation sum of parcel projected areas [m2]
+                inline DimensionedField<scalar, volMesh>& radAreaP();
+
+                //- Radiation sum of parcel projected areas [m2]
+                inline const DimensionedField<scalar, volMesh>&
+                    radAreaP() const;
+
+                //- Radiation sum of parcel temperature^4 [K4]
+                inline DimensionedField<scalar, volMesh>& radT4();
+
+                //- Radiation sum of parcel temperature^4 [K4]
+                inline const DimensionedField<scalar, volMesh>& radT4() const;
+
+                //- Radiation sum of parcel projected area*temperature^4 [m2K4]
+                inline DimensionedField<scalar, volMesh>& radAreaPT4();
+
+                //- Radiation sum of parcel temperature^4 [m2K4]
+                inline const DimensionedField<scalar, volMesh>&
+                    radAreaPT4() const;
+
 
             // Sources
 

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

@@ -89,6 +89,114 @@ inline bool Foam::ThermoCloud<CloudType>::radiation() const
 }
 
 
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaP()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaP()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaP_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaP() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaP()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaP_();
+}
+
+
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radT4()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radT4_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radT4() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radT4_();
+}
+
+
+template<class CloudType>
+inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaPT4()
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaPT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaPT4_();
+}
+
+
+template<class CloudType>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::ThermoCloud<CloudType>::radAreaPT4() const
+{
+    if (!radiation_)
+    {
+        FatalErrorIn
+        (
+            "inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
+            "Foam::ThermoCloud<CloudType>::radAreaPT4()"
+        )   << "Radiation field requested, but radiation model not active"
+            << abort(FatalError);
+    }
+
+    return radAreaPT4_();
+}
+
+
 template<class CloudType>
 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 Foam::ThermoCloud<CloudType>::hsTrans()
@@ -182,21 +290,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::Ep() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
+    if (radiation_)
-    if (radiation_ && this->solution().coupled())
     {
         scalarField& Ep = tEp().internalField();
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
+        const scalarField& sumAreaPT4 = radAreaPT4_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
+        Ep = sumAreaPT4*epsilon*physicoChemical::sigma.value()/V/dt;
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
-        }
-
-        Ep *= epsilon*physicoChemical::sigma.value()/V;
     }
 
     return tEp;
@@ -224,21 +326,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::ap() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
+    if (radiation_)
-    if (radiation_ && this->solution().coupled())
     {
         scalarField& ap = tap().internalField();
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
+        const scalarField& sumAreaP = radAreaP_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
+        ap = sumAreaP*epsilon/V/dt;
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            ap[cellI] += p.nParticle()*p.areaP();
-        }
-
-        ap *= epsilon/V;
     }
 
     return tap;
@@ -267,23 +363,16 @@ Foam::ThermoCloud<CloudType>::sigmap() const
         )
     );
 
-    // Need to check if coupled as field is created on-the-fly
+    if (radiation_)
-    if (radiation_ && this->solution().coupled())
     {
         scalarField& sigmap = tsigmap().internalField();
-
+        const scalar dt = this->db().time().deltaTValue();
         const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
         const scalar f = constProps_.f0();
+        const scalarField& sumAreaP = radAreaP_->field();
 
-        forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
+        sigmap *= sumAreaP*(1.0 - f)*(1.0 - epsilon)/V/dt;
-        {
-            const parcelType& p = iter();
-            const label cellI = p.cell();
-            sigmap[cellI] += p.nParticle()*p.areaP();
-        }
-
-        sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
     }
 
     return tsigmap;

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -363,7 +363,7 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
     // Invoke post-processing model
-    td.cloud().functions().postPatch(p, patchI, pp.whichFace(p.face()));
+    td.cloud().functions().postPatch(p, pp, trackFraction, tetIs);
 
     // Invoke surface film model
     if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -268,12 +268,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     (
         td,
         dt,
+        this->age_,
         Ts,
         d0,
         T0,
         mass0,
         this->mass0_,
         YMix[GAS]*YGas_,
+        YMix[LIQ]*YLiquid_,
+        YMix[SLD]*YSolid_,
         canCombust_,
         dMassDV,
         Sh,
@@ -466,6 +469,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 
@@ -476,12 +489,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 (
     TrackData& td,
     const scalar dt,
+    const scalar age,
     const scalar Ts,
     const scalar d,
     const scalar T,
     const scalar mass,
     const scalar mass0,
     const scalarField& YGasEff,
+    const scalarField& YLiquidEff,
+    const scalarField& YSolidEff,
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
@@ -510,10 +526,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     td.cloud().devolatilisation().calculate
     (
         dt,
+        age,
         mass0,
         mass,
         T,
         YGasEff,
+        YLiquidEff,
+        YSolidEff,
         canCombust,
         dMassDV
     );
@@ -638,7 +657,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType(p),
     YGas_(p.YGas_),
     YLiquid_(p.YLiquid_),
-    YSolid_(p.YSolid_)
+    YSolid_(p.YSolid_),
+    canCombust_(p.canCombust_)
 {}
 
 
@@ -652,7 +672,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
     ParcelType(p, mesh),
     YGas_(p.YGas_),
     YLiquid_(p.YLiquid_),
-    YSolid_(p.YSolid_)
+    YSolid_(p.YSolid_),
+    canCombust_(p.canCombust_)
 {}
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -196,12 +196,15 @@ protected:
         (
             TrackData& td,
             const scalar dt,           // timestep
-            const scalar Ts,           // Surface temperature
+            const scalar age,          // age
+            const scalar Ts,           // surface temperature
             const scalar d,            // diameter
             const scalar T,            // temperature
             const scalar mass,         // mass
             const scalar mass0,        // mass (initial on injection)
-            const scalarField& YGasEff,// Gas component mass fractions
+            const scalarField& YGasEff,// gas component mass fractions
+            const scalarField& YLiquidEff,// liquid component mass fractions
+            const scalarField& YSolidEff,// solid component mass fractions
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -458,6 +458,16 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Update sensible enthalpy transfer
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -144,7 +144,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
         Ts = td.cloud().constProps().TMin();
     }
 
-    // Assuming thermo props vary linearly with T for small dT
+    // Assuming thermo props vary linearly with T for small d(T)
     const scalar TRatio = Tc_/Ts;
 
     rhos = this->rhoc_*TRatio;
@@ -252,6 +252,16 @@ void Foam::ThermoParcel<ParcelType>::calc
 
         // Update sensible enthalpy coefficient
         td.cloud().hsCoeff()[cellI] += np0*Sph;
+
+        // Update radiation fields
+        if (td.cloud().radiation())
+        {
+            const scalar ap = this->areaP();
+            const scalar T4 = pow4(this->T_);
+            td.cloud().radAreaP()[cellI] += dt*np0*ap;
+            td.cloud().radT4()[cellI] += dt*np0*T4;
+            td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
+        }
     }
 }
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,8 +107,9 @@ template<class CloudType>
 void Foam::CloudFunctionObject<CloudType>::postPatch
 (
     const typename CloudType::parcelType&,
-    const label,
+    const polyPatch&,
-    const label
+    const scalar,
+    const tetIndices&
 )
 {
     // do nothing

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -141,8 +141,9 @@ public:
             virtual void postPatch
             (
                 const typename CloudType::parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& testIs
             );
 
             //- Post-face hook

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,13 +146,14 @@ template<class CloudType>
 void Foam::CloudFunctionObjectList<CloudType>::postPatch
 (
     const typename CloudType::parcelType& p,
-    const label patchI,
+    const polyPatch& pp,
-    const label patchFaceI
+    const scalar trackFraction,
+    const tetIndices& tetIs
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postPatch(p, patchI, patchFaceI);
+        this->operator[](i).postPatch(p, pp, trackFraction, tetIs);
     }
 }
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,8 +121,9 @@ public:
             virtual void postPatch
             (
                 const typename CloudType::parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 
             //- Post-face hook

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C

@@ -166,29 +166,32 @@ template<class CloudType>
 void Foam::ParticleErosion<CloudType>::postPatch
 (
     const parcelType& p,
-    const label patchI,
+    const polyPatch& pp,
-    const label patchFaceI
+    const scalar trackFraction,
+    const tetIndices& tetIs
 )
 {
+    const label patchI = pp.index();
+
     const label localPatchI = applyToPatch(patchI);
 
     if (localPatchI != -1)
     {
-        const fvMesh& mesh = this->owner().mesh();
+        vector nw;
+        vector Up;
 
         // patch-normal direction
-        vector nw = p.currentTetIndices().faceTri(mesh).normal();
+        this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
-        // particle direction of travel
+        // particle velocity reletive to patch
-        const vector& U = p.U();
+        const vector& U = p.U() - Up;
 
         // quick reject if particle travelling away from the patch
-        if ((-nw & U) < 0)
+        if ((nw & U) < 0)
         {
             return;
         }
 
-        nw /= mag(nw);
         const scalar magU = mag(U);
         const vector Udir = U/magU;
 
@@ -197,6 +200,7 @@ void Foam::ParticleErosion<CloudType>::postPatch
 
         const scalar coeff = p.nParticle()*p.mass()*sqr(magU)/(p_*psi_*K_);
 
+        const label patchFaceI = pp.whichFace(p.face());
         scalar& Q = QPtr_->boundaryField()[patchI][patchFaceI];
         if (tan(alpha) < K_/6.0)
         {

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,8 +126,9 @@ public:
             virtual void postPatch
             (
                 const parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 };
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -218,11 +218,14 @@ template<class CloudType>
 void Foam::PatchPostProcessing<CloudType>::postPatch
 (
     const parcelType& p,
-    const label patchI,
+    const polyPatch& pp,
-    const label
+    const scalar,
+    const tetIndices& tetIs
 )
 {
+    const label patchI = pp.index();
     const label localPatchI = applyToPatch(patchI);
+
     if (localPatchI != -1 && patchData_[localPatchI].size() < maxStoredParcels_)
     {
         times_[localPatchI].append(this->owner().time().value());

src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,8 +127,9 @@ public:
             virtual void postPatch
             (
                 const parcelType& p,
-                const label patchI,
+                const polyPatch& pp,
-                const label patchFaceI
+                const scalar trackFraction,
+                const tetIndices& tetIs
             );
 };
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C

@@ -306,8 +306,9 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                         this->owner().functions().postPatch
                         (
                             p,
-                            patchI,
+                            mesh.boundaryMesh()[patchI],
-                            patchFaceI
+                            1.0,
+                            p.currentTetIndices()
                         );
                      }
                 }

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

@@ -185,60 +185,7 @@ Foam::CellZoneInjection<CloudType>::CellZoneInjection
         )
     )
 {
-    const fvMesh& mesh = owner.mesh();
+    updateMesh();
-    const label zoneI = mesh.cellZones().findZoneID(cellZoneName_);
-
-    if (zoneI < 0)
-    {
-        FatalErrorIn
-        (
-            "Foam::CellZoneInjection<CloudType>::CellZoneInjection"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )   << "Unknown cell zone name: " << cellZoneName_
-            << ". Valid cell zones are: " << mesh.cellZones().names()
-            << nl << exit(FatalError);
-    }
-
-    const labelList& cellZoneCells = mesh.cellZones()[zoneI];
-    const label nCells = cellZoneCells.size();
-    const scalar nCellsTotal = returnReduce(nCells, sumOp<label>());
-    const scalar VCells = sum(scalarField(mesh.V(), cellZoneCells));
-    const scalar VCellsTotal = returnReduce(VCells, sumOp<scalar>());
-    Info<< "    cell zone size      = " << nCellsTotal << endl;
-    Info<< "    cell zone volume    = " << VCellsTotal << endl;
-
-    if ((nCells == 0) || (VCellsTotal*numberDensity_ < 1))
-    {
-        WarningIn
-        (
-            "Foam::CellZoneInjection<CloudType>::CellZoneInjection"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )   << "Number of particles to be added to cellZone " << cellZoneName_
-            << " is zero" << endl;
-    }
-    else
-    {
-        setPositions(cellZoneCells);
-
-        Info<< "    number density      = " << numberDensity_ << nl
-            << "    number of particles = " << positions_.size() << endl;
-
-        // Construct parcel diameters
-        diameters_.setSize(positions_.size());
-        forAll(diameters_, i)
-        {
-            diameters_[i] = sizeDistribution_->sample();
-        }
-    }
-
-    // Determine volume of particles to inject
-    this->volumeTotal_ = sum(pow3(diameters_))*constant::mathematical::pi/6.0;
 }
 
 
@@ -270,6 +217,55 @@ Foam::CellZoneInjection<CloudType>::~CellZoneInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::CellZoneInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    const fvMesh& mesh = this->owner().mesh();
+    const label zoneI = mesh.cellZones().findZoneID(cellZoneName_);
+
+    if (zoneI < 0)
+    {
+        FatalErrorIn("Foam::CellZoneInjection<CloudType>::updateMesh()")
+            << "Unknown cell zone name: " << cellZoneName_
+            << ". Valid cell zones are: " << mesh.cellZones().names()
+            << nl << exit(FatalError);
+    }
+
+    const labelList& cellZoneCells = mesh.cellZones()[zoneI];
+    const label nCells = cellZoneCells.size();
+    const scalar nCellsTotal = returnReduce(nCells, sumOp<label>());
+    const scalar VCells = sum(scalarField(mesh.V(), cellZoneCells));
+    const scalar VCellsTotal = returnReduce(VCells, sumOp<scalar>());
+    Info<< "    cell zone size      = " << nCellsTotal << endl;
+    Info<< "    cell zone volume    = " << VCellsTotal << endl;
+
+    if ((nCells == 0) || (VCellsTotal*numberDensity_ < 1))
+    {
+        WarningIn("Foam::CellZoneInjection<CloudType>::updateMesh()")
+            << "Number of particles to be added to cellZone " << cellZoneName_
+            << " is zero" << endl;
+    }
+    else
+    {
+        setPositions(cellZoneCells);
+
+        Info<< "    number density      = " << numberDensity_ << nl
+            << "    number of particles = " << positions_.size() << endl;
+
+        // Construct parcel diameters
+        diameters_.setSize(positions_.size());
+        forAll(diameters_, i)
+        {
+            diameters_[i] = sizeDistribution_->sample();
+        }
+    }
+
+    // Determine volume of particles to inject
+    this->volumeTotal_ = sum(pow3(diameters_))*constant::mathematical::pi/6.0;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::CellZoneInjection<CloudType>::timeEnd() const
 {

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

@@ -130,6 +130,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -126,17 +126,7 @@ Foam::ConeInjection<CloudType>::ConeInjection
     // Set total volume to inject
     this->volumeTotal_ = flowRateProfile_.integrate(0.0, duration_);
 
-    // Set/cache the injector cells
+    updateMesh();
-    forAll(positionAxis_, i)
-    {
-        this->findCellAtPosition
-        (
-            injectorCells_[i],
-            injectorTetFaces_[i],
-            injectorTetPts_[i],
-            positionAxis_[i].first()
-        );
-    }
 }
 
 
@@ -173,6 +163,23 @@ Foam::ConeInjection<CloudType>::~ConeInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ConeInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    forAll(positionAxis_, i)
+    {
+        this->findCellAtPosition
+        (
+            injectorCells_[i],
+            injectorTetFaces_[i],
+            injectorTetPts_[i],
+            positionAxis_[i].first()
+        );
+    }
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::ConeInjection<CloudType>::timeEnd() const
 {

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

@@ -144,6 +144,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -201,6 +201,8 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
 
     // Set total volume to inject
     this->volumeTotal_ = flowRateProfile_.integrate(0.0, duration_);
+
+    updateMesh();
 }
 
 
@@ -244,6 +246,30 @@ Foam::ConeNozzleInjection<CloudType>::~ConeNozzleInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    switch (injectionMethod_)
+    {
+        case imPoint:
+        {
+            this->findCellAtPosition
+            (
+                injectorCell_,
+                tetFaceI_,
+                tetPtI_,
+                position_
+            );
+        }
+        default:
+        {
+            // do nothing
+        }
+    }
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::ConeNozzleInjection<CloudType>::timeEnd() const
 {
@@ -342,6 +368,7 @@ void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell
                     "const label, "
                     "const scalar, "
                     "vector&, "
+                    "label&, "
                     "label&"
                 ")"
             )<< "Unknown injectionMethod type" << nl

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

@@ -212,6 +212,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -96,17 +96,7 @@ Foam::FieldActivatedInjection<CloudType>::FieldActivatedInjection
     this->volumeTotal_ =
         nParcelsPerInjector_*sum(pow3(diameters_))*pi/6.0;
 
-    // Set/cache the injector cells
+    updateMesh();
-    forAll(positions_, i)
-    {
-        this->findCellAtPosition
-        (
-            injectorCells_[i],
-            injectorTetFaces_[i],
-            injectorTetPts_[i],
-            positions_[i]
-        );
-    }
 }
 
 
@@ -142,6 +132,23 @@ Foam::FieldActivatedInjection<CloudType>::~FieldActivatedInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::FieldActivatedInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    forAll(positions_, i)
+    {
+        this->findCellAtPosition
+        (
+            injectorCells_[i],
+            injectorTetFaces_[i],
+            injectorTetPts_[i],
+            positions_[i]
+        );
+    }
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::FieldActivatedInjection<CloudType>::timeEnd() const
 {

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

@@ -150,6 +150,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -153,6 +153,13 @@ Foam::InflationInjection<CloudType>::~InflationInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::InflationInjection<CloudType>::updateMesh()
+{
+    // do nothing
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
 {
@@ -234,8 +241,8 @@ Foam::label Foam::InflationInjection<CloudType>::parcelsToInject
                 "Foam::label "
                 "Foam::InflationInjection<CloudType>::parcelsToInject"
                 "("
-                    "const scalar time0, "
+                    "const scalar, "
-                    "const scalar time1"
+                    "const scalar"
                 ")"
             )
                 << "Maximum particle split iterations ("

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

@@ -146,6 +146,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -284,6 +284,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
     SOI_(0.0),
     volumeTotal_(0.0),
     massTotal_(0.0),
+    massFlowRate_(owner.db().time(), "massFlowRate"),
     massInjected_(this->template getModelProperty<scalar>("massInjected")),
     nInjections_(this->template getModelProperty<label>("nInjections")),
     parcelsAddedTotal_
@@ -310,6 +311,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     SOI_(0.0),
     volumeTotal_(0.0),
     massTotal_(0.0),
+    massFlowRate_(owner.db().time(), "massFlowRate"),
     massInjected_(this->template getModelProperty<scalar>("massInjected")),
     nInjections_(this->template getModelProperty<scalar>("nInjections")),
     parcelsAddedTotal_
@@ -335,7 +337,8 @@ Foam::InjectionModel<CloudType>::InjectionModel
     }
     else
     {
-        this->coeffDict().lookup("massFlowRate") >> massTotal_;
+        massFlowRate_.reset(this->coeffDict());
+        massTotal_ = massFlowRate_.value(owner.db().time().value());
     }
 
     const word parcelBasisType = this->coeffDict().lookup("parcelBasisType");
@@ -384,6 +387,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
     SOI_(im.SOI_),
     volumeTotal_(im.volumeTotal_),
     massTotal_(im.massTotal_),
+    massFlowRate_(im.massFlowRate_),
     massInjected_(im.massInjected_),
     nInjections_(im.nInjections_),
     parcelsAddedTotal_(im.parcelsAddedTotal_),
@@ -403,6 +407,13 @@ Foam::InjectionModel<CloudType>::~InjectionModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::InjectionModel<CloudType>::updateMesh()
+{
+    // do nothing
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::InjectionModel<CloudType>::timeEnd() const
 {
@@ -458,7 +469,16 @@ Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
 template<class CloudType>
 Foam::scalar Foam::InjectionModel<CloudType>::averageParcelMass()
 {
-    label nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
+    label nTotal = 0.0;
+    if (this->owner().solution().transient())
+    {
+        nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
+    }
+    else
+    {
+        nTotal = parcelsToInject(0.0, 1.0);
+    }
+
     return massTotal_/nTotal;
 }
 
@@ -594,6 +614,8 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
     const polyMesh& mesh = this->owner().mesh();
     typename TrackData::cloudType& cloud = td.cloud();
 
+    massTotal_ = massFlowRate_.value(mesh.time().value());
+
     // Reset counters
     time0_ = 0.0;
     label parcelsAdded = 0;

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

@@ -53,6 +53,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "SubModelBase.H"
 #include "vector.H"
+#include "TimeDataEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -102,6 +103,9 @@ protected:
             //- Total mass to inject [kg]
             scalar massTotal_;
 
+            //- Mass flow rate profile for steady calculations
+            TimeDataEntry<scalar> massFlowRate_;
+
             //- Total mass injected to date [kg]
             scalar massInjected_;
 
@@ -246,6 +250,12 @@ public:
 
     // Member Functions
 
+        // Mapping
+
+            //- Update mesh
+            virtual void updateMesh();
+
+
         // Global information
 
             //- Return the start-of-injection time

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

@@ -168,6 +168,16 @@ Foam::scalar Foam::InjectionModelList<CloudType>::averageParcelMass()
 }
 
 
+template<class CloudType>
+void Foam::InjectionModelList<CloudType>::updateMesh()
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).updateMesh();
+    }
+}
+
+
 template<class CloudType>
 template<class TrackData>
 void Foam::InjectionModelList<CloudType>::inject(TrackData& td)

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

@@ -97,6 +97,12 @@ public:
             scalar averageParcelMass();
 
 
+        // Edit
+
+            //- Set injector locations when mesh is updated
+            void updateMesh();
+
+
         // Per-injection event functions
 
             //- Main injection loop

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

@@ -65,16 +65,7 @@ Foam::KinematicLookupTableInjection<CloudType>::KinematicLookupTableInjection
     injectorTetFaces_.setSize(injectors_.size());
     injectorTetPts_.setSize(injectors_.size());
 
-    forAll(injectors_, i)
+    updateMesh();
-    {
-        this->findCellAtPosition
-        (
-            injectorCells_[i],
-            injectorTetFaces_[i],
-            injectorTetPts_[i],
-            injectors_[i].x()
-        );
-    }
 
     // Determine volume of particles to inject
     this->volumeTotal_ = 0.0;
@@ -112,6 +103,23 @@ Foam::KinematicLookupTableInjection<CloudType>::~KinematicLookupTableInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::KinematicLookupTableInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    forAll(injectors_, i)
+    {
+        this->findCellAtPosition
+        (
+            injectorCells_[i],
+            injectorTetFaces_[i],
+            injectorTetPts_[i],
+            injectors_[i].x()
+        );
+    }
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::timeEnd() const
 {

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

@@ -129,6 +129,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -26,7 +26,6 @@ License
 #include "ManualInjection.H"
 #include "mathematicalConstants.H"
 #include "PackedBoolList.H"
-#include "Switch.H"
 
 using namespace Foam::constant::mathematical;
 
@@ -65,48 +64,13 @@ Foam::ManualInjection<CloudType>::ManualInjection
             this->coeffDict().subDict("sizeDistribution"),
             owner.rndGen()
         )
-    )
+    ),
-{
+    ignoreOutOfBounds_
-    Switch ignoreOutOfBounds
     (
         this->coeffDict().lookupOrDefault("ignoreOutOfBounds", false)
-    );
+    )
-
+{
-    label nRejected = 0;
+    updateMesh();
-
-    PackedBoolList keep(positions_.size(), true);
-
-    forAll(positions_, pI)
-    {
-        if
-        (
-            !this->findCellAtPosition
-            (
-                injectorCells_[pI],
-                injectorTetFaces_[pI],
-                injectorTetPts_[pI],
-                positions_[pI],
-                !ignoreOutOfBounds
-            )
-        )
-        {
-            keep[pI] = false;
-
-            nRejected++;
-        }
-    }
-
-    if (nRejected > 0)
-    {
-        inplaceSubset(keep, positions_);
-        inplaceSubset(keep, diameters_);
-        inplaceSubset(keep, injectorCells_);
-        inplaceSubset(keep, injectorTetFaces_);
-        inplaceSubset(keep, injectorTetPts_);
-
-        Info<< "    " << nRejected
-            << " particles ignored, out of bounds." << endl;
-    }
 
     // Construct parcel diameters
     forAll(diameters_, i)
@@ -133,7 +97,8 @@ Foam::ManualInjection<CloudType>::ManualInjection
     injectorTetFaces_(im.injectorTetFaces_),
     injectorTetPts_(im.injectorTetPts_),
     U0_(im.U0_),
-    sizeDistribution_(im.sizeDistribution_().clone().ptr())
+    sizeDistribution_(im.sizeDistribution_().clone().ptr()),
+    ignoreOutOfBounds_(im.ignoreOutOfBounds_)
 {}
 
 
@@ -146,6 +111,46 @@ Foam::ManualInjection<CloudType>::~ManualInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ManualInjection<CloudType>::updateMesh()
+{
+    label nRejected = 0;
+
+    PackedBoolList keep(positions_.size(), true);
+
+    forAll(positions_, pI)
+    {
+        if
+        (
+            !this->findCellAtPosition
+            (
+                injectorCells_[pI],
+                injectorTetFaces_[pI],
+                injectorTetPts_[pI],
+                positions_[pI],
+                !ignoreOutOfBounds_
+            )
+        )
+        {
+            keep[pI] = false;
+            nRejected++;
+        }
+    }
+
+    if (nRejected > 0)
+    {
+        inplaceSubset(keep, positions_);
+        inplaceSubset(keep, diameters_);
+        inplaceSubset(keep, injectorCells_);
+        inplaceSubset(keep, injectorTetFaces_);
+        inplaceSubset(keep, injectorTetPts_);
+
+        Info<< "    " << nRejected
+            << " particles ignored, out of bounds" << endl;
+    }
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::ManualInjection<CloudType>::timeEnd() const
 {

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

@@ -43,6 +43,7 @@ SourceFiles
 
 #include "InjectionModel.H"
 #include "distributionModel.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +51,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class ManualInjection Declaration
+                       Class ManualInjection Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -84,6 +85,9 @@ class ManualInjection
         //- Parcel size distribution model
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
+        //- Flag to suppress errors if particle injection site is out-of-bounds
+        Switch ignoreOutOfBounds_;
+
 
 public:
 
@@ -120,6 +124,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -82,7 +82,7 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
 
     duration_ = owner.db().time().userTimeToTime(duration_);
 
-    cellOwners_ = patch.faceCells();
+    updateMesh();
 
     // TODO: retrieve mean diameter from distrution model
     scalar pMeanDiameter =
@@ -139,6 +139,15 @@ Foam::PatchFlowRateInjection<CloudType>::~PatchFlowRateInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::PatchFlowRateInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    const polyPatch& patch = this->owner().mesh().boundaryMesh()[patchId_];
+    cellOwners_ = patch.faceCells();
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::PatchFlowRateInjection<CloudType>::timeEnd() const
 {

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

@@ -146,6 +146,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

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

@@ -80,11 +80,9 @@ Foam::PatchInjection<CloudType>::PatchInjection
             << nl << exit(FatalError);
     }
 
-    const polyPatch& patch = owner.mesh().boundaryMesh()[patchId_];
-
     duration_ = owner.db().time().userTimeToTime(duration_);
 
-    cellOwners_ = patch.faceCells();
+    updateMesh();
 
     label patchSize = cellOwners_.size();
     label totalPatchSize = patchSize;
@@ -125,6 +123,15 @@ Foam::PatchInjection<CloudType>::~PatchInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::PatchInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    const polyPatch& patch = this->owner().mesh().boundaryMesh()[patchId_];
+    cellOwners_ = patch.faceCells();
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::PatchInjection<CloudType>::timeEnd() const
 {

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

@@ -129,6 +129,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C

@@ -233,7 +233,7 @@ bool Foam::LocalInteraction<CloudType>::correct
                 vector nw;
                 vector Up;
 
-                this->patchData(p, pp, trackFraction, tetIs, nw, Up);
+                this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
                 // Calculate motion relative to patch velocity
                 U -= Up;

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -179,159 +179,6 @@ bool Foam::PatchInteractionModel<CloudType>::correct
 }
 
 
-template<class CloudType>
-void Foam::PatchInteractionModel<CloudType>::patchData
-(
-    typename CloudType::parcelType& p,
-    const polyPatch& pp,
-    const scalar trackFraction,
-    const tetIndices& tetIs,
-    vector& nw,
-    vector& Up
-) const
-{
-    const fvMesh& mesh = this->owner().mesh();
-
-    const volVectorField& Ufield =
-        mesh.objectRegistry::lookupObject<volVectorField>(UName_);
-
-    label patchI = pp.index();
-    label patchFaceI = pp.whichFace(p.face());
-
-    vector n = tetIs.faceTri(mesh).normal();
-    n /= mag(n);
-
-    vector U = Ufield.boundaryField()[patchI][patchFaceI];
-
-    // Unless the face is rotating, the required normal is n;
-    nw = n;
-
-    if (!mesh.moving())
-    {
-        // Only wall patches may have a non-zero wall velocity from
-        // the velocity field when the mesh is not moving.
-
-        if (isA<wallPolyPatch>(pp))
-        {
-            Up = U;
-        }
-        else
-        {
-            Up = vector::zero;
-        }
-    }
-    else
-    {
-        vector U00 = Ufield.oldTime().boundaryField()[patchI][patchFaceI];
-
-        vector n00 = tetIs.oldFaceTri(mesh).normal();
-
-        // Difference in normal over timestep
-        vector dn = vector::zero;
-
-        if (mag(n00) > SMALL)
-        {
-            // If the old normal is zero (for example in layer
-            // addition) then use the current normal, meaning that the
-            // motion can only be translational, and dn remains zero,
-            // otherwise, calculate dn:
-
-            n00 /= mag(n00);
-
-            dn = n - n00;
-        }
-
-        // Total fraction thought the timestep of the motion,
-        // including stepFraction before the current tracking step
-        // and the current trackFraction
-        // i.e.
-        // let s = stepFraction, t = trackFraction
-        // Motion of x in time:
-        // |-----------------|---------|---------|
-        // x00               x0        xi        x
-        //
-        // where xi is the correct value of x at the required
-        // tracking instant.
-        //
-        // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
-        //
-        // i.e. the motion covered by previous tracking portions
-        // within this timestep, and
-        //
-        // xi = x0 + t*(x - x0)
-        //    = t*x + (1 - t)*x0
-        //    = t*x + (1 - t)*(s*x + (1 - s)*x00)
-        //    = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
-        //
-        // let m = (s + t - s*t)
-        //
-        // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
-        //
-        // In the same form as before.
-
-        scalar m =
-            p.stepFraction()
-          + trackFraction
-          - (p.stepFraction()*trackFraction);
-
-        // When the mesh is moving, the velocity field on wall patches
-        // will contain the velocity associated with the motion of the
-        // mesh, in which case it is interpolated in time using m.
-        // For other patches the face velocity will need to be
-        // reconstructed from the face centre motion.
-
-        const vector& Cf = mesh.faceCentres()[p.face()];
-
-        vector Cf00 = mesh.faces()[p.face()].centre(mesh.oldPoints());
-
-        if (isA<wallPolyPatch>(pp))
-        {
-            Up = U00 + m*(U - U00);
-        }
-        else
-        {
-            Up = (Cf - Cf00)/this->owner().time().deltaTValue();
-        }
-
-        if (mag(dn) > SMALL)
-        {
-            // Rotational motion, nw requires interpolation and a
-            // rotational velocity around face centre correction to Up
-            // is required.
-
-            nw = n00 + m*dn;
-
-            // Cf at tracking instant
-            vector Cfi = Cf00 + m*(Cf - Cf00);
-
-            // Normal vector cross product
-            vector omega = (n00 ^ n);
-
-            scalar magOmega = mag(omega);
-
-            // magOmega = sin(angle between unit normals)
-            // Normalise omega vector by magOmega, then multiply by
-            // angle/dt to give the correct angular velocity vector.
-            omega *=
-                Foam::asin(magOmega)
-               /(magOmega*this->owner().time().deltaTValue());
-
-            // Project position onto face and calculate this position
-            // relative to the face centre.
-            vector facePos =
-                p.position()
-              - ((p.position() - Cfi) & nw)*nw
-              - Cfi;
-
-            Up += (omega ^ facePos);
-        }
-
-        // No further action is required if the motion is
-        // translational only, nw and Up have already been set.
-    }
-}
-
-
 template<class CloudType>
 void Foam::PatchInteractionModel<CloudType>::info(Ostream& os)
 {

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -164,18 +164,6 @@ public:
             const tetIndices& tetIs
         );
 
-        //- Calculate the patch normal and velocity to interact with,
-        //  accounting for patch motion if required.
-        void patchData
-        (
-            typename CloudType::parcelType& p,
-            const polyPatch& pp,
-            const scalar trackFraction,
-            const tetIndices& tetIs,
-            vector& normal,
-            vector& Up
-        ) const;
-
 
         // I-O
 

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,7 +74,7 @@ bool Foam::Rebound<CloudType>::correct
     vector nw;
     vector Up;
 
-    this->patchData(p, pp, trackFraction, tetIs, nw, Up);
+    this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
     // Calculate motion relative to patch velocity
     U -= Up;

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C

@@ -148,7 +148,7 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 vector nw;
                 vector Up;
 
-                this->patchData(p, pp, trackFraction, tetIs, nw, Up);
+                this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
 
                 // Calculate motion relative to patch velocity
                 U -= Up;

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C

@@ -156,18 +156,20 @@ Foam::CompositionModel<CloudType>::componentNames(const label phaseI) const
 template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::globalCarrierId
 (
-    const word& cmptName
+    const word& cmptName,
+    const bool allowNotFound
 ) const
 {
     label id = thermo_.carrierId(cmptName);
 
-    if (id < 0)
+    if (id < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
             "Foam::label Foam::CompositionModel<CloudType>::globalCarrierId"
             "("
-                "const word&"
+                "const word&, "
+                "const bool"
             ") const"
         )   << "Unable to determine global id for requested component "
             << cmptName << ". Available components are " << nl
@@ -182,19 +184,21 @@ template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::globalId
 (
     const label phaseI,
-    const word& cmptName
+    const word& cmptName,
+    const bool allowNotFound
 ) const
 {
     label id = phaseProps_[phaseI].globalId(cmptName);
 
-    if (id < 0)
+    if (id < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
             "Foam::label Foam::CompositionModel<CloudType>::globalId"
             "("
                 "const label, "
-                "const word&"
+                "const word&, "
+                "const bool"
             ") const"
         )   << "Unable to determine global id for requested component "
             << cmptName << abort(FatalError);
@@ -218,19 +222,21 @@ template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::localId
 (
     const label phaseI,
-    const word& cmptName
+    const word& cmptName,
+    const bool allowNotFound
 ) const
 {
     label id = phaseProps_[phaseI].id(cmptName);
 
-    if (id < 0)
+    if (id < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
             "Foam::label Foam::CompositionModel<CloudType>::localId"
             "("
                 "const label, "
-                "const word&"
+                "const word&, "
+                "const bool"
             ") const"
         )   << "Unable to determine local id for component " << cmptName
             << abort(FatalError);
@@ -244,12 +250,13 @@ template<class CloudType>
 Foam::label Foam::CompositionModel<CloudType>::localToGlobalCarrierId
 (
     const label phaseI,
-    const label id
+    const label id,
+    const bool allowNotFound
 ) const
 {
     label gid = phaseProps_[phaseI].globalCarrierIds()[id];
 
-    if (gid < 0)
+    if (gid < 0 && !allowNotFound)
     {
         FatalErrorIn
         (
@@ -257,7 +264,8 @@ Foam::label Foam::CompositionModel<CloudType>::localToGlobalCarrierId
             "Foam::CompositionModel<CloudType>::localToGlobalCarrierId"
             "("
                 "const label, "
-                "const label"
+                "const label, "
+                "const bool"
             ") const"
         )   << "Unable to determine global carrier id for phase "
             << phaseI << " with local id " << id

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -166,22 +166,37 @@ public:
                 const wordList& componentNames(const label phaseI) const;
 
                 //- Return global id of component cmptName in carrier thermo
-                label globalCarrierId(const word& cmptName) const;
+                label globalCarrierId
+                (
+                    const word& cmptName,
+                    const bool allowNotFound = false
+                ) const;
 
                 //- Return global id of component cmptName in phase phaseI
-                label globalId(const label phaseI, const word& cmptName) const;
+                label globalId
+                (
+                    const label phaseI,
+                    const word& cmptName,
+                    const bool allowNotFound = false
+                ) const;
 
                 //- Return global ids of for phase phaseI
                 const labelList& globalIds(const label phaseI) const;
 
                 //- Return local id of component cmptName in phase phaseI
-                label localId(const label phaseI, const word& cmptName) const;
+                label localId
+                (
+                    const label phaseI,
+                    const word& cmptName,
+                    const bool allowNotFound = false
+                ) const;
 
                 //- Return global carrier id of component given local id
                 label localToGlobalCarrierId
                 (
                     const label phaseI,
-                    const label id
+                    const label id,
+                    const bool allowNotFound = false
                 ) const;
 
                 //- Return the list of phase phaseI mass fractions

src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C

@@ -64,16 +64,7 @@ Foam::ReactingLookupTableInjection<CloudType>::ReactingLookupTableInjection
     injectorTetFaces_.setSize(injectors_.size());
     injectorTetPts_.setSize(injectors_.size());
 
-    forAll(injectors_, i)
+    updateMesh();
-    {
-        this->findCellAtPosition
-        (
-            injectorCells_[i],
-            injectorTetFaces_[i],
-            injectorTetPts_[i],
-            injectors_[i].x()
-        );
-    }
 
     // Determine volume of particles to inject
     this->volumeTotal_ = 0.0;
@@ -111,6 +102,23 @@ Foam::ReactingLookupTableInjection<CloudType>::~ReactingLookupTableInjection()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ReactingLookupTableInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    forAll(injectors_, i)
+    {
+        this->findCellAtPosition
+        (
+            injectorCells_[i],
+            injectorTetFaces_[i],
+            injectorTetPts_[i],
+            injectors_[i].x()
+        );
+    }
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::timeEnd() const
 {

src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H

@@ -132,6 +132,9 @@ public:
 
     // Member Functions
 
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
         //- Return the end-of-injection time
         scalar timeEnd() const;
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,10 +102,13 @@ template<class CloudType>
 void Foam::ConstantRateDevolatilisation<CloudType>::calculate
 (
     const scalar dt,
+    const scalar age,
     const scalar mass0,
     const scalar mass,
     const scalar T,
     const scalarField& YGasEff,
+    const scalarField& YLiquidEff,
+    const scalarField& YSolidEff,
     bool& canCombust,
     scalarField& dMassDV
 ) const

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/ConstantRateDevolatilisation/ConstantRateDevolatilisation.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,10 +104,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.C

@@ -78,6 +78,9 @@ void Foam::DevolatilisationModel<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
+    const scalar,
+    const scalarField&,
+    const scalarField&,
     const scalarField&,
     bool&,
     scalarField&
@@ -91,6 +94,9 @@ void Foam::DevolatilisationModel<CloudType>::calculate
             "const scalar, "
             "const scalar, "
             "const scalar, "
+            "const scalar, "
+            "const scalarField&, "
+            "const scalarField&, "
             "const scalarField&, "
             "bool&, "
             "scalarField&"

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/DevolatilisationModel/DevolatilisationModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class DevolatilisationModel Declaration
+                    Class DevolatilisationModel Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -126,10 +126,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,9 @@ void Foam::NoDevolatilisation<CloudType>::calculate
     const scalar,
     const scalar,
     const scalar,
+    const scalar,
+    const scalarField&,
+    const scalarField&,
     const scalarField&,
     bool& canCombust,
     scalarField&

src/lagrangian/intermediate/submodels/ReactingMultiphase/DevolatilisationModel/NoDevolatilisation/NoDevolatilisation.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,10 +85,13 @@ public:
         virtual void calculate
         (
             const scalar dt,
+            const scalar age,
             const scalar mass0,
             const scalar mass,
             const scalar T,
             const scalarField& YGasEff,
+            const scalarField& YLiquidEff,
+            const scalarField& YSolidEff,
             bool& canCombust,
             scalarField& dMassDV
         ) const;