Merge branch 'feature-extract-eulerian-particles' into 'develop'

New feature extract eulerian particles

New functionality to extract particle data from multiphase calculations and replay the data in lagrangian cases, using both the raw input particle data, and data processed into a (smaller) set of injection locations.

See merge request !82

applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,7 +69,6 @@ int main(int argc, char *argv[])
 
     bool writeAgglom = readBool(agglomDict.lookup("writeFacesAgglomeration"));
 
-
     const polyBoundaryMesh& boundary = mesh.boundaryMesh();
 
     labelListIOList finalAgglom
@@ -95,15 +94,20 @@ int main(int argc, char *argv[])
         {
             label patchi =  patchids[i];
             const polyPatch& pp = boundary[patchi];
+
             if (!pp.coupled())
             {
                 Info << "\nAgglomerating patch : " << pp.name() << endl;
+
                 pairPatchAgglomeration agglomObject
                 (
-                    pp,
+                    pp.localFaces(),
+                    pp.localPoints(),
                     agglomDict.subDict(pp.name())
                 );
+
                 agglomObject.agglomerate();
+
                 finalAgglom[patchi] =
                     agglomObject.restrictTopBottomAddressing();
 
@@ -116,42 +120,42 @@ int main(int argc, char *argv[])
     }
 
 
-    // - All patches which are not agglomarated are identity for finalAgglom
-    forAll(boundary, patchid)
+    // All patches which are not agglomerated are identity for finalAgglom
+    forAll(boundary, patchi)
     {
-        if (finalAgglom[patchid].size() == 0)
+        if (finalAgglom[patchi].size() == 0)
         {
-            finalAgglom[patchid] = identity(boundary[patchid].size());
+            finalAgglom[patchi] = identity(boundary[patchi].size());
         }
     }
 
     // Sync agglomeration across coupled patches
     labelList nbrAgglom(mesh.nFaces() - mesh.nInternalFaces(), -1);
 
-    forAll(boundary, patchid)
+    forAll(boundary, patchi)
     {
-        const polyPatch& pp = boundary[patchid];
+        const polyPatch& pp = boundary[patchi];
         if (pp.coupled())
         {
-            finalAgglom[patchid] = identity(pp.size());
+            finalAgglom[patchi] = identity(pp.size());
             forAll(pp, i)
             {
-                nbrAgglom[pp.start() - mesh.nInternalFaces() + i] =
-                    finalAgglom[patchid][i];
+                const label agglomi = pp.start() - mesh.nInternalFaces() + i;
+                nbrAgglom[agglomi] = finalAgglom[patchi][i];
             }
         }
     }
 
     syncTools::swapBoundaryFaceList(mesh, nbrAgglom);
-    forAll(boundary, patchid)
+    forAll(boundary, patchi)
     {
-        const polyPatch& pp = boundary[patchid];
+        const polyPatch& pp = boundary[patchi];
         if (pp.coupled() && !refCast<const coupledPolyPatch>(pp).owner())
         {
             forAll(pp, i)
             {
-                finalAgglom[patchid][i] =
-                    nbrAgglom[pp.start() - mesh.nInternalFaces() + i];
+                const label agglomi = pp.start() - mesh.nInternalFaces() + i;
+                finalAgglom[patchi][i] = nbrAgglom[agglomi];
             }
         }
     }
@@ -179,13 +183,13 @@ int main(int argc, char *argv[])
             facesAgglomeration.boundaryFieldRef();
 
         label coarsePatchIndex = 0;
-        forAll(boundary, patchid)
+        forAll(boundary, patchi)
         {
-            const polyPatch& pp = boundary[patchid];
+            const polyPatch& pp = boundary[patchi];
             if (pp.size() > 0)
             {
                 fvPatchScalarField& bFacesAgglomeration =
-                    facesAgglomerationBf[patchid];
+                    facesAgglomerationBf[patchi];
 
                 forAll(bFacesAgglomeration, j)
                 {
@@ -193,11 +197,11 @@ int main(int argc, char *argv[])
                         index.toGlobal
                         (
                             Pstream::myProcNo(),
-                            finalAgglom[patchid][j] + coarsePatchIndex
+                            finalAgglom[patchi][j] + coarsePatchIndex
                         );
                 }
 
-                coarsePatchIndex += max(finalAgglom[patchid]) + 1;
+                coarsePatchIndex += max(finalAgglom[patchi]) + 1;
             }
         }
 

src/OpenFOAM/containers/Lists/ListOps/ListOps.H

@@ -208,7 +208,7 @@ template<class ListType>
 label findMin(const ListType&, const label start=0);
 
 
-//- Find first occurence of given element in sorted list and return index,
+//- Find first occurrence of given element in sorted list and return index,
 //  return -1 if not found. Binary search.
 template<class ListType>
 label findSortedIndex

src/OpenFOAM/db/functionObjects/writeFile/writeFile.C

@@ -182,7 +182,7 @@ Foam::functionObjects::writeFile::~writeFile()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::functionObjects::writeFile::read(const dictionary& dict)
+bool Foam::functionObjects::writeFile::read(const dictionary& dict)
 {
     writePrecision_ =
         dict.lookupOrDefault("writePrecision", IOstream::defaultPrecision());
@@ -190,6 +190,8 @@ void Foam::functionObjects::writeFile::read(const dictionary& dict)
     // Only write on master process
     writeToFile_ = dict.lookupOrDefault("writeToFile", true);
     writeToFile_ = writeToFile_ && Pstream::master();
+
+    return true;
 }
 
 

src/OpenFOAM/db/functionObjects/writeFile/writeFile.H

@@ -148,28 +148,28 @@ public:
     // Member Functions
 
         //- Read
-        void read(const dictionary& dict);
+        virtual bool read(const dictionary& dict);
 
         //- Return access to the file (if only 1)
-        OFstream& file();
+        virtual OFstream& file();
 
         //- Flag to allow writing to file
-        bool writeToFile() const;
+        virtual bool writeToFile() const;
 
         //- Return width of character stream output
-        label charWidth() const;
+        virtual label charWidth() const;
 
         //- Write a commented string to stream
-        void writeCommented(Ostream& os, const string& str) const;
+        virtual void writeCommented(Ostream& os, const string& str) const;
 
         //- Write a tabbed string to stream
-        void writeTabbed(Ostream& os, const string& str) const;
+        virtual void writeTabbed(Ostream& os, const string& str) const;
 
         //- Write a commented header to stream
-        void writeHeader(Ostream& os, const string& str) const;
+        virtual void writeHeader(Ostream& os, const string& str) const;
 
         //- Write the current time to stream
-        void writeTime(Ostream& os) const;
+        virtual void writeTime(Ostream& os) const;
 
         //- Write a (commented) header property and value pair
         template<class Type>

src/OpenFOAM/fields/cloud/cloud.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,7 +36,6 @@ namespace Foam
     word cloud::defaultName("defaultCloud");
 }
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::cloud::cloud(const objectRegistry& obr, const word& cloudName)
@@ -70,4 +69,10 @@ void Foam::cloud::autoMap(const mapPolyMesh&)
 }
 
 
+void Foam::cloud::writeObjects(objectRegistry& obr) const
+{
+    NotImplemented;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/fields/cloud/cloud.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,6 +36,7 @@ SourceFiles
 #define cloud_H
 
 #include "objectRegistry.H"
+#include "IOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,6 +93,24 @@ public:
             //- Remap the cells of particles corresponding to the
             //  mesh topology change
             virtual void autoMap(const mapPolyMesh&);
+
+
+        // I-O
+
+            //- Read particle fields from objects in the obr registry
+            //virtual void readObjects(objectRegistry& obr);
+
+            //- Write particle fields as objects into the obr registry
+            virtual void writeObjects(objectRegistry& obr) const;
+
+            //- Helper to construct IOField on a supplied object registry
+            template<class Type>
+            static IOField<Type>& createIOField
+            (
+                const word& fieldName,
+                const label nParticle,
+                objectRegistry& obr
+            );
 };
 
 
