diff --git a/applications/utilities/preProcessing/createViewFactors/Allwmake b/applications/utilities/preProcessing/createViewFactors/Allwmake
new file mode 100755
index 0000000000000000000000000000000000000000..9de13a5b033475b4a416150305a0e8467984c408
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/Allwmake
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd "${0%/*}" || exit # Run from this directory
+
+ . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments
+
+#------------------------------------------------------------------------------
+
+wmake $targetType viewFactorModels
+wmake $targetType createViewFactors
+
+#------------------------------------------------------------------------------
diff --git a/applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/files b/applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..3f36ae84ca3c63749236d04150032f210a9ed0d2
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/files
@@ -0,0 +1,3 @@
+createViewFactors.C
+
+EXE = $(FOAM_APPBIN)/createViewFactors
diff --git a/applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/options b/applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..71f2b983ab772cd5be2b024557291d0631d3a301
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+ -I../viewFactorModels/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+ -lviewFactorModels \
+ -lfiniteVolume
diff --git a/applications/utilities/preProcessing/createViewFactors/createViewFactors/createViewFactors.C b/applications/utilities/preProcessing/createViewFactors/createViewFactors/createViewFactors.C
new file mode 100644
index 0000000000000000000000000000000000000000..84ad369a99c67a71488945ea16268015ef015a6b
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/createViewFactors/createViewFactors.C
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+ createViewFactors
+
+Group
+ grpPreProcessingUtilities
+
+Description
+ Creates view factors to be used in the view-factor radiation model.
+
+ Operands:
+ \table
+ Operand | Type | Location
+ input | dictionary | \<constant\>/viewFactorsDict
+ input | dictionary | \<constant\>/finalAgglom
+ output | scalarListList | \<constant\>/F
+ output | mapDistribute | \<constant\>/mapDist
+ output | labelListList | \<constant\>/globalFaceFaces
+ output | volScalarField | \<time\>/viewVectorField
+ output | OBJ | allVisibleFaces.obj
+ \endtable
+
+ where the dictionaries mean:
+ \table
+ Dictionary | Description
+ viewFactorsDict | Main-control dictionary
+ finalAgglom | (Optional) Agglomeration addressing (from faceAgglomerate)
+ F | View factors (matrix)
+ mapDist | Map used for parallel running
+ globalFaceFaces | Face addressing
+ viewVectorField | View factors as a volume field
+ allVisibleFaces.obj | The visualisation of the rays
+ \endtable
+
+Usage
+ Minimal example in \c <constant>/viewFactorsDict:
+ \verbatim
+ // Inherited entries
+ raySearchEngine <word>;
+ agglomerate <bool>;
+ nRayPerFace <label>;
+ writeViewFactors <bool>;
+ writeRays <bool>;
+ ...
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ raySearchEngine | Ray search engine type | word | yes | -
+ agglomerate | Flag to agglomeration | bool | yes | -
+ nRayPerFace | Number of rays issued per face | label | yes | -
+ writeViewFactors | Flag to write the view factor field | bool | yes |-
+ writeRays | Flag to write the ray geometry | bool | no | false
+ \endtable
+
+ Options for the \c raySearchEngine entry:
+ \verbatim
+ voxel | Ray search engine discretising space into uniform voxels
+ \endverbatim
+
+ The inherited entries are elaborated in:
+ - \link viewFactorModel.H \endlink
+ - \link raySearchEngine.H \endlink
+
+Note
+
+ - Participating patches must be in the \c vewFactorWall group, i.e. using the
+ \c inGroups entry of the "\<case\>/polyMesh/boundary" file.
+
+ \verbatim
+ myPatch
+ {
+ type wall;
+ inGroups 2(wall viewFactorWall);
+ ...
+ }
+ \endverbatim
+
+ Reads:
+
+ - <constant>/viewFactorsDict : main controls
+ - <constant>/finalAgglom : agglomeration addressing (from faceAgglomerate)
+
+
+ Generates:
+
+ - <constant>/F : view factors (matrix)
+ - <constant>/mapDist : map used for parallel running
+ - <constant>/globalFaceFaces : face addressing
+
+SeeAlso
+- Foam::VF::raySearchEngine
+- Foam::VE::viewFactorModel
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "viewFactorModel.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ #include "addRegionOption.H"
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createNamedMesh.H"
+
+ IOdictionary dict
+ (
+ IOobject
+ (
+ "viewFactorsDict",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ
+ )
+ );
+
+
+ // Calculate the view factors
+ auto modelPtr = VF::viewFactorModel::New(mesh, dict);
+
+ modelPtr->calculate();
+
+ Info<< nl;
+
+ runTime.printExecutionTime(Info);
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
\ No newline at end of file
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/files b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..eddbe68d672cb51ff132a2a792ebb6c1cc206627
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/files
@@ -0,0 +1,12 @@
+raySearchEngine/raySearchEngine/raySearchEngine.C
+raySearchEngine/raySearchEngine/raySearchEngineNew.C
+raySearchEngine/voxel/voxelRaySearchEngine.C
+/* raySearchEngine/nonUniformVoxel/nonUniformVoxelRaySearchEngine.C */
+/* raySearchEngine/allToAll/allToAllRaySearchEngine.C */
+viewFactorModel/viewFactorModel/viewFactorModel.C
+viewFactorModel/viewFactorModel/viewFactorModelNew.C
+viewFactorModel/viewFactor2AI/viewFactor2AI.C
+viewFactorModel/viewFactor2LI/viewFactor2LI.C
+viewFactorModel/viewFactorHottel/viewFactorHottel.C
+
+LIB = $(FOAM_LIBBIN)/libviewFactorModels
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/options b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..ddcc52fdc7bfa96f701e0225c7149bd8e5263eb9
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/surfMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/parallel/distributed/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lsurfMesh \
+ -lmeshTools \
+ -ldistributed \
+ -lradiationModels \
+ -lfileFormats
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.C
new file mode 100644
index 0000000000000000000000000000000000000000..f1a37bbc5f389adc3274f29c3adf730e272c5d76
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.C
@@ -0,0 +1,620 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "raySearchEngine.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "meshTools.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace VF
+{
+ defineTypeNameAndDebug(raySearchEngine, 0);
+ defineRunTimeSelectionTable(raySearchEngine, mesh);
+
+ const label raySearchEngine::maxDynListLength = 1000000000;
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+void Foam::VF::raySearchEngine::check
+(
+ const labelList& nVisibleFaceFaces
+)
+{
+ label nRay = 0;
+ label nFaceMin = labelMax;
+ label nFaceMax = labelMin;
+ for (const label n : nVisibleFaceFaces)
+ {
+ nFaceMin = min(nFaceMin, n);
+ nFaceMax = max(nFaceMax, n);
+ nRay += n;
+ }
+
+ const label nFace = nVisibleFaceFaces.size();
+ const label nGlobalRays = returnReduce(nRay, sumOp<label>());
+
+ if (nGlobalRays == 0)
+ {
+ FatalErrorInFunction
+ << "No rays identified - view factors will not be calculated"
+ << exit(FatalError);
+ }
+
+ const label globalNFacesMin = returnReduce(nFaceMin, minOp<label>());
+ const label globalNFacesMax = returnReduce(nFaceMax, maxOp<label>());
+ const label nGlobalFaces = returnReduce(nFace, sumOp<label>());
+ const scalar avgFaces = nGlobalRays/scalar(nGlobalFaces);
+
+ Info<< "\nRay summary:" << nl
+ << " Number of rays: " << nGlobalRays << nl
+ << " Number of rays-per-face (min, max, average): ("
+ << globalNFacesMin << ", "
+ << globalNFacesMax << ", "
+ << avgFaces << ")" << endl;
+}
+
+
+Foam::label Foam::VF::raySearchEngine::closestPointIndex
+(
+ const point& p0,
+ const List<point>& pts
+)
+{
+ label pointi = -1;
+
+ scalar distSqr = GREAT;
+ forAll(pts, pti)
+ {
+ const scalar d2 = magSqr(pts[pti] - p0);
+ if (d2 < distSqr)
+ {
+ pointi = pti;
+ distSqr = d2;
+ }
+ }
+
+ return pointi;
+}
+
+
+void Foam::VF::raySearchEngine::createAgglomeration(const IOobject& io)
+{
+ Info<< "\nFace agglomeration: active" << nl
+ << " Reading file " << io.name() << endl;
+
+ // Read agglomeration map
+ const labelListIOList finalAgglom(io);
+
+ Info<< " Creating coarse mesh" << nl;
+
+ agglomMeshPtr_.reset
+ (
+ new singleCellFvMesh
+ (
+ IOobject
+ (
+ IOobject::scopedName("agglom", mesh_.name()),
+ mesh_.time().timeName(),
+ mesh_.time(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ finalAgglom
+ )
+ );
+
+ const auto& coarseMesh = agglomMeshPtr_();
+
+
+ // Calculate total number of fine and coarse faces
+
+ nCoarseFace_ = 0;
+ nFace_ = 0;
+
+ const polyBoundaryMesh& finePatches = mesh_.boundaryMesh();
+ const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
+
+ for (const label patchi : patchIDs_)
+ {
+ nCoarseFace_ += coarsePatches[patchi].size();
+ nFace_ += finePatches[patchi].size();
+ }
+
+ Info<< "\nTotal number of coarse faces: "
+ << returnReduce(nCoarseFace_, sumOp<label>())
+ << endl;
+
+ Info<< "\nTotal number of fine faces: "
+ << returnReduce(nFace_, sumOp<label>())
+ << endl;
+
+ // Collect local Cf, Sf, agglom index on coarse mesh
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ DynamicList<point> localCf(nCoarseFace_);
+ DynamicList<vector> localSf(nCoarseFace_);
+ DynamicList<label> localAgg(nCoarseFace_);
+
+ for (const label patchi : patchIDs_)
+ {
+ const labelList& agglom = finalAgglom[patchi];
+
+ if (agglom.empty()) continue;
+
+ label nAgglom = max(agglom) + 1;
+ const labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
+ const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchi];
+
+ const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchi];
+ const vectorField& coarseSf = coarseMesh.Sf().boundaryField()[patchi];
+
+ const polyPatch& pp = finePatches[patchi];
+ patchAreas_[patchi] += sum(coarseMesh.magSf().boundaryField()[patchi]);
+
+ forAll(coarseCf, facei)
+ {
+ const label coarseFacei = coarsePatchFace[facei];
+ const label agglomi = coarseFacei + coarsePatches[patchi].start();
+
+ // Construct single coarse face
+ const labelList& fineFaces = coarseToFine[coarseFacei];
+ uindirectPrimitivePatch cf
+ (
+ UIndirectList<face>(pp, fineFaces),
+ pp.points()
+ );
+
+ // Collect all points (vertices, face centres)
+ const label nFaces = cf.faceCentres().size();
+ const label nPoints = cf.localPoints().size();
+ List<point> allPoints(nFaces + nPoints);
+ SubList<point>(allPoints, nFaces) = cf.faceCentres();
+ SubList<point>(allPoints, nPoints, nFaces) = cf.localPoints();
+
+ // Coarse face centre set to closest point
+ const label pti = closestPointIndex(coarseCf[facei], allPoints);
+
+ if (pti != -1)
+ {
+ localCf.push_back(allPoints[pti]);
+ localSf.push_back(coarseSf[facei]);
+ localAgg.push_back(agglomi);
+ }
+ }
+ }
+
+ Info<< "\nAssembled coarse patch data" << endl;
+
+ // Distribute local coarse Cf and Sf for shooting rays
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ allCf_[Pstream::myProcNo()].transfer(localCf);
+ allSf_[Pstream::myProcNo()].transfer(localSf);
+ allAgg_[Pstream::myProcNo()].transfer(localAgg);
+
+ Pstream::allGatherList(allCf_);
+ Pstream::allGatherList(allSf_);
+ Pstream::allGatherList(allAgg_);
+
+ Pstream::listCombineGather(patchAreas_, plusEqOp<scalar>());
+ Pstream::broadcast(patchAreas_);
+
+ globalNumbering_ = globalIndex(nCoarseFace_);
+}
+
+
+void Foam::VF::raySearchEngine::createGeometry()
+{
+ DynamicList<point> localCf(mesh_.nBoundaryFaces());
+ DynamicList<vector> localSf(mesh_.nBoundaryFaces());
+
+ const auto& pbm = mesh_.boundaryMesh();
+
+ for (const label patchi : patchIDs_)
+ {
+ localCf.push_back(pbm[patchi].faceCentres());
+ localSf.push_back(pbm[patchi].faceAreas());
+
+ patchAreas_[patchi] += sum(mesh_.magSf().boundaryField()[patchi]);
+ }
+
+ Info<< "\nAssembled patch data" << endl;
+
+ nFace_ = localCf.size();
+ nCoarseFace_ = -1;
+
+ allCf_[Pstream::myProcNo()].transfer(localCf);
+ allSf_[Pstream::myProcNo()].transfer(localSf);
+
+ Pstream::allGatherList(allCf_);
+ Pstream::allGatherList(allSf_);
+
+ Pstream::listCombineGather(patchAreas_, plusEqOp<scalar>());
+ Pstream::broadcast(patchAreas_);
+
+ globalNumbering_ = globalIndex(nFace_);
+}
+
+
+void Foam::VF::raySearchEngine::createParallelAddressing
+(
+ labelList& rayEndFace
+) const
+{
+ // Construct compact numbering
+ // - return map from remote to compact indices
+ // (per processor (!= myProcNo) a map from remote index to compact index)
+ // - construct distribute map
+ // - renumber rayEndFace into compact addressing
+
+ DebugInfo << "\nCreating map distribute" << endl;
+
+ List<Map<label>> compactMap(Pstream::nProcs());
+ mapPtr_.reset(new mapDistribute(globalNumbering_, rayEndFace, compactMap));
+
+ DebugInfo << "\nCreating compact-to-global addressing" << endl;
+
+ // Invert compactMap (from processor+localface to compact) to go
+ // from compact to processor+localface (expressed as a globalIndex)
+ compactToGlobal_.resize_nocopy(mapPtr_->constructSize());
+
+ // Local indices first
+ // Note: are not in compactMap
+ for (label i = 0; i < globalNumbering_.localSize(); ++i)
+ {
+ compactToGlobal_[i] = globalNumbering_.toGlobal(i);
+ }
+
+ forAll(compactMap, proci)
+ {
+ const Map<label>& localToCompactMap = compactMap[proci];
+
+ forAllConstIters(localToCompactMap, iter)
+ {
+ compactToGlobal_[*iter] =
+ globalNumbering_.toGlobal(proci, iter.key());
+ }
+ }
+}
+
+
+Foam::coordSystem::cartesian Foam::VF::raySearchEngine::createCoordSystem
+(
+ const point& origin,
+ const vector& dir
+) const
+{
+ vector axis(Zero);
+
+ for (direction d=0; d<3; ++d)
+ {
+ axis = dir^tensor::I.col(d);
+
+ // Remove empty direction for 2D
+ if (mesh_.nSolutionD() == 2)
+ {
+ meshTools::constrainDirection(mesh_, mesh_.solutionD(), axis);
+ }
+
+ if (magSqr(axis) > 0)
+ {
+ axis.normalise();
+ break;
+ }
+ }
+
+ return coordSystem::cartesian(origin, dir, axis);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::VF::raySearchEngine::createHemiPoints
+(
+ const label nRayPerFace
+) const
+{
+ auto themiPts = tmp<pointField>::New(nRayPerFace);
+ auto& hemiPts = themiPts.ref();
+
+ const label nPoints = hemiPts.size();
+
+ if (mesh_.nSolutionD() == 3)
+ {
+ // Point in range -1 < x < 1; -1 < y < 1; 0 < z 1
+
+ forAll(hemiPts, pointi)
+ {
+ const scalar phi = Foam::acos(1 - (pointi + 0.5)/nPoints);
+ const scalar theta =
+ mathematical::pi*(1 + Foam::sqrt(5.0))*(pointi + 0.5);
+
+ hemiPts[pointi] =
+ vector
+ (
+ Foam::cos(theta)*Foam::sin(phi),
+ Foam::sin(theta)*Foam::sin(phi),
+ Foam::cos(phi)
+ );
+ }
+ }
+ else if (mesh_.nSolutionD() == 2)
+ {
+ // Point in range -1 < x < 1; y = 0; 0 < z < 1; _\|/_
+
+ forAll(hemiPts, pointi)
+ {
+ const scalar theta = mathematical::pi*(pointi+0.5)/nPoints;
+ hemiPts[pointi] = vector(Foam::cos(theta), 0, Foam::sin(theta));
+ }
+ }
+
+ return themiPts;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::VF::raySearchEngine::raySearchEngine
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ mesh_(mesh),
+ mapPtr_(nullptr),
+ compactToGlobal_(),
+ globalNumbering_(),
+ patchGroup_(dict.getOrDefault<word>("patchGroup", "viewFactorWall")),
+ patchIDs_(mesh_.boundaryMesh().indices(patchGroup_)),
+ patchAreas_(mesh_.boundaryMesh().nNonProcessor(), Zero),
+ agglomerate_(dict.get<bool>("agglomerate")),
+ agglomMeshPtr_(nullptr),
+ nFace_(-1),
+ nCoarseFace_(-1),
+ allCf_(UPstream::nProcs()),
+ allSf_(UPstream::nProcs()),
+ allAgg_(UPstream::nProcs())
+{
+ Info<< "\nParticipating patches:" << endl;
+
+ forAll(patchIDs_, i)
+ {
+ const label patchi = patchIDs_[i];
+ Info<< " " << i << ": " << mesh_.boundaryMesh()[patchi].name()
+ << endl;
+ }
+
+ const word agglomName(dict.getOrDefault<word>("agglom", "finalAgglom"));
+
+ IOobject agglomIO
+ (
+ agglomName,
+ mesh_.facesInstance(),
+ mesh_,
+ IOobject::MUST_READ
+ );
+
+
+ if (agglomerate_)
+ {
+ // Sets allCf_, allSf_, allAgg_ based on coarse mesh
+ // Sets nFace_, nCoarseFace_
+ createAgglomeration(agglomIO);
+ }
+ else
+ {
+ // Check for presence of finalAgglom - will cause problems in later
+ // calculations with viewFactor radiation model
+ if (agglomIO.typeHeaderOk<labelListIOList>())
+ {
+ WarningInFunction
+ << "Found agglomeration file: " << agglomIO.objectPath() << nl
+ << " This is inconsistent with the view factor calculation "
+ << "and should be removed" << nl << endl;
+ }
+
+ // Sets allCf_, allSf_ based on fine mesh
+ // Sets nFace_; nCoarseFace_ remains unset (-1)
+ createGeometry();
+ }
+
+ globalNumbering_ =
+ nCoarseFace_ == -1 ? globalIndex(nFace_) : globalIndex(nCoarseFace_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::VF::raySearchEngine::correct
+(
+ labelListList& visibleFaceFaces
+) const
+{
+ labelList rayStartFace;
+ labelList rayEndFace;
+ shootRays(rayStartFace, rayEndFace);
+
+ const label nFace = nParticipatingFaces();
+
+ // Calculate number of visible faces from each local start face
+ labelList nVisibleFaceFaces(nFace, Zero);
+ for (const label facei : rayStartFace)
+ {
+ ++nVisibleFaceFaces[facei];
+ }
+
+ check(nVisibleFaceFaces);
+
+ createParallelAddressing(rayEndFace);
+
+ // Set visible face-faces
+
+ // visibleFaceFaces has:
+ // (local face, local viewed face) = compact viewed face
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ visibleFaceFaces.resize_nocopy(nFace);
+ forAll(nVisibleFaceFaces, facei)
+ {
+ visibleFaceFaces[facei].resize_nocopy(nVisibleFaceFaces[facei]);
+ }
+
+ nVisibleFaceFaces = 0;
+ forAll(rayStartFace, i)
+ {
+ const label facei = rayStartFace[i];
+ const label sloti = rayEndFace[i];
+ visibleFaceFaces[facei][nVisibleFaceFaces[facei]++] = sloti;
+ }
+}
+
+
+void Foam::VF::raySearchEngine::compactAddressing
+(
+ const mapDistribute& map,
+ pointField& compactCf,
+ vectorField& compactSf,
+ List<List<vector>>& compactFineSf,
+ List<List<point>>& compactFineCf,
+ DynamicList<List<point>>& compactPoints,
+ DynamicList<label>& compactPatchId
+) const
+{
+ compactCf.resize_nocopy(map.constructSize());
+ compactSf.resize_nocopy(map.constructSize());
+ compactFineSf.resize_nocopy(map.constructSize());
+ compactFineCf.resize_nocopy(map.constructSize());
+ compactPoints.setCapacity(map.constructSize());
+ compactPatchId.setCapacity(map.constructSize());
+
+ // Insert my local values area and centre values
+ if (agglomMeshPtr_)
+ {
+ SubList<vector>(compactSf, nCoarseFace_) = allSf_[Pstream::myProcNo()];
+ SubList<point>(compactCf, nCoarseFace_) = allCf_[Pstream::myProcNo()];
+
+ const auto& coarseMesh = agglomMeshPtr_();
+ const auto& coarsePatches = coarseMesh.boundaryMesh();
+ const auto& coarseFaces = coarseMesh.faces();
+ const auto& coarsePoints = coarseMesh.points();
+
+ const auto& finalAgglom = coarseMesh.patchFaceAgglomeration();
+
+ // Insert my fine local values per coarse face
+ label sloti = 0;
+ for (const label patchi : patchIDs_)
+ {
+ const auto& fineCfp = mesh_.Cf().boundaryField()[patchi];
+ const auto& fineSfp = mesh_.Sf().boundaryField()[patchi];
+ const labelList& agglom = finalAgglom[patchi];
+
+ if (agglom.empty()) continue;
+
+ const label nAgglom = max(agglom) + 1;
+ const labelListList coarseToFine = invertOneToMany(nAgglom, agglom);
+ const labelList& coarsePatchFace =
+ coarseMesh.patchFaceMap()[patchi];
+
+ const label coarseStart = coarsePatches[patchi].start();
+
+ forAll(coarseToFine, coarsei)
+ {
+ compactPatchId.push_back(patchi);
+
+ const vectorField localPoints
+ (
+ coarsePoints,
+ coarseFaces[coarseStart + coarsei]
+ );
+ compactPoints.push_back(localPoints);
+
+ const label coarseFacei = coarsePatchFace[coarsei];
+ const labelList& fineFaces = coarseToFine[coarseFacei];
+
+ List<point>& fineCf = compactFineCf[sloti];
+ fineCf.resize_nocopy(fineFaces.size());
+ fineCf = UIndirectList<point>(fineCfp, fineFaces);
+
+ List<vector>& fineSf = compactFineSf[sloti];
+ fineSf.resize_nocopy(fineFaces.size());
+ fineSf = UIndirectList<vector>(fineSfp, fineFaces);
+
+ ++sloti;
+ }
+ }
+ }
+ else
+ {
+ SubList<vector>(compactSf, nFace_) = allSf_[Pstream::myProcNo()];
+ SubList<point>(compactCf, nFace_) = allCf_[Pstream::myProcNo()];
+
+ const auto& patches = mesh_.boundaryMesh();
+ const faceList& faces = mesh_.faces();
+
+ label sloti = 0;
+
+ for (const label patchi : patchIDs_)
+ {
+ const auto& Sfp = mesh_.Sf().boundaryField()[patchi];
+ const auto& Cfp = mesh_.Cf().boundaryField()[patchi];
+
+ const polyPatch& pp = patches[patchi];
+
+ forAll(pp, facei)
+ {
+ compactPatchId.push_back(patchi);
+
+ const auto& fpts = faces[facei + pp.start()];
+ compactPoints.push_back(List<point>(mesh_.points(), fpts));
+
+ compactFineCf[sloti] = List<point>({Cfp[facei]});
+ compactFineSf[sloti] = List<vector>({Sfp[facei]});
+ ++sloti;
+ }
+ }
+ }
+
+
+ // Do all swapping
+ map.distribute(compactSf);
+ map.distribute(compactCf);
+ map.distribute(compactFineCf);
+ map.distribute(compactFineSf);
+ map.distribute(compactPoints);
+ map.distribute(compactPatchId);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.H
new file mode 100644
index 0000000000000000000000000000000000000000..31a6d6f8d95a08f8130dd51708c65a6b677a8ae7
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.H
@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::VF::raySearchEngine
+
+Description
+ Base class for ray search engines
+
+ Participating patches must be in the \c viewFactorWall group, i.e. using the
+ \c inGroups entry of the "\<case\>/polyMesh/boundary" file.
+
+ \verbatim
+ myPatch
+ {
+ type wall;
+ inGroups 2(wall viewFactorWall);
+ ...
+ }
+ \endverbatim
+
+ Face agglomeration can be employed, created using the \c faceAgglomerate
+ utility. The file name to be read can be user-defined:
+
+ \verbatim
+ // Name of agglomeration file; default = finalAgglom
+ agglom finalAgglom;
+ \endverbatim
+
+SourceFiles
+ raySearchEngine.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vf_raySearchEngine_H
+#define Foam_vf_raySearchEngine_H
+
+#include "cartesianCS.H"
+#include "mapDistribute.H"
+#include "singleCellFvMesh.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace VF
+{
+
+/*---------------------------------------------------------------------------*\
+ Class raySearchEngine Declaration
+\*---------------------------------------------------------------------------*/
+
+class raySearchEngine
+{
+protected:
+
+ // Protected Data
+
+ //- Reference to the mesh
+ const fvMesh& mesh_;
+
+ //- Parallel map
+ mutable autoPtr<mapDistribute> mapPtr_;
+
+ //- Compact to global addressing
+ mutable labelList compactToGlobal_;
+
+ //- Global numbering
+ globalIndex globalNumbering_;
+
+ //- Name of patch group to identify participating patches
+ const word patchGroup_;
+
+ //- List of participating patch IDs
+ labelList patchIDs_;
+
+ //- Patch areas
+ scalarList patchAreas_;
+
+ //- Agglomeration flag
+ bool agglomerate_;
+
+ //- Agglomerated mesh representation
+ autoPtr<singleCellFvMesh> agglomMeshPtr_;
+
+ //- Number of original faces
+ label nFace_;
+
+ //- Number of coarse faces
+ label nCoarseFace_;
+
+ //- List of all face centres per processor
+ List<pointField> allCf_;
+
+ //- List of all face areas per processor
+ List<vectorField> allSf_;
+
+ //- List of all face agglomeration index per processor
+ List<labelField> allAgg_;
+
+
+ // Protected Member Functions
+
+ static void check(const labelList& nVisibleFaceFaces);
+
+ static label closestPointIndex
+ (
+ const point& p0,
+ const List<point>& pts
+ );
+
+ //- Create patch geometry based on the original mesh
+ void createGeometry();
+
+ //- Create parallel addressing - map, compact-to-global
+ void createParallelAddressing(labelList& rayEndFace) const;
+
+ //- Create Cartesian co-ordinate system
+ coordSystem::cartesian createCoordSystem
+ (
+ const point& origin,
+ const vector& dir
+ ) const;
+
+ //- Create patch geometry based on the agglomerated mesh
+ void createAgglomeration(const IOobject& io);
+
+ //- Create a set of points describing a hemisphere
+ // Note: origin is (0 0 0)
+ tmp<pointField> createHemiPoints(const label nRayPerFace) const;
+
+
+public:
+
+ static const label maxDynListLength;
+
+ //- Run-time type information
+ TypeName("raySearchEngine");
+
+ //- Selection table
+ declareRunTimeSelectionTable
+ (
+ autoPtr,
+ raySearchEngine,
+ mesh,
+ (
+ const fvMesh& mesh,
+ const dictionary& dict
+ ),
+ (mesh, dict)
+ );
+
+ //- Selector
+ static autoPtr<raySearchEngine> New
+ (
+ const fvMesh& mesh,
+ const dictionary& dict
+ );
+
+
+ // Generated Methods
+
+ //- No copy construct
+ raySearchEngine(const raySearchEngine&) = delete;
+
+ //- No copy assignment
+ void operator=(const raySearchEngine&) = delete;
+
+
+ //- Constructor
+ raySearchEngine(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ virtual ~raySearchEngine() = default;
+
+
+ // Public Member Functions
+
+ // Access
+
+ //- Reference to the mesh
+ inline const fvMesh& mesh() const noexcept;
+
+ //- Parallel map
+ inline const mapDistribute& map() const;
+
+ //- Compact to global addressing
+ inline const labelList& compactToGlobal() const noexcept;
+
+ //- Global numbering
+ inline const globalIndex& globalNumbering() const noexcept;
+
+ //- List of participating patch IDs
+ inline const labelList& patchIDs() const noexcept;
+
+ //- Patch areas
+ inline const scalarList& patchAreas() const noexcept;
+
+ //- Number of participating faces
+ inline label nParticipatingFaces() const;
+
+ //- List of all face centres per processor
+ inline const List<pointField>& allCf() const noexcept;
+
+ //- List of all face areas per processor
+ inline const List<vectorField>& allSf() const noexcept;
+
+ //- List of all face agglomeration index per processor
+ inline const List<labelField>& allAgg() const noexcept;
+
+
+ // Main calculation functions
+
+ //- Shoot rays; returns lists of ray start and end faces
+ virtual void shootRays
+ (
+ labelList& rayStartFaceOut,
+ labelList& rayEndFaceOut
+ ) const = 0;
+
+ //- Correct
+ virtual void correct(labelListList& visibleFaceFaces) const;
+
+ //- Create compact addressing
+ void compactAddressing
+ (
+ const mapDistribute& map,
+ pointField& compactCf,
+ vectorField& compactSf,
+ List<List<vector>>& compactFineSf,
+ List<List<point>>& compactFineCf,
+ DynamicList<List<point>>& compactPoints,
+ DynamicList<label>& compactPatchId
+ ) const;
+
+ //- Interpolate field
+ template<class Type>
+ void interpolate
+ (
+ GeometricField<Type, fvPatchField, volMesh>& fld,
+ const List<List<Type>>& values
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace VF
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "raySearchEngineI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "raySearchEngineTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineI.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineI.H
new file mode 100644
index 0000000000000000000000000000000000000000..bc410e8b21b9548f582aeb7ee06ec421d00d9a48
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineI.H
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+const Foam::fvMesh& Foam::VF::raySearchEngine::mesh() const noexcept
+{
+ return mesh_;
+}
+
+
+const Foam::mapDistribute& Foam::VF::raySearchEngine::map() const
+{
+ if (!mapPtr_)
+ {
+ FatalErrorInFunction
+ << "mapDistribute has not been set"
+ << abort(FatalError);
+ }
+
+ return mapPtr_;
+}
+
+
+const Foam::labelList&
+Foam::VF::raySearchEngine::compactToGlobal() const noexcept
+{
+ return compactToGlobal_;
+}
+
+
+const Foam::globalIndex&
+Foam::VF::raySearchEngine::globalNumbering() const noexcept
+{
+ return globalNumbering_;
+}
+
+
+const Foam::labelList& Foam::VF::raySearchEngine::patchIDs() const noexcept
+{
+ return patchIDs_;
+}
+
+
+const Foam::scalarList& Foam::VF::raySearchEngine::patchAreas() const noexcept
+{
+ return patchAreas_;
+}
+
+
+Foam::label Foam::VF::raySearchEngine::nParticipatingFaces() const
+{
+ if (nCoarseFace_ == -1) return nFace_;
+ return nCoarseFace_;
+}
+
+
+const Foam::List<Foam::pointField>&
+Foam::VF::raySearchEngine::allCf() const noexcept
+{
+ return allCf_;
+}
+
+
+const Foam::List<Foam::vectorField>&
+Foam::VF::raySearchEngine::allSf() const noexcept
+{
+ return allSf_;
+}
+
+
+const Foam::List<Foam::labelField>&
+Foam::VF::raySearchEngine::allAgg() const noexcept
+{
+ return allAgg_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineNew.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..78af6adc5222ac4e4c5f7e222ef2663dd23523be
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineNew.C
@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "raySearchEngine.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::VF::raySearchEngine> Foam::VF::raySearchEngine::New
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+{
+ const word modelType(dict.get<word>("raySearchEngine"));
+
+ Info<< "Selecting " << typeName << ": " << modelType << endl;
+
+ auto* ctorPtr = meshConstructorTable(modelType);
+
+ if (!ctorPtr)
+ {
+ FatalIOErrorInLookup
+ (
+ dict,
+ typeName,
+ modelType,
+ *meshConstructorTablePtr_
+ ) << exit(FatalIOError);
+ }
+
+ return autoPtr<raySearchEngine>(ctorPtr(mesh, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineTemplates.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..b345451ada9f120c1cf74ea98f846f683ed7e420
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineTemplates.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+void Foam::VF::raySearchEngine::interpolate
+(
+ GeometricField<Type, fvPatchField, volMesh>& fld,
+ const List<List<Type>>& values
+) const
+{
+ label compacti = 0;
+
+ auto& vfbf = fld.boundaryFieldRef();
+
+ if (agglomMeshPtr_)
+ {
+ const auto& coarseMesh = agglomMeshPtr_();
+ const auto& finalAgglom = coarseMesh.patchFaceAgglomeration();
+
+ for (const label patchi : patchIDs_)
+ {
+ const labelList& agglom = finalAgglom[patchi];
+
+ if (agglom.empty()) continue;
+
+ label nAgglom = max(agglom) + 1;
+ const labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
+ const labelList& coarsePatchFace =
+ coarseMesh.patchFaceMap()[patchi];
+
+ forAll(coarseToFine, i)
+ {
+ const label coarseFacei = coarsePatchFace[i];
+ const labelList& fineFaces = coarseToFine[coarseFacei];
+ const Type sumValues = sum(values[compacti]);
+
+ for (const label fineFacei : fineFaces)
+ {
+ vfbf[patchi][fineFacei] = sumValues;
+ }
+
+ ++compacti;
+ }
+ }
+ }
+ else
+ {
+ label compacti = 0;
+ for (const label patchi : patchIDs_)
+ {
+ auto& vfp = vfbf[patchi];
+
+ for (Type& vfi : vfp)
+ {
+ vfi = sum(values[compacti++]);
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.C
new file mode 100644
index 0000000000000000000000000000000000000000..757f1620443708e5e5e3cd47e0ccb6cb5d9aa1bb
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.C
@@ -0,0 +1,948 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "voxelRaySearchEngine.H"
+#include "processorPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace VF
+{
+ defineTypeNameAndDebug(voxel, 0);
+ addToRunTimeSelectionTable(raySearchEngine, voxel, mesh);
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+void Foam::VF::voxel::setTriangulation(const fvMesh& mesh)
+{
+ Info<< "\nCreating triangulated surface" << endl;
+
+ // Storage for surfaceMesh. Size estimate.
+ DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
+ DynamicList<label> globalFaces(mesh.nBoundaryFaces());
+ label nFace = 0;
+
+ const auto& pbm = mesh.boundaryMesh();
+
+ forAll(patchIDs_, i)
+ {
+ const label patchi = patchIDs_[i];
+ const polyPatch& patch = pbm[patchi];
+ const pointField& points = patch.points();
+
+ for (const face& f : patch)
+ {
+ label nTri = 0;
+ faceList triFaces(f.nTriangles(points));
+ f.triangles(points, nTri, triFaces);
+
+ const label globalFacei =
+ globalNumbering_.toGlobal(Pstream::myProcNo(), nFace++);
+
+ for (const face& f : triFaces)
+ {
+ triangles.push_back(labelledTri(f[0], f[1], f[2], i));
+ globalFaces.push_back(globalFacei);
+ }
+ }
+ }
+
+ triToGlobalFace_.transfer(globalFaces);
+
+ Info<< " Total number of triangles: "
+ << returnReduce(triangles.size(), sumOp<label>())
+ << endl;
+
+ triangles.shrink();
+ surface_ = triSurface(triangles, mesh.points());
+ surface_.compactPoints();
+}
+
+
+Foam::labelList Foam::VF::voxel::setFaceVertexHits
+(
+ const fvMesh& mesh,
+ const labelList& patchIDs
+)
+{
+ labelList vertHits(mesh.points().size(), Zero);
+
+ if (mesh.nSolutionD() == 3)
+ {
+ const auto& pbm = mesh.boundaryMesh();
+
+ // Create a new triangulation based on the surface agglomeration
+ for (const label patchI : patchIDs)
+ {
+ const polyPatch& patch = pbm[patchI];
+ for (const face& f : patch)
+ {
+ for (const label fpi : f)
+ {
+ ++vertHits[fpi];
+ }
+ }
+ }
+
+ for (const auto& pp : pbm)
+ {
+ const labelList& meshPts = pp.meshPoints();
+
+ if (pp.size())
+ {
+ if (isA<processorPolyPatch>(pp))
+ {
+ // Add all processor patch points
+ for (const label pi : meshPts)
+ {
+ ++vertHits[pi];
+ }
+ }
+ else
+ {
+ // Add boundary points
+
+ const auto& bndyPts = pp.boundaryPoints();
+
+ for (const label pi : bndyPts)
+ {
+ ++vertHits[meshPts[pi]];
+ }
+ }
+ }
+ }
+ }
+
+ return vertHits;
+}
+
+
+void Foam::VF::voxel::setCoarseTriangulation(const fvMesh& mesh)
+{
+ Info<< "\nCreating triangulated surface" << endl;
+
+
+ // Filter out fine mesh points along coarse mesh faces
+
+
+ // Storage for surfaceMesh. Size estimate.
+ DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
+ DynamicList<label> globalFaces(mesh.nBoundaryFaces());
+ labelList vertHits = setFaceVertexHits(mesh, patchIDs_);
+
+
+ // Only simplifying edges for 3D
+ const label nVertMin = mesh.nSolutionD() == 3 ? 2 : 0;
+
+ label nInvalid = 0;
+ label nFace = 0;
+
+ const auto& pbm = mesh.boundaryMesh();
+
+ for (const label patchi : patchIDs_)
+ {
+ const polyPatch& patch = pbm[patchi];
+ const pointField& points = patch.points();
+
+ for (const face& f : patch)
+ {
+ DynamicList<label> faceVerts;
+ for (const label fpi : f)
+ {
+ if (vertHits[fpi] > nVertMin)
+ {
+ faceVerts.push_back(fpi);
+ }
+ }
+
+ if (faceVerts.size() < 3)
+ {
+ ++nInvalid;
+ continue;
+ }
+
+ label nTri = 0;
+ const face reducedFace(faceVerts);
+ faceList triFaces(reducedFace.nTriangles(points));
+ reducedFace.triangles(points, nTri, triFaces);
+
+ const label globalFacei =
+ globalNumbering_.toGlobal(Pstream::myProcNo(), nFace++);
+
+ for (const face& f : triFaces)
+ {
+ triangles.push_back(labelledTri(f[0], f[1], f[2], patchi));
+ globalFaces.push_back(globalFacei);
+ }
+ }
+ }
+
+ triToGlobalFace_.transfer(globalFaces);
+
+ Info<< " Total number of triangles: "
+ << returnReduce(triangles.size(), sumOp<label>())
+ << "\n Number of invalid (removed) triangles: "
+ << returnReduce(nInvalid, sumOp<label>())
+ << endl;
+
+ triangles.shrink();
+ surface_ = triSurface(triangles, mesh.points());
+ surface_.compactPoints();
+}
+
+
+void Foam::VF::voxel::broadcast()
+{
+ // Every processor has whole surface
+ const globalIndex globalFaceIdx(globalIndex::gatherOnly{}, surface_.size());
+ const globalIndex globalPointIdx
+ (
+ globalIndex::gatherOnly{},
+ surface_.points().size()
+ );
+
+ List<labelledTri> globalSurfaceTris(globalFaceIdx.gather(surface_));
+ pointField globalSurfacePoints(globalPointIdx.gather(surface_.points()));
+ List<label> globalTriToGlobalFace(globalFaceIdx.gather(triToGlobalFace_));
+
+
+ for (const label proci : globalPointIdx.allProcs())
+ {
+ const label offset = globalPointIdx.localStart(proci);
+
+ if (offset)
+ {
+ for
+ (
+ labelledTri& tri
+ : globalSurfaceTris.slice(globalFaceIdx.range(proci))
+ )
+ {
+ tri[0] += offset;
+ tri[1] += offset;
+ tri[2] += offset;
+ }
+ }
+ }
+
+ surface_ =
+ triSurface
+ (
+ std::move(globalSurfaceTris),
+ std::move(globalSurfacePoints)
+ );
+
+ Pstream::broadcast(surface_);
+
+ triToGlobalFace_ = std::move(globalTriToGlobalFace);
+ Pstream::broadcast(triToGlobalFace_);
+}
+
+
+Foam::pointHit Foam::VF::voxel::rayTriIntersect
+(
+ const label trii,
+ const point& origin,
+ const vector& direction
+) const
+{
+ // Note: origin was made relative to voxel mesh on entry to hit function
+ // - need to convert back into global position to be consistent with
+ // triangles for intersection tests
+
+ const auto& ind = surface_[trii];
+ const auto& pts = surface_.points();
+
+ // Note: flipped orientation of triangle (into domain) so that we can use
+ // visibility check when performing ray-triangle intersections
+ const triPointRef tri(pts[ind[0]], pts[ind[2]], pts[ind[1]]);
+
+ static scalar tolerance = 1e-6;
+
+ return
+ tri.intersection
+ (
+ globalPosition(origin),
+ direction,
+ intersection::VISIBLE,
+ tolerance
+ );
+}
+
+
+void Foam::VF::voxel::writeBox
+(
+ OBJstream& os,
+ bool lines,
+ const boundBox& bb
+) const
+{
+ os.write(treeBoundBox(bb), lines);
+}
+
+
+void Foam::VF::voxel::writeVoxels(const word& fName) const
+{
+ if (!UPstream::master()) return;
+
+ OBJstream os(fName);
+ Info<< "Writing voxels to " << os.name() << endl;
+
+ boundBox bb;
+ const bool lines = true;
+ for (label i = 0; i < nijk_[0]; ++i)
+ {
+ for (label j = 0; j < nijk_[1]; ++j)
+ {
+ for (label k = 0; k < nijk_[2]; ++k)
+ {
+ bb.min() = voxelMin(i, j, k);
+ bb.max() = voxelMax(i, j, k);
+ writeBox(os, lines, bb);
+ }
+ }
+ }
+
+ Info<< "- done" << endl;
+}
+
+
+void Foam::VF::voxel::writeTriBoundBoxes(const word& fName) const
+{
+ if (!UPstream::master()) return;
+
+ OBJstream os(fName);
+ Info<< "Writing triangle boundBoxes to " << os.name() << endl;
+
+ const bool lines = true;
+ for (const auto& voxeli : objects_)
+ {
+ for (const label trii : voxeli)
+ {
+ writeBox(os, lines, objectBbs_[trii]);
+ }
+ }
+
+ Info<< "- done" << endl;
+}
+
+
+void Foam::VF::voxel::writeHitObjects
+(
+ const label voxeli,
+ const point& origin,
+ const vector& dir
+) const
+{
+ OBJstream os("voxel_" + Foam::name(voxeli) + ".obj");
+
+ // Write voxel
+ labelVector ijk = this->ijk(voxeli);
+
+ boundBox voxelBb;
+ voxelBb.min() = voxelMin(ijk[0], ijk[1], ijk[2]);
+ voxelBb.max() = voxelMax(ijk[0], ijk[1], ijk[2]);
+
+ writeBox(os, true, voxelBb);
+
+ for (const auto& trii : objects_[voxeli])
+ {
+ writeBox(os, true, objectBbs_[trii]);
+
+ const auto& ind = surface_[trii];
+ const auto& pts = surface_.points();
+ const triPointRef tri(pts[ind[0]], pts[ind[1]], pts[ind[2]]);
+ os.write(tri, false);
+ }
+}
+
+
+Foam::pointIndexHit Foam::VF::voxel::hitObject
+(
+ const label voxeli,
+ const point& origin,
+ const vector& dir,
+ const vector& t,
+ const scalar minDistance
+) const
+{
+ if (objects_[voxeli].empty()) return pointIndexHit();
+
+ // boundBox rayBb;
+ // rayBb.add(origin);
+ // rayBb.add(origin + dir*(dir & t));
+
+ label triHiti = -1;
+ //rayBb.add(origin + dir);
+ //rayBb.inflate(0.01);
+
+ if (debug > 2)
+ {
+ writeHitObjects(voxeli, origin, dir);
+ }
+
+ // Determine all triangles that intersect with ray
+ // - only keep nearest hit
+
+ pointHit closestHit;
+ for (const auto& trii : objects_[voxeli])
+ {
+ // Only perform ray/tri intersection if bound boxes intersect
+ //if (objectBbs_[trii].overlaps(rayBb))
+ {
+ pointHit pi = rayTriIntersect(trii, origin, dir);
+
+ if
+ (
+ pi.hit()
+ && (
+ pi.distance() > minDistance
+ && pi.distance() < closestHit.distance()
+ )
+ )
+ {
+ triHiti = trii;
+ closestHit = pi;
+ }
+ }
+ }
+
+ return pointIndexHit(closestHit, triHiti);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::VF::voxel::voxel(const fvMesh& mesh, const dictionary& dict)
+:
+ raySearchEngine(mesh, dict),
+ bb0_(),
+ span0_(Zero),
+ nijk_(Zero),
+ dxyz_(Zero),
+ nRayPerFace_(dict.get<label>("nRayPerFace")),
+ nTriPerVoxelMax_(dict.getOrDefault<label>("nTriPerVoxelMax", 50)),
+ depthMax_(dict.getOrDefault<label>("depthMax", 5)),
+ objects_(),
+ objectBbs_()
+{
+ if (agglomMeshPtr_)
+ {
+ setCoarseTriangulation(*agglomMeshPtr_);
+ }
+ else
+ {
+ setTriangulation(mesh);
+ }
+
+ broadcast();
+
+ objectBbs_.resize_nocopy(surface_.size());
+
+ bb0_.add(surface_.points());
+ bb0_.inflate(0.01);
+ span0_ = bb0_.span();
+
+ {
+ scalar maxDx = span0_.x();
+ scalar maxDy = span0_.y();
+ scalar maxDz = span0_.z();
+
+ const label refDn = 8;
+ scalar maxDim = max(maxDx, max(maxDy, maxDz));
+
+ setVoxelDims
+ (
+ refDn*maxDx/maxDim,
+ refDn*maxDy/maxDim,
+ refDn*maxDz/maxDim
+ );
+
+ objects_.resize_nocopy(nVoxel());
+ }
+
+ label depth = 0;
+ label trii = 0;
+ voxelise(objects_, trii, depth);
+
+ Info<< "\nCreated voxel mesh: " << nijk_ << endl;
+
+ if ((debug > 3) && UPstream::master())
+ {
+ writeVoxels("voxels.obj");
+ writeTriBoundBoxes("triBoundBoxes.obj");
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::VF::voxel::refineObjects
+(
+ List<DynamicList<label>>& objects,
+ const label triMax
+)
+{
+ refineVoxelDims();
+
+ if (debug > 2) Pout<< "Refining voxels: n=" << nijk_ << endl;
+
+ List<DynamicList<label>> objectsNew(objects.size()*8);
+
+ for (label trii = 0; trii <= triMax; ++trii)
+ {
+ addBbToVoxels(objectBbs_[trii], trii, objectsNew);
+ }
+
+ objects.transfer(objectsNew);
+}
+
+
+Foam::label Foam::VF::voxel::addBbToVoxels
+(
+ const boundBox& bb,
+ const label trii,
+ List<DynamicList<label>>& objects
+) const
+{
+ //Pout<< "addBbToVoxels: trii=" << trii << " n=" << nijk_ << endl;
+
+ const point minbb(localPosition(bb.min()));
+ const label i0 = max(0, floor(minbb.x()/dxyz_[0]));
+ const label j0 = max(0, floor(minbb.y()/dxyz_[1]));
+ const label k0 = max(0, floor(minbb.z()/dxyz_[2]));
+
+ const point maxbb(localPosition(bb.max()));
+ const label i1 = min(nijk_[0], ceil(maxbb.x()/dxyz_[0]));
+ const label j1 = min(nijk_[1], ceil(maxbb.y()/dxyz_[1]));
+ const label k1 = min(nijk_[2], ceil(maxbb.z()/dxyz_[2]));
+
+ label nTriMax = 0;
+
+ for (label ii = i0; ii < i1; ++ii)
+ {
+ for (label jj = j0; jj < j1; ++jj)
+ {
+ for (label kk = k0; kk < k1; ++kk)
+ {
+ const label voxeli = this->voxeli(ii, jj, kk);
+
+ objects[voxeli].push_back(trii);
+ nTriMax = max(nTriMax, objects[voxeli].size());
+ }
+ }
+ }
+
+ return nTriMax;
+}
+
+
+void Foam::VF::voxel::voxelise
+(
+ List<DynamicList<label>>& objects,
+ const label trii0,
+ const label depth
+)
+{
+ if (debug > 2)
+ {
+ Pout<< "voxelise - start at tri=" << trii0
+ << " depth=" << depth
+ << endl;
+ }
+
+ const auto& points = surface_.points();
+
+ for (label trii = trii0; trii < surface_.size(); ++trii)
+ {
+ // Triangle bounding box
+ boundBox bb(points, surface_[trii]);
+ bb.inflate(0.01);
+ objectBbs_[trii] = bb;
+
+ const label nVoxelMax = addBbToVoxels(bb, trii, objects);
+
+ // Number of triangles per voxel - if exceed limit, refine voxels...
+ if (nVoxelMax > nTriPerVoxelMax_ && depth < depthMax_)
+ {
+ refineObjects(objects, trii);
+ voxelise(objects, trii + 1, depth + 1);
+ break;
+ }
+ }
+}
+
+
+Foam::pointIndexHit Foam::VF::voxel::hit
+(
+ const label tri0,
+ const vector& direction0
+) const
+{
+ if (tri0 > surface_.size()-1)
+ {
+ FatalErrorInFunction
+ << "Index tri0 out of bounds: " << tri0
+ << ". Surface size: " << surface_.size()
+ << abort(FatalError);
+ }
+
+ return hit(surface_.faceCentres()[tri0], direction0);
+}
+
+
+Foam::pointIndexHit Foam::VF::voxel::hit
+(
+ const point& origin0,
+ const vector& direction0
+) const
+{
+ // Initialise return value
+ pointIndexHit hitInfo;
+
+ const point origin(localPosition(origin0));
+
+ if (cmptMin(origin) < 0)
+ {
+ FatalErrorInFunction
+ << "Point located outside voxel mesh"
+ << " - possible coarsening problem?"
+ << abort(FatalError);
+ }
+
+ if (magSqr(direction0) < ROOTVSMALL)
+ {
+ WarningInFunction
+ << "Supplied direction has zero size"
+ << endl;
+
+ return hitInfo;
+ }
+
+ const vector dir(normalised(direction0));
+
+ labelVector ijk(Zero);
+ labelVector step(Zero);
+ vector tDelta(vector::max);
+ vector tMax(vector::max);
+
+ for (direction d = 0; d < 3; ++d)
+ {
+ ijk[d] = floor(origin[d]/dxyz_[d]);
+ step[d] = sign0(dir[d]);
+ if (step[d] != 0)
+ {
+ tDelta[d] = mag(dxyz_[d]/dir[d]);
+
+ scalar voxelMax = (1 + ijk[d] - neg(dir[d]))*dxyz_[d];
+ tMax[d] = (voxelMax - origin[d])/dir[d];
+ }
+ }
+
+ if (debug > 2)
+ {
+ Info<< "surfBb:" << boundBox(surface_.points())
+ << " bb:" << bb0_
+ << " origin" << origin0
+ << " voxel_origin:" << origin
+ << " ijk:" << ijk
+ << " step:" << step
+ << " dxyz:" << dxyz_
+ << " tDelta:" << tDelta
+ << " tMax:" << tMax
+ << endl;
+ }
+
+ auto traverse = [&](const label i)
+ {
+ ijk[i] += step[i];
+ if (outOfBounds(ijk, i)) return false;
+ tMax[i] += tDelta[i];
+ return true;
+ };
+
+
+ while (true)
+ {
+ const label voxeli = this->voxeli(ijk);
+
+ if (debug > 2)
+ {
+ Info<< "ijk:" << ijk
+ << " voxeli:" << voxeli
+ << " t:" << tMax
+ << " objs:" << objects_[voxeli].size()
+ << endl;
+ }
+
+ hitInfo = hitObject(voxeli, origin, dir, tMax);
+
+ if (hitInfo.hit())
+ {
+ // Valid hit
+ break;
+ }
+ else
+ {
+ if (tMax[0] < tMax[1] && tMax[0] < tMax[2])
+ {
+ if (!traverse(0)) break;
+ }
+ else if (tMax[1] < tMax[2])
+ {
+ if (!traverse(1)) break;
+ }
+ else
+ {
+ if (!traverse(2)) break;
+ }
+ }
+ }
+
+ return hitInfo;
+}
+
+
+void Foam::VF::voxel::shootRays
+(
+ labelList& rayStartFaceOut,
+ labelList& rayEndFaceOut
+) const
+{
+ (void)mesh_.time().cpuTimeIncrement();
+
+ const pointField& myFc = allCf_[Pstream::myProcNo()];
+ const vectorField& myArea = allSf_[Pstream::myProcNo()];
+
+ const label nTotalRays = myFc.size()*nRayPerFace_;
+ if (nTotalRays > maxDynListLength)
+ {
+ FatalErrorInFunction
+ << "Dynamic list needs more capacity to support " << nRayPerFace_
+ << " rays per face (" << nTotalRays << "). "
+ << "Limit set by maxDynListLength: " << maxDynListLength
+ << abort(FatalError);
+ }
+
+ DynamicList<label> rayStartFace(nTotalRays);
+ DynamicList<label> rayEndFace(rayStartFace.size());
+
+ DynamicList<label> endFacei(nTotalRays);
+ DynamicList<label> startFacei(nTotalRays);
+
+ const pointField hemi(createHemiPoints(nRayPerFace_));
+
+ Info<< "\nShooting rays" << endl;
+
+ const scalar nDiv = 10;
+ const label nStep = floor(myFc.size()/nDiv);
+ labelHashSet localFaceHits;
+
+ for (label stepi=0; stepi<nDiv; ++stepi)
+ {
+ const label step0 = stepi*nStep;
+ const label step1 = stepi == nDiv - 1 ? myFc.size() : (stepi+1)*nStep;
+
+ for (label i = step0; i < step1; ++i)
+ {
+ // Info<< "i/N = " << i << "/" << (myFc.size()-1)
+ // << " step0:" << step0 << " step1:" << step1
+ // << endl;
+
+ localFaceHits.clear();
+
+ const point& origin = myFc[i];
+ const vector dir(-normalised(myArea[i]));
+
+ const coordSystem::cartesian cs = createCoordSystem(origin, dir);
+
+ const vectorField pts(cs.transformPoint(hemi));
+
+ for (const auto& p : pts)
+ {
+ const pointIndexHit hitInfo = hit(origin, p-origin);
+
+ if (hitInfo.hit())
+ {
+ label hitFacei = triToGlobalFace_[hitInfo.index()];
+
+ if (localFaceHits.insert(hitFacei))
+ {
+ endFacei.push_back(hitFacei);
+ startFacei.push_back(i);
+ }
+ }
+ else
+ {
+ // Miss
+ }
+ }
+ }
+
+ const label globalStepi = returnReduce(stepi, minOp<label>()) + 1;
+
+ Info<< " " << globalStepi/nDiv*100 << "% complete" << endl;
+ }
+
+
+ // Ensure all rays are unique/filter out duplicates
+ // - add symmetric connections for non-self-intersecting rays
+
+ if (UPstream::parRun())
+ {
+ edgeHashSet uniqueRays;
+ List<DynamicList<label>> remoteStartFace(Pstream::nProcs());
+ List<DynamicList<label>> remoteEndFace(Pstream::nProcs());
+
+ const labelList globalStartFaceIDs
+ (
+ globalNumbering_.toGlobal(startFacei)
+ );
+
+ forAll(globalStartFaceIDs, rayi)
+ {
+ label start = globalStartFaceIDs[rayi];
+ label end = endFacei[rayi];
+
+ const label endProci = globalNumbering_.whichProcID(end);
+
+ if (start > end) Swap(start, end);
+
+ if (uniqueRays.insert(edge(start, end)))
+ {
+ if (endProci != UPstream::myProcNo())
+ {
+ // Convert local face into global face and vice-versa
+ remoteStartFace[endProci].push_back(start);
+ remoteEndFace[endProci].push_back(end);
+ }
+ }
+ }
+
+ PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+ // Send remote data
+ for (const int domain : Pstream::allProcs())
+ {
+ if (domain != Pstream::myProcNo())
+ {
+ UOPstream str(domain, pBufs);
+ str << remoteStartFace[domain]
+ << remoteEndFace[domain];
+ }
+ }
+
+ pBufs.finishedSends();
+
+ // Consume
+ for (const int domain : Pstream::allProcs())
+ {
+ if (domain != Pstream::myProcNo())
+ {
+ UIPstream is(domain, pBufs);
+ is >> remoteStartFace[domain]
+ >> remoteEndFace[domain];
+
+ forAll(remoteStartFace[domain], i)
+ {
+ const label start = remoteStartFace[domain][i];
+ const label end = remoteEndFace[domain][i];
+ uniqueRays.insert(edge(start, end));
+ }
+ }
+ }
+
+ // Populate ray addressing based on uniqueRays
+ for (const edge& e : uniqueRays)
+ {
+ const label start = e.first();
+ const label end = e.second();
+
+ bool sameFace = start == end;
+
+ if (globalNumbering_.isLocal(start))
+ {
+ // Ray originates from this processor
+ const label localStart = globalNumbering_.toLocal(start);
+
+ rayStartFace.append(localStart);
+ rayEndFace.append(end);
+ }
+
+ if (!sameFace && globalNumbering_.isLocal(end))
+ {
+ // Ray hits this processor
+ // - add symmetric ray from end->start
+ const label localEnd = globalNumbering_.toLocal(end);
+
+ rayStartFace.append(localEnd);
+ rayEndFace.append(start);
+ }
+ }
+ }
+ else
+ {
+ // Populate ray addressing based on uniqueRays
+
+ edgeHashSet uniqueRays;
+
+ forAll(startFacei, rayi)
+ {
+ label start = startFacei[rayi];
+ label end = endFacei[rayi];
+
+ if (start > end) Swap(start, end);
+
+ if (uniqueRays.insert(edge(start, end)))
+ {
+ rayStartFace.append(start);
+ rayEndFace.append(end);
+
+ if (start != end)
+ {
+ // Add symmetric ray from end->start
+ rayStartFace.append(end);
+ rayEndFace.append(start);
+ }
+ }
+ }
+ }
+
+ rayStartFaceOut.transfer(rayStartFace);
+ rayEndFaceOut.transfer(rayEndFace);
+
+ Info<< " Time taken: " << mesh_.time().cpuTimeIncrement() << " s"
+ << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.H
new file mode 100644
index 0000000000000000000000000000000000000000..1103e1a881ceea25283207420c60c96abd700720
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.H
@@ -0,0 +1,284 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::VF::voxel
+
+Description
+ Ray search engine based on discretising space into uniform voxels
+
+ Voxels are refined using 2x2x2 refinement.
+
+ The number of rays per face is supplied by the user, whereby rays are
+ issued uniformly across a hemisphere.
+
+Usage
+ Minimal example by using \c <constant>/viewFactorsDict:
+ \verbatim
+ // Mandatory entries
+ nRayPerFace <label>;
+
+ // Optional entries
+ nTriPerVoxelMax <label>;
+ depthMax <label>;
+
+ // Inherited entries
+ ...
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ nRayPerFace | Number of rays issued per face | label | yes | -
+ nRayPerFace | Target number of triangles per voxel | label | no | 50
+ depthMax | Maximum voxel refinement depth | label | no | 5
+ \endtable
+
+ The inherited entries are elaborated in:
+ - \link raySearchEngine.H \endlink
+
+SourceFiles
+ voxelRaySearchEngine.C
+
+SeeAlso
+- Foam::VF::raySearchEngine
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vf_voxelRaySearchEngine_H
+#define Foam_vf_voxelRaySearchEngine_H
+
+#include "raySearchEngine.H"
+#include "labelVector.H"
+#include "OBJstream.H"
+#include "pointIndexHit.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace VF
+{
+
+/*---------------------------------------------------------------------------*\
+ Class voxel Declaration
+\*---------------------------------------------------------------------------*/
+
+class voxel
+:
+ public raySearchEngine
+{
+ // Private Data
+
+ //- Triangulation of view factor patches
+ triSurface surface_;
+
+ //- Triangle to mesh face addressing
+ labelList triToGlobalFace_;
+
+ //- Surface bounding box
+ boundBox bb0_;
+
+ //- Span of surface bounding box
+ vector span0_;
+
+ //- Voxel discretisation
+ labelVector nijk_;
+
+ //- Voxel dimensions
+ vector dxyz_;
+
+ //- Number of rays issued per face
+ const label nRayPerFace_;
+
+ //- Maximum number of triangles per voxel (trigger to refine voxels)
+ const label nTriPerVoxelMax_;
+
+ //- Maximum depth of voxel refinement. Note: voxels are refined 2x2x2
+ const label depthMax_;
+
+ //- List of triangle per voxel
+ List<DynamicList<label>> objects_;
+
+ //- List of triangle bounding boxes
+ List<boundBox> objectBbs_;
+
+
+ // Private Member Functions
+
+ inline bool outOfBounds(const labelVector& ijk, const label dir) const;
+
+ inline point localPosition(const vector& globalPosition) const;
+
+ inline point globalPosition(const vector& localPosition) const;
+
+ inline void setVoxelDims(const label i, const label j, const label k);
+
+ inline void refineVoxelDims();
+
+ inline point voxelMin
+ (
+ const label i,
+ const label j,
+ const label k
+ ) const;
+
+ inline point voxelMax
+ (
+ const label i,
+ const label j,
+ const label k
+ ) const;
+
+ inline constexpr label sign0(const scalar x) const;
+
+ //- Set triangulation based on original mesh
+ void setTriangulation(const fvMesh& mesh);
+
+ //- Set the participating face vertices when simplifying edges
+ static labelList setFaceVertexHits
+ (
+ const fvMesh& mesh,
+ const labelList& patchIDs
+ );
+
+ //- Set triangulation based on agglomeration
+ void setCoarseTriangulation(const fvMesh& mesh);
+
+ //- Broadcast triangulation across all procs
+ void broadcast();
+
+ void refineObjects
+ (
+ List<DynamicList<label>>& objects,
+ const label triMax
+ );
+
+ label addBbToVoxels
+ (
+ const boundBox& bb,
+ const label trii,
+ List<DynamicList<label>>& objects
+ ) const;
+
+ void voxelise
+ (
+ List<DynamicList<label>>& objects,
+ const label trii0,
+ const label depth
+ );
+
+ pointHit rayTriIntersect
+ (
+ const label trii,
+ const point& origin,
+ const vector& direction
+ ) const;
+
+ pointIndexHit hitObject
+ (
+ const label voxeli,
+ const point& origin,
+ const vector& direction,
+ const vector& t,
+ const scalar minDistance = 1e-6
+ ) const;
+
+ void writeBox
+ (
+ OBJstream& os,
+ bool lines,
+ const boundBox& bb
+ ) const;
+
+ void writeVoxels(const word& fName) const;
+
+ void writeTriBoundBoxes(const word& fName) const;
+
+ void writeHitObjects
+ (
+ const label voxeli,
+ const point& origin,
+ const vector& dir
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("voxel");
+
+ //- Constructor
+ voxel(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ virtual ~voxel() = default;
+
+
+ // Public Member Functions
+
+ inline labelVector nijk() const noexcept;
+
+ inline label nVoxel() const noexcept;
+
+ inline label voxeli(const labelVector ijk) const noexcept;
+
+ inline label voxeli
+ (
+ const label i,
+ const label j,
+ const label k
+ ) const noexcept;
+
+ inline labelVector ijk(const label voxeli) const noexcept;
+
+ pointIndexHit hit(const point& origin, const vector& dir) const;
+
+ pointIndexHit hit(const label tri0, const vector& dir) const;
+
+ virtual void shootRays
+ (
+ labelList& rayStartFaceOut,
+ labelList& rayEndFaceOut
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace VF
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "voxelRaySearchEngineI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngineI.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngineI.H
new file mode 100644
index 0000000000000000000000000000000000000000..a46ca7e1da5fae77d5a8f03cae70ba02a742e945
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngineI.H
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+bool Foam::VF::voxel::outOfBounds
+(
+ const labelVector& ijk,
+ const label dir
+) const
+{
+ return (ijk[dir] < 0 || ijk[dir] >= nijk_[dir]);
+};
+
+
+Foam::point Foam::VF::voxel::localPosition(const vector& globalPosition) const
+{
+ return globalPosition - bb0_.min();
+}
+
+
+Foam::point Foam::VF::voxel::globalPosition(const vector& localPosition) const
+{
+ return bb0_.min() + localPosition;
+}
+
+
+void Foam::VF::voxel::setVoxelDims(const label i, const label j, const label k)
+{
+ nijk_[0] = max(1, i);
+ nijk_[1] = max(1, j);
+ nijk_[2] = max(1, k);
+
+ dxyz_[0] = span0_[0]/nijk_[0];
+ dxyz_[1] = span0_[1]/nijk_[1];
+ dxyz_[2] = span0_[2]/nijk_[2];
+}
+
+
+void Foam::VF::voxel::refineVoxelDims()
+{
+ nijk_ *= 2;
+
+ // Do not refine empty direction for 2D
+ const auto& solutionD = mesh_.solutionD();
+ for (direction d=0; d<3; ++d)
+ {
+ if (solutionD[d] == -1)
+ {
+ nijk_[d] = 1;
+ }
+ }
+
+ setVoxelDims(nijk_[0], nijk_[1], nijk_[2]);
+}
+
+
+Foam::point Foam::VF::voxel::voxelMin
+(
+ const label i,
+ const label j,
+ const label k
+) const
+{
+ return point(i*dxyz_[0], j*dxyz_[1], k*dxyz_[2]);
+}
+
+
+Foam::point Foam::VF::voxel::voxelMax
+(
+ const label i,
+ const label j,
+ const label k
+) const
+{
+ return voxelMin(i+1, j+1, k+1);
+}
+
+
+constexpr Foam::label Foam::VF::voxel::sign0(const scalar x) const
+{
+ if (x > 0) return 1;
+ if (x < 0) return -1;
+ return 0;
+};
+
+
+Foam::labelVector Foam::VF::voxel::nijk() const noexcept
+{
+ return nijk_;
+}
+
+
+Foam::label Foam::VF::voxel::nVoxel() const noexcept
+{
+ return nijk_[0]*nijk_[1]*nijk_[2];
+}
+
+
+Foam::label Foam::VF::voxel::voxeli
+(
+ const labelVector ijk
+) const noexcept
+{
+ return voxeli(ijk[0], ijk[1], ijk[2]);
+}
+
+
+Foam::label Foam::VF::voxel::voxeli
+(
+ const label i,
+ const label j,
+ const label k
+) const noexcept
+{
+ return i + (nijk_[0]*(j + (nijk_[1]*k)));
+}
+
+
+Foam::labelVector Foam::VF::voxel::ijk(const label voxeli) const noexcept
+{
+ const label nx = nijk_[0];
+ const label ny = nijk_[1];
+
+ const label i = voxeli%nx;
+ const label k = voxeli/nx/ny;
+ const label j = (voxeli/nx)%ny;
+
+ return labelVector(i, j, k);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.C
new file mode 100644
index 0000000000000000000000000000000000000000..2b7afe34bff644309a57a8532eaa415eb71ca37d
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.C
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "viewFactor2AI.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace VF
+{
+ defineTypeNameAndDebug(viewFactor2AI, 0);
+ addToRunTimeSelectionTable(viewFactorModel, viewFactor2AI, mesh);
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+Foam::scalar Foam::VF::viewFactor2AI::calculateFij
+(
+ const point& xi,
+ const point& xj,
+ const vector& dAi,
+ const vector& dAj
+)
+{
+ const vector r(xj - xi);
+ const scalar rMag = mag(r);
+ const scalar dAiMag(mag(dAi));
+ const scalar dAjMag(mag(dAj));
+
+ if (rMag > SMALL && dAiMag > SMALL && dAjMag > SMALL)
+ {
+ const vector nr(r/rMag);
+ const vector ni(dAi/dAiMag);
+ const vector nj(dAj/dAjMag);
+
+ const scalar Fij =
+ -(nr & ni)*(nr & nj)/sqr(rMag)*dAjMag*dAiMag/mathematical::pi;
+
+ if (Fij > 0) return Fij;
+ }
+
+ return 0;
+}
+
+
+Foam::scalarListList Foam::VF::viewFactor2AI::calculate
+(
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCf,
+ const vectorField& compactSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+) const
+{
+ // Fill local view factor matrix
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ scalarListList Fij(visibleFaceFaces.size());
+
+ forAll(visibleFaceFaces, facei)
+ {
+ if (debug > 1)
+ {
+ Pout<< "facei:" << facei << "/" << visibleFaceFaces.size()
+ << endl;
+ }
+
+ const labelList& visibleFaces = visibleFaceFaces[facei];
+
+ Fij[facei].resize(visibleFaces.size(), Zero);
+
+ const point& dCi = compactCf[facei];
+ const vector& Ai = compactSf[facei];
+ const scalar magAi = mag(Ai);
+
+ if (magAi < SMALL) continue;
+
+ forAll(visibleFaces, visibleFacei)
+ {
+ const label sloti = visibleFaces[visibleFacei];
+ const point& dCj = compactCf[sloti];
+ const vector& Aj = compactSf[sloti];
+
+ const scalar Fij2AI = calculateFij(dCi, dCj, Ai, Aj);
+
+ Fij[facei][visibleFacei] = Fij2AI/magAi;
+ }
+ }
+
+
+ return Fij;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::VF::viewFactor2AI::viewFactor2AI
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ viewFactorModel(mesh, dict)
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.H
new file mode 100644
index 0000000000000000000000000000000000000000..9ba377435b4d55193731e0d6ef6b1012e01888ca
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.H
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::VF::viewFactor2AI
+
+Description
+ Computes view factors according to the double area integral (2AI) method.
+
+Usage
+ Minimal example in \c <constant>/viewFactorsDict:
+ \verbatim
+ // Inherited entries
+ ...
+ \endverbatim
+
+ The inherited entries are elaborated in:
+ - \link viewFactorModel.H \endlink
+
+SourceFiles
+ viewFactorModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vf_viewFactor2AI_H
+#define Foam_vf_viewFactor2AI_H
+
+#include "viewFactorModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace VF
+{
+
+/*---------------------------------------------------------------------------*\
+ Class viewFactor2AI Declaration
+\*---------------------------------------------------------------------------*/
+
+class viewFactor2AI
+:
+ public viewFactorModel
+{
+
+protected:
+
+ // Protected Member Functions
+
+ //- Calculate view factor using the double-area integral
+ static scalar calculateFij
+ (
+ const point& xi,
+ const point& xj,
+ const vector& dAi,
+ const vector& dAj
+ );
+
+ //- Calculate
+ virtual scalarListList calculate
+ (
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCf,
+ const vectorField& compactSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("viewFactor2AI");
+
+ //- Constructor
+ viewFactor2AI(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ virtual ~viewFactor2AI() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace VF
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.C
new file mode 100644
index 0000000000000000000000000000000000000000..273b5c2ac5e3f22ccddd99c6de10960989cb8b85
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.C
@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "viewFactor2LI.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace VF
+{
+ defineTypeNameAndDebug(viewFactor2LI, 0);
+ addToRunTimeSelectionTable(viewFactorModel, viewFactor2LI, mesh);
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+Foam::scalar Foam::VF::viewFactor2LI::calculateFij
+(
+ const List<point>& lPoints,
+ const List<point>& rPoints,
+ const scalar alpha
+)
+{
+ scalar Fij = 0;
+
+ forAll(lPoints, i)
+ {
+ // Edge vector and centroid of edge i
+ const vector si(lPoints[i] - lPoints.rcValue(i));
+ const point ci(0.5*(lPoints[i] + lPoints.rcValue(i)));
+
+ forAll(rPoints, j)
+ {
+ // Edge vector and centroid of edge j
+ const vector sj(rPoints[j] - rPoints.rcValue(j));
+ const point cj(0.5*(rPoints[j] + rPoints.rcValue(j)));
+
+ vector r(ci - cj);
+ if (mag(r) < SMALL)
+ {
+ r = alpha*si;
+ }
+
+ Fij += (si & sj)*Foam::log(r & r);
+ }
+ }
+
+ return max(0, 0.25*Fij/mathematical::pi);
+}
+
+
+Foam::scalarListList Foam::VF::viewFactor2LI::calculate
+(
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCf,
+ const vectorField& compactSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+) const
+{
+ // Fill local view factor matrix
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ scalarListList Fij(visibleFaceFaces.size());
+
+ forAll(visibleFaceFaces, facei)
+ {
+ if (debug > 1)
+ {
+ Pout<< "facei:" << facei << "/" << visibleFaceFaces.size()
+ << endl;
+ }
+
+ const labelList& visibleFaces = visibleFaceFaces[facei];
+
+ Fij[facei].resize(visibleFaces.size(), Zero);
+
+ const vector& Ai = compactSf[facei];
+ const scalar magAi = mag(Ai);
+
+ forAll(visibleFaces, visibleFacei)
+ {
+ const label sloti = visibleFaces[visibleFacei];
+ const List<point>& lPoints = compactPoints[facei];
+ const List<point>& rPoints = compactPoints[sloti];
+
+ const scalar Fij2LI = calculateFij(lPoints, rPoints, alpha_);
+
+ Fij[facei][visibleFacei] = Fij2LI/magAi;
+ }
+ }
+
+ return Fij;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::VF::viewFactor2LI::viewFactor2LI
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ viewFactorModel(mesh, dict),
+ alpha_(dict.getOrDefault("alpha", 0.21))
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.H
new file mode 100644
index 0000000000000000000000000000000000000000..0d2c0d273d3403cc4000dc32a54d1d093aecfb12
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.H
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::VF::viewFactor2LI
+
+Description
+ Computes view factors according to the double line integral (2LI) method.
+
+Usage
+ Minimal example in \c <constant>/viewFactorsDict:
+ \verbatim
+ // Optional entries
+ alpha <scalar>;
+
+ // Inherited entries
+ ...
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ alpha | Perturbation for common edges | scalar | no | 0.21
+ \endtable
+
+ The inherited entries are elaborated in:
+ - \link viewFactorModel.H \endlink
+
+SourceFiles
+ viewFactorModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vf_viewFactor2LI_H
+#define Foam_vf_viewFactor2LI_H
+
+#include "viewFactorModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace VF
+{
+
+/*---------------------------------------------------------------------------*\
+ Class viewFactor2LI Declaration
+\*---------------------------------------------------------------------------*/
+
+class viewFactor2LI
+:
+ public viewFactorModel
+{
+protected:
+
+ // Protected Data
+
+ //- Perturbation for common edges; default = 0.21
+ const scalar alpha_;
+
+
+ // Protected Member Functions
+
+ //- Calculate view factor using the double-area integral
+ static scalar calculateFij
+ (
+ const List<point>& lPoints,
+ const List<point>& rPoints,
+ const scalar alpha
+ );
+
+ //- Calculate
+ virtual scalarListList calculate
+ (
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCf,
+ const vectorField& compactSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("viewFactor2LI");
+
+ //- Constructor
+ viewFactor2LI(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ virtual ~viewFactor2LI() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace VF
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.C
new file mode 100644
index 0000000000000000000000000000000000000000..3b352ab69f8c976833f9566fe2977a5c054ebf22
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.C
@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "viewFactorHottel.H"
+#include "mathematicalConstants.H"
+#include "fvMesh.H"
+#include "meshTools.H"
+//#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace VF
+{
+ defineTypeNameAndDebug(viewFactorHottel, 0);
+// addToRunTimeSelectionTable(viewFactorModel, viewFactorHottel, mesh);
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+Foam::scalar Foam::VF::viewFactorHottel::calculateFij
+(
+ const point& p0,
+ const point& p1,
+ const point& p2,
+ const point& p3
+)
+{
+ return 0.5*(mag(p2-p1) + mag(p3-p0) - mag(p2-p0) - mag(p3-p1));
+}
+
+
+Foam::scalarListList Foam::VF::viewFactorHottel::calculate
+(
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCf,
+ const vectorField& compactSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+) const
+{
+ // Fill local view factor matrix
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ scalarListList Fij(visibleFaceFaces.size());
+
+ forAll(visibleFaceFaces, facei)
+ {
+ if (debug > 1)
+ {
+ Pout<< "facei:" << facei << "/" << visibleFaceFaces.size()
+ << endl;
+ }
+
+ const labelList& visibleFaces = visibleFaceFaces[facei];
+
+ Fij[facei].resize_nocopy(visibleFaces.size());
+
+ const point& dCi = compactCf[facei];
+ const vector& Ai = compactSf[facei];
+ const scalar magAi = mag(Ai);
+
+ const vector d1((Ai/magAi) ^ emptyDir_);
+ const vector l1(0.5*magAi/w_*d1);
+ const point p0(dCi + l1);
+ const point p1(dCi - l1);
+
+ forAll(visibleFaces, visibleFacei)
+ {
+ const label sloti = visibleFaces[visibleFacei];
+
+ const point& dCj = compactCf[sloti];
+ const vector& Aj = compactSf[sloti];
+ const scalar magAj = mag(Aj);
+
+ const vector d2((Aj/magAj) ^ emptyDir_);
+ const vector l2(0.5*magAj/w_*d2);
+ const point p2(dCj - l2);
+ const point p3(dCj + l2);
+
+ const scalar FijH = calculateFij(p0, p1, p2, p3);
+
+ Fij[facei][visibleFacei] = FijH/(magAi/w_);
+ }
+ }
+
+
+ return Fij;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::VF::viewFactorHottel::viewFactorHottel
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ viewFactorModel(mesh, dict),
+ emptyDir_(vector::one),
+ w_(0)
+{
+ if (mesh.nSolutionD() != 2)
+ {
+ FatalErrorInFunction
+ << "Hottel crossed strings method only applicable to 2D cases"
+ << exit(FatalError);
+ }
+
+ meshTools::constrainDirection(mesh, mesh.solutionD(), emptyDir_);
+ emptyDir_ = vector::one - emptyDir_;
+ emptyDir_.normalise();
+
+ // 2D width - assume slab
+ // TODO: warn wedge/axisymmetric?
+ w_ = mesh.bounds().span() & emptyDir_;
+
+ Info<< "\nEmpty direction: " << emptyDir_
+ << "\nWidth: " << w_ << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.H
new file mode 100644
index 0000000000000000000000000000000000000000..da7bade587b3356ca40c9e4dfa406c362eee8cc8
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.H
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::VF::viewFactorHottel
+
+Description
+ Computes view factors according to Hottel's crossed strings method.
+
+ Reference:
+ \verbatim
+ Hottel, H. C., & Saforim, A. F. (1967).
+ Radiative transfer.
+ McGraw-Hill Book Company, New York.
+ \endverbatim
+
+Usage
+ Minimal example in \c <constant>/viewFactorsDict:
+ \verbatim
+ // Inherited entries
+ ...
+ \endverbatim
+
+ The inherited entries are elaborated in:
+ - \link viewFactorModel.H \endlink
+
+Note
+ Only applicable to 2D cases
+
+SourceFiles
+ viewFactorModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vf_viewFactorHottel_H
+#define Foam_vf_viewFactorHottel_H
+
+#include "viewFactorModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace VF
+{
+
+/*---------------------------------------------------------------------------*\
+ Class viewFactorHottel Declaration
+\*---------------------------------------------------------------------------*/
+
+class viewFactorHottel
+:
+ public viewFactorModel
+{
+ // Private Data
+
+ //- Mesh empty direction
+ vector emptyDir_;
+
+ //- Mesh 2D width
+ scalar w_;
+
+
+protected:
+
+ // Protected Member Functions
+
+ //- Calculate view factor using the double-area integral
+ static scalar calculateFij
+ (
+ const point& p0,
+ const point& p1,
+ const point& p2,
+ const point& p3
+ );
+
+ //- Calculate
+ virtual scalarListList calculate
+ (
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCf,
+ const vectorField& compactSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("viewFactorHottel");
+
+ //- Constructor
+ viewFactorHottel(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ virtual ~viewFactorHottel() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace VF
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..cc4e1442f6cc596342acc4f6f6011d86b6648a4e
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.C
@@ -0,0 +1,290 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "viewFactorModel.H"
+#include "raySearchEngine.H"
+#include "OBJstream.H"
+#include "volFields.H"
+#include "IOmapDistribute.H"
+#include "scalarListIOList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace VF
+{
+ defineTypeNameAndDebug(viewFactorModel, 0);
+ defineRunTimeSelectionTable(viewFactorModel, mesh);
+}
+}
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+void Foam::VF::viewFactorModel::writeRays
+(
+ const fileName& fName,
+ const pointField& compactCf,
+ const labelListList& visibleFaceFaces
+)
+{
+ OBJstream str(fName);
+
+ Pout<< "Writing rays to " << str.name() << endl;
+
+ forAll(visibleFaceFaces, facei)
+ {
+ const labelList& visibleSlots = visibleFaceFaces[facei];
+
+ for (const label sloti : visibleSlots)
+ {
+ str.write(linePointRef(compactCf[facei], compactCf[sloti]));
+ }
+ }
+
+ str.flush();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::VF::viewFactorModel::viewFactorModel
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+:
+ mesh_(mesh),
+ searchEnginePtr_(raySearchEngine::New(mesh, dict)),
+ writeViewFactors_(dict.get<bool>("writeViewFactors")),
+ writeRays_(dict.getOrDefault<bool>("writeRays", false))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::VF::viewFactorModel::~viewFactorModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::VF::viewFactorModel::calculate()
+{
+ const auto& searchEngine = *searchEnginePtr_;
+
+ labelListList visibleFaceFaces;
+ searchEngine.correct(visibleFaceFaces);
+ const auto& map = searchEngine.map();
+
+ // Construct data in compact addressing
+ // Compact addressing has all local data, followed by remote contributions
+
+ pointField compactCf;
+ vectorField compactSf;
+ List<List<vector>> compactFineSf;
+ List<List<point>> compactFineCf;
+ DynamicList<List<point>> compactPoints;
+ DynamicList<label> compactPatchId;
+
+ searchEngine.compactAddressing
+ (
+ map,
+ compactCf,
+ compactSf,
+ compactFineSf,
+ compactFineCf,
+ compactPoints,
+ compactPatchId
+ );
+
+ if (writeRays_)
+ {
+ // Write all rays between visible faces using .OBJ format
+
+ writeRays
+ (
+ mesh_.time().path()/"allVisibleFaces.obj",
+ compactCf,
+ visibleFaceFaces
+ );
+ }
+
+ (void)mesh_.time().cpuTimeIncrement();
+
+ Info<< "\nCalculating view factors" << endl;
+
+ scalarListIOList Fij
+ (
+ IOobject
+ (
+ "F",
+ mesh_.facesInstance(),
+ mesh_,
+ IOobject::NO_READ
+ ),
+ calculate
+ (
+ visibleFaceFaces,
+ compactCf,
+ compactSf,
+ compactFineSf,
+ compactFineCf,
+ compactPoints,
+ compactPatchId
+ )
+ );
+
+ const label totalPatches = mesh_.boundaryMesh().nNonProcessor();
+
+ // Matrix sum in j(Fij) for each i
+ // Note: if enclosure sum should = 1
+ scalarSquareMatrix viewFactorPatch(totalPatches, Zero);
+
+ forAll(visibleFaceFaces, startFacei)
+ {
+ const scalar magAi = mag(compactSf[startFacei]);
+
+ const labelList& visibleSlots = visibleFaceFaces[startFacei];
+ const label patchi = compactPatchId[startFacei];
+
+ forAll(visibleSlots, visSloti)
+ {
+ const label sloti = visibleSlots[visSloti];
+ const label patchj = compactPatchId[sloti];
+
+ viewFactorPatch[patchi][patchj] += Fij[startFacei][visSloti]*magAi;
+ }
+ }
+
+ reduce(viewFactorPatch, sumOp<scalarSquareMatrix>());
+
+ const scalarList patchArea = searchEngine.patchAreas();
+
+
+ if (UPstream::master())
+ {
+ Info<< "\nPatch view factor contributions:" << nl << endl;
+
+ scalar vfSum = 0;
+
+ const auto& patchIDs = searchEngine.patchIDs();
+ const auto& patches = mesh_.boundaryMesh();
+
+ for (const label patchi : patchIDs)
+ {
+ Info<< " Patch " << patchi << ": " << patches[patchi].name()
+ << endl;
+
+ scalar vfPatch = 0;
+ for (const label patchj: patchIDs)
+ {
+ scalar vf = viewFactorPatch[patchi][patchj]/patchArea[patchi];
+ vfPatch += vf;
+
+ Info<< " F" << patchi << patchj << ": " << vf << endl;
+ }
+
+ Info<< " Sum: " << vfPatch << nl << endl;
+
+ vfSum += vfPatch;
+ }
+
+ Info<< "Sum(all patches) = " << vfSum << endl;
+ }
+
+ // Write view factors matrix in list-list form
+ Info<< "\nWriting view factor matrix" << endl;
+ Fij.write();
+
+ if (writeViewFactors_)
+ {
+ Info<< "\nWriting view factor field" << endl;
+
+ auto tviewFactorField =
+ volScalarField::New
+ (
+ "viewFactorField",
+ mesh_,
+ dimensionedScalar(dimless, Zero)
+ );
+
+ auto& viewFactorField = tviewFactorField.ref();
+
+ searchEngine.interpolate(viewFactorField, Fij);
+
+ viewFactorField.write();
+ }
+
+
+ // Create globalFaceFaces needed to insert view factors
+ // in F to the global matrix Fmatrix
+ {
+ labelListIOList IOglobalFaceFaces
+ (
+ IOobject
+ (
+ "globalFaceFaces",
+ mesh_.facesInstance(),
+ mesh_,
+ IOobject::NO_READ
+ ),
+ visibleFaceFaces.size()
+ );
+
+ forAll(IOglobalFaceFaces, facei)
+ {
+ IOglobalFaceFaces[facei] = renumber
+ (
+ searchEngine.compactToGlobal(),
+ visibleFaceFaces[facei]
+ );
+ }
+
+ IOglobalFaceFaces.write();
+ }
+
+ // Write parallel map
+ {
+ IOmapDistribute IOmapDist
+ (
+ IOobject
+ (
+ "mapDist",
+ mesh_.facesInstance(),
+ mesh_,
+ IOobject::NO_READ
+ ),
+ std::move(map)
+ );
+
+ IOmapDist.write();
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.H b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..549c50978da3f15556a93d72b74903f79024eb84
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.H
@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+ Foam::VF
+
+Description
+ A namespace for various \c viewFactor model implementations.
+
+Class
+ Foam::VF::viewFactorModel
+
+Description
+ A base class for \c viewFactor models.
+
+Usage
+ Minimal example in \c <constant>/viewFactorsDict:
+ \verbatim
+ // Mandatory entries
+ writeViewFactors <bool>;
+
+ // Optional entries
+ writeRays <bool>;
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ writeViewFactors | Flag to write the view factor field | bool | yes | -
+ writeRays | Flag to write the ray geometry | bool | no | false
+ \endtable
+
+SourceFiles
+ viewFactorModel.C
+ viewFactorModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_vf_viewFactorModel_H
+#define Foam_vf_viewFactorModel_H
+
+#include "autoPtr.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class fvMesh;
+
+namespace VF
+{
+
+// Forward Declarations
+class raySearchEngine;
+
+/*---------------------------------------------------------------------------*\
+ Class viewFactorModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class viewFactorModel
+{
+protected:
+
+ // Protected Data
+
+ //- Reference to the mesh database
+ const fvMesh& mesh_;
+
+ //- Run-time selectable ray search engine
+ autoPtr<raySearchEngine> searchEnginePtr_;
+
+ //- Flag to write the view factor field
+ bool writeViewFactors_;
+
+ //- Flag to write the ray geometry
+ bool writeRays_;
+
+
+ // Protected Member Functions
+
+ //- Write ray geometry to file
+ static void writeRays
+ (
+ const fileName& fName,
+ const pointField& compactCf,
+ const labelListList& visibleFaceFaces
+ );
+
+ //- Calculate the view factors using run-time selectable model
+ virtual scalarListList calculate
+ (
+ const labelListList& visibleFaceFaces,
+ const pointField& compactCoarseCf,
+ const vectorField& compactCoarseSf,
+ const List<List<vector>>& compactFineSf,
+ const List<List<point>>& compactFineCf,
+ const DynamicList<List<point>>& compactPoints,
+ const DynamicList<label>& compactPatchId
+ ) const = 0;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("viewFactorModel");
+
+ //- Selection table
+ declareRunTimeSelectionTable
+ (
+ autoPtr,
+ viewFactorModel,
+ mesh,
+ (
+ const fvMesh& mesh,
+ const dictionary& dict
+ ),
+ (mesh, dict)
+ );
+
+ //- Selector
+ static autoPtr<viewFactorModel> New
+ (
+ const fvMesh& mesh,
+ const dictionary& dict
+ );
+
+
+ // Generated Methods
+
+ //- No copy construct
+ viewFactorModel(const viewFactorModel&) = delete;
+
+ //- No copy assignment
+ void operator=(const viewFactorModel&) = delete;
+
+
+ //- Constructor
+ viewFactorModel(const fvMesh& mesh, const dictionary& dict);
+
+ //- Destructor
+ virtual ~viewFactorModel();
+
+
+ // Public Member Functions
+
+ //- Calculate the view factors
+ virtual void calculate();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace VF
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModelNew.C b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..e87f7d459fd696e0a9a2371853e6ad05c902c4f2
--- /dev/null
+++ b/applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModelNew.C
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023-2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "viewFactorModel.H"
+#include "viewFactorHottel.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::VF::viewFactorModel> Foam::VF::viewFactorModel::New
+(
+ const fvMesh& mesh,
+ const dictionary& dict
+)
+{
+ // Intercept 2D cases
+ if (mesh.nSolutionD() == 2)
+ {
+ Info<< "Selecting " << typeName << ": " << viewFactorHottel::typeName
+ << " for 2D cases" << nl << endl;
+ return autoPtr<viewFactorModel>(new viewFactorHottel(mesh, dict));
+ }
+
+
+ // 3D cases - use run-time selection
+
+ const word modelType(dict.get<word>("viewFactorModel"));
+
+ Info<< "Selecting " << typeName << ": " << modelType << endl;
+
+ auto* ctorPtr = meshConstructorTable(modelType);
+
+ if (!ctorPtr)
+ {
+ FatalIOErrorInLookup
+ (
+ dict,
+ typeName,
+ modelType,
+ *meshConstructorTablePtr_
+ ) << exit(FatalIOError);
+ }
+
+ return autoPtr<viewFactorModel>(ctorPtr(mesh, dict));
+}
+
+
+// ************************************************************************* //