diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
index 1ea52cc6aa..cb068960e8 100644
--- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
+++ b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index ec1d09a00f..6c8cbe4aa1 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -390,7 +390,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     // slightly towards the cell-centre.
     if (trackFraction < SMALL)
     {
-        position_ += 1.0e-6*(mesh.cellCentres()[celli_] - position_);
+        position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
     }
 
     return trackFraction;
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 4c4bcd0a43..53177e9afe 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -37,6 +37,7 @@ 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>
@@ -113,87 +114,161 @@ void Foam::DsmcCloud<ParcelType>::initialise
 
     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)
+    forAll(mesh_.cells(), cell)
     {
-        const word& moleculeName(molecules[i]);
+        const vector& cC = mesh_.cellCentres()[cell];
+        const labelList& cellFaces = mesh_.cells()[cell];
+        const scalar cV = mesh_.cellVolumes()[cell];
 
-        label typeId(findIndex(typeIdList_, moleculeName));
+        label nTets = 0;
 
-        if (typeId == -1)
+        // Each face is split into nEdges (or nVertices) - 2 tets.
+        forAll(cellFaces, face)
         {
-            FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
-                << "typeId " << moleculeName << "not defined." << nl
-                << abort(FatalError);
+            nTets += mesh_.faces()[cellFaces[face]].size() - 2;
         }
 
-        const typename ParcelType::constantProperties& cP = constProps(typeId);
+        // Calculate the cumulative tet volumes circulating around the cell and
+        // record the vertex labels of each.
+        scalarList cTetVFracs(nTets, 0.0);
 
-        scalar numberDensity = numberDensities[i];
+        List<labelList> tetPtIs(nTets, labelList(3,-1));
 
-        scalar spacing = pow(numberDensity,-(1.0/3.0));
+        // Keep track of which tet this is.
+        label tet = 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++)
+        forAll(cellFaces, face)
         {
-            for (int j = 0; j < nj; j++)
+            const labelList& facePoints = mesh_.faces()[cellFaces[face]];
+
+            label pointI = 1;
+            while ((pointI + 1) < facePoints.size())
             {
-                for (int k = 0; k < nk; k++)
+
+                const vector& pA = mesh_.points()[facePoints[0]];
+                const vector& pB = mesh_.points()[facePoints[pointI]];
+                const vector& pC = mesh_.points()[facePoints[pointI + 1]];
+
+                cTetVFracs[tet] =
+                    mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
+                  + cTetVFracs[max((tet - 1),0)];
+
+                tetPtIs[tet][0] = facePoints[0];
+                tetPtIs[tet][1] = facePoints[pointI];
+                tetPtIs[tet][2] = facePoints[pointI + 1];
+
+                pointI++;
+                tet++;
+            }
+        }
+
+        // Force the last volume fraction value to 1.0 to avoid any
+        // rounding/non-flat face errors giving a value < 1.0
+        cTetVFracs[nTets - 1] = 1.0;
+
+        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];
+
+            // Calculate the number of particles required
+            scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
+
+            // Only integer numbers of particles can be inserted
+            label nParticlesToInsert = label(particlesRequired);
+
+            // Add another particle with a probability proportional to the
+            // remainder of taking the integer part of particlesRequired
+            if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
+            {
+                nParticlesToInsert++;
+            }
+
+            for (label pI = 0; pI < nParticlesToInsert; pI++)
+            {
+                // Choose a random point in a generic tetrahedron
+
+                scalar s = rndGen_.scalar01();
+                scalar t = rndGen_.scalar01();
+                scalar u = rndGen_.scalar01();
+
+                if (s + t > 1.0)
                 {
-                    point p
-                    (
-                        x0 + (i + 0.5)*delta.x(),
-                        y0 + (j + 0.5)*delta.y(),
-                        z0 + (k + 0.5)*delta.z()
-                    );
+                    s = 1.0 - s;
+                    t = 1.0 - t;
+                }
 
-                    p.x() += pert*(rndGen_.scalar01() - 0.5);
-                    p.y() += pert*(rndGen_.scalar01() - 0.5);
-                    p.z() += pert*(rndGen_.scalar01() - 0.5);
+                if (t + u > 1.0)
+                {
+                    scalar tmp = u;
+                    u = 1.0 - s - t;
+                    t = 1.0 - tmp;
+                }
+                else if (s + t + u > 1.0)
+                {
+                    scalar tmp = u;
+                    u = s + t + u - 1.0;
+                    s = 1.0 - t - tmp;
+                }
 
-                    label cell = mesh_.findCell(p);
+                // Choose a tetrahedron to insert in, based on their relative
+                // volumes
+                scalar tetSelection = rndGen_.scalar01();
 
-                    vector U = equipartitionLinearVelocity
-                    (
-                        temperature,
-                        cP.mass()
-                    );
+                // Selected tetrahedron
+                label sTet = -1;
 
-                    scalar Ei = equipartitionInternalEnergy
-                    (
-                        temperature,
-                        cP.internalDegreesOfFreedom()
-                    );
+                forAll(cTetVFracs, tet)
+                {
+                    sTet = tet;
 
-                    U += velocity;
-
-                    if (cell >= 0)
+                    if (cTetVFracs[tet] >= tetSelection)
                     {
-                        addNewParcel
-                        (
-                            p,
-                            U,
-                            Ei,
-                            cell,
-                            typeId
-                        );
+                        break;
                     }
                 }
+
+                vector p =
+                    (1 - s - t - u)*cC
+                  + s*mesh_.points()[tetPtIs[sTet][0]]
+                  + t*mesh_.points()[tetPtIs[sTet][1]]
+                  + u*mesh_.points()[tetPtIs[sTet][2]];
+
+                vector U = equipartitionLinearVelocity
+                (
+                    temperature,
+                    cP.mass()
+                );
+
+                scalar Ei = equipartitionInternalEnergy
+                (
+                    temperature,
+                    cP.internalDegreesOfFreedom()
+                );
+
+                U += velocity;
+
+                addNewParcel
+                (
+                    p,
+                    U,
+                    Ei,
+                    cell,
+                    typeId
+                );
             }
         }
     }
