ENH: Patch injection model - refactored and improved injection location functionality

src/lagrangian/intermediate/Make/files

@@ -54,13 +54,13 @@ KINEMATICINJECTION=submodels/Kinematic/InjectionModel
 $(KINEMATICINJECTION)/KinematicLookupTableInjection/kinematicParcelInjectionData.C
 $(KINEMATICINJECTION)/KinematicLookupTableInjection/kinematicParcelInjectionDataIO.C
 $(KINEMATICINJECTION)/KinematicLookupTableInjection/kinematicParcelInjectionDataIOList.C
+$(KINEMATICINJECTION)/PatchInjection/patchInjectionBase.C
 
 THERMOINJECTION=submodels/Thermodynamic/InjectionModel
 $(THERMOINJECTION)/ThermoLookupTableInjection/thermoParcelInjectionData.C
 $(THERMOINJECTION)/ThermoLookupTableInjection/thermoParcelInjectionDataIO.C
 $(THERMOINJECTION)/ThermoLookupTableInjection/thermoParcelInjectionDataIOList.C
 
-
 REACTINGINJECTION=submodels/Reacting/InjectionModel
 $(REACTINGINJECTION)/ReactingLookupTableInjection/reactingParcelInjectionData.C
 $(REACTINGINJECTION)/ReactingLookupTableInjection/reactingParcelInjectionDataIO.C

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,8 +38,7 @@ Foam::PatchInjection<CloudType>::PatchInjection
 )
 :
     InjectionModel<CloudType>(dict, owner, modelName, typeName),
-    patchName_(this->coeffDict().lookup("patchName")),
+    patchInjectionBase(owner.mesh(), this->coeffDict().lookup("patchName")),
-    patchId_(owner.mesh().boundaryMesh().findPatchID(patchName_)),
     duration_(readScalar(this->coeffDict().lookup("duration"))),
     parcelsPerSecond_
     (
@@ -62,32 +61,11 @@ Foam::PatchInjection<CloudType>::PatchInjection
             this->coeffDict().subDict("sizeDistribution"),
             owner.rndGen()
         )
-    ),
+    )
-    cellOwners_(),
-    fraction_(1.0)
 {
-    if (patchId_ < 0)
-    {
-        FatalErrorIn
-        (
-            "PatchInjection<CloudType>::PatchInjection"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )   << "Requested patch " << patchName_ << " not found" << nl
-            << "Available patches are: " << owner.mesh().boundaryMesh().names()
-            << nl << exit(FatalError);
-    }
-
     duration_ = owner.db().time().userTimeToTime(duration_);
 
-    updateMesh();
+    patchInjectionBase::updateMesh(owner.mesh());
-
-    label patchSize = cellOwners_.size();
-    label totalPatchSize = patchSize;
-    reduce(totalPatchSize, sumOp<label>());
-    fraction_ = scalar(patchSize)/totalPatchSize;
 
     // Set total volume/mass to inject
     this->volumeTotal_ = fraction_*flowRateProfile_.integrate(0.0, duration_);
@@ -102,15 +80,12 @@ Foam::PatchInjection<CloudType>::PatchInjection
 )
 :
     InjectionModel<CloudType>(im),
-    patchName_(im.patchName_),
+    patchInjectionBase(im),
-    patchId_(im.patchId_),
     duration_(im.duration_),
     parcelsPerSecond_(im.parcelsPerSecond_),
     U0_(im.U0_),
     flowRateProfile_(im.flowRateProfile_),
-    sizeDistribution_(im.sizeDistribution_().clone().ptr()),
+    sizeDistribution_(im.sizeDistribution_().clone().ptr())
-    cellOwners_(im.cellOwners_),
-    fraction_(im.fraction_)
 {}
 
 
@@ -126,9 +101,7 @@ Foam::PatchInjection<CloudType>::~PatchInjection()
 template<class CloudType>
 void Foam::PatchInjection<CloudType>::updateMesh()
 {
-    // Set/cache the injector cells
+    patchInjectionBase::updateMesh(this->owner().mesh());
-    const polyPatch& patch = this->owner().mesh().boundaryMesh()[patchId_];
-    cellOwners_ = patch.faceCells();
 }
 
 
