diff --git a/applications/solvers/dsmc/dsmcFoam/Make/files b/applications/solvers/dsmc/dsmcFoam/Make/files
new file mode 100755
index 0000000000..a6860410e7
--- /dev/null
+++ b/applications/solvers/dsmc/dsmcFoam/Make/files
@@ -0,0 +1,3 @@
+dsmcFoam.C
+
+EXE = $(FOAM_APPBIN)/dsmcFoam
diff --git a/applications/solvers/dsmc/dsmcFoam/Make/options b/applications/solvers/dsmc/dsmcFoam/Make/options
new file mode 100755
index 0000000000..bc99834af6
--- /dev/null
+++ b/applications/solvers/dsmc/dsmcFoam/Make/options
@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -ldsmc
+
diff --git a/applications/solvers/dsmc/dsmcFoam/createFields.H b/applications/solvers/dsmc/dsmcFoam/createFields.H
new file mode 100644
index 0000000000..d024bd2017
--- /dev/null
+++ b/applications/solvers/dsmc/dsmcFoam/createFields.H
@@ -0,0 +1,162 @@
+
+    Info<< nl << "Reading field boundaryT" << endl;
+    volScalarField boundaryT
+    (
+        IOobject
+        (
+            "boundaryT",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field boundaryU" << endl;
+    volVectorField boundaryU
+    (
+        IOobject
+        (
+            "boundaryU",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoN (number density)" << endl;
+    volScalarField rhoN
+    (
+        IOobject
+        (
+            "rhoN",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoM (mass density)" << endl;
+    volScalarField rhoM
+    (
+        IOobject
+        (
+            "rhoM",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field rhoNdsmc (dsmc particle density)" << endl;
+    volScalarField dsmcRhoN
+    (
+        IOobject
+        (
+            "dsmcRhoN",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field momentum (momentum density)" << endl;
+    volVectorField momentum
+    (
+        IOobject
+        (
+            "momentum",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field linearKE (linear kinetic energy density)"
+        << endl;
+
+    volScalarField linearKE
+    (
+        IOobject
+        (
+            "linearKE",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field internalE (internal energy density)" << endl;
+    volScalarField internalE
+    (
+        IOobject
+        (
+            "internalE",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field iDof (internal degree of freedom density)"
+        << endl;
+
+    volScalarField iDof
+    (
+        IOobject
+        (
+            "iDof",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field q (surface heat transfer)" << endl;
+    volScalarField q
+    (
+        IOobject
+        (
+            "q",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field fD (surface force density)" << endl;
+    volVectorField fD
+    (
+        IOobject
+        (
+            "fD",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Constructing dsmcCloud " << endl;
+
+    dsmcCloud dsmc("dsmc", boundaryT, boundaryU);
diff --git a/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
new file mode 100644
index 0000000000..1ea52cc6aa
--- /dev/null
+++ b/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    dsmcFoam
+
+Description
+    Direct Simulation Monte Carlo Solver for 3D, transient, multi-species flows
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Carry out dsmcCloud timestep
+
+        dsmc.evolve();
+
+        // Retrieve flow field data from dsmcCloud
+
+        rhoN = dsmc.rhoN();
+        rhoN.correctBoundaryConditions();
+
+        rhoM = dsmc.rhoM();
+        rhoM.correctBoundaryConditions();
+
+        dsmcRhoN = dsmc.dsmcRhoN();
+        dsmcRhoN.correctBoundaryConditions();
+
+        momentum = dsmc.momentum();
+        momentum.correctBoundaryConditions();
+
+        linearKE = dsmc.linearKE();
+        linearKE.correctBoundaryConditions();
+
+        internalE = dsmc.internalE();
+        internalE.correctBoundaryConditions();
+
+        iDof = dsmc.iDof();
+        iDof.correctBoundaryConditions();
+
+        // Retrieve surface field data from dsmcCloud
+
+        q = dsmc.q();
+
+        fD = dsmc.fD();
+
+        // Print status of dsmcCloud
+
+        dsmc.info();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/files b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/files
new file mode 100755
index 0000000000..5cba9c270d
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/files
@@ -0,0 +1,3 @@
+dsmcFieldsCalc.C
+
+EXE = $(FOAM_APPBIN)/dsmcFieldsCalc
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/options b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/options
new file mode 100755
index 0000000000..4e81bb55de
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/options
@@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(LIB_SRC)/postProcessing/postCalc \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/postProcessing/functionObjects/utilities/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    $(FOAM_LIBBIN)/postCalc.o \
+    -lmeshTools \
+    -lfiniteVolume \
+    -lutilityFunctionObjects \
+    -llagrangian \
+    -ldsmc
+
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C
new file mode 100644
index 0000000000..ee4fd2c50f
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    dsmcFields
+
+Description
+    Calculate intensive fields (U and T) from averaged extensive fields from a
+    DSMC calculation.
+
+\*---------------------------------------------------------------------------*/
+
+#include "calc.H"
+#include "fvc.H"
+#include "dsmcCloud.H"
+#include "dsmcFields.H"
+#include "IOobjectList.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    template<class Type>
+    bool addFieldsToList
+    (
+        const fvMesh& mesh,
+        PtrList<GeometricField<Type, fvPatchField, volMesh> >& list,
+        const wordList& fieldNames
+    )
+    {
+        typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+        label index = 0;
+        forAll(fieldNames, i)
+        {
+            IOobject obj
+            (
+                fieldNames[i],
+                mesh.time().timeName(),
+                mesh,
+                IOobject::MUST_READ
+            );
+
+            if (obj.headerOk() && obj.headerClassName() == fieldType::typeName)
+            {
+                list.set(index++, new fieldType(obj, mesh));
+            }
+            else
+            {
+                Info<< "Could not find " << fieldNames[i] << endl;
+
+                return false;
+            }
+        }
+
+        return true;
+    }
+}
+
+
+void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
+{
+    bool writeResults = !args.options().found("noWrite");
+
+    wordList extensiveVSFNames
+    (
+        IStringStream
+        (
+            "( \
+                rhoNMean \
+                rhoMMean \
+                linearKEMean \
+                internalEMean \
+                iDofMean \
+            )"
+        )()
+    );
+
+    PtrList<volScalarField> extensiveVSFs(extensiveVSFNames.size());
+
+    if
+    (
+        !addFieldsToList
+        (
+            mesh,
+            extensiveVSFs,
+            extensiveVSFNames
+        )
+    )
+    {
+        return;
+    }
+
+    wordList extensiveVVFNames(IStringStream ("(momentumMean)")());
+
+    PtrList<volVectorField> extensiveVVFs(extensiveVVFNames.size());
+
+    if
+    (
+        !addFieldsToList
+        (
+            mesh,
+            extensiveVVFs,
+            extensiveVVFNames
+        )
+    )
+    {
+        return;
+    }
+
+    dsmcFields dF
+    (
+        "dsmcFieldsUtility",
+        mesh,
+        dictionary::null,
+        false
+    );
+
+    if (writeResults)
+    {
+        dF.write();
+    }
+}
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/dsmcInitialise/Make/files b/applications/utilities/preProcessing/dsmcInitialise/Make/files
new file mode 100755
index 0000000000..7853cb2bd6
--- /dev/null
+++ b/applications/utilities/preProcessing/dsmcInitialise/Make/files
@@ -0,0 +1,3 @@
+dsmcInitialise.C
+
+EXE = $(FOAM_APPBIN)/dsmcInitialise
diff --git a/applications/utilities/preProcessing/dsmcInitialise/Make/options b/applications/utilities/preProcessing/dsmcInitialise/Make/options
new file mode 100755
index 0000000000..bc99834af6
--- /dev/null
+++ b/applications/utilities/preProcessing/dsmcInitialise/Make/options
@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -ldsmc
+
diff --git a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
new file mode 100644
index 0000000000..c2f797c17e
--- /dev/null
+++ b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    dsmcFoam
+
+Description
+    Initialise a case for dsmcFoam by reading the initialisation dictionary
+    system/dsmcInitialise
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "Initialising dsmc for Time = " << runTime.timeName() << nl << endl;
+
+    dsmcCloud dsmc("dsmc", mesh);
+
+    label totalMolecules = dsmc.size();
+
+    if (Pstream::parRun())
+    {
+        reduce(totalMolecules, sumOp<label>());
+    }
+
+    Info<< nl << "Total number of molecules added: " << totalMolecules
+        << nl << endl;
+
+    IOstream::defaultPrecision(15);
+
+    if (!mesh.write())
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing dsmcCloud."
+            << nl << exit(FatalError);
+    }
+
+    Info<< nl << "ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/Allwmake b/src/lagrangian/Allwmake
index 41b2fd9bf9..fea6fce12e 100755
--- a/src/lagrangian/Allwmake
+++ b/src/lagrangian/Allwmake
@@ -6,6 +6,7 @@ wmake libso basic
 wmake libso solidParticle
 wmake libso intermediate
 wmake libso dieselSpray
+wmake libso dsmc
 
 molecularDynamics/Allwmake
 
diff --git a/src/lagrangian/dsmc/Make/files b/src/lagrangian/dsmc/Make/files
new file mode 100644
index 0000000000..9f6d33d209
--- /dev/null
+++ b/src/lagrangian/dsmc/Make/files
@@ -0,0 +1,16 @@
+/* Parcels */
+parcels/derived/dsmcParcel/dsmcParcel.C
+
+/* Cloud base classes */
+clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C
+
+/* Clouds */
+clouds/derived/dsmcCloud/dsmcCloud.C
+
+/* submodels */
+parcels/derived/dsmcParcel/defineDsmcParcel.C
+parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C
+parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
+parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C
+
+LIB = $(FOAM_LIBBIN)/libdsmc
diff --git a/src/lagrangian/dsmc/Make/options b/src/lagrangian/dsmc/Make/options
new file mode 100644
index 0000000000..15874b7b55
--- /dev/null
+++ b/src/lagrangian/dsmc/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude
+
+LIB_LIBS = \
+    -llagrangian \
+    -lfiniteVolume
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
new file mode 100644
index 0000000000..4c4bcd0a43
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -0,0 +1,840 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "DsmcCloud.H"
+#include "BinaryCollisionModel.H"
+#include "WallInteractionModel.H"
+#include "InflowBoundaryModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
+
+template<class ParcelType>
+Foam::scalar Foam::DsmcCloud<ParcelType>::Tref = 273;
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::buildConstProps()
+{
+    Info<< nl << "Constructing constant properties for" << endl;
+    constProps_.setSize(typeIdList_.size());
+
+    dictionary moleculeProperties
+    (
+        particleProperties_.subDict("moleculeProperties")
+    );
+
+    forAll(typeIdList_, i)
+    {
+        const word& id(typeIdList_[i]);
+
+        Info<< "    " << id << endl;
+
+        const dictionary& molDict(moleculeProperties.subDict(id));
+
+        constProps_[i] =
+        typename ParcelType::constantProperties::constantProperties(molDict);
+    }
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::buildCellOccupancy()
+{
+    forAll(cellOccupancy_, cO)
+    {
+        cellOccupancy_[cO].clear();
+    }
+
+    forAllIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        cellOccupancy_[iter().cell()].append(&iter());
+    }
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::initialise
+(
+    const IOdictionary& dsmcInitialiseDict
+)
+{
+    Info<< nl << "Initialising particles" << endl;
+
+    const scalar temperature
+    (
+        readScalar(dsmcInitialiseDict.lookup("temperature"))
+    );
+
+    const vector velocity(dsmcInitialiseDict.lookup("velocity"));
+
+    const dictionary& numberDensitiesDict
+    (
+        dsmcInitialiseDict.subDict("numberDensities")
+    );
+
+    List<word> molecules(numberDensitiesDict.toc());
+
+    Field<scalar> numberDensities(molecules.size());
+
+    forAll(molecules, i)
+    {
+        numberDensities[i] = readScalar
+        (
+            numberDensitiesDict.lookup(molecules[i])
+        );
+    }
+
+    numberDensities /= nParticle_;
+
+    scalar x0 = mesh_.bounds().min().x();
+    scalar xR = mesh_.bounds().max().x() - x0;
+
+    scalar y0 = mesh_.bounds().min().y();
+    scalar yR = mesh_.bounds().max().y() - y0;
+
+    scalar z0 = mesh_.bounds().min().z();
+    scalar zR = mesh_.bounds().max().z() - z0;
+
+    forAll(molecules, i)
+    {
+        const word& moleculeName(molecules[i]);
+
+        label typeId(findIndex(typeIdList_, moleculeName));
+
+        if (typeId == -1)
+        {
+            FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
+                << "typeId " << moleculeName << "not defined." << nl
+                << abort(FatalError);
+        }
+
+        const typename ParcelType::constantProperties& cP = constProps(typeId);
+
+        scalar numberDensity = numberDensities[i];
+
+        scalar spacing = pow(numberDensity,-(1.0/3.0));
+
+        int ni = label(xR/spacing) + 1;
+        int nj = label(yR/spacing) + 1;
+        int nk = label(zR/spacing) + 1;
+
+        vector delta(xR/ni, yR/nj, zR/nk);
+
+        scalar pert = spacing;
+
+        for (int i = 0; i < ni; i++)
+        {
+            for (int j = 0; j < nj; j++)
+            {
+                for (int k = 0; k < nk; k++)
+                {
+                    point p
+                    (
+                        x0 + (i + 0.5)*delta.x(),
+                        y0 + (j + 0.5)*delta.y(),
+                        z0 + (k + 0.5)*delta.z()
+                    );
+
+                    p.x() += pert*(rndGen_.scalar01() - 0.5);
+                    p.y() += pert*(rndGen_.scalar01() - 0.5);
+                    p.z() += pert*(rndGen_.scalar01() - 0.5);
+
+                    label cell = mesh_.findCell(p);
+
+                    vector U = equipartitionLinearVelocity
+                    (
+                        temperature,
+                        cP.mass()
+                    );
+
+                    scalar Ei = equipartitionInternalEnergy
+                    (
+                        temperature,
+                        cP.internalDegreesOfFreedom()
+                    );
+
+                    U += velocity;
+
+                    if (cell >= 0)
+                    {
+                        addNewParcel
+                        (
+                            p,
+                            U,
+                            Ei,
+                            cell,
+                            typeId
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    // Initialise the sigmaTcRMax_ field to the product of the cross section of
+    // the most abundant species and the most probable thermal speed (Bird,
+    // p222-223)
+
+    label mostAbundantType(findMax(numberDensities));
+
+    const typename ParcelType::constantProperties& cP = constProps
+    (
+        mostAbundantType
+    );
+
+    sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
+    (
+        temperature,
+        cP.mass()
+    );
+
+    sigmaTcRMax_.correctBoundaryConditions();
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::collisions()
+{
+    buildCellOccupancy();
+
+    // Temporary storage for subCells
+    List<DynamicList<label> > subCells(8);
+
+    scalar deltaT = mesh_.time().deltaT().value();
+
+    label collisionCandidates = 0;
+
+    label collisions = 0;
+
+    forAll(cellOccupancy_, celli)
+    {
+        const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
+
+        label nC(cellParcels.size());
+
+        if (nC > 1)
+        {
+
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            // Assign particles to one of 8 Cartesian subCells
+
+            // Clear temporary lists
+            forAll(subCells, i)
+            {
+                subCells[i].clear();
+            }
+
+            // Inverse addressing specifying which subCell a parcel is in
+            List<label> whichSubCell(cellParcels.size());
+
+            const point& cC = mesh_.cellCentres()[celli];
+
+            forAll(cellParcels, i)
+            {
+                ParcelType* p = cellParcels[i];
+
+                vector relPos = p->position() - cC;
+
+                label subCell =
+                    pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
+
+                subCells[subCell].append(i);
+
+                whichSubCell[i] = subCell;
+            }
+
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            scalar sigmaTcRMax = sigmaTcRMax_[celli];
+
+            scalar selectedPairs =
+                collisionSelectionRemainder_[celli]
+              + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT
+               /mesh_.cellVolumes()[celli];
+
+            label nCandidates(selectedPairs);
+
+            collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
+
+            collisionCandidates += nCandidates;
+
+            for (label c = 0; c < nCandidates; c++)
+            {
+                // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+                // subCell candidate selection procedure
+
+                // Select the first collision candidate
+                label candidateP = rndGen_.integer(0, nC - 1);
+
+                // Declare the second collision candidate
+                label candidateQ = -1;
+
+                List<label> subCellPs = subCells[whichSubCell[candidateP]];
+
+                label nSC = subCellPs.size();
+
+                if (nSC > 1)
+                {
+                    // If there are two or more particle in a subCell, choose
+                    // another from the same cell.  If the same candidate is
+                    // chosen, choose again.
+
+                    do
+                    {
+                        candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
+
+                    } while(candidateP == candidateQ);
+                }
+                else
+                {
+                    // Select a possible second collision candidate from the
+                    // whole cell.  If the same candidate is chosen, choose
+                    // again.
+
+                    do
+                    {
+                        candidateQ = rndGen_.integer(0, nC - 1);
+
+                    } while(candidateP == candidateQ);
+                }
+
+                // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+                // uniform candidate selection procedure
+
+                // // Select the first collision candidate
+                // label candidateP = rndGen_.integer(0, nC-1);
+
+                // // Select a possible second collision candidate
+                // label candidateQ = rndGen_.integer(0, nC-1);
+
+                // // If the same candidate is chosen, choose again
+                // while(candidateP == candidateQ)
+                // {
+                //     candidateQ = rndGen_.integer(0, nC-1);
+                // }
+
+                // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+                ParcelType* parcelP = cellParcels[candidateP];
+                ParcelType* parcelQ = cellParcels[candidateQ];
+
+                label typeIdP = parcelP->typeId();
+                label typeIdQ = parcelQ->typeId();
+
+                scalar sigmaTcR = binaryCollision().sigmaTcR
+                (
+                    typeIdP,
+                    typeIdQ,
+                    parcelP->U(),
+                    parcelQ->U()
+                );
+
+                // Update the maximum value of sigmaTcR stored, but use the
+                // initial value in the acceptance-rejection criteria because
+                // the number of collision candidates selected was based on this
+
+                if (sigmaTcR > sigmaTcRMax_[celli])
+                {
+                    sigmaTcRMax_[celli] = sigmaTcR;
+                }
+
+                if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
+                {
+                    binaryCollision().collide
+                    (
+                        typeIdP,
+                        typeIdQ,
+                        parcelP->U(),
+                        parcelQ->U(),
+                        parcelP->Ei(),
+                        parcelQ->Ei()
+                    );
+
+                    collisions++;
+                }
+            }
+        }
+    }
+
+    reduce(collisions, sumOp<label>());
+
+    reduce(collisionCandidates, sumOp<label>());
+
+    if (collisionCandidates)
+    {
+        Info<< "    Collisions                      = "
+            << collisions << nl
+            << "    Acceptance rate                 = "
+            << scalar(collisions)/scalar(collisionCandidates) << nl
+            << endl;
+    }
+    else
+    {
+        Info<< "    No collisions" << endl;
+    }
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::resetSurfaceDataFields()
+{
+    volScalarField::GeometricBoundaryField& qBF(q_.boundaryField());
+
+    forAll(qBF, p)
+    {
+        qBF[p] = 0.0;
+    }
+
+    volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField());
+
+    forAll(fDBF, p)
+    {
+        fDBF[p] = vector::zero;
+    }
+}
+
+
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const vector& U,
+    const scalar Ei,
+    const label cellId,
+    const label typeId
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        position,
+        U,
+        Ei,
+        cellId,
+        typeId
+    );
+
+    addParticle(pPtr);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::DsmcCloud<ParcelType>::DsmcCloud
+(
+    const word& cloudType,
+    const volScalarField& T,
+    const volVectorField& U
+)
+:
+    Cloud<ParcelType>(T.mesh(), cloudType, false),
+    DsmcBaseCloud(),
+    cloudType_(cloudType),
+    mesh_(T.mesh()),
+    particleProperties_
+    (
+        IOobject
+        (
+            cloudType + "Properties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    typeIdList_(particleProperties_.lookup("typeIdList")),
+    nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
+    cellOccupancy_(mesh_.nCells()),
+    sigmaTcRMax_
+    (
+        IOobject
+        (
+            this->name() + "SigmaTcRMax",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    collisionSelectionRemainder_(mesh_.nCells(), 0),
+    q_
+    (
+        IOobject
+        (
+            this->name() + "q_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
+    ),
+    fD_
+    (
+        IOobject
+        (
+            this->name() + "fD_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            dimensionSet(1, -1, -2, 0, 0),
+            vector::zero
+        )
+    ),
+    constProps_(),
+    rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
+    T_(T),
+    U_(U),
+    binaryCollisionModel_
+    (
+        BinaryCollisionModel<DsmcCloud<ParcelType> >::New
+        (
+            particleProperties_,
+            *this
+        )
+    ),
+    wallInteractionModel_
+    (
+        WallInteractionModel<DsmcCloud<ParcelType> >::New
+        (
+            particleProperties_,
+            *this
+        )
+    ),
+    inflowBoundaryModel_
+    (
+        InflowBoundaryModel<DsmcCloud<ParcelType> >::New
+        (
+            particleProperties_,
+            *this
+        )
+    )
+{
+    buildConstProps();
+
+    buildCellOccupancy();
+}
+
+
+template<class ParcelType>
+Foam::DsmcCloud<ParcelType>::DsmcCloud
+(
+    const word& cloudType,
+    const fvMesh& mesh
+)
+    :
+    Cloud<ParcelType>(mesh, cloudType, false),
+    DsmcBaseCloud(),
+    cloudType_(cloudType),
+    mesh_(mesh),
+    particleProperties_
+    (
+        IOobject
+        (
+            cloudType + "Properties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    typeIdList_(particleProperties_.lookup("typeIdList")),
+    nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
+    cellOccupancy_(),
+    sigmaTcRMax_
+    (
+        IOobject
+        (
+            this->name() + "SigmaTcRMax",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, 3, -1, 0, 0), 0.0)
+    ),
+    collisionSelectionRemainder_(),
+    q_
+    (
+        IOobject
+        (
+            this->name() + "q_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
+    ),
+    fD_
+    (
+        IOobject
+        (
+            this->name() + "fD_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            dimensionSet(1, -1, -2, 0, 0),
+            vector::zero
+        )
+    ),
+    constProps_(),
+    rndGen_(label(971501) + 1526*Pstream::myProcNo()),
+    T_
+    (
+        volScalarField
+        (
+            IOobject
+            (
+                "T",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
+        )
+    ),
+    U_
+    (
+        volVectorField
+        (
+            IOobject
+            (
+                "U",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector
+            (
+                "zero",
+                dimensionSet(0, 1, -1, 0, 0),
+                vector::zero
+            )
+        )
+    ),
+    binaryCollisionModel_(),
+    wallInteractionModel_(),
+    inflowBoundaryModel_()
+{
+    clear();
+
+    buildConstProps();
+
+    IOdictionary dsmcInitialiseDict
+    (
+        IOobject
+        (
+            "dsmcInitialiseDict",
+            mesh_.time().system(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    initialise(dsmcInitialiseDict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::DsmcCloud<ParcelType>::~DsmcCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::evolve()
+{
+    typename ParcelType::trackData td(*this);
+
+    // Reset the surface data collection fields
+    resetSurfaceDataFields();
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
+
+    // Insert new particles from the inflow boundary
+    this->inflowBoundary().inflow();
+
+    // Move the particles ballistically with their current velocities
+    Cloud<ParcelType>::move(td);
+
+    // Calculate new velocities via stochastic collisions
+    collisions();
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::info() const
+{
+    label nDsmcParticles = this->size();
+    reduce(nDsmcParticles, sumOp<label>());
+
+    scalar nMol = nDsmcParticles*nParticle_;
+
+    vector linearMomentum = linearMomentumOfSystem();
+    reduce(linearMomentum, sumOp<vector>());
+
+    scalar linearKineticEnergy = linearKineticEnergyOfSystem();
+    reduce(linearKineticEnergy, sumOp<scalar>());
+
+    scalar internalEnergy = internalEnergyOfSystem();
+    reduce(internalEnergy, sumOp<scalar>());
+
+    Info<< "Cloud name: " << this->name() << nl
+        << "    Number of dsmc particles        = "
+        << nDsmcParticles << nl
+        << "    Number of molecules             = "
+        << nMol << nl
+        << "    Mass in system                  = "
+        << returnReduce(massInSystem(), sumOp<scalar>()) << nl
+        << "    Average linear momentum         = "
+        << linearMomentum/nMol << nl
+        << "   |Average linear momentum|        = "
+        << mag(linearMomentum)/nMol << nl
+        << "    Average linear kinetic energy   = "
+        << linearKineticEnergy/nMol << nl
+        << "    Average internal energy         = "
+        << internalEnergy/nMol << nl
+        << "    Average total energy            = "
+        << (internalEnergy + linearKineticEnergy)/nMol << nl
+        << endl;
+}
+
+
+template<class ParcelType>
+Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
+(
+    scalar temperature,
+    scalar mass
+)
+{
+    return
+        sqrt(kb*temperature/mass)
+       *vector
+        (
+            rndGen_.GaussNormal(),
+            rndGen_.GaussNormal(),
+            rndGen_.GaussNormal()
+        );
+}
+
+
+template<class ParcelType>
+Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
+(
+    scalar temperature,
+    scalar iDof
+)
+{
+    scalar Ei = 0.0;
+
+    if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
+    {
+        // Special case for iDof = 2, i.e. diatomics;
+        Ei = -log(rndGen_.scalar01())*kb*temperature;
+    }
+    else
+    {
+        scalar a = 0.5*iDof - 1;
+
+        scalar energyRatio;
+
+        scalar P;
+
+        do
+        {
+            energyRatio = 10*rndGen_.scalar01();
+
+            P = pow((energyRatio/a), a)*exp(a - energyRatio);
+
+        } while (P < rndGen_.scalar01());
+
+        Ei = energyRatio*kb*temperature;
+    }
+
+    return Ei;
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::dumpParticlePositions() const
+{
+    OFstream pObj
+    (
+        this->db().time().path()/"parcelPositions_"
+      + this->name() + "_"
+      + this->db().time().timeName() + ".obj"
+    );
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+
+        pObj<< "v " << p.position().x()
+            << " "  << p.position().y()
+            << " "  << p.position().z()
+            << nl;
+    }
+
+    pObj.flush();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
new file mode 100644
index 0000000000..b246b04990
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -0,0 +1,434 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::DsmcCloud
+
+Description
+    Templated base class for dsmc cloud
+
+SourceFiles
+    DsmcCloudI.H
+    DsmcCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DsmcCloud_H
+#define DsmcCloud_H
+
+#include "Cloud.H"
+#include "DsmcBaseCloud.H"
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "Random.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "scalarIOField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+template<class CloudType>
+class BinaryCollisionModel;
+
+template<class CloudType>
+class WallInteractionModel;
+
+template<class CloudType>
+class InflowBoundaryModel;
+
+/*---------------------------------------------------------------------------*\
+                         Class DsmcCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class DsmcCloud
+:
+    public Cloud<ParcelType>,
+    public DsmcBaseCloud
+{
+    // Private data
+
+        //- Cloud type - used to set the name of the parcel properties
+        //  dictionary by appending "Properties"
+        const word cloudType_;
+
+        //- References to the mesh and time databases
+        const fvMesh& mesh_;
+
+        //- Dictionary of particle properties
+        IOdictionary particleProperties_;
+
+        //- A list of unique instances of molecule types in the simulation.
+        // The position of an entry in the list maps to the label identifying
+        // the typeId, i.e. where typeIdList_ = (N2 O2 CO2)
+        // N2 has typeId label = 0, O2 = 1, CO2 = 2.
+        List<word> typeIdList_;
+
+        //- Number of real atoms/molecules represented by a parcel
+        scalar nParticle_;
+
+        //- A data structure holding which particles are in which cell
+        List<DynamicList<ParcelType*> > cellOccupancy_;
+
+        //- An IOField holding the value of (sigmaT * cR)max for each cell (see
+        // Bird p220). Initialised with the parcels, updated as required, and
+        // read in on start/restart.
+        volScalarField sigmaTcRMax_;
+
+        //- A field holding the remainder from the previous collision selections
+        scalarField collisionSelectionRemainder_;
+
+        //- Heat flux at surface field
+        volScalarField q_;
+
+        //- Force density at surface field
+        volVectorField fD_;
+
+        //- Parcel constant properties - one for each type
+        List<typename ParcelType::constantProperties> constProps_;
+
+        //- Random number generator
+        Random rndGen_;
+
+
+        // References to the macroscopic fields
+
+            //- Temperature
+            const volScalarField& T_;
+
+            //- Velocity
+            const volVectorField& U_;
+
+
+        // References to the cloud sub-models
+
+            //- Binary collision model
+            autoPtr<BinaryCollisionModel<DsmcCloud<ParcelType> > >
+                binaryCollisionModel_;
+
+            //- Wall interaction model
+            autoPtr<WallInteractionModel<DsmcCloud<ParcelType> > >
+                wallInteractionModel_;
+
+            //- Inflow boundary model
+            autoPtr<InflowBoundaryModel<DsmcCloud<ParcelType> > >
+                inflowBoundaryModel_;
+
+
+    // Private Member Functions
+
+        //- Build the constant properties for all of the species
+        void buildConstProps();
+
+        //- Record which particles are in which cell
+        void buildCellOccupancy();
+
+        //- Initialise the system
+        void initialise(const IOdictionary& dsmcInitialiseDict);
+
+        //- Calculate collisions between molecules
+        void collisions();
+
+        //- Reset the surface data accumulation field values
+        void resetSurfaceDataFields();
+
+        //- Disallow default bitwise copy construct
+        DsmcCloud(const DsmcCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const DsmcCloud&);
+
+
+public:
+
+    // Static data members
+
+        //- Boltzmann constant
+        static scalar kb;
+
+        //- Reference temperature for all models
+        static scalar Tref;
+
+
+    // Constructors
+
+        //- Construct given name and mesh, will read Parcels from file
+        DsmcCloud
+        (
+            const word& cloudType,
+            const volScalarField& T,
+            const volVectorField& U
+        );
+
+        //- Construct given name and mesh.  Used to initialise.
+        DsmcCloud
+        (
+            const word& cloudType,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~DsmcCloud();
+
+
+    // Member Functions
+
+        // Access
+
+            // References to the mesh and databases
+
+                //- Return the cloud type
+                inline const word& cloudType() const;
+
+                //- Return refernce to the mesh
+                inline const fvMesh& mesh() const;
+
+
+            // References to the dsmc specific data
+
+                //- Return particle properties dictionary
+                inline const IOdictionary& particleProperties() const;
+
+                //- Return the idList
+                inline const List<word>& typeIdList() const;
+
+                //- Return the number of real particles represented by one
+                //  parcel
+                inline scalar nParticle() const;
+
+                //- Return the cell occupancy addressing
+                inline const List<DynamicList<ParcelType*> >&
+                    cellOccupancy() const;
+
+                //- Return the sigmaTcRMax field.  non-const access to allow
+                // updating.
+                inline volScalarField& sigmaTcRMax();
+
+                //- Return the collision selection remainder field.  non-const
+                // access to allow updating.
+                inline scalarField& collisionSelectionRemainder();
+
+                //- Return all of the constant properties
+                inline const List<typename ParcelType::constantProperties>&
+                    constProps() const;
+
+                //- Return the constant properties of the given typeId
+                inline const typename ParcelType::constantProperties&
+                    constProps(label typeId) const;
+
+                //- Return refernce to the random object
+                inline Random& rndGen();
+
+
+            // References to the surface data collection fields
+
+                //- Return heat flux at surface field
+                inline const volScalarField& q() const;
+
+                //- Return non-const heat flux at surface field
+                inline volScalarField& q();
+
+                //- Return force density at surface field
+                inline const volVectorField& fD() const;
+
+                //- Return non-const force density at surface field
+                inline volVectorField& fD();
+
+
+            // References to the macroscopic fields
+
+                //- Return macroscopic temperature
+                inline const volScalarField& T() const;
+
+                //- Return macroscopic velocity
+                inline const volVectorField& U() const;
+
+
+            // Kinetic theory helper functions
+
+                //- Generate a random velocity sampled from the Maxwellian speed
+                // distribution
+                vector equipartitionLinearVelocity
+                (
+                    scalar temperature,
+                    scalar mass
+                );
+
+                //- Generate a random internal energy, sampled from the
+                // equilibrium distribution (Bird eqn 11.22 and 11.23 and
+                // adapting code from DSMC3.FOR)
+                scalar equipartitionInternalEnergy
+                (
+                    scalar temperature,
+                    scalar internalDegreesOfFreedom
+                );
+
+
+                // From the Maxwellian distribution:
+                //- Average particle speed
+                inline scalar maxwellianAverageSpeed
+                (
+                    scalar temperature,
+                    scalar mass
+                ) const;
+
+                //- RMS particle speed
+                inline scalar maxwellianRMSSpeed
+                (
+                    scalar temperature,
+                    scalar mass
+                ) const;
+
+                //- Most probable speed
+                inline scalar maxwellianMostProbableSpeed
+                (
+                    scalar temperature,
+                    scalar mass
+                ) const;
+
+
+            // Sub-models
+
+                //- Return reference to binary elastic collision model
+                inline const BinaryCollisionModel<DsmcCloud<ParcelType> >&
+                    binaryCollision() const;
+
+                //- Return non-const reference to binary elastic collision model
+                inline BinaryCollisionModel<DsmcCloud<ParcelType> >&
+                    binaryCollision();
+
+                //- Return reference to wall interaction model
+                inline const WallInteractionModel<DsmcCloud<ParcelType> >&
+                    wallInteraction() const;
+
+                //- Return non-const reference to wall interaction model
+                inline WallInteractionModel<DsmcCloud<ParcelType> >&
+                    wallInteraction();
+
+                //- Return reference to wall interaction model
+                inline const InflowBoundaryModel<DsmcCloud<ParcelType> >&
+                    inflowBoundary() const;
+
+                //- Return non-const reference to wall interaction model
+                inline InflowBoundaryModel<DsmcCloud<ParcelType> >&
+                    inflowBoundary();
+
+
+        // Check
+
+            //- Total mass injected
+            inline scalar massInjected() const;
+
+            //- Total mass in system
+            inline scalar massInSystem() const;
+
+            //- Total linear momentum of the system
+            inline vector linearMomentumOfSystem() const;
+
+            //- Total linear kinetic energy in the system
+            inline scalar linearKineticEnergyOfSystem() const;
+
+            //- Total internal energy in the system
+            inline scalar internalEnergyOfSystem() const;
+
+            //- Print cloud information
+            void info() const;
+
+            //- Dump particle positions to .obj file
+            void dumpParticlePositions() const;
+
+
+            // Fields
+
+                //- Return the real particle number density field
+                inline const tmp<volScalarField> rhoN() const;
+
+                //- Return the particle mass density field
+                inline const tmp<volScalarField> rhoM() const;
+
+                //- Return the field of number of DSMC particles
+                inline const tmp<volScalarField> dsmcRhoN() const;
+
+                //- Return the momentum density field
+                inline const tmp<volVectorField> momentum() const;
+
+                //- Return the total linear kinetic energy (translational and
+                // thermal density field
+                inline const tmp<volScalarField> linearKE() const;
+
+                //- Return the internal energy density field
+                inline const tmp<volScalarField> internalE() const;
+
+                //- Return the average internal degrees of freedom  field
+                inline const tmp<volScalarField> iDof() const;
+
+
+        // Cloud evolution functions
+
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const vector& U,
+                const scalar Ei,
+                const label cellId,
+                const label typeId
+            );
+
+            //- Evolve the cloud (move, collide)
+            void evolve();
+
+
+            //- Clear the Cloud
+            inline void clear();
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "DsmcCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DsmcCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
new file mode 100644
index 0000000000..73899636e8
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -0,0 +1,598 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline const Foam::word& Foam::DsmcCloud<ParcelType>::cloudType() const
+{
+    return cloudType_;
+}
+
+
+template<class ParcelType>
+inline const Foam::fvMesh& Foam::DsmcCloud<ParcelType>::mesh() const
+{
+    return mesh_;
+}
+
+
+template<class ParcelType>
+inline const Foam::IOdictionary&
+Foam::DsmcCloud<ParcelType>::particleProperties() const
+{
+    return particleProperties_;
+}
+
+
+template<class ParcelType>
+inline const Foam::List<Foam::word>&
+Foam::DsmcCloud<ParcelType>::typeIdList() const
+{
+    return typeIdList_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::nParticle() const
+{
+    return nParticle_;
+}
+
+
+template<class ParcelType>
+inline const Foam::List<Foam::DynamicList<ParcelType*> >&
+Foam::DsmcCloud<ParcelType>::cellOccupancy() const
+{
+    return cellOccupancy_;
+}
+
+
+template<class ParcelType>
+inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::sigmaTcRMax()
+{
+    return sigmaTcRMax_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalarField&
+Foam::DsmcCloud<ParcelType>::collisionSelectionRemainder()
+{
+    return collisionSelectionRemainder_;
+}
+
+
+template<class ParcelType>
+inline const Foam::List<typename ParcelType::constantProperties>&
+Foam::DsmcCloud<ParcelType>::constProps() const
+{
+    return constProps_;
+}
+
+
+template<class ParcelType>
+inline const typename ParcelType::constantProperties&
+Foam::DsmcCloud<ParcelType>::constProps
+(
+    label typeId
+) const
+{
+    if (typeId < 0 || typeId >= constProps_.size())
+    {
+        FatalErrorIn("Foam::DsmcCloud<ParcelType>::constProps(label typeId)")
+            << "constantProperties for requested typeId index "
+            << typeId << " do not exist" << nl
+            << abort(FatalError);
+    }
+
+    return constProps_[typeId];
+}
+
+
+template<class ParcelType>
+inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
+{
+    return rndGen_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
+{
+    return q_;
+}
+
+
+template<class ParcelType>
+inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q()
+{
+    return q_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
+{
+    return fD_;
+}
+
+
+template<class ParcelType>
+inline Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD()
+{
+    return fD_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::T() const
+{
+    return T_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::U() const
+{
+    return U_;
+}
+
+
+template<class ParcelType>
+inline const Foam::BinaryCollisionModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::binaryCollision() const
+{
+    return binaryCollisionModel_;
+}
+
+
+template<class ParcelType>
+inline Foam::BinaryCollisionModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::binaryCollision()
+{
+    return binaryCollisionModel_();
+}
+
+
+template<class ParcelType>
+inline const Foam::WallInteractionModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::wallInteraction() const
+{
+    return wallInteractionModel_;
+}
+
+
+template<class ParcelType>
+inline Foam::WallInteractionModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::wallInteraction()
+{
+    return wallInteractionModel_();
+}
+
+
+template<class ParcelType>
+inline const Foam::InflowBoundaryModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::inflowBoundary() const
+{
+    return inflowBoundaryModel_;
+}
+
+
+template<class ParcelType>
+inline Foam::InflowBoundaryModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::inflowBoundary()
+{
+    return inflowBoundaryModel_();
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::massInSystem() const
+{
+    scalar sysMass = 0.0;
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+
+        const typename ParcelType::constantProperties& cP = constProps
+        (
+            p.typeId()
+        );
+
+        sysMass += cP.mass();
+    }
+
+    return nParticle_*sysMass;
+}
+
+
+template<class ParcelType>
+inline Foam::vector Foam::DsmcCloud<ParcelType>::linearMomentumOfSystem() const
+{
+    vector linearMomentum(vector::zero);
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+
+        const typename ParcelType::constantProperties& cP = constProps
+        (
+            p.typeId()
+        );
+
+        linearMomentum += cP.mass()*p.U();
+    }
+
+    return nParticle_*linearMomentum;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::DsmcCloud<ParcelType>::linearKineticEnergyOfSystem() const
+{
+    scalar linearKineticEnergy = 0.0;
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+
+        const typename ParcelType::constantProperties& cP = constProps
+        (
+            p.typeId()
+        );
+
+        linearKineticEnergy += 0.5*cP.mass()*(p.U() & p.U());
+    }
+
+    return nParticle_*linearKineticEnergy;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::DsmcCloud<ParcelType>::internalEnergyOfSystem() const
+{
+    scalar internalEnergy = 0.0;
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+
+        internalEnergy += p.Ei();
+    }
+
+    return nParticle_*internalEnergy;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
+(
+    scalar temperature,
+    scalar mass
+) const
+{
+    return 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass));
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
+(
+    scalar temperature,
+    scalar mass
+) const
+{
+    return sqrt(3.0*kb*temperature/mass);
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
+(
+    scalar temperature,
+    scalar mass
+) const
+{
+    return sqrt(2.0*kb*temperature/mass);
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::rhoN() const
+{
+    tmp<volScalarField> trhoN
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "rhoN",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& rhoN = trhoN().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        rhoN[cellI]++;
+    }
+
+    rhoN *= nParticle_/mesh().cellVolumes();
+
+    return trhoN;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::rhoM() const
+{
+    tmp<volScalarField> trhoM
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "rhoM",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& rhoM = trhoM().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        rhoM[cellI] += constProps(p.typeId()).mass();
+    }
+
+    rhoM *= nParticle_/mesh().cellVolumes();
+
+    return trhoM;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::dsmcRhoN() const
+{
+    tmp<volScalarField> tdsmcRhoN
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "dsmcRhoN",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& dsmcRhoN = tdsmcRhoN().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        dsmcRhoN[cellI]++;
+    }
+
+    return tdsmcRhoN;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volVectorField>
+Foam::DsmcCloud<ParcelType>::momentum() const
+{
+    tmp<volVectorField> tmomentum
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                this->name() + "momentum",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedVector
+            (
+                "zero",
+                dimensionSet(1, -2, -1, 0, 0),
+                vector::zero
+            )
+        )
+    );
+
+    vectorField& momentum = tmomentum().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        momentum[cellI] += constProps(p.typeId()).mass()*p.U();
+    }
+
+    momentum *= nParticle_/mesh().cellVolumes();
+
+    return tmomentum;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::linearKE() const
+{
+    tmp<volScalarField> tlinearKE
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "linearKE",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& linearKE = tlinearKE().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
+    }
+
+    linearKE *= nParticle_/mesh().cellVolumes();
+
+    return tlinearKE;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::internalE() const
+{
+    tmp<volScalarField> tinternalE
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "internalE",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& internalE = tinternalE().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        internalE[cellI] += p.Ei();
+    }
+
+    internalE *= nParticle_/mesh().cellVolumes();
+
+    return tinternalE;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::iDof() const
+{
+    tmp<volScalarField> tiDof
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "iDof",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& iDof = tiDof().internalField();
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
+    }
+
+    iDof *= nParticle_/mesh().cellVolumes();
+
+    return tiDof;
+}
+
+
+template<class ParcelType>
+inline void Foam::DsmcCloud<ParcelType>::clear()
+{
+    return IDLList<ParcelType>::clear();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C b/src/lagrangian/dsmc/clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C
new file mode 100644
index 0000000000..1e0c395ca1
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "DsmcBaseCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(DsmcBaseCloud, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::DsmcBaseCloud::DsmcBaseCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
+
+Foam::DsmcBaseCloud::~DsmcBaseCloud()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.H b/src/lagrangian/dsmc/clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.H
new file mode 100644
index 0000000000..df49f8e148
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.H
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::DsmcBaseCloud
+
+Description
+    Virtual abstract base class for templated DsmcCloud
+
+SourceFiles
+    DsmcBaseCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DsmcBaseCloud_H
+#define DsmcBaseCloud_H
+
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class DsmcBaseCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class DsmcBaseCloud
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        DsmcBaseCloud(const DsmcBaseCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const DsmcBaseCloud&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DsmcBaseCloud");
+
+    // Constructors
+
+        //- Null constructor
+        DsmcBaseCloud();
+
+    //- Destructor
+    virtual ~DsmcBaseCloud();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
new file mode 100644
index 0000000000..000ffcd53e
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dsmcCloud, 0);
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dsmcCloud::dsmcCloud
+(
+    const word& cloudType,
+    const volScalarField& T,
+    const volVectorField& U
+)
+:
+    DsmcCloud<dsmcParcel>(cloudType, T, U)
+{
+    dsmcParcel::readFields(*this);
+}
+
+
+Foam::dsmcCloud::dsmcCloud
+(
+    const word& cloudType,
+    const fvMesh& mesh
+)
+:
+    DsmcCloud<dsmcParcel>(cloudType, mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dsmcCloud::~dsmcCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::dsmcCloud::writeFields() const
+{
+    dsmcParcel::writeFields(*this);
+}
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
new file mode 100644
index 0000000000..2f88bbaad5
--- /dev/null
+++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::dsmcCloud
+
+Description
+    Cloud class to simulate dsmc parcels
+
+SourceFiles
+    dsmcCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dsmcCloud_H
+#define dsmcCloud_H
+
+#include "DsmcCloud.H"
+#include "dsmcParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class dsmcCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class dsmcCloud
+:
+    public DsmcCloud<dsmcParcel>
+{
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        dsmcCloud(const dsmcCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dsmcCloud&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dsmcCloud");
+
+
+    // Constructors
+
+        //- Construct from components
+        dsmcCloud
+        (
+            const word& cloudType,
+            const volScalarField& T,
+            const volVectorField& U
+        );
+
+        //- Construct from name and mesh, used to initialise.
+        dsmcCloud
+        (
+            const word& cloudType,
+            const fvMesh& mesh
+        );
+
+    //- Destructor
+    ~dsmcCloud();
+
+
+    // Member functions
+
+        //- Write fields
+        void writeFields() const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
new file mode 100644
index 0000000000..f9c1357189
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "DsmcParcel.H"
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackData>
+bool Foam::DsmcParcel<ParcelType>::move
+(
+    TrackData& td
+)
+{
+    ParcelType& p = static_cast<ParcelType&>(*this);
+
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    const polyMesh& mesh = td.cloud().pMesh();
+    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
+
+    const scalar deltaT = mesh.time().deltaT().value();
+    scalar tEnd = (1.0 - p.stepFraction())*deltaT;
+    const scalar dtMax = tEnd;
+
+    while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
+    {
+        // Set the Lagrangian time-step
+        scalar dt = min(dtMax, tEnd);
+
+        dt *= p.trackToFace(p.position() + dt*U_, td);
+
+        tEnd -= dt;
+
+        p.stepFraction() = 1.0 - tEnd/deltaT;
+
+        if (p.onBoundary() && td.keepParticle)
+        {
+            if (p.face() > -1)
+            {
+                if
+                (
+                    isType<processorPolyPatch>(pbMesh[p.patch(p.face())])
+                )
+                {
+                    td.switchProcessor = true;
+                }
+            }
+        }
+    }
+
+    return td.keepParticle;
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::DsmcParcel<ParcelType>::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    TrackData& td
+)
+{
+    td.switchProcessor = true;
+}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    int&
+)
+{}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::DsmcParcel<ParcelType>::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    TrackData& td
+)
+{
+    const constantProperties& constProps(td.cloud().constProps(typeId_));
+
+    scalar m = constProps.mass();
+
+    // pre-interaction energy
+    scalar preIE = 0.5*m*(U_ & U_) + Ei_;
+
+    // pre-interaction momentum
+    vector preIMom = m*U_;
+
+    td.cloud().wallInteraction().correct
+    (
+        wpp,
+        this->face(),
+        U_,
+        Ei_,
+        typeId_
+    );
+
+    // post-interaction energy
+    scalar postIE = 0.5*m*(U_ & U_) + Ei_;
+
+    // post-interaction momentum
+    vector postIMom = m*U_;
+
+    label wppIndex = wpp.index();
+
+    label wppLocalFace = wpp.whichFace(this->face());
+
+    const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
+
+    const scalar deltaT = td.cloud().mesh().time().deltaT().value();
+
+    scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);
+
+    vector deltaFD = td.cloud().nParticle()*(preIMom - postIMom)/(deltaT*fA);
+
+    td.cloud().q().boundaryField()[wppIndex][wppLocalFace] += deltaQ;
+
+    td.cloud().fD().boundaryField()[wppIndex][wppLocalFace] += deltaFD;
+}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::hitWallPatch
+(
+    const wallPolyPatch&,
+    int&
+)
+{}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::DsmcParcel<ParcelType>::hitPatch
+(
+    const polyPatch&,
+    TrackData& td
+)
+{
+    td.keepParticle = false;
+}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::hitPatch
+(
+    const polyPatch&,
+    int&
+)
+{}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::transformProperties
+(
+    const tensor& T
+)
+{
+    Particle<ParcelType>::transformProperties(T);
+    U_ = transform(T, U_);
+}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::transformProperties
+(
+    const vector& separation
+)
+{
+    Particle<ParcelType>::transformProperties(separation);
+}
+
+
+// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
+
+#include "DsmcParcelIO.C"
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
new file mode 100644
index 0000000000..bbc2fd1539
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
@@ -0,0 +1,347 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::DsmcParcel
+
+Description
+    DSMC parcel class
+
+SourceFiles
+    DsmcParcelI.H
+    DsmcParcel.C
+    DsmcParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DsmcParcel_H
+#define DsmcParcel_H
+
+#include "Particle.H"
+#include "IOstream.H"
+#include "autoPtr.H"
+#include "contiguous.H"
+#include "mathematicalConstants.H"
+
+#include "DsmcCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class ParcelType>
+class DsmcParcel;
+
+// Forward declaration of friend functions
+
+template<class ParcelType>
+Ostream& operator<<
+(
+    Ostream&,
+    const DsmcParcel<ParcelType>&
+);
+
+/*---------------------------------------------------------------------------*\
+                         Class DsmcParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class DsmcParcel
+:
+    public Particle<ParcelType>
+{
+public:
+
+    //- Class to hold DSMC particle constant properties
+    class constantProperties
+    {
+        // Private data
+
+            //- Particle mass [kg] (constant)
+            scalar mass_;
+
+            //- Particle hard sphere diameter [m] (constant)
+            scalar d_;
+
+            //- Internal degrees of freedom
+            scalar internalDegreesOfFreedom_;
+
+            //- Viscosity index
+            scalar omega_;
+
+
+    public:
+
+        // Constrcutors
+
+            //- Null constructor, allows List of constantProperties to be
+            //  created before the contents is initialised
+            inline constantProperties();
+
+            //- Constructor from dictionary
+            inline constantProperties(const dictionary& dict);
+
+
+        // Member functions
+
+            //- Return const access to the particle density
+            inline scalar mass() const;
+
+            //- Return const access to the minimum particle mass
+            inline scalar d() const;
+
+            //- Return the reference total collision cross section
+            inline scalar sigmaT() const;
+
+            //- Return the internalDegreesOfFreedom
+            inline scalar internalDegreesOfFreedom() const;
+
+            //- Return the viscosity index
+            inline scalar omega() const;
+
+    };
+
+
+    //- Class used to pass kinematic tracking data to the trackToFace function
+    class trackData
+    :
+        public Particle<ParcelType>::trackData
+    {
+        // Private data
+
+            //- Reference to the cloud containing this particle
+            DsmcCloud<ParcelType>& cloud_;
+
+
+    public:
+
+        // Constructors
+
+            //- Construct from components
+            inline trackData
+            (
+                DsmcCloud<ParcelType>& cloud
+            );
+
+
+        // Member functions
+
+            //- Return access to the owner cloud
+            inline DsmcCloud<ParcelType>& cloud();
+    };
+
+
+protected:
+
+    // Protected member data
+
+        // Parcel properties
+
+            //- Velocity of Parcel [m/s]
+            vector U_;
+
+            //- Internal energy of the Parcel, covering all non-translational
+            // degrees of freedom [J]
+            scalar Ei_;
+
+            //- Parcel type id
+            label typeId_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DsmcParcel");
+
+    friend class Cloud<ParcelType>;
+
+
+    // Constructors
+
+        //- Construct from components
+        inline DsmcParcel
+        (
+            DsmcCloud<ParcelType>& owner,
+            const vector& position,
+            const vector& U,
+            const scalar Ei,
+            const label celli,
+            const label typeId
+        );
+
+        //- Construct from Istream
+        DsmcParcel
+        (
+            const Cloud<ParcelType>& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<ParcelType> clone() const
+        {
+            return autoPtr<ParcelType>(new DsmcParcel<ParcelType>(*this));
+        }
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return type id
+            inline label typeId() const;
+
+            //- Return const access to velocity
+            inline const vector& U() const;
+
+            //- Return const access to internal energy
+            inline scalar Ei() const;
+
+        // Edit
+
+            //- Return access to velocity
+            inline vector& U();
+
+            //- Return access to internal energy
+            inline scalar& Ei();
+
+
+        // Main calculation loop
+
+            // Tracking
+
+                //- Move the parcel
+                template<class TrackData>
+                bool move(TrackData& td);
+
+
+        // Patch interactions
+
+            //- Overridable function to handle the particle hitting a
+            //  processorPatch
+            template<class TrackData>
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                TrackData& td
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //  processorPatch without trackData
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a wallPatch
+            template<class TrackData>
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                TrackData& td
+            );
+
+            //- Overridable function to handle the particle hitting a wallPatch
+            //  without trackData
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a polyPatch
+            template<class TrackData>
+            void hitPatch
+            (
+                const polyPatch&,
+                TrackData& td
+            );
+
+            //- Overridable function to handle the particle hitting a polyPatch
+            //- without trackData
+            void hitPatch
+            (
+                const polyPatch&,
+                int&
+            );
+
+            //- Transform the physical properties of the particle
+            //  according to the given transformation tensor
+            void transformProperties(const tensor& T);
+
+            //- Transform the physical properties of the particle
+            //  according to the given separation vector
+            void transformProperties(const vector& separation);
+
+
+        // I-O
+
+            static void readFields(DsmcCloud<ParcelType>& c);
+
+            static void writeFields(const DsmcCloud<ParcelType>& c);
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParcelType>
+        (
+            Ostream&,
+            const DsmcParcel<ParcelType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "DsmcParcelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define defineParcelTypeNameAndDebug(Type, DebugSwitch)                       \
+    template<>                                                                \
+    const Foam::word DsmcParcel<Type>::typeName(#Type);                       \
+    template<>                                                                \
+    int DsmcParcel<Type>::debug                                               \
+    (                                                                         \
+        Foam::debug::debugSwitch(#Type, DebugSwitch)                          \
+    );
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "DsmcParcel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
new file mode 100644
index 0000000000..ff86419186
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class ParcelType>
+inline Foam::DsmcParcel<ParcelType>::constantProperties::constantProperties()
+:
+    mass_(0),
+    d_(0)
+{}
+
+
+template <class ParcelType>
+inline Foam::DsmcParcel<ParcelType>::constantProperties::constantProperties
+(
+    const dictionary& dict
+)
+:
+    mass_(readScalar(dict.lookup("mass"))),
+    d_(readScalar(dict.lookup("diameter"))),
+    internalDegreesOfFreedom_
+    (
+        readScalar(dict.lookup("internalDegreesOfFreedom"))
+    ),
+    omega_(readScalar(dict.lookup("omega")))
+{}
+
+
+template <class ParcelType>
+inline Foam::DsmcParcel<ParcelType>::trackData::trackData
+(
+    DsmcCloud<ParcelType>& cloud
+)
+:
+    Particle<ParcelType>::trackData(cloud),
+    cloud_(cloud)
+{}
+
+
+template <class ParcelType>
+inline Foam::DsmcParcel<ParcelType>::DsmcParcel
+(
+    DsmcCloud<ParcelType>& owner,
+    const vector& position,
+    const vector& U,
+    const scalar Ei,
+    const label celli,
+    const label typeId
+)
+:
+    Particle<ParcelType>(owner, position, celli),
+    U_(U),
+    Ei_(Ei),
+    typeId_(typeId)
+{}
+
+
+// * * * * * * * * * constantProperties Member Functions * * * * * * * * * * //
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::mass() const
+{
+    return mass_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::d() const
+{
+    return d_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::sigmaT() const
+{
+    return mathematicalConstant::pi*d_*d_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::internalDegreesOfFreedom()
+const
+{
+    return internalDegreesOfFreedom_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::omega() const
+{
+    return omega_;
+}
+
+
+// * * * * * * * * * * * trackData Member Functions  * * * * * * * * * * * * //
+
+template <class ParcelType>
+inline Foam::DsmcCloud<ParcelType>&
+Foam::DsmcParcel<ParcelType>::trackData::cloud()
+{
+    return cloud_;
+}
+
+
+// * * * * * * * * * * DsmcParcel Member Functions  * * * * * * * * * * //
+
+template <class ParcelType>
+inline Foam::label Foam::DsmcParcel<ParcelType>::typeId() const
+{
+    return typeId_;
+}
+
+
+template <class ParcelType>
+inline const Foam::vector& Foam::DsmcParcel<ParcelType>::U() const
+{
+    return U_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar Foam::DsmcParcel<ParcelType>::Ei() const
+{
+    return Ei_;
+}
+
+
+template <class ParcelType>
+inline Foam::vector& Foam::DsmcParcel<ParcelType>::U()
+{
+    return U_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar& Foam::DsmcParcel<ParcelType>::Ei()
+{
+    return Ei_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
new file mode 100644
index 0000000000..b8bc688c09
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "DsmcParcel.H"
+#include "IOstreams.H"
+#include "IOField.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class ParcelType>
+Foam::DsmcParcel<ParcelType>::DsmcParcel
+(
+    const Cloud<ParcelType>& cloud,
+    Istream& is,
+    bool readFields
+)
+:
+    Particle<ParcelType>(cloud, is, readFields),
+    U_(vector::zero),
+    Ei_(0.0),
+    typeId_(-1)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            is >> U_;
+            Ei_ = readScalar(is);
+            typeId_ = readLabel(is);
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&U_),
+                sizeof(U_)
+              + sizeof(Ei_)
+              + sizeof(typeId_)
+            );
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "DsmcParcel<ParcelType>::DsmcParcel"
+        "(const Cloud<ParcelType>&, Istream&, bool)"
+    );
+}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::readFields
+(
+    DsmcCloud<ParcelType>& c
+)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    IOField<vector> U(c.fieldIOobject("U", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, U);
+
+    IOField<scalar> Ei(c.fieldIOobject("Ei", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, Ei);
+
+    IOField<label> typeId(c.fieldIOobject("typeId", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, typeId);
+
+    label i = 0;
+    forAllIter(typename Cloud<ParcelType>, c, iter)
+    {
+        ParcelType& p = iter();
+
+        p.U_ = U[i];
+        p.Ei_ = Ei[i];
+        p.typeId_ = typeId[i];
+        i++;
+    }
+}
+
+
+template<class ParcelType>
+void Foam::DsmcParcel<ParcelType>::writeFields
+(
+    const DsmcCloud<ParcelType>& c
+)
+{
+    Particle<ParcelType>::writeFields(c);
+
+    label np =  c.size();
+
+    IOField<vector> U(c.fieldIOobject("U", IOobject::NO_READ), np);
+    IOField<scalar> Ei(c.fieldIOobject("Ei", IOobject::NO_READ), np);
+    IOField<label> typeId(c.fieldIOobject("typeId", IOobject::NO_READ), np);
+
+    label i = 0;
+    forAllConstIter(typename Cloud<ParcelType>, c, iter)
+    {
+        const DsmcParcel<ParcelType>& p = iter();
+
+        U[i] = p.U();
+        Ei[i] = p.Ei();
+        typeId[i] = p.typeId();
+        i++;
+    }
+
+    U.write();
+    Ei.write();
+    typeId.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const DsmcParcel<ParcelType>& p
+)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const Particle<ParcelType>& >(p)
+            << token::SPACE << p.U()
+            << token::SPACE << p.Ei()
+            << token::SPACE << p.typeId();
+    }
+    else
+    {
+        os  << static_cast<const Particle<ParcelType>& >(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.U_),
+            sizeof(p.U())
+          + sizeof(p.Ei())
+          + sizeof(p.typeId())
+        );
+    }
+
+    // Check state of Ostream
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, const DsmcParcel<ParcelType>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/defineDsmcParcel.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/defineDsmcParcel.C
new file mode 100644
index 0000000000..d3b71359da
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/defineDsmcParcel.C
@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcParcel.H"
+#include "DsmcCloud.H"
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(Cloud<dsmcParcel>, 0);
+
+    defineParcelTypeNameAndDebug(DsmcCloud<dsmcParcel>, 0);
+};
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C
new file mode 100644
index 0000000000..8f52fcdabf
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C
@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcParcel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dsmcParcel, 0);
+    defineParticleTypeNameAndDebug(dsmcParcel, 0);
+    defineParcelTypeNameAndDebug(dsmcParcel, 0);
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dsmcParcel::dsmcParcel
+(
+    DsmcCloud<dsmcParcel>& owner,
+    const vector& position,
+    const vector& U,
+    const scalar Ei,
+    const label celli,
+    const label typeId
+)
+:
+    DsmcParcel<dsmcParcel>
+    (
+        owner,
+        position,
+        U,
+        Ei,
+        celli,
+        typeId
+    )
+{}
+
+
+Foam::dsmcParcel::dsmcParcel
+(
+    const Cloud<dsmcParcel>& cloud,
+    Istream& is,
+    bool readFields
+)
+:
+    DsmcParcel<dsmcParcel>(cloud, is, readFields)
+{}
+
+
+// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
+
+Foam::dsmcParcel::~dsmcParcel()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H
new file mode 100644
index 0000000000..24d9c9ada4
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::dsmcParcel
+
+Description
+
+
+SourceFiles
+    dsmcParcel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dsmcParcel_H
+#define dsmcParcel_H
+
+#include "DsmcParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class dsmcParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+class dsmcParcel
+:
+    public DsmcParcel<dsmcParcel>
+{
+
+public:
+
+    //- Run-time type information
+    TypeName("dsmcParcel");
+
+
+    // Constructors
+
+        //- Construct from components
+        dsmcParcel
+        (
+            DsmcCloud<dsmcParcel>& owner,
+            const vector& position,
+            const vector& U,
+            const scalar Ei,
+            const label celli,
+            const label typeId
+        );
+
+        //- Construct from Istream
+        dsmcParcel
+        (
+            const Cloud<dsmcParcel>& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<dsmcParcel> clone() const
+        {
+            return autoPtr<dsmcParcel>(new dsmcParcel(*this));
+        }
+
+
+    //- Destructor
+    virtual ~dsmcParcel();
+};
+
+
+template<>
+inline bool contiguous<dsmcParcel>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C
new file mode 100644
index 0000000000..de2444687b
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcParcel.H"
+#include "DsmcCloud.H"
+#include "VariableHardSphere.H"
+#include "LarsenBorgnakkeVariableHardSphere.H"
+
+namespace Foam
+{
+    makeBinaryCollisionModel(DsmcCloud<dsmcParcel>);
+
+    // Add instances of collision model to the table
+    makeBinaryCollisionModelType
+    (
+        VariableHardSphere,
+        DsmcCloud,
+        dsmcParcel
+    );
+    makeBinaryCollisionModelType
+    (
+        LarsenBorgnakkeVariableHardSphere,
+        DsmcCloud,
+        dsmcParcel
+    );
+};
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C
new file mode 100644
index 0000000000..c3ee121ff0
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcParcel.H"
+#include "DsmcCloud.H"
+#include "FreeStream.H"
+#include "NoInflow.H"
+
+namespace Foam
+{
+    makeInflowBoundaryModel(DsmcCloud<dsmcParcel>);
+
+    // Add instances of inflow boundary model to the table
+    makeInflowBoundaryModelType
+    (
+        FreeStream,
+        DsmcCloud,
+        dsmcParcel
+    );
+    makeInflowBoundaryModelType
+    (
+        NoInflow,
+        DsmcCloud,
+        dsmcParcel
+    );
+};
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
new file mode 100644
index 0000000000..874c787d6a
--- /dev/null
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcParcel.H"
+#include "DsmcCloud.H"
+#include "MaxwellianThermal.H"
+#include "SpecularReflection.H"
+
+namespace Foam
+{
+    makeWallInteractionModel(DsmcCloud<dsmcParcel>);
+
+    // Add instances of wall interaction model to the table
+    makeWallInteractionModelType
+    (
+        MaxwellianThermal,
+        DsmcCloud,
+        dsmcParcel
+    );
+    makeWallInteractionModelType
+    (
+        SpecularReflection,
+        DsmcCloud,
+        dsmcParcel
+    );
+};
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C
new file mode 100644
index 0000000000..0118014f9c
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "BinaryCollisionModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::BinaryCollisionModel<CloudType>::BinaryCollisionModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::BinaryCollisionModel<CloudType>::~BinaryCollisionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+const CloudType&
+Foam::BinaryCollisionModel<CloudType>::owner() const
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+CloudType&
+Foam::BinaryCollisionModel<CloudType>::owner()
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary&
+Foam::BinaryCollisionModel<CloudType>::dict() const
+{
+    return dict_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary&
+Foam::BinaryCollisionModel<CloudType>::coeffDict() const
+{
+    return coeffDict_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NewBinaryCollisionModel.C"
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
new file mode 100644
index 0000000000..771115c984
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::BinaryCollisionModel
+
+Description
+    Templated DSMC particle collision class
+
+SourceFiles
+    BinaryCollisionModel.C
+    NewBinaryCollisionModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BinaryCollisionModel_H
+#define BinaryCollisionModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class BinaryCollisionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class BinaryCollisionModel
+{
+    // Private data
+
+        //- The cloud dictionary
+        const dictionary& dict_;
+
+        // reference to the owner cloud class
+        CloudType& owner_;
+
+        //- The coefficients dictionary
+        const dictionary coeffDict_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("BinaryCollisionModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        BinaryCollisionModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct from components
+        BinaryCollisionModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+
+    // Destructor
+    virtual ~BinaryCollisionModel();
+
+
+    //- Selector
+    static autoPtr<BinaryCollisionModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Access
+
+        //- Return the owner cloud object
+        const CloudType& owner() const;
+
+        //- Return non-const access to the owner cloud object
+        CloudType& owner();
+
+        //- Return the dictionary
+        const dictionary& dict() const;
+
+        //- Return the coefficients dictionary
+        const dictionary& coeffDict() const;
+
+
+    // Member Functions
+
+        //- Return the collision cross section * relative velocity product
+        virtual scalar sigmaTcR
+        (
+            label typeIdP,
+            label typeIdQ,
+            const vector& UP,
+            const vector& UQ
+        ) const = 0;
+
+        //- Apply collision
+        virtual void collide
+        (
+            label typeIdP,
+            label typeIdQ,
+            vector& UP,
+            vector& UQ,
+            scalar& EiP,
+            scalar& EiQ
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeBinaryCollisionModel(CloudType)                            \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        BinaryCollisionModel<CloudType>,                               \
+        0                                                                     \
+    );                                                                        \
+                                                                              \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        BinaryCollisionModel<CloudType>,                               \
+        dictionary                                                            \
+    );
+
+
+#define makeBinaryCollisionModelType(SS, CloudType, ParcelType)        \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(SS<CloudType<ParcelType> >, 0);       \
+                                                                              \
+    BinaryCollisionModel<CloudType<ParcelType> >::                     \
+        adddictionaryConstructorToTable<SS<CloudType<ParcelType> > >          \
+            add##SS##CloudType##ParcelType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BinaryCollisionModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/NewBinaryCollisionModel.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/NewBinaryCollisionModel.C
new file mode 100644
index 0000000000..8878164cdb
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/NewBinaryCollisionModel.C
@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "BinaryCollisionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::BinaryCollisionModel<CloudType> >
+Foam::BinaryCollisionModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word BinaryCollisionModelType(dict.lookup("BinaryCollisionModel"));
+
+    Info<< "Selecting BinaryCollisionModel "
+        << BinaryCollisionModelType
+        << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(BinaryCollisionModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "BinaryCollisionModel<CloudType>::New"
+            "(const dictionary&, CloudType&)"
+        )
+            << "Unknown BinaryCollisionModelType type "
+            << BinaryCollisionModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid BinaryCollisionModel types are:" << nl
+            << dictionaryConstructorTablePtr_->toc() << exit(FatalError);
+    }
+
+    return autoPtr<BinaryCollisionModel<CloudType> >
+    (
+        cstrIter()(dict, owner)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
new file mode 100644
index 0000000000..b6deaf1f00
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
@@ -0,0 +1,271 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "LarsenBorgnakkeVariableHardSphere.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
+(
+    scalar ChiA,
+    scalar ChiB
+)
+{
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    scalar ChiAMinusOne = ChiA - 1;
+
+    scalar ChiBMinusOne = ChiB - 1;
+
+    if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
+    {
+        return rndGen.scalar01();
+    }
+
+    scalar energyRatio;
+
+    scalar P;
+
+    do
+    {
+        P = 0;
+
+        energyRatio = rndGen.scalar01();
+
+        if (ChiAMinusOne < SMALL)
+        {
+            P = 1.0 - pow(energyRatio, ChiB);
+        }
+        else if (ChiBMinusOne < SMALL)
+        {
+            P = 1.0 - pow(energyRatio, ChiA);
+        }
+        else
+        {
+            P =
+                pow
+                (
+                    (ChiAMinusOne + ChiBMinusOne)*energyRatio/ChiAMinusOne,
+                    ChiAMinusOne
+                )
+               *pow
+                (
+                    (ChiAMinusOne + ChiBMinusOne)*(1 - energyRatio)
+                    /ChiBMinusOne,
+                    ChiBMinusOne
+                );
+        }
+    } while (P < rndGen.scalar01());
+
+    return energyRatio;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::LarsenBorgnakkeVariableHardSphere
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    BinaryCollisionModel<CloudType>(dict, cloud, typeName),
+    relaxationCollisionNumber_
+    (
+        readScalar(this->coeffDict().lookup("relaxationCollisionNumber"))
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::
+~LarsenBorgnakkeVariableHardSphere()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+template <class CloudType>
+Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::sigmaTcR
+(
+    label typeIdP,
+    label typeIdQ,
+    const vector& UP,
+    const vector& UQ
+) const
+{
+    const CloudType& cloud(this->owner());
+
+    scalar dPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).d()
+          + cloud.constProps(typeIdQ).d()
+        );
+
+    scalar omegaPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).omega()
+          + cloud.constProps(typeIdQ).omega()
+        );
+
+    scalar cR = mag(UP - UQ);
+
+    if (cR < VSMALL)
+    {
+        return 0;
+    }
+
+    scalar mP = cloud.constProps(typeIdP).mass();
+
+    scalar mQ = cloud.constProps(typeIdQ).mass();
+
+    scalar mR = mP*mQ/(mP + mQ);
+
+    // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
+    scalar sigmaTPQ =
+        mathematicalConstant::pi*dPQ*dPQ
+       *pow(2.0*CloudType::kb*CloudType::Tref/(mR*cR*cR), omegaPQ - 0.5)
+       /exp(Foam::lgamma(2.5 - omegaPQ));
+
+    return sigmaTPQ*cR;
+}
+
+
+template <class CloudType>
+void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
+(
+    label typeIdP,
+    label typeIdQ,
+    vector& UP,
+    vector& UQ,
+    scalar& EiP,
+    scalar& EiQ
+)
+{
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
+
+    // Larsen Borgnakke internal energy redistribution part.  Using the serial
+    // application of the LB method, as per the INELRS subroutine in Bird's
+    // DSMC0R.FOR
+
+    scalar preCollisionEiP = EiP;
+
+    scalar preCollisionEiQ = EiQ;
+
+    scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
+
+    scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
+
+    scalar omegaPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).omega()
+          + cloud.constProps(typeIdQ).omega()
+        );
+
+    scalar mP = cloud.constProps(typeIdP).mass();
+
+    scalar mQ = cloud.constProps(typeIdQ).mass();
+
+    scalar mR = mP*mQ/(mP + mQ);
+
+    vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
+
+    scalar cRsqr = magSqr(UP - UQ);
+
+    scalar availableEnergy = 0.5*mR*cRsqr;
+
+    scalar ChiB = 2.5 - omegaPQ;
+
+    if (iDofP > 0)
+    {
+        if (inverseCollisionNumber > rndGen.scalar01())
+        {
+            availableEnergy += preCollisionEiP;
+
+            scalar ChiA = 0.5*iDofP;
+
+            EiP = energyRatio(ChiA, ChiB)*availableEnergy;
+
+            availableEnergy -= EiP;
+        }
+    }
+
+    if (iDofQ > 0)
+    {
+        if (inverseCollisionNumber > rndGen.scalar01())
+        {
+            availableEnergy += preCollisionEiQ;
+
+            // Change to general LB ratio calculation
+            scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
+
+            EiQ = energyRatio*availableEnergy;
+
+            availableEnergy -= EiQ;
+        }
+    }
+
+    // Rescale the translational energy
+    scalar cR = sqrt(2.0*availableEnergy/mR);
+
+    // Variable Hard Sphere collision part
+
+    scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
+
+    scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
+
+    scalar phi = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+
+    vector postCollisionRelU =
+        cR
+       *vector
+        (
+            cosTheta,
+            sinTheta*cos(phi),
+            sinTheta*sin(phi)
+        );
+
+    UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
+
+    UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H
new file mode 100644
index 0000000000..c95d70f24c
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+Class
+    Foam::LarsenBorgnakkeVariableHardSphere
+
+Description
+    Variable Hard Sphere BinaryCollision Model with Larsen Borgnakke internal
+    energy redistribution.  Based on the INELRS subroutine in Bird's DSMC0R.FOR
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LarsenBorgnakkeVariableHardSphere_H
+#define LarsenBorgnakkeVariableHardSphere_H
+
+#include "BinaryCollisionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                      Class LarsenBorgnakkeVariableHardSphere Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class LarsenBorgnakkeVariableHardSphere
+:
+    public BinaryCollisionModel<CloudType>
+{
+    // Private data
+
+        //- Relaxation collision number
+        const scalar relaxationCollisionNumber_;
+
+
+    // Private Member Functions
+
+        //- Calculate the energy ratio for distribution to internal degrees of
+        // freedom
+        scalar energyRatio
+        (
+            scalar ChiA,
+            scalar ChiB
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LarsenBorgnakkeVariableHardSphere");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        LarsenBorgnakkeVariableHardSphere
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~LarsenBorgnakkeVariableHardSphere();
+
+
+    // Member Functions
+
+        //- Return the collision cross section * relative velocity product
+        virtual scalar sigmaTcR
+        (
+            label typeIdP,
+            label typeIdQ,
+            const vector& UP,
+            const vector& UQ
+        ) const;
+
+        //- Apply collision
+        virtual void collide
+        (
+            label typeIdP,
+            label typeIdQ,
+            vector& UP,
+            vector& UQ,
+            scalar& EiP,
+            scalar& EiQ
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LarsenBorgnakkeVariableHardSphere.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
new file mode 100644
index 0000000000..de47f1550d
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "VariableHardSphere.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::VariableHardSphere<CloudType>::VariableHardSphere
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    BinaryCollisionModel<CloudType>(dict, cloud, typeName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::VariableHardSphere<CloudType>::~VariableHardSphere()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+template <class CloudType>
+Foam::scalar Foam::VariableHardSphere<CloudType>::sigmaTcR
+(
+    label typeIdP,
+    label typeIdQ,
+    const vector& UP,
+    const vector& UQ
+) const
+{
+    const CloudType& cloud(this->owner());
+
+    scalar dPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).d()
+          + cloud.constProps(typeIdQ).d()
+        );
+
+    scalar omegaPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).omega()
+          + cloud.constProps(typeIdQ).omega()
+        );
+
+    scalar cR = mag(UP - UQ);
+
+    if (cR < VSMALL)
+    {
+        return 0;
+    }
+
+    scalar mP = cloud.constProps(typeIdP).mass();
+
+    scalar mQ = cloud.constProps(typeIdQ).mass();
+
+    scalar mR = mP*mQ/(mP + mQ);
+
+    // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
+    scalar sigmaTPQ =
+        mathematicalConstant::pi*dPQ*dPQ
+       *pow(2.0*CloudType::kb*CloudType::Tref/(mR*cR*cR), omegaPQ - 0.5)
+       /exp(Foam::lgamma(2.5 - omegaPQ));
+
+    return sigmaTPQ*cR;
+}
+
+
+template <class CloudType>
+void Foam::VariableHardSphere<CloudType>::collide
+(
+    label typeIdP,
+    label typeIdQ,
+    vector& UP,
+    vector& UQ,
+    scalar& EiP,
+    scalar& EiQ
+)
+{
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    scalar mP = cloud.constProps(typeIdP).mass();
+
+    scalar mQ = cloud.constProps(typeIdQ).mass();
+
+    vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
+
+    scalar cR = mag(UP - UQ);
+
+    scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
+
+    scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
+
+    scalar phi = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+
+    vector postCollisionRelU =
+        cR
+       *vector
+        (
+            cosTheta,
+            sinTheta*cos(phi),
+            sinTheta*sin(phi)
+        );
+
+    UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
+
+    UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
new file mode 100644
index 0000000000..a1dbdb9a27
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+Class
+    Foam::VariableHardSphere
+
+Description
+    Variable Hard Sphere BinaryCollision Model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VariableHardSphere_H
+#define VariableHardSphere_H
+
+#include "BinaryCollisionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                      Class VariableHardSphere Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class VariableHardSphere
+:
+    public BinaryCollisionModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("VariableHardSphere");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        VariableHardSphere
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~VariableHardSphere();
+
+
+    // Member Functions
+
+        //- Return the collision cross section * relative velocity product
+        virtual scalar sigmaTcR
+        (
+            label typeIdP,
+            label typeIdQ,
+            const vector& UP,
+            const vector& UQ
+        ) const;
+
+        //- Apply collision
+        virtual void collide
+        (
+            label typeIdP,
+            label typeIdQ,
+            vector& UP,
+            vector& UQ,
+            scalar& EiP,
+            scalar& EiQ
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "VariableHardSphere.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
new file mode 100644
index 0000000000..f4203f4d66
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -0,0 +1,280 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "FreeStream.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::FreeStream<CloudType>::FreeStream
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    InflowBoundaryModel<CloudType>(dict, cloud, typeName),
+    patchIndex_(),
+    temperature_(readScalar(this->coeffDict().lookup("temperature"))),
+    velocity_(this->coeffDict().lookup("velocity")),
+    moleculeTypeIds_(),
+    numberDensities_(),
+    particleFluxAccumulators_()
+{
+    word patchName = this->coeffDict().lookup("patch");
+
+    patchIndex_ = cloud.mesh().boundaryMesh().findPatchID(patchName);
+
+    const polyPatch& patch = cloud.mesh().boundaryMesh()[patchIndex_];
+
+    if (patchIndex_ == -1)
+    {
+        FatalErrorIn
+        (
+            "Foam::FreeStream<CloudType>::FreeStream"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "patch " << patchName << " not found." << nl
+            << abort(FatalError);
+    }
+
+    const dictionary& numberDensitiesDict
+    (
+        this->coeffDict().subDict("numberDensities")
+    );
+
+    List<word> molecules(numberDensitiesDict.toc());
+
+    numberDensities_.setSize(molecules.size());
+
+    moleculeTypeIds_.setSize(molecules.size());
+
+    forAll(molecules, i)
+    {
+        numberDensities_[i] = readScalar
+        (
+            numberDensitiesDict.lookup(molecules[i])
+        );
+
+        moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]);
+
+        if (moleculeTypeIds_[i] == -1)
+        {
+            FatalErrorIn
+            (
+                "Foam::FreeStream<CloudType>::FreeStream"
+                "("
+                    "const dictionary&, "
+                    "CloudType&"
+                ")"
+            )   << "typeId " << molecules[i] << "not defined in cloud." << nl
+                << abort(FatalError);
+        }
+    }
+
+    numberDensities_ /= cloud.nParticle();
+
+    particleFluxAccumulators_.setSize
+    (
+        molecules.size(),
+        Field<scalar>(patch.size(), 0)
+    );
+}
+
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::FreeStream<CloudType>::~FreeStream()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::FreeStream<CloudType>::inflow()
+{
+    CloudType& cloud(this->owner());
+
+    const polyMesh& mesh(cloud.mesh());
+
+    const scalar deltaT = mesh.time().deltaT().value();
+
+    Random& rndGen(cloud.rndGen());
+
+    const polyPatch& patch = mesh.boundaryMesh()[patchIndex_];
+
+    label particlesInserted = 0;
+
+    // Add mass to the accumulators.  negative face area dotted with the
+    // velocity to point flux into the domain.
+
+    forAll(particleFluxAccumulators_, i)
+    {
+        particleFluxAccumulators_[i] +=
+           -patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
+    }
+
+    forAll(patch, f)
+    {
+        // Loop over all faces as the outer loop to avoid calculating
+        // geometrical properties multiple times.
+
+        labelList faceVertices = patch[f];
+
+        label nVertices = faceVertices.size();
+
+        label globalFaceIndex = f + patch.start();
+
+        label cell = mesh.faceOwner()[globalFaceIndex];
+
+        const vector& fC = patch.faceCentres()[f];
+
+        scalar fA = mag(patch.faceAreas()[f]);
+
+        // Cummulative triangle area fractions
+        List<scalar> cTriAFracs(nVertices);
+
+        for (label v = 0; v < nVertices - 1; v++)
+        {
+            const point& vA = mesh.points()[faceVertices[v]];
+
+            const point& vB = mesh.points()[faceVertices[(v + 1)]];
+
+            cTriAFracs[v] =
+                0.5*mag((vA - fC)^(vB - fC))/fA
+              + cTriAFracs[max((v - 1), 0)];
+        }
+
+        // Force the last area fraction value to 1.0 to avoid any
+        // rounding/non-flat face errors giving a value < 1.0
+        cTriAFracs[nVertices - 1] = 1.0;
+
+        // Normal unit vector *negative* so normal is pointing into the
+        // domain
+        vector nw = patch.faceAreas()[f];
+        nw /= -mag(nw);
+
+        // Wall tangential unit vector. Use the direction between the
+        // face centre and the first vertex in the list
+        vector tw1 = fC - (mesh.points()[faceVertices[0]]);
+        tw1 /= mag(tw1);
+
+        // Other tangential unit vector.  Rescaling in case face is not
+        // flat and nw and tw1 aren't perfectly orthogonal
+        vector tw2 = nw^tw1;
+        tw2 /= mag(tw2);
+
+        forAll(particleFluxAccumulators_, i)
+        {
+            scalar& faceAccumulator = particleFluxAccumulators_[i][f];
+
+            // Number of particles to insert
+            label nI = max(label(faceAccumulator), 0);
+
+            faceAccumulator -= nI;
+
+            label typeId = moleculeTypeIds_[i];
+
+            scalar mass = cloud.constProps(typeId).mass();
+
+            for (label n = 0; n < nI; n++)
+            {
+                // Choose a triangle to insert on, based on their relative area
+
+                scalar triSelection = rndGen.scalar01();
+
+                // Selected triangle
+                label sTri = -1;
+
+                forAll(cTriAFracs, tri)
+                {
+                    sTri = tri;
+
+                    if (cTriAFracs[tri] >= triSelection)
+                    {
+                        break;
+                    }
+                }
+
+                // Randomly distribute the points on the triangle, using the
+                // method from:
+                // Generating Random Points in Triangles
+                // by Greg Turk
+                // from "Graphics Gems", Academic Press, 1990
+                // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
+
+                const point& A = fC;
+                const point& B = mesh.points()[faceVertices[sTri]];
+                const point& C =
+                    mesh.points()[faceVertices[(sTri + 1) % nVertices]];
+
+                scalar s = rndGen.scalar01();
+                scalar t = sqrt(rndGen.scalar01());
+
+                point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
+
+                vector U =
+                    sqrt(CloudType::kb*temperature_/mass)
+                   *(
+                        rndGen.GaussNormal()*tw1
+                      + rndGen.GaussNormal()*tw2
+                      - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+                    );
+
+                U += velocity_;
+
+                scalar Ei = cloud.equipartitionInternalEnergy
+                (
+                    temperature_,
+                    cloud.constProps(typeId).internalDegreesOfFreedom()
+                );
+
+                cloud.addNewParcel
+                (
+                    p,
+                    U,
+                    Ei,
+                    cell,
+                    typeId
+                );
+
+                particlesInserted++;
+            }
+        }
+    }
+
+    reduce(particlesInserted, sumOp<label>());
+
+    Info<< "    Particles inserted              = "
+        << particlesInserted << endl;
+
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
new file mode 100644
index 0000000000..f1bd4cc321
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::FreeStream
+
+Description
+    Inserting new particles across the faces of a specified patch for a free
+    stream.  Uniform values of temperature, velocity and number densities
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FreeStream_H
+#define FreeStream_H
+
+#include "InflowBoundaryModel.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                        Class FreeStream Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class FreeStream
+:
+    public InflowBoundaryModel<CloudType>
+{
+    // Private data
+
+        //- Index of patch to introduce particles across
+        label patchIndex_;
+
+        //- Temperature of the free stream
+        scalar temperature_;
+
+        //- Velocity of the free stream
+        vector velocity_;
+
+        //- The molecule types to be introduced
+        List<label> moleculeTypeIds_;
+
+        //- The number density of the species in the inflow
+        Field<scalar> numberDensities_;
+
+        //- A List of Fields, one Field for every species to be introduced, each
+        // field entry corresponding to a face on the patch to be injected
+        // across.
+        List<Field<scalar> > particleFluxAccumulators_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("FreeStream");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        FreeStream
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~FreeStream();
+
+
+    // Member Functions
+
+        //- Introduce particles
+        virtual void inflow();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "FreeStream.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/InflowBoundaryModel.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/InflowBoundaryModel.C
new file mode 100644
index 0000000000..8e99e3f43f
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/InflowBoundaryModel.C
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "InflowBoundaryModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InflowBoundaryModel<CloudType>::InflowBoundaryModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InflowBoundaryModel<CloudType>::~InflowBoundaryModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+const CloudType& Foam::InflowBoundaryModel<CloudType>::owner() const
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+CloudType& Foam::InflowBoundaryModel<CloudType>::owner()
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary& Foam::InflowBoundaryModel<CloudType>::dict() const
+{
+    return dict_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary& Foam::InflowBoundaryModel<CloudType>::coeffDict() const
+{
+    return coeffDict_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NewInflowBoundaryModel.C"
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/InflowBoundaryModel.H b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/InflowBoundaryModel.H
new file mode 100644
index 0000000000..6b01a82f40
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/InflowBoundaryModel.H
@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::InflowBoundaryModel
+
+
+Description
+    Templated inflow boundary model class
+
+SourceFiles
+    InflowBoundaryModel.C
+    NewInflowBoundaryModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef InflowBoundaryModel_H
+#define InflowBoundaryModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class InflowBoundaryModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class InflowBoundaryModel
+{
+    // Private data
+
+        //- The cloud dictionary
+        const dictionary& dict_;
+
+        // Reference to the owner cloud class
+        CloudType& owner_;
+
+        //- The coefficients dictionary
+        const dictionary coeffDict_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("InflowBoundaryModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        InflowBoundaryModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        InflowBoundaryModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+
+    //- Destructor
+    virtual ~InflowBoundaryModel();
+
+
+    //- Selector
+    static autoPtr<InflowBoundaryModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Access
+
+        //- Return const access the owner cloud object
+        inline const CloudType& owner() const;
+
+        //- Return non-const access the owner cloud object for manipulation
+        inline CloudType& owner();
+
+        //- Return the owner cloud dictionary
+        inline const dictionary& dict() const;
+
+        //- Return the coefficients dictionary
+        inline const dictionary& coeffDict() const;
+
+    //- Introduce particles
+    virtual void inflow() = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeInflowBoundaryModel(CloudType)                                    \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(InflowBoundaryModel<CloudType>, 0);   \
+                                                                              \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        InflowBoundaryModel<CloudType>,                                       \
+        dictionary                                                            \
+    );
+
+
+#define makeInflowBoundaryModelType(SS, CloudType, ParcelType)                \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(SS<CloudType<ParcelType> >, 0);       \
+                                                                              \
+    InflowBoundaryModel<CloudType<ParcelType> >::                             \
+        adddictionaryConstructorToTable<SS<CloudType<ParcelType> > >          \
+            add##SS##CloudType##ParcelType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "InflowBoundaryModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/NewInflowBoundaryModel.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/NewInflowBoundaryModel.C
new file mode 100644
index 0000000000..3abc5b2936
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/InflowBoundaryModel/NewInflowBoundaryModel.C
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "InflowBoundaryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::InflowBoundaryModel<CloudType> >
+Foam::InflowBoundaryModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word InflowBoundaryModelType(dict.lookup("InflowBoundaryModel"));
+
+    Info<< "Selecting InflowBoundaryModel " << InflowBoundaryModelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(InflowBoundaryModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "InflowBoundaryModel<CloudType>::New"
+            "(const dictionary&, CloudType&)"
+        )   << "Unknown InflowBoundaryModelType type "
+            << InflowBoundaryModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid InflowBoundaryModel types are:" << nl
+            << dictionaryConstructorTablePtr_->toc() << exit(FatalError);
+    }
+
+    return autoPtr<InflowBoundaryModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/NoInflow/NoInflow.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/NoInflow/NoInflow.C
new file mode 100644
index 0000000000..d04db500da
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/NoInflow/NoInflow.C
@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "NoInflow.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoInflow<CloudType>::NoInflow
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    InflowBoundaryModel<CloudType>(dict, cloud, typeName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoInflow<CloudType>::~NoInflow()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::NoInflow<CloudType>::inflow()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/NoInflow/NoInflow.H b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/NoInflow/NoInflow.H
new file mode 100644
index 0000000000..e3f46622e9
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/NoInflow/NoInflow.H
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::NoInflow
+
+Description
+    Not inserting any particles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoInflow_H
+#define NoInflow_H
+
+#include "InflowBoundaryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                          Class NoInflow Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoInflow
+:
+    public InflowBoundaryModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("NoInflow");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        NoInflow
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~NoInflow();
+
+
+    // Member Functions
+
+        //- Introduce particles (none in this case)
+        virtual void inflow();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoInflow.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
new file mode 100644
index 0000000000..c2c12ce5ca
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "MaxwellianThermal.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::MaxwellianThermal<CloudType>::MaxwellianThermal
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    WallInteractionModel<CloudType>(dict, cloud, typeName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::MaxwellianThermal<CloudType>::~MaxwellianThermal()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::MaxwellianThermal<CloudType>::correct
+(
+    const wallPolyPatch& wpp,
+    const label faceId,
+    vector& U,
+    scalar& Ei,
+    label typeId
+)
+{
+    label wppIndex = wpp.index();
+
+    label wppLocalFace = wpp.whichFace(faceId);
+
+    vector nw = wpp.faceAreas()[wppLocalFace];
+
+    // Normal unit vector
+    nw /= mag(nw);
+
+    // Normal velocity magnitude
+    scalar magUn = U & nw;
+
+    // Wall tangential velocity (flow direction)
+    vector Ut = U - magUn*nw;
+
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    while (mag(Ut) < SMALL)
+    {
+        // If the incident velocity is parallel to the face normal, no
+        // tangential direction can be chosen.  Add a perturbation to the
+        // incoming velocity and recalculate.
+
+        U = vector
+        (
+            U.x()*(0.8 + 0.2*rndGen.scalar01()),
+            U.y()*(0.8 + 0.2*rndGen.scalar01()),
+            U.z()*(0.8 + 0.2*rndGen.scalar01())
+        );
+
+        magUn = U & nw;
+
+        Ut = U - magUn*nw;
+    }
+
+    // Wall tangential unit vector
+    vector tw1 = Ut/mag(Ut);
+
+    // Other tangential unit vector
+    vector tw2 = nw^tw1;
+
+    scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
+
+    scalar mass = cloud.constProps(typeId).mass();
+
+    scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
+
+    U =
+        sqrt(CloudType::kb*T/mass)
+       *(
+            rndGen.GaussNormal()*tw1
+          + rndGen.GaussNormal()*tw2
+          - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+        );
+
+    U += cloud.U().boundaryField()[wppIndex][wppLocalFace];
+
+    Ei = cloud.equipartitionInternalEnergy(T, iDof);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H
new file mode 100644
index 0000000000..2b71beaba0
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::MaxwellianThermal
+
+Description
+    Wall interaction setting microscopic velocity to a random one drawn from a
+    Maxwellian distribution corresponding to a specified temperature
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MaxwellianThermal_H
+#define MaxwellianThermal_H
+
+#include "WallInteractionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                      Class MaxwellianThermal Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class MaxwellianThermal
+:
+    public WallInteractionModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("MaxwellianThermal");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        MaxwellianThermal
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~MaxwellianThermal();
+
+
+    // Member Functions
+
+        //- Apply wall correction
+        virtual void correct
+        (
+            const wallPolyPatch& wpp,
+            const label faceId,
+            vector& U,
+            scalar& Ei,
+            label typeId
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MaxwellianThermal.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
new file mode 100644
index 0000000000..4de5507bf0
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "SpecularReflection.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::SpecularReflection<CloudType>::SpecularReflection
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    WallInteractionModel<CloudType>(dict, cloud, typeName)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::SpecularReflection<CloudType>::~SpecularReflection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::SpecularReflection<CloudType>::correct
+(
+    const wallPolyPatch& wpp,
+    const label faceId,
+    vector& U,
+    scalar& Ei,
+    label typeId
+)
+{
+    vector nw = wpp.faceAreas()[wpp.whichFace(faceId)];
+    nw /= mag(nw);
+
+    scalar magUn = U & nw;
+
+    if (magUn > 0.0)
+    {
+        U -= 2.0*magUn*nw;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H
new file mode 100644
index 0000000000..d2208bd350
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::SpecularReflection
+
+Description
+    Reversing the wall-normal component of the particle velocity
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SpecularReflection_H
+#define SpecularReflection_H
+
+#include "WallInteractionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                      Class SpecularReflection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class SpecularReflection
+:
+    public WallInteractionModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("SpecularReflection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        SpecularReflection
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~SpecularReflection();
+
+
+    // Member Functions
+
+        //- Apply wall correction
+        virtual void correct
+        (
+            const wallPolyPatch& wpp,
+            const label faceId,
+            vector& U,
+            scalar& Ei,
+            label typeId
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SpecularReflection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/NewWallInteractionModel.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/NewWallInteractionModel.C
new file mode 100644
index 0000000000..85652f2379
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/NewWallInteractionModel.C
@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "WallInteractionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::WallInteractionModel<CloudType> >
+Foam::WallInteractionModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word WallInteractionModelType(dict.lookup("WallInteractionModel"));
+
+    Info<< "Selecting WallInteractionModel " << WallInteractionModelType
+        << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(WallInteractionModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "WallInteractionModel<CloudType>::New"
+            "(const dictionary&, CloudType&)"
+        )
+            << "Unknown WallInteractionModelType type "
+            << WallInteractionModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid WallInteractionModel types are:" << nl
+            << dictionaryConstructorTablePtr_->toc() << exit(FatalError);
+    }
+
+    return autoPtr<WallInteractionModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.C
new file mode 100644
index 0000000000..060d5dd76d
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.C
@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "WallInteractionModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::WallInteractionModel<CloudType>::WallInteractionModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::WallInteractionModel<CloudType>::~WallInteractionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+const CloudType& Foam::WallInteractionModel<CloudType>::owner() const
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+CloudType& Foam::WallInteractionModel<CloudType>::owner()
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary& Foam::WallInteractionModel<CloudType>::dict() const
+{
+    return dict_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary&
+Foam::WallInteractionModel<CloudType>::coeffDict() const
+{
+    return coeffDict_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NewWallInteractionModel.C"
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H
new file mode 100644
index 0000000000..582bf2b106
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H
@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::WallInteractionModel
+
+Description
+    Templated wall interaction model class
+
+SourceFiles
+    WallInteractionModel.C
+    NewWallInteractionModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef WallInteractionModel_H
+#define WallInteractionModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class WallInteractionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class WallInteractionModel
+{
+    // Private data
+
+        //- The cloud dictionary
+        const dictionary& dict_;
+
+        // reference to the owner cloud class
+        CloudType& owner_;
+
+        //- The coefficients dictionary
+        const dictionary coeffDict_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("WallInteractionModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        WallInteractionModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct from components
+        WallInteractionModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+
+    // Destructor
+    virtual ~WallInteractionModel();
+
+
+    //- Selector
+    static autoPtr<WallInteractionModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Access
+
+        //- Return the owner cloud object
+        const CloudType& owner() const;
+
+        //- Return non-const access to the owner cloud object
+        CloudType& owner();
+
+        //- Return the dictionary
+        const dictionary& dict() const;
+
+        //- Return the coefficients dictionary
+        const dictionary& coeffDict() const;
+
+
+    // Member Functions
+
+        //- Apply wall correction
+        virtual void correct
+        (
+            const wallPolyPatch& wpp,
+            const label faceId,
+            vector& U,
+            scalar& Ei,
+            label typeId
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeWallInteractionModel(CloudType)                                   \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(WallInteractionModel<CloudType>, 0);  \
+                                                                              \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        WallInteractionModel<CloudType>,                                      \
+        dictionary                                                            \
+    );
+
+
+#define makeWallInteractionModelType(SS, CloudType, ParcelType)               \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(SS<CloudType<ParcelType> >, 0);       \
+                                                                              \
+    WallInteractionModel<CloudType<ParcelType> >::                            \
+        adddictionaryConstructorToTable<SS<CloudType<ParcelType> > >          \
+            add##SS##CloudType##ParcelType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "WallInteractionModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
index 8f047d1d41..24c3fd53d0 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
@@ -308,6 +308,7 @@ Foam::fieldAverage::fieldAverage
     active_(true),
     prevTimeIndex_(-1),
     cleanRestart_(false),
+    resetOnOutput_(false),
     faItems_(),
     meanScalarFields_(),
     meanVectorFields_(),
@@ -355,6 +356,7 @@ void Foam::fieldAverage::read(const dictionary& dict)
     if (active_)
     {
         dict.readIfPresent("cleanRestart", cleanRestart_);
+        dict.readIfPresent("resetOnOutput", resetOnOutput_);
         dict.lookup("fields") >> faItems_;
 
         initialize();
@@ -387,6 +389,17 @@ void Foam::fieldAverage::write()
         calcAverages();
         writeAverages();
         writeAveragingProperties();
+
+        if (resetOnOutput_)
+        {
+            Info<< "fieldAverage: restarting averaging at time "
+                << obr_.time().timeName() << nl << endl;
+
+            initialize();
+
+            // ensure first averaging works unconditionally
+            prevTimeIndex_ = -1;
+        }
     }
 }
 
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
index bffb23b2e0..29ac4c2a3b 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
@@ -41,6 +41,11 @@ Description
         // averaging info if available
         cleanRestart true;
 
+        // Whether to reset the averaged fields after they have been written.
+        // Used to average over only the preceding write interval for transient
+        // cases.
+        resetOnOutput true;
+
         // Fields to be averaged. runTime modifiable!
         fields
         (
@@ -132,6 +137,9 @@ protected:
         //- Clean restart flag
         Switch cleanRestart_;
 
+        //- resetOnOutput flag
+        Switch resetOnOutput_;
+
         //- List of field average items, describing what averages to be
         //  calculated and output
         List<fieldAverageItem> faItems_;
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index 3fb9cc135a..1ffc36cab0 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.C
@@ -151,6 +151,8 @@ Foam::forces::forces
     patchSet_(),
     pName_(""),
     UName_(""),
+    directForceDensity_(false),
+    fDName_(""),
     rhoRef_(0),
     CofR_(vector::zero),
     forcesFilePtr_(NULL)
@@ -161,7 +163,13 @@ Foam::forces::forces
         active_ = false;
         WarningIn
         (
-            "forces::forces(const objectRegistry& obr, const dictionary& dict)"
+            "Foam::forces::forces"
+            "("
+                "const word&, "
+                "const objectRegistry&, "
+                "const dictionary&, "
+                "const bool"
+            ")"
         )   << "No fvMesh available, deactivating."
             << endl;
     }
@@ -189,27 +197,50 @@ void Foam::forces::read(const dictionary& dict)
         patchSet_ =
             mesh.boundaryMesh().patchSet(wordList(dict.lookup("patches")));
 
-        // Optional entries U and p
-        pName_ = dict.lookupOrDefault<word>("pName", "p");
-        UName_ = dict.lookupOrDefault<word>("UName", "U");
+        dict.readIfPresent("directForceDensity", directForceDensity_);
 
-        // Check whether UName and pName exists, if not deactivate forces
-        if
-        (
-            !obr_.foundObject<volVectorField>(UName_)
-         || !obr_.foundObject<volScalarField>(pName_)
-        )
+        if (directForceDensity_)
         {
-            active_ = false;
-            WarningIn("void forces::read(const dictionary& dict)")
-                << "Could not find " << UName_ << " or "
-                << pName_ << " in database." << nl
-                << "    De-activating forces."
-                << endl;
-        }
+            // Optional entry for fDName
+            fDName_ = dict.lookupOrDefault<word>("fDName", "fD");
 
-        // Reference density needed for incompressible calculations
-        rhoRef_ = readScalar(dict.lookup("rhoInf"));
+            // Check whether fDName exists, if not deactivate forces
+            if
+            (
+                !obr_.foundObject<volVectorField>(fDName_)
+            )
+            {
+                active_ = false;
+                WarningIn("void forces::read(const dictionary& dict)")
+                << "Could not find " << fDName_ << " in database." << nl
+                    << "    De-activating forces."
+                    << endl;
+            }
+        }
+        else
+        {
+            // Optional entries U and p
+            pName_ = dict.lookupOrDefault<word>("pName", "p");
+            UName_ = dict.lookupOrDefault<word>("UName", "U");
+
+            // Check whether UName and pName exists, if not deactivate forces
+            if
+            (
+                !obr_.foundObject<volVectorField>(UName_)
+             || !obr_.foundObject<volScalarField>(pName_)
+            )
+            {
+                active_ = false;
+                WarningIn("void forces::read(const dictionary& dict)")
+                << "Could not find " << UName_ << " or "
+                    << pName_ << " in database." << nl
+                    << "    De-activating forces."
+                    << endl;
+            }
+
+            // Reference density needed for incompressible calculations
+            rhoRef_ = readScalar(dict.lookup("rhoInf"));
+        }
 
         // Centre of rotation for moment calculations
         CofR_ = dict.lookup("CofR");
@@ -307,40 +338,76 @@ void Foam::forces::write()
 
 Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
 {
-    const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
-    const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
-
-    const fvMesh& mesh = U.mesh();
-
     forcesMoments fm
     (
         pressureViscous(vector::zero, vector::zero),
         pressureViscous(vector::zero, vector::zero)
     );
 
-    const surfaceVectorField::GeometricBoundaryField& Sfb =
-        mesh.Sf().boundaryField();
-
-    tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
-    const volSymmTensorField::GeometricBoundaryField& devRhoReffb
-        = tdevRhoReff().boundaryField();
-
-    forAllConstIter(labelHashSet, patchSet_, iter)
+    if (directForceDensity_)
     {
-        label patchi = iter.key();
+        const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
 
-        vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
+        const fvMesh& mesh = fD.mesh();
 
-        vectorField pf =
-            mesh.Sf().boundaryField()[patchi]*p.boundaryField()[patchi];
+        const surfaceVectorField::GeometricBoundaryField& Sfb =
+            mesh.Sf().boundaryField();
 
-        fm.first().first() += rho(p)*sum(pf);
-        fm.second().first() += rho(p)*sum(Md ^ pf);
+        forAllConstIter(labelHashSet, patchSet_, iter)
+        {
+            label patchi = iter.key();
 
-        vectorField vf = Sfb[patchi] & devRhoReffb[patchi];
+            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
 
-        fm.first().second() += sum(vf);
-        fm.second().second() += sum(Md ^ vf);
+            scalarField sA = mag(Sfb[patchi]);
+
+            // Normal force = surfaceUnitNormal * (surfaceNormal & forceDensity)
+            vectorField fN =
+                Sfb[patchi]/sA
+               *(
+                    Sfb[patchi] & fD.boundaryField()[patchi]
+                );
+
+            fm.first().first() += sum(fN);
+            fm.second().first() += sum(Md ^ fN);
+
+            // Tangential force (total force minus normal fN)
+            vectorField fT = sA*fD.boundaryField()[patchi] - fN;
+
+            fm.first().second() += sum(fT);
+            fm.second().second() += sum(Md ^ fT);
+        }
+    }
+    else
+    {
+        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
+        const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
+
+        const fvMesh& mesh = U.mesh();
+
+        const surfaceVectorField::GeometricBoundaryField& Sfb =
+            mesh.Sf().boundaryField();
+
+        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
+        const volSymmTensorField::GeometricBoundaryField& devRhoReffb
+            = tdevRhoReff().boundaryField();
+
+        forAllConstIter(labelHashSet, patchSet_, iter)
+        {
+            label patchi = iter.key();
+
+            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
+
+            vectorField pf = Sfb[patchi]*p.boundaryField()[patchi];
+
+            fm.first().first() += rho(p)*sum(pf);
+            fm.second().first() += rho(p)*sum(Md ^ pf);
+
+            vectorField vf = Sfb[patchi] & devRhoReffb[patchi];
+
+            fm.first().second() += sum(vf);
+            fm.second().second() += sum(Md ^ vf);
+        }
     }
 
     reduce(fm, sumOp());
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.H b/src/postProcessing/functionObjects/forces/forces/forces.H
index 73edae9da1..96651af97e 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.H
+++ b/src/postProcessing/functionObjects/forces/forces/forces.H
@@ -121,7 +121,7 @@ protected:
         //- Switch to send output to Info as well as to file
         Switch log_;
 
-        // Read from dictonary
+        // Read from dictionary
 
             //- Patches to integrate forces over
             labelHashSet patchSet_;
@@ -132,6 +132,12 @@ protected:
             //- Name of velocity field
             word UName_;
 
+            //- Is the force density being supplied directly?
+            Switch directForceDensity_;
+
+            //- The name of the force density (fD) field
+            word fDName_;
+
             //- Reference density needed for incompressible calculations
             scalar rhoRef_;
 
diff --git a/src/postProcessing/functionObjects/utilities/Make/files b/src/postProcessing/functionObjects/utilities/Make/files
index ae95721ac8..2d37c7c8f3 100644
--- a/src/postProcessing/functionObjects/utilities/Make/files
+++ b/src/postProcessing/functionObjects/utilities/Make/files
@@ -1,4 +1,7 @@
 dynamicPressure/dynamicPressure.C
 dynamicPressure/dynamicPressureFunctionObject.C
 
+dsmcFields/dsmcFields.C
+dsmcFields/dsmcFieldsFunctionObject.C
+
 LIB = $(FOAM_LIBBIN)/libutilityFunctionObjects
diff --git a/src/postProcessing/functionObjects/utilities/Make/options b/src/postProcessing/functionObjects/utilities/Make/options
index 5166bcc9e3..623013bc21 100644
--- a/src/postProcessing/functionObjects/utilities/Make/options
+++ b/src/postProcessing/functionObjects/utilities/Make/options
@@ -1,9 +1,13 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -lsampling
+    -lsampling \
+    -llagrangian \
+    -ldsmc
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/IOdsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/IOdsmcFields.C
new file mode 100644
index 0000000000..32627acb3c
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/IOdsmcFields.C
@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Typedef
+    Foam::IOdsmcFields
+
+Description
+    Instance of the generic IOOutputFilter for dsmcFields.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOdsmcFields_H
+#define IOdsmcFields_H
+
+#include "dsmcFields.H"
+#include "IOOutputFilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOOutputFilter<dsmcFields> IOdsmcFields;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
new file mode 100644
index 0000000000..3aaa9c514a
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcFields.H"
+#include "volFields.H"
+#include "dictionary.H"
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dsmcFields, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dsmcFields::dsmcFields
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    name_(name),
+    obr_(obr),
+    active_(true)
+{
+    // Check if the available mesh is an fvMesh, otherwise deactivate
+    if (!isA<fvMesh>(obr_))
+    {
+        active_ = false;
+        WarningIn
+        (
+            "dsmcFields::dsmcFields"
+            "(const objectRegistry&, const dictionary&)"
+        )   << "No fvMesh available, deactivating." << nl
+            << endl;
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dsmcFields::~dsmcFields()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::dsmcFields::read(const dictionary& dict)
+{
+    if (active_)
+    {
+
+    }
+}
+
+
+void Foam::dsmcFields::execute()
+{
+    // Do nothing - only valid on write
+}
+
+
+void Foam::dsmcFields::end()
+{
+    // Do nothing - only valid on write
+}
+
+
+void Foam::dsmcFields::write()
+{
+    if (active_)
+    {
+        word rhoNMeanName = "rhoNMean";
+        word rhoMMeanName = "rhoMMean";
+        word momentumMeanName = "momentumMean";
+        word linearKEMeanName = "linearKEMean";
+        word internalEMeanName = "internalEMean";
+        word iDofMeanName = "iDofMean";
+
+        const volScalarField& rhoNMean = obr_.lookupObject<volScalarField>
+        (
+            rhoNMeanName
+        );
+
+        const volScalarField& rhoMMean = obr_.lookupObject<volScalarField>
+        (
+            rhoMMeanName
+        );
+
+        const volVectorField& momentumMean = obr_.lookupObject<volVectorField>
+        (
+            momentumMeanName
+        );
+
+        const volScalarField& linearKEMean = obr_.lookupObject<volScalarField>
+        (
+            linearKEMeanName
+        );
+
+        const volScalarField& internalEMean = obr_.lookupObject<volScalarField>
+        (
+            internalEMeanName
+        );
+
+        const volScalarField& iDofMean = obr_.lookupObject<volScalarField>
+        (
+            iDofMeanName
+        );
+
+        if (min(rhoNMean).value() > VSMALL)
+        {
+            Info<< "Calculating dsmcFields." << endl;
+
+            Info<< "    Calculating UMean field." << endl;
+            volVectorField UMean
+            (
+                IOobject
+                (
+                    "UMean",
+                    obr_.time().timeName(),
+                    obr_,
+                    IOobject::NO_READ
+                ),
+                momentumMean/rhoMMean
+            );
+
+            Info<< "    Calculating translationalT field." << endl;
+            volScalarField translationalT
+            (
+                IOobject
+                (
+                    "translationalT",
+                    obr_.time().timeName(),
+                    obr_,
+                    IOobject::NO_READ
+                ),
+                2.0/(3.0*dsmcCloud::kb*rhoNMean)
+                *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
+            );
+
+            Info<< "    Calculating internalT field." << endl;
+            volScalarField internalT
+            (
+                IOobject
+                (
+                    "internalT",
+                    obr_.time().timeName(),
+                    obr_,
+                    IOobject::NO_READ
+                ),
+                2.0/(dsmcCloud::kb*iDofMean)*internalEMean
+            );
+
+            Info<< "    Calculating overallT field." << endl;
+            volScalarField overallT
+            (
+                IOobject
+                (
+                    "overallT",
+                    obr_.time().timeName(),
+                    obr_,
+                    IOobject::NO_READ
+                ),
+                2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean))
+                *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
+            );
+
+            Info<< "    mag(UMean) max/min : "
+                << max(mag(UMean)).value() << " "
+                << min(mag(UMean)).value() << endl;
+
+            Info<< "    translationalT max/min : "
+                << max(translationalT).value() << " "
+                << min(translationalT).value() << endl;
+
+            Info<< "    internalT max/min : "
+                << max(internalT).value() << " "
+                << min(internalT).value() << endl;
+
+            Info<< "    overallT max/min : "
+                << max(overallT).value() << " "
+                << min(overallT).value() << endl;
+
+            UMean.write();
+
+            translationalT.write();
+
+            internalT.write();
+
+            overallT.write();
+
+            Info<< "dsmcFields written." << nl << endl;
+        }
+        else
+        {
+            Info<< "Small or negative value (" << min(rhoNMean)
+                << ") found in rhoNMean field. "
+                << "Not calculating dsmcFields to avoid division by zero "
+                << "or invalid results."
+                << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.H b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.H
new file mode 100644
index 0000000000..b53dda6700
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.H
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::dsmcFields
+
+Description
+    Calculate intensive fields:
+    - UMean
+    - translationalT
+    - internalT
+    - overallT
+    from averaged extensive fields from a DSMC calculation.
+
+SourceFiles
+    dsmcFields.C
+    IOdsmcFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dsmcFields_H
+#define dsmcFields_H
+
+#include "pointFieldFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class dsmcFields Declaration
+\*---------------------------------------------------------------------------*/
+
+class dsmcFields
+{
+    // Private data
+
+        //- Name of this set of dsmcFields objects
+        word name_;
+
+        const objectRegistry& obr_;
+
+        //- on/off switch
+        bool active_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        dsmcFields(const dsmcFields&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dsmcFields&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dsmcFields");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        dsmcFields
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    // Destructor
+
+        virtual ~dsmcFields();
+
+
+    // Member Functions
+
+        //- Return name of the set of dsmcFields
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the dsmcFields data
+        virtual void read(const dictionary&);
+
+        //- Execute, currently does nothing
+        virtual void execute();
+
+        //- Execute at the final time-loop, currently does nothing
+        virtual void end();
+
+        //- Calculate the dsmcFields and write
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&)
+        {}
+
+        //- Update for changes of mesh
+        virtual void movePoints(const pointField&)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFieldsFunctionObject.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFieldsFunctionObject.C
new file mode 100644
index 0000000000..cd5324f25b
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFieldsFunctionObject.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "dsmcFieldsFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(dsmcFieldsFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        dsmcFieldsFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFieldsFunctionObject.H b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFieldsFunctionObject.H
new file mode 100644
index 0000000000..cf34cdaa2e
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFieldsFunctionObject.H
@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Typedef
+    Foam::dsmcFieldsFunctionObject
+
+Description
+    FunctionObject wrapper around dsmcFields to allow it to be created via
+    the functions list within controlDict.
+
+SourceFiles
+    dsmcFieldsFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dsmcFieldsFunctionObject_H
+#define dsmcFieldsFunctionObject_H
+
+#include "dsmcFields.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<dsmcFields>
+        dsmcFieldsFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //