From 5160b81e74ad22334674ac2a71c80c41fadc67c9 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 9 Jan 2009 16:25:17 +0000
Subject: [PATCH] restructuring so that injection is handled by injection model
 and not cloud

---
 .../Templates/KinematicCloud/KinematicCloud.C | 249 ++------------
 .../Templates/KinematicCloud/KinematicCloud.H |  95 +-----
 .../KinematicCloud/KinematicCloudI.H          |  35 --
 .../Templates/ReactingCloud/ReactingCloud.C   | 162 +++-------
 .../Templates/ReactingCloud/ReactingCloud.H   |  17 +-
 .../Templates/ThermoCloud/ThermoCloud.C       |  51 ++-
 .../Templates/ThermoCloud/ThermoCloud.H       |  20 +-
 .../InjectionModel/InjectionModel.C           | 304 ++++++++++++------
 .../InjectionModel/InjectionModel.H           | 158 ++++++---
 9 files changed, 470 insertions(+), 621 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 1a59476437..8ed5e0c262 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -36,43 +36,32 @@ License
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>
-Foam::scalar Foam::KinematicCloud<ParcelType>::setNumberOfParticles
+void Foam::KinematicCloud<ParcelType>::addNewParcel
 (
-    const label nParcels,
-    const scalar pDiameter,
-    const scalar pVolumeFraction,
-    const scalar pRho,
-    const scalar pVolume
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
 )
 {
-    scalar nP = 0.0;
-    switch (parcelBasis_)
-    {
-        case pbMass:
-        {
-            nP = pVolumeFraction*massTotal_/nParcels
-               /(pRho*mathematicalConstant::pi/6.0*pow(pDiameter, 3));
-            break;
-        }
-        case pbNumber:
-        {
-            nP = pVolumeFraction*massTotal_/(pRho*pVolume);
-            break;
-        }
-        default:
-        {
-            nP = 0.0;
-            FatalErrorIn
-            (
-                "Foam::KinematicCloud<ParcelType>::setNumberOfParticles"
-                "(const label, const scalar, const scalar, const scalar, "
-                "const scalar)"
-            )<< "Unknown parcelBasis type" << nl
-             << exit(FatalError);
-        }
-    }
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        parcelTypeId_,
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        constProps_
+    );
 
-    return nP;
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
 }
 
 
@@ -107,14 +96,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
     coupled_(particleProperties_.lookup("coupled")),
     rndGen_(label(0)),
-    time0_(this->db().time().value()),
-    parcelBasisType_(particleProperties_.lookup("parcelBasisType")),
-    parcelBasis_(pbNumber),
-    massTotal_
-    (
-        dimensionedScalar(particleProperties_.lookup("massTotal")).value()
-    ),
-    massInjected_(0.0),
     rho_(rho),
     U_(U),
     mu_(mu),
@@ -160,9 +141,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             particleProperties_.subDict("integrationSchemes")
         )
     ),
-    nInjections_(0),
-    nParcelsAdded_(0),
-    nParcelsAddedTotal_(0),
     UTrans_
     (
         IOobject
@@ -191,27 +169,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
         mesh_,
         dimensionedScalar("zero",  dimensionSet(1, 0, -1, 0, 0), 0.0)
     )
