Merge commit 'OpenCFD/master' into olesenm

applications/solvers/dsmc/dsmcFoam/Make/files

@@ -0,0 +1,3 @@
+dsmcFoam.C
+
+EXE = $(FOAM_APPBIN)/dsmcFoam

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
+

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);

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);
+}
+
+
+// ************************************************************************* //

applications/test/List/ListTest.C

@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
     List<vector> list(IStringStream("1 ((0 1 2))")());
     Info<< list << endl;
 
-    List<vector> list2(IStringStream("((0 1 2) (3 4 5) (3 4 5))")());
+    List<vector> list2(IStringStream("((0 1 2) (3 4 5) (6 7 8))")());
     Info<< list2 << endl;
 
     Info<< findIndex(list2, vector(3, 4, 5)) << endl;
@@ -59,6 +59,13 @@ int main(int argc, char *argv[])
     Info<< list2 << nl
         << list3 << endl;
 
+
+    // Subset
+    const labelList map(IStringStream("2 (0 2)")());
+    List<vector> subList3(list3, map);
+    Info<< "Elements " << map << " out of " << list3
+        << " : " << subList3 << endl;
+
     return 0;
 }
 

applications/test/UIndirectListTest/Make/files

@@ -0,0 +1,3 @@
+UIndirectListTest.C
+
+EXE = $(FOAM_USER_APPBIN)/UIndirectListTest

applications/test/UIndirectListTest/UIndirectListTest.C

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "UIndirectList.H"
+#include "IOstreams.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    List<double> completeList(10);
+
+    forAll(completeList, i)
+    {
+        completeList[i] = 0.1*i;
+    }
+
+    List<label> addresses(5);
+    addresses[0] = 1;
+    addresses[1] = 0;
+    addresses[2] = 7;
+    addresses[3] = 8;
+    addresses[4] = 5;
+
+    UIndirectList<double> idl(completeList, addresses);
+
+    forAll(idl, i)
+    {
+        Info<< idl[i] << token::SPACE;
+    }
+
+    Info<< endl;
+
+    idl[1] = -666;
+
+    Info<< "idl[1] changed:" << idl() << endl;
+
+    idl = -999;
+
+    Info<< "idl changed:" << idl() << endl;
+
+    UIndirectList<double> idl2(idl);
+
+    Info<< "idl2:" << idl2() << endl;
+
+    idl = idl2();
+
+    Info<< "idl assigned from UList:" << idl() << endl;
+
+
+    List<double> realList = UIndirectList<double>(completeList, addresses);
+
+    Info<< "realList:" << realList << endl;
+
+    List<double> realList2(UIndirectList<double>(completeList, addresses));
+
+    Info<< "realList2:" << realList2 << endl;
+
+
+    Info << "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C

@@ -177,7 +177,7 @@ label mergePatchFaces
         List<faceList> allFaceSetsFaces(allFaceSets.size());
         forAll(allFaceSets, setI)
         {
-            allFaceSetsFaces[setI] = IndirectList<face>
+            allFaceSetsFaces[setI] = UIndirectList<face>
             (
                 mesh.faces(),
                 allFaceSets[setI]

applications/utilities/mesh/conversion/polyDualMesh/meshDualiser.C

@@ -76,7 +76,7 @@ void Foam::meshDualiser::checkPolyTopoChange(const polyTopoChange& meshMod)
                     "meshDualiser::checkPolyTopoChange(const polyTopoChange&)"
                 )   << "duplicate verts:" << newToOld[newI]
                     << " coords:"
-                    << IndirectList<point>(points, newToOld[newI])()
+                    << UIndirectList<point>(points, newToOld[newI])()
                     << abort(FatalError);
             }
         }
@@ -226,10 +226,7 @@ Foam::label Foam::meshDualiser::addInternalFace
 
     if (debug)
     {
-        pointField facePoints
-        (
-            IndirectList<point>(meshMod.points(), newFace)()
-        );
+        pointField facePoints(meshMod.points(), newFace);
 
         labelList oldToNew;
         pointField newPoints;
@@ -289,7 +286,7 @@ Foam::label Foam::meshDualiser::addInternalFace
         //n /= mag(n);
         //Pout<< "Generated internal dualFace:" << dualFaceI
         //    << " verts:" << newFace
-        //    << " points:" << IndirectList<point>(meshMod.points(), newFace)()
+        //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
         //    << " n:" << n
         //    << " between dualowner:" << dualCell0
         //    << " dualneigbour:" << dualCell1
@@ -316,7 +313,7 @@ Foam::label Foam::meshDualiser::addInternalFace
         //n /= mag(n);
         //Pout<< "Generated internal dualFace:" << dualFaceI
         //    << " verts:" << newFace
-        //    << " points:" << IndirectList<point>(meshMod.points(), newFace)()
+        //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
         //    << " n:" << n
         //    << " between dualowner:" << dualCell1
         //    << " dualneigbour:" << dualCell0
@@ -373,7 +370,7 @@ Foam::label Foam::meshDualiser::addBoundaryFace
     //n /= mag(n);
     //Pout<< "Generated boundary dualFace:" << dualFaceI
     //    << " verts:" << newFace
-    //    << " points:" << IndirectList<point>(meshMod.points(), newFace)()
+    //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)()
     //    << " n:" << n
     //    << " on dualowner:" << dualCellI
     //    << endl;
@@ -568,7 +565,7 @@ void Foam::meshDualiser::createFaceFromInternalFace
 
     //Pout<< "createFaceFromInternalFace : At face:" << faceI
     //    << " verts:" << f
-    //    << " points:" << IndirectList<point>(mesh_.points(), f)()
+    //    << " points:" << UIndirectList<point>(mesh_.points(), f)()
     //    << " started walking at edge:" << fEdges[fp]
     //    << " verts:" << mesh_.edges()[fEdges[fp]]
     //    << endl;
@@ -617,7 +614,7 @@ void Foam::meshDualiser::createFaceFromInternalFace
             {
                 FatalErrorIn("createFacesFromInternalFace(..)")
                     << "face:" << faceI << " verts:" << f
-                    << " points:" << IndirectList<point>(mesh_.points(), f)()
+                    << " points:" << UIndirectList<point>(mesh_.points(), f)()
                     << " no feature edge between " << f[fp]
                     << " and " << f[nextFp] << " although have different"
                     << " dual cells." << endl

applications/utilities/mesh/manipulation/setSet/writePointSet.C

@@ -70,7 +70,7 @@ void writePointSet
 
     labelList pointLabels(set.toc());
 
-    pointField setPoints(IndirectList<point>(mesh.points(), pointLabels)());
+    pointField setPoints(mesh.points(), pointLabels);
 
     // Write points
 

applications/utilities/mesh/manipulation/splitMesh/regionSide.C

@@ -22,14 +22,11 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "regionSide.H"
 #include "meshTools.H"
 #include "primitiveMesh.H"
-#include "IndirectList.H"
 
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -113,7 +110,7 @@ Foam::label Foam::regionSide::otherEdge
     )   << "Cannot find other edge on face " << faceI << " that uses point "
         << pointI << " but not point " << freePointI << endl
         << "Edges on face:" << fEdges
-        << " verts:" << IndirectList<edge>(mesh.edges(), fEdges)()
+        << " verts:" << UIndirectList<edge>(mesh.edges(), fEdges)()
         << " Vertices on face:"
         << mesh.faces()[faceI]
         << " Vertices on original edge:" << e << abort(FatalError);

applications/utilities/mesh/manipulation/splitMesh/splitMesh.C

@@ -50,15 +50,12 @@ Description
 #include "attachDetach.H"
 #include "attachPolyTopoChanger.H"
 #include "regionSide.H"
+#include "primitiveFacePatch.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Calculation engine for set of faces in a mesh
-typedef PrimitivePatch<face, List, const pointField&> facePatch;
-
-
 // Find edge between points v0 and v1.
 label findEdge(const primitiveMesh& mesh, const label v0, const label v1)
 {
@@ -163,10 +160,16 @@ int main(int argc, char *argv[])
     // set of edges on side of this region. Use PrimitivePatch to find these.
     //
 
-    IndirectList<face> zoneFaces(mesh.faces(), faces);
-
     // Addressing on faces only in mesh vertices.
-    facePatch fPatch(zoneFaces(), mesh.points());
+    primitiveFacePatch fPatch
+    (
+        UIndirectList<face>
+        (
+            mesh.faces(),
+            faces
+        ),
+        mesh.points()
+    );
 
     const labelList& meshPoints = fPatch.meshPoints();
 

applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C

@@ -60,7 +60,6 @@ Description
 #include "polyTopoChanger.H"
 #include "mapPolyMesh.H"
 #include "ListOps.H"
-#include "IndirectList.H"
 #include "slidingInterface.H"
 #include "perfectInterface.H"
 #include "IOobjectList.H"

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -35,14 +35,22 @@ Description
     -rotate (vector vector)
         Rotates the points from the first vector to the second,
 
+     or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
+     or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
+
     -scale vector
         Scales the points by the given vector.
 
     The any or all of the three options may be specified and are processed
     in the above order.
 
-    With -rotateFields (in combination with -rotate) it will also
-    read & transform vector & tensor fields.
+    With -rotateFields (in combination with -rotate/yawPitchRoll/rollPitchYaw)
+    it will also read & transform vector & tensor fields.
+
+    Note:
+    yaw (rotation about z)
+    pitch (rotation about y)
+    roll (rotation about x)
 
 \*---------------------------------------------------------------------------*/
 
@@ -58,6 +66,7 @@ Description
 #include "IStringStream.H"
 
 using namespace Foam;
+using namespace Foam::mathematicalConstant;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -131,6 +140,8 @@ int main(int argc, char *argv[])
 {
     argList::validOptions.insert("translate", "vector");
     argList::validOptions.insert("rotate", "(vector vector)");
+    argList::validOptions.insert("rollPitchYaw", "(roll pitch yaw)");
+    argList::validOptions.insert("yawPitchRoll", "(yaw pitch roll)");
     argList::validOptions.insert("rotateFields", "");
     argList::validOptions.insert("scale", "vector");
  
@@ -185,6 +196,58 @@ int main(int argc, char *argv[])
             rotateFields(runTime, T);
         }
     }
+    else if (args.options().found("rollPitchYaw"))
+    {
+        vector v(IStringStream(args.options()["rollPitchYaw"])());
+
+        Info<< "Rotating points by" << nl
+            << "    roll  " << v.x() << nl
+            << "    pitch " << v.y() << nl
+            << "    yaw   " << v.z() << endl;
+
+
+        // Convert to radians
+        v *= pi/180.0;
+
+        quaternion R(v.x(), v.y(), v.z());
+
+        Info<< "Rotating points by quaternion " << R << endl;
+        points = transform(R, points);
+
+        if (args.options().found("rotateFields"))
+        {
+            rotateFields(runTime, R.R());
+        }
+    }
+    else if (args.options().found("yawPitchRoll"))
+    {
+        vector v(IStringStream(args.options()["yawPitchRoll"])());
+
+        Info<< "Rotating points by" << nl
+            << "    yaw   " << v.x() << nl
+            << "    pitch " << v.y() << nl
+            << "    roll  " << v.z() << endl;
+
+
+        // Convert to radians
+        v *= pi/180.0;
+
+        scalar yaw = v.x();
+        scalar pitch = v.y();
+        scalar roll = v.z();
+
+        quaternion R = quaternion(vector(0, 0, 1), yaw);
+        R *= quaternion(vector(0, 1, 0), pitch);
+        R *= quaternion(vector(1, 0, 0), roll);
+
+        Info<< "Rotating points by quaternion " << R << endl;
+        points = transform(R, points);
+
+        if (args.options().found("rotateFields"))
+        {
+            rotateFields(runTime, R.R());
+        }
+    }
 
     if (args.options().found("scale"))
     {

applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "writePointSet.H"
@@ -76,7 +74,7 @@ void writePointSet
 
     writeFuns::insert
     (
-        IndirectList<point>(vMesh.mesh().points(), set.toc())(),
+        UIndirectList<point>(vMesh.mesh().points(), set.toc())(),
         ptField
     );
 

applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/Make/files

@@ -0,0 +1,3 @@
+dsmcFieldsCalc.C
+
+EXE = $(FOAM_APPBIN)/dsmcFieldsCalc

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
+

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();
+    }
+}
+
+// ************************************************************************* //

