From ab5f6dbf4145a0fc955e7533938af9e007f05595 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Fri, 14 Jun 2024 16:53:28 +0100
Subject: [PATCH] ENH: Added new createViewFactors utility

Creates view factors for the view factor radiation model.

User-selectable models:

- raySearchEngine: model to generate rays, i.e. face-to-face connections
- viewFactorModel: model to compute the view factors

For visualisation, use:

- Write the view factors as a volume field

    writeViewFactors    yes;

- Write the rays using OBJ format:

    writeRays       yes; // default = no

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
---
 .../preProcessing/createViewFactors/Allwmake  |  11 +
 .../createViewFactors/Make/files              |   3 +
 .../createViewFactors/Make/options            |   7 +
 .../createViewFactors/createViewFactors.C     | 165 +++
 .../viewFactorModels/Make/files               |  12 +
 .../viewFactorModels/Make/options             |  14 +
 .../raySearchEngine/raySearchEngine.C         | 620 ++++++++++++
 .../raySearchEngine/raySearchEngine.H         | 289 ++++++
 .../raySearchEngine/raySearchEngineI.H        | 101 ++
 .../raySearchEngine/raySearchEngineNew.C      |  59 ++
 .../raySearchEngineTemplates.C                |  86 ++
 .../voxel/voxelRaySearchEngine.C              | 948 ++++++++++++++++++
 .../voxel/voxelRaySearchEngine.H              | 284 ++++++
 .../voxel/voxelRaySearchEngineI.H             | 155 +++
 .../viewFactor2AI/viewFactor2AI.C             | 138 +++
 .../viewFactor2AI/viewFactor2AI.H             | 117 +++
 .../viewFactor2LI/viewFactor2LI.C             | 141 +++
 .../viewFactor2LI/viewFactor2LI.H             | 130 +++
 .../viewFactorHottel/viewFactorHottel.C       | 153 +++
 .../viewFactorHottel/viewFactorHottel.H       | 135 +++
 .../viewFactorModel/viewFactorModel.C         | 290 ++++++
 .../viewFactorModel/viewFactorModel.H         | 185 ++++
 .../viewFactorModel/viewFactorModelNew.C      |  72 ++
 23 files changed, 4115 insertions(+)
 create mode 100755 applications/utilities/preProcessing/createViewFactors/Allwmake
 create mode 100644 applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/files
 create mode 100644 applications/utilities/preProcessing/createViewFactors/createViewFactors/Make/options
 create mode 100644 applications/utilities/preProcessing/createViewFactors/createViewFactors/createViewFactors.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/files
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/Make/options
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngine.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineI.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineNew.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/raySearchEngine/raySearchEngineTemplates.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngine.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/raySearchEngine/voxel/voxelRaySearchEngineI.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2AI/viewFactor2AI.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactor2LI/viewFactor2LI.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorHottel/viewFactorHottel.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.C
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModel.H
 create mode 100644 applications/utilities/preProcessing/createViewFactors/viewFactorModels/viewFactorModel/viewFactorModel/viewFactorModelNew.C

diff --git a/applications/utilities/preProcessing/createViewFactors/Allwmake b/applications/utilities/preProcessing/createViewFactors/Allwmake
new file mode 100755
index 0000000000..9de13a5b03
--- /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 0000000000..3f36ae84ca
--- /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 0000000000..71f2b983ab
--- /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 0000000000..84ad369a99
--- /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 0000000000..eddbe68d67
--- /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 0000000000..ddcc52fdc7
--- /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 0000000000..f1a37bbc5f
--- /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 0000000000..31a6d6f8d9
--- /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 0000000000..bc410e8b21
--- /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 0000000000..78af6adc52
--- /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 0000000000..b345451ada
--- /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 0000000000..757f162044
--- /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 0000000000..1103e1a881
--- /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 0000000000..a46ca7e1da
--- /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 0000000000..2b7afe34bf
--- /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 0000000000..9ba377435b
--- /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 0000000000..273b5c2ac5
--- /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 0000000000..0d2c0d273d
--- /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 0000000000..3b352ab69f
--- /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 0000000000..da7bade587
--- /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 0000000000..cc4e1442f6
--- /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 0000000000..549c50978d
--- /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 0000000000..e87f7d459f
--- /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));
+}
+
+
+// ************************************************************************* //