Merge branch 'dsmc' into particleInteractions

applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H

@@ -1,162 +0,0 @@
-
-    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/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C

@@ -41,53 +41,21 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createFields.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    Info<< nl << "Constructing dsmcCloud " << endl;
+
+    dsmcCloud dsmc("dsmc", mesh);
+
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (runTime.loop())
     {
-        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();

applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C

@@ -113,7 +113,16 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
         return;
     }
 
-    wordList extensiveVVFNames(IStringStream ("(momentumMean)")());
+    wordList extensiveVVFNames
+    (
+        IStringStream
+        (
+            "( \
+                momentumMean \
+                fDMean \
+            )"
+        )()
+    );
 
     PtrList<volVectorField> extensiveVVFs(extensiveVVFNames.size());
 

applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C

@@ -44,9 +44,21 @@ int main(int argc, char *argv[])
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    IOdictionary dsmcInitialiseDict
+    (
+        IOobject
+        (
+            "dsmcInitialiseDict",
+            mesh.time().system(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
     Info<< "Initialising dsmc for Time = " << runTime.timeName() << nl << endl;
 
-    dsmcCloud dsmc("dsmc", mesh);
+    dsmcCloud dsmc("dsmc", mesh, dsmcInitialiseDict);
 
     label totalMolecules = dsmc.size();
 

etc/controlDict

@@ -377,6 +377,7 @@ DebugSwitches
     displacementLaplacian 0;
     displacementSBRStress 0;
     distanceSurface     0;
+    distribution        0;
     downwind            0;
     dragModel           0;
     duplicatePoints     0;

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C

@@ -475,21 +475,92 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
 
 template<class ParcelType>
-void Foam::DsmcCloud<ParcelType>::resetSurfaceDataFields()
+void Foam::DsmcCloud<ParcelType>::resetFields()
 {
-    volScalarField::GeometricBoundaryField& qBF(q_.boundaryField());
+    q_ = dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0);
 
-    forAll(qBF, p)
+    fD_ = dimensionedVector
+    (
+        "zero",
+        dimensionSet(1, -1, -2, 0, 0),
+        vector::zero
+    );
+
+    rhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
+
+    rhoM_ =  dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL);
+
+    dsmcRhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0);
+
+    linearKE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
+
+    internalE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
+
+    iDof_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
+
+    momentum_ = dimensionedVector
+    (
+        "zero",
+        dimensionSet(1, -2, -1, 0, 0),
+        vector::zero
+    );
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::calculateFields()
+{
+    scalarField& rhoN = rhoN_.internalField();
+
+    scalarField& rhoM = rhoM_.internalField();
+
+    scalarField& dsmcRhoN = dsmcRhoN_.internalField();
+
+    scalarField& linearKE = linearKE_.internalField();
+
+    scalarField& internalE = internalE_.internalField();
+
+    scalarField& iDof = iDof_.internalField();
+
+    vectorField& momentum = momentum_.internalField();
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
     {
-        qBF[p] = 0.0;
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        rhoN[cellI]++;
+
+        rhoM[cellI] += constProps(p.typeId()).mass();
+
+        dsmcRhoN[cellI]++;
+
+        linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
+
+        internalE[cellI] += p.Ei();
+
+        iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
+
+        momentum[cellI] += constProps(p.typeId()).mass()*p.U();
     }
 
-    volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField());
+    rhoN *= nParticle_/mesh().cellVolumes();
+    rhoN_.correctBoundaryConditions();
 
-    forAll(fDBF, p)
+    rhoM *= nParticle_/mesh().cellVolumes();
-    {
+    rhoM_.correctBoundaryConditions();
-        fDBF[p] = vector::zero;
+
-    }
+    linearKE *= nParticle_/mesh().cellVolumes();
+    linearKE_.correctBoundaryConditions();
+
+    internalE *= nParticle_/mesh().cellVolumes();
+    internalE_.correctBoundaryConditions();
+
+    iDof *= nParticle_/mesh().cellVolumes();
+    iDof_.correctBoundaryConditions();
+
+    momentum *= nParticle_/mesh().cellVolumes();
+    momentum_.correctBoundaryConditions();
 }
 
 
@@ -525,14 +596,13 @@ template<class ParcelType>
 Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudName,
-    const volScalarField& T,
+    const fvMesh& mesh
-    const volVectorField& U
 )
 :