@@ -148,7 +121,7 @@ Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
 {
     if ((time0 >= 0.0) && (time0 < duration_))
     {
-        scalar nParcels = fraction_*(time1 - time0)*parcelsPerSecond_;
+        scalar nParcels = this->fraction_*(time1 - time0)*parcelsPerSecond_;
 
         cachedRandom& rnd = this->owner().rndGen();
 
@@ -186,7 +159,7 @@ Foam::scalar Foam::PatchInjection<CloudType>::volumeToInject
 {
     if ((time0 >= 0.0) && (time0 < duration_))
     {
-        return fraction_*flowRateProfile_.integrate(time0, time1);
+        return this->fraction_*flowRateProfile_.integrate(time0, time1);
     }
     else
     {
@@ -207,39 +180,15 @@ void Foam::PatchInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    if (cellOwners_.size() > 0)
+    patchInjectionBase::setPositionAndCell
-    {
+    (
-        cachedRandom& rnd = this->owner().rndGen();
+        this->owner().mesh(),
-
+        this->owner().rndGen(),
-        label cellI = rnd.position<label>(0, cellOwners_.size() - 1);
+        position,
-
+        cellOwner,
-        cellOwner = cellOwners_[cellI];
+        tetFaceI,
-
+        tetPtI
-        // The position is between the face and cell centre, which could be
+    );
-        // in any tet of the decomposed cell, so arbitrarily choose the first
-        // face of the cell as the tetFace and the first point after the base
-        // point on the face as the tetPt.  The tracking will pick the cell
-        // consistent with the motion in the firsttracking step.
-        tetFaceI = this->owner().mesh().cells()[cellOwner][0];
-        tetPtI = 1;
-
-        // position perturbed between cell and patch face centres
-        const vector& pc = this->owner().mesh().C()[cellOwner];
-        const vector& pf =
-            this->owner().mesh().Cf().boundaryField()[patchId_][cellI];
-
-        const scalar a = rnd.sample01<scalar>();
-        const vector d = pf - pc;
-        position = pc + 0.5*a*d;
-    }
-    else
-    {
-        cellOwner = -1;
-        tetFaceI = -1;
-        tetPtI = -1;
-        // dummy position
-        position = pTraits<vector>::max;
-    }
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,7 @@ Description
       - Initial parcel velocity
       - Injection volume flow rate
     - Parcel diameters obtained by distribution model
-    - Parcels injected at cell centres adjacent to patch
+    - Parcels injected randomly across the patch
 
 SourceFiles
     PatchInjection.C
@@ -44,6 +44,7 @@ SourceFiles
 #define PatchInjection_H
 
 #include "InjectionModel.H"
+#include "patchInjectionBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,16 +63,11 @@ class distributionModel;
 template<class CloudType>
 class PatchInjection
 :
-    public InjectionModel<CloudType>
+    public InjectionModel<CloudType>,
+    public patchInjectionBase
 {
     // Private data
 
-        //- Name of patch
-        const word patchName_;
-
-        //- Id of patch
-        const label patchId_;
-
         //- Injection duration [s]
         scalar duration_;
 
@@ -87,12 +83,6 @@ class PatchInjection
         //- Parcel size distribution model
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
-        //- List of cell labels corresponding to injector positions
-        labelList cellOwners_;
-
-        //- Fraction of injection controlled by this processor
-        scalar fraction_;
-
 
 public:
 
@@ -133,7 +123,7 @@ public:
         virtual void updateMesh();
 
         //- Return the end-of-injection time
-        scalar timeEnd() const;
+        virtual scalar timeEnd() const;
 
         //- Number of parcels to introduce relative to SOI
         virtual label parcelsToInject(const scalar time0, const scalar time1);

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.C

@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchInjectionBase.H"
+#include "polyMesh.H"
+#include "SubField.H"
+#include "cachedRandom.H"
+#include "triPointRef.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchInjectionBase::patchInjectionBase
+(
+    const polyMesh& mesh,
+    const word& patchName
+)
+:
+    patchName_(patchName),
+    patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
+    patchArea_(0.0),
+    patchNormal_(),
+    cellOwners_(),
+    fraction_(1.0),
+    triFace_(),
+    triToFace_(),
+    triCumulativeMagSf_(),
+    sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
+{
+    if (patchId_ < 0)
+    {
+        FatalErrorIn
+        (
+            "Foam::patchInjectionBase::patchInjectionBase"
+            "("
+                "const polyMesh&, "
+                "const word&"
+            ")"
+        )   << "Requested patch " << patchName_ << " not found" << nl
+            << "Available patches are: " << mesh.boundaryMesh().names()
+            << nl << exit(FatalError);
+    }
+
+    updateMesh(mesh);
+}
+
+
+Foam::patchInjectionBase::patchInjectionBase(const patchInjectionBase& pib)
+:
+    patchName_(pib.patchName_),
+    patchId_(pib.patchId_),
+    patchArea_(pib.patchArea_),
+    patchNormal_(pib.patchNormal_),
+    cellOwners_(pib.cellOwners_),
+    fraction_(pib.fraction_),
+    triFace_(pib.triFace_),
+    triToFace_(pib.triToFace_),
+    triCumulativeMagSf_(pib.triCumulativeMagSf_),
+    sumTriMagSf_(pib.sumTriMagSf_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchInjectionBase::~patchInjectionBase()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
+{
+    // Set/cache the injector cells
+    const polyPatch& patch = mesh.boundaryMesh()[patchId_];
+    const pointField& points = patch.points();
+
+    cellOwners_ = patch.faceCells();
+
+    // triangulate the patch faces and create addressing
+    DynamicList<label> triToFace(2*patch.size());
+    DynamicList<scalar> triMagSf(2*patch.size());
+    DynamicList<face> triFace(2*patch.size());
+    DynamicList<face> tris(5);
+
+    forAll(patch, faceI)
+    {
+        const face& f = patch[faceI];
+
+        tris.clear();
+        f.triangles(points, tris);
+
+        forAll(tris, i)
+        {
+            triToFace.append(faceI);
+            triFace.append(tris[i]);
+            triMagSf.append(tris[i].mag(points));
+        }
+    }
+
+    forAll(sumTriMagSf_, i)
+    {
+        sumTriMagSf_[i] = 0.0;
+    }
+
+    sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
+
+    Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
+    Pstream::listCombineScatter(sumTriMagSf_);
+
+    for (label i = 1; i < triMagSf.size(); i++)
+    {
+        triMagSf[i] += triMagSf[i-1];
+    }
+
+    // transfer to persistent storage
+    triFace_.transfer(triFace);
+    triToFace_.transfer(triToFace);
+    triCumulativeMagSf_.transfer(triMagSf);
+
+    // fraction of injection volume to be injected by this patch
+    fraction_ = sumTriMagSf_[Pstream::myProcNo() + 1]/sum(sumTriMagSf_);
+
+    // convert sumTriMagSf_ into cumulative sum of areas per proc
+    for (label i = 1; i < sumTriMagSf_.size(); i++)
+    {
+        sumTriMagSf_[i] += sumTriMagSf_[i - 1];
+    }
+
+    const scalarField magSf(mag(patch.faceAreas()));
+    patchArea_ = sum(magSf);
+    patchNormal_ = patch.faceAreas()/magSf;
+    reduce(patchArea_, sumOp<scalar>());
+}
+
+
+void Foam::patchInjectionBase::setPositionAndCell
+(
+    const polyMesh& mesh,
+    cachedRandom& rnd,
+    vector& position,
+    label& cellOwner,
+    label& tetFaceI,
+    label& tetPtI
+)
+{
+    if (cellOwners_.size() > 0)
+    {
+        // determine which processor to inject from
+        scalar areaFraction = rnd.position<scalar>(0, patchArea_);
+
+        label procI = 0;
+        forAllReverse(sumTriMagSf_, i)
+        {
+            if (areaFraction > sumTriMagSf_[i])
+            {
+                procI = i;
+                break;
+            }
+        }
+
+        if (Pstream::myProcNo() == procI)
+        {
+            // find corresponding decomposed face triangle
+            label triI = 0;
+            scalar offset = sumTriMagSf_[procI];
+            forAllReverse(triCumulativeMagSf_, i)
+            {
+                if (areaFraction > triCumulativeMagSf_[i] + offset)
+                {
+                    triI = i;
+                    break;
+                }
+            }
+
+            // set cellOwner
+            label faceI = triToFace_[triI];
+            cellOwner = cellOwners_[faceI];
+
+            // find random point in triangle
+            const polyPatch& patch = mesh.boundaryMesh()[patchId_];
+            const pointField& points = patch.points();
+            const face& tf = triFace_[triI];
+            const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
+            const point pf(tri.randomPoint(rnd));
+
+            // position perturbed away from face (into domain)
+            const scalar a = rnd.position(0.1, 0.5);
+            const vector& pc = mesh.cellCentres()[cellOwner];
+            const vector d = mag(pf - pc)*patchNormal_[faceI];
+
+            position = pf - a*d;
+
+            // The position is between the face and cell centre, which could
+            // be in any tet of the decomposed cell, so arbitrarily choose the
+            // first face of the cell as the tetFace and the first point after
+            // the base point on the face as the tetPt.  The tracking will pick
+            // the cell consistent with the motion in the first tracking step
+            tetFaceI = mesh.cells()[cellOwner][0];
+            tetPtI = 1;
+        }
+        else
+        {
+            cellOwner = -1;
+            tetFaceI = -1;
+            tetPtI = -1;
+
+            // dummy position
+            position = pTraits<vector>::max;
+        }
+    }
+    else
+    {
+        cellOwner = -1;
+        tetFaceI = -1;
+        tetPtI = -1;
+
+        // dummy position
+        position = pTraits<vector>::max;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/patchInjectionBase.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PatchInjectionBase
+
+Description
+    Base class for patch-based injection models.
+
+    Class handles injecting at a random point adjacent to the patch faces to
+    provide a more stochsatic view of the injection process.  Patch faces are
+    triangulated, and area fractions accumulated.  The fractional areas are
+    then applied to determine across which face a parcel is to be injected.
+
+SourceFiles
+    PatchInjectionBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchInjectionBase_H
+#define patchInjectionBase_H
+
+#include "word.H"
+#include "labelList.H"
+#include "scalarList.H"
+#include "vectorList.H"
+#include "faceList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward class declarations
+class polyMesh;
+class cachedRandom;
+
+/*---------------------------------------------------------------------------*\
+                     Class patchInjectionBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchInjectionBase
+{
+protected:
+
+    // Protected data
+
+        //- Patch name
+        const word patchName_;
+
+        //- Patch ID
+        const label patchId_;
+
+        //- Patch area - total across all processors
+        scalar patchArea_;
+
+        //- Patch face normal directions
+        vectorList patchNormal_;
+
+        //- List of cell labels corresponding to injector positions
+        labelList cellOwners_;
+
+        //- Fraction of injection controlled by this processor
+        scalar fraction_;
+
+        //- Decomposed patch faces as a list of triangles
+        faceList triFace_;
+
+        //- Addressing from per triangle to patch face
+        labelList triToFace_;
+
+        //- Cumulative triangle area per triangle face
+        scalarList triCumulativeMagSf_;
+
+        //- Cumulative area fractions per processor
+        scalarList sumTriMagSf_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh and patch name
+        patchInjectionBase(const polyMesh& mesh, const word& patchName);
+
+        //- Copy constructor
+        patchInjectionBase(const patchInjectionBase& pib);
+
+
+    //- Destructor
+    virtual ~patchInjectionBase();
+
+
+    // Member Functions
+
+        //- Update patch geometry and derived info for injection locations
+        virtual void updateMesh(const polyMesh& mesh);
+
+        //- Set the injection position and owner cell, tetFace and tetPt
+        virtual void setPositionAndCell
+        (
+            const polyMesh& mesh,
+            cachedRandom& rnd,
+            vector& position,
+            label& cellOwner,
+            label& tetFaceI,
+            label& tetPtI
+        );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // end namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //