lagrangian: Move penetration function from KinematicCloud to SprayCloud

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=992

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -400,7 +400,6 @@ public:
 
 
             //- Optional particle forces
-//            inline const typename parcelType::forceType& forces() const;
             inline const forceType& forces() const;
 
             //- Return the optional particle forces
@@ -499,9 +498,6 @@ public:
             //- Total rotational kinetic energy in the system
             inline scalar rotationalKineticEnergyOfSystem() const;
 
-            //- Penetration for fraction [0-1] of the current total mass
-            inline scalar penetration(const scalar fraction) const;
-
             //- Mean diameter Dij
             inline scalar Dij(const label i, const label j) const;
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -355,142 +355,6 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 }
 
 
-template<class CloudType>
-inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
-(
-    const scalar fraction
-) const
-{
-    if ((fraction < 0) || (fraction > 1))
-    {
-        FatalErrorInFunction
-            << "fraction should be in the range 0 < fraction < 1"
-            << exit(FatalError);
-    }
-
-    scalar distance = 0.0;
-
-    const label nParcel = this->size();
-    globalIndex globalParcels(nParcel);
-    const label nParcelSum = globalParcels.size();
-
-    if (nParcelSum == 0)
-    {
-        return distance;
-    }
-
-    // lists of parcels mass and distance from initial injection point
-    List<List<scalar>> procMass(Pstream::nProcs());
-    List<List<scalar>> procDist(Pstream::nProcs());
-
-    List<scalar>& mass = procMass[Pstream::myProcNo()];
-    List<scalar>& dist = procDist[Pstream::myProcNo()];
-
-    mass.setSize(nParcel);
-    dist.setSize(nParcel);
-
-    label i = 0;
-    scalar mSum = 0.0;
-    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
-    {
-        const parcelType& p = iter();
-        scalar m = p.nParticle()*p.mass();
-        scalar d = mag(p.position() - p.position0());
-        mSum += m;
-
-        mass[i] = m;
-        dist[i] = d;
-
-        i++;
-    }
-
-    // calculate total mass across all processors
-    reduce(mSum, sumOp<scalar>());
-    Pstream::gatherList(procMass);
-    Pstream::gatherList(procDist);
-
-    if (Pstream::master())
-    {
-        // flatten the mass lists
-        List<scalar> allMass(nParcelSum, 0.0);
-        SortableList<scalar> allDist(nParcelSum, 0.0);
-        for (label procI = 0; procI < Pstream::nProcs(); procI++)
-        {
-            SubList<scalar>
-            (
-                allMass,
-                globalParcels.localSize(procI),
-                globalParcels.offset(procI)
-            ).assign(procMass[procI]);
-
-            // flatten the distance list
-            SubList<scalar>
-            (
-                allDist,
-                globalParcels.localSize(procI),
-                globalParcels.offset(procI)
-            ).assign(procDist[procI]);
-        }
-
-        // sort allDist distances into ascending order
-        // note: allMass masses are left unsorted
-        allDist.sort();
-
-        if (nParcelSum > 1)
-        {
-            const scalar mLimit = fraction*mSum;
-            const labelList& indices = allDist.indices();
-
-            if (mLimit > (mSum - allMass[indices.last()]))
-            {
-                distance = allDist.last();
-            }
-            else
-            {
-                // assuming that 'fraction' is generally closer to 1 than 0,
-                // loop through in reverse distance order
-                const scalar mThreshold = (1.0 - fraction)*mSum;
-                scalar mCurrent = 0.0;
-                label i0 = 0;
-
-                forAllReverse(indices, i)
-                {
-                    label indI = indices[i];
-
-                    mCurrent += allMass[indI];
-
-                    if (mCurrent > mThreshold)
-                    {
-                        i0 = i;
-                        break;
-                    }
-                }
-
-                if (i0 == indices.size() - 1)
-                {
-                    distance = allDist.last();
-                }
-                else
-                {
-                    // linearly interpolate to determine distance
-                    scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
-                    distance =
-                        allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
-                }
-            }
-        }
-        else
-        {
-            distance = allDist.first();
-        }
-    }
-
-    Pstream::scatter(distance);
-
-    return distance;
-}
-
-
 template<class CloudType>
 inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen()
 {

src/lagrangian/intermediate/clouds/baseClasses/kinematicCloud/kinematicCloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,12 +85,6 @@ public:
             //- Total linear kinetic energy in the system
             virtual scalar linearKineticEnergyOfSystem() const = 0;
 
-            //- Total rotational kinetic energy in the system
-//            virtual scalar rotationalKineticEnergyOfSystem() const = 0;
-
-            //- Penetration for percentage of the current total mass
-//            virtual scalar penetration(const scalar& fraction) const = 0;
-
             //- Mean diameter Dij
             virtual scalar Dij(const label i, const label j) const = 0;
 

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H

@@ -183,6 +183,12 @@ public:
             inline scalar averageParcelMass() const;
 
 
+        // Check
+
+            //- Penetration for fraction [0-1] of the current total mass
+            inline scalar penetration(const scalar fraction) const;
+
+
             // Sub-models
 
                 //- Return const-access to the atomization model

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H

@@ -72,4 +72,140 @@ inline Foam::scalar Foam::SprayCloud<CloudType>::averageParcelMass() const
 }
 
 