applications/utilities/preProcessing/dsmcInitialise/Make/files

@@ -0,0 +1,3 @@
+dsmcInitialise.C
+
+EXE = $(FOAM_APPBIN)/dsmcInitialise

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
+

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);
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C

@@ -306,7 +306,7 @@ label findEdge
 
     FatalErrorIn("findEdge") << "Cannot find edge with labels " << v0
         << ' ' << v1 << " in candidates " << edgeLabels
-        << " with vertices:" << IndirectList<edge>(surf.edges(), edgeLabels)()
+        << " with vertices:" << UIndirectList<edge>(surf.edges(), edgeLabels)()
         << abort(FatalError);
 
     return -1;
@@ -346,7 +346,7 @@ label otherEdge
     FatalErrorIn("otherEdge") << "Cannot find other edge on face " << faceI
         << " verts:" << surf.localPoints()[faceI]
         << " connected to point " << pointI
-        << " faceEdges:" << IndirectList<edge>(surf.edges(), fEdges)()
+        << " faceEdges:" << UIndirectList<edge>(surf.edges(), fEdges)()
         << abort(FatalError);
 
     return -1;

applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C

@@ -23,9 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Transform (scale/rotate) a surface. Like transforPoints but then for
+    Transform (scale/rotate) a surface. Like transformPoints but then for
     surfaces.
 
+    The rollPitchYaw option takes three angles (degrees):
+    - roll (rotation about x) followed by
+    - pitch (rotation about y) followed by
+    - yaw (rotation about z)
+
+    The yawPitchRoll does yaw followed by pitch followed by roll.
+
 \*---------------------------------------------------------------------------*/
 
 #include "triSurface.H"
@@ -35,8 +42,10 @@ Description
 #include "boundBox.H"
 #include "transformField.H"
 #include "Pair.H"
+#include "quaternion.H"
 
 using namespace Foam;
+using namespace Foam::mathematicalConstant;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,6 +61,8 @@ int main(int argc, char *argv[])
     argList::validOptions.insert("translate", "vector");
     argList::validOptions.insert("rotate", "(vector vector)");
     argList::validOptions.insert("scale", "vector");
+    argList::validOptions.insert("rollPitchYaw", "(roll pitch yaw)");
+    argList::validOptions.insert("yawPitchRoll", "(yaw pitch roll)");
     argList args(argc, argv);
 
     fileName surfFileName(args.additionalArgs()[0]);
@@ -96,6 +107,48 @@ int main(int argc, char *argv[])
 
         points = transform(T, points);
     }
+    else if (args.options().found("rollPitchYaw"))
+    {
+        vector v(IStringStream(args.options()["rollPitchYaw"])());
+
+        Info<< "Rotating points by" << nl
+            << "    roll  " << v.x() << nl
+            << "    pitch " << v.y() << nl
+            << "    yaw   " << v.z() << endl;
+
+
+        // Convert to radians
+        v *= pi/180.0;
+
+        quaternion R(v.x(), v.y(), v.z());
+
+        Info<< "Rotating points by quaternion " << R << endl;
+        points = transform(R, points);
+    }
+    else if (args.options().found("yawPitchRoll"))
+    {
+        vector v(IStringStream(args.options()["yawPitchRoll"])());
+
+        Info<< "Rotating points by" << nl
+            << "    yaw   " << v.x() << nl
+            << "    pitch " << v.y() << nl
+            << "    roll  " << v.z() << endl;
+
+
+        // Convert to radians
+        v *= pi/180.0;
+
+        scalar yaw = v.x();
+        scalar pitch = v.y();
+        scalar roll = v.z();
+
+        quaternion R = quaternion(vector(0, 0, 1), yaw);
+        R *= quaternion(vector(0, 1, 0), pitch);
+        R *= quaternion(vector(1, 0, 0), roll);
+
+        Info<< "Rotating points by quaternion " << R << endl;
+        points = transform(R, points);
+    }
 
     if (args.options().found("scale"))
     {

src/OpenFOAM/containers/Lists/List/List.C

@@ -31,6 +31,7 @@ License
 #include "PtrList.H"
 #include "SLList.H"
 #include "IndirectList.H"
+#include "UIndirectList.H"
 #include "BiIndirectList.H"
 #include "contiguous.H"
 
@@ -171,6 +172,29 @@ Foam::List<T>::List(List<T>& a, bool reUse)
 }
 
 