-{
-    if (parcelBasisType_ == "mass")
-    {
-        parcelBasis_ = pbMass;
-    }
-    else if (parcelBasisType_ == "number")
-    {
-        parcelBasis_ = pbNumber;
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "Foam::KinematicCloud<ParcelType>::KinematicCloud"
-            "(const word&, const volScalarField&"
-            ", const volVectorField&, const volScalarField&, const "
-            "dimensionedVector&)"
-        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
-         << exit(FatalError);
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -265,7 +223,12 @@ void Foam::KinematicCloud<ParcelType>::evolve()
         g_.value()
     );
 
-    inject();
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (coupled_)
     {
@@ -276,160 +239,16 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::inject()
-{
-    scalar time = this->db().time().value();
-
-    scalar pRho = constProps_.rho0();
-
-    this->injection().prepareForNextTimeStep(time0_, time);
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcels();
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume();
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction = this->injection().volumeFraction();
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->db().time().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - time0_
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - time0_
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->injection().findInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                *this,
-                parcelTypeId_,
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                constProps_
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->time().deltaT().value();
-
-            this->injectParcel(pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-}
-
-
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::injectParcel(ParcelType* p)
-{
-    addParticle(p);
-    nParcelsAdded_++;
-    nParcelsAddedTotal_++;
-    massInjected_ += p->mass()*p->nParticle();
-}
-
-
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::postInjectCheck()
-{
-    if (nParcelsAdded_)
-    {
-        Pout<< "\n--> Cloud: " << this->name() << nl
-            << "    Added " << nParcelsAdded_
-            <<  " new parcels" << nl << endl;
-    }
-
-    // Reset parcel counters
-    nParcelsAdded_ = 0;
-
-    // Set time for start of next injection
-    time0_ = this->db().time().value();
-
-    // Increment number of injections
-    nInjections_++;
-}
-
-
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::info() const
 {
     Info<< "Cloud name: " << this->name() << nl
         << "    Parcels added during this run   = "
-        << returnReduce(nParcelsAddedTotal_, sumOp<label>()) << nl
+        << returnReduce(this->injection().nParcelsAddedTotal(), sumOp<label>())
+            << nl
         << "    Mass introduced during this run = "
-        << returnReduce(massInjected_, sumOp<scalar>()) << nl
+        << returnReduce(this->injection().massInjected(), sumOp<scalar>())
+            << nl
         << "    Current number of parcels       = "
         << returnReduce(this->size(), sumOp<label>()) << nl
         << "    Current mass in system          = "
@@ -445,7 +264,7 @@ void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
     (
         this->db().time().path()/"parcelPositions_"
       + this->name() + "_"
-      + name(this->nInjections_) + ".obj"
+      + name(this->injection().nInjections()) + ".obj"
     );
 
     forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 0d3781413c..556dbc3546 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -85,19 +85,6 @@ class KinematicCloud
     public kinematicCloud
 {
 
-public:
-
-    // Enumerations
-
-        //- Parcel basis representation options
-        //  i.e constant number of particles OR constant mass per parcel
-        enum parcelBasis
-        {
-            pbNumber,
-            pbMass
-        };
-
-
 private:
 
     // Private data
@@ -126,22 +113,6 @@ private:
         //- Random number generator - used by some injection routines
         Random rndGen_;
 
-        //- Time at beginning of timestep
-        scalar time0_;
-
-
-        // Injection properties
-
-            //- Parcel basis
-            const word parcelBasisType_;
-            parcelBasis parcelBasis_;
-
-            //- Total mass to inject [kg]
-            scalar massTotal_;
-
-            //- Total mass injected to date [kg]
-            scalar massInjected_;
-
 
         // References to the carrier gas fields
 
@@ -161,9 +132,8 @@ private:
             const dimensionedVector& g_;
 
 
-        // Interpolation
-
-            dictionary interpolationSchemes_;
+        //- Interpolation schemes dictionary
+        dictionary interpolationSchemes_;
 
 
         // References to the cloud sub-models
@@ -190,17 +160,6 @@ private:
             autoPtr<vectorIntegrationScheme> UIntegrator_;
 
 
-        // Counters
-
-            //- Number of injections counter
-            label nInjections_;
-
-            //- Running counters describing parcels added during each
-            //  injection
-            label nParcelsAdded_;
-            label nParcelsAddedTotal_;
-
-
         // Sources
 
             //- Momentum
@@ -219,30 +178,6 @@ private:
         void operator=(const KinematicCloud&);
 
 
-protected:
-
-    // Protected member functions
-
-        //- Set the number of particles per parcel
-        scalar setNumberOfParticles
-        (
-            const label nParcels,
-            const scalar pDiameter,
-            const scalar pVolumeFraction,
-            const scalar pRho,
-            const scalar pVolume
-        );
-
-        //- Inject more parcels
-        void inject();
-
-        //- Inject parcel if it is valid - delete otherwise
-        void injectParcel(ParcelType* p);
-
-        //- Post-injection checks
-        void postInjectCheck();
-
-
 public:
 
     // Constructors
@@ -286,12 +221,6 @@ public:
             //- Return refernce to the random object
             inline Random& rndGen();
 
-            //- Return the start of injection interval time
-            inline scalar time0() const;
-
-            //- Return a reference to the mass of particles to introduce
-            inline scalar massTotal() const;
-
 
             // References to the carrier gas fields
 
@@ -380,15 +309,6 @@ public:
             void dumpParticlePositions() const;
 
 
-            // Counters
-
-                //- Return the number of injections
-                inline label nInjections() const;
-
-                //- Return the total number parcels added
-                inline label nParcelsAddedTotal() const;
-
-
             // Fields
 
                 //- Return the particle volume fraction field
@@ -402,6 +322,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 73a891dcb3..7124410731 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -155,41 +155,6 @@ Foam::KinematicCloud<ParcelType>::UIntegrator() const
 }
 
 
-template<class ParcelType>
-inline Foam::label Foam::KinematicCloud<ParcelType>::nInjections() const
-{
-    return nInjections_;
-}
-
-
-template<class ParcelType>
-inline Foam::label Foam::KinematicCloud<ParcelType>::nParcelsAddedTotal() const
-{
-    return nParcelsAddedTotal_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::time0() const
-{
-    return time0_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::massTotal() const
-{
-    return massTotal_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::massInjected() const
-{
-    return massInjected_;
-}
-
-
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicCloud<ParcelType>::massInSystem() const
 {
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
index 4bad764272..95f2e9810e 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C
@@ -29,6 +29,42 @@ License
 #include "MassTransferModel.H"
 #include "SurfaceReactionModel.H"
 
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::ReactingCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        this->parcelTypeId(),
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        composition().YGas0(),
+        composition().YLiquid0(),
+        composition().YSolid0(),
+        composition().YMixture0(),
+        constProps_
+    );
+
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -174,7 +210,12 @@ void Foam::ReactingCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject();
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (this->coupled())
     {
@@ -185,123 +226,4 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-void Foam::ReactingCloud<ParcelType>::inject()
-{
-    scalar time = this->db().time().value();
-
-    scalar pRho = this->constProps().rho0();
-
-    this->injection().prepareForNextTimeStep(this->time0(), time);
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcels();
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume();
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction = this->injection().volumeFraction();
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->db().time().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - this->time0()
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - this->time0()
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = this->time0() + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = this->setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->injection().findInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                *this,
-                this->parcelTypeId(),
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                composition().YGas0(),
-                composition().YLiquid0(),
-                composition().YSolid0(),
-                composition().YMixture0(),
-                this->constProps()
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->db().time().deltaT().value();
-
-            this->injectParcel(pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
index 9ce5c98b1a..e257d04484 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H
@@ -114,12 +114,6 @@ class ReactingCloud
         void operator=(const ReactingCloud&);
 
 
-protected:
-
-    //- Inject more parcels
-    void inject();
-
-
 public:
 
     //- Runtime type information
@@ -199,6 +193,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
index c585d22fd3..184d6d8aae 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@@ -30,6 +30,38 @@ License
 #include "interpolationCellPoint.H"
 #include "ThermoParcel.H"
 
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::ThermoCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        this->parcelTypeId(),
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        constProps_
+    );
+
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -167,7 +199,12 @@ void Foam::ThermoCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject(td);
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (this->coupled())
     {
@@ -178,16 +215,4 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoCloud<ParcelType>::inject
-(
-    TrackingData& td
-)
-{
-    // Injection is same as for KinematicCloud<ParcelType>
-    KinematicCloud<ParcelType>::inject(td);
-}
-
-
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
index ba7930af30..a3a0044c47 100644
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@@ -116,13 +116,6 @@ class ThermoCloud
         void operator=(const ThermoCloud&);
 
 
-protected:
-
-    //- Inject more parcels
-    template<class TrackingData>
-    void inject(TrackingData& td);
-
-
 public:
 
     //- Runtime type information
@@ -173,7 +166,7 @@ public:
 
             // Modelling options
 
-                //- Radiation flag
+                 //- Radiation flag
                 inline bool radiation() const;
 
 
@@ -208,6 +201,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index f23dc19c0a..f3c1caaa31 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -25,109 +25,22 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "InjectionModel.H"
+#include "mathematicalConstants.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InjectionModel<CloudType>::InjectionModel
-(
-    const dictionary& dict,
-    CloudType& owner,
-    const word& type
-)
-:   dict_(dict),
-    owner_(owner),
-    coeffDict_(dict.subDict(type + "Coeffs")),
-    SOI_(readScalar(coeffDict_.lookup("SOI"))),
-    volumeTotal_(0.0),
-    timeStep0_(0.0),
-    nParcels_(0),
-    volume_(0.0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InjectionModel<CloudType>::~InjectionModel()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class CloudType>
-const CloudType& Foam::InjectionModel<CloudType>::owner() const
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-CloudType& Foam::InjectionModel<CloudType>::owner()
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::InjectionModel<CloudType>::dict() const
-{
-    return dict_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::InjectionModel<CloudType>::coeffDict() const
-{
-    return coeffDict_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::timeStart() const
-{
-    return SOI_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeTotal() const
-{
-    return volumeTotal_;
-}
-
-
-template<class CloudType>
-Foam::label Foam::InjectionModel<CloudType>::nParcels() const
-{
-    return nParcels_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volume() const
-{
-    return volume_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeFraction() const
-{
-    return volume_/volumeTotal_;
-}
-
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class CloudType>
 void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 (
     const scalar time0,
-    const scalar time1
+    const scalar time1,
+    label& nParcels,
+    scalar& volume
 )
 {
     // Initialise values
-    nParcels_ = 0;
-    volume_ = 0.0;
+    nParcels = 0;
+    volume = 0.0;
 
     // Return if not started injection event
     if (time1 < SOI_)
@@ -141,13 +54,13 @@ void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
     scalar t1 = time1 - SOI_;
 
     // Number of parcels to inject
-    nParcels_ = nParcelsToInject(t0, t1);
+    nParcels = nParcelsToInject(t0, t1);
 
     // Volume of parcels to inject
-    volume_ = volumeToInject(t0, t1);
+    volume = volumeToInject(t0, t1);
 
     // Hold previous time if no parcels, but non-zero volume fraction
-    if ((nParcels_ == 0) && (volume_ > 0.0))
+    if ((nParcels == 0) && (volume > 0.0))
     {
         // hold value of timeStep0_
     }
@@ -210,6 +123,203 @@ void Foam::InjectionModel<CloudType>::findInjectorCellAndPosition
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
+(
+    const label nParcels,
+    const scalar diameter,
+    const scalar volumeFraction,
+    const scalar rho,
+    const scalar volume
+)
+{
+    scalar nP = 0.0;
+    switch (parcelBasis_)
+    {
+        case pbMass:
+        {
+            nP = volumeFraction*massTotal_/nParcels
+               /(rho*mathematicalConstant::pi/6.0*pow3(diameter));
+            break;
+        }
+        case pbNumber:
+        {
+            nP = volumeFraction*massTotal_/(rho*volume);
+            break;
+        }
+        default:
+        {
+            nP = 0.0;
+            FatalErrorIn
+            (
+                "void Foam::InjectionModel<CloudType>::setNumberOfParticles"
+                "(const label, const scalar, const scalar, const scalar, "
+                "const scalar)"
+            )<< "Unknown parcelBasis type" << nl
+             << exit(FatalError);
+        }
+    }
+
+    return nP;
+}
+
+
+template<class CloudType>
+void Foam::InjectionModel<CloudType>::postInjectCheck()
+{
+    if (nParcelsAdded_ > 0)
+    {
+        Pout<< "\n--> Cloud: " << owner_.name() << nl
+            << "    Added " << nParcelsAdded_
+            <<  " new parcels" << nl << endl;
+    }
+
+    // Increment total number of parcels added
+    nParcelsAddedTotal_ += nParcelsAdded_;
+
+    // Reset parcel counters
+    nParcelsAdded_ = 0;
+
+    // Update time for start of next injection
+    time0_ = owner_.db().time().value();
+
+    // Increment number of injections
+    nInjections_++;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectionModel<CloudType>::InjectionModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:   dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs")),
+    SOI_(readScalar(coeffDict_.lookup("SOI"))),
+    volumeTotal_(0.0),
+    massTotal_(dimensionedScalar(coeffDict_.lookup("massTotal")).value()),
+    massInjected_(0.0),
+    nInjections_(0),
+    nParcelsAdded_(0),
+    nParcelsAddedTotal_(0),
+    parcelBasisType_(coeffDict_.lookup("parcelBasisType")),
+    parcelBasis_(pbNumber),
+    time0_(owner.db().time().value()),
+    timeStep0_(0.0)
+{
+    if (parcelBasisType_ == "mass")
+    {
+        parcelBasis_ = pbMass;
+    }
+    else if (parcelBasisType_ == "number")
+    {
+        parcelBasis_ = pbNumber;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::InjectionModel<CloudType>::InjectionModel"
+            "(const dictionary&, CloudType&, const word&)"
+        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
+         << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectionModel<CloudType>::~InjectionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+template<class TrackData>
+void Foam::InjectionModel<CloudType>::inject(TrackData& td)
+{
+    const scalar time = owner_.db().time().value();
+    const scalar continuousDt = owner_.db().time().deltaT().value();
+
+    // Prepare for next time step
+    nParcelsAdded_ = 0;
+    label nParcels = 0;
+    scalar volume = 0.0;
+    prepareForNextTimeStep(time0_, time, nParcels, volume);
+
+    // Return if no parcels are required
+    if (nParcels == 0)
+    {
+        postInjectCheck();
+        return;
+    }
+
+    // Particle density given by constant properties
+    const scalar rho = td.constProps().rho0();
+
+    // Volume fraction to introduce during this timestep
+    const scalar volFraction = volumeFraction(volume);
+
+    // Duration of injection period during this timestep
+    const scalar deltaT = min
+    (
+        continuousDt,
+        min(time - SOI_, timeEnd() - time0_)
+    );
+
+    // Pad injection time if injection starts during this timestep
+    const scalar padTime = max(0.0, SOI_ - time0_);
+
+    // Introduce new parcels linearly with time
+    for (label iParcel=0; iParcel<nParcels; iParcel++)
+    {
+        // Calculate the pseudo time of injection for parcel 'iParcel'
+        scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
+
+        // Determine injected parcel properties
+        vector pos = position(iParcel, timeInj, owner_.meshInfo());
+
+        // Diameter of parcels
+        scalar d = d0(iParcel, timeInj);
+
+        // Number of particles per parcel
+        scalar nP = setNumberOfParticles
+        (
+            nParcels,
+            d,
+            volFraction,
+            rho,
+            volume
+        );
+
+        // Velocity of parcels
+        vector U = velocity(iParcel, timeInj, owner_.meshInfo());
+
+        // Determine the injection cell
+        label cellI = -1;
+        findInjectorCellAndPosition(cellI, pos);
+
+        if (cellI >= 0)
+        {
+            scalar dt = time - timeInj;
+            td.cloud().addNewParcel(pos, cellI, d, U, nP, dt);
+
+            massInjected_ += nP*rho*mathematicalConstant::pi*pow3(d)/6.0;
+            nParcelsAdded_++;
+        }
+    }
+
+    postInjectCheck();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewInjectionModel.C"
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
index 0b00ebeda9..0df2a5ef6d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H
@@ -55,6 +55,21 @@ template<class CloudType>
 class InjectionModel
 {
 
+public:
+
+    // Enumerations
+
+        //- Parcel basis representation options
+        //  i.e constant number of particles OR constant mass per parcel
+        enum parcelBasis
+        {
+            pbNumber,
+            pbMass
+        };
+
+
+private:
+
     // Private data
 
         //- The cloud dictionary
@@ -74,24 +89,47 @@ protected:
         // Global injection properties
 
             //- Start of injection [s]
-            scalar SOI_;
+            const scalar SOI_;
 
             //- Total volume of parcels to introduce [m^3]
             //  Initialised in the individual injection models
             scalar volumeTotal_;
 
+            //- Total mass to inject [kg]
+            const scalar massTotal_;
+
+            //- Total mass injected to date [kg]
+            scalar massInjected_;
+
+
+        // Counters
+
+            //- Number of injections counter
+            label nInjections_;
+
+            //- Running counter of parcels added during each injection
+            label nParcelsAdded_;
+
+            //- Running counter of total number of parcels added
+            label nParcelsAddedTotal_;
+
 
         // Injection properties per Lagrangian time step
 
+            // Parcel basis
+
+                //- Parcel basis name
+                const word parcelBasisType_;
+
+                //- Parcel basis enumeration
+                parcelBasis parcelBasis_;
+
+            //- Continuous phase time at start of injection time step [s]
+            scalar time0_;
+
             //- Time at start of injection time step [s]
             scalar timeStep0_;
 
-            //- Number of parcels to introduce []
-            label nParcels_;
-
-            //- Volume of parcels to introduce [m^3]
-            scalar volume_;
-
 
     // Protected member functions
 
@@ -110,6 +148,37 @@ protected:
         ) const = 0;
 
 
+        //- Determine properties for next time step/injection interval
+        void prepareForNextTimeStep
+        (
+            const scalar time0,
+            const scalar time1,
+            label& nParcels,
+            scalar& volume
+        );
+
+        //- Find the cell that contains the injector position
+        //  Will modify position slightly towards the owner cell centroid
+        virtual void findInjectorCellAndPosition
+        (
+            label& cellI,
+            vector& position
+        );
+
+        //- Set number of particles to inject given parcel properties
+        scalar setNumberOfParticles
+        (
+            const label nParcels,
+            const scalar diameter,
+            const scalar volumeFraction,
+            const scalar rho,
+            const scalar volume
+        );
+
+        //- Post injection checks
+        void postInjectCheck();
+
+
 public:
 
     //- Runtime type information
@@ -140,33 +209,31 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~InjectionModel();
+    //- Destructor
+    virtual ~InjectionModel();
 
 
-    // Selector
-
-        static autoPtr<InjectionModel<CloudType> > New
-        (
-            const dictionary& dict,
-            CloudType& owner
-        );
+    //- Selector
+    static autoPtr<InjectionModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
 
 
     // Access
 
+        //- Return the owner cloud dictionary
+        inline const dictionary& dict() const;
+
         //- Return const access the owner cloud object
-        const CloudType& owner() const;
+        inline const CloudType& owner() const;
 
         //- Return non-const access the owner cloud object for manipulation
-        CloudType& owner();
-
-        //- Return the dictionary
-        const dictionary& dict() const;
+        inline CloudType& owner();
 
         //- Return the coefficients dictionary
-        const dictionary& coeffDict() const;
+        inline const dictionary& coeffDict() const;
 
 
     // Member Functions
@@ -178,44 +245,41 @@ public:
         // Global information
 
             //- Return the start-of-injection time
-            scalar timeStart() const;
+            inline scalar timeStart() const;
 
             //- Return the total volume to be injected across the event
-            scalar volumeTotal() const;
+            inline scalar volumeTotal() const;
+
+            //- Return mass of particles to introduce
+            inline scalar massTotal() const;
+
+            //- Return mass of particles injected (cummulative)
+            inline scalar massInjected() const;
 
             //- Return the end-of-injection time
             virtual scalar timeEnd() const = 0;
 
+            // Counters
 
-        // Per Lagrangian time step properties
+                //- Return the number of injections
+                inline label nInjections() const;
 
-            //- Determine properties for next time step/injection interval
-            void prepareForNextTimeStep
-            (
-                const scalar time0,
-                const scalar time1
-            );
+                //- Return the total number parcels added
+                inline label nParcelsAddedTotal() const;
 
-            //- Return the number of parcels to introduce
-            label nParcels() const;
 
-            //- Return the volume of parcels to introduce
-            scalar volume() const;
+        // Per-injection event functions
+
+            //- Main injection loop
+            template<class TrackData>
+            void inject(TrackData& td);
 
             //- Return the volume fraction to introduce
-            scalar volumeFraction() const;
+            inline scalar volumeFraction(const scalar volume) const;
 
 
         // Injection geometry
 
-            //- Find the cell that contains the injector position
-            //  Will modify position slightly towards the owner cell centroid
-            virtual void findInjectorCellAndPosition
-            (
-                label& cellI,
-                vector& position
-            );
-
             //- Return the injection position
             virtual vector position
             (
@@ -241,6 +305,10 @@ public:
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "InjectionModelI.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam