diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 4c4bcd0a43dc965a3bd1daae9b06801b0d359c77..ac3b625638d165c82577210fa16cfefba12430f8 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -113,87 +113,161 @@ void Foam::DsmcCloud<ParcelType>::initialise
numberDensities /= nParticle_;
- scalar x0 = mesh_.bounds().min().x();
- scalar xR = mesh_.bounds().max().x() - x0;
+ forAll(mesh_.cells(), cell)
+ {
+ const vector& cC = mesh_.cellCentres()[cell];
+ const labelList& cellFaces = mesh_.cells()[cell];
+ const scalar cV = mesh_.cellVolumes()[cell];
- scalar y0 = mesh_.bounds().min().y();
- scalar yR = mesh_.bounds().max().y() - y0;
+ label nTets = 0;
- scalar z0 = mesh_.bounds().min().z();
- scalar zR = mesh_.bounds().max().z() - z0;
+ // Each face is split into nEdges (or nVertices) - 2 tets.
+ forAll(cellFaces, face)
+ {
+ nTets += mesh_.faces()[cellFaces[face]].size() - 2;
+ }
- forAll(molecules, i)
- {
- const word& moleculeName(molecules[i]);
+ // Calculate the cumulative tet volumes circulating around the cell and
+ // record the vertex labels of each.
+ scalarList cTetVFracs(nTets, 0.0);
+
+ List<labelList> tetPtIs(nTets, labelList(3,-1));
- label typeId(findIndex(typeIdList_, moleculeName));
+ // Keep track of which tet this is.
+ label tet = 0;
- if (typeId == -1)
+ forAll(cellFaces, face)
{
- FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
- << "typeId " << moleculeName << "not defined." << nl
- << abort(FatalError);
- }
+ const labelList& facePoints = mesh_.faces()[cellFaces[face]];
- const typename ParcelType::constantProperties& cP = constProps(typeId);
+ label pointI = 1;
+ while ((pointI + 1) < facePoints.size())
+ {
- scalar numberDensity = numberDensities[i];
+ const vector& pA = mesh_.points()[facePoints[0]];
+ const vector& pB = mesh_.points()[facePoints[pointI]];
+ const vector& pC = mesh_.points()[facePoints[pointI + 1]];
- scalar spacing = pow(numberDensity,-(1.0/3.0));
+ cTetVFracs[tet] =
+ mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
+ + cTetVFracs[max((tet - 1),0)];
- int ni = label(xR/spacing) + 1;
- int nj = label(yR/spacing) + 1;
- int nk = label(zR/spacing) + 1;
+ tetPtIs[tet][0] = facePoints[0];
+ tetPtIs[tet][1] = facePoints[pointI];
+ tetPtIs[tet][2] = facePoints[pointI + 1];
- vector delta(xR/ni, yR/nj, zR/nk);
+ pointI++;
+ tet++;
+ }
+ }
- scalar pert = spacing;
+ // 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;
- for (int i = 0; i < ni; i++)
+ forAll(molecules, i)
{
- for (int j = 0; j < nj; j++)
+ const word& moleculeName(molecules[i]);
+
+ label typeId(findIndex(typeIdList_, moleculeName));
+
+ if (typeId == -1)
{
- for (int k = 0; k < nk; k++)
- {
- point p
- (
- x0 + (i + 0.5)*delta.x(),
- y0 + (j + 0.5)*delta.y(),
- z0 + (k + 0.5)*delta.z()
- );
+ FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
+ << "typeId " << moleculeName << "not defined." << nl
+ << abort(FatalError);
+ }
- p.x() += pert*(rndGen_.scalar01() - 0.5);
- p.y() += pert*(rndGen_.scalar01() - 0.5);
- p.z() += pert*(rndGen_.scalar01() - 0.5);
+ const typename ParcelType::constantProperties& cP =
+ constProps(typeId);
- label cell = mesh_.findCell(p);
+ scalar numberDensity = numberDensities[i];
- vector U = equipartitionLinearVelocity
- (
- temperature,
- cP.mass()
- );
+ // Calculate the number of particles required
+ scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
- scalar Ei = equipartitionInternalEnergy
- (
- temperature,
- cP.internalDegreesOfFreedom()
- );
+ // 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++;
+ }
- U += velocity;
+ 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 (cell >= 0)
+ if (s + t > 1.0)
+ {
+ s = 1.0 - s;
+ t = 1.0 - t;
+ }
+
+ 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;
+ }
+
+ // Choose a tetrahedron to insert in, based on their relative
+ // volumes
+ scalar tetSelection = rndGen_.scalar01();
+
+ // Selected tetrahedron
+ label sTet = -1;
+
+ forAll(cTetVFracs, tet)
+ {
+ sTet = tet;
+
+ 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
+ );
}
}
}
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 73899636e8f6cf2a7e47754a15c8cf8aeb711712..ede1e3559bb9dbb76619f45f05bdbd5e17a4515a 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -343,7 +343,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 +380,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 +568,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 f4203f4d66f7e85eaf335eee7c5c3b80521dd4e0..02ed6d28b5852cb35ef1eff08e94a1cc044d8a7a 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -157,7 +157,7 @@ void Foam::FreeStream<CloudType>::inflow()
scalar fA = mag(patch.faceAreas()[f]);
- // Cummulative triangle area fractions
+ // Cumulative triangle area fractions
List<scalar> cTriAFracs(nVertices);
for (label v = 0; v < nVertices - 1; v++)