ENH: Adding patchProbes to named Patch and Qr and Qe to fvDOM

src/sampling/probes/patchProbes.C

@@ -28,7 +28,9 @@ License
 #include "IOmanip.H"
 // For 'nearInfo' helper class only
 #include "directMappedPatchBase.H"
-#include "meshSearch.H"
+//#include "meshSearch.H"
+#include "treeBoundBox.H"
+#include "treeDataFace.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,27 +43,116 @@ namespace Foam
 
 void Foam::patchProbes::findElements(const fvMesh& mesh)
 {
-    elementList_.clear();
-    elementList_.setSize(size());
+
+    const polyBoundaryMesh& bm = mesh.boundaryMesh();
+
+    label patchI = bm.findPatchID(patchName_);
+
+    if (patchI == -1)
+    {
+        FatalErrorIn
+        (
+            " Foam::patchProbes::findElements(const fvMesh&)"
+        )   << " Unknown patch name "
+            << patchName_ << endl
+            << exit(FatalError);
+    }
+
      // All the info for nearest. Construct to miss
     List<directMappedPatchBase::nearInfo> nearest(this->size());
 
-    // Octree based search engine
-    meshSearch meshSearchEngine(mesh, false);
+    const polyPatch& pp = bm[patchI];
 
-    forAll(*this, probeI)
+    if (pp.size() > 0)
     {
-        const vector& sample = operator[](probeI);
-        label faceI = meshSearchEngine.findNearestBoundaryFace(sample);
-        const point& fc = mesh.faceCentres()[faceI];
-        nearest[probeI].first() = pointIndexHit
+        labelList bndFaces(pp.size());
+        forAll(bndFaces, i)
+        {
+            bndFaces[i] =  pp.start() + i;
+        }
+
+        treeBoundBox overallBb(pp.points());
+        Random rndGen(123456);
+        overallBb = overallBb.extend(rndGen, 1E-4);
+        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        const indexedOctree<treeDataFace> boundaryTree
         (
-            true,
-            fc,
-            faceI
+            treeDataFace    // all information needed to search faces
+            (
+                false,                      // do not cache bb
+                mesh,
+                bndFaces                    // patch faces only
+            ),
+            overallBb,                      // overall search domain
+            8,                              // maxLevel
+            10,                             // leafsize
+            3.0                             // duplicity
         );
-        nearest[probeI].second().first() = magSqr(fc-sample);
-        nearest[probeI].second().second() = Pstream::myProcNo();
+
+
+        if (elementList_.empty())
+        {
+            elementList_.setSize(probeLocations().size());
+
+            forAll(probeLocations(), probeI)
+            {
+                const point sample = probeLocations()[probeI];
+
+                scalar span = boundaryTree.bb().mag();
+
+                pointIndexHit info = boundaryTree.findNearest
+                (
+                    sample,
+                    Foam::sqr(span)
+                );
+
+                if (!info.hit())
+                {
+                    info = boundaryTree.findNearest
+                    (
+                        sample,
+                        Foam::sqr(GREAT)
+                    );
+                }
+
+                label faceI = boundaryTree.shapes().faceLabels()[info.index()];
+
+                const label patchi = bm.whichPatch(faceI);
+
+                if (isA<emptyPolyPatch>(bm[patchi]))
+                {
+                    WarningIn
+                    (
+                        " Foam::patchProbes::findElements(const fvMesh&)"
+                    )
+                    << " The sample point: " << sample
+                    << " belongs to " << patchi
+                    << " which is an empty patch. This is not permitted. "
+                    << " This sample will not be included "
+                    << endl;
+                }
+                else
+                {
+                    const point& fc = mesh.faceCentres()[faceI];
+
+                    directMappedPatchBase::nearInfo sampleInfo;
+
+                    sampleInfo.first() = pointIndexHit
+                    (
+                        true,
+                        fc,
+                        faceI
+                    );
+
+                    sampleInfo.second().first() = magSqr(fc-sample);
+                    sampleInfo.second().second() = Pstream::myProcNo();
+
+                    nearest[probeI]= sampleInfo;
+                }
+            }
+        }
     }
 
 
@@ -93,6 +184,11 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
             localI = nearest[sampleI].first().index();
         }
 
+        if (elementList_.empty())
+        {
+             elementList_.setSize(probeLocations().size());
+        }
+
         elementList_[sampleI] = localI;
     }
 }

src/sampling/probes/patchProbes.H

@@ -58,6 +58,11 @@ class patchProbes
 :
     public probes
 {
+    // Private data
+
+        //- Patch name
+        word patchName_;
+
 
     // Private Member Functions
 

src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -221,6 +221,10 @@ updateCoeffs()
                   + emissivity()()[faceI]*physicoChemical::sigma.value()
                   * pow4(Tp[faceI])
                 )/pi;