@@ -101,6 +120,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "cloudTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/fields/cloud/cloudTemplates.C

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::IOField<Type>& Foam::cloud::createIOField
+(
+    const word& fieldName,
+    const label nParticle,
+    objectRegistry& obr
+)
+{
+    IOField<Type>* fieldPtr
+    (
+        new IOField<Type>
+        (
+            IOobject
+            (
+                fieldName,
+                obr.time().timeName(),
+                obr,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            nParticle
+        )
+    );
+
+    fieldPtr->store();
+
+    return *fieldPtr;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/subModelBase/subModelBase.C

@@ -167,6 +167,34 @@ bool Foam::subModelBase::writeTime() const
 }
 
 
+bool Foam::subModelBase::getModelDict
+(
+    const word& entryName,
+    dictionary& dict
+) const
+{
+    if (properties_.found(baseName_))
+    {
+        const dictionary& baseDict = properties_.subDict(baseName_);
+
+        if (inLine() && baseDict.found(modelName_))
+        {
+            const dictionary& modelDict = baseDict.subDict(modelName_);
+            dict = modelDict.subOrEmptyDict(entryName);
+            return true;
+        }
+        else if (baseDict.found(modelType_))
+        {
+            const dictionary& modelDict = baseDict.subDict(modelType_);
+            dict = modelDict.subOrEmptyDict(entryName);
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
 void Foam::subModelBase::write(Ostream& os) const
 {
     os  << coeffDict_;

src/OpenFOAM/primitives/subModelBase/subModelBase.H

@@ -179,6 +179,13 @@ public:
 
             // Model properties
 
+                //- Retrieve dictionary, return true if set
+                bool getModelDict
+                (
+                    const word& entryName,
+                    dictionary& dict
+                ) const;
+
                 //- Retrieve generic property from the sub-model
                 template<class Type>
                 void getModelProperty(const word& entryName, Type& value) const;

src/functionObjects/field/Make/files

@@ -69,13 +69,19 @@ streamFunction/streamFunction.C
 valueAverage/valueAverage.C
 fluxSummary/fluxSummary.C
 mapFields/mapFields.C
-reactionSensitivityAnalysis/reactionsSensitivityAnalysisObjects.C 
+reactionSensitivityAnalysis/reactionsSensitivityAnalysisObjects.C
 DESModelRegions/DESModelRegions.C
 
 externalCoupled/externalCoupled.C
 externalCoupled/externalCoupledMixed/externalCoupledMixedFvPatchFields.C
 externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C
 
+extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
+extractEulerianParticles/eulerianParticle/eulerianParticle.C
+
+particleDistribution/particleDistribution.C
+
+
 ddt2/ddt2.C
 zeroGradient/zeroGradient.C
 

src/functionObjects/field/Make/options

@@ -2,6 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
@@ -19,18 +20,22 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/fvAgglomerationMethods/pairPatchAgglomeration/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
+    -lmeshTools \
+    -llagrangian \
+    -ldistributionModels \
+    -lsampling \
+    -lsurfMesh \
     -lfluidThermophysicalModels \
     -lincompressibleTransportModels \
     -lturbulenceModels \
     -lcompressibleTransportModels \
     -lincompressibleTurbulenceModels \
     -lcompressibleTurbulenceModels \
-    -lmeshTools \
-    -lsampling \
-    -lsurfMesh \
     -lchemistryModel \
-    -lreactionThermophysicalModels
+    -lreactionThermophysicalModels \
+    -lpairPatchAgglomeration

src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.C

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "eulerianParticle.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::eulerianParticle::eulerianParticle()
+:
+    faceIHit(-1),
+    VC(vector::zero),
+    VU(vector::zero),
+    V(0),
+    time(0)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const functionObjects::eulerianParticle& p
+)
+{
+    os  << p.faceIHit << token::SPACE
+        << p.VC << token::SPACE
+        << p.VU << token::SPACE
+        << p.V << token::SPACE
+        << p.time;
+
+    return os;
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    functionObjects::eulerianParticle& p
+)
+{
+    is  >> p.faceIHit
+        >> p.VC
+        >> p.VU
+        >> p.V
+        >> p.time;
+
+    return is;
+}
+
+
+void Foam::functionObjects::eulerianParticle::write(Ostream& os) const
+{
+    scalar pDiameter = cbrt(6*V/constant::mathematical::pi);
+    vector U = VU/(V + ROOTVSMALL);
+    vector C = VC/(V + ROOTVSMALL);
+
+    os  << time << token::SPACE
+        << faceIHit << token::SPACE
+        << C << token::SPACE
+        << pDiameter << token::SPACE
+        << U << token::SPACE
+        << endl;
+}
+
+
+Foam::dictionary Foam::functionObjects::eulerianParticle::writeDict() const
+{
+    scalar pDiameter = cbrt(6*V/constant::mathematical::pi);
+    vector U = VU/(V + ROOTVSMALL);
+    vector C = VC/(V + ROOTVSMALL);
+
+    dictionary dict;
+    dict.add("time", time);
+    dict.add("meshFace", faceIHit);
+    dict.add("position", C);
+    dict.add("diameter", pDiameter);
+    dict.add("U", U);
+
+    return dict;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::eulerianParticle
+
+Description
+    Lightweight class to store particle data derived from VOF calculations,
+    with special handling for input, output and parallel reduction.
+
+SourceFiles
+    eulerianParticle.H
+    eulerianParticle.C
+    eulerianParticleTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_eulerianParticle_H
+#define functionObjects_eulerianParticle_H
+
+#include "label.H"
+#include "scalar.H"
+#include "vector.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+namespace functionObjects
+{
+    class eulerianParticle;
+}
+
+// Forward declaration of friend functions and operators
+Istream& operator>>(Istream&, functionObjects::eulerianParticle&);
+Ostream& operator<<(Ostream&, const functionObjects::eulerianParticle&);
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class eulerianParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class eulerianParticle
+{
+
+public:
+
+    // Public data
+
+        //- Index of face in faceZone that this particle hits.  Also used to
+        //  identify the index of the coarse face of the surface agglomeration
+        //  Note: value of -1 used to indicate that the particle has not
+        //        been initialised
+        label faceIHit;
+
+        //- Volume multiplied by face centres [m4]
+        vector VC;
+
+        //- Volume multiplied by velocity [m4/s]
+        vector VU;
+
+        //- Volume [m3]
+        scalar V;
+
+        //- Injection time - set at collection [s]
+        scalar time;
+
+
+    //- Constructor
+    eulerianParticle();
+
+
+    // Public Member Functions
+
+        //- Write to stream
+        void write(Ostream& os) const;
+
+        //- Write to dictionary
+        Foam::dictionary writeDict() const;
+
+
+    // Operators
+
+        friend bool operator==
+        (
+            const eulerianParticle& a,
+            const eulerianParticle& b
+        )
+        {
+            return
+                a.faceIHit == b.faceIHit
+             && a.VC == b.VC
+             && a.VU == b.VU
+             && a.V == b.V
+             && a.time == b.time;
+        }
+
+        friend bool operator!=
+        (
+            const eulerianParticle& a,
+            const eulerianParticle& b
+        )
+        {
+            return !(a == b);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "eulerianParticleTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticleTemplates.C

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class VOFParticle>
+class sumParticleOp
+{
+    public:
+    eulerianParticle operator()
+    (
+        const eulerianParticle& p0,
+        const eulerianParticle& p1
+    ) const
+    {
+        if ((p0.faceIHit != -1) && (p1.faceIHit == -1))
+        {
+            return p0;
+        }
+        else if ((p0.faceIHit == -1) && (p1.faceIHit != -1))
+        {
+            return p1;
+        }
+        else if ((p0.faceIHit != -1) && (p1.faceIHit != -1))
+        {
+            // Choose particle with the largest collected volume and
+            // accumulate total volume
+            if (p0.V > p1.V)
+            {
+                eulerianParticle p = p0;
+                p.V = p0.V + p1.V;
+                p.VC = p0.VC + p1.VC;
+                p.VU = p0.VU + p1.VU;
+                return p;
+            }
+            else
+            {
+                eulerianParticle p = p1;
+                p.V = p0.V + p1.V;
+                p.VC = p0.VC + p1.VC;
+                p.VU = p0.VU + p1.VU;
+                return p;
+            }
+        }
+        else
+        {
+            eulerianParticle p;
+            return p;
+        }
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //

src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C

@@ -0,0 +1,681 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extractEulerianParticles.H"
+#include "regionSplit2D.H"
+#include "mathematicalConstants.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "surfaceInterpolate.H"
+#include "pairPatchAgglomeration.H"
+#include "emptyPolyPatch.H"
+#include "coupledPolyPatch.H"
+#include "binned.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(extractEulerianParticles, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        extractEulerianParticles,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::fileName
+Foam::functionObjects::extractEulerianParticles::dictBaseFileDir() const
+{
+    fileName baseDir(".."); //  = mesh_.time().path();
+
+    if (Pstream::parRun())
+    {
+        // Put in undecomposed case (Note: gives problems for
+        // distributed data running)
+        baseDir = baseDir/".."/"postProcessing";
+    }
+    else
+    {
+        baseDir = baseDir/"postProcessing";
+    }
+
+    return baseDir;
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::checkFaceZone()
+{
+    DebugInFunction << endl;
+
+    zoneID_ = mesh_.faceZones().findZoneID(faceZoneName_);
+    if (zoneID_ == -1)
+    {
+        FatalErrorInFunction
+            << "Unable to find faceZone " << faceZoneName_
+            << ".  Available faceZones are: " << mesh_.faceZones().names()
+            << exit(FatalError);
+    }
+
+    const faceZone& fz = mesh_.faceZones()[zoneID_];
+    const label nFaces = fz.size();
+    const label allFaces = returnReduce(nFaces, sumOp<label>());
+
+    if (allFaces < nInjectorLocations_)
+    {
+        FatalErrorInFunction
+            << "faceZone " << faceZoneName_
+            << ": Number of faceZone faces (" << allFaces
+            << ") is less than the number of requested locations ("
+            << nInjectorLocations_ << ")."
+            << exit(FatalError);
+    }
+
+    Info<< type() << " " << name() << " output:" << nl
+        << "    faceZone : " << faceZoneName_ << nl
+        << "    faces    : " << allFaces << nl
+        << endl;
+
+    // Initialise old iteration blocked faces
+    // Note: for restart, this info needs to be written/read
+    regions0_.setSize(fz.size(), -1);
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::initialiseBins()
+{
+    DebugInFunction << endl;
+
+    if (!nInjectorLocations_)
+    {
+        return;
+    }
+
+    const faceZone& fz = mesh_.faceZones()[zoneID_];
+
+    // Agglomerate faceZone faces into nInjectorLocations_ global locations
+    const indirectPrimitivePatch patch
+    (
+        IndirectList<face>(mesh_.faces(), fz),
+        mesh_.points()
+    );
+
+    const label nFaces = fz.size();
+    label nLocations = nInjectorLocations_;
+
+    if (Pstream::parRun())
+    {
+        label nGlobalFaces = returnReduce(nFaces, sumOp<label>());
+        scalar fraction = scalar(nFaces)/scalar(nGlobalFaces);
+        nLocations = ceil(fraction*nInjectorLocations_);
+        if (debug)
+        {
+            Pout<< "nFaces:" << nFaces
+                << ", nGlobalFaces:" << nGlobalFaces
+                << ", fraction:" << fraction
+                << ", nLocations:" << nLocations
+                << endl;
+        }
+    }
+
+    pairPatchAgglomeration ppa
+    (
+        patch.localFaces(),
+        patch.localPoints(),
+        10,
+        50,
+        nLocations,
+        labelMax,
+        180
+    );
+
+    ppa.agglomerate();
+
+    label nCoarseFaces = 0;
+    if (nFaces != 0)
+    {
+        fineToCoarseAddr_ = ppa.restrictTopBottomAddressing();
+        nCoarseFaces = max(fineToCoarseAddr_) + 1;
+    }
+
+    globalCoarseFaces_ = globalIndex(nCoarseFaces);
+
+    Info<< "Created " << returnReduce(nCoarseFaces, sumOp<label>())
+        << " coarse faces" << endl;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::functionObjects::extractEulerianParticles::phiU() const
+{
+    DebugInFunction << endl;
+
+    const surfaceScalarField& phi
+    (
+        mesh_.lookupObject<surfaceScalarField>(phiName_)
+    );
+
+    if (phi.dimensions() == dimMass/dimTime)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        return phi/fvc::interpolate(rho);
+    }
+
+    return phi;
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::setBlockedFaces
+(
+    const surfaceScalarField& alphaf,
+    const faceZone& fz,
+    boolList& blockedFaces
+)
+{
+    DebugInFunction << endl;
+
+    // Initialise storage for patch and patch-face indices where faceZone
+    // intersects mesh patch(es)
+    patchIDs_.setSize(fz.size(), -1);
+    patchFaceIDs_.setSize(fz.size(), -1);
+
+    label nBlockedFaces = 0;
+    forAll(fz, localFacei)
+    {
+        const label facei = fz[localFacei];
+
+        if (mesh_.isInternalFace(facei))
+        {
+            if (alphaf[facei] > alphaThreshold_)
+            {
+                blockedFaces[localFacei] = true;
+            }
+        }
+        else
+        {
+            label patchi = mesh_.boundaryMesh().whichPatch(facei);
+            label patchFacei = -1;
+
+            const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+            const scalarField& alphafp = alphaf.boundaryField()[patchi];
+
+            if (isA<coupledPolyPatch>(pp))
+            {
+                const coupledPolyPatch& cpp =
+                    refCast<const coupledPolyPatch>(pp);
+
+                if (cpp.owner())
+                {
+                    patchFacei = cpp.whichFace(facei);
+                }
+            }
+            else if (!isA<emptyPolyPatch>(pp))
+            {
+                patchFacei = pp.whichFace(facei);
+            }
+
+            if (patchFacei == -1)
+            {
+                patchi = -1;
+            }
+            else if (alphafp[patchFacei] > alphaThreshold_)
+            {
+                blockedFaces[localFacei] = true;
+            }
+
+            patchIDs_[localFacei] = patchi;
+            patchFaceIDs_[localFacei] = patchFacei;
+        }
+    }
+
+    DebugInFunction << "Number of blocked faces: " << nBlockedFaces << endl;
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::collectParticles
+(
+    const scalar time,
+    const boolList& collectParticle
+)
+{
+    DebugInFunction << "collectParticle: " << collectParticle << endl;
+
+    // Collect particles on local processors that have passed through faceZone
+    forAll(collectParticle, regioni)
+    {
+        if (!collectParticle[regioni])
+        {
+            continue;
+        }
+
+        Map<label>::const_iterator iter = regionToParticleMap_.find(regioni);
+        eulerianParticle p = particles_[iter()];
+
+        if (p.faceIHit != -1 && nInjectorLocations_)
+        {
+            // Use coarse face index for tag output
+            label coarseFacei = fineToCoarseAddr_[p.faceIHit];
+            p.faceIHit = globalCoarseFaces_.toGlobal(coarseFacei);
+        }
+
+        reduce(p, sumParticleOp<eulerianParticle>());
+
+        const scalar pDiameter = cbrt(6.0*p.V/constant::mathematical::pi);
+
+        if ((pDiameter > minDiameter_) && (pDiameter < maxDiameter_))
+        {
+            if (Pstream::master())
+            {
+                const scalar d = cbrt(6*p.V/constant::mathematical::pi);
+                const point position = p.VC/(p.V + ROOTVSMALL);
+                const vector U = p.VU/(p.V + ROOTVSMALL);
+                label tag = -1;
+                if (nInjectorLocations_)
+                {
+                    tag = p.faceIHit;
+                }
+
+                injectedParticle* ip = new injectedParticle
+                (
+                    mesh_,
+                    position,
+                    tag,
+                    time,
+                    d,
+                    U
+                );
+
+                cloud_.addParticle(ip);
+            }
+
+            nCollectedParticles_++;
+        }
+        else
+        {
+            // Discard particles over/under diameter threshold
+            nDiscardedParticles_++;
+            discardedVolume_ += p.V;
+        }
+    }
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::calculateAddressing
+(
+    const label nRegionsOld,
+    const label nRegionsNew,
+    const scalar time,
+    labelList& regionFaceIDs
+)
+{
+    DebugInFunction << endl;
+
+    // New region can only point to one old region
+    // Old region can only point to one new region.  If old region intersects
+    // multiple new regions, select max of new region indices.
+
+    labelList oldToNewRegion(nRegionsOld, -1);
+    labelList newToOldRegion(nRegionsNew, -1);
+
+    forAll(regionFaceIDs, facei)
+    {
+        label newRegioni = regionFaceIDs[facei];
+        label oldRegioni = regions0_[facei];
+
+        if (newRegioni != -1)
+        {
+            newToOldRegion[newRegioni] = oldRegioni;
+
+            if (oldRegioni != -1)
+            {
+                // New region linked to old (so can accumulate particle data)
+                // Old region might already be mapped to a new region
+                oldToNewRegion[oldRegioni] =
+                    max(newRegioni, oldToNewRegion[oldRegioni]);
+            }
+        }
+    }
+
+    // Need to re-number the new region indices based on knowledge of which
+    // old region they intersect.  After operations, there should be a
+    // one-to-one match between the old and new regions.
+
+    // Ensure all old regions point to the same new regions (if present)
+    Pstream::listCombineGather(oldToNewRegion, maxEqOp<label>());
+    Pstream::listCombineScatter(oldToNewRegion);
+
+    // Any new region that points to an old region should be renumbered to the
+    // new region specified by the oldToNewRegion index
+
+    if (oldToNewRegion.size())
+    {
+        // Create corrected new to new addressing
+        labelList newToNewRegionCorr(newToOldRegion.size(), -1);
+        forAll(newToOldRegion, newRegioni)
+        {
+            label oldRegioni = newToOldRegion[newRegioni];
+            if (oldRegioni != -1)
+            {
+                label newRegionICorr = oldToNewRegion[oldRegioni];
+                newToNewRegionCorr[newRegioni] = newRegionICorr;
+                newToOldRegion[newRegioni] = oldRegioni;
+            }
+        }
+
+        // Renumber the new (current) face region IDs
+        forAll(regionFaceIDs, facei)
+        {
+            label newRegioni = regionFaceIDs[facei];
+
+            if (newRegioni != -1)
+            {
+                label newRegionICorr = newToNewRegionCorr[newRegioni];
+
+                // Note: newRegionICorr can be -1 if the region is new, since
+                // the address corrections are based on inverting the
+                // old-to-new addressing
+                if (newRegionICorr != -1)
+                {
+                    regionFaceIDs[facei] = newRegionICorr;
+                }
+            }
+        }
+
+
+        boolList collectParticleFlag(nRegionsOld, true);
+        forAll(oldToNewRegion, oldRegioni)
+        {
+            label newRegioni = oldToNewRegion[oldRegioni];
+            if (newRegioni != -1)
+            {
+                collectParticleFlag[oldRegioni] = false;
+            }
+        }
+
+        // Collect particles whose IDs are no longer active
+        collectParticles(time, collectParticleFlag);
+    }
+
+
+    // Re-order collection bins
+
+    Map<label> newRegionToParticleMap(nRegionsNew);
+    List<eulerianParticle> newParticles(nRegionsNew);
+    label particlei = 0;
+
+    forAll(newToOldRegion, newRegioni)
+    {
+        label oldRegioni = newToOldRegion[newRegioni];
+        if (oldRegioni == -1)
+        {
+            // No mapping from old to new - this is a new particle
+            newRegionToParticleMap.insert(newRegioni, particlei);
+            particlei++;
+        }
+        else
+        {
+            // Update addressing for old to new regions
+            label oldParticlei = regionToParticleMap_[oldRegioni];
+            if (newRegionToParticleMap.insert(newRegioni, particlei))
+            {
+                // First time this particle has been seen
+                newParticles[particlei] = particles_[oldParticlei];
+                particlei++;
+            }
+            else
+            {
+                // Combine with existing contribution
+                label newParticlei = newRegionToParticleMap[newRegioni];
+                newParticles[newParticlei] =
+                    sumParticleOp<eulerianParticle>()
+                    (
+                        newParticles[newParticlei],
+                        particles_[oldParticlei]
+                    );
+            }
+        }
+    }
+
+    // Reset the particles list and addressing for latest available info
+    particles_.transfer(newParticles);
+    regionToParticleMap_ = newRegionToParticleMap;
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::accumulateParticleInfo
+(
+    const surfaceScalarField& alphaf,
+    const surfaceScalarField& phi,
+    const labelList& regionFaceIDs,
+    const faceZone& fz
+)
+{
+    DebugInFunction << endl;
+
+    const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
+    const surfaceVectorField Uf(fvc::interpolate(U));
+
+    const scalar deltaT = mesh_.time().deltaTValue();
+    const pointField& faceCentres = mesh_.faceCentres();
+
+    forAll(regionFaceIDs, localFacei)
+    {
+        const label newRegioni = regionFaceIDs[localFacei];
+
+        if (newRegioni != -1)
+        {
+            const label particlei = regionToParticleMap_[newRegioni];
+            const label meshFacei = fz[localFacei];
+            eulerianParticle& p = particles_[particlei];
+
+            if (p.faceIHit < 0)
+            {
+                // New particle - does not exist in particles_ list
+                p.faceIHit = localFacei;
+                p.V = 0;
+                p.VC = vector::zero;
+                p.VU = vector::zero;
+            }
+
+            // Accumulate particle properties
+            scalar magPhii = mag(faceValue(phi, localFacei, meshFacei));
+            vector Ufi = faceValue(Uf, localFacei, meshFacei);
+            scalar dV = magPhii*deltaT;
+            p.V += dV;
+            p.VC += dV*faceCentres[meshFacei];
+            p.VU += dV*Ufi;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::extractEulerianParticles::extractEulerianParticles
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(runTime, name),
+    cloud_(mesh_, "eulerianParticleCloud"),
+    faceZoneName_(word::null),
+    zoneID_(-1),
+    patchIDs_(),
+    patchFaceIDs_(),
+    alphaName_("alpha"),
+    alphaThreshold_(0.1),
+    UName_("U"),
+    rhoName_("rho"),
+    phiName_("phi"),
+    nInjectorLocations_(0),
+    fineToCoarseAddr_(),
+    globalCoarseFaces_(),
+    regions0_(),
+    nRegions0_(0),
+    particles_(),
+    regionToParticleMap_(),
+    minDiameter_(ROOTVSMALL),
+    maxDiameter_(GREAT),
+    nCollectedParticles_(0),
+    nDiscardedParticles_(0),
+    discardedVolume_(0)
+{
+    if (mesh_.nSolutionD() != 3)
+    {
+        FatalErrorInFunction
+            << name << " function object only applicable to 3-D cases"
+            << exit(FatalError);
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::extractEulerianParticles::~extractEulerianParticles()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::extractEulerianParticles::read
+(
+    const dictionary& dict
+)
+{
+    DebugInFunction << endl;
+
+    if (fvMeshFunctionObject::read(dict) && writeFile::read(dict))
+    {
+        dict.lookup("faceZone") >> faceZoneName_;
+        dict.lookup("alpha") >> alphaName_;
+
+        dict.readIfPresent("alphaThreshold", alphaThreshold_);
+        dict.readIfPresent("U", UName_);
+        dict.readIfPresent("rho", rhoName_);
+        dict.readIfPresent("phi", phiName_);
+        dict.readIfPresent("nLocations", nInjectorLocations_);
+        dict.readIfPresent("minDiameter", minDiameter_);
+        dict.readIfPresent("maxDiameter", maxDiameter_);
+
+        checkFaceZone();
+
+        if (nInjectorLocations_)
+        {
+            initialiseBins();
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::extractEulerianParticles::execute()
+{
+    DebugInFunction << endl;
+
+    Log << type() << " " << name() << " output:" << nl;
+
+    const volScalarField& alpha =
+        mesh_.lookupObject<volScalarField>(alphaName_);
+
+    const surfaceScalarField alphaf
+    (
+        typeName + ":alphaf",
+        fvc::interpolate(alpha)
+    );
+
+    const faceZone& fz = mesh_.faceZones()[zoneID_];
+    const indirectPrimitivePatch patch
+    (
+        IndirectList<face>(mesh_.faces(), fz),
+        mesh_.points()
+    );
+
+    // Set the blocked faces, i.e. where alpha > alpha threshold value
+    boolList blockedFaces(fz.size(), false);
+    setBlockedFaces(alphaf, fz, blockedFaces);
+
+    // Split the  faceZone according to the blockedFaces
+    // - Returns a list of (disconnected) region index per face zone face
+    regionSplit2D regionFaceIDs(mesh_, patch, blockedFaces);
+
+    // Global number of regions
+    const label nRegionsNew = regionFaceIDs.nRegions();
+
+    // Calculate the addressing between the old and new region information
+    // Also collects particles that have traversed the faceZone
+    calculateAddressing
+    (
+        nRegions0_,
+        nRegionsNew,
+        mesh_.time().value(),
+        regionFaceIDs
+    );
+
+    // Process latest region information
+    tmp<surfaceScalarField> tphi = phiU();
+    accumulateParticleInfo(alphaf, tphi(), regionFaceIDs, fz);
+
+    // Reset the blocked faces for the next integration step
+    nRegions0_ = nRegionsNew;
+    regions0_ = regionFaceIDs;
+
+    Log << "    Collected particles: " << nCollectedParticles_ << nl
+        << "    Discarded particles: " << nDiscardedParticles_ << nl
+        << "    Discarded volume: " << discardedVolume_ << nl
+        << endl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::extractEulerianParticles::write()
+{
+    DebugInFunction << endl;
+
+    cloud_.write();
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H

@@ -0,0 +1,295 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::extractEulerianParticles
+Group
+    grpFieldFunctionObjects
+
+Description
+    Generates particle size information from Eulerian calculations, e.g. VoF.
+
+Usage
+    extractEulerianParticles1
+    {
+        type            extractEulerianParticles;
+        libs            ("libfieldFunctionObjects.so");
+        ...
+        faceZone        f0;
+        nLocations      10;
+        alphaName       alpha.water;
+        UName           U;
+        rhoName         rho;
+        phiName         phi;
+    }
+    \endverbatim
+
+    where the entries comprise:
+    \table
+        Property     | Description             | Required | Default value
+        type         | type name: extractEulerianParticles | yes |
+        faceZone     | Name of faceZone used as collection surface | yes |
+        nLocations   | Number of injection bins to generate  | yes |
+        aplhaName    | Name of phase indicator field | yes   |
+        rhoName      | Name of density field   |      yes |
+        phiNane      | Name of flux field      |      yes |
+    \endtable
+
+SourceFiles
+    extractEulerianParticles.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_extractEulerianParticles_H
+#define functionObjects_extractEulerianParticles_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "runTimeSelectionTables.H"
+#include "polyMesh.H"
+#include "surfaceFieldsFwd.H"
+#include "vectorList.H"
+#include "globalIndex.H"
+#include "cachedRandom.H"
+#include "eulerianParticle.H"
+#include "IOdictionary.H"
+#include "injectedParticleCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+/*---------------------------------------------------------------------------*\
+           Class extractEulerianParticlesFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class extractEulerianParticles
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+
+protected:
+
+    // Protected data
+
+        //- Storage for collected particles
+        injectedParticleCloud cloud_;
+
+
+        // faceZone info
+
+            //- Name of faceZone to sample
+            word faceZoneName_;
+
+            //- Index of the faceZone
+            label zoneID_;
+
+            //- Patch indices where faceZone face intersect patch
+            labelList patchIDs_;
+
+            //- Patch face indices where faceZone face intersect patch
+            labelList patchFaceIDs_;
+
+
+        // Field names
+
+            //- Name of phase fraction field
+            word alphaName_;
+
+            //- Value of phase fraction used to identify particle boundaries
+            scalar alphaThreshold_;
+
+            //- Name of the velocity field, default = U
+            word UName_;
+
+            //- Name of the density field, default = rho
+            word rhoName_;
+
+            //- Name of the flux field, default ="rho"
+            word phiName_;
+
+
+        // Agglomeration
+
+            //- Number of sample locations to generate
+            label nInjectorLocations_;
+
+            //- Agglomeration addressing from fine to coarse (local proc only)
+            labelList fineToCoarseAddr_;
+
+            //- Global coarse face addressing
+            globalIndex globalCoarseFaces_;
+
+
+        // Particle collection info
+
+            //- Region indices in faceZone faces from last iteration
+            labelList regions0_;
+
+            //- Number of regions from last iteration
+            label nRegions0_;
+
+            //- Particle properties (partial, being accumulated)
+            List<eulerianParticle> particles_;
+
+            //- Map from region to index in particles_ list
+            Map<label> regionToParticleMap_;
+
+            //- Minimum diameter (optional)
+            //  Can be used to filter out 'small' particles
+            scalar minDiameter_;
+
+            //- Maximum diameter (optional)
+            //  Can be used to filter out 'large' particles
+            scalar maxDiameter_;
+
+
+        // Statistics
+
+            //- Total number of collected particles
+            label nCollectedParticles_;
+
+            //- Total number of discarded particles
+            label nDiscardedParticles_;
+
+            //- Total discarded volume [m3]
+            scalar discardedVolume_;
+
+
+    // Protected Member Functions
+
+        //- Return the base directory for dictionary output
+        virtual fileName dictBaseFileDir() const;
+
+        //- Check that the faceZone is valid
+        virtual void checkFaceZone();
+
+        //- Initialise the particle collection bins
+        virtual void initialiseBins();
+
+        //- Return the volumetric flux
+        virtual tmp<surfaceScalarField> phiU() const;
+
+        //- Set the blocked faces, i.e. where alpha > alpha threshold value
+        virtual void setBlockedFaces
+        (
+            const surfaceScalarField& alphaf,
+            const faceZone& fz,
+            boolList& blockedFaces
+        );
+
+        //- Calculate the addressing between regions between iterations
+        virtual void calculateAddressing
+        (
+            const label nRegionsOld,
+            const label nRegionsNew,
+            const scalar time,
+            labelList& regionFaceIDs
+        );
+
+        //- Collect particles that have passed through the faceZone
+        virtual void collectParticles
+        (
+            const scalar time,
+            const boolList& collectParticle
+        );
+
+        //- Process latest region information
+        virtual void accumulateParticleInfo
+        (
+            const surfaceScalarField& alphaf,
+            const surfaceScalarField& phi,
+            const labelList& regionFaceIDs,
+            const faceZone& fz
+        );
+
+        template<class Type>
+        inline Type faceValue
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
+            const label localFaceI,
+            const label globalFaceI
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        extractEulerianParticles(const extractEulerianParticles&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extractEulerianParticles&);
+
+
+public:
+
+    // Static data members
+
+        //- Runtime type information
+        TypeName("extractEulerianParticles");
+
+
+    // Constructors
+
+        //- Construct from components
+        extractEulerianParticles
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~extractEulerianParticles();
+
+
+    // Member Functions
+
+        //- Read the field min/max data
+        virtual bool read(const dictionary&);
+
+        //- Execute
+        virtual bool execute();
+
+        //- Write
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "extractEulerianParticlesTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticlesTemplates.C

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMesh.H"
+
+template<class Type>
+Type Foam::functionObjects::extractEulerianParticles::faceValue
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& field,
+    const label localFacei,
+    const label globalFacei
+) const
+{
+    if (mesh_.isInternalFace(globalFacei))
+    {
+        return field[globalFacei];
+    }
+    else
+    {
+        label patchi = patchIDs_[localFacei];
+        label pFacei = patchFaceIDs_[localFacei];
+        if (patchi != 0)
+        {
+            return field.boundaryField()[patchi][pFacei];
+        }
+        else
+        {
+            return pTraits<Type>::zero;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/fluxSummary/fluxSummary.H

@@ -41,7 +41,7 @@ Usage
     fluxSummary1
     {
         type        fluxSummary;
-        libs        ("libutilityFunctionObjects.so");
+        libs        ("libfieldFunctionObjects.so");
         ...
         write       yes;
         log         yes;
@@ -240,7 +240,7 @@ public:
 
     // Member Functions
 
-        //- Read the field min/max data
+        //- Read the field fluxSummary data
         virtual bool read(const dictionary&);
 
         //- Execute, currently does nothing

src/functionObjects/field/particleDistribution/particleDistribution.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "particleDistribution.H"
+#include "addToRunTimeSelectionTable.H"
+#include "general.H"
+#include "fvMesh.H"
+#include "cloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(particleDistribution, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        particleDistribution,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::particleDistribution::particleDistribution
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(runTime, name, typeName, dict),
+    cloudName_("unknown-cloudName"),
+    fieldNames_(),
+    tagFieldName_("unknown-tagFieldName"),
+    distributionBinWidth_(0),
+    rndGen_(1234, -1)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::particleDistribution::~particleDistribution()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::particleDistribution::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict) && writeFile::read(dict))
+    {
+        dict.lookup("cloud") >> cloudName_;
+        dict.lookup("fields") >> fieldNames_;
+        dict.lookup("tagField") >> tagFieldName_;
+        dict.lookup("distributionBinWidth") >> distributionBinWidth_;
+
+        Info<< type() << " " << name() << " output:"
+            << "    Processing cloud : " << cloudName_ << nl
+            << endl;
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::particleDistribution::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::particleDistribution::write()
+{
+    Log << type() << " " << name() << " output:" << endl;
+
+    if (!mesh_.foundObject<cloud>(cloudName_))
+    {
+        WarningInFunction
+            << "Unable to find cloud " << cloudName_
+            << " in the mesh database" << endl;
+
+        return false;
+    }
+
+    const cloud& c = mesh_.lookupObject<cloud>(cloudName_);
+
+    objectRegistry cloudObr
+    (
+        IOobject
+        (
+            name() & "CloudRegistry",
+            mesh_.time().timeName(),
+            cloud::prefix,
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    c.writeObjects(cloudObr);
+
+    List<DynamicList<label>> tagAddr;
+    if
+    (
+        tagFieldName_ != "none"
+     && cloudObr.foundObject<IOField<scalar>>(tagFieldName_)
+    )
+    {
+        // Tag field present - generate distribution per tag
+        const IOField<label>& tag =
+            cloudObr.lookupObject<IOField<label>>(tagFieldName_);
+        const HashSet<label> tagMap(tag);
+        const label tagMax = tagMap.size();
+
+        List<DynamicList<label>> tagAddr(tagMax);
+        forAll(tag, i)
+        {
+            label newTag = tagMap[tag[i]];
+            tagAddr[newTag].append(i);
+        }
+    }
+
+    bool ok = false;
+    forAll(fieldNames_, i)
+    {
+        const word fName = fieldNames_[i];
+        ok = ok || processField<scalar>(cloudObr, fName, tagAddr);
+        ok = ok || processField<vector>(cloudObr, fName, tagAddr);
+        ok = ok || processField<tensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<sphericalTensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<symmTensor>(cloudObr, fName, tagAddr);
+        ok = ok || processField<tensor>(cloudObr, fName, tagAddr);
+
+        if (log && !ok)
+        {
+            WarningInFunction
+                << "Unable to find field " << fName << " in the "
+                << cloudName_ << " cloud database" << endl;
+        }
+    }
+
+
+    return true;
+}
+
+
+void Foam::functionObjects::particleDistribution::generateDistribution
+(
+    const word& fieldName,
+    const scalarField& field,
+    const label tag
+)
+{
+    Ostream& os = file();
+
+    if (tag != -1)
+    {
+        os  << tag << token::TAB;
+    }
+
+    distributionModels::general distribution
+    (
+        field,
+        distributionBinWidth_,
+        rndGen_
+    );
+
+    distribution.writeData(os);
+
+    os  << endl;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/particleDistribution/particleDistribution.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::particleDistribution
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Generates a particle distrbution for lagrangian data.
+
+Usage
+    \verbatim
+    particleDistribution1
+    {
+        type        particleDistribution;
+        libs        ("libfieldFunctionObjects.so");
+        ...
+        cloud       "myCloud";
+        field       "d"
+        distributionBinWidth 0.1;
+    }
+    \endverbatim
+
+    Where the entries comprise:
+    \table
+        Property    | Description              |  Required  | Default value
+        type        | Type name: particleDistribution | yes |
+        cloud       | Name of cloud to process | Yes        |
+        field       | Name of cloud field to process | Yes  |
+        distributionBinWidth | Width of distribution bins | yes |
+    \endtable
+
+See also
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObjects::writeFile
+    Foam::functionObjects::timeControl
+
+SourceFiles
+    particleDistribution.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_particleDistribution_H
+#define functionObjects_particleDistribution_H
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "scalarField.H"
+#include "cachedRandom.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class particleDistribution Declaration
+\*---------------------------------------------------------------------------*/
+
+class particleDistribution
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+protected:
+
+    // Protected data
+
+        //- Cloud name
+        word cloudName_;
+
+        //- Field name
+        wordList fieldNames_;
+
+        //- Tag field name - used to filter the particles into groups
+        word tagFieldName_;
+
+        //- Distribution bin width
+        scalar distributionBinWidth_;
+
+        //- Random number generator - used by distribution models
+        cachedRandom rndGen_;
+
+
+    // Private Member Functions
+
+        //- Generate the distribution
+        void generateDistribution
+        (
+            const word& fieldName,
+            const scalarField& field,
+            const label tag = -1
+        );
+
+        // Helper function to retrieve the particle data
+        template<class Type>
+        bool processField
+        (
+            const objectRegistry& obr,
+            const word& fieldName,
+            const List<DynamicList<label>>& addr
+        );
+
+        //- Disallow default bitwise copy construct
+        particleDistribution(const particleDistribution&) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const particleDistribution&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("particleDistribution");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        particleDistribution
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~particleDistribution();
+
+
+    // Member Functions
+
+        //- Read the particleDistribution data
+        virtual bool read(const dictionary&);
+
+        //- Execute, currently does nothing
+        virtual bool execute();
+
+        //- Write the particleDistribution
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "particleDistributionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/particleDistribution/particleDistributionTemplates.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "IOField.H"
+
+template<class Type>
+bool Foam::functionObjects::particleDistribution::processField
+(
+    const objectRegistry& obr,
+    const word& fieldName,
+    const List<DynamicList<label>>& addr
+)
+{
+    if (obr.foundObject<IOField<Type>>(fieldName))
+    {
+        const IOField<Type>& field =
+            obr.lookupObject<IOField<Type>>(fieldName);
+
+        if (addr.size())
+        {
+            forAll(addr, i)
+            {
+                const Field<Type> subField(field, addr[i]);
+                for (direction d = 0; d < pTraits<Type>::nComponents; ++d)
+                {
+                    generateDistribution(fieldName, subField.component(d), i);
+                }
+            }
+        }
+        else
+        {
+            for (direction d = 0; d < pTraits<Type>::nComponents; ++d)
+            {
+                generateDistribution(fieldName, field.component(d));
+            }
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,17 +39,28 @@ void Foam::pairPatchAgglomeration::compactLevels(const label nCreatedLevels)
 
 bool Foam::pairPatchAgglomeration::continueAgglomerating
 (
-    const label nCoarseFaces
+    const label nLocal,
+    const label nLocalOld
 )
 {
-    // Check the need for further agglomeration on all processors
-    label localnCoarseFaces = nCoarseFaces;
-    bool contAgg = localnCoarseFaces >= nFacesInCoarsestLevel_;
-    return contAgg;
+    // Keep agglomerating
+    // - if global number of faces is still changing
+    // - and if local number of faces still too large (on any processor)
+    //       or if global number of faces still too large
+
+    label nGlobal = returnReduce(nLocal, sumOp<label>());
+    label nGlobalOld = returnReduce(nLocalOld, sumOp<label>());
+
+    return
+    (
+        returnReduce(nLocal > nFacesInCoarsestLevel_, orOp<bool>())
+     || nGlobal > nGlobalFacesInCoarsestLevel_
+    )
+    && nGlobal != nGlobalOld;
 }
 
 
-void Foam::pairPatchAgglomeration::setBasedEdgeWeights()
+void Foam::pairPatchAgglomeration::setLevel0EdgeWeights()
 {
     const bPatch& coarsePatch = patchLevels_[0];
     forAll(coarsePatch.edges(), i)
@@ -78,7 +89,7 @@ void Foam::pairPatchAgglomeration::setBasedEdgeWeights()
                     facePairWeight_.insert(edgeCommon, edgeLength);
                 }
 
-                if (cosI < Foam::cos(degToRad(featureAngle_)))
+                if (mag(cosI) < Foam::cos(degToRad(featureAngle_)))
                 {
                     facePairWeight_[edgeCommon] = -1.0;
                 }
@@ -107,14 +118,12 @@ void Foam::pairPatchAgglomeration::setEdgeWeights
     const label fineLevelIndex
 )
 {
-
     const bPatch& coarsePatch = patchLevels_[fineLevelIndex];
-
     const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
     const label nCoarseI =  max(fineToCoarse) + 1;
     labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));
 
-    HashSet<edge, Hash<edge>> fineFeaturedFaces(coarsePatch.nEdges()/10);
+    HashSet<edge, Hash<edge> > fineFeaturedFaces(coarsePatch.nEdges()/10);
 
     // Map fine faces with featured edge into coarse faces
     forAllConstIter(EdgeMap<scalar>, facePairWeight_, iter)
@@ -147,7 +156,6 @@ void Foam::pairPatchAgglomeration::setEdgeWeights
             if (eFaces.size() == 2)
             {
                 const edge edgeCommon = edge(eFaces[0], eFaces[1]);
-
                 if (facePairWeight_.found(edgeCommon))
                 {
                     facePairWeight_[edgeCommon] += edgeLength;
@@ -156,7 +164,6 @@ void Foam::pairPatchAgglomeration::setEdgeWeights
                 {
                     facePairWeight_.insert(edgeCommon, edgeLength);
                 }
-
                 // If the fine 'pair' faces was featured edge so it is
                 // the coarse 'pair'
                 if (fineFeaturedFaces.found(edgeCommon))
@@ -188,9 +195,9 @@ void Foam::pairPatchAgglomeration::setEdgeWeights
 
 Foam::pairPatchAgglomeration::pairPatchAgglomeration
 (
-    const polyPatch& patch,
-    const dictionary& controlDict,
-    const bool additionalWeights
+    const faceList& faces,
+    const pointField& points,
+    const dictionary& controlDict
 )
 :
     mergeLevels_
@@ -202,32 +209,61 @@ Foam::pairPatchAgglomeration::pairPatchAgglomeration
     (
         readLabel(controlDict.lookup("nFacesInCoarsestLevel"))
     ),
+    nGlobalFacesInCoarsestLevel_(labelMax),
+    //(
+    //    readLabel(controlDict.lookup("nGlobalFacesInCoarsestLevel"))
+    //),
     featureAngle_
     (
         controlDict.lookupOrDefault<scalar>("featureAngle", 0)
     ),
     nFaces_(maxLevels_),
     restrictAddressing_(maxLevels_),
-    restrictTopBottomAddressing_(identity(patch.size())),
+    restrictTopBottomAddressing_(identity(faces.size())),
     patchLevels_(maxLevels_),
-    facePairWeight_(patch.size())
+    facePairWeight_(faces.size())
 {
     // Set base fine patch
-    patchLevels_.set
-    (
-        0,
-        new bPatch
-        (
-            patch.localFaces(),
-            patch.localPoints()
-        )
-    );
+    patchLevels_.set(0, new bPatch(faces, points));
 
     // Set number of faces for the base patch
-    nFaces_[0] = patch.size();
+    nFaces_[0] = faces.size();
 
     // Set edge weights for level 0
-    setBasedEdgeWeights();
+    setLevel0EdgeWeights();
+}
+
+
+Foam::pairPatchAgglomeration::pairPatchAgglomeration
+(
+    const faceList& faces,
+    const pointField& points,
+    const label mergeLevels,
+    const label maxLevels,
+    const label nFacesInCoarsestLevel,          // local number of cells
+    const label nGlobalFacesInCoarsestLevel,    // global number of cells
+    const scalar featureAngle
+)
+:
+    mergeLevels_(mergeLevels),
+    maxLevels_(maxLevels),
+    nFacesInCoarsestLevel_(nFacesInCoarsestLevel),
+    nGlobalFacesInCoarsestLevel_(nGlobalFacesInCoarsestLevel),
+    featureAngle_(featureAngle),
+    nFaces_(maxLevels_),
+    restrictAddressing_(maxLevels_),
+    restrictTopBottomAddressing_(identity(faces.size())),
+    patchLevels_(maxLevels_),
+    facePairWeight_(faces.size())
+{
+    // Set base fine patch
+    patchLevels_.set(0, new bPatch(faces, points));
+
+    // Set number of faces for the base patch
+    nFaces_[0] = faces.size();
+
+    // Set edge weights for level 0
+    setLevel0EdgeWeights();
 }
 
 
@@ -239,7 +275,8 @@ Foam::pairPatchAgglomeration::~pairPatchAgglomeration()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::bPatch& Foam::pairPatchAgglomeration::patchLevel
+const Foam::pairPatchAgglomeration::bPatch&
+Foam::pairPatchAgglomeration::patchLevel
 (
     const label i
 ) const
@@ -277,7 +314,7 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
 
     if (fineToCoarse.size() == 0)
     {
-        return true;
+        return false;
     }
 
     if (fineToCoarse.size() != patch.size())
@@ -289,7 +326,7 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
             << abort(FatalError);
     }
 
-    const label nCoarseI =  max(fineToCoarse)+1;
+    const label nCoarseI =  max(fineToCoarse) + 1;
     List<face> patchFaces(nCoarseI);
 
 
@@ -323,7 +360,6 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
                 facePairWeight_[e1] = -1.0;
                 facePairWeight_[e2] = -1.0;
             }
-
             return false;
         }
 
@@ -355,71 +391,67 @@ void Foam::pairPatchAgglomeration::agglomerate()
 {
     label nPairLevels = 0;
     label nCreatedLevels = 1; // 0 level is the base patch
+
     label nCoarseFaces = 0;
     label nCoarseFacesOld = 0;
 
     while (nCreatedLevels < maxLevels_)
     {
         const bPatch& patch = patchLevels_[nCreatedLevels - 1];
-        tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
-        bool agglomOK = false;
 
-        do
+        // Agglomerate locally
+        tmp<labelField> tfinalAgglom;
+
+        bool createdLevel = false;
+        while (!createdLevel)
         {
-            label nCoarseFacesPrev = nCoarseFaces;
+            // Agglomerate locally using edge weights
+            // - calculates nCoarseFaces; returns fine to coarse addressing
+            tfinalAgglom = agglomerateOneLevel(nCoarseFaces, patch);
 
-            finalAgglomPtr = agglomerateOneLevel
-            (
-                nCoarseFaces,
-                patch
-            );
-
-            if (nCoarseFaces > 0 && nCoarseFaces != nCoarseFacesPrev)
+            if (nCoarseFaces == 0)
             {
-                if
-                (
-                    (
-                        agglomOK = agglomeratePatch
-                        (
-                            patch,
-                            finalAgglomPtr,
-                            nCreatedLevels
-                        )
-                    )
-                )
-                {
-                    restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
-                    mapBaseToTopAgglom(nCreatedLevels);
-                    setEdgeWeights(nCreatedLevels);
-
-                    if (nPairLevels % mergeLevels_)
-                    {
-                        combineLevels(nCreatedLevels);
-                    }
-                    else
-                    {
-                        nCreatedLevels++;
-                    }
-
-                    nPairLevels++;
-                }
+                break;
             }
             else
             {
-                agglomOK = true;
+                // Attempt to create coarse face addressing
+                // - returns true if successful; otherwise resets edge weights
+                //   and assume try again...
+                createdLevel = agglomeratePatch
+                (
+                    patch,
+                    tfinalAgglom,
+                    nCreatedLevels
+                );
+            }
+        }
+
+        if (createdLevel)
+        {
+            restrictAddressing_.set(nCreatedLevels, tfinalAgglom);
+
+            mapBaseToTopAgglom(nCreatedLevels);
+
+            setEdgeWeights(nCreatedLevels);
+
+            if (nPairLevels % mergeLevels_)
+            {
+                combineLevels(nCreatedLevels);
+            }
+            else
+            {
+                nCreatedLevels++;
             }
 
-            reduce(nCoarseFaces, sumOp<label>());
+            nPairLevels++;
 
-        } while (!agglomOK);
+            nFaces_[nCreatedLevels] = nCoarseFaces;
+        }
 
-        nFaces_[nCreatedLevels] = nCoarseFaces;
-
-        if
-        (
-            !continueAgglomerating(nCoarseFaces)
-          || (nCoarseFacesOld ==  nCoarseFaces)
-        )
+        // Check to see if we need to continue agglomerating
+        // - Note: performs parallel reductions
+        if (!continueAgglomerating(nCoarseFaces, nCoarseFacesOld))
         {
             break;
         }
@@ -496,8 +528,8 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
                     if
                     (
                         facePairWeight_[edgeCommon] > clusterMaxFaceCoeff
-                        && facePairWeight_[edgeCommon] != -1.0
-                        && coarseCellMap[faceNeig] > 0
+                     && facePairWeight_[edgeCommon] != -1.0
+                     && coarseCellMap[faceNeig] >= 0
                     )
                     {
                         clusterMatchFaceNo = faceNeig;
@@ -505,14 +537,14 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
                     }
                 }
 
-                if (clusterMatchFaceNo >= 0)
+                if (clusterMatchFaceNo > 0)
                 {
                     // Add the cell to the best cluster
                     coarseCellMap[facei] = coarseCellMap[clusterMatchFaceNo];
                 }
                 else
                 {
-                    // if not create single-cell "clusters" for each
+                    // If not create single-cell "clusters" for each
                     coarseCellMap[facei] = nCoarseFaces;
                     nCoarseFaces++;
                 }
@@ -521,15 +553,14 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
     }
 
     // Check that all faces are part of clusters,
-
     for (label facei=0; facei<nFineFaces; facei++)
     {
         if (coarseCellMap[facei] < 0)
         {
             FatalErrorInFunction
-              << " face " << facei
-            << " is not part of a cluster"
-            << exit(FatalError);
+                << " face " << facei
+                << " is not part of a cluster"
+                << exit(FatalError);
         }
     }
 

src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,16 +36,10 @@ SourceFiles
 #define pairPatchAgglomeration_H
 
 #include "mathematicalConstants.H"
-#include "polyPatch.H"
 #include "indirectPrimitivePatch.H"
 #include "List.H"
 #include "EdgeMap.H"
 
-namespace Foam
-{
-    typedef PrimitivePatch<face, List, const pointField> bPatch;
-}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -59,6 +53,8 @@ class pairPatchAgglomeration
 {
 protected:
 
+    typedef PrimitivePatch<face, List, const pointField> bPatch;
+
     // Protected data
 
         //- Number of levels to merge, 1 = don't merge, 2 = merge pairs etc.
@@ -70,6 +66,9 @@ protected:
         //- Number of faces in coarsest level
         label nFacesInCoarsestLevel_;
 
+        //- Global number of faces in coarsest level
+        label nGlobalFacesInCoarsestLevel_;
+
         //- Feature angle
         scalar featureAngle_;
 
@@ -116,13 +115,17 @@ private:
         void compactLevels(const label fineLevelIndex);
 
         //- Check the need for further agglomeration
-        bool continueAgglomerating(const label fineLevelIndex);
+        bool continueAgglomerating
+        (
+            const label nLocal,
+            const label nLocalOld
+        );
 
         //- Set edge weights
         void setEdgeWeights(const label indexLevel);
 
         //- Set base patch edge weights
-        void setBasedEdgeWeights();
+        void setLevel0EdgeWeights();
 
         //- Maps current level with base patch
         void mapBaseToTopAgglom(const label fineLevelIndex);
@@ -138,12 +141,24 @@ public:
 
     // Constructors
 
-        //- Construct given mesh and controls
+        //- Construct given faces, points and control dictionary
         pairPatchAgglomeration
         (
-            const polyPatch& patch,
-            const dictionary& controlDict,
-            const bool additionalWeights = false
+            const faceList& faces,
+            const pointField& points,
+            const dictionary& controlDict
+        );
+
+        //- Construct from components
+        pairPatchAgglomeration
+        (
+            const faceList& faces,
+            const pointField& points,
+            const label mergeLevels,
+            const label maxLevels,
+            const label nFacesInCoarsestLevel,
+            const label nGlobalFacesInCoarsestLevel,
+            const scalar featureAngle = 0
         );
 
 

src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomerationTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,7 +46,7 @@ void Foam::pairPatchAgglomeration::restrictField
             << abort(FatalError);
     }
 
-    cf = Zero;
+    cf = pTraits<Type>::zero;
 
     forAll(ff, i)
     {

src/lagrangian/basic/Make/files

@@ -3,6 +3,10 @@ particle/particleIO.C
 passiveParticle/passiveParticleCloud.C
 indexedParticle/indexedParticleCloud.C
 
+injectedParticle/injectedParticle.C
+injectedParticle/injectedParticleIO.C
+injectedParticle/injectedParticleCloud.C
+
 InteractionLists/referredWallFace/referredWallFace.C
 
 LIB = $(FOAM_LIBBIN)/liblagrangian

src/lagrangian/basic/injectedParticle/injectedParticle.C

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectedParticle.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(injectedParticle, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::injectedParticle::injectedParticle(const injectedParticle& p)
+:
+    particle(p),
+    tag_(p.tag_),
+    soi_(p.soi_),
+    d_(p.d_),
+    U_(p.U_)
+{}
+
+
+Foam::injectedParticle::injectedParticle
+(
+    const injectedParticle& p,
+    const polyMesh& mesh
+)
+:
+    particle(p, mesh),
+    tag_(p.tag_),
+    soi_(p.soi_),
+    d_(p.d_),
+    U_(p.U_)
+{}
+
+
+// ************************************************************************* //

src/lagrangian/basic/injectedParticle/injectedParticle.H

@@ -0,0 +1,250 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::injectedParticle
+
+Description
+
+SourceFiles
+    injectedParticle.C
+    injectedParticleIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef injectedParticle_H
+#define injectedParticle_H
+
+#include "particle.H"
+#include "IOstream.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class injectedParticle;
+
+// Forward declaration of friend functions
+
+Ostream& operator<<
+(
+    Ostream&,
+    const injectedParticle&
+);
+
+/*---------------------------------------------------------------------------*\
+                         Class injectedParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class injectedParticle
+:
+    public particle
+{
+public:
+
+    //- Size in bytes of the fields
+    static const std::size_t sizeofFields;
+
+
+protected:
+
+    // Protected data
+
+        // Particle properties
+
+            //- Tag
+            label tag_;
+
+            //- Start of injection [s]
+            scalar soi_;
+
+            //- Diameter [m]
+            scalar d_;
+
+            //- Velocity [m/s]
+            vector U_;
+
+
+public:
+
+    // Static data members
+
+        //- Runtime type information
+        TypeName("injectedParticle");
+
+        //- String representation of properties
+        AddToPropertyList
+        (
+            particle,
+            " tag"
+          + " soi"
+          + " d"
+          + " (Ux Uy Uz)";
+        );
+
+        //- String representation of property types
+        AddToPropertyTypes
+        (
+            particle,
+            "{label"
+          + " scalar"
+          + " scalar"
+          + " vector}"
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        inline injectedParticle
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label tag,
+            const scalar soi,
+            const scalar d,
+            const vector& U
+        );
+
+        //- Construct from Istream
+        injectedParticle
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct as a copy
+        injectedParticle(const injectedParticle& p);
+
+        //- Construct as a copy
+        injectedParticle(const injectedParticle& p, const polyMesh& mesh);
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new injectedParticle(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone(const polyMesh& mesh) const
+        {
+            return autoPtr<particle>(new injectedParticle(*this, mesh));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<injectedParticle> operator()(Istream& is) const
+            {
+                return autoPtr<injectedParticle>
+                (
+                    new injectedParticle(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the tag
+            inline label tag() const;
+
+            //- Return const access to the start of injection
+            inline scalar soi() const;
+
+            //- Return const access to diameter
+            inline scalar d() const;
+
+            //- Return const access to velocity
+            inline const vector& U() const;
+
+
+        // Edit
+
+            //- Return the tag
+            inline label& tag();
+
+            //- Return the start of injection
+            inline scalar& soi();
+
+            //- Return access to diameter
+            inline scalar& d();
+
+            //- Return access to velocity
+            inline vector& U();
+
+
+        // I-O
+
+            //- Read
+            static void readFields(Cloud<injectedParticle>& c);
+
+            //- Write
+            static void writeFields(const Cloud<injectedParticle>& c);
+
+            //- Write particle fields as objects into the obr registry
+            static void writeObjects
+            (
+                const Cloud<injectedParticle>& c,
+                objectRegistry& obr
+            );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const injectedParticle&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "injectedParticleI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/basic/injectedParticle/injectedParticleCloud.C

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectedParticleCloud.H"
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(Cloud<injectedParticle>, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::injectedParticleCloud::injectedParticleCloud
+(
+    const polyMesh& mesh,
+    const word& cloudName,
+    bool readFields
+)
+:
+    Cloud<injectedParticle>(mesh, cloudName, false)
+{
+    if (readFields)
+    {
+        injectedParticle::readFields(*this);
+    }
+}
+
+
+Foam::injectedParticleCloud::injectedParticleCloud
+(
+    const injectedParticleCloud& c,
+    const word& cloudName
+)
+:
+    Cloud<injectedParticle>(c.pMesh(), cloudName, c)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::injectedParticleCloud::~injectedParticleCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::injectedParticleCloud::writeObjects(objectRegistry& obr) const
+{
+    injectedParticle::writeObjects(*this, obr);
+}
+
+
+// ************************************************************************* //

src/lagrangian/basic/injectedParticle/injectedParticleCloud.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::injectedParticleCloud
+
+Group
+
+Description
+
+SourceFiles
+    injectedParticleCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef injectedParticleCloud_H
+#define injectedParticleCloud_H
+
+#include "Cloud.H"
+#include "injectedParticle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+/*---------------------------------------------------------------------------*\
+                   Class injectedParticleCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class injectedParticleCloud
+:
+    public Cloud<injectedParticle>
+{
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        injectedParticleCloud(const injectedParticleCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const injectedParticleCloud&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh and cloud name
+        injectedParticleCloud
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        injectedParticleCloud(const injectedParticleCloud& c, const word& name);
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<injectedParticleCloud> clone(const word& name) const
+        {
+            return autoPtr<injectedParticleCloud>
+            (
+                new injectedParticleCloud(*this, name)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~injectedParticleCloud();
+
+
+    // Member Functions
+
+        // I-O
+
+            //- Write particle fields as objects into the obr registry
+            virtual void writeObjects(objectRegistry& obr) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/basic/injectedParticle/injectedParticleI.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::injectedParticle::injectedParticle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label tag,
+    const scalar soi,
+    const scalar d,
+    const vector& U
+)
+:
+    particle(mesh, position, -1, false),
+    tag_(tag),
+    soi_(soi),
+    d_(d),
+    U_(U)
+{}
+
+
+inline Foam::label Foam::injectedParticle::tag() const
+{
+    return tag_;
+}
+
+
+inline Foam::scalar Foam::injectedParticle::soi() const
+{
+    return soi_;
+}
+
+
+inline Foam::scalar Foam::injectedParticle::d() const
+{
+    return d_;
+}
+
+
+inline const Foam::vector& Foam::injectedParticle::U() const
+{
+    return U_;
+}
+
+
+inline Foam::label& Foam::injectedParticle::tag()
+{
+    return tag_;
+}
+
+
+inline Foam::scalar& Foam::injectedParticle::soi()
+{
+    return soi_;
+}
+
+
+inline Foam::scalar& Foam::injectedParticle::d()
+{
+    return d_;
+}
+
+
+inline Foam::vector& Foam::injectedParticle::U()
+{
+    return U_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/basic/injectedParticle/injectedParticleIO.C

@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "injectedParticle.H"
+#include "IOstreams.H"
+#include "IOField.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::string Foam::injectedParticle::propertyList_ =
+    Foam::injectedParticle::propertyList();
+
+Foam::string Foam::injectedParticle::propertyTypes_ =
+    Foam::injectedParticle::propertyTypes();
+
+const std::size_t Foam::injectedParticle::sizeofFields
+(
+    sizeof(label) + sizeof(scalar) + sizeof(scalar) + sizeof(vector)
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::injectedParticle::injectedParticle
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    particle(mesh, is, readFields),
+    tag_(-1),
+    soi_(0.0),
+    d_(0.0),
+    U_(Zero)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            tag_ = readLabel(is);
+            soi_ = readScalar(is);
+            d_ = readScalar(is);
+            is >> U_;
+        }
+        else
+        {
+            is.read(reinterpret_cast<char*>(&soi_), sizeofFields);
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "injectedParticle::injectedParticle"
+        "(const polyMesh&, Istream&, bool)"
+    );
+}
+
+
+void Foam::injectedParticle::readFields(Cloud<injectedParticle>& c)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    particle::readFields(c);
+
+    IOField<label> tag(c.fieldIOobject("tag", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, tag);
+
+    IOField<scalar> soi(c.fieldIOobject("soi", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, soi);
+
+    IOField<scalar> d(c.fieldIOobject("d", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, d);
+
+    IOField<vector> U(c.fieldIOobject("U", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, U);
+
+    label i = 0;
+
+    forAllIter(Cloud<injectedParticle>, c, iter)
+    {
+        injectedParticle& p = iter();
+
+        p.tag_ = tag[i];
+        p.soi_ = soi[i];
+        p.d_ = d[i];
+        p.U_ = U[i];
+
+        i++;
+    }
+}
+
+
+void Foam::injectedParticle::writeFields(const Cloud<injectedParticle>& c)
+{
+    particle::writeFields(c);
+
+    label np =  c.size();
+
+    IOField<label> tag(c.fieldIOobject("tag", IOobject::NO_READ), np);
+    IOField<scalar> soi(c.fieldIOobject("soi", IOobject::NO_READ), np);
+    IOField<scalar> d(c.fieldIOobject("d", IOobject::NO_READ), np);
+    IOField<vector> U(c.fieldIOobject("U", IOobject::NO_READ), np);
+
+    label i = 0;
+
+    forAllConstIter(Cloud<injectedParticle>, c, iter)
+    {
+        const injectedParticle& p = iter();
+
+        tag[i] = p.tag();
+        soi[i] = p.soi();
+        d[i] = p.d();
+        U[i] = p.U();
+
+        i++;
+    }
+
+    tag.write();
+    soi.write();
+    d.write();
+    U.write();
+}
+
+
+void Foam::injectedParticle::writeObjects
+(
+    const Cloud<injectedParticle>& c,
+    objectRegistry& obr
+)
+{
+    particle::writeObjects(c, obr);
+
+    label np = c.size();
+
+    IOField<label>& tag(cloud::createIOField<label>("tag", np, obr));
+    IOField<scalar>& soi(cloud::createIOField<scalar>("soi", np, obr));
+    IOField<scalar>& d(cloud::createIOField<scalar>("d", np, obr));
+    IOField<vector>& U(cloud::createIOField<vector>("U", np, obr));
+
+    label i = 0;
+
+    forAllConstIter(Cloud<injectedParticle>, c, iter)
+    {
+        const injectedParticle& p = iter();
+
+        tag[i] = p.tag();
+        soi[i] = p.soi();
+        d[i] = p.d();
+        U[i] = p.U();
+
+        i++;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const injectedParticle& p
+)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const particle&>(p)
+            << token::SPACE << p.tag()
+            << token::SPACE << p.soi()
+            << token::SPACE << p.d()
+            << token::SPACE << p.U();
+    }
+    else
+    {
+        os  << static_cast<const particle&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.tag_),
+            injectedParticle::sizeofFields
+        );
+    }
+
+    // Check state of Ostream
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, const injectedParticle&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/basic/particle/particle.H

@@ -574,6 +574,10 @@ public:
         template<class CloudType>
         static void writeFields(const CloudType& c);
 
+        //- Write particle fields as objects into the obr registry
+        template<class CloudType>
+        static void writeObjects(const CloudType& c, objectRegistry& obr);
+
         //- Write the particle position and cell
         void writePosition(Ostream&) const;
 

src/lagrangian/basic/particle/particleTemplates.C

@@ -162,7 +162,7 @@ void Foam::particle::writeFields(const CloudType& c)
     IOPosition<CloudType> ioP(c);
     ioP.write();
 
-    label np =  c.size();
+    label np = c.size();
 
     IOField<label> origProc
     (
@@ -184,6 +184,29 @@ void Foam::particle::writeFields(const CloudType& c)
 }
 
 
+template<class CloudType>
+void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
+{
+    label np = c.size();
+
+    IOField<vector>& position
+    (
+        cloud::createIOField<vector>("position", np, obr)
+    );
+    IOField<label>& origProc(cloud::createIOField<label>("origProc", np, obr));
+    IOField<label>& origId(cloud::createIOField<label>("origId", np, obr));
+
+    label i = 0;
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        position[i] = iter().position_;
+        origProc[i] = iter().origProc_;
+        origId[i] = iter().origId_;
+        i++;
+    }
+}
+
+
 template<class TrackData>
 Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)
 {

src/lagrangian/distributionModels/Make/files

@@ -1,6 +1,7 @@
 distributionModel/distributionModel.C
 distributionModel/distributionModelNew.C
 
+binned/binned.C
 exponential/exponential.C
 fixedValue/fixedValue.C
 general/general.C

src/lagrangian/distributionModels/RosinRammler/RosinRammler.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::RosinRammler
+    Foam::distributionModels::RosinRammler
 
 Description
     Rosin-Rammler distributionModel
@@ -39,8 +39,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef RosinRammler_H
-#define RosinRammler_H
+#ifndef distributionModels_RosinRammler_H
+#define distributionModels_RosinRammler_H
 
 #include "distributionModel.H"
 

src/lagrangian/distributionModels/binned/binned.C

@@ -0,0 +1,282 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "binned.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace distributionModels
+    {
+        defineTypeNameAndDebug(binned, 0);
+        addToRunTimeSelectionTable(distributionModel, binned, dictionary);
+    }
+}
+
+
+const char* Foam::distributionModels::binned::header =
+    "(bin probability)";
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::distributionModels::binned::initialise()
+{
+    const label nSample(xy_.size());
+    forAll(xy_, bini)
+    {
+        xy_[bini][1] /= scalar(nSample);
+    }
+
+    // Convert values to integral values
+    for (label bini = 1; bini < nSample; ++bini)
+    {
+        xy_[bini][1] += xy_[bini - 1][1];
+    }
+
+    // Calculate the mean value
+    label bini = 0;
+    forAll(xy_, i)
+    {
+        if (xy_[i][1] > 0.5)
+        {
+            bini = i;
+            break;
+        }
+    }
+
+    meanValue_ = xy_[bini][1];
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributionModels::binned::binned
+(
+    const dictionary& dict,
+    cachedRandom& rndGen
+)
+:
+    distributionModel(typeName, dict, rndGen),
+    xy_(distributionModelDict_.lookup("distribution")),
+    meanValue_(0)
+{
+    if (maxValue() < minValue())
+    {
+        FatalErrorInFunction
+            << "Maximum value is smaller than the minimum value:"
+            << "    maxValue = " << maxValue()
+            << ", minValue = " << minValue()
+            << exit(FatalError);
+    }
+
+    initialise();
+}
+
+
+Foam::distributionModels::binned::binned
+(
+    const UList<scalar>& sampleData,
+    const scalar binWidth,
+    cachedRandom& rndGen
+)
+:
+    distributionModel(typeName, dictionary::null, rndGen),
+    xy_(),
+    meanValue_(0)
+{
+    scalar minValue = GREAT;
+    scalar maxValue = -GREAT;
+    forAll(sampleData, i)
+    {
+        minValue = min(minValue, sampleData[i]);
+        maxValue = max(maxValue, sampleData[i]);
+    }
+
+    const label bin0 = floor(minValue/binWidth);
+    const label bin1 = ceil(maxValue/binWidth);
+    const label nBin = bin1 - bin0;
+
+    if (nBin == 0)
+    {
+        WarningInFunction
+            << "Data cannot be binned - zero bins generated" << nl
+            << "   Bin width   : " << binWidth << nl
+            << "   Sample data : " << sampleData
+            << endl;
+
+        return;
+    }
+
+    // Populate bin boundaries and initialise occurrences
+    xy_.setSize(nBin);
+    forAll(xy_, bini)
+    {
+        xy_[bini][0] = (bin0 + bini)*binWidth;
+        xy_[bini][1] = 0;
+    }
+
+    // Bin the data
+    forAll(sampleData, i)
+    {
+        // Choose the nearest bin
+        label bini = floor(sampleData[i]/binWidth) - bin0;
+        label binii = min(bini + 1, nBin - 1);
+
+        scalar d1 = mag(sampleData[i] - xy_[bini][0]);
+        scalar d2 = mag(xy_[binii][0] - sampleData[i]);
+
+        if (d1 < d2)
+        {
+            xy_[bini][1]++;
+        }
+        else
+        {
+            xy_[binii][1]++;
+        }
+    }
+
+    initialise();
+}
+
+
+Foam::distributionModels::binned::binned(const binned& p)
+:
+    distributionModel(p),
+    xy_(p.xy_),
+    meanValue_(p.meanValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::distributionModels::binned::~binned()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::distributionModels::binned::sample() const
+{
+    scalar y = rndGen_.sample01<scalar>();
+
+    for (label i = 0; i < xy_.size() - 1; ++i)
+    {
+        if (xy_[i][1] > y)
+        {
+            scalar d1 = y - xy_[i][1];
+            scalar d2 = xy_[i+1][1] - y;
+            if (d1 < d2)
+            {
+                return xy_[i][0];
+            }
+            else
+            {
+                return xy_[i+1][0];
+            }
+        }
+    }
+
+    return xy_.last()[0];
+}
+
+
+Foam::scalar Foam::distributionModels::binned::minValue() const
+{
+    return xy_.first()[0];
+}
+
+
+Foam::scalar Foam::distributionModels::binned::maxValue() const
+{
+    return xy_.last()[0];
+}
+
+
+Foam::scalar Foam::distributionModels::binned::meanValue() const
+{
+    return meanValue_;
+}
+
+
+void Foam::distributionModels::binned::readData(Istream& is)
+{
+//    distributionModel::readData(is);
+    is  >> xy_;
+}
+
+
+void Foam::distributionModels::binned::writeData(Ostream& os) const
+{
+//    distributionModel::writeData(os);
+    os  << xy_ ;
+}
+
+
+Foam::dictionary Foam::distributionModels::binned::writeDict
+(
+    const word& dictName
+) const
+{
+//    dictionary dict = distributionModel::writeDict(dictName);
+    dictionary dict(dictName);
+    dict.add("distribution", xy_);
+
+    return dict;
+}
+
+
+void Foam::distributionModels::binned::readDict(const dictionary& dict)
+{
+//    distributionModel::readDict(dict);
+    dict.lookup("distribution") >> xy_;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const distributionModels::binned& b
+)
+{
+    os.check(FUNCTION_NAME);
+
+    b.writeData(os);
+    return os;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, distributionModels::binned& b)
+{
+    is.check(FUNCTION_NAME);
+
+    b.readData(is);
+    return is;
+}
+
+
+// ************************************************************************* //

src/lagrangian/distributionModels/binned/binned.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::distributionModels::binned
+
+Description
+    A binned distribution model where the random sample is taken from a
+    discrete set of bins
+
+SourceFiles
+    binned.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributionModels_binned_H
+#define distributionModels_binned_H
+
+#include "distributionModel.H"
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+namespace distributionModels
+{
+    class binned;
+}
+
+// Forward declaration of friend functions and operators
+Istream& operator>>(Istream&, distributionModels::binned&);
+Ostream& operator<<(Ostream&, const distributionModels::binned&);
+
+namespace distributionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class binned Declaration
+\*---------------------------------------------------------------------------*/
+
+class binned
+:
+    public distributionModel
+{
+    // Private data
+
+        typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
+
+        // List of (bin probability)
+        List<pair> xy_;
+
+        //- Distribution mean value
+        scalar meanValue_;
+
+
+    // Private member functions
+
+        //- Initialise the distribution parameters
+        void initialise();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("binned");
+
+    static const char* header;
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        binned(const dictionary& dict, cachedRandom& rndGen);
+
+        //- Construct from components
+        binned
+        (
+            const UList<scalar>& sampleData,
+            const scalar binWidth,
+            cachedRandom& rndGen
+        );
+
+        //- Construct copy
+        binned(const binned& p);
+
+        //- Construct and return a clone
+        virtual autoPtr<distributionModel> clone() const
+        {
+            return autoPtr<distributionModel>(new binned(*this));
+        }
+
+
+    //- Destructor
+    virtual ~binned();
+
+
+    // Member Functions
+
+        //- Sample the distributionModel
+        virtual scalar sample() const;
+
+        //- Return the minimum value
+        virtual scalar minValue() const;
+
+        //- Return the maximum value
+        virtual scalar maxValue() const;
+
+        //- Return the mean value
+        virtual scalar meanValue() const;
+
+        //- Write data to stream
+        virtual void writeData(Ostream& os) const;
+
+        //- Read data from stream
+        virtual void readData(Istream& os);
+
+        //- Write data in dictionary format
+        virtual dictionary writeDict(const word& dictName) const;
+
+        //- Read data from dictionary
+        virtual void readDict(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace distributionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/distributionModels/distributionModel/distributionModel.C

@@ -69,7 +69,7 @@ Foam::distributionModels::distributionModel::distributionModel
     cachedRandom& rndGen
 )
 :
-    distributionModelDict_(dict.subDict(name + "Distribution")),
+    distributionModelDict_(dict),
     rndGen_(rndGen)
 {}
 

src/lagrangian/distributionModels/distributionModel/distributionModelNew.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,7 +50,10 @@ Foam::distributionModels::distributionModel::New
             << exit(FatalError);
     }
 
-    return autoPtr<distributionModel>(cstrIter()(dict, rndGen));
+    const dictionary distributionDict =
+        dict.subOrEmptyDict(modelType & "Distribution");
+
+    return autoPtr<distributionModel>(cstrIter()(distributionDict, rndGen));
 }
 
 

src/lagrangian/distributionModels/exponential/exponential.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::exponential
+    Foam::distributionModels::exponential
 
 Description
     exponential distribution model
@@ -32,8 +32,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef exponential_H
-#define exponential_H
+#ifndef distributionModels_exponential_H
+#define distributionModels_exponential_H
 
 #include "distributionModel.H"
 

src/lagrangian/distributionModels/fixedValue/fixedValue.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::fixedValue
+    Foam::distributionModels::fixedValue
 
 Description
     Returns a fixed value
@@ -32,8 +32,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef distributionModelFixedValue_H
-#define distributionModelFixedValue_H
+#ifndef distributionModels_fixedValue_H
+#define distributionModels_fixedValue_H
 
 #include "distributionModel.H"
 
@@ -44,7 +44,7 @@ namespace Foam
 namespace distributionModels
 {
 /*---------------------------------------------------------------------------*\
-                       Class fixedValue Declaration
+                         Class fixedValue Declaration
 \*---------------------------------------------------------------------------*/
 
 class fixedValue

src/lagrangian/distributionModels/general/general.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,30 +37,19 @@ namespace distributionModels
 }
 }
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::distributionModels::general::general
-(
-    const dictionary& dict,
-    cachedRandom& rndGen
-)
-:
-    distributionModel(typeName, dict, rndGen),
-    xy_(distributionModelDict_.lookup("distribution")),
-    nEntries_(xy_.size()),
-    minValue_(xy_[0][0]),
-    maxValue_(xy_[nEntries_-1][0]),
-    meanValue_(0.0),
-    integral_(nEntries_)
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::distributionModels::general::initialise()
 {
-    check();
+    const label nEntries = xy_.size();
 
-    // normalize the cumulative distributionModel
+    integral_.setSize(nEntries);
 
+    // Normalize the cumulative distribution
     integral_[0] = 0.0;
-    for (label i=1; i<nEntries_; i++)
+    for (label i = 1; i < nEntries; i++)
     {
-
         scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
         scalar d = xy_[i-1][1] - k*xy_[i-1][0];
         scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
@@ -72,14 +61,87 @@ Foam::distributionModels::general::general
 
     scalar sumArea = integral_.last();
 
-    meanValue_ = sumArea/(maxValue_ - minValue_);
+    meanValue_ = sumArea/(maxValue() - minValue());
 
-    for (label i=0; i<nEntries_; i++)
+    for (label i=0; i < nEntries; i++)
     {
         xy_[i][1] /= sumArea;
         integral_[i] /= sumArea;
     }
+}
 
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributionModels::general::general
+(
+    const dictionary& dict,
+    cachedRandom& rndGen
+)
+:
+    distributionModel(typeName, dict, rndGen),
+    xy_(distributionModelDict_.lookup("distribution")),
+    meanValue_(0.0),
+    integral_()
+{
+    check();
+
+    initialise();
+}
+
+
+Foam::distributionModels::general::general
+(
+    const UList<scalar>& sampleData,
+    const scalar binWidth,
+    cachedRandom& rndGen
+)
+:
+    distributionModel(typeName, dictionary::null, rndGen),
+    xy_(),
+    meanValue_(0.0),
+    integral_()
+{
+    scalar minValue = GREAT;
+    scalar maxValue = -GREAT;
+    forAll(sampleData, i)
+    {
+        minValue = min(minValue, sampleData[i]);
+        maxValue = max(maxValue, sampleData[i]);
+    }
+
+    label bin0 = floor(minValue/binWidth);
+    label bin1 = ceil(maxValue/binWidth);
+    label nEntries = bin1 - bin0;
+
+    if (nEntries == 0)
+    {
+        WarningInFunction
+            << "Data cannot be binned - zero bins generated" << nl
+            << "   Bin width   : " << binWidth << nl
+            << "   Sample data : " << sampleData
+            << endl;
+
+        return;
+    }
+
+    xy_.setSize(nEntries);
+
+    // Populate bin boundaries and initialise occurrences
+    for (label bini = 0; bini < nEntries; ++bini)
+    {
+        xy_[bini][0] = (bin0 + bini)*binWidth;
+        xy_[bini][1] = 0;
+    }
+
+    // Populate occurrences
+    forAll(sampleData, i)
+    {
+        label bini = floor(sampleData[i]/binWidth) - bin0;
+        xy_[bini][1]++;
+    }
+
+    initialise();
 }
 
 
@@ -87,9 +149,6 @@ Foam::distributionModels::general::general(const general& p)
 :
     distributionModel(p),
     xy_(p.xy_),
-    nEntries_(p.nEntries_),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_),
     integral_(p.integral_)
 {}
 
@@ -106,8 +165,8 @@ Foam::scalar Foam::distributionModels::general::sample() const
 {
     scalar y = rndGen_.sample01<scalar>();
 
-    // find the interval where y is in the table
-    label n=1;
+    // Find the interval where y is in the table
+    label n = 1;
     while (integral_[n] <= y)
     {
         n++;
@@ -119,7 +178,7 @@ Foam::scalar Foam::distributionModels::general::sample() const
     scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
     scalar x = 0.0;
 
-    // if k is small it is a linear equation, otherwise it is of second order
+    // If k is small it is a linear equation, otherwise it is of second order
     if (mag(k) > SMALL)
     {
         scalar p = 2.0*d/k;
@@ -148,13 +207,13 @@ Foam::scalar Foam::distributionModels::general::sample() const
 
 Foam::scalar Foam::distributionModels::general::minValue() const
 {
-    return minValue_;
+    return xy_.first()[0];
 }
 
 
 Foam::scalar Foam::distributionModels::general::maxValue() const
 {
-    return maxValue_;
+    return xy_.last()[0];
 }
 
 
@@ -164,4 +223,62 @@ Foam::scalar Foam::distributionModels::general::meanValue() const
 }
 
 
+void Foam::distributionModels::general::readData(Istream& is)
+{
+//    distributionModel::readData(is);
+    is  >> xy_;
+    initialise();
+}
+
+
+void Foam::distributionModels::general::writeData(Ostream& os) const
+{
+//    distributionModel::writeData(os);
+    os  << xy_;
+}
+
+
+Foam::dictionary Foam::distributionModels::general::writeDict
+(
+    const word& dictName
+) const
+{
+//    dictionary dict = distributionModel::writeDict(dictName);
+    dictionary dict(dictName);
+    dict.add("distribution", xy_);
+
+    return dict;
+}
+
+
+void Foam::distributionModels::general::readDict(const dictionary& dict)
+{
+//    distributionModel::readDict(dict);
+    dict.lookup("distribution") >> xy_;
+    initialise();
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const distributionModels::general& b
+)
+{
+    os.check(FUNCTION_NAME);
+
+    b.writeData(os);
+    return os;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, distributionModels::general& b)
+{
+    is.check(FUNCTION_NAME);
+
+    b.readData(is);
+    return is;
+}
+
+
 // ************************************************************************* //

src/lagrangian/distributionModels/general/general.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::general
+    Foam::distributionModels::general
 
 Description
     general distribution model
@@ -32,8 +32,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef general_H
-#define general_H
+#ifndef distributionModels_general_H
+#define distributionModels_general_H
 
 #include "distributionModel.H"
 #include "Vector.H"
@@ -43,6 +43,20 @@ SourceFiles
 
 namespace Foam
 {
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+namespace distributionModels
+{
+    class general;
+}
+
+// Forward declaration of friend functions and operators
+Istream& operator>>(Istream&, distributionModels::general&);
+Ostream& operator<<(Ostream&, const distributionModels::general&);
+
 namespace distributionModels
 {
 
@@ -58,19 +72,20 @@ class general
 
         typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
 
+        // List of (bin probability)
         List<pair> xy_;
 
-        label nEntries_;
-
-        //- Min and max values of the distribution
-        scalar minValue_;
-        scalar maxValue_;
-
         scalar meanValue_;
 
         List<scalar> integral_;
 
 
+    // Private member functions
+
+        //- Initialise the distribution parameters
+        void initialise();
+
+
 public:
 
     //- Runtime type information
@@ -82,6 +97,14 @@ public:
         //- Construct from components
         general(const dictionary& dict, cachedRandom& rndGen);
 
+        //- Construct from components
+        general
+        (
+            const UList<scalar>& sampleData,
+            const scalar binWidth,
+            cachedRandom& rndGen
+        );
+
         //- Construct copy
         general(const general& p);
 
@@ -109,6 +132,18 @@ public:
 
         //- Return the mean value
         virtual scalar meanValue() const;
+
+        //- Write data to stream
+        virtual void writeData(Ostream& os) const;
+
+        //- Read data from stream
+        virtual void readData(Istream& os);
+
+        //- Write data in dictionary format
+        virtual dictionary writeDict(const word& dictName) const;
+
+        //- Read data from dictionary
+        virtual void readDict(const dictionary& dict);
 };
 
 

src/lagrangian/distributionModels/multiNormal/multiNormal.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::multiNormal
+    Foam::distributionModels::multiNormal
 
 Description
     A multiNormal distribution model
@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef multiNormal_H
-#define multiNormal_H
+#ifndef distributionModels_multiNormal_H
+#define distributionModels_multiNormal_H
 
 #include "distributionModel.H"
 
@@ -50,7 +50,7 @@ namespace distributionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class multiNormal Declaration
+                         Class multiNormal Declaration
 \*---------------------------------------------------------------------------*/
 
 class multiNormal

src/lagrangian/distributionModels/normal/normal.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::normal
+    Foam::distributionModels::normal
 
 Description
     A normal distribution model
@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef normal_H
-#define normal_H
+#ifndef distributionModels_normal_H
+#define distributionModels_normal_H
 
 #include "distributionModel.H"
 

src/lagrangian/distributionModels/uniform/uniform.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::uniform
+    Foam::distributionModels::uniform
 
 Description
     Uniform/equally-weighted distribution model
@@ -32,8 +32,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef uniform_H
-#define uniform_H
+#ifndef distributionModels_uniform_H
+#define distributionModels_uniform_H
 
 #include "distributionModel.H"
 

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

@@ -895,4 +895,11 @@ void Foam::KinematicCloud<CloudType>::info()
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::writeObjects(objectRegistry& obr) const
+{
+    parcelType::writeObjects(*this, obr);
+}
+
+
 // ************************************************************************* //

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

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -612,6 +612,9 @@ public:
 
             //- Print cloud information
             void info();
+
+            //- Write particle fields as objects into the obr registry
+            virtual void writeObjects(objectRegistry& obr) const;
 };
 
 

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

@@ -362,4 +362,11 @@ void Foam::ReactingCloud<CloudType>::writeFields() const
 }
 
 
+template<class CloudType>
+void Foam::ReactingCloud<CloudType>::writeObjects(objectRegistry& obr) const
+{
+    CloudType::particleType::writeObjects(*this, this->composition(), obr);
+}
+
+
 // ************************************************************************* //

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

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -313,6 +313,9 @@ public:
 
             //- Write the field data for the cloud
             virtual void writeFields() const;
+
+            //- Write particle fields as objects into the obr registry
+            virtual void writeObjects(objectRegistry& obr) const;
 };
 
 

src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H

@@ -63,11 +63,15 @@ public:
     //- Runtime type information
     TypeName("kinematicCloud");
 
+
     // Constructors
 
         //- Null constructor
         kinematicCloud();
 
+    //- Destructor
+    virtual ~kinematicCloud();
+
 
     // Member functions
 
@@ -92,26 +96,22 @@ public:
             virtual scalar Dmax() const = 0;
 
 
-            // Fields
+        // Fields
 
-                //- Volume swept rate of parcels per cell
-                virtual const tmp<volScalarField> vDotSweep() const = 0;
+            //- Volume swept rate of parcels per cell
+            virtual const tmp<volScalarField> vDotSweep() const = 0;
 
-                //- Return the particle volume fraction field
-                //  Note: for particles belonging to this cloud only
-                virtual const tmp<volScalarField> theta() const = 0;
+            //- Return the particle volume fraction field
+            //  Note: for particles belonging to this cloud only
+            virtual const tmp<volScalarField> theta() const = 0;
 
-                //- Return the particle mass fraction field
-                //  Note: for particles belonging to this cloud only
-                virtual const tmp<volScalarField> alpha() const = 0;
+            //- Return the particle mass fraction field
+            //  Note: for particles belonging to this cloud only
+            virtual const tmp<volScalarField> alpha() const = 0;
 
-                //- Return the particle effective density field
-                //  Note: for particles belonging to this cloud only
-                virtual const tmp<volScalarField> rhoEff() const = 0;
-
-
-    //- Destructor
-    virtual ~kinematicCloud();
+            //- Return the particle effective density field
+            //  Note: for particles belonging to this cloud only
+            virtual const tmp<volScalarField> rhoEff() const = 0;
 };
 
 

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

@@ -307,6 +307,10 @@ public:
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType>
+            static void writeObjects(const CloudType& c, objectRegistry& obr);
+
 
     // Ostream Operator
 

src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcelIO.C

@@ -274,6 +274,39 @@ void Foam::CollidingParcel<ParcelType>::writeFields(const CloudType& c)
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::CollidingParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    label np = c.size();
+
+    IOField<vector>& f(cloud::createIOField<vector>("f", np, obr));
+    IOField<vector>& angularMomentum
+    (
+        cloud::createIOField<vector>("angularMomentum", np, obr)
+    );
+    IOField<vector>& torque(cloud::createIOField<vector>("torque", np, obr));
+
+    label i = 0;
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        const CollidingParcel<ParcelType>& p = iter();
+
+        f[i] = p.f();
+        angularMomentum[i] = p.angularMomentum();
+        torque[i] = p.torque();
+
+        i++;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

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

@@ -660,6 +660,10 @@ public:
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType>
+            static void writeObjects(const CloudType& c, objectRegistry& obr);
+
 
     // Ostream Operator
 

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

@@ -171,7 +171,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
 {
     ParcelType::writeFields(c);
 
-    label np =  c.size();
+    label np = c.size();
 
     IOField<label> active(c.fieldIOobject("active", IOobject::NO_READ), np);
     IOField<label> typeId(c.fieldIOobject("typeId", IOobject::NO_READ), np);
@@ -221,6 +221,54 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::KinematicParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    label np = c.size();
+
+    IOField<label>& active(cloud::createIOField<label>("active", np, obr));
+    IOField<label>& typeId(cloud::createIOField<label>("typeId", np, obr));
+    IOField<scalar>& nParticle
+    (
+        cloud::createIOField<scalar>("nParticle", np, obr)
+    );
+    IOField<scalar>& d(cloud::createIOField<scalar>("d", np, obr));
+    IOField<scalar>& dTarget(cloud::createIOField<scalar>("dTarget", np, obr));
+    IOField<vector>& U(cloud::createIOField<vector>("U", np, obr));
+    IOField<scalar>& rho(cloud::createIOField<scalar>("rho", np, obr));
+    IOField<scalar>& age(cloud::createIOField<scalar>("age", np, obr));
+    IOField<scalar>& tTurb(cloud::createIOField<scalar>("tTurb", np, obr));
+    IOField<vector>& UTurb(cloud::createIOField<vector>("UTurb", np, obr));
+
+    label i = 0;
+
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        const KinematicParcel<ParcelType>& p = iter();
+
+        active[i] = p.active();
+        typeId[i] = p.typeId();
+        nParticle[i] = p.nParticle();
+        d[i] = p.d();
+        dTarget[i] = p.dTarget();
+        U[i] = p.U();
+        rho[i] = p.rho();
+        age[i] = p.age();
+        tTurb[i] = p.tTurb();
+        UTurb[i] = p.UTurb();
+
+        i++;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

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

@@ -297,6 +297,10 @@ public:
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType>
+            static void writeObjects(const CloudType& c, objectRegistry& obr);
+
 
         // Ostream operator
 

src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelIO.C

@@ -110,7 +110,7 @@ void Foam::MPPICParcel<ParcelType>::writeFields(const CloudType& c)
 {
     ParcelType::writeFields(c);
 
-    label np =  c.size();
+    label np = c.size();
 
     IOField<vector>
         UCorrect(c.fieldIOobject("UCorrect", IOobject::NO_READ), np);
@@ -130,6 +130,34 @@ void Foam::MPPICParcel<ParcelType>::writeFields(const CloudType& c)
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::MPPICParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    label np = c.size();
+
+    IOField<vector>&
+        UCorrect(cloud::createIOField<vector>("UCorrect", np, obr));
+
+    label i = 0;
+
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        const MPPICParcel<ParcelType>& p = iter();
+
+        UCorrect[i] = p.UCorrect();
+
+        i++;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

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

@@ -433,7 +433,7 @@ public:
 
         // I-O
 
-            //- Read
+            //- Read - composition supplied
             template<class CloudType, class CompositionType>
             static void readFields
             (
@@ -445,7 +445,7 @@ public:
             template<class CloudType>
             static void readFields(CloudType& c);
 
-            //- Write
+            //- Write - composition supplied
             template<class CloudType, class CompositionType>
             static void writeFields
             (
@@ -453,10 +453,28 @@ public:
                 const CompositionType& compModel
             );
 
-            //- Read - composition supplied
+            //- Read - no composition
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
+            //- Write particle fields as objects into the obr registry
+            //  - no composition
+            template<class CloudType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                objectRegistry& obr
+            );
+
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType, class CompositionType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                const CompositionType& compModel,
+                objectRegistry& obr
+            );
+
 
     // Ostream Operator
 

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

@@ -327,6 +327,108 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::writeFields
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::ReactingMultiphaseParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ReactingMultiphaseParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    const CompositionType& compModel,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    label np = c.size();
+
+    // Write the composition fractions
+    if (np > 0)
+    {
+        const wordList& stateLabels = compModel.stateLabels();
+
+        const label idGas = compModel.idGas();
+        const wordList& gasNames = compModel.componentNames(idGas);
+        forAll(gasNames, j)
+        {
+            const word fieldName = "Y" + gasNames[j] + stateLabels[idGas];
+            IOField<scalar>& YGas
+            (
+                cloud::createIOField<scalar>(fieldName, np, obr)
+            );
+
+            label i = 0;
+            forAllConstIter
+            (
+                typename Cloud<ReactingMultiphaseParcel<ParcelType>>,
+                c,
+                iter
+            )
+            {
+                const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
+                YGas[i++] = p0.YGas()[j]*p0.Y()[GAS];
+            }
+        }
+
+        const label idLiquid = compModel.idLiquid();
+        const wordList& liquidNames = compModel.componentNames(idLiquid);
+        forAll(liquidNames, j)
+        {
+            const word fieldName = "Y" + liquidNames[j] + stateLabels[idLiquid];
+            IOField<scalar>& YLiquid
+            (
+                cloud::createIOField<scalar>(fieldName, np, obr)
+            );
+
+            label i = 0;
+            forAllConstIter
+            (
+                typename Cloud<ReactingMultiphaseParcel<ParcelType>>,
+                c,
+                iter
+            )
+            {
+                const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
+                YLiquid[i++] = p0.YLiquid()[j]*p0.Y()[LIQ];
+            }
+        }
+
+        const label idSolid = compModel.idSolid();
+        const wordList& solidNames = compModel.componentNames(idSolid);
+        forAll(solidNames, j)
+        {
+            const word fieldName = "Y" + solidNames[j] + stateLabels[idSolid];
+            IOField<scalar>& YSolid
+            (
+                cloud::createIOField<scalar>(fieldName, np, obr)
+            );
+
+            label i = 0;
+            forAllConstIter
+            (
+                typename Cloud<ReactingMultiphaseParcel<ParcelType>>,
+                c,
+                iter
+            )
+            {
+                const ReactingMultiphaseParcel<ParcelType>& p0 = iter();
+                YSolid[i++] = p0.YSolid()[j]*p0.Y()[SLD];
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

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

@@ -391,7 +391,7 @@ public:
 
         // I-O
 
-            //- Read
+            //- Read - composition supplied
             template<class CloudType, class CompositionType>
             static void readFields
             (
@@ -403,7 +403,7 @@ public:
             template<class CloudType>
             static void readFields(CloudType& c);
 
-            //- Write
+            //- Write - composition supplied
             template<class CloudType, class CompositionType>
             static void writeFields
             (
@@ -411,11 +411,29 @@ public:
                 const CompositionType& compModel
             );
 
-            //- Write - composition supplied
+            //- Write - no composition
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
 
+            //- Write particle fields as objects into the obr registry
+            //  - no composition
+            template<class CloudType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                objectRegistry& obr
+            );
+
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType, class CompositionType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                const CompositionType& compModel,
+                objectRegistry& obr
+            );
+
     // Ostream Operator
 
         friend Ostream& operator<< <ParcelType>

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

@@ -234,6 +234,74 @@ void Foam::ReactingParcel<ParcelType>::writeFields
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::ReactingParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ReactingParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    const CompositionType& compModel,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    label np = c.size();
+
+    if (np > 0)
+    {
+        IOField<scalar>& mass0(cloud::createIOField<scalar>("mass0", np, obr));
+
+        label i = 0;
+        forAllConstIter(typename Cloud<ReactingParcel<ParcelType>>, c, iter)
+        {
+            const ReactingParcel<ParcelType>& p = iter();
+            mass0[i++] = p.mass0_;
+        }
+
+        // Write the composition fractions
+        const wordList& phaseTypes = compModel.phaseTypes();
+        wordList stateLabels(phaseTypes.size(), "");
+        if (compModel.nPhase() == 1)
+        {
+            stateLabels = compModel.stateLabels()[0];
+        }
+
+        forAll(phaseTypes, j)
+        {
+            const word fieldName = "Y" + phaseTypes[j] + stateLabels[j];
+            IOField<scalar>& Y
+            (
+                cloud::createIOField<scalar>(fieldName, np, obr)
+            );
+
+            label i = 0;
+            forAllConstIter
+            (
+                typename Cloud<ReactingParcel<ParcelType>>,
+                c,
+                iter
+            )
+            {
+                const ReactingParcel<ParcelType>& p = iter();
+                Y[i++] = p.Y()[j];
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

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

@@ -450,6 +450,10 @@ public:
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType>
+            static void writeObjects(const CloudType& c, objectRegistry& obr);
+
 
     // Ostream Operator
 

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

@@ -137,6 +137,33 @@ void Foam::ThermoParcel<ParcelType>::writeFields(const CloudType& c)
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::ThermoParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    label np = c.size();
+
+    IOField<scalar>& T(cloud::createIOField<scalar>("T", np, obr));
+    IOField<scalar>& Cp(cloud::createIOField<scalar>("Cp", np, obr));
+
+    label i = 0;
+    forAllConstIter(typename Cloud<ThermoParcel<ParcelType>>, c, iter)
+    {
+        const ThermoParcel<ParcelType>& p = iter();
+
+        T[i] = p.T_;
+        Cp[i] = p.Cp_;
+        i++;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

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

@@ -32,6 +32,8 @@ License
 #include "ConeInjection.H"
 #include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
+#include "InjectedParticleDistributionInjection.H"
+#include "InjectedParticleInjection.H"
 #include "InflationInjection.H"
 #include "KinematicLookupTableInjection.H"
 #include "ManualInjection.H"
@@ -49,6 +51,8 @@ License
     makeInjectionModelType(ConeInjection, CloudType);                          \
     makeInjectionModelType(ConeNozzleInjection, CloudType);                    \
     makeInjectionModelType(FieldActivatedInjection, CloudType);                \
+    makeInjectionModelType(InjectedParticleDistributionInjection, CloudType);  \
+    makeInjectionModelType(InjectedParticleInjection, CloudType);              \
     makeInjectionModelType(InflationInjection, CloudType);                     \
     makeInjectionModelType(KinematicLookupTableInjection, CloudType);          \
     makeInjectionModelType(ManualInjection, CloudType);                        \

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

@@ -0,0 +1,525 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "InjectedParticleDistributionInjection.H"
+#include "mathematicalConstants.H"
+#include "PackedBoolList.H"
+#include "injectedParticleCloud.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::InjectedParticleDistributionInjection<CloudType>::initialise()
+{
+    injectedParticleCloud ipCloud(this->owner().mesh(), cloudName_);
+
+    List<label> tag(ipCloud.size());
+    List<point> position(ipCloud.size());
+    List<vector> U(ipCloud.size());
+    List<scalar> soi(ipCloud.size());
+    List<scalar> d(ipCloud.size());
+
+    // Flatten all data
+    label particlei = 0;
+    forAllConstIter(injectedParticleCloud, ipCloud, iter)
+    {
+        const injectedParticle& p = iter();
+        tag[particlei] = p.tag();
+        position[particlei] = p.position();
+        U[particlei] = p.U();
+        soi[particlei] = p.soi();
+        d[particlei] = p.d();
+        particlei++;
+    }
+
+    // Combine all proc data
+    if (Pstream::parRun())
+    {
+        List<List<label>> procTag(Pstream::nProcs());
+        procTag[Pstream::myProcNo()].transfer(tag);
+        Pstream::gatherList(procTag);
+        Pstream::scatterList(procTag);
+        tag =
+            ListListOps::combine<List<label>>
+            (
+                procTag, accessOp<List<label>>()
+            );
+
+        List<List<point>> procPosition(Pstream::nProcs());
+        procPosition[Pstream::myProcNo()].transfer(position);
+        Pstream::gatherList(procPosition);
+        Pstream::scatterList(procPosition);
+        position =
+            ListListOps::combine<List<point>>
+            (
+                procPosition, accessOp<List<point>>()
+            );
+
+        List<List<vector>> procU(Pstream::nProcs());
+        procU[Pstream::myProcNo()].transfer(U);
+        Pstream::gatherList(procU);
+        Pstream::scatterList(procU);
+        U =
+            ListListOps::combine<List<vector>>
+            (
+                procU, accessOp<List<vector>>()
+            );
+
+        List<List<scalar>> procSOI(Pstream::nProcs());
+        procSOI[Pstream::myProcNo()].transfer(soi);
+        Pstream::gatherList(procSOI);
+        Pstream::scatterList(procSOI);
+        soi =
+            ListListOps::combine<List<scalar>>
+            (
+                procSOI, accessOp<List<scalar>>()
+            );
+
+        List<List<scalar>> procD(Pstream::nProcs());
+        procD[Pstream::myProcNo()].transfer(d);
+        Pstream::gatherList(procD);
+        Pstream::scatterList(procD);
+        d =
+            ListListOps::combine<List<scalar>>
+            (
+                procD, accessOp<List<scalar>>()
+            );
+    }
+
+    label maxTag = -1;
+    forAll(tag, particlei)
+    {
+        maxTag = max(maxTag, tag[particlei]);
+    }
+
+    label nInjectors = maxTag + 1;
+    List<scalar> injStartTime(nInjectors, GREAT);
+    List<scalar> injEndTime(nInjectors, -GREAT);
+    List<DynamicList<point>> injPosition(nInjectors);
+    List<DynamicList<vector>> injU(nInjectors);
+    List<DynamicList<scalar>> injDiameter(nInjectors);
+
+    // Cache the particle information per tag
+    forAll(tag, i)
+    {
+        const label tagi = tag[i];
+        const scalar t = soi[i];
+        injStartTime[tagi] = min(t, injStartTime[tagi]);
+        injEndTime[tagi] = max(t, injEndTime[tagi]);
+        injPosition[tagi].append(position[i]);
+        injU[tagi].append(U[i]);
+        injDiameter[tagi].append(d[i]);
+    }
+
+    // Remove single particles and injectors where injection interval is 0
+    // - cannot generate a volume flow rate
+    scalar sumVolume = 0;
+    startTime_.setSize(nInjectors, 0);
+    endTime_.setSize(nInjectors, 0);
+    sizeDistribution_.setSize(nInjectors);
+    position_.setSize(nInjectors);
+    U_.setSize(nInjectors);
+    volumeFlowRate_.setSize(nInjectors, 0);
+
+    scalar minTime = GREAT;
+
+    // Populate injector properties, filtering out invalid entries
+    cachedRandom& rnd = this->owner().rndGen();
+    label injectori = 0;
+    forAll(injDiameter, i)
+    {
+        const DynamicList<scalar>& diameters = injDiameter[i];
+        const label nParticle = diameters.size();
+        const scalar dTime = injEndTime[i] - injStartTime[i];
+
+        if ((nParticle > 1) && (dTime > ROOTVSMALL))
+        {
+            minTime = min(minTime, injStartTime[i]);
+
+            startTime_[injectori] = injStartTime[i];
+            endTime_[injectori] = injEndTime[i];
+
+            // Re-sample the cloud data
+            position_[injectori].setSize(resampleSize_);
+            U_[injectori].setSize(resampleSize_);
+            List<point>& positioni = position_[injectori];
+            List<vector>& Ui = U_[injectori];
+
+            for (label samplei = 0; samplei < resampleSize_; ++samplei)
+            {
+                label posi = rnd.globalPosition<label>(0, nParticle - 1);
+                positioni[samplei] = injPosition[i][posi] + positionOffset_;
+                Ui[samplei] = injU[i][posi];
+            }
+
+            // Calculate the volume flow rate
+            scalar sumPow3 = 0;
+            forAll(diameters, particlei)
+            {
+                sumPow3 += pow3(diameters[particlei]);
+            }
+
+            const scalar volume = sumPow3*mathematical::pi/16.0;
+            sumVolume += volume;
+            volumeFlowRate_[injectori] = volume/dTime;
+
+            // Create the size distribution using the user-specified bin width
+            sizeDistribution_.set
+            (
+                injectori,
+                new distributionModels::general
+                (
+                    diameters,
+                    binWidth_,
+                    this->owner().rndGen()
+                )
+            );
+
+            injectori++;
+        }
+    }
+
+    // Resize
+    startTime_.setSize(injectori);
+    endTime_.setSize(injectori);
+    position_.setSize(injectori);
+    U_.setSize(injectori);
+    volumeFlowRate_.setSize(injectori);
+    sizeDistribution_.setSize(injectori);
+
+    // Reset start time to zero
+    forAll(startTime_, injectori)
+    {
+        startTime_[injectori] -= minTime;
+        endTime_[injectori] -= minTime;
+    }
+
+
+    // Set the volume of parcels to inject
+    this->volumeTotal_ = sumVolume;
+
+    // Provide some feedback
+    Info<< "    Read " << position_.size() << " injectors with "
+        << tag.size() << " total particles" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectedParticleDistributionInjection<CloudType>::
+InjectedParticleDistributionInjection
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    InjectionModel<CloudType>(dict, owner, modelName, typeName),
+    cloudName_(this->coeffDict().lookup("cloud")),
+    startTime_(this->template getModelProperty<scalarList>("startTime")),
+    endTime_(this->template getModelProperty<scalarList>("endTime")),
+    position_(this->template getModelProperty<List<vectorList>>("position")),
+    positionOffset_(this->coeffDict().lookup("positionOffset")),
+    volumeFlowRate_
+    (
+        this->template getModelProperty<scalarList>("volumeFlowRate")
+    ),
+    U_(this->template getModelProperty<List<vectorList>>("U")),
+    binWidth_(readScalar(this->coeffDict().lookup("binWidth"))),
+    sizeDistribution_(),
+    parcelsPerInjector_
+    (
+        ceil(readScalar(this->coeffDict().lookup("parcelsPerInjector")))
+    ),
+    resampleSize_
+    (
+        this->coeffDict().template lookupOrDefault<label>("resampleSize", 100)
+    ),
+    applyDistributionMassTotal_
+    (
+        readBool(dict.lookup("applyDistributionMassTotal"))
+    ),
+    ignoreOutOfBounds_
+    (
+        this->coeffDict().lookupOrDefault("ignoreOutOfBounds", false)
+    ),
+    nParcelsInjected_(this->parcelsAddedTotal()),
+    nParcelsInjected0_(0),
+    currentInjectori_(0),
+    currentSamplei_(0)
+{
+    if (startTime_.size())
+    {
+        // Restart
+        sizeDistribution_.setSize(startTime_.size());
+        forAll(sizeDistribution_, i)
+        {
+            const word dictName("distribution" + Foam::name(i));
+            dictionary dict;
+            this->getModelDict(dictName, dict);
+
+            sizeDistribution_.set
+            (
+                i,
+                new distributionModels::general(dict, this->owner().rndGen())
+            );
+        }
+    }
+    else
+    {
+        // Clean start
+        initialise();
+    }
+
+    // Set the mass of parcels to inject from distribution if requested
+    if (applyDistributionMassTotal_)
+    {
+        this->massTotal_ = this->volumeTotal_*this->owner().constProps().rho0();
+        Info<< "    Set mass to inject from distribution: "
+            << this->massTotal_ << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::InjectedParticleDistributionInjection<CloudType>::
+InjectedParticleDistributionInjection
+(
+    const InjectedParticleDistributionInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    cloudName_(im.cloudName_),
+    startTime_(im.startTime_),
+    endTime_(im.endTime_),
+    position_(im.position_),
+    positionOffset_(im.positionOffset_),
+    volumeFlowRate_(im.volumeFlowRate_),
+    U_(im.U_),
+    binWidth_(im.binWidth_),
+    sizeDistribution_(im.sizeDistribution_.size()),
+    parcelsPerInjector_(im.parcelsPerInjector_),
+    resampleSize_(im.resampleSize_),
+    applyDistributionMassTotal_(im.applyDistributionMassTotal_),
+    ignoreOutOfBounds_(im.ignoreOutOfBounds_),
+    nParcelsInjected_(im.nParcelsInjected_),
+    nParcelsInjected0_(im.nParcelsInjected0_),
+    currentInjectori_(0),
+    currentSamplei_(0)
+{
+    forAll(sizeDistribution_, injectori)
+    {
+        if (sizeDistribution_.set(injectori))
+        {
+            sizeDistribution_.set
+            (
+                injectori,
+                new distributionModels::general
+                (
+                    im.sizeDistribution_[injectori]
+                )
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectedParticleDistributionInjection<CloudType>::
+~InjectedParticleDistributionInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::InjectedParticleDistributionInjection<CloudType>::updateMesh()
+{}
+
+
+template<class CloudType>
+Foam::scalar
+Foam::InjectedParticleDistributionInjection<CloudType>::timeEnd() const
+{
+    return max(endTime_);
+}
+
+
+template<class CloudType>
+Foam::label
+Foam::InjectedParticleDistributionInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    // Ensure all procs know the latest parcel count
+    nParcelsInjected_ += returnReduce(nParcelsInjected0_, sumOp<label>());
+    nParcelsInjected0_ = 0;
+
+    scalar targetVolume = 0;
+    forAll(startTime_, injectori)
+    {
+        if (time1 > startTime_[injectori])
+        {
+            scalar totalDuration =
+                min(time1, endTime_[injectori]) - startTime_[injectori];
+
+            targetVolume += volumeFlowRate_[injectori]*totalDuration;
+        }
+    }
+
+    const label targetParcels =
+        round
+        (
+            scalar(startTime_.size()*parcelsPerInjector_)
+           *targetVolume/this->volumeTotal_
+        );
+
+    const label nParcels = targetParcels - nParcelsInjected_;
+
+    return nParcels;
+}
+
+
+template<class CloudType>
+Foam::scalar
+Foam::InjectedParticleDistributionInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0;
+    forAll(startTime_, injectori)
+    {
+        if ((time1 > startTime_[injectori]) && (time1 <= endTime_[injectori]))
+        {
+            scalar duration = min(time1, endTime_[injectori]) - time0;
+            volume += volumeFlowRate_[injectori]*duration;
+        }
+    }
+
+    return volume;
+}
+
+
+template<class CloudType>
+void Foam::InjectedParticleDistributionInjection<CloudType>::setPositionAndCell
+(
+    const label parcelI,
+    const label nParcels,
+    const scalar time,
+    vector& position,
+    label& cellOwner,
+    label& tetFaceI,
+    label& tetPtI
+)
+{
+    cachedRandom& rnd = this->owner().rndGen();
+    currentInjectori_ = rnd.globalPosition<label>(0, position_.size() - 1);
+    currentSamplei_ = rnd.globalPosition<label>(0, resampleSize_ - 1);
+
+    position = position_[currentInjectori_][currentSamplei_];
+
+    // Cache all mesh props for each position?
+    this->findCellAtPosition
+    (
+        cellOwner,
+        tetFaceI,
+        tetPtI,
+        position
+    );
+}
+
+
+template<class CloudType>
+void Foam::InjectedParticleDistributionInjection<CloudType>::setProperties
+(
+    const label parcelI,
+    const label,
+    const scalar,
+    typename CloudType::parcelType& parcel
+)
+{
+    // Set particle velocity
+    parcel.U() = U_[currentInjectori_][currentSamplei_];
+
+    // Set particle diameter
+    parcel.d() = sizeDistribution_[currentInjectori_].sample();
+
+    // Increment number of particles injected
+    // Note: local processor only!
+    nParcelsInjected0_++;
+}
+
+
+template<class CloudType>
+bool
+Foam::InjectedParticleDistributionInjection<CloudType>::fullyDescribed() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::InjectedParticleDistributionInjection<CloudType>::validInjection
+(
+    const label
+)
+{
+    return true;
+}
+
+
+template<class CloudType>
+void Foam::InjectedParticleDistributionInjection<CloudType>::info(Ostream& os)
+{
+    InjectionModel<CloudType>::info(os);
+
+    if (this->writeTime())
+    {
+        this->setModelProperty("startTime", startTime_);
+        this->setModelProperty("endTime", endTime_);
+        this->setModelProperty("position", position_);
+        this->setModelProperty("volumeFlowRate", volumeFlowRate_);
+        this->setModelProperty("U", U_);
+        forAll(sizeDistribution_, i)
+        {
+            const distributionModels::general& dist = sizeDistribution_[i];
+            const word dictName("distribution" + Foam::name(i));
+            dictionary dict(dist.writeDict(dictName));
+            this->setModelProperty(dictName, dict);
+        }
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::InjectedParticleDistributionInjection
+
+Description
+    Interrogates an injectedParticleCloud to convert the raw particle data
+    into a set of 'binned' injectors.
+
+    The bins are set according to the particle \c tag property, from which:
+    - diameters are converted into \c general distributions with a
+      user-specified bin width
+    - raw velocity and diameter data are resampled and stored to provide
+      variations per injector
+
+    The mass to inject can be set according to the raw input data mass total
+    by using the \c applyDistributionMassTotal switch
+
+
+Usage
+    \verbatim
+    model1
+    {
+        type            injectedParticleDistributionInjection;
+        SOI             0;
+        parcelBasisType mass;
+        cloud           eulerianParticleCloud;
+        positionOffset  (-0.025 2 -0.025);
+        binWidth        0.1e-3;
+        parcelsPerInjector 500;
+        resampleSize    100;                    // optional
+        applyDistributionMassTotal yes;
+
+        // Placeholder only when using applyDistributionMassTotal
+        massTotal       0;
+    }
+    \endverbatim
+
+Note
+    The each injector location is assumed to be operating under steady
+    conditions, i.e. using a constant flow rate profile
+
+SourceFiles
+    InjectedParticleDistributionInjection.C
+
+See also
+    Foam::injectedParticle
+    Foam::injectedParticleCloud
+    Foam::functionObjects::extractEulerianParticles
+    Foam::distributionModels::general
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef InjectedParticleDistributionInjection_H
+#define InjectedParticleDistributionInjection_H
+
+#include "InjectionModel.H"
+#include "general.H"
+#include "Switch.H"
+#include "vectorList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+           Class InjectedParticleDistributionInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class InjectedParticleDistributionInjection
+:
+    public InjectionModel<CloudType>
+{
+protected:
+
+    // Protected data
+
+        //- Name of cloud used to seed the new particles
+        const word cloudName_;
+
+        //- List of start time per injector
+        scalarList startTime_;
+
+        //- List of end time per injector
+        scalarList endTime_;
+
+        //- List of position per injector
+        List<vectorList> position_;
+
+        //- Position offset to apply to input positions
+        vector positionOffset_;
+
+        //- List of volume flow rate per injector [m3/s]
+        scalarList volumeFlowRate_;
+
+        //- List of parcel velocity per injector
+        List<vectorList> U_;
+
+        //- Bin width when generating particle distributions
+        scalar binWidth_;
+
+        //- List of size distribution model per injector
+        PtrList<distributionModels::general> sizeDistribution_;
+
+        //- Target number of parcels to inject per injector
+        scalar parcelsPerInjector_;
+
+        //- Resample size
+        label resampleSize_;
+
+        //- Flag to apply mass calculated from distribution instead of
+        //  InjectionModel massTotal_
+        bool applyDistributionMassTotal_;
+
+        //- Flag to suppress errors if particle injection site is out-of-bounds
+        Switch ignoreOutOfBounds_;
+
+        //- Running total of number of parcels injected
+        label nParcelsInjected_;
+
+        //- Number of parcels injected in last step (local proc only)
+        label nParcelsInjected0_;
+
+        //- Current injector
+        label currentInjectori_;
+
+        //- Current sample
+        label currentSamplei_;
+
+
+    // Protected Member Functions
+
+        //- Initialise injectors
+        void initialise();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("injectedParticleDistributionInjection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        InjectedParticleDistributionInjection
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& modelName
+        );
+
+        //- Construct copy
+        InjectedParticleDistributionInjection
+        (
+            const InjectedParticleDistributionInjection<CloudType>& im
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<InjectionModel<CloudType> > clone() const
+        {
+            return autoPtr<InjectionModel<CloudType> >
+            (
+                new InjectedParticleDistributionInjection<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~InjectedParticleDistributionInjection();
+
+
+    // Member Functions
+
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
+        //- Return the end-of-injection time
+        scalar timeEnd() const;
+
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
+        // Injection geometry
+
+            //- Set the injection position and owner cell, tetFace and tetPt
+            virtual void setPositionAndCell
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                vector& position,
+                label& cellOwner,
+                label& tetFaceI,
+                label& tetPtI
+            );
+
+            //- Set the parcel properties
+            virtual void setProperties
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                typename CloudType::parcelType& parcel
+            );
+
+            //- Flag to identify whether model fully describes the parcel
+            virtual bool fullyDescribed() const;
+
+            //- Return flag to identify whether or not injection of parcelI is
+            //  permitted
+            virtual bool validInjection(const label parcelI);
+
+
+        // I-O
+
+            //- Write injection info to stream
+            void info(Ostream& os);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "InjectedParticleDistributionInjection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,395 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "InjectedParticleInjection.H"
+#include "mathematicalConstants.H"
+#include "PackedBoolList.H"
+#include "SortableList.H"
+#include "injectedParticleCloud.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::InjectedParticleInjection<CloudType>::initialise()
+{
+    const injectedParticleCloud cloud(this->owner().mesh(), cloudName_);
+
+    label nParticles = cloud.size();
+    List<scalar> time(nParticles);
+    List<vector> position(nParticles);
+    List<scalar> diameter(nParticles);
+    List<vector> U(nParticles);
+
+    label particlei = 0;
+
+    forAllConstIter(injectedParticleCloud, cloud, iter)
+    {
+        const injectedParticle& p = iter();
+
+        time[particlei] = p.soi();
+        position[particlei] = p.position() + positionOffset_;
+        diameter[particlei] = p.d();
+        U[particlei] = p.U();
+
+        particlei++;
+    }
+
+    // Combine all proc data
+    if (Pstream::parRun())
+    {
+        List<List<scalar>> procTime(Pstream::nProcs());
+        procTime[Pstream::myProcNo()].transfer(time);
+        Pstream::gatherList(procTime);
+        Pstream::scatterList(procTime);
+        time =
+            ListListOps::combine<List<scalar>>
+            (
+                procTime, accessOp<List<scalar>>()
+            );
+
+        List<List<point>> procPosition(Pstream::nProcs());
+        procPosition[Pstream::myProcNo()].transfer(position);
+        Pstream::gatherList(procPosition);
+        Pstream::scatterList(procPosition);
+        position =
+            ListListOps::combine<List<point>>
+            (
+                procPosition, accessOp<List<point>>()
+            );
+
+        List<List<scalar>> procD(Pstream::nProcs());
+        procD[Pstream::myProcNo()].transfer(diameter);
+        Pstream::gatherList(procD);
+        Pstream::scatterList(procD);
+        diameter =
+            ListListOps::combine<List<scalar>>
+            (
+                procD, accessOp<List<scalar>>()
+            );
+
+        List<List<vector>> procU(Pstream::nProcs());
+        procU[Pstream::myProcNo()].transfer(U);
+        Pstream::gatherList(procU);
+        Pstream::scatterList(procU);
+        U =
+            ListListOps::combine<List<vector>>
+            (
+                procU, accessOp<List<vector>>()
+            );
+    }
+
+    nParticles = time.size();
+
+    // Reset SOI according to user selection
+    scalar minTime = min(time);
+    forAll(time, i)
+    {
+        time[i] -= minTime;
+    }
+
+    // Sort and renumber to ensure lists in ascending time
+    labelList sortedIndices;
+    Foam::sortedOrder(time, sortedIndices);
+    time_ = UIndirectList<scalar>(time, sortedIndices);
+    position_ = UIndirectList<point>(position, sortedIndices);
+    diameter_ = UIndirectList<scalar>(diameter, sortedIndices);
+    U_ = UIndirectList<vector>(U, sortedIndices);
+
+    // Pre-calculate injected particle volumes
+    List<scalar> volume(nParticles);
+    scalar sumVolume = 0;
+    forAll(volume, i)
+    {
+        scalar vol = pow3(diameter_[i])*mathematical::pi/16.0;
+        volume[i] = vol;
+        sumVolume += vol;
+    }
+    volume_.transfer(volume);
+
+    // Set the volume of particles to inject
+    this->volumeTotal_ = sumVolume;
+
+    // Provide some feedback
+    Info<< "    Read " << nParticles << " particles" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectedParticleInjection<CloudType>::InjectedParticleInjection
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    InjectionModel<CloudType>(dict, owner, modelName, typeName),
+    cloudName_(this->coeffDict().lookup("cloud")),
+    injectorCells_(),
+    injectorTetFaces_(),
+    injectorTetPts_(),
+    time_(this->template getModelProperty<scalarList>("time")),
+    position_(this->template getModelProperty<vectorList>("position")),
+    positionOffset_(this->coeffDict().lookup("positionOffset")),
+    diameter_(this->template getModelProperty<scalarList>("diameter")),
+    U_(this->template getModelProperty<vectorList>("U")),
+    volume_(this->template getModelProperty<scalarList>("volume")),
+    ignoreOutOfBounds_
+    (
+        this->coeffDict().lookupOrDefault("ignoreOutOfBounds", false)
+    ),
+    currentParticlei_
+    (
+        this->template getModelProperty<label>
+        (
+            "currentParticlei",
+            -1
+        )
+    )
+{
+    if (this->parcelBasis_ != InjectionModel<CloudType>::pbFixed)
+    {
+        FatalErrorInFunction
+            << "Injector model: " << this->modelName()
+            << " Parcel basis must be set to fixed"
+            << exit(FatalError);
+    }
+
+    if (!time_.size())
+    {
+        // Clean start
+        initialise();
+    }
+
+    injectorCells_.setSize(position_.size());
+    injectorTetFaces_.setSize(position_.size());
+    injectorTetPts_.setSize(position_.size());
+
+    updateMesh();
+
+    this->massTotal_ = this->volumeTotal_*this->owner().constProps().rho0();
+}
+
+
+template<class CloudType>
+Foam::InjectedParticleInjection<CloudType>::InjectedParticleInjection
+(
+    const InjectedParticleInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    cloudName_(im.cloudName_),
+    injectorCells_(im.injectorCells_),
+    injectorTetFaces_(im.injectorTetFaces_),
+    injectorTetPts_(im.injectorTetPts_),
+    time_(im.time_),
+    position_(im.position_),
+    positionOffset_(im.positionOffset_),
+    diameter_(im.diameter_),
+    U_(im.U_),
+    volume_(im.volume_),
+    ignoreOutOfBounds_(im.ignoreOutOfBounds_),
+    currentParticlei_(im.currentParticlei_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectedParticleInjection<CloudType>::~InjectedParticleInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::InjectedParticleInjection<CloudType>::updateMesh()
+{
+    label nRejected = 0;
+
+    PackedBoolList keep(position_.size(), true);
+
+    forAll(position_, particlei)
+    {
+        if
+        (
+            !this->findCellAtPosition
+            (
+                injectorCells_[particlei],
+                injectorTetFaces_[particlei],
+                injectorTetPts_[particlei],
+                position_[particlei],
+                !ignoreOutOfBounds_
+            )
+        )
+        {
+            keep[particlei] = false;
+            nRejected++;
+        }
+    }
+
+    if (nRejected > 0)
+    {
+        inplaceSubset(keep, time_);
+        inplaceSubset(keep, position_);
+        inplaceSubset(keep, diameter_);
+        inplaceSubset(keep, U_);
+        inplaceSubset(keep, volume_);
+        inplaceSubset(keep, injectorCells_);
+        inplaceSubset(keep, injectorTetFaces_);
+        inplaceSubset(keep, injectorTetPts_);
+
+        Info<< "    " << nRejected
+            << " particles ignored, out of bounds" << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectedParticleInjection<CloudType>::timeEnd() const
+{
+    return max(time_);
+}
+
+
+template<class CloudType>
+Foam::label Foam::InjectedParticleInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    label nParticles = 0;
+    forAll(time_, particlei)
+    {
+        if ((time_[particlei] >= time0) && (time_[particlei] < time1))
+        {
+            nParticles++;
+        }
+    }
+
+    return nParticles;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectedParticleInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar sumVolume = 0;
+    forAll(time_, particlei)
+    {
+        if ((time_[particlei] >= time0) && (time_[particlei] < time1))
+        {
+            sumVolume += volume_[particlei];
+        }
+    }
+
+    return sumVolume;
+}
+
+
+template<class CloudType>
+void Foam::InjectedParticleInjection<CloudType>::setPositionAndCell
+(
+    const label parceli,
+    const label nParcels,
+    const scalar time,
+    vector& position,
+    label& cellOwner,
+    label& tetFacei,
+    label& tetPti
+)
+{
+    // Note: optimisation - consume particle from lists to reduce storage
+    // as injection proceeds
+
+    currentParticlei_++;
+
+    position = position_[currentParticlei_];
+    cellOwner = injectorCells_[currentParticlei_];
+    tetFacei = injectorTetFaces_[currentParticlei_];
+    tetPti = injectorTetPts_[currentParticlei_];
+}
+
+
+template<class CloudType>
+void Foam::InjectedParticleInjection<CloudType>::setProperties
+(
+    const label parceli,
+    const label,
+    const scalar,
+    typename CloudType::parcelType& parcel
+)
+{
+    // Set particle velocity
+    parcel.U() = U_[currentParticlei_];
+
+    // Set particle diameter
+    parcel.d() = diameter_[currentParticlei_];
+}
+
+
+template<class CloudType>
+bool Foam::InjectedParticleInjection<CloudType>::fullyDescribed() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::InjectedParticleInjection<CloudType>::validInjection
+(
+    const label
+)
+{
+    return true;
+}
+
+
+template<class CloudType>
+void Foam::InjectedParticleInjection<CloudType>::info(Ostream& os)
+{
+    InjectionModel<CloudType>::info(os);
+
+    if (this->writeTime())
+    {
+        this->setModelProperty("currentParticlei", currentParticlei_);
+        this->setModelProperty("time", time_);
+        this->setModelProperty("position", position_);
+        this->setModelProperty("diameter", diameter_);
+        this->setModelProperty("U", U_);
+        this->setModelProperty("volume", volume_);
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::InjectedParticleInjection
+
+Description
+    Replays an set of particle data based on an injectedParticleCloud, using
+    the assumption of one particle per parcel.
+
+Usage
+    \verbatim
+    model1
+    {
+        type            injectedParticleInjection;
+        SOI             0;
+        massTotal       0; // Place holder only
+        parcelBasisType fixed;
+        nParticle       1; // 1 particle per parcel
+        cloud           eulerianParticleCloud;
+        positionOffset  (-0.025 2 -0.025);
+    }
+    \endverbatim
+
+SourceFiles
+    InjectedParticleInjection.C
+
+See also
+    Foam::injectedParticle
+    Foam::injectedParticleCloud
+    Foam::functionObjects::extractEulerianParticles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef InjectedParticleInjection_H
+#define InjectedParticleInjection_H
+
+#include "InjectionModel.H"
+#include "Switch.H"
+#include "vectorList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class injectedParticleCloud;
+
+/*---------------------------------------------------------------------------*\
+                 Class InjectedParticleInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class InjectedParticleInjection
+:
+    public InjectionModel<CloudType>
+{
+protected:
+
+    // Protected data
+
+        //- Name of cloud used to seed the new particles
+        const word cloudName_;
+
+        //- List of cell label per injector
+        labelList injectorCells_;
+
+        //- List of tetFace label per injector
+        labelList injectorTetFaces_;
+
+        //- List of tetPt label per injector
+        labelList injectorTetPts_;
+
+        //- List of injection time per particle [s]
+        scalarList time_;
+
+        //- List of position per particle [m]
+        vectorList position_;
+
+        //- Position offset to apply to input positions
+        vector positionOffset_;
+
+        //- List of diameter per particle [m]
+        scalarList diameter_;
+
+        //- List of velocity per particle [m/s]
+        vectorList U_;
+
+        //- List of volume per particle [m3]
+        scalarList volume_;
+
+        //- Flag to suppress errors if particle injection site is out-of-bounds
+        Switch ignoreOutOfBounds_;
+
+        //- Index of current particle
+        label currentParticlei_;
+
+
+        // Protected Member Functions
+
+            //- Initialise injectors
+            void initialise();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("injectedParticleInjection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        InjectedParticleInjection
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& modelName
+        );
+
+        //- Construct copy
+        InjectedParticleInjection
+        (
+            const InjectedParticleInjection<CloudType>& im
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<InjectionModel<CloudType> > clone() const
+        {
+            return autoPtr<InjectionModel<CloudType> >
+            (
+                new InjectedParticleInjection<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~InjectedParticleInjection();
+
+
+    // Member Functions
+
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
+        //- Return the end-of-injection time
+        scalar timeEnd() const;
+
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
+        // Injection geometry
+
+            //- Set the injection position and owner cell, tetFace and tetPt
+            virtual void setPositionAndCell
+            (
+                const label parceli,
+                const label nParcels,
+                const scalar time,
+                vector& position,
+                label& cellOwner,
+                label& tetFacei,
+                label& tetPti
+            );
+
+            //- Set the parcel properties
+            virtual void setProperties
+            (
+                const label parceli,
+                const label nParcels,
+                const scalar time,
+                typename CloudType::parcelType& parcel
+            );
+
+            //- Flag to identify whether model fully describes the parcel
+            virtual bool fullyDescribed() const;
+
+            //- Return flag to identify whether or not injection of parcelI is
+            //  permitted
+            virtual bool validInjection(const label parceli);
+
+
+        // I-O
+
+            //- Write injection info to stream
+            void info(Ostream& os);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "InjectedParticleInjection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -262,8 +262,8 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
 :
     CloudSubModelBase<CloudType>(owner),
     SOI_(0.0),
-    volumeTotal_(0.0),
-    massTotal_(0.0),
+    volumeTotal_(this->template getModelProperty<scalar>("volumeTotal")),
+    massTotal_(0),
     massFlowRate_(owner.db().time(), "massFlowRate"),
     massInjected_(this->template getModelProperty<scalar>("massInjected")),
     nInjections_(this->template getModelProperty<label>("nInjections")),
@@ -291,8 +291,8 @@ Foam::InjectionModel<CloudType>::InjectionModel
 :
     CloudSubModelBase<CloudType>(modelName, owner, dict, typeName, modelType),
     SOI_(0.0),
-    volumeTotal_(0.0),
-    massTotal_(0.0),
+    volumeTotal_(this->template getModelProperty<scalar>("volumeTotal")),
+    massTotal_(0),
     massFlowRate_(owner.db().time(), "massFlowRate"),
     massInjected_(this->template getModelProperty<scalar>("massInjected")),
     nInjections_(this->template getModelProperty<scalar>("nInjections")),
@@ -438,7 +438,6 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
 
     if (prepareForNextTimeStep(time, newParcels, newVolumeFraction))
     {
-
         const scalar trackTime = this->owner().solution().trackTime();
         const polyMesh& mesh = this->owner().mesh();
         typename TrackData::cloudType& cloud = td.cloud();
@@ -659,6 +658,7 @@ void Foam::InjectionModel<CloudType>::info(Ostream& os)
 
     if (this->writeTime())
     {
+        this->setModelProperty("volumeTotal", volumeTotal_);
         this->setModelProperty("massInjected", massInjected_);
         this->setModelProperty("nInjections", nInjections_);
         this->setModelProperty("parcelsAddedTotal", parcelsAddedTotal_);

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H

@@ -506,6 +506,24 @@ public:
             template<class CloudType>
             static void writeFields(const CloudType& c);
 
+            //- Write particle fields as objects into the obr registry
+            //  - no composition
+            template<class CloudType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                objectRegistry& obr
+            );
+
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType, class CompositionType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                const CompositionType& compModel,
+                objectRegistry& obr
+            );
+
 
     // Ostream Operator
 

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C

@@ -283,6 +283,76 @@ void Foam::SprayParcel<ParcelType>::writeFields
 }
 
 
+template<class ParcelType>
+template<class CloudType>
+void Foam::SprayParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::SprayParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    const CompositionType& compModel,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, compModel, obr);
+
+    label np = c.size();
+
+    IOField<scalar>& d0(cloud::createIOField<scalar>("d0", np, obr));
+    IOField<vector>& position0
+    (
+        cloud::createIOField<vector>("position0", np, obr)
+    );
+    IOField<scalar>& sigma(cloud::createIOField<scalar>("sigma", np, obr));
+    IOField<scalar>& mu(cloud::createIOField<scalar>("mu", np, obr));
+    IOField<scalar>& liquidCore
+    (
+        cloud::createIOField<scalar>("liquidCore", np, obr)
+    );
+    IOField<scalar>& KHindex(cloud::createIOField<scalar>("KHindex", np, obr));
+    IOField<scalar>& y(cloud::createIOField<scalar>("y", np, obr));
+    IOField<scalar>& yDot(cloud::createIOField<scalar>("yDot", np, obr));
+    IOField<scalar>& tc(cloud::createIOField<scalar>("tc", np, obr));
+    IOField<scalar>& ms(cloud::createIOField<scalar>("ms", np, obr));
+    IOField<scalar>& injector
+    (
+        cloud::createIOField<scalar>("injector", np, obr)
+    );
+    IOField<scalar>& tMom(cloud::createIOField<scalar>("tMom", np, obr));
+    IOField<scalar>& user(cloud::createIOField<scalar>("user", np, obr));
+
+    label i = 0;
+    forAllConstIter(typename Cloud<SprayParcel<ParcelType>>, c, iter)
+    {
+        const SprayParcel<ParcelType>& p = iter();
+        d0[i] = p.d0_;
+        position0[i] = p.position0_;
+        sigma[i] = p.sigma_;
+        mu[i] = p.mu_;
+        liquidCore[i] = p.liquidCore_;
+        KHindex[i] = p.KHindex_;
+        y[i] = p.y_;
+        yDot[i] = p.yDot_;
+        tc[i] = p.tc_;
+        ms[i] = p.ms_;
+        injector[i] = p.injector_;
+        tMom[i] = p.tMom_;
+        user[i] = p.user_;
+        i++;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class ParcelType>

src/lagrangian/spray/parcels/include/makeSprayParcelInjectionModels.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,8 @@ License
 #include "ConeInjection.H"
 #include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
+#include "InjectedParticleDistributionInjection.H"
+#include "InjectedParticleInjection.H"
 #include "InflationInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"
@@ -48,6 +50,8 @@ License
     makeInjectionModelType(ConeInjection, CloudType);                          \
     makeInjectionModelType(ConeNozzleInjection, CloudType);                    \
     makeInjectionModelType(FieldActivatedInjection, CloudType);                \
+    makeInjectionModelType(InjectedParticleDistributionInjection, CloudType);  \
+    makeInjectionModelType(InjectedParticleInjection, CloudType);              \
     makeInjectionModelType(InflationInjection, CloudType);                     \
     makeInjectionModelType(ManualInjection, CloudType);                        \
     makeInjectionModelType(NoInjection, CloudType);                            \

src/meshTools/Make/files

@@ -55,6 +55,8 @@ $(meshWave)/FaceCellWaveName.C
 regionSplit/regionSplit.C
 regionSplit/localPointRegion.C
 
+regionSplit2D/regionSplit2D.C
+
 indexedOctree/treeDataEdge.C
 indexedOctree/treeDataFace.C
 indexedOctree/treeDataPoint.C

src/meshTools/regionSplit2D/regionSplit2D.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSplit2D.H"
+#include "polyMesh.H"
+#include "PatchEdgeFaceWave.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+Foam::regionSplit2D::regionSplit2D
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& patch,
+    const boolList& blockedFaces,
+    const label offset
+)
+:
+    labelList(patch.size(), -1),
+    nRegions_(0)
+{
+    globalIndex globalFaces(blockedFaces.size());
+    label regionI = globalFaces.toGlobal(0);
+    List<patchEdgeFaceRegion> allEdgeInfo(patch.nEdges());
+    List<patchEdgeFaceRegion> allFaceInfo(patch.size());
+    DynamicList<label> changedEdges;
+    DynamicList<patchEdgeFaceRegion> changedRegions;
+    label nBlockedFaces = 0;
+    forAll(blockedFaces, faceI)
+    {
+        if (blockedFaces[faceI])
+        {
+            const labelList& fEdges = patch.faceEdges()[faceI];
+            forAll(fEdges, feI)
+            {
+                changedEdges.append(fEdges[feI]);
+
+                // Append globally unique value
+                changedRegions.append(regionI);
+            }
+            nBlockedFaces++;
+            regionI++;
+        }
+        else
+        {
+            // Block all non-seeded faces from the walk
+            allFaceInfo[faceI] = -2;
+        }
+    }
+
+    // Early exit if there are no blocked faces
+    if (returnReduce(nBlockedFaces, sumOp<label>()) == 0)
+    {
+        return;
+    }
+
+    PatchEdgeFaceWave
+    <
+        indirectPrimitivePatch,
+        patchEdgeFaceRegion
+    >
+    (
+        mesh,
+        patch,
+        changedEdges,
+        changedRegions,
+        allEdgeInfo,
+        allFaceInfo,
+        returnReduce(patch.nEdges(), sumOp<label>())
+    );
+
+
+    // Map from regions to local compact indexing
+    // - only for regions that originate from this processor
+    Map<label> regionToCompactAddr(changedRegions.size());
+    label compactRegionI = 0;
+    forAll(allFaceInfo, faceI)
+    {
+        label regionI = allFaceInfo[faceI].region();
+        if
+        (
+            globalFaces.isLocal(regionI)
+         && regionToCompactAddr.insert(regionI, compactRegionI)
+        )
+        {
+            compactRegionI++;
+        }
+    }
+
+    // In-place renumber the local regionI to global (compact) regionI
+    globalIndex giCompact(compactRegionI);
+    forAllIter(Map<label>, regionToCompactAddr, iter)
+    {
+        label compactRegionI = iter();
+        iter() = giCompact.toGlobal(compactRegionI);
+    }
+
+    // Ensure regionToCompactAddr consistent across all processors
+    // - not concerned about the op (keys are unique)
+    // - map size will be the number of regions in the set of faces
+    Pstream::mapCombineGather(regionToCompactAddr, minEqOp<label>());
+    Pstream::mapCombineScatter(regionToCompactAddr);
+
+    nRegions_ = regionToCompactAddr.size();
+
+    // Set the region index per face
+    forAll(allFaceInfo, faceI)
+    {
+        label regionI = allFaceInfo[faceI].region();
+        if (regionI >= 0)
+        {
+            this->operator[](faceI) = regionToCompactAddr[regionI] + offset;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionSplit2D::~regionSplit2D()
+{}
+
+
+// ************************************************************************* //

src/meshTools/regionSplit2D/regionSplit2D.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionSplit2D
+
+Description
+    Splits a patch into regions based on a mask field.  Result is a globally
+    consistent label list of region index per patch face.
+
+SourceFiles
+    regionSplit2D.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionSplit2D_H
+#define regionSplit2D_H
+
+#include "DynamicList.H"
+#include "boolList.H"
+#include "labelList.H"
+#include "indirectPrimitivePatch.H"
+#include "patchEdgeFaceRegion.H"
+#include "globalIndex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class faceZone;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                        Class regionSplit2D Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionSplit2D
+:
+    public labelList
+{
+    // Private data
+
+        //- Number of regions
+        label nRegions_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        regionSplit2D(const regionSplit2D&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const regionSplit2D&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh and list of blocked faces
+        regionSplit2D
+        (
+            const polyMesh& mesh,
+            const indirectPrimitivePatch& patch,
+            const boolList& blockedFaces,
+            const label offset = 0
+        );
+
+
+    //- Destructor
+    ~regionSplit2D();
+
+
+    // Member Functions
+
+        //- Return the global number of regions
+        label nRegions() const
+        {
+            return nRegions_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/Allrun

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+(cd eulerianInjection && ./Allrun)
+(cd lagrangianParticleInjection && ./Allrun)
+(cd lagrangianDistributionInjection && ./Allrun)

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/0.orig/U

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus-develop                          |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    top
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    bottom
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    sides
+    {
+        type            pressureInletOutletVelocity;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/0.orig/alpha.water

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus-develop                          |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    top
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    sides
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/0.orig/p_rgh

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus-develop                          |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    top
+    {
+        type            fixedFluxPressure;
+        gradient        uniform 0;
+        value           uniform 0;
+    }
+    bottom
+    {
+        type            fixedFluxPressure;
+        gradient        uniform 0;
+        value           uniform 0;
+    }
+    sides
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/Allclean

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+\rm -rf 0
+
+#------------------------------------------------------------------------------

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/Allrun

@@ -0,0 +1,20 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication  -s createBlockage topoSet -dict system/topoSetDict.createBlockage
+runApplication subsetMesh -overwrite blockage -patch bottom
+runApplication -s createPatch topoSet -dict system/topoSetDict.createPatch
+runApplication -s createCollector topoSet -dict system/topoSetDict.createCollector
+
+restore0Dir
+
+runApplication decomposePar
+
+runParallel setFields
+runParallel $(getApplication)
+
+#------------------------------------------------------------------------------

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/constant/transportProperties

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+    transportModel  Newtonian;
+    nu              [0 2 -1 0 0 0 0] 1e-06;
+    rho             [1 -3 0 0 0 0 0] 1000;
+}
+
+air
+{
+    transportModel  Newtonian;
+    nu              [0 2 -1 0 0 0 0] 1.48e-05;
+    rho             [1 -3 0 0 0 0 0] 1;
+}
+
+sigma           [1 0 -2 0 0 0 0] 0.07;
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/constant/turbulenceProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/blockMeshDict

@@ -0,0 +1,75 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.05;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 1)
+    (1 0 1)
+    (1 1 1)
+    (0 1 1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (75 75 75) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    top
+    {
+        type        patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    bottom
+    {
+        type        wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    sides
+    {
+        type        patch;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/controlDict

@@ -0,0 +1,70 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     interFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.2;
+
+deltaT          1e-4;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.001;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           0.5;
+maxAlphaCo      1;
+
+maxDeltaT       1;
+
+functions
+{
+    extractEulerianParticles1
+    {
+        type            extractEulerianParticles;
+        libs            ("libfieldFunctionObjects.so");
+        writeControl    writeTime;
+        faceZone        collector;
+        alpha           alpha.water;
+        nLocations      20;
+    }
+}
+
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/decomposeParDict

@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method          simple;
+
+simpleCoeffs
+{
+    n           (2 1 2);
+    delta       0.001;
+}
+
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/fvSchemes

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(rhoPhi,U)   Gauss limitedLinear 0.2;
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss linear;
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         orthogonal;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/fvSolution

@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    alpha.water
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 5;
+        cAlpha          2;
+    }
+
+    pcorr
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+        tolerance       1e-10;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+        tolerance       1e-07;
+        relTol          0.05;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance       1e-07;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor no;
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/setFieldsDict

@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.water 0
+);
+
+regions
+(
+    faceToFace
+    {
+        set         injectorFaces;
+
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+            volVectorFieldValue U (0 -15 0)
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/topoSetDict.createBlockage

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // Create blockage cellSet
+    {
+        name    blockage;
+        type    cellSet;
+        action  new;
+//        source  boxToCell;
+        source  cylinderToCell;
+        sourceInfo
+        {
+//            box     (0.015 0.010 0.015) (0.035 0.015 0.035);
+            p1      (0.025 0.02 0.025);
+            p2      (0.025 0.022 0.025);
+            radius  0.01;
+        }
+    }
+    {
+        name    blockage;
+        type    cellSet;
+        action  invert;
+    }
+);
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/topoSetDict.createCollector

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // Get all faces in cellSet
+    {
+        name    f0;
+        type    faceSet;
+        action  new;
+        source  patchToFace;
+        sourceInfo
+        {
+            name    sides;
+        }
+    }
+
+    {
+        name    collector;
+        type    faceZoneSet;
+        action  new;
+        source  setToFaceZone;
+        sourceInfo
+        {
+            faceSet f0;
+        }
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/eulerianInjection/system/topoSetDict.createPatch

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // Load initial cellSet
+    {
+        name    c0;
+        type    cellSet;
+        action  new;
+        source  cylinderToCell;
+        sourceInfo
+        {
+            p1          (0.025 1 0.025);
+            p2          (0.025 0.049 0.025);
+            radius      0.0015;
+        }
+    }
+
+    // Get all faces in cellSet and assign to injectorFaces
+    {
+        name    injectorFaces;
+        type    faceSet;
+        action  new;
+        source  cellToFace;
+        sourceInfo
+        {
+            set     c0;
+            option  all;
+        }
+    }
+
+    // Keep in injectorFaces all faces in boundary faces
+    {
+        name    injectorFaces;
+        type    faceSet;
+        action  subset;
+        source  boundaryToFace;
+        sourceInfo
+        {
+        }
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/N2

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      N2;
+}
+// ************************************************************************* //
+
+dimensions          [0 0 0 0 0 0 0];
+
+internalField       uniform 0.766;
+
+boundaryField
+{
+    walls
+    {
+        type                zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/O2

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      O2;
+}
+// ************************************************************************* //
+
+dimensions          [0 0 0 0 0 0 0];
+
+internalField       uniform 0.234;
+
+boundaryField
+{
+    walls
+    {
+        type                zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/T

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// ************************************************************************* //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 293;
+
+boundaryField
+{
+    walls
+    {
+        type                zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/U

@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// ************************************************************************* //
+
+dimensions          [0 1 -1 0 0 0 0];
+
+internalField       uniform (0 0 0);
+
+boundaryField
+{
+    walls
+    {
+        type                fixedValue;
+        value               uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/Ydefault

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Ydefault;
+}
+// ************************************************************************* //
+
+dimensions          [0 0 0 0 0 0 0];
+
+internalField       uniform 0;
+
+boundaryField
+{
+    walls
+    {
+        type                zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/alphat

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/interFoam/laminar/vofToLagrangian/lagrangianDistributionInjection/0.orig/epsilon

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 90;
+
+boundaryField
+{
+    walls
+    {
+        type            epsilonWallFunction;
+        value           uniform 90;
+    }
+}
+
+
+// ************************************************************************* //