From dc7c3ae7af2a94ee014e4d275e1443891c655dee Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Wed, 4 Mar 2009 14:17:42 +0000
Subject: [PATCH] Properly random distribution of particles across inflow patch
 faces.  Modified info() reporting to average per-molecule energy and momentum
 - more useful.

---
 .../clouds/Templates/DsmcCloud/DsmcCloud.C    |  33 +++--
 .../clouds/Templates/DsmcCloud/DsmcCloud.H    |   2 +-
 .../FreeStream/FreeStream.C                   | 116 +++++++++++++-----
 .../MaxwellianThermal/MaxwellianThermal.C     |   4 +-
 4 files changed, 104 insertions(+), 51 deletions(-)

diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index b6d20bd630..c588bb0cb7 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -555,6 +555,11 @@ void Foam::DsmcCloud<ParcelType>::evolve()
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::info() const
 {
+    label nDsmcParticles = this->size();
+    reduce(nDsmcParticles, sumOp<label>());
+
+    scalar nMol = nDsmcParticles*nParticle_;
+
     vector linearMomentum = linearMomentumOfSystem();
     reduce(linearMomentum, sumOp<vector>());
 
@@ -565,20 +570,22 @@ void Foam::DsmcCloud<ParcelType>::info() const
     reduce(internalEnergy, sumOp<scalar>());
 
     Info<< "Cloud name: " << this->name() << nl
-        << "    Current number of parcels       = "
-        << returnReduce(this->size(), sumOp<label>()) << nl
-        << "    Current mass in system          = "
+        << "    Number of dsmc particles        = "
+        << nDsmcParticles << nl
+        << "    Number of molecules             = "
+        << nMol << nl
+        << "    Mass in system                  = "
         << returnReduce(massInSystem(), sumOp<scalar>()) << nl
-        << "    Linear momentum                 = "
-        << linearMomentum << nl
-        << "    Linear momentum magnitude       = "
-        << mag(linearMomentum) << nl
-        << "    Linear kinetic energy           = "
-        << linearKineticEnergy << nl
-        << "    Internal energy                 = "
-        << internalEnergy << nl
-        << "    Total energy                    = "
-        << internalEnergy + linearKineticEnergy << nl
+        << "    Average linear momentum         = "
+        << linearMomentum/nMol << nl
+        << "   |Average linear momentum|        = "
+        << mag(linearMomentum)/nMol << nl
+        << "    Average linear kinetic energy   = "
+        << linearKineticEnergy/nMol << nl
+        << "    Average internal energy         = "
+        << internalEnergy/nMol << nl
+        << "    Average total energy            = "
+        << (internalEnergy + linearKineticEnergy)/nMol << nl
         << endl;
 }
 
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index 5d9edc9c06..d0bb26bc9c 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -97,7 +97,7 @@ class DsmcCloud
         List<DynamicList<ParcelType*> > cellOccupancy_;
 
         //- An IOField holding the value of (sigmaT * cR)max for each cell (see
-        // Bird p220). Initialised with the parcelsm updated as required, and
+        // Bird p220). Initialised with the parcels, updated as required, and
         // read in on start/restart.
         volScalarField sigmaTcRMax_;
 
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index ed55544077..a69f74aa99 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -126,8 +126,59 @@ void Foam::FreeStream<CloudType>::inflow()
     {
         particleFluxAccumulators_[i] +=
            -patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
+    }
 
-        forAll(particleFluxAccumulators_[i], f)
+    forAll(patch, f)
+    {
+        // Loop over all faces as the outer loop to avoid calculating
+        // geometrical properties multiple times.
+
+        labelList faceVertices = patch[f];
+
+        label nVertices = faceVertices.size();
+
+        label globalFaceIndex = f + patch.start();
+
+        label cell = mesh.faceOwner()[globalFaceIndex];
+
+        const vector& fC = patch.faceCentres()[f];
+
+        scalar fA = mag(patch.faceAreas()[f]);
+
+        // Cummulative triangle area fractions
+        List<scalar> cTriAFracs(nVertices);
+
+        for (label v = 0; v < nVertices - 1; v++)
+        {
+            const point& vA = mesh.points()[faceVertices[v]];
+
+            const point& vB = mesh.points()[faceVertices[(v + 1)]];
+
+            cTriAFracs[v] =
+                0.5*mag((vA - fC)^(vB - fC))/fA
+              + cTriAFracs[max((v - 1), 0)];
+        }
+
+        // Force the last area fraction value to 1.0 to avoid any
+        // rounding/non-flat face errors giving a value < 1.0
+        cTriAFracs[nVertices - 1] = 1.0;
+
+        // Normal unit vector *negative* so normal is pointing into the
+        // domain
+        vector nw = patch.faceAreas()[f];
+        nw /= -mag(nw);
+
+        // Wall tangential unit vector. Use the direction between the
+        // face centre and the first vertex in the list
+        vector tw1 = fC - (mesh.points()[faceVertices[0]]);
+        tw1 /= mag(tw1);
+
+        // Other tangential unit vector.  Rescaling in case face is not
+        // flat and nw and tw1 aren't perfectly orthogonal
+        vector tw2 = nw ^ tw1;
+        tw2 /= mag(tw2);
+
+        forAll(particleFluxAccumulators_, i)
         {
             scalar& faceAccumulator = particleFluxAccumulators_[i][f];
 
@@ -140,39 +191,44 @@ void Foam::FreeStream<CloudType>::inflow()
 
             scalar mass = cloud.constProps(typeId).mass();
 
-            labelList faceVertices = patch[f];
-
-            label globalFaceIndex = f + patch.start();
-
-            label cell = mesh.faceOwner()[globalFaceIndex];
-
-            const vector& fC = patch.faceCentres()[f];
-
             for (label n = 0; n < nI; n++)
             {
-                // Temporarily insert particles half way between the face and
-                // cell centres
-                vector p = 0.5*(fC + mesh.cellCentres()[cell]);
+                // Choose a triangle to insert on, based on their relative area
 
-                // Normal unit vector *negative* so normal is pointing into the
-                // domain
-                vector nw = patch.faceAreas()[f];
-                nw /= -mag(nw);
+                scalar triSelection = rndGen.scalar01();
 
-                // 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);
+                // Selected triangle
+                label sTri = -1;
 
-                // 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(cTriAFracs, tri)
+                {
+                    sTri = tri;
 
-                scalar C = sqrt(CloudType::kb*temperature_/mass);
+                    if (cTriAFracs[tri] >= triSelection)
+                    {
+                        break;
+                    }
+                }
+
+                // Randomly distribute the points on the triangle, using the
+                // method from:
+                // Generating Random Points in Triangles
+                // by Greg Turk
+                // from "Graphics Gems", Academic Press, 1990
+                // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
+
+                const point& A = fC;
+                const point& B = mesh.points()[faceVertices[sTri]];
+                const point& C =
+                    mesh.points()[faceVertices[(sTri + 1) % nVertices]];
+
+                scalar s = rndGen.scalar01();
+                scalar t = sqrt(rndGen.scalar01());
+
+                point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
 
                 vector U =
-                    C
+                    sqrt(CloudType::kb*temperature_/mass)
                    *(
                         rndGen.GaussNormal()*tw1
                       + rndGen.GaussNormal()*tw2
@@ -204,14 +260,6 @@ void Foam::FreeStream<CloudType>::inflow()
     Info<< "    Particles inserted              = "
         << particlesInserted << endl;
 
-    // Info<< "insert particles now! " << nl
-    //     << temperature_ << nl
-    //     << velocity_ << nl
-    //     << moleculeTypeIds_ << nl
-    //     << numberDensities_ << nl
-    //     << particleFluxAccumulators_ << nl
-    //     << patch
-    //     << endl;
 }
 
 
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 23cd0bed34..a6ed8e5718 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -102,10 +102,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
 
-    scalar C = sqrt(CloudType::kb*T/mass);
-
     U =
-        C
+        sqrt(CloudType::kb*T/mass)
        *(
             rndGen.GaussNormal()*tw1
           + rndGen.GaussNormal()*tw2