+
+             // Emmited heat flux from this ray direction
+            ray.Qem().boundaryField()[patchI][faceI] =
+                refValue()[faceI]*(n[faceI] & ray.dAve());
         }
         else
         {
@@ -228,6 +232,10 @@ updateCoeffs()
             valueFraction()[faceI] = 0.0;
             refGrad()[faceI] = 0.0;
             refValue()[faceI] = 0.0; //not used
+
+             // Incident heat flux on this ray direction
+            ray.Qin().boundaryField()[patchI][faceI] =
+                Iw[faceI]*(n[faceI] & ray.dAve());
         }
     }
 

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.C

@@ -82,6 +82,32 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         mesh_,
         dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
     ),
+    Qem_
+    (
+        IOobject
+        (
+            "Qem",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+    ),
+    Qin_
+    (
+        IOobject
+        (
+            "Qin",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+    ),
     a_
     (
         IOobject
@@ -366,13 +392,16 @@ void Foam::radiation::fvDOM::updateG()
 {
     G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
     Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
+    Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
+    Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 
     forAll(IRay_, rayI)
     {
         IRay_[rayI].addIntensity();
         G_ += IRay_[rayI].I()*IRay_[rayI].omega();
-        //Qr_ += IRay_[rayI].Qr();
         Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
+        Qem_.boundaryField() += IRay_[rayI].Qem().boundaryField();
+        Qin_.boundaryField() += IRay_[rayI].Qin().boundaryField();
     }
 }
 

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.H

@@ -83,6 +83,12 @@ class fvDOM
         //- Total radiative heat flux [W/m2]
         volScalarField Qr_;
 
+         //- Emmited radiative heat flux [W/m2]
+        volScalarField Qem_;
+
+        //- Incidet radiative heat flux [W/m2]
+        volScalarField Qin_;
+
         //- Total absorption coefficient [1/m]
         volScalarField a_;
 
@@ -213,6 +219,12 @@ public:
             //- Const access to total radiative heat flux field
             inline const volScalarField& Qr() const;
 
+            //- Const access to incident radiative heat flux field
+            inline const volScalarField& Qin() const;
+
+            //- Const access to emitted radiative heat flux field
+            inline const volScalarField& Qem() const;
+
             //- Const access to black body
             inline const blackBodyEmission& blackBody() const;
 };

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOMI.H

@@ -91,6 +91,16 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::Qr() const
     return Qr_;
 }
 
+inline const Foam::volScalarField& Foam::radiation::fvDOM::Qin() const
+{
+    return Qin_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::fvDOM::Qem() const
+{
+    return Qem_;
+}
 
 inline const Foam::radiation::blackBodyEmission&
 Foam::radiation::fvDOM::blackBody() const

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -81,6 +81,32 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         mesh_,
         dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
     ),
+    Qin_
+    (
+        IOobject
+        (
+            "Qin" + name(rayId),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
+    ),
+    Qem_
+    (
+        IOobject
+        (
+            "Qem" + name(rayId),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
+    ),
     d_(vector::zero),
     dAve_(vector::zero),
     theta_(theta),

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -81,6 +81,12 @@ private:
         //- Total radiative heat flux on boundary
         volScalarField Qr_;
 
+        //- Incident radiative heat flux on boundary
+        volScalarField Qin_;
+
+        //- Emitted radiative heat flux on boundary
+        volScalarField Qem_;
+
         //- Direction
         vector d_;
 
@@ -171,6 +177,18 @@ public:
             //- Return non-const access to the boundary heat flux
             inline volScalarField& Qr();
 
+            //- Return non-const access to the boundary incident heat flux
+            inline volScalarField& Qin();
+
+            //- Return non-const access to the boundary emmited heat flux
+            inline volScalarField& Qem();
+
+            //- Return const access to the boundary incident heat flux
+            inline const volScalarField& Qin() const;
+
+            //- Return const access to the boundary emmited heat flux
+            inline const volScalarField& Qem() const;
+
             //- Return direction
             inline const vector& d() const;
 

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H

@@ -42,6 +42,31 @@ inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qr()
     return Qr_;
 }
 
+inline const Foam::volScalarField& Foam::radiation::
+radiativeIntensityRay::Qin() const
+{
+    return Qin_;
+}
+
+
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qin()
+{
+    return Qin_;
+}
+
+
+inline const Foam::volScalarField& Foam::radiation::
+radiativeIntensityRay::Qem() const
+{
+    return Qem_;
+}
+
+
+inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qem()
+{
+    return Qem_;
+}
+
 
 inline const Foam::vector& Foam::radiation::radiativeIntensityRay::d() const
 {