@@ -276,7 +351,7 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
             scalar selectedPairs =
                 collisionSelectionRemainder_[celli]
-              + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT
+              + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
                /mesh_.cellVolumes()[celli];
 
             label nCandidates(selectedPairs);
@@ -551,6 +626,13 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     buildConstProps();
 
     buildCellOccupancy();
+
+    // Initialise the collision selection remainder to a random value between 0
+    // and 1.
+    forAll(collisionSelectionRemainder_, i)
+    {
+        collisionSelectionRemainder_[i] = rndGen_.scalar01();
+    }
 }
 
 
@@ -739,22 +821,27 @@ void Foam::DsmcCloud<ParcelType>::info() const
 
     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
+        << nDsmcParticles
         << endl;
+
+    if (nDsmcParticles)
+    {
+        Info<< "    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
+            << endl;
+    }
 }
 
 
@@ -785,7 +872,11 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
 {
     scalar Ei = 0.0;
 
-    if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
+    if (iDof < SMALL)
+    {
+        return Ei;
+    }
+    else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
     {
         // Special case for iDof = 2, i.e. diatomics;
         Ei = -log(rndGen_.scalar01())*kb*temperature;
@@ -796,7 +887,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
 
         scalar energyRatio;
 
-        scalar P;
+        scalar P = -1;
 
         do
         {
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index b246b04990..20bed1ee7b 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -299,6 +299,12 @@ public:
                     scalar mass
                 ) const;
 
+                inline scalarField maxwellianAverageSpeed
+                (
+                    scalarField temperature,
+                    scalar mass
+                ) const;
+
                 //- RMS particle speed
                 inline scalar maxwellianRMSSpeed
                 (
@@ -306,6 +312,12 @@ public:
                     scalar mass
                 ) const;
 
+                inline scalarField maxwellianRMSSpeed
+                (
+                    scalarField temperature,
+                    scalar mass
+                ) const;
+
                 //- Most probable speed
                 inline scalar maxwellianMostProbableSpeed
                 (
@@ -313,6 +325,11 @@ public:
                     scalar mass
                 ) const;
 
+                inline scalarField maxwellianMostProbableSpeed
+                (
+                    scalarField temperature,
+                    scalar mass
+                ) const;
 
             // Sub-models
 
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 73899636e8..d39de4954f 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -302,6 +302,17 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
+(
+    scalarField 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
 (
@@ -313,6 +324,17 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
+(
+    scalarField temperature,
+    scalar mass
+) const
+{
+    return sqrt(3.0*kb*temperature/mass);
+}
+
+
 template<class ParcelType>
 inline Foam::scalar
 Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
@@ -325,6 +347,18 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField
+Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
+(
+    scalarField 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
@@ -343,7 +377,7 @@ Foam::DsmcCloud<ParcelType>::rhoN() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
         )
     );
 
@@ -380,7 +414,7 @@ Foam::DsmcCloud<ParcelType>::rhoM() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL)
         )
     );
 