+// Construct as subset
+template<class T>
+Foam::List<T>::List(const UList<T>& a, const unallocLabelList& map)
+:
+    UList<T>(NULL, map.size())
+{
+    if (this->size_)
+    {
+        this->v_ = new T[this->size_];
+
+        List_ACCESS(T, (*this), vp);
+        List_CONST_ACCESS(T, a, ap);
+        List_FOR_ALL(map, i)
+            List_ELEM((*this), vp, i) = List_ELEM(a, ap, (map[i]));
+        List_END_FOR_ALL
+    }
+    else
+    {
+        this->v_ = 0;
+    }
+}
+
+
 // Construct given start and end iterators.
 template<class T>
 template<class InputIterator>
@@ -298,6 +322,28 @@ Foam::List<T>::List(const IndirectList<T>& lst)
 }
 
 
+// Construct as copy of UIndirectList<T>
+template<class T>
+Foam::List<T>::List(const UIndirectList<T>& lst)
+:
+    UList<T>(NULL, lst.size())
+{
+    if (this->size_)
+    {
+        this->v_ = new T[this->size_];
+
+        forAll(*this, i)
+        {
+            this->operator[](i) = lst[i];
+        }
+    }
+    else
+    {
+        this->v_ = 0;
+    }
+}
+
+
 // Construct as copy of BiIndirectList<T>
 template<class T>
 Foam::List<T>::List(const BiIndirectList<T>& lst)
@@ -567,6 +613,28 @@ void Foam::List<T>::operator=(const IndirectList<T>& lst)
 }
 
 
+// Assignment operator. Takes linear time.
+template<class T>
+void Foam::List<T>::operator=(const UIndirectList<T>& lst)
+{
+    if (lst.size() != this->size_)
+    {
+        if (this->v_) delete[] this->v_;
+        this->v_ = 0;
+        this->size_ = lst.size();
+        if (this->size_) this->v_ = new T[this->size_];
+    }
+
+    if (this->size_)
+    {
+        forAll(*this, i)
+        {
+            this->operator[](i) = lst[i];
+        }
+    }
+}
+
+
 // Assignment operator. Takes linear time.
 template<class T>
 void Foam::List<T>::operator=(const BiIndirectList<T>& lst)

src/OpenFOAM/containers/Lists/List/List.H

@@ -66,8 +66,10 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
     class DynamicList;
 template<class T> class SortableList;
 template<class T> class IndirectList;
+template<class T> class UIndirectList;
 template<class T> class BiIndirectList;
 
