diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index e8e4f79e4e..b4bbed1442 100644
--- a/src/functionObjects/field/Make/files
+++ b/src/functionObjects/field/Make/files
@@ -72,6 +72,8 @@ externalCoupled/externalCoupled.C
 externalCoupled/externalCoupledMixed/externalCoupledMixedFvPatchFields.C
 externalCoupled/externalCoupledTemperatureMixed/externalCoupledTemperatureMixedFvPatchScalarField.C
 
+extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
+
 ddt2/ddt2.C
 zeroGradient/zeroGradient.C
 
diff --git a/src/functionObjects/field/Make/options b/src/functionObjects/field/Make/options
index ea32eb29ff..efaea7c1b8 100644
--- a/src/functionObjects/field/Make/options
+++ b/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
diff --git a/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.C b/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.C
new file mode 100644
index 0000000000..8bc18a3ce3
--- /dev/null
+++ b/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.C
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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()
+:
+    globalFaceIHit(-1),
+    VC(vector::zero),
+    VU(vector::zero),
+    V(0),
+    time(0),
+    timeIndex(0)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const eulerianParticle& p)
+{
+    os  << p.globalFaceIHit << token::SPACE
+        << p.VC << token::SPACE
+        << p.VU << token::SPACE
+        << p.V << token::SPACE
+        << p.time << token::SPACE
+        << p.timeIndex;
+
+    return os;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, eulerianParticle& p)
+{
+    is  >> p.globalFaceIHit
+        >> p.VC
+        >> p.VU
+        >> p.V
+        >> p.time
+        >> p.timeIndex;
+
+    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
+        << globalFaceIHit << 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", globalFaceIHit);
+    dict.add("position", C);
+    dict.add("diameter", pDiameter);
+    dict.add("U", U);
+
+    return dict;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.H b/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.H
new file mode 100644
index 0000000000..e5d4bfc2b4
--- /dev/null
+++ b/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticle.H
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 globalFaceIHit;
+
+        //- 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;
+
+        //- Index of last output time
+        label timeIndex;
+
+
+    //- 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.globalFaceIHit == b.globalFaceIHit
+             && a.VC == b.VC
+             && a.VU == b.VU
+             && a.V == b.V
+             && a.time == b.time
+             && a.timeIndex == b.timeIndex;
+        }
+
+        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
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticleTemplates.C b/src/functionObjects/field/extractEulerianParticles/eulerianParticle/eulerianParticleTemplates.C
new file mode 100644
index 0000000000..d0fbee15fe
--- /dev/null
+++ b/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.globalFaceIHit != -1) && (p1.globalFaceIHit == -1))
+        {
+            return p0;
+        }
+        else if ((p0.globalFaceIHit == -1) && (p1.globalFaceIHit != -1))
+        {
+            return p1;
+        }
+        else if ((p0.globalFaceIHit != -1) && (p1.globalFaceIHit != -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
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
new file mode 100644
index 0000000000..60821e8092
--- /dev/null
+++ b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.C
@@ -0,0 +1,893 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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(".."); //  = obr_.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: "
+            << 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);
+    }
+
+    Log << 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);
+
+    // Create global addressing for faceZone
+    globalFaces_ = globalIndex(fz.size());
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::initialiseBins()
+{
+    DebugInFunction << endl;
+
+    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, 10, 50, nLocations, labelMax, 180);
+    ppa.agglomerate();
+
+    label nCoarseFaces = 0;
+    if (nFaces != 0)
+    {
+        fineToCoarseAddr_ = ppa.restrictTopBottomAddressing();
+        nCoarseFaces = max(fineToCoarseAddr_) + 1;
+        coarseToFineAddr_ = invertOneToMany(nCoarseFaces, fineToCoarseAddr_);
+
+        // Set coarse face centres as area average of fine face centres
+        const vectorField& faceCentres = mesh_.faceCentres();
+        const vectorField& faceAreas = mesh_.faceAreas();
+        coarsePosition_.setSize(coarseToFineAddr_.size());
+        forAll(coarseToFineAddr_, coarsei)
+        {
+            const labelList& fineFaces = coarseToFineAddr_[coarsei];
+
+            scalar sumArea = 0;
+            vector averagePosition(vector::zero);
+            forAll(fineFaces, i)
+            {
+                label facei = fz[fineFaces[i]];
+                scalar magSf = mag(faceAreas[facei]);
+                sumArea += magSf;
+                averagePosition += magSf*faceCentres[facei];
+            }
+            coarsePosition_[coarsei] = averagePosition/sumArea;
+        }
+    }
+
+    // Create global addressing for coarse face addressing
+    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 =
+            obr_.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;
+        if (iter != Map<label>::end())
+        {
+            // Particle on local processor
+            p = particles_[iter()];
+        }
+        reduce(p, sumParticleOp<eulerianParticle>());
+
+        const scalar pDiameter = cbrt(6.0*p.V/constant::mathematical::pi);
+
+        if ((pDiameter > minDiameter_) && (pDiameter < maxDiameter_))
+        {
+            p.time = time;
+            p.timeIndex = outputTimes_.size() - 1;
+
+            if (Pstream::master() && writeRawData_)
+            {
+                rawParticlesDictPtr_->add
+                (
+                    word("particle" + Foam::name(nCollectedParticles_)),
+                    p.writeDict()
+                );
+            }
+
+            if (globalFaces_.isLocal(p.globalFaceIHit))
+            {
+                collectedParticles_.append(p);
+            }
+
+            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.globalFaceIHit < 0)
+            {
+                // New particle - does not exist in particles_ list
+                p.globalFaceIHit = globalFaces_.toGlobal(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;
+        }
+    }
+}
+
+
+void Foam::functionObjects::extractEulerianParticles::writeBinnedParticleData()
+{
+    DebugInFunction << endl;
+
+    if (!createDistribution_)
+    {
+        // No need to store collected particles if not creating a distribution
+        collectedParticles_.clear();
+        return;
+    }
+
+    // Gather particles ready for collection from all procs
+    List<List<eulerianParticle> > allProcParticles(Pstream::nProcs());
+    allProcParticles[Pstream::myProcNo()] = collectedParticles_;
+    Pstream::gatherList(allProcParticles);
+    Pstream::scatterList(allProcParticles);
+    List<eulerianParticle> allParticles =
+        ListListOps::combine<List<eulerianParticle> >
+        (
+            allProcParticles,
+            accessOp<List<eulerianParticle> >()
+        );
+
+
+    // Determine coarse face index (global) and position for each particle
+    label nCoarseFaces = globalCoarseFaces_.size();
+    List<label> particleCoarseFacei(allParticles.size(), -1);
+    List<point> particleCoarseFacePosition(nCoarseFaces, point::min);
+
+    forAll(allParticles, particlei)
+    {
+        const eulerianParticle& p = allParticles[particlei];
+        label globalFaceHiti = p.globalFaceIHit;
+
+        if (globalFaces_.isLocal(globalFaceHiti))
+        {
+            label localFacei = globalFaces_.toLocal(globalFaceHiti);
+            label coarseFacei = fineToCoarseAddr_[localFacei];
+            label globalCoarseFacei = globalCoarseFaces_.toGlobal(coarseFacei);
+
+            particleCoarseFacei[particlei] = globalCoarseFacei;
+            particleCoarseFacePosition[globalCoarseFacei] =
+                coarsePosition_[coarseFacei];
+        }
+    }
+    Pstream::listCombineGather(particleCoarseFacei, maxEqOp<label>());
+    Pstream::listCombineGather(particleCoarseFacePosition, maxEqOp<point>());
+
+    // Write the agglomerated particle data to file
+    DynamicList<label> processedCoarseFaces;
+    if (Pstream::master())
+    {
+        fileName baseDir(dictBaseFileDir()/name());
+
+        IOdictionary dict
+        (
+            IOobject
+            (
+                "particleDistribution",
+                obr_.time().timeName(),
+                baseDir,
+                obr_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            )
+        );
+
+        labelListList coarseFaceToParticle =
+            invertOneToMany(nCoarseFaces, particleCoarseFacei);
+
+        // Process the allParticles per coarse face
+        forAll(coarseFaceToParticle, globalCoarseFacei)
+        {
+            const List<label>& particleIDs =
+                coarseFaceToParticle[globalCoarseFacei];
+
+            const label nParticle = particleIDs.size();
+
+            if (nParticle == 0)
+            {
+                continue;
+            }
+
+            Field<scalar> pd(particleIDs.size());
+            scalar sumV = 0;
+            vector sumVU = vector::zero;
+            scalar startTime = GREAT;
+            scalar endTime = -GREAT;
+            forAll(particleIDs, i)
+            {
+                const label particlei = particleIDs[i];
+                const eulerianParticle& p = allParticles[particlei];
+                scalar pDiameter = cbrt(6*p.V/constant::mathematical::pi);
+                pd[i] = pDiameter;
+                sumV += p.V;
+                sumVU += p.VU;
+
+                startTime = min(startTime, outputTimes_[p.timeIndex]);
+                endTime = max(endTime, outputTimes_[p.timeIndex + 1]);
+            }
+
+            if (sumV < ROOTVSMALL)
+            {
+                // Started collecting particle info, but not accumulated any
+                // volume yet
+                continue;
+            }
+
+
+            distributionModels::binned binnedDiameters
+            (
+                pd,
+                distributionBinWidth_,
+                rndGen_
+            );
+
+            // Velocity info hard-coded to volume average
+            vector Uave = sumVU/sumV;
+
+            dictionary particleDict;
+            particleDict.add("startTime", startTime);
+            particleDict.add("endTime", endTime);
+            particleDict.add("nParticle", nParticle);
+            particleDict.add
+            (
+                "position",
+                particleCoarseFacePosition[globalCoarseFacei]
+            );
+            particleDict.add("volume", sumV);
+            particleDict.add("U", Uave);
+            particleDict.add
+            (
+                "binnedDistribution",
+                binnedDiameters.writeDict("distribution")
+            );
+            dict.add
+            (
+                word("sample" + Foam::name(globalCoarseFacei)),
+                particleDict
+            );
+
+            processedCoarseFaces.append(globalCoarseFacei);
+        }
+
+        dict.regIOobject::write();
+    }
+
+
+    if (resetDistributionOnWrite_)
+    {
+        // Remove processed coarse faces from collectedParticles_
+        Pstream::scatter(processedCoarseFaces);
+        labelHashSet processedFaces(processedCoarseFaces);
+        DynamicList<eulerianParticle> nonProcessedParticles;
+        forAll(collectedParticles_, particlei)
+        {
+            const eulerianParticle& p = collectedParticles_[particlei];
+            label localFacei = globalFaces_.toLocal(p.globalFaceIHit);
+            label coarseFacei = fineToCoarseAddr_[localFacei];
+            label globalCoarseFacei = globalCoarseFaces_.toGlobal(coarseFacei);
+            if (!processedFaces.found(globalCoarseFacei))
+            {
+                nonProcessedParticles.append(p);
+            }
+        }
+        collectedParticles_.transfer(nonProcessedParticles);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::extractEulerianParticles::extractEulerianParticles
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(runTime, name),
+    writeRawData_(false),
+    rawParticlesDictPtr_(NULL),
+    outputTimes_(100),
+    faceZoneName_(word::null),
+    zoneID_(-1),
+    patchIDs_(),
+    patchFaceIDs_(),
+    alphaName_("alpha"),
+    alphaThreshold_(0.1),
+    UName_("U"),
+    rhoName_("rho"),
+    phiName_("phi"),
+    nInjectorLocations_(0),
+    fineToCoarseAddr_(),
+    coarseToFineAddr_(),
+    coarsePosition_(),
+    globalCoarseFaces_(),
+    regions0_(),
+    nRegions0_(0),
+    particles_(),
+    regionToParticleMap_(),
+    collectedParticles_(),
+    minDiameter_(ROOTVSMALL),
+    maxDiameter_(GREAT),
+    createDistribution_(false),
+    resetDistributionOnWrite_(false),
+    distributionBinWidth_(0),
+    rndGen_(1234, -1),
+    nCollectedParticles_(0),
+    nDiscardedParticles_(0),
+    discardedVolume_(0)
+{
+    if (mesh_.nSolutionD() != 3)
+    {
+        FatalErrorInFunction
+            << name << " function object only applicable to 3-D cases"
+            << exit(FatalError);
+    }
+
+    read(dict);
+
+    if (Pstream::master() && writeRawData_)
+    {
+        fileName baseDir(dictBaseFileDir()/name);
+
+        rawParticlesDictPtr_.reset
+        (
+            new IOdictionary
+            (
+                IOobject
+                (
+                    "rawParticlesDict",
+                    mesh_.time().timeName(),
+                    baseDir,
+                    mesh_,
+                    IOobject::NO_READ, // _IF_PRESENT,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
+DebugVar(rawParticlesDictPtr_->objectPath());
+DebugVar(rawParticlesDictPtr_->path());
+    }
+
+    outputTimes_.append(mesh_.time().value());
+}
+
+
+// * * * * * * * * * * * * * * * * 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("writeRawData") >> writeRawData_;
+
+        dict.lookup("faceZone") >> faceZoneName_;
+        dict.lookup("nLocations") >> nInjectorLocations_;
+        dict.lookup("alphaName") >> alphaName_;
+        dict.readIfPresent("alphaThreshold", alphaThreshold_);
+        dict.lookup("UName") >> UName_;
+        dict.lookup("rhoName") >> rhoName_;
+        dict.lookup("phiName") >> phiName_;
+
+        dict.readIfPresent("minDiameter", minDiameter_);
+        dict.readIfPresent("maxDiameter", maxDiameter_);
+        dict.lookup("createDistribution") >> createDistribution_;
+        dict.lookup("distributionBinWidth") >> distributionBinWidth_;
+        dict.lookup("resetDistributionOnWrite") >> resetDistributionOnWrite_;
+
+        checkFaceZone();
+
+        initialiseBins();
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::extractEulerianParticles::execute()
+{
+    DebugInFunction << endl;
+
+    Log << type() << " " << name() << " output:" << nl;
+
+    const volScalarField& alpha =
+        obr_.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);
+    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 fot 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;
+
+    outputTimes_.append(obr_.time().value());
+
+    writeBinnedParticleData();
+
+    if (Pstream::master() && writeRawData_)
+    {
+        rawParticlesDictPtr_->regIOobject::write();
+    }
+
+    return true;
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H
new file mode 100644
index 0000000000..5cf7d1cd69
--- /dev/null
+++ b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H
@@ -0,0 +1,346 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    Particle data is written to the directory:
+
+    \verbatim
+        $FOAM_CASE/postProcessing/<name>/rawParticlesDict
+        $FOAM_CASE/postProcessing/<name>/particleDistribution
+    \endverbatim
+
+Usage
+    extractEulerianParticles1
+    {
+        type            extractEulerianParticles;
+        libs            ("libfieldFunctionObjects.so");
+        ...
+        faceZone        f0;
+        nLocations      10;
+        alphaName       alpha.water;
+        UName           U;
+        rhoName         rho;
+        phiName         phi;
+
+        createDistribution yes;
+        distributionBinWidth 1e-4;
+        resetDistributionOnWrite no;
+
+        writeRawData    yes;
+    }
+    \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 |
+        createDistribution | Flag to create a binned distribution | yes |
+        distributionBinWidth | Binned distribution bin width| yes |
+        writeRawData | Flag to write raw particle data | 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"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+/*---------------------------------------------------------------------------*\
+         Class extractEulerianParticlesFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class extractEulerianParticles
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+
+protected:
+
+    // Protected data
+
+        //- Flag to write the raw collected particle data
+        bool writeRawData_;
+
+        //- Raw particle output dictionary
+        autoPtr<IOdictionary> rawParticlesDictPtr_;
+
+        //- List to keep track of output times
+        DynamicList<scalar> outputTimes_;
+
+
+        // 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_;
+
+            //- Global face addressing
+            globalIndex globalFaces_;
+
+
+        // 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, defauls = 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
+            labelList fineToCoarseAddr_;
+
+            //- Agglomeration addressing from coarse to fine
+            labelListList coarseToFineAddr_;
+
+            //- Coarse face positions
+            vectorList coarsePosition_;
+
+            //- Global coarse face addressing
+            globalIndex globalCoarseFaces_;
+
+
+        // Particle collection info
+
+            //- Time of last write step
+            scalar lastOutputTime_;
+
+            //- 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_;
+
+            // Collected particles
+            DynamicList<eulerianParticle> collectedParticles_;
+
+            //- 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_;
+
+            //- Flag to create a distribution
+            bool createDistribution_;
+
+            //- Flag to reset the distribution on each write
+            bool resetDistributionOnWrite_;
+
+            //- Diameter distribution bin width
+            scalar distributionBinWidth_;
+
+            //- Random class needed by distribution models
+            cachedRandom rndGen_;
+
+
+        // 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
+        );
+
+        //- Write agglomerated particle data to stream
+        virtual void writeBinnedParticleData();
+
+        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
+
+        //- Static data staticData
+        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
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticlesTemplates.C b/src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticlesTemplates.C
new file mode 100644
index 0000000000..e232098178
--- /dev/null
+++ b/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;
+        }
+    }
+}
+
+
+// ************************************************************************* //