@@ -568,7 +602,7 @@ Foam::DsmcCloud<ParcelType>::iDof() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
         )
     );
 
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index f4203f4d66..cd18d9064d 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -36,32 +36,27 @@ Foam::FreeStream<CloudType>::FreeStream
 )
 :
     InflowBoundaryModel<CloudType>(dict, cloud, typeName),
-    patchIndex_(),
-    temperature_(readScalar(this->coeffDict().lookup("temperature"))),
-    velocity_(this->coeffDict().lookup("velocity")),
+    patches_(),
     moleculeTypeIds_(),
     numberDensities_(),
     particleFluxAccumulators_()
 {
-    word patchName = this->coeffDict().lookup("patch");
+    // Identify which patches to use
 
-    patchIndex_ = cloud.mesh().boundaryMesh().findPatchID(patchName);
+    DynamicList<label> patches;
 
-    const polyPatch& patch = cloud.mesh().boundaryMesh()[patchIndex_];
-
-    if (patchIndex_ == -1)
+    forAll(cloud.mesh().boundaryMesh(), p)
     {
-        FatalErrorIn
-        (
-            "Foam::FreeStream<CloudType>::FreeStream"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )   << "patch " << patchName << " not found." << nl
-            << abort(FatalError);
+        const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
+
+        if (patch.type() == polyPatch::typeName)
+        {
+            patches.append(p);
+        }
     }
 
+    patches_.transfer(patches);
+
     const dictionary& numberDensitiesDict
     (
         this->coeffDict().subDict("numberDensities")
@@ -69,10 +64,24 @@ Foam::FreeStream<CloudType>::FreeStream
 
     List<word> molecules(numberDensitiesDict.toc());
 
-    numberDensities_.setSize(molecules.size());
+    // Initialise the particleFluxAccumulators_
+    particleFluxAccumulators_.setSize(patches_.size());
+
+    forAll(patches_, p)
+    {
+        const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
+
+        particleFluxAccumulators_[p] = List<Field<scalar> >
+        (
+            molecules.size(),
+            Field<scalar>(patch.size(), 0.0)
+        );
+    }
 
     moleculeTypeIds_.setSize(molecules.size());
 
+    numberDensities_.setSize(molecules.size());
+
     forAll(molecules, i)
     {
         numberDensities_[i] = readScalar
@@ -97,12 +106,6 @@ Foam::FreeStream<CloudType>::FreeStream
     }
 
     numberDensities_ /= cloud.nParticle();
-
-    particleFluxAccumulators_.setSize
-    (
-        molecules.size(),
-        Field<scalar>(patch.size(), 0)
-    );
 }
 
 
@@ -127,144 +130,262 @@ void Foam::FreeStream<CloudType>::inflow()
 
     Random& rndGen(cloud.rndGen());
 
-    const polyPatch& patch = mesh.boundaryMesh()[patchIndex_];
+    scalar sqrtPi = sqrt(mathematicalConstant::pi);
 
     label particlesInserted = 0;
 
-    // Add mass to the accumulators.  negative face area dotted with the
-    // velocity to point flux into the domain.
+    const volScalarField::GeometricBoundaryField& boundaryT
+    (
+        cloud.T().boundaryField()
+    );
 
-    forAll(particleFluxAccumulators_, i)
+    const volVectorField::GeometricBoundaryField& boundaryU
+    (
+        cloud.U().boundaryField()
+    );
+
+    forAll(patches_, p)
     {
-        particleFluxAccumulators_[i] +=
-           -patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
-    }
+        label patchI = patches_[p];
 
-    forAll(patch, f)
-    {
-        // Loop over all faces as the outer loop to avoid calculating
-        // geometrical properties multiple times.
+        const polyPatch& patch = mesh.boundaryMesh()[patchI];
 
-        labelList faceVertices = patch[f];
+        // Add mass to the accumulators.  negative face area dotted with the
+        // velocity to point flux into the domain.
 
-        label nVertices = faceVertices.size();
+        // Take a reference to the particleFluxAccumulator for this patch
+        List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
 
-        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++)
+        forAll(pFA, i)
         {
-            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++)
+            scalarField mostProbableSpeed
+            (
+                cloud.maxwellianMostProbableSpeed
+                (
+                    boundaryT[patchI],
+                    mass
+                )
+            );
+
+            // Dotting boundary velocity with the face unit normal (which points
+            // out of the domain, so it must be negated), dividing by the most
+            // probable speed to form molecularSpeedRatio * cosTheta
+
+            scalarField sCosTheta =
+                boundaryU[patchI]
+              & -patch.faceAreas()/mag(patch.faceAreas())
+               /mostProbableSpeed;
+
+            // From Bird eqn 4.22
+
+            pFA[i] +=
+                mag(patch.faceAreas())*numberDensities_[i]*deltaT
+               *mostProbableSpeed
+               *(
+                   exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
+                )
+               /(2.0*sqrtPi);
+        }
+
+        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]);
+
+            // Cumulative triangle area fractions
+            List<scalar> cTriAFracs(nVertices);
+
+            for (label v = 0; v < nVertices - 1; v++)
             {
-                // Choose a triangle to insert on, based on their relative area
+                const point& vA = mesh.points()[faceVertices[v]];
 
-                scalar triSelection = rndGen.scalar01();
+                const point& vB = mesh.points()[faceVertices[(v + 1)]];
 
-                // Selected triangle
-                label sTri = -1;
+                cTriAFracs[v] =
+                    0.5*mag((vA - fC)^(vB - fC))/fA
+                  + cTriAFracs[max((v - 1), 0)];
+            }
 