+typedef UList<label> unallocLabelList;
 
 /*---------------------------------------------------------------------------*\
                            Class List Declaration
@@ -112,6 +114,9 @@ public:
         //- Construct as copy or re-use as specified.
         List(List<T>&, bool reUse);
 
+        //- Construct as subset.
+        List(const UList<T>&, const unallocLabelList& mapAddressing);
+
         //- Construct given start and end iterators.
         template<class InputIterator>
         List(InputIterator first, InputIterator last);
@@ -129,6 +134,9 @@ public:
         //- Construct as copy of IndirectList<T>
         List(const IndirectList<T>&);
 
+        //- Construct as copy of UIndirectList<T>
+        List(const UIndirectList<T>&);
+
         //- Construct as copy of BiIndirectList<T>
         List(const BiIndirectList<T>&);
 
@@ -206,6 +214,9 @@ public:
         //- Assignment from IndirectList operator. Takes linear time.
         void operator=(const IndirectList<T>&);
 
+        //- Assignment from UIndirectList operator. Takes linear time.
+        void operator=(const UIndirectList<T>&);
+
         //- Assignment from BiIndirectList operator. Takes linear time.
         void operator=(const BiIndirectList<T>&);
 

src/OpenFOAM/containers/Lists/UIndirectList/UIndirectList.H

@@ -23,17 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::IndirectList
+    Foam::UIndirectList
 
 Description
+    A List with indirect addressing. Like IndirectList but does not store
+    addressing.
 
 SourceFiles
-    IndirectListI.H
+    UIndirectListI.H
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef IndirectList_H
-#define IndirectList_H
+#ifndef UIndirectList_H
+#define UIndirectList_H
 
 #include "List.H"
 
@@ -43,22 +45,16 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class IndirectList Declaration
+                           Class UIndirectList Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class T>
-class IndirectList
+class UIndirectList
 {
     // Private data
 
-        const UList<T>& completeList_;
-        List<label> addresses_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise assignment
-        void operator=(const IndirectList<T>&);
+        UList<T>& completeList_;
+        const UList<label>& addressing_;
 
 
 public:
@@ -66,7 +62,7 @@ public:
     // Constructors
 
         //- Construct given the complete list and the addressing array
-        inline IndirectList(const UList<T>&, const List<label>&);
+        inline UIndirectList(const UList<T>&, const UList<label>&);
 
 
     // Member Functions
@@ -74,14 +70,28 @@ public:
         // Access
 
             inline label size() const;
+            inline bool  empty() const;
+
             inline const UList<T>& completeList() const;
-            inline const List<label>& addresses() const;
+            inline const List<label>& addressing() const;
 
 
         // Member Operators
 
+            //- Return the addressed elements as a List
+            inline List<T> operator()() const;
+
+            //- Return non-const access to an element
+            inline T& operator[](const label);
+
             //- Return const access to an element
             inline const T& operator[](const label) const;
+
+            //- Assignment from UList of addressed elements
+            inline void operator=(const UList<T>&);
+
+            //- Assignment of all entries to the given value
+            inline void operator=(const T&);
 };
 
 
@@ -91,7 +101,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "IndirectListI.H"
+#include "UIndirectListI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/containers/Lists/UIndirectList/UIndirectListI.H

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 T>
+inline Foam::UIndirectList<T>::UIndirectList
+(
+    const UList<T>& completeList,
+    const UList<label>& addr
+)
+:
+    completeList_(const_cast<UList<T>&>(completeList)),
+    addressing_(addr)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class T>
+inline Foam::label Foam::UIndirectList<T>::size() const
+{
+    return addressing_.size();
+}
+
+
+template<class T>
+inline bool Foam::UIndirectList<T>::empty() const
+{
+    return addressing_.empty();
+}
+
+
+template<class T>
+inline const Foam::UList<T>& Foam::UIndirectList<T>::completeList() const
+{
+    return completeList_;
+}
+
+
+template<class T>
+inline const Foam::List<Foam::label>& Foam::UIndirectList<T>::addressing() const
+{
+    return addressing_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class T>
+inline Foam::List<T> Foam::UIndirectList<T>::operator()() const
+{
+    List<T> result(size());
+
+    forAll(*this, i)
+    {
+        result[i] = operator[](i);
+    }
+
+    return result;
+}
+
+
+template<class T>
+inline T& Foam::UIndirectList<T>::operator[](const label i)
+{
+    return completeList_[addressing_[i]];
+}
+
+
+template<class T>
+inline const T& Foam::UIndirectList<T>::operator[](const label i) const
+{
+    return completeList_[addressing_[i]];
+}
+
+
+template<class T>
+inline void Foam::UIndirectList<T>::operator=(const UList<T>& ae)
+{
+    if (addressing_.size() != ae.size())
+    {
+        FatalErrorIn("UIndirectList<T>::operator=(const UList<T>&)")
+            << "Addressing and list of addressed elements "
+               "have different sizes: "
+            << addressing_.size() << " " << ae.size()
+            << abort(FatalError);
+    }
+
+    forAll(addressing_, i)
+    {
+        completeList_[addressing_[i]] = ae[i];
+    }
+}
+
+
+template<class T>
+inline void Foam::UIndirectList<T>::operator=(const T& t)
+{
+    forAll(addressing_, i)
+    {
+        completeList_[addressing_[i]] = t;
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/db/objectRegistry/objectRegistry.C

@@ -50,7 +50,8 @@ Foam::objectRegistry::objectRegistry
             IOobject::NO_READ,
             IOobject::AUTO_WRITE,
             false
-        )
+        ),
+        true    // to flag that this is the top-level regIOobject
     ),
     HashTable<regIOobject*>(nIoObjects),
     time_(t),

src/OpenFOAM/db/regIOobject/regIOobject.C

@@ -41,13 +41,18 @@ int Foam::regIOobject::fileModificationSkew
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from IOobject
-Foam::regIOobject::regIOobject(const IOobject& io)
+Foam::regIOobject::regIOobject(const IOobject& io, const bool isTime)
 :
     IOobject(io),
     registered_(false),
     ownedByRegistry_(false),
     lastModified_(0),
-    eventNo_(db().getEvent()),
+    eventNo_                // Do not get event for top level Time database
+    (
+        isTime
+      ? 0
+      : db().getEvent()
+    ),
     isPtr_(NULL)
 {
     // Register with objectRegistry if requested
@@ -156,7 +161,7 @@ bool Foam::regIOobject::checkOut()
 }
 
 
-bool Foam::regIOobject::uptodate(const word& a) const
+bool Foam::regIOobject::upToDate(const word& a) const
 {
     if (db().lookupObject<regIOobject>(a).eventNo() >= eventNo_)
     {
@@ -169,7 +174,7 @@ bool Foam::regIOobject::uptodate(const word& a) const
 }
 
 
-bool Foam::regIOobject::uptodate(const word& a, const word& b) const
+bool Foam::regIOobject::upToDate(const word& a, const word& b) const
 {
     if
     (
@@ -186,7 +191,7 @@ bool Foam::regIOobject::uptodate(const word& a, const word& b) const
 }
 
 
-bool Foam::regIOobject::uptodate
+bool Foam::regIOobject::upToDate
 (
     const word& a,
     const word& b,
@@ -209,7 +214,7 @@ bool Foam::regIOobject::uptodate
 }
 
 
-bool Foam::regIOobject::uptodate
+bool Foam::regIOobject::upToDate
 (
     const word& a,
     const word& b,
@@ -235,7 +240,7 @@ bool Foam::regIOobject::uptodate
 
 
 //- Flag me as up to date
-void Foam::regIOobject::setUptodate()
+void Foam::regIOobject::setUpToDate()
 {
     eventNo_ = db().getEvent();
 }

src/OpenFOAM/db/regIOobject/regIOobject.H

@@ -99,8 +99,9 @@ public:
 
     // Constructors
 
-        //- Construct from IOobject
-        regIOobject(const IOobject&);
+        //- Construct from IOobject. Optional flag for if IOobject is the
+        //  top level regIOobject.
+        regIOobject(const IOobject&, const bool isTime = false);
 
         //- Construct as copy
         regIOobject(const regIOobject&);
@@ -154,14 +155,20 @@ public:
             inline label& eventNo();
 
             //- Am I uptodate with respect to other regIOobjects
-            bool uptodate(const word&) const;
-            bool uptodate(const word&, const word&) const;
-            bool uptodate(const word&, const word&, const word&) const;
-            bool uptodate(const word&, const word&, const word&, const word&)
-             const;
+            bool upToDate(const word&) const;
+            bool upToDate(const word&, const word&) const;
+            bool upToDate(const word&, const word&, const word&) const;
+            bool upToDate
+            (
+                const word&,
+                const word&,
+                const word&,
+                const word&
+            ) const;
 
             //- Flag me as up to date
-            void setUptodate();
+            void setUpToDate();
+
 
         // Edit
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -649,6 +649,7 @@ typename
 Foam::GeometricField<Type, PatchField, GeoMesh>::DimensionedInternalField&
 Foam::GeometricField<Type, PatchField, GeoMesh>::dimensionedInternalField()
 {
+    this->setUpToDate();
     storeOldTimes();
     return *this;
 }
@@ -659,6 +660,7 @@ typename
 Foam::GeometricField<Type, PatchField, GeoMesh>::InternalField&
 Foam::GeometricField<Type, PatchField, GeoMesh>::internalField()
 {
+    this->setUpToDate();
     storeOldTimes();
     return *this;
 }
@@ -670,6 +672,7 @@ typename
 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField&
 Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField()
 {
+    this->setUpToDate();
     storeOldTimes();
     return boundaryField_;
 }
@@ -822,6 +825,7 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::
 correctBoundaryConditions()
 {
+    this->setUpToDate();
     storeOldTimes();
     boundaryField_.evaluate();
 }

src/OpenFOAM/matrices/solution/solution.C

@@ -27,7 +27,7 @@ License
 #include "solution.H"
 #include "Time.H"
 
-// these are for old syntax compatibility:
+// These are for old syntax compatibility:
 #include "BICCG.H"
 #include "ICCG.H"
 #include "IStringStream.H"
@@ -36,7 +36,7 @@ License
 
 int Foam::solution::debug(::Foam::debug::debugSwitch("solution", 0));
 
-// list of sub-dictionaries to rewrite
+// List of sub-dictionaries to rewrite
 //! @cond localScope
 static const Foam::List<Foam::word> subDictNames
 (
@@ -59,7 +59,12 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
             IOobject::NO_WRITE
         )
     ),
-    relaxationFactors_(ITstream("relaxationFactors", tokenList())()),
+    relaxationFactors_
+    (
+        ITstream("relaxationFactors",
+        tokenList())()
+    ),
+    defaultRelaxationFactor_(0),
     solvers_(ITstream("solvers", tokenList())())
 {
     read();
@@ -157,6 +162,11 @@ bool Foam::solution::read()
             relaxationFactors_ = dict.subDict("relaxationFactors");
         }
 
+        if (relaxationFactors_.found("default"))
+        {
+            relaxationFactors_.lookup("default") >> defaultRelaxationFactor_;
+        }
+
         if (dict.found("solvers"))
         {
             solvers_ = dict.subDict("solvers");
@@ -192,7 +202,9 @@ bool Foam::solution::relax(const word& name) const
         Info<< "Find relax for " << name << endl;
     }
 
-    return relaxationFactors_.found(name);
+    return
+        relaxationFactors_.found(name)
+     || relaxationFactors_.found("default");
 }
 
 
@@ -203,8 +215,27 @@ Foam::scalar Foam::solution::relaxationFactor(const word& name) const
         Info<< "Lookup relaxationFactor for " << name << endl;
     }
 
+    if (relaxationFactors_.found(name))
+    {
         return readScalar(relaxationFactors_.lookup(name));
     }
+    else if (defaultRelaxationFactor_ > SMALL)
+    {
+        return defaultRelaxationFactor_;
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "Foam::solution::relaxationFactor(const word& name)",
+            relaxationFactors_
+        )   << "Cannot find relaxationFactor for '" << name
+            << "' or a suitable default value."
+            << exit(FatalIOError);
+
+        return 0;
+    }
+}
 
 
 const Foam::dictionary& Foam::solution::solverDict(const word& name) const

src/OpenFOAM/matrices/solution/solution.H

@@ -56,6 +56,9 @@ class solution
         //- Dictionary of relaxation factors for all the fields
         dictionary relaxationFactors_;
 
+        //- Optional default relaxation factor for all the fields
+        scalar defaultRelaxationFactor_;
+
         //- Dictionary of solver parameters for all the fields
         dictionary solvers_;
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -629,7 +629,7 @@ Foam::pointField Foam::globalMeshData::sharedPoints() const
 
             toMaster
                 << sharedPointAddr_ 
-                << IndirectList<point>(mesh_.points(), sharedPointLabels_)();
+                << UIndirectList<point>(mesh_.points(), sharedPointLabels_)();
         }
 
         // Receive sharedPoints

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

@@ -653,8 +653,8 @@ void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
             }
 
             toNeighbour
-                << IndirectList<label>(sharedPointAddr_, changedIndices)()
-                << IndirectList<label>(sharedPointLabels_, changedIndices)();
+                << UIndirectList<label>(sharedPointAddr_, changedIndices)()
+                << UIndirectList<label>(sharedPointLabels_, changedIndices)();
         }
     }
 }

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C

@@ -124,7 +124,7 @@ Foam::List<Foam::labelPair> Foam::mapDistribute::schedule
     );
 
     // Processors involved in my schedule
-    return IndirectList<labelPair>(allComms, mySchedule);
+    return UIndirectList<labelPair>(allComms, mySchedule);
 
 
     //if (debug)

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -188,11 +188,11 @@ void Foam::cyclicPolyPatch::calcTransforms()
                      << endl
                     << "Mesh face:" << start()+facei
                     << " vertices:"
-                    << IndirectList<point>(points, operator[](facei))()
+                    << UIndirectList<point>(points, operator[](facei))()
                     << endl
                     << "Neighbour face:" << start()+nbrFacei
                     << " vertices:"
-                    << IndirectList<point>(points, operator[](nbrFacei))()
+                    << UIndirectList<point>(points, operator[](nbrFacei))()
                     << endl
                     << "Rerun with cyclic debug flag set"
                     << " for more information." << exit(FatalError);
@@ -403,12 +403,12 @@ bool Foam::cyclicPolyPatch::getGeometricHalves
             fileName nm0(casePath/name()+"_half0_faces.obj");
             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
                 << " faces to OBJ file " << nm0 << endl;
-            writeOBJ(nm0, IndirectList<face>(pp, half0ToPatch)(), pp.points());
+            writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
 
             fileName nm1(casePath/name()+"_half1_faces.obj");
             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
                 << " faces to OBJ file " << nm1 << endl;
-            writeOBJ(nm1, IndirectList<face>(pp, half1ToPatch)(), pp.points());
+            writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
         }
 
         // Dump face centres
@@ -672,7 +672,7 @@ bool Foam::cyclicPolyPatch::matchAnchors
                 )   << "Patch:" << name() << " : "
                     << "Cannot find point on face " << f
                     << " with vertices:"
-                    << IndirectList<point>(pp.points(), f)()
+                    << UIndirectList<point>(pp.points(), f)()
                     << " that matches point " << wantedAnchor
                     << " when matching the halves of cyclic patch " << name()
                     << endl
@@ -1133,8 +1133,8 @@ bool Foam::cyclicPolyPatch::order
     half1ToPatch = half0ToPatch + halfSize;
 
     // Get faces
-    faceList half0Faces(IndirectList<face>(pp, half0ToPatch));
-    faceList half1Faces(IndirectList<face>(pp, half1ToPatch));
+    faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
+    faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
 
     // Get geometric quantities
     pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
@@ -1221,8 +1221,8 @@ bool Foam::cyclicPolyPatch::order
         }
 
         // And redo all matching
-        half0Faces = IndirectList<face>(pp, half0ToPatch);
-        half1Faces = IndirectList<face>(pp, half1ToPatch);
+        half0Faces = UIndirectList<face>(pp, half0ToPatch);
+        half1Faces = UIndirectList<face>(pp, half1ToPatch);
 
         getCentresAndAnchors
         (
@@ -1334,8 +1334,8 @@ bool Foam::cyclicPolyPatch::order
         if (baffleI == halfSize)
         {
             // And redo all matching
-            half0Faces = IndirectList<face>(pp, half0ToPatch);
-            half1Faces = IndirectList<face>(pp, half1ToPatch);
+            half0Faces = UIndirectList<face>(pp, half0ToPatch);
+            half1Faces = UIndirectList<face>(pp, half1ToPatch);
 
             getCentresAndAnchors
             (
@@ -1420,8 +1420,8 @@ bool Foam::cyclicPolyPatch::order
         }
 
         // And redo all matching
-        half0Faces = IndirectList<face>(pp, half0ToPatch);
-        half1Faces = IndirectList<face>(pp, half1ToPatch);
+        half0Faces = UIndirectList<face>(pp, half0ToPatch);
+        half1Faces = UIndirectList<face>(pp, half1ToPatch);
 
         getCentresAndAnchors
         (

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C

@@ -207,7 +207,7 @@ void Foam::processorPolyPatch::calcGeometry()
                     << endl
                     << "Mesh face:" << start()+facei
                     << " vertices:"
-                    << IndirectList<point>(points(), operator[](facei))()
+                    << UIndirectList<point>(points(), operator[](facei))()
                     << endl
                     << "Rerun with processor debug flag set for"
                     << " more information." << exit(FatalError);
@@ -731,7 +731,7 @@ bool Foam::processorPolyPatch::order
                     << " : "
                     << "Cannot find point on face " << pp[oldFaceI]
                     << " with vertices "
-                    << IndirectList<point>(pp.points(), pp[oldFaceI])()
+                    << UIndirectList<point>(pp.points(), pp[oldFaceI])()
                     << " that matches point " << wantedAnchor
                     << " when matching the halves of processor patch " << name()
                     << "Continuing with incorrect face ordering from now on!"

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C

@@ -87,8 +87,7 @@ Foam::PackedBoolList Foam::syncTools::getMasterPoints(const polyMesh& mesh)
 
     labelList minProc(mesh.globalData().nGlobalPoints(), labelMax);
 
-    IndirectList<label>(minProc, sharedPointAddr) =
-        Pstream::myProcNo();
+    UIndirectList<label>(minProc, sharedPointAddr) = Pstream::myProcNo();
 
     Pstream::listCombineGather(minProc, minEqOp<label>());
     Pstream::listCombineScatter(minProc);
@@ -207,8 +206,7 @@ Foam::PackedBoolList Foam::syncTools::getMasterEdges(const polyMesh& mesh)
 
     labelList minProc(mesh.globalData().nGlobalEdges(), labelMax);
 
-    IndirectList<label>(minProc, sharedEdgeAddr) =
-        Pstream::myProcNo();
+    UIndirectList<label>(minProc, sharedEdgeAddr) = Pstream::myProcNo();
 
     Pstream::listCombineGather(minProc, minEqOp<label>());
     Pstream::listCombineScatter(minProc);

src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H

@@ -37,7 +37,6 @@ SourceFiles
 #define ZoneMesh_H
 
 #include "List.H"
-#include "IndirectList.H"
 #include "regIOobject.H"
 #include "HashSet.H"
 #include "pointFieldsFwd.H"

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsEdgeOwner.C

@@ -81,7 +81,7 @@ Foam::PatchTools::edgeOwner
                     << "Edge " << edgeI << " vertices:" << edges[edgeI]
                     << " is used by faces " << nbrFaces
                     << " vertices:"
-                    << IndirectList<Face>(localFaces, nbrFaces)()
+                    << UIndirectList<Face>(localFaces, nbrFaces)()
                     << " none of which use the edge vertices in the same order"
                     << nl << "I give up" << abort(FatalError);
             }

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortEdges.C

@@ -108,7 +108,7 @@ Foam::PatchTools::sortedEdgeFaces
 
             faceAngles.sort();
 
-            sortedEdgeFaces[edgeI] = IndirectList<label>
+            sortedEdgeFaces[edgeI] = UIndirectList<label>
             (
                 faceNbs,
                 faceAngles.indices()

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -1238,10 +1238,7 @@ void Foam::autoSnapDriver::smoothDisplacement
         magDisp().write();
 
         Pout<< "Writing actual patch displacement ..." << endl;
-        vectorField actualPatchDisp
-        (
-            IndirectList<point>(disp, pp.meshPoints())()
-        );
+        vectorField actualPatchDisp(disp, pp.meshPoints());
         dumpMove
         (
             mesh.time().path()/"actualPatchDisplacement.obj",

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -922,7 +922,7 @@ void Foam::refinementSurfaces::findNearest
     List<pointIndexHit>& hitInfo
 ) const
 {
-    labelList geometries(IndirectList<label>(surfaces_, surfacesToTest));
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
 
     // Do the tests. Note that findNearest returns index in geometries.
     searchableSurfacesQueries::findNearest
@@ -955,7 +955,7 @@ void Foam::refinementSurfaces::findNearestRegion
     labelList& hitRegion
 ) const
 {
-    labelList geometries(IndirectList<label>(surfaces_, surfacesToTest));
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
 
     // Do the tests. Note that findNearest returns index in geometries.
     List<pointIndexHit> hitInfo;
@@ -991,7 +991,7 @@ void Foam::refinementSurfaces::findNearestRegion
 
         List<pointIndexHit> localHits
         (
-            IndirectList<pointIndexHit>
+            UIndirectList<pointIndexHit>
             (
                 hitInfo,
                 localIndices

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

@@ -718,7 +718,7 @@ void Foam::fvMeshDistribute::getNeighbourData
             // Which processor they will end up on
             const labelList newProc
             (
-                IndirectList<label>(distribution, pp.faceCells())
+                UIndirectList<label>(distribution, pp.faceCells())
             );
 
             OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
@@ -1192,7 +1192,7 @@ void Foam::fvMeshDistribute::sendMesh
     //
     //    forAll(cellZones, zoneI)
     //    {
-    //        IndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
+    //        UIndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
     //    }
     //}
 

src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C

@@ -269,7 +269,7 @@ bool Foam::hexCellLooper::cut
         {
             FatalErrorIn("hexCellLooper::walkHex") << "Face:" << faceVerts
                 << " on points:" << facePoints << endl
-                << IndirectList<point>(facePoints, faceVerts)()
+                << UIndirectList<point>(facePoints, faceVerts)()
                 << abort(FatalError);
         }
     }

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -103,8 +103,8 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
 )
 {
     // Clear the fields for accumulation
-    IndirectList<vector>(cellCentres_, changedCells) = vector::zero;
-    IndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
+    UIndirectList<vector>(cellCentres_, changedCells) = vector::zero;
+    UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
 
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
@@ -112,9 +112,9 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
     // first estimate the approximate cell centre as the average of face centres
 
     vectorField cEst(mesh_.nCells());
-    IndirectList<vector>(cEst, changedCells) = vector::zero;
+    UIndirectList<vector>(cEst, changedCells) = vector::zero;
     scalarField nCellFaces(mesh_.nCells());
-    IndirectList<scalar>(nCellFaces, changedCells) = 0.0;
+    UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
 
     forAll(changedFaces, i)
     {

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

@@ -32,7 +32,6 @@ License
 #include "treeDataFace.H"
 #include "indexedOctree.H"
 #include "OFstream.H"
-#include "IndirectList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -153,11 +152,7 @@ void Foam::faceCoupleInfo::writePointsFaces() const
         (
             "cutToMasterPoints.obj",
             m.localPoints(),
-            pointField
-            (
-                IndirectList<point>(c.localPoints(), masterToCutPoints_)()
-            )
-        );
+            pointField(c.localPoints(), masterToCutPoints_));
     }
     {
         Pout<< "Writing cutToSlavePoints to cutToSlavePoints.obj" << endl;
@@ -166,10 +161,7 @@ void Foam::faceCoupleInfo::writePointsFaces() const
         (
             "cutToSlavePoints.obj",
             s.localPoints(),
-            pointField
-            (
-                IndirectList<point>(c.localPoints(), slaveToCutPoints_)()
-            )
+            pointField(c.localPoints(), slaveToCutPoints_)
         );
     }
 
@@ -405,7 +397,7 @@ Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
     if (report)
     {
         Pout<< "mostAlignedEdge : finding nearest edge among "
-            << IndirectList<edge>(cutFaces().edges(), pEdges)()
+            << UIndirectList<edge>(cutFaces().edges(), pEdges)()
             << " connected to point " << pointI
             << " coord:" << localPoints[pointI]
             << " running between " << edgeStart << " coord:"
@@ -623,7 +615,7 @@ void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
                         "(const labelList&)"
                     )   << " unsplitEdge:" << unsplitEdge
                         << " does not correspond to split edges "
-                        << IndirectList<edge>(cutEdges, stringedEdges)()
+                        << UIndirectList<edge>(cutEdges, stringedEdges)()
                         << abort(FatalError);
                 }
             }
@@ -631,7 +623,7 @@ void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
             //Pout<< "For master edge:"
             //    << unsplitEdge
             //    << " Found stringed points "
-            //    <<  IndirectList<point>
+            //    <<  UIndirectList<point>
             //        (
             //            cutFaces().localPoints(),
             //            splitPoints.shrink()
@@ -664,9 +656,9 @@ Foam::label Foam::faceCoupleInfo::matchFaces
             "(const scalar, const face&, const pointField&"
             ", const face&, const pointField&)"
         )   << "Different sizes for supposedly matching faces." << nl
-            << "f0:" << f0 << " coords:" << IndirectList<point>(points0, f0)()
+            << "f0:" << f0 << " coords:" << UIndirectList<point>(points0, f0)()
             << nl
-            << "f1:" << f1 << " coords:" << IndirectList<point>(points1, f1)()
+            << "f1:" << f1 << " coords:" << UIndirectList<point>(points1, f1)()
             << abort(FatalError);
     }
 
@@ -721,10 +713,9 @@ Foam::label Foam::faceCoupleInfo::matchFaces
             ", const face&, const pointField&)"
         )   << "No unique match between two faces" << nl
             << "Face " << f0 << " coords "
-            << IndirectList<point>(points0, f0)()
-            << nl
+            << UIndirectList<point>(points0, f0)() << nl
             << "Face " << f1 << " coords "
-            << IndirectList<point>(points1, f1)()
+            << UIndirectList<point>(points1, f1)()
             << "when using tolerance " << absTol
             << " and forwardMatching:" << sameOrientation
             << abort(FatalError);
@@ -1575,7 +1566,7 @@ void Foam::faceCoupleInfo::perfectPointMatch
 
 
     // Use compaction lists to renumber cutPoints.
-    cutPoints_ = IndirectList<point>(cutPoints_, compactToCut)();
+    cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
     {
         const faceList& cutLocalFaces = cutFaces().localFaces();
 
@@ -1770,11 +1761,11 @@ void Foam::faceCoupleInfo::subDivisionMatch
                 writeOBJ
                 (
                     "errorEdges.obj",
-                    IndirectList<edge>
+                    UIndirectList<edge>
                     (
                         cutFaces().edges(),
                         cutFaces().pointEdges()[cutPointI]
-                    )(),
+                    ),
                     cutFaces().localPoints(),
                     false
                 );
@@ -1894,7 +1885,7 @@ void Foam::faceCoupleInfo::subDivisionMatch
                 "(const polyMesh&, const bool, const scalar)"
             )   << "Did not match all of cutFaces to a master face" << nl
                 << "First unmatched cut face:" << cutFaceI << " with points:"
-                << IndirectList<point>(cutFaces().points(), cutF)()
+                << UIndirectList<point>(cutFaces().points(), cutF)()
                 << nl
                 << "This usually means that the slave patch is not a"
                 << " subdivision of the master patch"

src/dynamicMesh/polyMeshAdder/polyMeshAdder.C

@@ -1929,11 +1929,8 @@ Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
     //(
     //    pointField
     //    (
-    //        IndirectList<point>
-    //        (
     //        mesh.points(),
     //        sharedPointLabels
-    //        )()
     //    ),
     //    mergeDist,
     //    false,

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -110,7 +110,7 @@ Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces
     );
 
     // Extract pp part
-    return IndirectList<labelList>(globalEdgeFaces, meshEdges)();
+    return UIndirectList<labelList>(globalEdgeFaces, meshEdges);
 }
 
 
@@ -629,7 +629,7 @@ void Foam::addPatchCellLayer::setRefinement
 
         {
             labelList n(mesh_.nPoints(), 0);
-            IndirectList<label>(n, meshPoints) = nPointLayers;
+            UIndirectList<label>(n, meshPoints) = nPointLayers;
             syncTools::syncPointList(mesh_, n, maxEqOp<label>(), 0, false);
 
             // Non-synced
@@ -706,7 +706,7 @@ void Foam::addPatchCellLayer::setRefinement
 
         {
             pointField d(mesh_.nPoints(), wallPoint::greatPoint);
-            IndirectList<point>(d, meshPoints) = firstLayerDisp;
+            UIndirectList<point>(d, meshPoints) = firstLayerDisp;
             syncTools::syncPointList
             (
                 mesh_,

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -565,7 +565,7 @@ void Foam::combineFaces::setRefinement
             const labelList& setFaces = faceSets[setI];
 
             masterFace_[setI] = setFaces[0];
-            faceSetsVertices_[setI] = IndirectList<face>
+            faceSetsVertices_[setI] = UIndirectList<face>
             (
                 mesh_.faces(),
                 setFaces

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -28,7 +28,6 @@ License
 #include "polyMesh.H"
 #include "polyTopoChange.H"
 #include "ListOps.H"
-#include "IndirectList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -746,7 +746,7 @@ Foam::label Foam::hexRef8::findLevel
 
             FatalErrorIn("hexRef8::findLevel(..)")
                 << "face:" << f
-                << " level:" << IndirectList<label>(pointLevel_, f)()
+                << " level:" << UIndirectList<label>(pointLevel_, f)()
                 << " startFp:" << startFp
                 << " wantedLevel:" << wantedLevel
                 << abort(FatalError);
@@ -774,7 +774,7 @@ Foam::label Foam::hexRef8::findLevel
 
     FatalErrorIn("hexRef8::findLevel(..)")
         << "face:" << f
-        << " level:" << IndirectList<label>(pointLevel_, f)()
+        << " level:" << UIndirectList<label>(pointLevel_, f)()
         << " startFp:" << startFp
         << " wantedLevel:" << wantedLevel
         << abort(FatalError);
@@ -823,7 +823,7 @@ void Foam::hexRef8::checkInternalOrientation
 )
 {
     face compactFace(identity(newFace.size()));
-    pointField compactPoints(IndirectList<point>(meshMod.points(), newFace)());
+    pointField compactPoints(meshMod.points(), newFace);
 
     vector n(compactFace.normal(compactPoints));
 
@@ -869,7 +869,7 @@ void Foam::hexRef8::checkBoundaryOrientation
 )
 {
     face compactFace(identity(newFace.size()));
-    pointField compactPoints(IndirectList<point>(meshMod.points(), newFace)());
+    pointField compactPoints(meshMod.points(), newFace);
 
     vector n(compactFace.normal(compactPoints));
 
@@ -1308,11 +1308,11 @@ void Foam::hexRef8::createInternalFaces
                             << "cell:" << cellI << " cLevel:" << cLevel
                             << " cell points:" << cPoints
                             << " pointLevel:"
-                            << IndirectList<label>(pointLevel_, cPoints)()
+                            << UIndirectList<label>(pointLevel_, cPoints)()
                             << " face:" << faceI
                             << " f:" << f
                             << " pointLevel:"
-                            << IndirectList<label>(pointLevel_, f)()
+                            << UIndirectList<label>(pointLevel_, f)()
                             << " faceAnchorLevel:" << faceAnchorLevel[faceI]
                             << " faceMidPoint:" << faceMidPoint[faceI]
                             << " faceMidPointI:" << faceMidPointI
@@ -1381,11 +1381,11 @@ void Foam::hexRef8::createInternalFaces
                             << "cell:" << cellI << " cLevel:" << cLevel
                             << " cell points:" << cPoints
                             << " pointLevel:"
-                            << IndirectList<label>(pointLevel_, cPoints)()
+                            << UIndirectList<label>(pointLevel_, cPoints)()
                             << " face:" << faceI
                             << " f:" << f
                             << " pointLevel:"
-                            << IndirectList<label>(pointLevel_, f)()
+                            << UIndirectList<label>(pointLevel_, f)()
                             << " faceAnchorLevel:" << faceAnchorLevel[faceI]
                             << " faceMidPoint:" << faceMidPoint[faceI]
                             << " faceMidPointI:" << faceMidPointI
@@ -3511,7 +3511,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
                         << " lower level" << endl
                         << "cellPoints:" << cPoints << endl
                         << "pointLevels:"
-                        << IndirectList<label>(pointLevel_, cPoints)() << endl
+                        << UIndirectList<label>(pointLevel_, cPoints)() << endl
                         << abort(FatalError);
                 }
             }
@@ -4473,7 +4473,7 @@ void Foam::hexRef8::checkMesh() const
                     << "Coupled face " << faceI
                     << " on patch " << patchI
                     << " " << mesh_.boundaryMesh()[patchI].name()
-                    << " coords:" << IndirectList<point>(mesh_.points(), f)()
+                    << " coords:" << UIndirectList<point>(mesh_.points(), f)()
                     << " has face area:" << magArea
                     << " (coupled) neighbour face area differs:"
                     << neiFaceAreas[i]
@@ -4515,7 +4515,7 @@ void Foam::hexRef8::checkMesh() const
                     << "Coupled face " << faceI
                     << " on patch " << patchI
                     << " " << mesh_.boundaryMesh()[patchI].name()
-                    << " coords:" << IndirectList<point>(mesh_.points(), f)()
+                    << " coords:" << UIndirectList<point>(mesh_.points(), f)()
                     << " has size:" << f.size()
                     << " (coupled) neighbour face has size:"
                     << nVerts[i]
@@ -4565,7 +4565,7 @@ void Foam::hexRef8::checkMesh() const
                     << "Coupled face " << faceI
                     << " on patch " << patchI
                     << " " << mesh_.boundaryMesh()[patchI].name()
-                    << " coords:" << IndirectList<point>(mesh_.points(), f)()
+                    << " coords:" << UIndirectList<point>(mesh_.points(), f)()
                     << " has anchor vector:" << anchorVec
                     << " (coupled) neighbour face anchor vector differs:"
                     << anchorPoints[i]

src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C

@@ -375,7 +375,7 @@ void Foam::removeFaces::mergeFaces
     //{
     //    Pout<< "Modifying masterface " << faceI
     //        << " from faces:" << faceLabels
-    //        << " old verts:" << IndirectList<face>(mesh_.faces(), faceLabels)
+    //        << " old verts:" << UIndirectList<face>(mesh_.faces(), faceLabels)
     //        << " for new verts:"
     //        << mergedFace
     //        << " possibly new owner " << own

src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C

@@ -353,7 +353,7 @@ void Foam::removePoints::setRefinement
                 // Points from the mesh
                 List<point> meshPoints
                 (
-                    IndirectList<point>
+                    UIndirectList<point>
                     (
                         mesh_.points(),
                         mesh_.faces()[savedFaceLabels_[saveI]]  // mesh face

src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C

@@ -202,7 +202,7 @@ void Foam::enrichedPatch::calcCutFaces() const
                 // Grab the next point options
 //                 Pout << "curPointLabel: " << mp[curPointLabel] << endl;
                 const labelList& nextPoints = pp[curPointLabel];
-//                 Pout << "nextPoints: " << IndirectList<label>(mp, nextPoints) << endl;
+//                 Pout << "nextPoints: " << UIndirectList<label>(mp, nextPoints) << endl;
                 // Get the vector along the edge and the right vector
                 vector ahead = curPoint - lp[prevPointLabel];
                 ahead -= normal*(normal & ahead);

src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C

@@ -165,7 +165,7 @@ void directMappedFixedValueFvPatchField<Type>::getNewValues
         if (Pstream::myProcNo() == sendProc)
         {
             OPstream toProc(Pstream::scheduled, recvProc);
-            toProc<< IndirectList<Type>(sendValues, sendLabels[recvProc])();
+            toProc<< UIndirectList<Type>(sendValues, sendLabels[recvProc])();
         }
         else
         {
@@ -188,7 +188,11 @@ void directMappedFixedValueFvPatchField<Type>::getNewValues
 
     // Do data from myself
     {
-        IndirectList<Type> fromFld(sendValues, sendLabels[Pstream::myProcNo()]);
+        UIndirectList<Type> fromFld
+        (
+            sendValues,
+            sendLabels[Pstream::myProcNo()]
+        );
 
         // Destination faces
         const labelList& faceLabels = receiveFaceLabels[Pstream::myProcNo()];

src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C

@@ -162,8 +162,8 @@ void directMappedVelocityFluxFixedValueFvPatchField::getNewValues
         if (Pstream::myProcNo() == sendProc)
         {
             OPstream toProc(Pstream::scheduled, recvProc);
-            toProc<< IndirectList<vector>(sendUValues, sendLabels[recvProc])();
-            toProc<< IndirectList<scalar>
+            toProc<< UIndirectList<vector>(sendUValues, sendLabels[recvProc])();
+            toProc<< UIndirectList<scalar>
                 (
                     sendPhiValues,
                     sendLabels[recvProc]
@@ -192,13 +192,13 @@ void directMappedVelocityFluxFixedValueFvPatchField::getNewValues
 
     // Do data from myself
     {
-        IndirectList<vector> fromUFld
+        UIndirectList<vector> fromUFld
             (
                 sendUValues,
                 sendLabels[Pstream::myProcNo()]
             );
 
-        IndirectList<scalar> fromPhiFld
+        UIndirectList<scalar> fromPhiFld
             (
                 sendPhiValues,
                 sendLabels[Pstream::myProcNo()]

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C

@@ -29,14 +29,11 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-int fvSchemes::debug(Foam::debug::debugSwitch("fvSchemes", false));
+int Foam::fvSchemes::debug(Foam::debug::debugSwitch("fvSchemes", false));
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-fvSchemes::fvSchemes(const objectRegistry& obr)
+Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
 :
     IOdictionary
     (
@@ -124,7 +121,12 @@ fvSchemes::fvSchemes(const objectRegistry& obr)
         ITstream(objectPath() + "::fluxRequired",
         tokenList())()
     ),
-    defaultFluxRequired_(false)
+    defaultFluxRequired_(false),
+    cacheFields_
+    (
+        ITstream(objectPath() + "::cacheFields",
+        tokenList())()
+    )
 {
     read();
 }
@@ -132,7 +134,7 @@ fvSchemes::fvSchemes(const objectRegistry& obr)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool fvSchemes::read()
+bool Foam::fvSchemes::read()
 {
     if (regIOobject::read())
     {
@@ -329,6 +331,11 @@ bool fvSchemes::read()
             }
         }
 
+        if (dict.found("cacheFields"))
+        {
+            cacheFields_ = dict.subDict("cacheFields");
+        }
+
         return true;
     }
     else
@@ -338,7 +345,7 @@ bool fvSchemes::read()
 }
 
 
-const dictionary& fvSchemes::schemesDict() const
+const Foam::dictionary& Foam::fvSchemes::schemesDict() const
 {
     if (found("select"))
     {
@@ -351,7 +358,7 @@ const dictionary& fvSchemes::schemesDict() const
 }
 
 
-ITstream& fvSchemes::ddtScheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::ddtScheme(const word& name) const
 {
     if (debug)
     {
@@ -370,7 +377,7 @@ ITstream& fvSchemes::ddtScheme(const word& name) const
 }
 
 
-ITstream& fvSchemes::d2dt2Scheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::d2dt2Scheme(const word& name) const
 {
     if (debug)
     {
@@ -389,7 +396,7 @@ ITstream& fvSchemes::d2dt2Scheme(const word& name) const
 }
 
 
-ITstream& fvSchemes::interpolationScheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::interpolationScheme(const word& name) const
 {
     if (debug)
     {
@@ -412,7 +419,7 @@ ITstream& fvSchemes::interpolationScheme(const word& name) const
 }
 
 
-ITstream& fvSchemes::divScheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::divScheme(const word& name) const
 {
     if (debug)
     {
@@ -431,7 +438,7 @@ ITstream& fvSchemes::divScheme(const word& name) const
 }
 
 
-ITstream& fvSchemes::gradScheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::gradScheme(const word& name) const
 {
     if (debug)
     {
@@ -450,7 +457,7 @@ ITstream& fvSchemes::gradScheme(const word& name) const
 }
 
 
-ITstream& fvSchemes::snGradScheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::snGradScheme(const word& name) const
 {
     if (debug)
     {
@@ -469,7 +476,7 @@ ITstream& fvSchemes::snGradScheme(const word& name) const
 }
 
 
-ITstream& fvSchemes::laplacianScheme(const word& name) const
+Foam::ITstream& Foam::fvSchemes::laplacianScheme(const word& name) const
 {
     if (debug)
     {
@@ -488,7 +495,7 @@ ITstream& fvSchemes::laplacianScheme(const word& name) const
 }
 
 
-bool fvSchemes::fluxRequired(const word& name) const
+bool Foam::fvSchemes::fluxRequired(const word& name) const
 {
     if (debug)
     {
@@ -506,8 +513,22 @@ bool fvSchemes::fluxRequired(const word& name) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+bool Foam::fvSchemes::cache(const word& name) const
+{
+    if (debug)
+    {
+        Info<< "Lookup cache for " << name << endl;
+    }
+
+    if (cacheFields_.found(name))
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H

@@ -79,6 +79,8 @@ class fvSchemes
         dictionary fluxRequired_;
         bool defaultFluxRequired_;
 
+        dictionary cacheFields_;
+
 
     // Private Member Functions
 
@@ -121,6 +123,8 @@ public:
 
             bool fluxRequired(const word& name) const;
 
+            bool cache(const word& name) const;
+
 
         // Read
 

src/lagrangian/Allwmake

@@ -6,6 +6,7 @@ wmake libso basic
 wmake libso solidParticle
 wmake libso intermediate
 wmake libso dieselSpray
+wmake libso dsmc
 
 molecularDynamics/Allwmake
 

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

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

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();
+}
+
+
+// ************************************************************************* //

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

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();
+}
+
+
+// ************************************************************************* //

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()
+{}
+
+
+// ************************************************************************* //

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

src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C

@@ -22,56 +22,54 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dsmcCloud, 0);
+};
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//- Construct given size
-template<class T>
-inline Foam::IndirectList<T>::IndirectList
+Foam::dsmcCloud::dsmcCloud
 (
-    const Foam::UList<T>& completeList,
-    const Foam::List<label>& addresses
+    const word& cloudType,
+    const volScalarField& T,
+    const volVectorField& U
 )
 :
-    completeList_(completeList),
-    addresses_(addresses)
+    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  * * * * * * * * * * * * * //
 
-template<class T>
-inline Foam::label Foam::IndirectList<T>::size() const
+void Foam::dsmcCloud::writeFields() const
 {
-    return addresses_.size();
+    dsmcParcel::writeFields(*this);
 }
 
-
-template<class T>
-inline const Foam::UList<T>& Foam::IndirectList<T>::
-completeList() const
-{
-    return completeList_;
-}
-
-
-template<class T>
-inline const Foam::List<Foam::label>& Foam::IndirectList<T>::addresses() const
-{
-    return addresses_;
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-template<class T>
-inline const T& Foam::IndirectList<T>::operator[](const Foam::label i) const
-{
-    return completeList_[addresses_[i]];
-}
-
-
 // ************************************************************************* //

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

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

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

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_;
+}
+
+
+// ************************************************************************* //

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

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);
+};
+
+
+// ************************************************************************* //

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()
+{}
+
+
+// ************************************************************************* //

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

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
+    );
+};
+
+
+// ************************************************************************* //

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
+    );
+};
+
+
+// ************************************************************************* //

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
+    );
+};
+
+
+// ************************************************************************* //

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

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

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)
+    );
+}
+
+
+// ************************************************************************* //

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);
+}
+
+
+// ************************************************************************* //

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

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);
+}
+
+
+// ************************************************************************* //

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

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

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

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

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

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));
+}
+
+
+// ************************************************************************* //

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()
+{}
+
+
+// ************************************************************************* //

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

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);
+}
+
+
+// ************************************************************************* //

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

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

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

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));
+}
+
+
+// ************************************************************************* //