Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -5,7 +5,7 @@
       + fvm::div(rhoPhi, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
-            fvc::div(fvc::absolute(phi, U), p)
+            divU*p
           + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
         )
        *(

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -82,12 +82,12 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
 
-        {
-            // Store divU from the previous mesh so that it can be mapped
-            // and used in correctPhi to ensure the corrected phi has the
-            // same divergence
-            volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
+        // Store divU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
 
+        {
             #include "CourantNo.H"
             #include "setDeltaT.H"
 
@@ -110,8 +110,12 @@ int main(int argc, char *argv[])
                 ghf = (g & mesh.Cf()) - ghRef;
             }
 
-            if (mesh.changing() && correctPhi)
+            if ((correctPhi && mesh.changing()) || mesh.topoChanging())
             {
+                // Calculate absolute flux from the mapped surface velocity
+                // SAF: temporary fix until mapped Uf is assessed
+                Uf = fvc::interpolate(U);
+
                 // Calculate absolute flux from the mapped surface velocity
                 phi = mesh.Sf() & Uf;
 
@@ -119,6 +123,8 @@ int main(int argc, char *argv[])
 
                 // Make the fluxes relative to the mesh motion
                 fvc::makeRelative(phi, U);
+
+                mesh.topoChanging(false);
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -87,15 +87,6 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
-
-            dgdt =
-            (
-                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
-              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
-            );
-
             phi = phiHbyA + p_rghEqnIncomp.flux();
 
             U = HbyA
@@ -116,10 +107,16 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
-    // Correct p_rgh for consistency with p and the updated densities
+    p = max(p_rgh + rho*gh, pMin);
     p_rgh = p - rho*gh;
     p_rgh.correctBoundaryConditions();
 
+    dgdt =
+    (
+        pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+      - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+    );
+
     K = 0.5*magSqr(U);
 
     Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -98,6 +98,7 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
+            volScalarField divU(fvc::div(fvc::absolute(phi, U)));
             #include "TEqn.H"
 
             // --- Pressure corrector loop

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
+
+// ************************************************************************* //

src/topoChangerFvMesh/Make/files

@@ -1,5 +1,6 @@
 topoChangerFvMesh/topoChangerFvMesh.C
 rawTopoChangerFvMesh/rawTopoChangerFvMesh.C
+dynamicMotionSolverTopoFvMesh/dynamicMotionSolverTopoFvMesh.C
 
 /*
 linearValveFvMesh/linearValveFvMesh.C

src/topoChangerFvMesh/Make/options

@@ -2,10 +2,12 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -ldynamicFvMesh \
     -ldynamicMesh \
-    -lmeshTools
+    -lmeshTools \
+    -lsurfMesh

src/topoChangerFvMesh/dynamicMotionSolverTopoFvMesh/dynamicMotionSolverTopoFvMesh.C

@@ -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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "addToRunTimeSelectionTable.H"
+#include "dynamicMotionSolverTopoFvMesh.H"
+#include "mapPolyMesh.H"
+#include "OBJstream.H"
+#include "Time.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dynamicMotionSolverTopoFvMesh, 0);
+
+    addToRunTimeSelectionTable
+    (
+        dynamicFvMesh,
+        dynamicMotionSolverTopoFvMesh,
+        IOobject
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dynamicMotionSolverTopoFvMesh::dynamicMotionSolverTopoFvMesh
+(
+    const IOobject& io
+)
+:
+    topoChangerFvMesh(io),
+    motionPtr_(motionSolver::New(*this))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dynamicMotionSolverTopoFvMesh::~dynamicMotionSolverTopoFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::dynamicMotionSolverTopoFvMesh::update()
+{
+    // Do mesh changes (not using inflation - points added directly into mesh)
+    autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(false);
+
+    if (topoChangeMap.valid())
+    {
+        Info << "Executing mesh topology update" << endl;
+        motionPtr_->updateMesh(topoChangeMap());
+
+        setV0() = V();
+
+        pointField newPoints(motionPtr_->newPoints());
+        movePoints(newPoints);
+
+        if (debug)
+        {
+            OBJstream osOld("oldPts_" + time().timeName() + ".obj");
+            const pointField& oldPts = oldPoints();
+            forAll(oldPts, i)
+            {
+                osOld.write(oldPts[i]);
+            }
+
+            OBJstream osNew("newPts_" + time().timeName() + ".obj");
+            forAll(points(), i)
+            {
+                osNew.write(points()[i]);
+            }
+        }
+    }
+    else
+    {
+        // Calculate the new point positions using the motion solver
+        pointField newPoints(motionPtr_->newPoints());
+
+        // The mesh now contains the cells with zero volume
+        Info << "Executing mesh motion" << endl;
+        movePoints(newPoints);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/topoChangerFvMesh/dynamicMotionSolverTopoFvMesh/dynamicMotionSolverTopoFvMesh.H

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::dynamicMotionSolverTopoFvMesh
+
+Description
+    Sample topoChangerFvMesh that moves an object using a motion solver
+
+SourceFiles
+    dynamicMotionSolverTopoFvMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicMotionSolverTopoFvMesh_H
+#define dynamicMotionSolverTopoFvMesh_H
+
+#include "topoChangerFvMesh.H"
+#include "motionSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                      Class dynamicMotionSolverTopoFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicMotionSolverTopoFvMesh
+:
+    public topoChangerFvMesh
+{
+    // Private data
+
+        //- Motion dictionary
+        dictionary motionDict_;
+
+        //- Motion solver
+        autoPtr<motionSolver> motionPtr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        dynamicMotionSolverTopoFvMesh(const dynamicMotionSolverTopoFvMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dynamicMotionSolverTopoFvMesh&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicMotionSolverTopoFvMesh");
+
+
+    // Constructors
+
+        //- Construct from database
+        explicit dynamicMotionSolverTopoFvMesh(const IOobject& io);
+
+
+    //- Destructor
+    virtual ~dynamicMotionSolverTopoFvMesh();
+
+
+    // Member Functions
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/T

@@ -0,0 +1,63 @@
+/*--------------------------------*- 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;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 300;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 300;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 300;
+        value           uniform 300;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    farFieldWall
+    {
+        type            zeroGradient;
+    }
+    sphere
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/T.air

@@ -0,0 +1,63 @@
+/*--------------------------------*- 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;
+    location    "0";
+    object      T.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 300;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 300;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 300;
+        value           uniform 300;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    farFieldWall
+    {
+        type            zeroGradient;
+    }
+    sphere
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/T.water

@@ -0,0 +1,63 @@
+/*--------------------------------*- 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;
+    location    "0";
+    object      T.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 300;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 300;
+    }
+    top
+    {
+        type            inletOutlet;
+        inletValue      uniform 300;
+        value           uniform 300;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    farFieldWall
+    {
+        type            zeroGradient;
+    }
+    sphere
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/U

@@ -0,0 +1,64 @@
+/*--------------------------------*- 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;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform (0 0 0);
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform (0 0 0);
+    }
+    top
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    farFieldWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    sphere
+    {
+        type            movingWallVelocity;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/alpha.water

@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus.master                           |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.walter;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 0;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 0;
+    }
+    top
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    farFieldWall
+    {
+        type            zeroGradient;
+    }
+    sphere
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/p

@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus.master                           |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 101325;
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 101325;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 101325;
+    }
+    top
+    {
+        type            calculated;
+        value           uniform 101325;
+    }
+    bottom
+    {
+        type            calculated;
+        value           uniform 101325;
+    }
+    farFieldWall
+    {
+        type            calculated;
+        value           uniform 101325;
+    }
+    sphere
+    {
+        type            calculated;
+        value           uniform 101325;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/p_rgh

@@ -0,0 +1,68 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus.master                           |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 101325;
+
+boundaryField
+{
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 101325;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+        value           uniform 101325;
+    }
+    top
+    {
+        type            fixedValue;
+        value           uniform 101325;
+    }
+    bottom
+    {
+        type            fixedFluxPressure;
+        gradient        uniform 0;
+        value           uniform 101325;
+    }
+    farFieldWall
+    {
+        type            fixedFluxPressure;
+        gradient        uniform 0;
+        value           uniform 101325;
+    }
+    sphere
+    {
+        type            fixedFluxPressure;
+        gradient        uniform 0;
+        value           uniform 101325;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/0.org/pointDisplacement

@@ -0,0 +1,99 @@
+/*--------------------------------*- 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       pointVectorField;
+    location    "0";
+    object      pointDisplacement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "$FOAM_CASE/constant/dynamicMeshDict"
+#include "$FOAM_CASE/system/blockMeshDict"
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    farFieldWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    bottom
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    top
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    sphere
+    {
+        type            sixDoFRigidBodyDisplacement;
+        centreOfMass    (0 $:sphereCentreHeight 0);
+        momentOfInertia (2e-3 2e-3 2e-3);
+        velocity        (0 -3.77 0);
+        accelerationRelaxation  1;
+        accelerationDamping     1;
+
+        mass            2.08e-3;  // sphere mass 0.15 kg, 5 degree wedge
+//      rhoInf          -1;       // needed only for solvers solving for
+                                  // kinematic pressure
+        report          on;
+
+        constraints
+        {
+            maxIterations   500;
+
+            fixedLine1
+            {
+                sixDoFRigidBodyMotionConstraint line;
+                centreOfRotation (0 0.147 0);
+                direction       (0 1 0);
+            }
+
+            fixedOrientation1
+            {
+                sixDoFRigidBodyMotionConstraint orientation;
+            }
+        }
+
+        solver
+        {
+            type Newmark;
+        }
+
+        value           uniform ( 0 0 0 );
+    }
+    innerAMI
+    {
+        type            cyclicAMI;
+    }
+    outerAMI
+    {
+        type            cyclicAMI;
+    }
+    back
+    {
+        type            wedge;
+    }
+    front
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/Allclean

@@ -0,0 +1,7 @@
+#!/bin/sh
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm -rf 0/
+cleanCase

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/Allrun

@@ -0,0 +1,25 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication topoSet
+
+# remove empty patches
+runApplication createPatch -overwrite
+
+# split mesh into 2 regions at the AMI coupled patch
+runApplication mergeOrSplitBaffles -overwrite
+
+echo "Copying 0.org to 0"
+cp -r 0.org 0
+
+runApplication setFields
+
+(cd constant/polyMesh && ln -s ../../include/meshModifiers .)
+
+runApplication $(getApplication)
+
+# ----------------------------------------------------------------- end-of-file

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/dynamicMeshDict

@@ -0,0 +1,108 @@
+/*--------------------------------*- 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       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicMotionSolverTopoFvMesh;
+
+motionSolverLibs ("libfvMotionSolvers.so");
+
+solver          displacementLayeredMotion;
+
+displacementLayeredMotionCoeffs
+{
+    // In addition to boundary conditions on pointDisplacement
+    // we can also have boundary conditions on faceZones bordering cellZones.
+    regions
+    {
+        LARTopCells
+        {
+            interpolationScheme oneSided;
+
+            boundaryField
+            {
+                top // patch
+                {
+                    type            uniformFollow; // use pointDisplacement bc
+                    patch           sphere;          // fixed/static
+                }
+
+                LARTopFaces
+                {
+                    type            slip;
+                }
+            }
+        }
+
+        topBlock
+        {
+            interpolationScheme oneSided;
+
+            boundaryField
+            {
+                LARTopFaces
+                {
+                    type            uniformFollow;
+                    patch           sphere;
+                }
+
+                centralTopFaces
+                {
+                    type            slip;
+                }
+
+            }
+        }
+
+        centralBlock
+        {
+            interpolationScheme oneSided;
+
+            boundaryField
+            {
+                centralTopFaces
+                {
+                    type            uniformFollow;
+                    patch           sphere;
+                }
+
+                centralBottomFaces
+                {
+                    type            slip;
+                }
+            }
+        }
+
+        bottomBlock
+        {
+            interpolationScheme oneSided;
+
+            boundaryField
+            {
+                centralBottomFaces
+                {
+                    type            uniformFollow;
+                    patch           sphere;
+                }
+                bottom
+                {
+                    type            slip;
+                }
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/g

@@ -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       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0 );
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/thermophysicalProperties

@@ -0,0 +1,24 @@
+/*--------------------------------*- 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      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+pMin            pMin [1 -1 -2 0 0 0 0] 1000;
+
+sigma           sigma [1 0 -2 0 0 0 0] 0.07;
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/thermophysicalProperties.air

@@ -0,0 +1,55 @@
+/*--------------------------------*- 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      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    equationOfState
+    {
+        rho         1.27;
+        R           90000;
+        rho0        1.27;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/thermophysicalProperties.water

@@ -0,0 +1,53 @@
+/*--------------------------------*- 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      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18.0;
+    }
+    equationOfState
+    {
+        rho         1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/transportProperties

@@ -0,0 +1,33 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+phases (water air);
+
+water
+{
+    transportModel  Newtonian;
+    nu              [0 2 -1 0 0 0 0] 1e-06;
+    rho             [1 -3 0 0 0 0 0] 998.2;
+}
+
+air
+{
+    transportModel  Newtonian;
+    nu              [0 2 -1 0 0 0 0] 1.48e-05;
+    rho             [1 -3 0 0 0 0 0] 1;
+}
+
+// ************************************************************************* //

tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/constant/turbulenceProperties

@@ -0,0 +1,28 @@
+/*--------------------------------*- 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;
+
+RAS
+{
+    RASModel        laminar;
+
+    turbulence      on;
+
+    printCoeffs     off;
+}
+// ************************************************************************* //