-                forAll(cTriAFracs, tri)
+            // 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 n = patch.faceAreas()[f];
+            n /= -mag(n);
+
+            // Wall tangential unit vector. Use the direction between the
+            // face centre and the first vertex in the list
+            vector t1 = fC - (mesh.points()[faceVertices[0]]);
+            t1 /= mag(t1);
+
+            // Other tangential unit vector.  Rescaling in case face is not
+            // flat and n and t1 aren't perfectly orthogonal
+            vector t2 = n^t1;
+            t2 /= mag(t2);
+
+            scalar faceTemperature = boundaryT[patchI][f];
+
+            const vector& faceVelocity = boundaryU[patchI][f];
+
+            forAll(pFA, i)
+            {
+                scalar& faceAccumulator = pFA[i][f];
+
+                // Number of whole particles to insert
+                label nI = max(label(faceAccumulator), 0);
+
+                // Add another particle with a probability proportional to the
+                // remainder of taking the integer part of faceAccumulator
+                if ((faceAccumulator - nI) > rndGen.scalar01())
                 {
-                    sTri = tri;
-
-                    if (cTriAFracs[tri] >= triSelection)
-                    {
-                        break;
-                    }
+                    nI++;
                 }
 
-                // 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
+                faceAccumulator -= nI;
 
-                const point& A = fC;
-                const point& B = mesh.points()[faceVertices[sTri]];
-                const point& C =
+                label typeId = moleculeTypeIds_[i];
+
+                scalar mass = cloud.constProps(typeId).mass();
+
+                for (label i = 0; i < nI; i++)
+                {
+                    // 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());
+                    scalar s = rndGen.scalar01();
+                    scalar t = sqrt(rndGen.scalar01());
 
-                point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
+                    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
+                    // Velocity generation
+
+                    scalar mostProbableSpeed
+                    (
+                        cloud.maxwellianMostProbableSpeed
+                        (
+                            faceTemperature,
+                            mass
+                        )
                     );
 
-                U += velocity_;
+                    scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
 
-                scalar Ei = cloud.equipartitionInternalEnergy
-                (
-                    temperature_,
-                    cloud.constProps(typeId).internalDegreesOfFreedom()
-                );
+                    // Coefficients required for Bird eqn 12.5
+                    scalar uNormProbCoeffA =
+                        sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
 
-                cloud.addNewParcel
-                (
-                    p,
-                    U,
-                    Ei,
-                    cell,
-                    typeId
-                );
+                    scalar uNormProbCoeffB =
+                        0.5*
+                        (
+                            1.0
+                          + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
+                        );
 
-                particlesInserted++;
+                    // Equivalent to the QA value in Bird's DSMC3.FOR
+                    scalar randomScaling = 3.0;
+
+                    if (sCosTheta < -3)
+                    {
+                        randomScaling = mag(sCosTheta) + 1;
+                    }
+
+                    scalar P = -1;
+
+                    // Normalised candidates for the normal direction velocity
+                    // component
+                    scalar uNormal;
+                    scalar uNormalThermal;
+
+                    // Select a velocity using Bird eqn 12.5
+                    do
+                    {
+                        uNormalThermal =
+                            randomScaling*(2.0*rndGen.scalar01() - 1);
+
+                        uNormal = uNormalThermal + sCosTheta;
+
+                        if (uNormal < 0.0)
+                        {
+                            P = -1;
+                        }
+                        else
+                        {
+                            P = 2.0*uNormal/uNormProbCoeffA
+                               *exp(uNormProbCoeffB - sqr(uNormalThermal));
+                        }
+
+                    } while (P < rndGen.scalar01());
+
+                    vector U =
+                        sqrt(CloudType::kb*faceTemperature/mass)
+                       *(
+                            rndGen.GaussNormal()*t1
+                          + rndGen.GaussNormal()*t2
+                        )
+                      + mostProbableSpeed*uNormal*n;
+
+                    scalar Ei = cloud.equipartitionInternalEnergy
+                    (
+                        faceTemperature,
+                        cloud.constProps(typeId).internalDegreesOfFreedom()
+                    );
+
+                    cloud.addNewParcel
+                    (
+                        p,
+                        U,
+                        Ei,
+                        cell,
+                        typeId
+                    );
+
+                    particlesInserted++;
+                }
             }
         }
     }
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
index f1bd4cc321..03c2fe2abb 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
@@ -26,8 +26,10 @@ 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
+    Inserting new particles across the faces of a all patched of type
+    "patch" for a free stream.  Uniform values number density, temperature
+    and velocity sourced face-by-face from the boundaryT and boundaryU fields
+    of the cloud.
 
 \*---------------------------------------------------------------------------*/
 