+template<class CloudType>
+inline Foam::scalar Foam::SprayCloud<CloudType>::penetration
+(
+    const scalar fraction
+) const
+{
+    if ((fraction < 0) || (fraction > 1))
+    {
+        FatalErrorInFunction
+            << "fraction should be in the range 0 < fraction < 1"
+            << exit(FatalError);
+    }
+
+    scalar distance = 0.0;
+
+    const label nParcel = this->size();
+    globalIndex globalParcels(nParcel);
+    const label nParcelSum = globalParcels.size();
+
+    if (nParcelSum == 0)
+    {
+        return distance;
+    }
+
+    // lists of parcels mass and distance from initial injection point
+    List<List<scalar>> procMass(Pstream::nProcs());
+    List<List<scalar>> procDist(Pstream::nProcs());
+
+    List<scalar>& mass = procMass[Pstream::myProcNo()];
+    List<scalar>& dist = procDist[Pstream::myProcNo()];
+
+    mass.setSize(nParcel);
+    dist.setSize(nParcel);
+
+    label i = 0;
+    scalar mSum = 0.0;
+    forAllConstIter(typename SprayCloud<CloudType>, *this, iter)
+    {
+        const parcelType& p = iter();
+        scalar m = p.nParticle()*p.mass();
+        scalar d = mag(p.position() - p.position0());
+        mSum += m;
+
+        mass[i] = m;
+        dist[i] = d;
+
+        i++;
+    }
+
+    // calculate total mass across all processors
+    reduce(mSum, sumOp<scalar>());
+    Pstream::gatherList(procMass);
+    Pstream::gatherList(procDist);
+
+    if (Pstream::master())
+    {
+        // flatten the mass lists
+        List<scalar> allMass(nParcelSum, 0.0);
+        SortableList<scalar> allDist(nParcelSum, 0.0);
+        for (label procI = 0; procI < Pstream::nProcs(); procI++)
+        {
+            SubList<scalar>
+            (
+                allMass,
+                globalParcels.localSize(procI),
+                globalParcels.offset(procI)
+            ).assign(procMass[procI]);
+
+            // flatten the distance list
+            SubList<scalar>
+            (
+                allDist,
+                globalParcels.localSize(procI),
+                globalParcels.offset(procI)
+            ).assign(procDist[procI]);
+        }
+
+        // sort allDist distances into ascending order
+        // note: allMass masses are left unsorted
+        allDist.sort();
+
+        if (nParcelSum > 1)
+        {
+            const scalar mLimit = fraction*mSum;
+            const labelList& indices = allDist.indices();
+
+            if (mLimit > (mSum - allMass[indices.last()]))
+            {
+                distance = allDist.last();
+            }
+            else
+            {
+                // assuming that 'fraction' is generally closer to 1 than 0,
+                // loop through in reverse distance order
+                const scalar mThreshold = (1.0 - fraction)*mSum;
+                scalar mCurrent = 0.0;
+                label i0 = 0;
+
+                forAllReverse(indices, i)
+                {
+                    label indI = indices[i];
+
+                    mCurrent += allMass[indI];
+
+                    if (mCurrent > mThreshold)
+                    {
+                        i0 = i;
+                        break;
+                    }
+                }
+
+                if (i0 == indices.size() - 1)
+                {
+                    distance = allDist.last();
+                }
+                else
+                {
+                    // linearly interpolate to determine distance
+                    scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
+                    distance =
+                        allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
+                }
+            }
+        }
+        else
+        {
+            distance = allDist.first();
+        }
+    }
+
+    Pstream::scatter(distance);
+
+    return distance;
+}
+
+
 // ************************************************************************* //