-    Cloud<ParcelType>(T.mesh(), cloudName, false),
+    Cloud<ParcelType>(mesh, cloudName, false),
     DsmcBaseCloud(),
     cloudName_(cloudName),
-    mesh_(T.mesh()),
+    mesh_(mesh),
     particleProperties_
     (
         IOobject
@@ -564,37 +634,142 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     (
         IOobject
         (
-            this->name() + "q_",
+            "q",
             mesh_.time().timeName(),
             mesh_,
-            IOobject::NO_READ,
+            IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        mesh_,
+        mesh_
-        dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
     ),
     fD_
     (
         IOobject
         (
-            this->name() + "fD_",
+            "fD",
             mesh_.time().timeName(),
             mesh_,
-            IOobject::NO_READ,
+            IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        mesh_,
+        mesh_
-        dimensionedVector
+    ),
+    rhoN_
+    (
+        IOobject
         (
-            "zero",
+            "rhoN",
-            dimensionSet(1, -1, -2, 0, 0),
+            mesh_.time().timeName(),
-            vector::zero
+            mesh_,
-        )
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    rhoM_
+    (
+        IOobject
+        (
+            "rhoM",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    dsmcRhoN_
+    (
+        IOobject
+        (
+            "dsmcRhoN",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    linearKE_
+    (
+        IOobject
+        (
+            "linearKE",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    internalE_
+    (
+        IOobject
+        (
+            "internalE",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    iDof_
+    (
+        IOobject
+        (
+            "iDof",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    momentum_
+    (
+        IOobject
+        (
+            "momentum",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
     ),
     constProps_(),
     rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
-    T_(T),
+    boundaryT_
-    U_(U),
+    (
+        volScalarField
+        (
+            IOobject
+            (
+                "boundaryT",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        )
+    ),
+    boundaryU_
+    (
+        volVectorField
+        (
+            IOobject
+            (
+                "boundaryU",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        )
+    ),
     binaryCollisionModel_
     (
         BinaryCollisionModel<DsmcCloud<ParcelType> >::New
@@ -637,7 +812,8 @@ template<class ParcelType>
 Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudName,
-    const fvMesh& mesh
+    const fvMesh& mesh,
+    const IOdictionary& dsmcInitialiseDict
 )
     :
     Cloud<ParcelType>(mesh, cloudName, false),
@@ -703,15 +879,111 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
             vector::zero
         )
     ),
+    rhoN_
+    (
+        IOobject
+        (
+            this->name() + "rhoN_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
+    ),
+    rhoM_
+    (
+        IOobject
+        (
+            this->name() + "rhoM_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL)
+    ),
+    dsmcRhoN_
+    (
+        IOobject
+        (
+            this->name() + "dsmcRhoN_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0)
+    ),
+    linearKE_
+    (
+        IOobject
+        (
+            this->name() + "linearKE_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
+    ),
+    internalE_
+    (
+        IOobject
+        (
+            this->name() + "internalE_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
+    ),
+    iDof_
+    (
+        IOobject
+        (
+            this->name() + "iDof_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
+    ),
+    momentum_
+    (
+        IOobject
+        (
+            this->name() + "momentum_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            dimensionSet(1, -2, -1, 0, 0),
+            vector::zero
+        )
+    ),
     constProps_(),
     rndGen_(label(971501) + 1526*Pstream::myProcNo()),
-    T_
+    boundaryT_
     (
         volScalarField
         (
             IOobject
             (
-                "T",
+                "boundaryT",
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
@@ -721,13 +993,13 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
             dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
         )
     ),
-    U_
+    boundaryU_
     (
         volVectorField
         (
             IOobject
             (
-                "U",
+                "boundaryU",
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
@@ -750,18 +1022,6 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
 
     buildConstProps();
 
-    IOdictionary dsmcInitialiseDict
-    (
-        IOobject
-        (
-            "dsmcInitialiseDict",
-            mesh_.time().system(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
     initialise(dsmcInitialiseDict);
 }
 
@@ -783,8 +1043,8 @@ void Foam::DsmcCloud<ParcelType>::evolve()
 
     typename ParcelType::trackData td(*this);
 
-    // Reset the surface data collection fields
+    // Reset the data collection fields
-    resetSurfaceDataFields();
+    resetFields();
 
     if (debug)
     {
@@ -799,6 +1059,9 @@ void Foam::DsmcCloud<ParcelType>::evolve()
 
     // Calculate new velocities via stochastic collisions
     collisions();
+
+    // Calculate the volume field data
+    calculateFields();
 }
 
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H

@@ -110,6 +110,27 @@ class DsmcCloud
         //- Force density at surface field
         volVectorField fD_;
 
+        //- number density field
+        volScalarField rhoN_;
+
+        //- Mass density field
+        volScalarField rhoM_;
+
+        //- dsmc particle density field
+        volScalarField dsmcRhoN_;
+
+        //- linear kinetic energy density field
+        volScalarField linearKE_;
+
+        //- Internal energy density field
+        volScalarField internalE_;
+
+        // Internal degree of freedom density field
+        volScalarField iDof_;
+
+        //- Momentum density field
+        volVectorField momentum_;
+
         //- Parcel constant properties - one for each type
         List<typename ParcelType::constantProperties> constProps_;
 
@@ -121,13 +142,13 @@ class DsmcCloud
         scalar cachedDeltaT_;
 
 
-        // References to the macroscopic fields
+        // boundary value fields
 
-            //- Temperature
+            //- boundary temperature
-            const volScalarField& T_;
+            volScalarField boundaryT_;
 
-            //- Velocity
+            //- boundary velocity
-            const volVectorField& U_;
+            volVectorField boundaryU_;
 
 
         // References to the cloud sub-models
@@ -159,8 +180,11 @@ class DsmcCloud
         //- Calculate collisions between molecules
         void collisions();
 
-        //- Reset the surface data accumulation field values
+        //- Reset the data accumulation field values to zero
-        void resetSurfaceDataFields();
+        void resetFields();
+
+        //- Calculate the volume field data
+        void calculateFields();
 
         //- Disallow default bitwise copy construct
         DsmcCloud(const DsmcCloud&);
@@ -179,21 +203,22 @@ public:
 
     // Constructors
 
-        //- Construct given name and mesh, will read Parcels from file
+        //- Construct given name and mesh, will read Parcels and fields from
-        DsmcCloud
+        //  file
-        (
-            const word& cloudName,
-            const volScalarField& T,
-            const volVectorField& U
-        );
-
-        //- Construct given name and mesh.  Used to initialise.
         DsmcCloud
         (
             const word& cloudName,
             const fvMesh& mesh
         );
 
+        //- Construct given name, mesh and initialisation dictionary.
+        DsmcCloud
+        (
+            const word& cloudName,
+            const fvMesh& mesh,
+            const IOdictionary& dsmcInitialiseDict
+        );
+
 
     //- Destructor
     virtual ~DsmcCloud();
@@ -254,28 +279,71 @@ public:
                 inline scalar cachedDeltaT() const;
 
 
-            // References to the surface data collection fields
+            // References to the boundary fields for surface data collection
 
-                //- Return heat flux at surface field
+                //- Return non-const heat flux boundary field reference
-                inline const volScalarField& q() const;
+                inline volScalarField::GeometricBoundaryField& qBF();
 
-                //- Return non-const heat flux at surface field
+                //- Return non-const force density at boundary field reference
-                inline volScalarField& q();
+                inline volVectorField::GeometricBoundaryField& fDBF();
 
-                //- Return force density at surface field
+                //- Return non-const number density boundary field reference
-                inline const volVectorField& fD() const;
+                inline volScalarField::GeometricBoundaryField& rhoNBF();
 
-                //- Return non-const force density at surface field
+                //- Return non-const mass density boundary field reference
-                inline volVectorField& fD();
+                inline volScalarField::GeometricBoundaryField& rhoMBF();
+
+                //- Return non-const linear kinetic energy density boundary
+                //  field reference
+                inline volScalarField::GeometricBoundaryField& linearKEBF();
+
+                //- Return non-const internal energy density boundary field
+                // reference
+                inline volScalarField::GeometricBoundaryField& internalEBF();
+
+                //- Return non-const internal degree of freedom density boundary
+                //  field reference
+                inline volScalarField::GeometricBoundaryField& iDofBF();
+
+                //- Return non-const momentum density boundary field reference
+                inline volVectorField::GeometricBoundaryField& momentumBF();
 
 
             // References to the macroscopic fields
 
                 //- Return macroscopic temperature
-                inline const volScalarField& T() const;
+                inline const volScalarField& boundaryT() const;
 
                 //- Return macroscopic velocity
-                inline const volVectorField& U() const;
+                inline const volVectorField& boundaryU() const;
+
+                //- Return heat flux at surface field
+                inline const volScalarField& q() const;
+
+                //- Return force density at surface field
+                inline const volVectorField& fD() const;
+
+                //- Return the real particle number density field
+                inline const volScalarField& rhoN() const;
+
+                //- Return the particle mass density field
+                inline const volScalarField& rhoM() const;
+
+                //- Return the field of number of DSMC particles
+                inline const volScalarField& dsmcRhoN() const;
+
+                //- Return the total linear kinetic energy (translational and
+                // thermal density field
+                inline const volScalarField& linearKE() const;
+
+                //- Return the internal energy density field
+                inline const volScalarField& internalE() const;
+
+                //- Return the average internal degrees of freedom  field
+                inline const volScalarField& iDof() const;
+
+                //- Return the momentum density field
+                inline const volVectorField& momentum() const;
 
 
             // Kinetic theory helper functions
@@ -389,29 +457,6 @@ public:
             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
@@ -429,7 +474,6 @@ public:
             //- Evolve the cloud (move, collide)
             void evolve();
 
-
             //- Clear the Cloud
             inline void clear();
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H

@@ -135,44 +135,82 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::cachedDeltaT() const
 
 
 template<class ParcelType>
-inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::qBF()
 {
-    return q_;
+    return q_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q()
+inline Foam::volVectorField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::fDBF()
 {
-    return q_;
+    return fD_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::rhoNBF()
 {
-    return fD_;
+    return rhoN_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD()
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::rhoMBF()
 {
-    return fD_;
+    return rhoM_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::T() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::linearKEBF()
 {
-    return T_;
+    return linearKE_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::U() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::internalEBF()
 {
-    return U_;
+    return internalE_.boundaryField();
+}
+
+
+template<class ParcelType>
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::iDofBF()
+{
+    return iDof_.boundaryField();
+}
+
+
+template<class ParcelType>
+inline Foam::volVectorField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::momentumBF()
+{
+    return momentum_.boundaryField();
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField&
+Foam::DsmcCloud<ParcelType>::boundaryT() const
+{
+    return boundaryT_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField&
+Foam::DsmcCloud<ParcelType>::boundaryU() const
+{
+    return boundaryU_;
 }
 
 
@@ -374,265 +412,70 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
+{
+    return q_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
+{
+    return fD_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::rhoN() const
 {
-    tmp<volScalarField> trhoN
+    return rhoN_;
-    (
-        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), VSMALL)
-        )
-    );
-
-    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>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::rhoM() const
-Foam::DsmcCloud<ParcelType>::rhoM() const
 {
-    tmp<volScalarField> trhoM
+    return rhoM_;
-    (
-        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), VSMALL)
-        )
-    );
-
-    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>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::dsmcRhoN() const
 {
-    tmp<volScalarField> tdsmcRhoN
+    return dsmcRhoN_;
-    (
-        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>
+inline const Foam::volScalarField&
-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
+    return linearKE_;
-    (
-        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>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::internalE() const
 {
-    tmp<volScalarField> tinternalE
+    return internalE_;
-    (
-        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>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::iDof() const
 {
-    tmp<volScalarField> tiDof
+    return iDof_;
-    (
+}
-        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), VSMALL)
-        )
-    );
 
-    scalarField& iDof = tiDof().internalField();
 
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+template<class ParcelType>
-    {
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::momentum() const
-        const ParcelType& p = iter();
+{
-        const label cellI = p.cell();
+    return momentum_;
-
-        iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
-    }
-
-    iDof *= nParticle_/mesh().cellVolumes();
-
-    return tiDof;
 }
 
 

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

@@ -39,11 +39,10 @@ namespace Foam
 Foam::dsmcCloud::dsmcCloud
 (
     const word& cloudName,
-    const volScalarField& T,
+    const fvMesh& mesh
-    const volVectorField& U
 )
 :
-    DsmcCloud<dsmcParcel>(cloudName, T, U)
+    DsmcCloud<dsmcParcel>(cloudName, mesh)
 {
     dsmcParcel::readFields(*this);
 }
@@ -52,10 +51,11 @@ Foam::dsmcCloud::dsmcCloud
 Foam::dsmcCloud::dsmcCloud
 (
     const word& cloudName,
-    const fvMesh& mesh
+    const fvMesh& mesh,
+    const IOdictionary& dsmcInitialiseDict
 )
 :
-    DsmcCloud<dsmcParcel>(cloudName, mesh)
+    DsmcCloud<dsmcParcel>(cloudName, mesh, dsmcInitialiseDict)
 {}
 
 

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

@@ -69,21 +69,22 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct given name and mesh, will read Parcels and fields from
-        dsmcCloud
+        //  file
-        (
-            const word& cloudName,
-            const volScalarField& T,
-            const volVectorField& U
-        );
-
-        //- Construct from name and mesh, used to initialise.
         dsmcCloud
         (
             const word& cloudName,
             const fvMesh& mesh
         );
 
+        //- Construct given name, mesh and initialisation dictionary.
+        dsmcCloud
+        (
+            const word& cloudName,
+            const fvMesh& mesh,
+            const IOdictionary& dsmcInitialiseDict
+        );
+
     //- Destructor
     ~dsmcCloud();
 

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C

@@ -114,10 +114,41 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
     TrackData& td
 )
 {
+    label wppIndex = wpp.index();
+
+    label wppLocalFace = wpp.whichFace(this->face());
+
+    const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
+
+    const scalar deltaT = td.cloud().cachedDeltaT();
+
     const constantProperties& constProps(td.cloud().constProps(typeId_));
 
     scalar m = constProps.mass();
 
+    vector nw = wpp.faceAreas()[wppLocalFace];
+    nw /= mag(nw);
+
+    scalar U_dot_nw = U_ & nw;
+
+    vector Ut = U_ - U_dot_nw*nw;
+
+    scalar invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);
+
+    td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
+
+    td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
+
+    td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
+        0.5*m*(U_ & U_)*invMagUnfA;
+
+    td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
+
+    td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
+        constProps.internalDegreesOfFreedom()*invMagUnfA;
+
+    td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
+
     // pre-interaction energy
     scalar preIE = 0.5*m*(U_ & U_) + Ei_;
 
@@ -133,27 +164,40 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
         typeId_
     );
 
+    U_dot_nw = U_ & nw;
+
+    Ut = U_ - U_dot_nw*nw;
+
+    invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);
+
+    td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
+
+    td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
+
+    td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
+        0.5*m*(U_ & U_)*invMagUnfA;
+
+    td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
+
+    td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
+        constProps.internalDegreesOfFreedom()*invMagUnfA;
+
+    td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
+
     // 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().cachedDeltaT();
-
     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().qBF()[wppIndex][wppLocalFace] += deltaQ;
+
+    td.cloud().fDBF()[wppIndex][wppLocalFace] += deltaFD;
 
-    td.cloud().fD().boundaryField()[wppIndex][wppLocalFace] += deltaFD;
 }
 
 

src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C

@@ -28,6 +28,7 @@ License
 #include "DsmcCloud.H"
 #include "MaxwellianThermal.H"
 #include "SpecularReflection.H"
+#include "MixedDiffuseSpecular.H"
 
 namespace Foam
 {
@@ -46,6 +47,12 @@ namespace Foam
         DsmcCloud,
         dsmcParcel
     );
+    makeWallInteractionModelType
+    (
+        MixedDiffuseSpecular,
+        DsmcCloud,
+        dsmcParcel
+    );
 };
 
 

src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

@@ -136,12 +136,12 @@ void Foam::FreeStream<CloudType>::inflow()
 
     const volScalarField::GeometricBoundaryField& boundaryT
     (
-        cloud.T().boundaryField()
+        cloud.boundaryT().boundaryField()
     );
 
     const volVectorField::GeometricBoundaryField& boundaryU
     (
-        cloud.U().boundaryField()
+        cloud.boundaryU().boundaryField()
     );
 
     forAll(patches_, p)
@@ -165,7 +165,8 @@ void Foam::FreeStream<CloudType>::inflow()
             if (min(boundaryT[patchI]) < SMALL)
             {
                 FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
-                    << "Zero boundary temperature detected, check boundaryT condition." << nl
+                    << "Zero boundary temperature detected, check boundaryT "
+                    << "condition." << nl
                     << nl << abort(FatalError);
             }
 

src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C

@@ -68,10 +68,10 @@ void Foam::MaxwellianThermal<CloudType>::correct
     nw /= mag(nw);
 
     // Normal velocity magnitude
-    scalar magUn = U & nw;
+    scalar U_dot_nw = U & nw;
 
     // Wall tangential velocity (flow direction)
-    vector Ut = U - magUn*nw;
+    vector Ut = U - U_dot_nw*nw;
 
     CloudType& cloud(this->owner());
 
@@ -90,9 +90,9 @@ void Foam::MaxwellianThermal<CloudType>::correct
             U.z()*(0.8 + 0.2*rndGen.scalar01())
         );
 
-        magUn = U & nw;
+        U_dot_nw = U & nw;
 
-        Ut = U - magUn*nw;
+        Ut = U - U_dot_nw*nw;
     }
 
     // Wall tangential unit vector
@@ -101,7 +101,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
     // Other tangential unit vector
     vector tw2 = nw^tw1;
 
-    scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
+    scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
 
     scalar mass = cloud.constProps(typeId).mass();
 
@@ -115,7 +115,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
           - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
         );
 
-    U += cloud.U().boundaryField()[wppIndex][wppLocalFace];
+    U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
 
     Ei = cloud.equipartitionInternalEnergy(T, iDof);
 }

src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "MixedDiffuseSpecular.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::MixedDiffuseSpecular<CloudType>::MixedDiffuseSpecular
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    WallInteractionModel<CloudType>(dict, cloud, typeName),
+    diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::MixedDiffuseSpecular<CloudType>::~MixedDiffuseSpecular()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::MixedDiffuseSpecular<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 U_dot_nw = U & nw;
+
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    if (diffuseFraction_ > rndGen.scalar01())
+    {
+        // Diffuse reflection
+
+        // Wall tangential velocity (flow direction)
+        vector Ut = U - U_dot_nw*nw;
+
+        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())
+            );
+
+            U_dot_nw = U & nw;
+
+            Ut = U - U_dot_nw*nw;
+        }
+
+        // Wall tangential unit vector
+        vector tw1 = Ut/mag(Ut);
+
+        // Other tangential unit vector
+        vector tw2 = nw^tw1;
+
+        scalar T = cloud.boundaryT().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.boundaryU().boundaryField()[wppIndex][wppLocalFace];
+
+        Ei = cloud.equipartitionInternalEnergy(T, iDof);
+    }
+    else
+    {
+        // Specular reflection
+
+        if (U_dot_nw > 0.0)
+        {
+            U -= 2.0*U_dot_nw*nw;
+        }
+    }
+
+}
+
+
+// ************************************************************************* //

src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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::MixedDiffuseSpecular
+
+Description
+    Wall interaction setting microscopic velocity to a random one drawn from a
+    Maxwellian distribution corresponding to a specified temperature
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MixedDiffuseSpecular_H
+#define MixedDiffuseSpecular_H
+
+#include "WallInteractionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                    Class MixedDiffuseSpecular Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class MixedDiffuseSpecular
+:
+    public WallInteractionModel<CloudType>
+{
+    // Private data
+
+        //- Fraction of wall interactions that are diffuse
+        scalar diffuseFraction_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("MixedDiffuseSpecular");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        MixedDiffuseSpecular
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~MixedDiffuseSpecular();
+
+
+    // 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 "MixedDiffuseSpecular.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C

@@ -63,11 +63,11 @@ void Foam::SpecularReflection<CloudType>::correct
     vector nw = wpp.faceAreas()[wpp.whichFace(faceId)];
     nw /= mag(nw);
 
-    scalar magUn = U & nw;
+    scalar U_dot_nw = U & nw;
 
-    if (magUn > 0.0)
+    if (U_dot_nw > 0.0)
     {
-        U -= 2.0*magUn*nw;
+        U -= 2.0*U_dot_nw*nw;
     }
 }
 

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C

@@ -25,27 +25,49 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "distribution.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    defineTypeNameAndDebug(distribution, 0);
+}
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::distribution::write
+(
+    const fileName& file,
+    const List<Pair<scalar> >& pairs
+)
+{
+    OFstream os(file);
+
+    forAll(pairs, i)
+    {
+        os  << pairs[i].first() << ' ' << pairs[i].second() << nl;
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-distribution::distribution()
+Foam::distribution::distribution()
 :
     Map<label>(),
     binWidth_(1)
 {}
 
 
-distribution::distribution(const scalar binWidth)
+Foam::distribution::distribution(const scalar binWidth)
 :
     Map<label>(),
     binWidth_(binWidth)
 {}
 
 
-distribution::distribution(const distribution& d)
+Foam::distribution::distribution(const distribution& d)
 :
     Map<label>(static_cast< Map<label> >(d)),
     binWidth_(d.binWidth())
@@ -54,13 +76,13 @@ distribution::distribution(const distribution& d)
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-distribution::~distribution()
+Foam::distribution::~distribution()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-label distribution::totalEntries() const
+Foam::label Foam::distribution::totalEntries() const
 {
     label sumOfEntries = 0;
 
@@ -88,7 +110,7 @@ label distribution::totalEntries() const
 }
 
 
-scalar distribution::approxTotalEntries() const
+Foam::scalar Foam::distribution::approxTotalEntries() const
 {
     scalar sumOfEntries = 0;
 
@@ -101,7 +123,7 @@ scalar distribution::approxTotalEntries() const
 }
 
 
-scalar distribution::mean() const
+Foam::scalar Foam::distribution::mean() const
 {
     scalar runningSum = 0;
 
@@ -124,7 +146,7 @@ scalar distribution::mean() const
 }
 
 
-scalar distribution::median()
+Foam::scalar Foam::distribution::median()
 {
     // From:
     // http://mathworld.wolfram.com/StatisticalMedian.html
@@ -188,7 +210,7 @@ scalar distribution::median()
 }
 
 
-void distribution::add(const scalar valueToAdd)
+void Foam::distribution::add(const scalar valueToAdd)
 {
     iterator iter(this->begin());
 
@@ -218,13 +240,13 @@ void distribution::add(const scalar valueToAdd)
 }
 
 
-void distribution::add(const label valueToAdd)
+void Foam::distribution::add(const label valueToAdd)
 {
     add(scalar(valueToAdd));
 }
 
 
-void distribution::insertMissingKeys()
+void Foam::distribution::insertMissingKeys()
 {
     iterator iter(this->begin());
 
@@ -247,7 +269,7 @@ void distribution::insertMissingKeys()
 }
 
 
-List< Pair<scalar> > distribution::normalised()
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalised()
 {
     scalar totEnt = approxTotalEntries();
 
@@ -268,17 +290,25 @@ List< Pair<scalar> > distribution::normalised()
         normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_;
     }
 
+    if (debug)
+    {
+        Info<< "totEnt: " << totEnt << endl;
+    }
+
     return normDist;
 }
 
 
-List< Pair<scalar> > distribution::normalisedMinusMean()
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedMinusMean()
 {
     return normalisedShifted(mean());
 }
 
 
-List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedShifted
+(
+    scalar shiftValue
+)
 {
     List<Pair<scalar> > oldDist(normalised());
 
@@ -301,20 +331,23 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
 
     label newKey = lowestNewKey;
 
-//     Info << shiftValue
+    if (debug)
-//         << nl << lowestOldBin
+    {
-//         << nl << lowestNewKey
+        Info<< shiftValue
-//         << nl << interpolationStartDirection
+            << nl << lowestOldBin
-//         << endl;
+            << nl << lowestNewKey
+            << nl << interpolationStartDirection
+            << endl;
 
-//     scalar checkNormalisation = 0;
+        scalar checkNormalisation = 0;
 
-//     forAll (oldDist, oD)
+        forAll (oldDist, oD)
-//     {
+        {
-//         checkNormalisation += oldDist[oD].second()*binWidth_;
+            checkNormalisation += oldDist[oD].second()*binWidth_;
-//     }
+        }
 
-//     Info << "Initial normalisation = " << checkNormalisation << endl;
+        Info<< "Initial normalisation = " << checkNormalisation << endl;
+    }
 
     forAll(oldDist,u)
     {
@@ -368,20 +401,23 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
         newKey++;
     }
 
-//     checkNormalisation = 0;
+    if (debug)
+    {
+        scalar checkNormalisation = 0;
 
-//     forAll (newDist, nD)
+        forAll (newDist, nD)
-//     {
+        {
-//         checkNormalisation += newDist[nD].second()*binWidth_;
+            checkNormalisation += newDist[nD].second()*binWidth_;
-//     }
+        }
 
-//     Info << "Shifted normalisation = " << checkNormalisation << endl;
+        Info<< "Shifted normalisation = " << checkNormalisation << endl;
+    }
 
     return newDist;
 }
 
 
-List<Pair<scalar> > distribution::raw()
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::raw()
 {
     insertMissingKeys();
 
@@ -406,7 +442,7 @@ List<Pair<scalar> > distribution::raw()
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void distribution::operator=(const distribution& rhs)
+void Foam::distribution::operator=(const distribution& rhs)
 {
     // Check for assignment to self
     if (this == &rhs)
@@ -424,7 +460,7 @@ void distribution::operator=(const distribution& rhs)
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const distribution& d)
+Foam::Ostream& Foam::operator<<(Ostream& os, const distribution& d)
 {
     os  << d.binWidth_
         << static_cast<const Map<label>&>(d);
@@ -440,8 +476,4 @@ Ostream& operator<<(Ostream& os, const distribution& d)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H

@@ -26,6 +26,8 @@ Class
     Foam::distribution
 
 Description
+    Accumulating histogram of values.  Specified bin resolution
+    automatic generation of bins.
 
 SourceFiles
     distributionI.H
@@ -59,6 +61,21 @@ class distribution
 
 public:
 
+    //- Runtime type information
+
+        TypeName("distribution");
+
+    // Static functions
+
+        //- write to file
+
+            static void write
+            (
+                const fileName& file,
+                const List<Pair<scalar> >& pairs
+            );
+
+
     // Constructors
 
         //- Construct null
@@ -73,7 +90,7 @@ public:
 
     // Destructor
 
-        ~distribution();
+        virtual ~distribution();
 
 
     // Member Functions
@@ -97,7 +114,7 @@ public:
 
         List<Pair<scalar> > normalisedMinusMean();
 
-        List<Pair<scalar> > normalisedShifted(const scalar shiftValue);
+        List<Pair<scalar> > normalisedShifted(scalar shiftValue);
 
         List<Pair<scalar> > raw();
 

src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C

@@ -106,6 +106,7 @@ void Foam::dsmcFields::write()
         word linearKEMeanName = "linearKEMean";
         word internalEMeanName = "internalEMean";
         word iDofMeanName = "iDofMean";
+        word fDMeanName = "fDMean";
 
         const volScalarField& rhoNMean = obr_.lookupObject<volScalarField>
         (
@@ -137,6 +138,11 @@ void Foam::dsmcFields::write()
             iDofMeanName
         );
 
+        volVectorField fDMean =  obr_.lookupObject<volVectorField>
+        (
+            fDMeanName
+        );
+
         if (min(mag(rhoNMean)).value() > VSMALL)
         {
             Info<< "Calculating dsmcFields." << endl;
@@ -165,7 +171,7 @@ void Foam::dsmcFields::write()
                     IOobject::NO_READ
                 ),
                 2.0/(3.0*dsmcCloud::kb*rhoNMean)
-                *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
             );
 
             Info<< "    Calculating internalT field." << endl;
@@ -178,7 +184,7 @@ void Foam::dsmcFields::write()
                     obr_,
                     IOobject::NO_READ
                 ),
-                2.0/(dsmcCloud::kb*iDofMean)*internalEMean
+                (2.0/dsmcCloud::kb)*(internalEMean/iDofMean)
             );
 
             Info<< "    Calculating overallT field." << endl;
@@ -192,9 +198,36 @@ void Foam::dsmcFields::write()
                     IOobject::NO_READ
                 ),
                 2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean))
-                *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
             );
 
+            Info<< "    Calculating pressure field." << endl;
+            volScalarField p
+            (
+                IOobject
+                (
+                    "p",
+                    obr_.time().timeName(),
+                    obr_,
+                    IOobject::NO_READ
+                ),
+                dsmcCloud::kb*rhoNMean*translationalT
+            );
+
+            const fvMesh& mesh = fDMean.mesh();
+
+            forAll(mesh.boundaryMesh(), i)
+            {
+                const polyPatch& patch = mesh.boundaryMesh()[i];
+
+                if (isA<wallPolyPatch>(patch))
+                {
+                    p.boundaryField()[i] =
+                        fDMean.boundaryField()[i]
+                      & (patch.faceAreas()/mag(patch.faceAreas()));
+                }
+            }
+
             Info<< "    mag(UMean) max/min : "
                 << max(mag(UMean)).value() << " "
                 << min(mag(UMean)).value() << endl;
@@ -211,6 +244,10 @@ void Foam::dsmcFields::write()
                 << max(overallT).value() << " "
                 << min(overallT).value() << endl;
 
+            Info<< "    p max/min : "
+                << max(p).value() << " "
+                << min(p).value() << endl;
+
             UMean.write();
 
             translationalT.write();
@@ -219,6 +256,8 @@ void Foam::dsmcFields::write()
 
             overallT.write();
 
+            p.write();
+
             Info<< "dsmcFields written." << nl << endl;
         }
         else