@@ -52,14 +54,8 @@ class FreeStream
 {
     // 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 indices of patches to introduce molecules across
+        labelList patches_;
 
         //- The molecule types to be introduced
         List<label> moleculeTypeIds_;
@@ -67,10 +63,13 @@ class FreeStream
         //- 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_;
+        //- A List of Lists of Fields specifying carry-over of mass flux from
+        // one timestep to the next
+        // + Outer List - one inner List for each patch
+        // + Inner List - one Field for every species to be introduced
+        // + Each field entry corresponding to a face to be injected across
+        //   with a particular species
+        List<List<Field<scalar> > > particleFluxAccumulators_;
 
 
 public:
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
index 3aaa9c514a..273c88f369 100644
--- a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
@@ -137,7 +137,7 @@ void Foam::dsmcFields::write()
             iDofMeanName
         );
 
-        if (min(rhoNMean).value() > VSMALL)
+        if (min(mag(rhoNMean)).value() > VSMALL)
         {
             Info<< "Calculating dsmcFields." << endl;
 
@@ -223,10 +223,9 @@ void Foam::dsmcFields::write()
         }
         else
         {
-            Info<< "Small or negative value (" << min(rhoNMean)
+            Info<< "Small value (" << min(mag(rhoNMean))
                 << ") found in rhoNMean field. "
-                << "Not calculating dsmcFields to avoid division by zero "
-                << "or invalid results."
+                << "Not calculating dsmcFields to avoid division by zero."
                 << endl;
         }
     }