From b9507c21f98372c90ebad377213759a2f20ba532 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Fri, 8 Apr 2022 10:30:25 -0700
Subject: [PATCH 01/10] ENH: Adding viewFactorsGenExt which uses pbrt for ray
 tracing

---
 .../preProcessing/viewFactorsGen/shootRays.H  |    1 +
 .../viewFactorsGenExt/Make/files              |    3 +
 .../viewFactorsGenExt/Make/options            |   39 +
 .../viewFactorsGenExt/searchingEngineExt.H    |  121 ++
 .../viewFactorsGenExt/shootRaysExt.H          |  196 +++
 .../viewFactorsGenExt/viewFactorsGenExt.C     | 1224 +++++++++++++++++
 .../radiationModels/viewFactor/viewFactor.C   |    2 +
 7 files changed, 1586 insertions(+)
 create mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/Make/files
 create mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/Make/options
 create mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
 create mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
 create mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C

diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
index 1df6e108f0..e35b2d4268 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -49,6 +49,7 @@ for (const int proci : Pstream::allProcs())
 
                     const vector d(remFc - fc);
 
+
                     if (((d & fA) < 0.) && ((d & remA) > 0))
                     {
                         start.append(fc + 0.001*d);
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/Make/files b/applications/utilities/preProcessing/viewFactorsGenExt/Make/files
new file mode 100644
index 0000000000..6d18b01b2c
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/Make/files
@@ -0,0 +1,3 @@
+viewFactorsGenExt.C
+
+EXE = $(FOAM_APPBIN)/viewFactorsGenExt
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/Make/options b/applications/utilities/preProcessing/viewFactorsGenExt/Make/options
new file mode 100644
index 0000000000..2c9476be08
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/Make/options
@@ -0,0 +1,39 @@
+PBRT = /home/sergio/pub/pbrt
+PBRTSRC = $(PBRT)/pbrt-v3/src
+EXTSRC = $(PBRT)/pbrt-v3/src/ext
+GLOGSRC = $(EXTSRC)/glog/src
+
+PBRTBUILD = /home/sergio/pub/pbrt/pbrt-v3-build
+EXTBUILD = $(PBRTBUILD)/src/ext
+GLOGBUILD = $(EXTBUILD)/glog
+
+EXE_INC = \
+    -Wno-old-style-cast \
+    -I$(PBRTSRC)/core \
+    -I$(PBRTSRC)/accelerators \
+    -I$(PBRTSRC)/shapes \
+    -I$(PBRTSRC) \
+    -I$(EXTSRC) \
+    -I$(GLOGSRC) \
+    -I$(GLOGBUILD) \
+    -I./rays \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/distributed/lnInclude \
+
+EXE_LIBS = \
+    -L$(GLOGBUILD) -lglog \
+    -L$(PBRTBUILD) -lpbrt \
+    -L$(EXTBUILD)/openexr/OpenEXR/IlmImf -lIlmImf \
+    -L$(EXTBUILD)/openexr/IlmBase/Imath -lImath \
+    -L$(EXTBUILD)/openexr/IlmBase/Half -lHalf \
+    -L$(EXTBUILD)/openexr/IlmBase/Iex -lIex \
+    -L$(EXTBUILD)/openexr/IlmBase/IexMath -lIexMath \
+    -L$(EXTBUILD)/openexr/IlmBase/IlmThread -lIlmThread \
+    -L$(EXTBUILD)/ptex/src/ptex -lPtex \
+    -lfiniteVolume \
+    -lsurfMesh \
+    -lmeshTools \
+    -ldistributed \
+    -lradiationModels
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
new file mode 100644
index 0000000000..678899fc7d
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
@@ -0,0 +1,121 @@
+Random rndGen(653213);
+
+// Determine mesh bounding boxes:
+List<treeBoundBox> meshBb
+(
+    1,
+    treeBoundBox
+    (
+        boundBox(coarseMesh.points(), false)
+    ).extend(rndGen, 1e-3)
+);
+
+// Dummy bounds dictionary
+dictionary dict;
+dict.add("bounds", meshBb);
+dict.add
+(
+    "distributionType",
+    distributedTriSurfaceMesh::distributionTypeNames_
+    [
+        distributedTriSurfaceMesh::FROZEN
+    ]
+);
+dict.add("mergeDistance", SMALL);
+
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+const triSurface localSurface = triangulate
+(
+    patches,
+    includePatches,
+    finalAgglom,
+    triSurfaceToAgglom,
+    globalNumbering,
+    coarsePatches
+);
+
+
+distributedTriSurfaceMesh surfacesMesh
+(
+    IOobject
+    (
+        "wallSurface.stl",
+        runTime.constant(),     // directory
+        "triSurface",           // instance
+        runTime,                // registry
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    localSurface,
+    dict
+);
+
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
+// pbrt surface
+std::vector<Point3f> vertices;
+
+const pointField& pts = surfacesMesh.localPoints();
+for (const auto& pt : pts)
+{
+    vertices.push_back(Point3f(pt[0], pt[1], pt[2]));
+}
+
+std::vector<int> indices;
+for (const auto& tri : surfacesMesh.localFaces())
+{
+    indices.push_back(tri[0]);
+    indices.push_back(tri[1]);
+    indices.push_back(tri[2]);
+}
+
+std::vector<int> faceIndices;
+faceIndices.reserve(surfacesMesh.localFaces().size());
+
+for (label faceI = 0; faceI < surfacesMesh.localFaces().size(); faceI++)
+{
+    faceIndices.push_back(faceI);
+}
+
+Transform o2w;
+Transform w2o;
+
+std::vector<std::shared_ptr<Shape>> surfacesMesh_pbrt
+(
+    CreateTriangleMesh
+    (
+        &o2w,
+        &w2o,
+        false,                  // bool reverseOrientation
+        surfacesMesh.size(),    // int nTriangles
+        &indices[0],            // int *vertexIndices
+        vertices.size(),
+        &vertices[0],
+        nullptr,
+        nullptr,
+        nullptr,
+        nullptr,
+        nullptr,
+        &faceIndices[0]
+    )
+);
+
+std::vector<std::shared_ptr<Primitive>> prims;
+
+prims.reserve(surfacesMesh_pbrt.size());
+for (auto s : surfacesMesh_pbrt)
+{
+    prims.push_back
+    (
+        std::make_shared<GeometricPrimitive>
+            (s, nullptr, nullptr, nullptr)
+    );
+}
+
+std::shared_ptr<Primitive> accel;
+
+ParamSet paramSet;
+
+accel = CreateBVHAccelerator(std::move(prims), paramSet);
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
new file mode 100644
index 0000000000..29ace62fe6
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
@@ -0,0 +1,196 @@
+// All rays expressed as start face (local) index end end face (global)
+// Pre-size by assuming a certain percentage is visible.
+
+// Maximum length for dynamicList
+const label maxDynListLength
+(
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+);
+
+for (const int proci : Pstream::allProcs())
+{
+
+    std::vector<pbrt::Point3f> start;
+    start.reserve(coarseMesh.nFaces());
+    std::vector<pbrt::Vector3f> dir;
+    dir.reserve(start.size());
+
+    std::vector<pbrt::Point3f> end(start.size());
+    end.reserve(start.size());
+
+    DynamicList<label> startIndex(start.size());
+    DynamicList<label> endIndex(start.size());
+
+    DynamicList<label> startAgg(start.size());
+    DynamicList<label> endAgg(start.size());
+
+    const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
+    const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
+    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
+
+    const pointField& remoteArea = remoteCoarseSf[proci];
+    const pointField& remoteFc = remoteCoarseCf[proci];
+    const labelField& remoteAgg = remoteCoarseAgg[proci];
+
+    pbrt::Vector3f smallV(SMALL, SMALL, SMALL);
+
+    label i = 0;
+    label j = 0;
+    do
+    {
+        for (; i < myFc.size(); i++)
+        {
+            const point& fc = myFc[i];
+            const vector& fA = myArea[i];
+            const label& fAgg = myAgg[i];
+
+            for (; j < remoteFc.size(); j++)
+            {
+                if (proci != Pstream::myProcNo() || i != j)
+                {
+                    const point& remFc = remoteFc[j];
+                    const vector& remA = remoteArea[j];
+                    const label& remAgg = remoteAgg[j];
+
+                    const vector d = remFc - fc;
+
+                    const vector nd = d/mag(d);
+                    const vector nfA = fA/mag(fA);
+                    const vector nremA = remA/mag(remA);
+
+                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
+                    {
+                        pbrt::Vector3f direction(d[0], d[1], d[2]);
+                        direction += smallV;
+                        pbrt::Point3f s(fc[0], fc[1], fc[2]);
+                        end.push_back(s + direction);
+
+                        s += pbrt::Point3f(0.01*d[0], 0.01*d[1], 0.01*d[2]);
+
+                        start.push_back(s);
+                        startIndex.append(i);
+                        startAgg.append(globalNumbering.toGlobal(proci, fAgg));
+
+                        dir.push_back(direction);
+
+                        label globalI = globalNumbering.toGlobal(proci, j);
+                        endIndex.append(globalI);
+                        endAgg.append(globalNumbering.toGlobal(proci, remAgg));
+                        if (startIndex.size() > maxDynListLength)
+                        {
+                            FatalErrorInFunction
+                                << "Dynamic list need from capacity."
+                                << "Actual size maxDynListLength : "
+                                <<  maxDynListLength
+                                << abort(FatalError);
+                        }
+                    }
+                }
+            }
+
+            if (j == remoteFc.size())
+            {
+                j = 0;
+            }
+        }
+
+    } while (i < myFc.size());
+
+    label totalRays(startIndex.size());
+    reduce(totalRays, sumOp<scalar>());
+
+    if (Pstream::master() && debug)
+    {
+        Info<< "Number of rays :" << totalRays << endl;
+    }
+
+    DynamicList<label> dRayIs;
+
+    std::vector<pbrt::Ray> raysInt;
+    raysInt.reserve(start.size());
+
+    std::vector<int> raysTriIndex;
+    raysTriIndex.reserve(start.size());
+
+    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
+    {
+        pbrt::Ray ray(start[rayI], dir[rayI]);
+
+        pbrt::SurfaceInteraction isect;
+        bool intersects = accel->Intersect(ray, &isect);
+        const Triangle* tri = dynamic_cast<const Triangle *>(isect.shape);
+
+        if (intersects)
+        {
+            const int index = tri->faceIndex;
+
+            const Vector3f delta(ray(ray.tMax) - end[rayI]);
+
+            if (delta.Length() < 1e-4)
+            {
+                rayStartFace.append(startIndex[rayI]);
+                rayEndFace.append(endIndex[rayI]);
+            }
+            else if (triSurfaceToAgglom[index] == startAgg[rayI])
+            {
+                dRayIs.append(rayI);
+                raysInt.push_back(ray);
+                raysTriIndex.push_back(index);
+            }
+        }
+    }
+    start.clear();
+
+    labelField rayIs;
+    rayIs.transfer(dRayIs);
+    dRayIs.clear();
+
+    boolList converged(rayIs.size(), false);
+    label nConverged = 0;
+    label iter = 0;
+    do
+    {
+        forAll(rayIs, rayI)
+        {
+            const label rayID = rayIs[rayI];
+            pbrt::Ray& rayHit = raysInt[rayI];
+            const int index =  raysTriIndex[rayI];
+
+            if (!converged[rayI])
+            {
+                if (triSurfaceToAgglom[index] == startAgg[rayID])
+                {
+                    rayHit.tMax *= 1.1;
+                    rayHit.o = rayHit(rayHit.tMax);
+
+                    pbrt::SurfaceInteraction isect;
+                    bool intersects = accel->Intersect(rayHit, &isect);
+                    if (intersects)
+                    {
+                        const Triangle* tri =
+                            dynamic_cast<const Triangle *>(isect.shape);
+                        const int index = tri->faceIndex;
+
+                        raysTriIndex[rayI] = index;
+                    }
+                }
+                else if (triSurfaceToAgglom[index] == endAgg[rayID])
+                {
+                    converged[rayI] = true;
+                    nConverged++;
+                    rayStartFace.append(startIndex[rayID]);
+                    rayEndFace.append(endIndex[rayID]);
+                }
+            }
+        }
+
+        iter ++;
+
+    } while (nConverged < converged.size() && iter < 10);
+
+    startIndex.clear();
+    end.clear();
+    endIndex.clear();
+    startAgg.clear();
+    endAgg.clear();
+}
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C b/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
new file mode 100644
index 0000000000..ed63ff5975
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
@@ -0,0 +1,1224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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
+    viewFactorsGenExt
+
+Group
+    grpPreProcessingUtilities
+
+Description
+    This view factors generation application uses a combined approach of
+    double area integral (2AI) and double linear integral (2LI). 2AI is used
+    when the two surfaces are 'far' apart and 2LI whenre they are 'close'.
+    2LI is integrated along edges using Gaussian quadrature.
+    The distance between faces is calculating a ratio between averaged areas
+    and the distance between face centres.
+
+    The input from viewFactorsDict are:
+
+        tolGaussQuad              0.1;      // GaussQuad  error
+        distTol                   8;        // R/Average(rm)
+        alpha                     0.22;     // Use for common edges for 2LI
+
+
+    For debugging purposes, the following entries can be set in viewFactorsDict:
+
+        writeViewFactorMatrix     true;
+        writeFacesAgglomeration   false;
+        dumpRays		          false;
+
+
+        writeViewFactorMatrix   writes the sum of the VF on each face.
+        writeFacesAgglomeration writes the agglomeration
+        dumpRays                dumps rays
+
+
+    In order to specify the participants patches in the VF calculation the
+    keywaord viewFactorWall should be added to the boundary file.
+
+    floor
+    {
+        type            wall;
+        inGroups        2(wall viewFactorWall);
+        nFaces          100;
+        startFace       3100;
+    }
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "distributedTriSurfaceMesh.H"
+#include "meshTools.H"
+#include "constants.H"
+
+#include "uindirectPrimitivePatch.H"
+#include "DynamicField.H"
+#include "unitConversion.H"
+
+#include "scalarMatrices.H"
+#include "labelListIOList.H"
+#include "scalarListIOList.H"
+
+#include "singleCellFvMesh.H"
+
+#include "IOmapDistribute.H"
+
+#define PBRT_CONSTEXPR constexpr
+#define PBRT_THREAD_LOCAL thread_local
+
+#include "pbrt.h"
+#include "shape.h"
+#include "triangle.h"
+#include "geometry.h"
+#include "paramset.h"
+
+#include "accelerators/bvh.h"
+
+using namespace Foam;
+using namespace Foam::constant;
+using namespace Foam::constant::mathematical;
+using namespace pbrt;
+
+
+triSurface triangulate
+(
+    const polyBoundaryMesh& bMesh,
+    const labelHashSet& includePatches,
+    const labelListIOList& finalAgglom,
+    labelList& triSurfaceToAgglom,
+    const globalIndex& globalNumbering,
+    const polyBoundaryMesh& coarsePatches
+)
+{
+    const polyMesh& mesh = bMesh.mesh();
+
+    // Storage for surfaceMesh. Size estimate.
+    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
+
+    label newPatchI = 0;
+    label localTriFaceI = 0;
+
+    for (const label patchI : includePatches)
+    {
+        const polyPatch& patch = bMesh[patchI];
+        const pointField& points = patch.points();
+
+        label nTriTotal = 0;
+
+        forAll(patch, patchFaceI)
+        {
+            const face& f = patch[patchFaceI];
+
+            faceList triFaces(f.nTriangles(points));
+
+            label nTri = 0;
+
+            f.triangles(points, nTri, triFaces);
+
+            forAll(triFaces, triFaceI)
+            {
+                const face& f = triFaces[triFaceI];
+
+                triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
+
+                nTriTotal++;
+
+                triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
+                (
+                    Pstream::myProcNo(),
+                    finalAgglom[patchI][patchFaceI]
+                  + coarsePatches[patchI].start()
+                );
+            }
+        }
+
+        newPatchI++;
+    }
+
+    //triSurfaceToAgglom.resize(localTriFaceI-1);
+
+    triangles.shrink();
+
+    triSurface rawSurface(triangles, mesh.points());
+
+    triSurface surface
+    (
+        rawSurface.localFaces(),
+        rawSurface.localPoints()
+    );
+
+    // Combine the triSurfaces across all processors
+    if (Pstream::parRun())
+    {
+        List<List<labelledTri>> surfaceProcTris(Pstream::nProcs());
+        List<pointField> surfaceProcPoints(Pstream::nProcs());
+
+        surfaceProcTris[Pstream::myProcNo()] = surface;
+        surfaceProcPoints[Pstream::myProcNo()] = surface.points();
+
+        Pstream::gatherList(surfaceProcTris);
+        Pstream::scatterList(surfaceProcTris);
+        Pstream::gatherList(surfaceProcPoints);
+        Pstream::scatterList(surfaceProcPoints);
+
+        label nTris = 0;
+        forAll(surfaceProcTris, i)
+        {
+            nTris += surfaceProcTris[i].size();
+        }
+
+        List<labelledTri> globalSurfaceTris(nTris);
+        label trii = 0;
+        label offset = 0;
+        forAll(surfaceProcTris, i)
+        {
+            forAll(surfaceProcTris[i], j)
+            {
+                globalSurfaceTris[trii] = surfaceProcTris[i][j];
+                globalSurfaceTris[trii][0] += offset;
+                globalSurfaceTris[trii][1] += offset;
+                globalSurfaceTris[trii][2] += offset;
+                trii++;
+            }
+            offset += surfaceProcPoints[i].size();
+        }
+
+        label nPoints = 0;
+        forAll(surfaceProcPoints, i)
+        {
+            nPoints += surfaceProcPoints[i].size();
+        }
+
+        pointField globalSurfacePoints(nPoints);
+
+        label pointi = 0;
+        forAll(surfaceProcPoints, i)
+        {
+            forAll(surfaceProcPoints[i], j)
+            {
+                globalSurfacePoints[pointi++] = surfaceProcPoints[i][j];
+            }
+        }
+
+        surface = triSurface(globalSurfaceTris, globalSurfacePoints);
+    }
+
+    // Add patch names to surface
+    surface.patches().setSize(newPatchI);
+
+    newPatchI = 0;
+
+    for (const label patchI : includePatches)
+    {
+        const polyPatch& patch = bMesh[patchI];
+
+        surface.patches()[newPatchI].index() = patchI;
+        surface.patches()[newPatchI].name() = patch.name();
+        surface.patches()[newPatchI].geometricType() = patch.type();
+
+        newPatchI++;
+    }
+
+    return surface;
+}
+
+
+void writeRays
+(
+    const fileName& fName,
+    const pointField& compactCf,
+    const pointField& myFc,
+    const labelListList& visibleFaceFaces
+)
+{
+    OFstream str(fName);
+    label vertI = 0;
+
+    Pout<< "Dumping rays to " << str.name() << endl;
+
+    forAll(myFc, faceI)
+    {
+        const labelList visFaces = visibleFaceFaces[faceI];
+        forAll(visFaces, faceRemote)
+        {
+            label compactI = visFaces[faceRemote];
+            const point& remoteFc = compactCf[compactI];
+
+            meshTools::writeOBJ(str, myFc[faceI]);
+            vertI++;
+            meshTools::writeOBJ(str, remoteFc);
+            vertI++;
+            str << "l " << vertI-1 << ' ' << vertI << nl;
+        }
+    }
+    str.flush();
+}
+
+
+scalar calculateViewFactorFij2AI
+(
+    const vector& i,
+    const vector& j,
+    const vector& dAi,
+    const vector& dAj
+)
+{
+    vector r = i - j;
+    scalar rMag = mag(r);
+
+    if (rMag > SMALL)
+    {
+        scalar dAiMag = mag(dAi);
+        scalar dAjMag = mag(dAj);
+
+        vector ni = dAi/dAiMag;
+        vector nj = dAj/dAjMag;
+        scalar cosThetaJ = mag(nj & r)/rMag;
+        scalar cosThetaI = mag(ni & r)/rMag;
+
+        return
+        (
+            (cosThetaI*cosThetaJ*dAjMag*dAiMag)
+           /(sqr(rMag)*constant::mathematical::pi)
+        );
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+void insertMatrixElements
+(
+    const globalIndex& globalNumbering,
+    const label fromProcI,
+    const labelListList& globalFaceFaces,
+    const scalarListList& viewFactors,
+    scalarSquareMatrix& matrix
+)
+{
+    forAll(viewFactors, faceI)
+    {
+        const scalarList& vf = viewFactors[faceI];
+        const labelList& globalFaces = globalFaceFaces[faceI];
+
+        label globalI = globalNumbering.toGlobal(fromProcI, faceI);
+        forAll(globalFaces, i)
+        {
+            matrix[globalI][globalFaces[i]] = vf[i];
+        }
+    }
+}
+
+scalar GaussQuad
+(
+
+    const scalarList& w,
+    const scalarList& p,
+    const scalar& magSi,
+    const scalar& magSj,
+    const vector& di,
+    const vector& dj,
+    const vector& ci,
+    const vector& cj,
+    const scalar cosij,
+    const scalar alpha,
+    label gi
+)
+{
+    scalar dIntFij = 0;
+    if (gi == 0)
+    {
+        vector r(ci - cj);
+        if (mag(r) < SMALL)
+        {
+            r = (alpha*magSi)*di;
+        }
+        dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
+    }
+    else
+    {
+        List<vector> pi(w.size());
+        forAll (pi, i)
+        {
+            pi[i] = ci + p[i]*(magSi/2)*di;
+        }
+
+
+        List<vector> pj(w.size());
+        forAll (pj, i)
+        {
+            pj[i] = cj + p[i]*(magSj/2)*dj;
+        }
+
+        forAll (w, i)
+        {
+            forAll (w, j)
+            {
+                vector r(pi[i] - pj[j]);
+
+                if (mag(r) < SMALL)
+                {
+                    r = (alpha*magSi)*di;
+                     dIntFij +=
+                        cosij*w[i]*w[j]*Foam::log(r&r);
+                }
+                else
+                {
+                    dIntFij +=
+                        cosij*w[i]*w[j]*Foam::log(r&r);
+                }
+
+            }
+        }
+
+        dIntFij *= (magSi/2) * (magSj/2);
+
+    }
+    return dIntFij;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Calculate view factors from face agglomeration array."
+        " The finalAgglom generated by faceAgglomerate utility."
+    );
+
+    #include "addRegionOption.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createNamedMesh.H"
+
+    // Read view factor dictionary
+    IOdictionary viewFactorDict
+    (
+       IOobject
+       (
+            "viewFactorsDict",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+       )
+    );
+
+    const word viewFactorWall("viewFactorWall");
+
+    const bool writeViewFactors =
+        viewFactorDict.getOrDefault("writeViewFactorMatrix", false);
+
+    const bool dumpRays =
+        viewFactorDict.getOrDefault("dumpRays", false);
+
+    const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
+
+    const scalar tolGaussQuad =
+        viewFactorDict.getOrDefault<scalar>("tolGaussQuad", 0.01);
+
+    const scalar distTol =
+         viewFactorDict.getOrDefault<scalar>("distTol", 8);
+
+    const scalar alpha =
+         viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
+
+    // Read agglomeration map
+    labelListIOList finalAgglom
+    (
+        IOobject
+        (
+            "finalAgglom",
+            mesh.facesInstance(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    // Create the coarse mesh  using agglomeration
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (debug)
+    {
+        Pout << "\nCreating single cell mesh..." << endl;
+    }
+
+    singleCellFvMesh coarseMesh
+    (
+        IOobject
+        (
+            "coarse:" + mesh.name(),
+            runTime.timeName(),
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        finalAgglom
+    );
+
+    if (debug)
+    {
+        Pout << "\nCreated single cell mesh..." << endl;
+    }
+
+
+    // Calculate total number of fine and coarse faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    label nCoarseFaces = 0;      //total number of coarse faces
+    label nFineFaces = 0;        //total number of fine faces
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
+
+    labelList viewFactorsPatches(patches.indices(viewFactorWall));
+    for (const label patchi : viewFactorsPatches)
+    {
+        nCoarseFaces += coarsePatches[patchi].size();
+        nFineFaces += patches[patchi].size();
+    }
+
+    // total number of coarse faces
+    label totalNCoarseFaces = nCoarseFaces;
+
+    reduce(totalNCoarseFaces, sumOp<label>());
+
+    if (Pstream::master())
+    {
+        Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
+    }
+
+    if (Pstream::master() && debug)
+    {
+        Pout << "\nView factor patches included in the calculation : "
+             << viewFactorsPatches << endl;
+    }
+
+    // Collect local Cf and Sf on coarse mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    DynamicList<point> localCoarseCf(nCoarseFaces);
+    DynamicList<point> localCoarseSf(nCoarseFaces);
+    DynamicList<label> localAgg(nCoarseFaces);
+    labelHashSet includePatches;
+
+    for (const label patchID : viewFactorsPatches)
+    {
+        const polyPatch& pp = patches[patchID];
+        const labelList& agglom = finalAgglom[patchID];
+
+        includePatches.insert(patchID);
+
+        if (agglom.size() > 0)
+        {
+            label nAgglom = max(agglom)+1;
+            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
+            const labelList& coarsePatchFace =
+                coarseMesh.patchFaceMap()[patchID];
+
+            const pointField& coarseCf =
+                coarseMesh.Cf().boundaryField()[patchID];
+            const pointField& coarseSf =
+                coarseMesh.Sf().boundaryField()[patchID];
+
+            forAll(coarseCf, faceI)
+            {
+                point cf = coarseCf[faceI];
+
+                const label coarseFaceI = coarsePatchFace[faceI];
+                const labelList& fineFaces = coarseToFine[coarseFaceI];
+                const label agglomI =
+                    agglom[fineFaces[0]] + coarsePatches[patchID].start();
+
+                // Construct single face
+                uindirectPrimitivePatch upp
+                (
+                    UIndirectList<face>(pp, fineFaces),
+                    pp.points()
+                );
+
+                List<point> availablePoints
+                (
+                    upp.faceCentres().size()
+                  + upp.localPoints().size()
+                );
+
+                SubList<point>
+                (
+                    availablePoints,
+                    upp.faceCentres().size()
+                ) = upp.faceCentres();
+
+                SubList<point>
+                (
+                    availablePoints,
+                    upp.localPoints().size(),
+                    upp.faceCentres().size()
+                ) = upp.localPoints();
+
+                point cfo = cf;
+                scalar dist = GREAT;
+                forAll(availablePoints, iPoint)
+                {
+                    point cfFine = availablePoints[iPoint];
+                    if (mag(cfFine-cfo) < dist)
+                    {
+                        dist = mag(cfFine-cfo);
+                        cf = cfFine;
+                    }
+                }
+
+                point sf = coarseSf[faceI];
+                localCoarseCf.append(cf);
+                localCoarseSf.append(sf);
+                localAgg.append(agglomI);
+
+            }
+        }
+    }
+
+    // Distribute local coarse Cf and Sf for shooting rays
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<pointField> remoteCoarseCf(Pstream::nProcs());
+    List<pointField> remoteCoarseSf(Pstream::nProcs());
+    List<labelField> remoteCoarseAgg(Pstream::nProcs());
+
+    remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
+    remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
+    remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
+
+    Pstream::gatherList(remoteCoarseCf);
+    Pstream::scatterList(remoteCoarseCf);
+    Pstream::gatherList(remoteCoarseSf);
+    Pstream::scatterList(remoteCoarseSf);
+    Pstream::gatherList(remoteCoarseAgg);
+    Pstream::scatterList(remoteCoarseAgg);
+
+
+    globalIndex globalNumbering(nCoarseFaces);
+
+    // Set up searching engine for obstacles
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    #include "searchingEngineExt.H"
+
+    // Determine rays between coarse face centres
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
+
+    DynamicList<label> rayEndFace(rayStartFace.size());
+
+    // Return rayStartFace in local index and rayEndFace in global index
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+     #include "shootRaysExt.H"
+
+    // Calculate number of visible faces from local index
+    labelList nVisibleFaceFaces(nCoarseFaces, Zero);
+
+    forAll(rayStartFace, i)
+    {
+        nVisibleFaceFaces[rayStartFace[i]]++;
+    }
+
+    labelListList visibleFaceFaces(nCoarseFaces);
+
+    label nViewFactors = 0;
+    forAll(nVisibleFaceFaces, faceI)
+    {
+        visibleFaceFaces[faceI].setSize(nVisibleFaceFaces[faceI]);
+        nViewFactors += nVisibleFaceFaces[faceI];
+    }
+
+    // - 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
+
+    List<Map<label>> compactMap(Pstream::nProcs());
+
+    mapDistribute map(globalNumbering, rayEndFace, compactMap);
+
+    // visibleFaceFaces has:
+    //    (local face, local viewed face) = compact viewed face
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    nVisibleFaceFaces = 0;
+    forAll(rayStartFace, i)
+    {
+        label faceI = rayStartFace[i];
+        label compactI = rayEndFace[i];
+        visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
+    }
+
+    // Construct data in compact addressing
+    // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2LI) need edges (li)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    pointField compactCoarseCf(map.constructSize(), Zero);
+    pointField compactCoarseSf(map.constructSize(), Zero);
+    List<List<point>> compactFineSf(map.constructSize());
+    List<List<point>> compactFineCf(map.constructSize());
+
+    DynamicList<List<point>> compactPoints(map.constructSize());
+
+    DynamicList<label> compactPatchId(map.constructSize());
+
+    // Insert my coarse local values
+    SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
+    SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
+
+    const faceList& faces = mesh.faces();
+
+    // Insert my fine local values
+    label compactI = 0;
+    forAll(viewFactorsPatches, i)
+    {
+        label patchID = viewFactorsPatches[i];
+
+        const labelList& agglom = finalAgglom[patchID];
+        if (agglom.size() > 0)
+        {
+            label nAgglom = max(agglom)+1;
+            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
+            const labelList& coarsePatchFace =
+                coarseMesh.patchFaceMap()[patchID];
+
+            const polyPatch& pp = patches[patchID];
+
+            forAll(coarseToFine, coarseI)
+            {
+                compactPatchId.append(patchID);
+                List<point>& fineCf = compactFineCf[compactI];
+                List<point>& fineSf = compactFineSf[compactI];
+
+                label startFace = pp.start();
+
+                const vectorField locPoints
+                (
+                    mesh.points(),
+                    faces[coarseI + startFace]
+                );
+
+                const label coarseFaceI = coarsePatchFace[coarseI];
+                const labelList& fineFaces = coarseToFine[coarseFaceI];
+
+                fineCf.setSize(fineFaces.size());
+                fineSf.setSize(fineFaces.size());
+
+                compactPoints.append(locPoints);
+
+                fineCf = UIndirectList<point>
+                (
+                    mesh.Cf().boundaryField()[patchID],
+                    coarseToFine[coarseFaceI]
+                );
+                fineSf = UIndirectList<point>
+                (
+                    mesh.Sf().boundaryField()[patchID],
+                    coarseToFine[coarseFaceI]
+                );
+
+                compactI++;
+            }
+        }
+    }
+
+    if (Pstream::master())
+    {
+        Info<< "map distribute..."  << endl;
+    }
+
+    // Do all swapping
+    map.distribute(compactCoarseSf);
+    map.distribute(compactCoarseCf);
+    map.distribute(compactFineCf);
+    map.distribute(compactFineSf);
+
+    map.distribute(compactPoints);
+
+    map.distribute(compactPatchId);
+
+    // Plot all rays between visible faces.
+    if (dumpRays)
+    {
+        writeRays
+        (
+            runTime.path()/"allVisibleFaces.obj",
+            compactCoarseCf,
+            remoteCoarseCf[Pstream::myProcNo()],
+            visibleFaceFaces
+        );
+    }
+
+
+    // Fill local view factor matrix
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    scalarListIOList F2LI
+    (
+        IOobject
+        (
+            "F",
+            mesh.facesInstance(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        nCoarseFaces
+    );
+
+    label totalPatches = coarsePatches.size();
+    reduce(totalPatches, maxOp<label>());
+
+    // Matrix sum in j(Fij) for each i (if enclosure sum = 1)
+    scalarSquareMatrix sumViewFactorPatch
+    (
+        totalPatches,
+        0.0
+    );
+
+    scalarList patchArea(totalPatches, Zero);
+
+    if (Pstream::master())
+    {
+        Info<< "\nCalculating view factors..." << endl;
+    }
+
+    List<scalarList> GaussPoints(5);
+    GaussPoints[0].setSize(1);
+    GaussPoints[0] = 0;
+
+    GaussPoints[1].setSize(2);
+    GaussPoints[1][0] =  1/std::sqrt(3);
+    GaussPoints[1][1] = -1/std::sqrt(3);
+
+    GaussPoints[2].setSize(3);
+    GaussPoints[2][0] =  0;
+    GaussPoints[2][1] =  0.774597;
+    GaussPoints[2][2] = -0.774597;
+
+    GaussPoints[3].setSize(4);
+    GaussPoints[3][0] =  0.339981;
+    GaussPoints[3][1] = -0.339981;
+    GaussPoints[3][2] = 0.861136;
+    GaussPoints[3][3] = -0.861136;
+
+    GaussPoints[4].setSize(5);
+    GaussPoints[4][0] =  0;
+    GaussPoints[4][1] = 0.538469;
+    GaussPoints[4][2] = -0.538469;
+    GaussPoints[4][3] = 0.90618;
+    GaussPoints[4][4] = -0.90618;
+
+
+    List<scalarList> GaussWeights(5);
+    GaussWeights[0].setSize(1);
+    GaussWeights[0] = 2;
+
+    GaussWeights[1].setSize(2);
+    GaussWeights[1][0] =  1;
+    GaussWeights[1][1] =  1;
+
+    GaussWeights[2].setSize(3);
+    GaussWeights[2][0] =  0.888889;
+    GaussWeights[2][1] =  0.555556;
+    GaussWeights[2][2] =  0.555556;
+
+    GaussWeights[3].setSize(4);
+    GaussWeights[3][0] =  0.652145;
+    GaussWeights[3][1] =  0.652145;
+    GaussWeights[3][2] =  0.347855;
+    GaussWeights[3][3] =  0.347855;
+
+    GaussWeights[4].setSize(5);
+    GaussWeights[4][0] =  0.568889;
+    GaussWeights[4][1] =  0.478629;
+    GaussWeights[4][2] =  0.478629;
+    GaussWeights[4][3] =  0.236927;
+    GaussWeights[4][4] =  0.236927;
+
+    label maxQuadOrder = 5;
+
+    if (mesh.nSolutionD() == 3)
+    {
+        forAll(localCoarseSf, coarseFaceI)
+        {
+            const List<point>& localFineSf = compactFineSf[coarseFaceI];
+            const vector Ai = sum(localFineSf);
+            const List<point>& localFineCf = compactFineCf[coarseFaceI];
+            const label fromPatchId = compactPatchId[coarseFaceI];
+
+            const List<point>& lPoints = compactPoints[coarseFaceI];
+
+            patchArea[fromPatchId] += mag(Ai);
+
+            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
+
+            forAll(visCoarseFaces, visCoarseFaceI)
+            {
+                //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
+                F2LI[coarseFaceI].setSize(visCoarseFaces.size());
+                label compactJ = visCoarseFaces[visCoarseFaceI];
+                const List<point>& remoteFineSj = compactFineSf[compactJ];
+                const List<point>& remoteFineCj = compactFineCf[compactJ];
+
+                const List<point>& rPointsCj = compactPoints[compactJ];
+
+                const label toPatchId = compactPatchId[compactJ];
+
+                bool far(false);
+                // Relative distance
+                forAll(localFineSf, i)
+                {
+                    const scalar dAi =
+                        Foam::sqrt
+                        (
+                            mag(localFineSf[i])/constant::mathematical::pi
+                        );
+                    const vector& dCi = localFineCf[i];
+
+                    forAll(remoteFineSj, j)
+                    {
+                        const scalar dAj =
+                            Foam::sqrt
+                            (
+                                mag(remoteFineSj[j])/constant::mathematical::pi
+                            );
+                        const vector& dCj = remoteFineCj[j];
+
+                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
+
+                        if (dist > distTol)
+                        {
+                            far = true;
+                        }
+                    }
+                }
+
+                if (far)
+                {
+                    // 2AI method
+                    scalar F2AIij = 0;
+
+                    forAll(localFineSf, i)
+                    {
+                        const vector& dAi = localFineSf[i];
+                        const vector& dCi = localFineCf[i];
+
+                        forAll(remoteFineSj, j)
+                        {
+                            const vector& dAj = remoteFineSj[j];
+                            const vector& dCj = remoteFineCj[j];
+
+                            scalar dIntFij = calculateViewFactorFij2AI
+                            (
+                                dCi,
+                                dCj,
+                                dAi,
+                                dAj
+                            );
+
+                            F2AIij += dIntFij;
+                        }
+                    }
+                    F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
+                }
+                else
+                {
+                    // 2LI method
+                    label nLocal = lPoints.size();
+                    label nRemote = rPointsCj.size();
+
+                    // Using sub-divisions (quadrature)
+                    scalar oldEToeInt = 0;
+                    for (label gi=0; gi < maxQuadOrder; gi++)
+                    {
+                        scalar F2LIij = 0;
+                        for(label i=0; i<nLocal; i++)
+                        {
+                            vector si;
+                            vector ci;
+
+                            vector sj;
+                            vector cj;
+
+                            if (i == 0)
+                            {
+                                si = lPoints[i] - lPoints[nLocal-1];
+                                ci = (lPoints[i] + lPoints[nLocal-1])/2;
+                            }
+                            else
+                            {
+                                si = lPoints[i] - lPoints[i-1];
+                                ci = (lPoints[i] + lPoints[i-1])/2;
+                            }
+
+                            for(label j=0; j<nRemote; j++)
+                            {
+                                if (j == 0)
+                                {
+                                    sj = rPointsCj[j]-rPointsCj[nRemote-1];
+                                    cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
+                                }
+                                else
+                                {
+                                    sj = rPointsCj[j] - rPointsCj[j-1];
+                                    cj = (rPointsCj[j] + rPointsCj[j-1])/2;
+                                }
+
+
+                                scalar magSi = mag(si);
+                                scalar magSj = mag(sj);
+                                scalar cosij = (si & sj)/(magSi * magSj);
+
+                                vector di = si/magSi;
+                                vector dj = sj/magSj;
+
+                                label quadOrder = gi;
+                                const vector r(ci - cj);
+                                // Common edges use n = 0
+                                if (mag(r) < SMALL)
+                                {
+                                    quadOrder = 0;
+                                }
+
+                                scalar dIntFij =
+                                    GaussQuad
+                                    (
+                                        GaussWeights[gi],
+                                        GaussPoints[gi],
+                                        magSi,
+                                        magSj,
+                                        di,
+                                        dj,
+                                        ci,
+                                        cj,
+                                        cosij,
+                                        alpha,
+                                        quadOrder
+                                    );
+
+                                F2LIij += dIntFij;
+                            }
+                        }
+
+                        scalar err = (F2LIij-oldEToeInt)/F2LIij;
+
+                        if
+                        (
+                            (mag(err) < tolGaussQuad && gi > 0)
+                         || gi == maxQuadOrder-1
+                        )
+                        {
+                            //DebugVar(gi)
+                            F2LI[coarseFaceI][visCoarseFaceI] =
+                                F2LIij/mag(Ai)/4/constant::mathematical::pi;
+                            break;
+                        }
+                        else
+                        {
+                            oldEToeInt = F2LIij;
+                        }
+                    }
+                }
+
+                sumViewFactorPatch[fromPatchId][toPatchId] +=
+                    F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
+            }
+        }
+    }
+    else
+    {
+         FatalErrorInFunction
+            << " View factors are not available in 2D "
+            << exit(FatalError);
+    }
+
+
+
+    // Write view factors matrix in listlist form
+    F2LI.write();
+
+    reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
+    reduce(patchArea, sumOp<scalarList>());
+
+    if (Pstream::master() && debug)
+    {
+        forAll(viewFactorsPatches, i)
+        {
+            label patchI =  viewFactorsPatches[i];
+            for (label j=i; j<viewFactorsPatches.size(); j++)
+            {
+                label patchJ =  viewFactorsPatches[j];
+
+                Info << "F" << patchI << patchJ << ": "
+                     << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
+                     << endl;
+            }
+        }
+    }
+
+    if (writeViewFactors)
+    {
+        if (Pstream::master())
+        {
+            Info << "Writing view factor matrix..." << endl;
+        }
+
+        volScalarField viewFactorField
+        (
+            IOobject
+            (
+                "viewFactorField",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimless, Zero)
+        );
+
+        label compactI = 0;
+
+        volScalarField::Boundary& vfbf = viewFactorField.boundaryFieldRef();
+        forAll(viewFactorsPatches, i)
+        {
+            label patchID = viewFactorsPatches[i];
+            const labelList& agglom = finalAgglom[patchID];
+            if (agglom.size() > 0)
+            {
+                label nAgglom = max(agglom)+1;
+                labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
+                const labelList& coarsePatchFace =
+                    coarseMesh.patchFaceMap()[patchID];
+
+                forAll(coarseToFine, coarseI)
+                {
+                    scalar FiSum = sum(F2LI[compactI]);
+
+                    const label coarseFaceID = coarsePatchFace[coarseI];
+                    const labelList& fineFaces = coarseToFine[coarseFaceID];
+                    forAll(fineFaces, fineId)
+                    {
+                        const label faceID = fineFaces[fineId];
+                        vfbf[patchID][faceID] = FiSum;
+                    }
+                    compactI++;
+                }
+            }
+        }
+        viewFactorField.write();
+    }
+
+
+    // Invert compactMap (from processor+localface to compact) to go
+    // from compact to processor+localface (expressed as a globalIndex)
+    // globalIndex globalCoarFaceNum(coarseMesh.nFaces());
+    labelList compactToGlobal(map.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()
+            );
+        }
+    }
+
+
+    labelListList globalFaceFaces(visibleFaceFaces.size());
+
+    // Create globalFaceFaces needed to insert view factors
+    // in F to the global matrix Fmatrix
+    forAll(globalFaceFaces, faceI)
+    {
+        globalFaceFaces[faceI] = renumber
+        (
+            compactToGlobal,
+            visibleFaceFaces[faceI]
+        );
+    }
+
+    labelListIOList IOglobalFaceFaces
+    (
+        IOobject
+        (
+            "globalFaceFaces",
+            mesh.facesInstance(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        std::move(globalFaceFaces)
+    );
+    IOglobalFaceFaces.write();
+
+
+    IOmapDistribute IOmapDist
+    (
+        IOobject
+        (
+            "mapDist",
+            mesh.facesInstance(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        std::move(map)
+    );
+
+    IOmapDist.write();
+
+    Info<< "End\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index e1f8e4c8cf..cbfa9005be 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -35,6 +35,8 @@ License
 #include "boundaryRadiationProperties.H"
 #include "lduCalculatedProcessorField.H"
 
+
+
 using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

From 565c68f454f046dd08dc04c2aa1c6810f6822749 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Wed, 20 Apr 2022 07:46:08 -0700
Subject: [PATCH 02/10] ENH: Adding intersection margen

---
 .../utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H   | 2 +-
 .../preProcessing/viewFactorsGenExt/viewFactorsGenExt.C        | 3 +++
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
index 29ace62fe6..0cc8240d93 100644
--- a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
@@ -126,7 +126,7 @@ for (const int proci : Pstream::allProcs())
 
             const Vector3f delta(ray(ray.tMax) - end[rayI]);
 
-            if (delta.Length() < 1e-4)
+            if (delta.Length() < intTol)
             {
                 rayStartFace.append(startIndex[rayI]);
                 rayEndFace.append(endIndex[rayI]);
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C b/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
index ed63ff5975..2f5f9ba496 100644
--- a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
+++ b/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
@@ -451,6 +451,9 @@ int main(int argc, char *argv[])
     const scalar alpha =
          viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
 
+    const scalar intTol =
+         viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
+
     // Read agglomeration map
     labelListIOList finalAgglom
     (

From 457979a7b78ce1e6d84d3c2a1c587ddd42f9fb65 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 28 Jun 2022 08:54:53 -0700
Subject: [PATCH 03/10] ENH: Adding new viewFactor utility using CGAL

---
 .../preProcessing/viewFactorsGen/Make/options |    8 +-
 .../viewFactorsGen/searchingEngine.H          |   33 +-
 .../preProcessing/viewFactorsGen/shootRays.H  |  137 +-
 .../viewFactorsGen/viewFactorsGen.C           |  573 ++++++--
 .../viewFactorsGenExt/Make/files              |    3 -
 .../viewFactorsGenExt/Make/options            |   39 -
 .../viewFactorsGenExt/searchingEngineExt.H    |  121 --
 .../viewFactorsGenExt/shootRaysExt.H          |  196 ---
 .../viewFactorsGenExt/viewFactorsGenExt.C     | 1227 -----------------
 9 files changed, 548 insertions(+), 1789 deletions(-)
 delete mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/Make/files
 delete mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/Make/options
 delete mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
 delete mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
 delete mode 100644 applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C

diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/options b/applications/utilities/preProcessing/viewFactorsGen/Make/options
index 1a9826a229..298501f447 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@@ -1,10 +1,16 @@
+include $(GENERAL_RULES)/cgal
+
 EXE_INC = \
+    -Wno-old-style-cast \
+    ${CGAL_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)/parallel/distributed/lnInclude
+
 
 EXE_LIBS = \
+    ${CGAL_LIBS} \
     -lfiniteVolume \
     -lsurfMesh \
     -lmeshTools \
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index df38326ded..9b89b848c9 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -25,7 +25,7 @@ dict.add("mergeDistance", SMALL);
 
 labelList triSurfaceToAgglom(5*nFineFaces);
 
-const triSurface localSurface = triangulate
+triSurface localSurface = triangulate
 (
     patches,
     includePatches,
@@ -36,6 +36,8 @@ const triSurface localSurface = triangulate
 );
 
 
+// CGAL surface
+
 distributedTriSurfaceMesh surfacesMesh
 (
     IOobject
@@ -45,13 +47,38 @@ distributedTriSurfaceMesh surfacesMesh
         "triSurface",           // instance
         runTime,                // registry
         IOobject::NO_READ,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     localSurface,
     dict
 );
 
-
 triSurfaceToAgglom.resize(surfacesMesh.size());
 
 surfacesMesh.setField(triSurfaceToAgglom);
+
+const pointField& pts = surfacesMesh.localPoints();
+
+std::list<Triangle> triangles;
+
+for (const auto& triLocal : surfacesMesh.localFaces())
+{
+    point p1l = pts[triLocal[0]];
+    point p2l = pts[triLocal[1]];
+    point p3l = pts[triLocal[2]];
+
+    Point p1(p1l[0], p1l[1], p1l[2]);
+    Point p2(p2l[0], p2l[1], p2l[2]);
+    Point p3(p3l[0], p3l[1], p3l[2]);
+
+    Triangle tri(p1, p2, p3);
+
+    if (tri.is_degenerate())
+    {
+        std::cout << tri << std::endl;
+    }
+    triangles.push_back(tri);
+}
+
+// constructs AABB tree
+Tree tree(triangles.begin(), triangles.end());
diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
index e35b2d4268..c738784b79 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -1,26 +1,23 @@
-// All rays expressed as start face (local) index end end face (global)
-// Pre-size by assuming a certain percentage is visible.
-
 // Maximum length for dynamicList
 const label maxDynListLength
 (
-    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000)
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
 );
 
 for (const int proci : Pstream::allProcs())
 {
-    // Shoot rays from me to proci. Note that even if processor has
-    // 0 faces we still need to call findLine to keep calls synced.
+    std::vector<Point> start;
+    start.reserve(coarseMesh.nFaces());
+
+    std::vector<Point> end(start.size());
+    end.reserve(start.size());
 
-    DynamicField<point> start(coarseMesh.nFaces());
-    DynamicField<point> end(start.size());
     DynamicList<label> startIndex(start.size());
     DynamicList<label> endIndex(start.size());
 
     DynamicList<label> startAgg(start.size());
     DynamicList<label> endAgg(start.size());
 
-
     const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
     const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
     const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
@@ -50,15 +47,30 @@ for (const int proci : Pstream::allProcs())
                     const vector d(remFc - fc);
 
 
-                    if (((d & fA) < 0.) && ((d & remA) > 0))
+                    const vector nd = d/mag(d);
+                    const vector nfA = fA/mag(fA);
+                    const vector nremA = remA/mag(remA);
+
+                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
                     {
-                        start.append(fc + 0.001*d);
+                        vector direction(d[0], d[1], d[2]);
+                        vector s(fc[0], fc[1], fc[2]);
+                        vector rayEnd(s + (1-intTol)*direction);
+                        end.push_back(Point(rayEnd[0], rayEnd[1], rayEnd[2]));
+
+                        s += vector(intTol*d[0], intTol*d[1], intTol*d[2]);
+
+                        start.push_back(Point(s[0], s[1], s[2]));
                         startIndex.append(i);
-                        startAgg.append(globalNumbering.toGlobal(proci, fAgg));
-                        end.append(fc + 0.999*d);
+                        if (useAgglomeration)
+                        {
+                            startAgg.append(globalNumbering.toGlobal(proci, fAgg));
+                            endAgg.append(globalNumbering.toGlobal(proci, remAgg));
+                        }
+
                         label globalI = globalNumbering.toGlobal(proci, j);
                         endIndex.append(globalI);
-                        endAgg.append(globalNumbering.toGlobal(proci, remAgg));
+
                         if (startIndex.size() > maxDynListLength)
                         {
                             FatalErrorInFunction
@@ -77,98 +89,33 @@ for (const int proci : Pstream::allProcs())
             }
         }
 
-    } while (returnReduce(i < myFc.size(), orOp<bool>()));
+    } while (i < myFc.size());
 
-    List<pointIndexHit> hitInfo(startIndex.size());
-    surfacesMesh.findLine(start, end, hitInfo);
+    label totalRays(startIndex.size());
+    reduce(totalRays, sumOp<scalar>());
 
-    // Return hit global agglo index
-    labelList aggHitIndex;
-    surfacesMesh.getField(hitInfo, aggHitIndex);
-
-    DynamicList<label> dRayIs;
-
-    // Collect the rays which have no obstacle in between rayStartFace
-    // and rayEndFace. If the ray hit itself, it gets stored in dRayIs
-    forAll(hitInfo, rayI)
+    if (debug)
     {
-        if (!hitInfo[rayI].hit())
+        Pout<< "Number of rays :" << totalRays << endl;
+    }
+
+    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
+    {
+        Segment ray(start[rayI], end[rayI]);
+        Segment_intersection intersects = tree.any_intersection(ray);
+
+        if (!intersects)
         {
             rayStartFace.append(startIndex[rayI]);
             rayEndFace.append(endIndex[rayI]);
         }
-        else if (aggHitIndex[rayI] == startAgg[rayI])
-        {
-            dRayIs.append(rayI);
-        }
     }
-
     start.clear();
 
-
-    // Continue rays which hit themself. If they hit the target
-    // agglomeration are added to rayStartFace and rayEndFace
-
-    bool firstLoop = true;
-    DynamicField<point> startHitItself;
-    DynamicField<point> endHitItself;
-    label iter = 0;
-
-    do
+    if (debug)
     {
-        labelField rayIs;
-        rayIs.transfer(dRayIs);
-        dRayIs.clear();
-        forAll(rayIs, rayI)
-        {
-            const label rayID = rayIs[rayI];
-            label hitIndex = -1;
-
-            if (firstLoop)
-            {
-                hitIndex = rayIs[rayI];
-            }
-            else
-            {
-                hitIndex = rayI;
-            }
-
-            if (hitInfo[hitIndex].hit())
-            {
-                if (aggHitIndex[hitIndex] == startAgg[rayID])
-                {
-                    const vector& endP = end[rayID];
-                    const vector& startP = hitInfo[hitIndex].hitPoint();
-                    const vector& d = endP - startP;
-
-                    startHitItself.append(startP + 0.01*d);
-                    endHitItself.append(startP + 1.01*d);
-
-                    dRayIs.append(rayID);
-                }
-                else if (aggHitIndex[hitIndex] == endAgg[rayID])
-                {
-                    rayStartFace.append(startIndex[rayID]);
-                    rayEndFace.append(endIndex[rayID]);
-                }
-
-            }
-        }
-
-        hitInfo.clear();
-        hitInfo.resize(dRayIs.size());
-
-        surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
-
-        surfacesMesh.getField(hitInfo, aggHitIndex);
-
-
-        endHitItself.clear();
-        startHitItself.clear();
-        firstLoop = false;
-        iter ++;
-
-    } while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
+        Pout << "hits : " <<  rayStartFace.size()<< endl;
+    }
 
     startIndex.clear();
     end.clear();
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index fafdfb811b..368fc8df1c 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -4,7 +4,7 @@
    \\    /   O peration     |
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
--------------------------------------------------------------------------------
+// -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
     Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
@@ -25,20 +25,40 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    viewFactorsGen
+    viewFactorsGenExt
 
 Group
     grpPreProcessingUtilities
 
 Description
-    View factors are calculated based on a face agglomeration array
-    (finalAgglom generated by faceAgglomerate utility).
+    This view factors generation application uses a combined approach of
+    double area integral (2AI) and double linear integral (2LI). 2AI is used
+    when the two surfaces are 'far' apart and 2LI whenre they are 'close'.
+    2LI is integrated along edges using Gaussian quadrature.
+    The distance between faces is calculating a ratio between averaged areas
+    and the distance between face centres.
 
-    Each view factor between the agglomerated faces i and j (Fij) is calculated
-    using a double integral of the sub-areas composing the agglomeration.
+    The input from viewFactorsDict are:
 
-    The patches involved in the view factor calculation are taken from the
-    boundary file and should be part on the group viewFactorWall. ie.:
+        tolGaussQuad              0.1;      // GaussQuad  error
+        distTol                   8;        // R/Average(rm)
+        alpha                     0.22;     // Use for common edges for 2LI
+
+
+    For debugging purposes, the following entries can be set in viewFactorsDict:
+
+        writeViewFactorMatrix     true;
+        writeFacesAgglomeration   false;
+        dumpRays		          false;
+
+
+        writeViewFactorMatrix   writes the sum of the VF on each face.
+        writeFacesAgglomeration writes the agglomeration
+        dumpRays                dumps rays
+
+
+    In order to specify the participants patches in the VF calculation the
+    keywaord viewFactorWall should be added to the boundary file.
 
     floor
     {
@@ -56,21 +76,55 @@ Description
 #include "surfaceFields.H"
 #include "distributedTriSurfaceMesh.H"
 #include "meshTools.H"
+#include "constants.H"
 
 #include "indirectPrimitivePatch.H"
 #include "DynamicField.H"
-#include "unitConversion.H"
+//#include "unitConversion.H"
 
 #include "scalarMatrices.H"
 #include "labelListIOList.H"
 #include "scalarListIOList.H"
 
 #include "singleCellFvMesh.H"
-
 #include "IOmapDistribute.H"
 
-using namespace Foam;
+#ifndef NO_CGAL
 
+// Silence boost bind deprecation warnings (before CGAL-5.2.1)
+#include "CGAL/version.h"
+#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000)
+#define BOOST_BIND_GLOBAL_PLACEHOLDERS
+#endif
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+#include <CGAL/Surface_mesh.h>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point;
+typedef K::Direction_3 Vector3C;
+typedef K::Triangle_3 Triangle;
+typedef K::Segment_3 Segment;
+
+typedef std::list<Triangle>::iterator Iterator;
+typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
+typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+typedef boost::optional
+<
+    Tree::Intersection_and_primitive_id<Segment>::Type
+> Segment_intersection;
+
+#endif // NO_CGAL
+
+using namespace Foam;
+using namespace Foam::constant;
+using namespace Foam::constant::mathematical;
+
+//using namespace pbrt;
 
 triSurface triangulate
 (
@@ -127,20 +181,74 @@ triSurface triangulate
         newPatchI++;
     }
 
-    //striSurfaceToAgglom.resize(localTriFaceI-1);
+    //triSurfaceToAgglom.resize(localTriFaceI-1);
 
     triangles.shrink();
 
-    // Create globally numbered tri surface
     triSurface rawSurface(triangles, mesh.points());
 
-    // Create locally numbered tri surface
     triSurface surface
     (
         rawSurface.localFaces(),
         rawSurface.localPoints()
     );
 
+    // Combine the triSurfaces across all processors
+    if (Pstream::parRun())
+    {
+        List<List<labelledTri>> surfaceProcTris(Pstream::nProcs());
+        List<pointField> surfaceProcPoints(Pstream::nProcs());
+
+        surfaceProcTris[Pstream::myProcNo()] = surface;
+        surfaceProcPoints[Pstream::myProcNo()] = surface.points();
+
+        Pstream::gatherList(surfaceProcTris);
+        Pstream::scatterList(surfaceProcTris);
+        Pstream::gatherList(surfaceProcPoints);
+        Pstream::scatterList(surfaceProcPoints);
+
+        label nTris = 0;
+        forAll(surfaceProcTris, i)
+        {
+            nTris += surfaceProcTris[i].size();
+        }
+
+        List<labelledTri> globalSurfaceTris(nTris);
+        label trii = 0;
+        label offset = 0;
+        forAll(surfaceProcTris, i)
+        {
+            forAll(surfaceProcTris[i], j)
+            {
+                globalSurfaceTris[trii] = surfaceProcTris[i][j];
+                globalSurfaceTris[trii][0] += offset;
+                globalSurfaceTris[trii][1] += offset;
+                globalSurfaceTris[trii][2] += offset;
+                trii++;
+            }
+            offset += surfaceProcPoints[i].size();
+        }
+
+        label nPoints = 0;
+        forAll(surfaceProcPoints, i)
+        {
+            nPoints += surfaceProcPoints[i].size();
+        }
+
+        pointField globalSurfacePoints(nPoints);
+
+        label pointi = 0;
+        forAll(surfaceProcPoints, i)
+        {
+            forAll(surfaceProcPoints[i], j)
+            {
+                globalSurfacePoints[pointi++] = surfaceProcPoints[i][j];
+            }
+        }
+
+        surface = triSurface(globalSurfaceTris, globalSurfacePoints);
+    }
+
     // Add patch names to surface
     surface.patches().setSize(newPatchI);
 
@@ -193,7 +301,7 @@ void writeRays
 }
 
 
-scalar calculateViewFactorFij
+scalar calculateViewFactorFij2AI
 (
     const vector& i,
     const vector& j,
@@ -249,7 +357,72 @@ void insertMatrixElements
     }
 }
 
+scalar GaussQuad
+(
 
+    const scalarList& w,
+    const scalarList& p,
+    const scalar& magSi,
+    const scalar& magSj,
+    const vector& di,
+    const vector& dj,
+    const vector& ci,
+    const vector& cj,
+    const scalar cosij,
+    const scalar alpha,
+    label gi
+)
+{
+    scalar dIntFij = 0;
+    if (gi == 0)
+    {
+        vector r(ci - cj);
+        if (mag(r) < SMALL)
+        {
+            r = (alpha*magSi)*di;
+        }
+        dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
+    }
+    else
+    {
+        List<vector> pi(w.size());
+        forAll (pi, i)
+        {
+            pi[i] = ci + p[i]*(magSi/2)*di;
+        }
+
+        List<vector> pj(w.size());
+        forAll (pj, i)
+        {
+            pj[i] = cj + p[i]*(magSj/2)*dj;
+        }
+
+        forAll (w, i)
+        {
+            forAll (w, j)
+            {
+                vector r(pi[i] - pj[j]);
+
+                if (mag(r) < SMALL)
+                {
+                    r = (alpha*magSi)*di;
+                     dIntFij +=
+                        cosij*w[i]*w[j]*Foam::log(r&r);
+                }
+                else
+                {
+                    dIntFij +=
+                        cosij*w[i]*w[j]*Foam::log(r&r);
+                }
+
+            }
+        }
+
+        dIntFij *= (magSi/2) * (magSj/2);
+
+    }
+    return dIntFij;
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -273,7 +446,7 @@ int main(int argc, char *argv[])
             "viewFactorsDict",
             runTime.constant(),
             mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::MUST_READ,
             IOobject::NO_WRITE
        )
     );
@@ -288,6 +461,23 @@ int main(int argc, char *argv[])
 
     const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
 
+    const scalar tolGaussQuad =
+        viewFactorDict.getOrDefault<scalar>("tolGaussQuad", 0.01);
+
+    const scalar distTol =
+         viewFactorDict.getOrDefault<scalar>("distTol", 8);
+
+    const scalar alpha =
+         viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
+
+    const scalar intTol =
+         viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
+
+    bool useAgglomeration(true);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const labelList viewFactorsPatches(patches.indices(viewFactorWall));
+
     // Read agglomeration map
     labelListIOList finalAgglom
     (
@@ -296,12 +486,22 @@ int main(int argc, char *argv[])
             "finalAgglom",
             mesh.facesInstance(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE,
             false
         )
     );
 
+    if (!finalAgglom.typeHeaderOk<labelListIOList>())
+    {
+        finalAgglom.setSize(patches.size());
+        for (label patchi=0;  patchi < patches.size(); patchi++)
+        {
+            finalAgglom[patchi] = identity(patches[patchi].size());
+        }
+        useAgglomeration = false;
+    }
+
     // Create the coarse mesh  using agglomeration
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -336,10 +536,8 @@ int main(int argc, char *argv[])
     label nCoarseFaces = 0;      //total number of coarse faces
     label nFineFaces = 0;        //total number of fine faces
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
     const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
 
-    labelList viewFactorsPatches(patches.indices(viewFactorWall));
     for (const label patchi : viewFactorsPatches)
     {
         nCoarseFaces += coarsePatches[patchi].size();
@@ -473,10 +671,8 @@ int main(int argc, char *argv[])
 
     DynamicList<label> rayEndFace(rayStartFace.size());
 
-
     // Return rayStartFace in local index and rayEndFace in global index
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
     #include "shootRays.H"
 
     // Calculate number of visible faces from local index
@@ -518,9 +714,9 @@ int main(int argc, char *argv[])
         visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
     }
 
-
     // Construct data in compact addressing
-    // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2LI) need edges (li)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     pointField compactCoarseCf(map.constructSize(), Zero);
@@ -528,12 +724,16 @@ int main(int argc, char *argv[])
     List<List<point>> compactFineSf(map.constructSize());
     List<List<point>> compactFineCf(map.constructSize());
 
+    DynamicList<List<point>> compactPoints(map.constructSize());
+
     DynamicList<label> compactPatchId(map.constructSize());
 
     // Insert my coarse local values
     SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
     SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
 
+    const faceList& faces = mesh.faces();
+
     // Insert my fine local values
     label compactI = 0;
     forAll(viewFactorsPatches, i)
@@ -548,11 +748,21 @@ int main(int argc, char *argv[])
             const labelList& coarsePatchFace =
                 coarseMesh.patchFaceMap()[patchID];
 
+            const polyPatch& pp = patches[patchID];
+
             forAll(coarseToFine, coarseI)
             {
                 compactPatchId.append(patchID);
                 List<point>& fineCf = compactFineCf[compactI];
-                List<point>& fineSf = compactFineSf[compactI++];
+                List<point>& fineSf = compactFineSf[compactI];
+
+                label startFace = pp.start();
+
+                const vectorField locPoints
+                (
+                    mesh.points(),
+                    faces[coarseI + startFace]
+                );
 
                 const label coarseFaceI = coarsePatchFace[coarseI];
                 const labelList& fineFaces = coarseToFine[coarseFaceI];
@@ -560,6 +770,8 @@ int main(int argc, char *argv[])
                 fineCf.setSize(fineFaces.size());
                 fineSf.setSize(fineFaces.size());
 
+                compactPoints.append(locPoints);
+
                 fineCf = UIndirectList<point>
                 (
                     mesh.Cf().boundaryField()[patchID],
@@ -570,19 +782,25 @@ int main(int argc, char *argv[])
                     mesh.Sf().boundaryField()[patchID],
                     coarseToFine[coarseFaceI]
                 );
+
+                compactI++;
             }
         }
     }
 
+    if (Pstream::master() && debug)
+    {
+        Info<< "map distribute..."  << endl;
+    }
+
     // Do all swapping
     map.distribute(compactCoarseSf);
     map.distribute(compactCoarseCf);
     map.distribute(compactFineCf);
     map.distribute(compactFineSf);
-
+    map.distribute(compactPoints);
     map.distribute(compactPatchId);
 
-
     // Plot all rays between visible faces.
     if (dumpRays)
     {
@@ -598,8 +816,7 @@ int main(int argc, char *argv[])
 
     // Fill local view factor matrix
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    scalarListIOList F
+    scalarListIOList F2LI
     (
         IOobject
         (
@@ -630,6 +847,61 @@ int main(int argc, char *argv[])
         Info<< "\nCalculating view factors..." << endl;
     }
 
+    FixedList<scalarList, 5> GaussPoints;
+    GaussPoints[0].setSize(1);
+    GaussPoints[0] = 0;
+
+    GaussPoints[1].setSize(2);
+    GaussPoints[1][0] =  1/std::sqrt(3);
+    GaussPoints[1][1] = -1/std::sqrt(3);
+
+    GaussPoints[2].setSize(3);
+    GaussPoints[2][0] =  0;
+    GaussPoints[2][1] =  0.774597;
+    GaussPoints[2][2] = -0.774597;
+
+    GaussPoints[3].setSize(4);
+    GaussPoints[3][0] =  0.339981;
+    GaussPoints[3][1] = -0.339981;
+    GaussPoints[3][2] = 0.861136;
+    GaussPoints[3][3] = -0.861136;
+
+    GaussPoints[4].setSize(5);
+    GaussPoints[4][0] =  0;
+    GaussPoints[4][1] = 0.538469;
+    GaussPoints[4][2] = -0.538469;
+    GaussPoints[4][3] = 0.90618;
+    GaussPoints[4][4] = -0.90618;
+
+
+    List<scalarList> GaussWeights(5);
+    GaussWeights[0].setSize(1);
+    GaussWeights[0] = 2;
+
+    GaussWeights[1].setSize(2);
+    GaussWeights[1][0] =  1;
+    GaussWeights[1][1] =  1;
+
+    GaussWeights[2].setSize(3);
+    GaussWeights[2][0] =  0.888889;
+    GaussWeights[2][1] =  0.555556;
+    GaussWeights[2][2] =  0.555556;
+
+    GaussWeights[3].setSize(4);
+    GaussWeights[3][0] =  0.652145;
+    GaussWeights[3][1] =  0.652145;
+    GaussWeights[3][2] =  0.347855;
+    GaussWeights[3][3] =  0.347855;
+
+    GaussWeights[4].setSize(5);
+    GaussWeights[4][0] =  0.568889;
+    GaussWeights[4][1] =  0.478629;
+    GaussWeights[4][2] =  0.478629;
+    GaussWeights[4][3] =  0.236927;
+    GaussWeights[4][4] =  0.236927;
+
+    label maxQuadOrder = 5;
+
     if (mesh.nSolutionD() == 3)
     {
         forAll(localCoarseSf, coarseFaceI)
@@ -638,119 +910,207 @@ int main(int argc, char *argv[])
             const vector Ai = sum(localFineSf);
             const List<point>& localFineCf = compactFineCf[coarseFaceI];
             const label fromPatchId = compactPatchId[coarseFaceI];
+
+            const List<point>& lPoints = compactPoints[coarseFaceI];
+
             patchArea[fromPatchId] += mag(Ai);
 
             const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 
             forAll(visCoarseFaces, visCoarseFaceI)
             {
-                F[coarseFaceI].setSize(visCoarseFaces.size());
+                //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
+                F2LI[coarseFaceI].setSize(visCoarseFaces.size());
                 label compactJ = visCoarseFaces[visCoarseFaceI];
                 const List<point>& remoteFineSj = compactFineSf[compactJ];
                 const List<point>& remoteFineCj = compactFineCf[compactJ];
 
+                const List<point>& rPointsCj = compactPoints[compactJ];
+
                 const label toPatchId = compactPatchId[compactJ];
 
-                scalar Fij = 0;
+                bool far(false);
+                // Relative distance
                 forAll(localFineSf, i)
                 {
-                    const vector& dAi = localFineSf[i];
+                    const scalar dAi =
+                        Foam::sqrt
+                        (
+                            mag(localFineSf[i])/constant::mathematical::pi
+                        );
                     const vector& dCi = localFineCf[i];
 
                     forAll(remoteFineSj, j)
                     {
-                        const vector& dAj = remoteFineSj[j];
+                        const scalar dAj =
+                            Foam::sqrt
+                            (
+                                mag(remoteFineSj[j])/constant::mathematical::pi
+                            );
                         const vector& dCj = remoteFineCj[j];
 
-                        scalar dIntFij = calculateViewFactorFij
-                        (
-                            dCi,
-                            dCj,
-                            dAi,
-                            dAj
-                        );
+                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
 
-                        Fij += dIntFij;
+                        if (dist > distTol)
+                        {
+                            far = true;
+                        }
                     }
                 }
-                F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
-                sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
+
+                if (far)
+                {
+                    // 2AI method
+                    scalar F2AIij = 0;
+
+                    forAll(localFineSf, i)
+                    {
+                        const vector& dAi = localFineSf[i];
+                        const vector& dCi = localFineCf[i];
+
+                        forAll(remoteFineSj, j)
+                        {
+                            const vector& dAj = remoteFineSj[j];
+                            const vector& dCj = remoteFineCj[j];
+
+                            scalar dIntFij = calculateViewFactorFij2AI
+                            (
+                                dCi,
+                                dCj,
+                                dAi,
+                                dAj
+                            );
+
+                            F2AIij += dIntFij;
+                        }
+                    }
+                    F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
+                }
+                else
+                {
+                    // 2LI method
+                    label nLocal = lPoints.size();
+                    label nRemote = rPointsCj.size();
+
+                    // Using sub-divisions (quadrature)
+                    scalar oldEToeInt = 0;
+                    for (label gi=0; gi < maxQuadOrder; gi++)
+                    {
+                        scalar F2LIij = 0;
+                        for(label i=0; i<nLocal; i++)
+                        {
+                            vector si;
+                            vector ci;
+
+                            vector sj;
+                            vector cj;
+
+                            if (i == 0)
+                            {
+                                si = lPoints[i] - lPoints[nLocal-1];
+                                ci = (lPoints[i] + lPoints[nLocal-1])/2;
+                            }
+                            else
+                            {
+                                si = lPoints[i] - lPoints[i-1];
+                                ci = (lPoints[i] + lPoints[i-1])/2;
+                            }
+
+                            for(label j=0; j<nRemote; j++)
+                            {
+                                if (j == 0)
+                                {
+                                    sj = rPointsCj[j]-rPointsCj[nRemote-1];
+                                    cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
+                                }
+                                else
+                                {
+                                    sj = rPointsCj[j] - rPointsCj[j-1];
+                                    cj = (rPointsCj[j] + rPointsCj[j-1])/2;
+                                }
+
+
+                                scalar magSi = mag(si);
+                                scalar magSj = mag(sj);
+                                scalar cosij = (si & sj)/(magSi * magSj);
+
+                                vector di = si/magSi;
+                                vector dj = sj/magSj;
+
+                                label quadOrder = gi;
+                                const vector r(ci - cj);
+                                // Common edges use n = 0
+                                if (mag(r) < SMALL)
+                                {
+                                    quadOrder = 0;
+                                }
+
+                                scalar dIntFij =
+                                    GaussQuad
+                                    (
+                                        GaussWeights[gi],
+                                        GaussPoints[gi],
+                                        magSi,
+                                        magSj,
+                                        di,
+                                        dj,
+                                        ci,
+                                        cj,
+                                        cosij,
+                                        alpha,
+                                        quadOrder
+                                    );
+
+                                F2LIij += dIntFij;
+                            }
+                        }
+
+                        scalar err = (F2LIij-oldEToeInt)/F2LIij;
+
+                        if
+                        (
+                            (mag(err) < tolGaussQuad && gi > 0)
+                         || gi == maxQuadOrder-1
+                        )
+                        {
+                            F2LI[coarseFaceI][visCoarseFaceI] =
+                                F2LIij/mag(Ai)/4/constant::mathematical::pi;
+                            break;
+                        }
+                        else
+                        {
+                            oldEToeInt = F2LIij;
+                        }
+                    }
+                }
+
+                sumViewFactorPatch[fromPatchId][toPatchId] +=
+                    F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
             }
         }
     }
-    else if (mesh.nSolutionD() == 2)
+    else
     {
-        const boundBox& box = mesh.bounds();
-        const Vector<label>& dirs = mesh.geometricD();
-        vector emptyDir = Zero;
-        forAll(dirs, i)
-        {
-            if (dirs[i] == -1)
-            {
-                emptyDir[i] = 1.0;
-            }
-        }
-
-        scalar wideBy2 = (box.span() & emptyDir)*2.0;
-
-        forAll(localCoarseSf, coarseFaceI)
-        {
-            const vector& Ai = localCoarseSf[coarseFaceI];
-            const vector& Ci = localCoarseCf[coarseFaceI];
-            vector Ain = Ai/mag(Ai);
-            vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
-            vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
-
-            const label fromPatchId = compactPatchId[coarseFaceI];
-            patchArea[fromPatchId] += mag(Ai);
-
-            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
-            forAll(visCoarseFaces, visCoarseFaceI)
-            {
-                F[coarseFaceI].setSize(visCoarseFaces.size());
-                label compactJ = visCoarseFaces[visCoarseFaceI];
-                const vector& Aj = compactCoarseSf[compactJ];
-                const vector& Cj = compactCoarseCf[compactJ];
-
-                const label toPatchId = compactPatchId[compactJ];
-
-                vector Ajn = Aj/mag(Aj);
-                vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
-                vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
-
-                scalar d1 = mag(R1i - R2j);
-                scalar d2 = mag(R2i - R1j);
-                scalar s1 = mag(R1i - R1j);
-                scalar s2 = mag(R2i - R2j);
-
-                scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
-
-                F[coarseFaceI][visCoarseFaceI] = Fij;
-                sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
-            }
-        }
-    }
-
-    if (Pstream::master())
-    {
-        Info << "Writing view factor matrix..." << endl;
+         FatalErrorInFunction
+            << " View factors are not available in 2D "
+            << exit(FatalError);
     }
 
     // Write view factors matrix in listlist form
-    F.write();
+    F2LI.write();
 
     reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
     reduce(patchArea, sumOp<scalarList>());
 
-
     if (Pstream::master() && debug)
     {
         forAll(viewFactorsPatches, i)
         {
             label patchI =  viewFactorsPatches[i];
-            forAll(viewFactorsPatches, i)
+            for (label j=i; j<viewFactorsPatches.size(); j++)
             {
-                label patchJ =  viewFactorsPatches[i];
+                label patchJ =  viewFactorsPatches[j];
+
                 Info << "F" << patchI << patchJ << ": "
                      << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
                      << endl;
@@ -758,9 +1118,13 @@ int main(int argc, char *argv[])
         }
     }
 
-
     if (writeViewFactors)
     {
+        if (Pstream::master())
+        {
+            Info << "Writing view factor matrix..." << endl;
+        }
+
         volScalarField viewFactorField
         (
             IOobject
@@ -791,13 +1155,14 @@ int main(int argc, char *argv[])
 
                 forAll(coarseToFine, coarseI)
                 {
-                    const scalar Fij = sum(F[compactI]);
+                    scalar FiSum = sum(F2LI[compactI]);
+
                     const label coarseFaceID = coarsePatchFace[coarseI];
                     const labelList& fineFaces = coarseToFine[coarseFaceID];
                     forAll(fineFaces, fineId)
                     {
                         const label faceID = fineFaces[fineId];
-                        vfbf[patchID][faceID] = Fij;
+                        vfbf[patchID][faceID] = FiSum;
                     }
                     compactI++;
                 }
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/Make/files b/applications/utilities/preProcessing/viewFactorsGenExt/Make/files
deleted file mode 100644
index 6d18b01b2c..0000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-viewFactorsGenExt.C
-
-EXE = $(FOAM_APPBIN)/viewFactorsGenExt
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/Make/options b/applications/utilities/preProcessing/viewFactorsGenExt/Make/options
deleted file mode 100644
index 2c9476be08..0000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/Make/options
+++ /dev/null
@@ -1,39 +0,0 @@
-PBRT = /home/sergio/pub/pbrt
-PBRTSRC = $(PBRT)/pbrt-v3/src
-EXTSRC = $(PBRT)/pbrt-v3/src/ext
-GLOGSRC = $(EXTSRC)/glog/src
-
-PBRTBUILD = /home/sergio/pub/pbrt/pbrt-v3-build
-EXTBUILD = $(PBRTBUILD)/src/ext
-GLOGBUILD = $(EXTBUILD)/glog
-
-EXE_INC = \
-    -Wno-old-style-cast \
-    -I$(PBRTSRC)/core \
-    -I$(PBRTSRC)/accelerators \
-    -I$(PBRTSRC)/shapes \
-    -I$(PBRTSRC) \
-    -I$(EXTSRC) \
-    -I$(GLOGSRC) \
-    -I$(GLOGBUILD) \
-    -I./rays \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/parallel/distributed/lnInclude \
-
-EXE_LIBS = \
-    -L$(GLOGBUILD) -lglog \
-    -L$(PBRTBUILD) -lpbrt \
-    -L$(EXTBUILD)/openexr/OpenEXR/IlmImf -lIlmImf \
-    -L$(EXTBUILD)/openexr/IlmBase/Imath -lImath \
-    -L$(EXTBUILD)/openexr/IlmBase/Half -lHalf \
-    -L$(EXTBUILD)/openexr/IlmBase/Iex -lIex \
-    -L$(EXTBUILD)/openexr/IlmBase/IexMath -lIexMath \
-    -L$(EXTBUILD)/openexr/IlmBase/IlmThread -lIlmThread \
-    -L$(EXTBUILD)/ptex/src/ptex -lPtex \
-    -lfiniteVolume \
-    -lsurfMesh \
-    -lmeshTools \
-    -ldistributed \
-    -lradiationModels
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
deleted file mode 100644
index 678899fc7d..0000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H
+++ /dev/null
@@ -1,121 +0,0 @@
-Random rndGen(653213);
-
-// Determine mesh bounding boxes:
-List<treeBoundBox> meshBb
-(
-    1,
-    treeBoundBox
-    (
-        boundBox(coarseMesh.points(), false)
-    ).extend(rndGen, 1e-3)
-);
-
-// Dummy bounds dictionary
-dictionary dict;
-dict.add("bounds", meshBb);
-dict.add
-(
-    "distributionType",
-    distributedTriSurfaceMesh::distributionTypeNames_
-    [
-        distributedTriSurfaceMesh::FROZEN
-    ]
-);
-dict.add("mergeDistance", SMALL);
-
-labelList triSurfaceToAgglom(5*nFineFaces);
-
-const triSurface localSurface = triangulate
-(
-    patches,
-    includePatches,
-    finalAgglom,
-    triSurfaceToAgglom,
-    globalNumbering,
-    coarsePatches
-);
-
-
-distributedTriSurfaceMesh surfacesMesh
-(
-    IOobject
-    (
-        "wallSurface.stl",
-        runTime.constant(),     // directory
-        "triSurface",           // instance
-        runTime,                // registry
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    localSurface,
-    dict
-);
-
-
-triSurfaceToAgglom.resize(surfacesMesh.size());
-
-// pbrt surface
-std::vector<Point3f> vertices;
-
-const pointField& pts = surfacesMesh.localPoints();
-for (const auto& pt : pts)
-{
-    vertices.push_back(Point3f(pt[0], pt[1], pt[2]));
-}
-
-std::vector<int> indices;
-for (const auto& tri : surfacesMesh.localFaces())
-{
-    indices.push_back(tri[0]);
-    indices.push_back(tri[1]);
-    indices.push_back(tri[2]);
-}
-
-std::vector<int> faceIndices;
-faceIndices.reserve(surfacesMesh.localFaces().size());
-
-for (label faceI = 0; faceI < surfacesMesh.localFaces().size(); faceI++)
-{
-    faceIndices.push_back(faceI);
-}
-
-Transform o2w;
-Transform w2o;
-
-std::vector<std::shared_ptr<Shape>> surfacesMesh_pbrt
-(
-    CreateTriangleMesh
-    (
-        &o2w,
-        &w2o,
-        false,                  // bool reverseOrientation
-        surfacesMesh.size(),    // int nTriangles
-        &indices[0],            // int *vertexIndices
-        vertices.size(),
-        &vertices[0],
-        nullptr,
-        nullptr,
-        nullptr,
-        nullptr,
-        nullptr,
-        &faceIndices[0]
-    )
-);
-
-std::vector<std::shared_ptr<Primitive>> prims;
-
-prims.reserve(surfacesMesh_pbrt.size());
-for (auto s : surfacesMesh_pbrt)
-{
-    prims.push_back
-    (
-        std::make_shared<GeometricPrimitive>
-            (s, nullptr, nullptr, nullptr)
-    );
-}
-
-std::shared_ptr<Primitive> accel;
-
-ParamSet paramSet;
-
-accel = CreateBVHAccelerator(std::move(prims), paramSet);
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H b/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
deleted file mode 100644
index 0cc8240d93..0000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H
+++ /dev/null
@@ -1,196 +0,0 @@
-// All rays expressed as start face (local) index end end face (global)
-// Pre-size by assuming a certain percentage is visible.
-
-// Maximum length for dynamicList
-const label maxDynListLength
-(
-    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
-);
-
-for (const int proci : Pstream::allProcs())
-{
-
-    std::vector<pbrt::Point3f> start;
-    start.reserve(coarseMesh.nFaces());
-    std::vector<pbrt::Vector3f> dir;
-    dir.reserve(start.size());
-
-    std::vector<pbrt::Point3f> end(start.size());
-    end.reserve(start.size());
-
-    DynamicList<label> startIndex(start.size());
-    DynamicList<label> endIndex(start.size());
-
-    DynamicList<label> startAgg(start.size());
-    DynamicList<label> endAgg(start.size());
-
-    const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
-    const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
-    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
-
-    const pointField& remoteArea = remoteCoarseSf[proci];
-    const pointField& remoteFc = remoteCoarseCf[proci];
-    const labelField& remoteAgg = remoteCoarseAgg[proci];
-
-    pbrt::Vector3f smallV(SMALL, SMALL, SMALL);
-
-    label i = 0;
-    label j = 0;
-    do
-    {
-        for (; i < myFc.size(); i++)
-        {
-            const point& fc = myFc[i];
-            const vector& fA = myArea[i];
-            const label& fAgg = myAgg[i];
-
-            for (; j < remoteFc.size(); j++)
-            {
-                if (proci != Pstream::myProcNo() || i != j)
-                {
-                    const point& remFc = remoteFc[j];
-                    const vector& remA = remoteArea[j];
-                    const label& remAgg = remoteAgg[j];
-
-                    const vector d = remFc - fc;
-
-                    const vector nd = d/mag(d);
-                    const vector nfA = fA/mag(fA);
-                    const vector nremA = remA/mag(remA);
-
-                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
-                    {
-                        pbrt::Vector3f direction(d[0], d[1], d[2]);
-                        direction += smallV;
-                        pbrt::Point3f s(fc[0], fc[1], fc[2]);
-                        end.push_back(s + direction);
-
-                        s += pbrt::Point3f(0.01*d[0], 0.01*d[1], 0.01*d[2]);
-
-                        start.push_back(s);
-                        startIndex.append(i);
-                        startAgg.append(globalNumbering.toGlobal(proci, fAgg));
-
-                        dir.push_back(direction);
-
-                        label globalI = globalNumbering.toGlobal(proci, j);
-                        endIndex.append(globalI);
-                        endAgg.append(globalNumbering.toGlobal(proci, remAgg));
-                        if (startIndex.size() > maxDynListLength)
-                        {
-                            FatalErrorInFunction
-                                << "Dynamic list need from capacity."
-                                << "Actual size maxDynListLength : "
-                                <<  maxDynListLength
-                                << abort(FatalError);
-                        }
-                    }
-                }
-            }
-
-            if (j == remoteFc.size())
-            {
-                j = 0;
-            }
-        }
-
-    } while (i < myFc.size());
-
-    label totalRays(startIndex.size());
-    reduce(totalRays, sumOp<scalar>());
-
-    if (Pstream::master() && debug)
-    {
-        Info<< "Number of rays :" << totalRays << endl;
-    }
-
-    DynamicList<label> dRayIs;
-
-    std::vector<pbrt::Ray> raysInt;
-    raysInt.reserve(start.size());
-
-    std::vector<int> raysTriIndex;
-    raysTriIndex.reserve(start.size());
-
-    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
-    {
-        pbrt::Ray ray(start[rayI], dir[rayI]);
-
-        pbrt::SurfaceInteraction isect;
-        bool intersects = accel->Intersect(ray, &isect);
-        const Triangle* tri = dynamic_cast<const Triangle *>(isect.shape);
-
-        if (intersects)
-        {
-            const int index = tri->faceIndex;
-
-            const Vector3f delta(ray(ray.tMax) - end[rayI]);
-
-            if (delta.Length() < intTol)
-            {
-                rayStartFace.append(startIndex[rayI]);
-                rayEndFace.append(endIndex[rayI]);
-            }
-            else if (triSurfaceToAgglom[index] == startAgg[rayI])
-            {
-                dRayIs.append(rayI);
-                raysInt.push_back(ray);
-                raysTriIndex.push_back(index);
-            }
-        }
-    }
-    start.clear();
-
-    labelField rayIs;
-    rayIs.transfer(dRayIs);
-    dRayIs.clear();
-
-    boolList converged(rayIs.size(), false);
-    label nConverged = 0;
-    label iter = 0;
-    do
-    {
-        forAll(rayIs, rayI)
-        {
-            const label rayID = rayIs[rayI];
-            pbrt::Ray& rayHit = raysInt[rayI];
-            const int index =  raysTriIndex[rayI];
-
-            if (!converged[rayI])
-            {
-                if (triSurfaceToAgglom[index] == startAgg[rayID])
-                {
-                    rayHit.tMax *= 1.1;
-                    rayHit.o = rayHit(rayHit.tMax);
-
-                    pbrt::SurfaceInteraction isect;
-                    bool intersects = accel->Intersect(rayHit, &isect);
-                    if (intersects)
-                    {
-                        const Triangle* tri =
-                            dynamic_cast<const Triangle *>(isect.shape);
-                        const int index = tri->faceIndex;
-
-                        raysTriIndex[rayI] = index;
-                    }
-                }
-                else if (triSurfaceToAgglom[index] == endAgg[rayID])
-                {
-                    converged[rayI] = true;
-                    nConverged++;
-                    rayStartFace.append(startIndex[rayID]);
-                    rayEndFace.append(endIndex[rayID]);
-                }
-            }
-        }
-
-        iter ++;
-
-    } while (nConverged < converged.size() && iter < 10);
-
-    startIndex.clear();
-    end.clear();
-    endIndex.clear();
-    startAgg.clear();
-    endAgg.clear();
-}
diff --git a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C b/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
deleted file mode 100644
index 2f5f9ba496..0000000000
--- a/applications/utilities/preProcessing/viewFactorsGenExt/viewFactorsGenExt.C
+++ /dev/null
@@ -1,1227 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 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
-    viewFactorsGenExt
-
-Group
-    grpPreProcessingUtilities
-
-Description
-    This view factors generation application uses a combined approach of
-    double area integral (2AI) and double linear integral (2LI). 2AI is used
-    when the two surfaces are 'far' apart and 2LI whenre they are 'close'.
-    2LI is integrated along edges using Gaussian quadrature.
-    The distance between faces is calculating a ratio between averaged areas
-    and the distance between face centres.
-
-    The input from viewFactorsDict are:
-
-        tolGaussQuad              0.1;      // GaussQuad  error
-        distTol                   8;        // R/Average(rm)
-        alpha                     0.22;     // Use for common edges for 2LI
-
-
-    For debugging purposes, the following entries can be set in viewFactorsDict:
-
-        writeViewFactorMatrix     true;
-        writeFacesAgglomeration   false;
-        dumpRays		          false;
-
-
-        writeViewFactorMatrix   writes the sum of the VF on each face.
-        writeFacesAgglomeration writes the agglomeration
-        dumpRays                dumps rays
-
-
-    In order to specify the participants patches in the VF calculation the
-    keywaord viewFactorWall should be added to the boundary file.
-
-    floor
-    {
-        type            wall;
-        inGroups        2(wall viewFactorWall);
-        nFaces          100;
-        startFace       3100;
-    }
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "fvMesh.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "distributedTriSurfaceMesh.H"
-#include "meshTools.H"
-#include "constants.H"
-
-#include "uindirectPrimitivePatch.H"
-#include "DynamicField.H"
-#include "unitConversion.H"
-
-#include "scalarMatrices.H"
-#include "labelListIOList.H"
-#include "scalarListIOList.H"
-
-#include "singleCellFvMesh.H"
-
-#include "IOmapDistribute.H"
-
-#define PBRT_CONSTEXPR constexpr
-#define PBRT_THREAD_LOCAL thread_local
-
-#include "pbrt.h"
-#include "shape.h"
-#include "triangle.h"
-#include "geometry.h"
-#include "paramset.h"
-
-#include "accelerators/bvh.h"
-
-using namespace Foam;
-using namespace Foam::constant;
-using namespace Foam::constant::mathematical;
-using namespace pbrt;
-
-
-triSurface triangulate
-(
-    const polyBoundaryMesh& bMesh,
-    const labelHashSet& includePatches,
-    const labelListIOList& finalAgglom,
-    labelList& triSurfaceToAgglom,
-    const globalIndex& globalNumbering,
-    const polyBoundaryMesh& coarsePatches
-)
-{
-    const polyMesh& mesh = bMesh.mesh();
-
-    // Storage for surfaceMesh. Size estimate.
-    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
-
-    label newPatchI = 0;
-    label localTriFaceI = 0;
-
-    for (const label patchI : includePatches)
-    {
-        const polyPatch& patch = bMesh[patchI];
-        const pointField& points = patch.points();
-
-        label nTriTotal = 0;
-
-        forAll(patch, patchFaceI)
-        {
-            const face& f = patch[patchFaceI];
-
-            faceList triFaces(f.nTriangles(points));
-
-            label nTri = 0;
-
-            f.triangles(points, nTri, triFaces);
-
-            forAll(triFaces, triFaceI)
-            {
-                const face& f = triFaces[triFaceI];
-
-                triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
-
-                nTriTotal++;
-
-                triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
-                (
-                    Pstream::myProcNo(),
-                    finalAgglom[patchI][patchFaceI]
-                  + coarsePatches[patchI].start()
-                );
-            }
-        }
-
-        newPatchI++;
-    }
-
-    //triSurfaceToAgglom.resize(localTriFaceI-1);
-
-    triangles.shrink();
-
-    triSurface rawSurface(triangles, mesh.points());
-
-    triSurface surface
-    (
-        rawSurface.localFaces(),
-        rawSurface.localPoints()
-    );
-
-    // Combine the triSurfaces across all processors
-    if (Pstream::parRun())
-    {
-        List<List<labelledTri>> surfaceProcTris(Pstream::nProcs());
-        List<pointField> surfaceProcPoints(Pstream::nProcs());
-
-        surfaceProcTris[Pstream::myProcNo()] = surface;
-        surfaceProcPoints[Pstream::myProcNo()] = surface.points();
-
-        Pstream::gatherList(surfaceProcTris);
-        Pstream::scatterList(surfaceProcTris);
-        Pstream::gatherList(surfaceProcPoints);
-        Pstream::scatterList(surfaceProcPoints);
-
-        label nTris = 0;
-        forAll(surfaceProcTris, i)
-        {
-            nTris += surfaceProcTris[i].size();
-        }
-
-        List<labelledTri> globalSurfaceTris(nTris);
-        label trii = 0;
-        label offset = 0;
-        forAll(surfaceProcTris, i)
-        {
-            forAll(surfaceProcTris[i], j)
-            {
-                globalSurfaceTris[trii] = surfaceProcTris[i][j];
-                globalSurfaceTris[trii][0] += offset;
-                globalSurfaceTris[trii][1] += offset;
-                globalSurfaceTris[trii][2] += offset;
-                trii++;
-            }
-            offset += surfaceProcPoints[i].size();
-        }
-
-        label nPoints = 0;
-        forAll(surfaceProcPoints, i)
-        {
-            nPoints += surfaceProcPoints[i].size();
-        }
-
-        pointField globalSurfacePoints(nPoints);
-
-        label pointi = 0;
-        forAll(surfaceProcPoints, i)
-        {
-            forAll(surfaceProcPoints[i], j)
-            {
-                globalSurfacePoints[pointi++] = surfaceProcPoints[i][j];
-            }
-        }
-
-        surface = triSurface(globalSurfaceTris, globalSurfacePoints);
-    }
-
-    // Add patch names to surface
-    surface.patches().setSize(newPatchI);
-
-    newPatchI = 0;
-
-    for (const label patchI : includePatches)
-    {
-        const polyPatch& patch = bMesh[patchI];
-
-        surface.patches()[newPatchI].index() = patchI;
-        surface.patches()[newPatchI].name() = patch.name();
-        surface.patches()[newPatchI].geometricType() = patch.type();
-
-        newPatchI++;
-    }
-
-    return surface;
-}
-
-
-void writeRays
-(
-    const fileName& fName,
-    const pointField& compactCf,
-    const pointField& myFc,
-    const labelListList& visibleFaceFaces
-)
-{
-    OFstream str(fName);
-    label vertI = 0;
-
-    Pout<< "Dumping rays to " << str.name() << endl;
-
-    forAll(myFc, faceI)
-    {
-        const labelList visFaces = visibleFaceFaces[faceI];
-        forAll(visFaces, faceRemote)
-        {
-            label compactI = visFaces[faceRemote];
-            const point& remoteFc = compactCf[compactI];
-
-            meshTools::writeOBJ(str, myFc[faceI]);
-            vertI++;
-            meshTools::writeOBJ(str, remoteFc);
-            vertI++;
-            str << "l " << vertI-1 << ' ' << vertI << nl;
-        }
-    }
-    str.flush();
-}
-
-
-scalar calculateViewFactorFij2AI
-(
-    const vector& i,
-    const vector& j,
-    const vector& dAi,
-    const vector& dAj
-)
-{
-    vector r = i - j;
-    scalar rMag = mag(r);
-
-    if (rMag > SMALL)
-    {
-        scalar dAiMag = mag(dAi);
-        scalar dAjMag = mag(dAj);
-
-        vector ni = dAi/dAiMag;
-        vector nj = dAj/dAjMag;
-        scalar cosThetaJ = mag(nj & r)/rMag;
-        scalar cosThetaI = mag(ni & r)/rMag;
-
-        return
-        (
-            (cosThetaI*cosThetaJ*dAjMag*dAiMag)
-           /(sqr(rMag)*constant::mathematical::pi)
-        );
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-void insertMatrixElements
-(
-    const globalIndex& globalNumbering,
-    const label fromProcI,
-    const labelListList& globalFaceFaces,
-    const scalarListList& viewFactors,
-    scalarSquareMatrix& matrix
-)
-{
-    forAll(viewFactors, faceI)
-    {
-        const scalarList& vf = viewFactors[faceI];
-        const labelList& globalFaces = globalFaceFaces[faceI];
-
-        label globalI = globalNumbering.toGlobal(fromProcI, faceI);
-        forAll(globalFaces, i)
-        {
-            matrix[globalI][globalFaces[i]] = vf[i];
-        }
-    }
-}
-
-scalar GaussQuad
-(
-
-    const scalarList& w,
-    const scalarList& p,
-    const scalar& magSi,
-    const scalar& magSj,
-    const vector& di,
-    const vector& dj,
-    const vector& ci,
-    const vector& cj,
-    const scalar cosij,
-    const scalar alpha,
-    label gi
-)
-{
-    scalar dIntFij = 0;
-    if (gi == 0)
-    {
-        vector r(ci - cj);
-        if (mag(r) < SMALL)
-        {
-            r = (alpha*magSi)*di;
-        }
-        dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
-    }
-    else
-    {
-        List<vector> pi(w.size());
-        forAll (pi, i)
-        {
-            pi[i] = ci + p[i]*(magSi/2)*di;
-        }
-
-
-        List<vector> pj(w.size());
-        forAll (pj, i)
-        {
-            pj[i] = cj + p[i]*(magSj/2)*dj;
-        }
-
-        forAll (w, i)
-        {
-            forAll (w, j)
-            {
-                vector r(pi[i] - pj[j]);
-
-                if (mag(r) < SMALL)
-                {
-                    r = (alpha*magSi)*di;
-                     dIntFij +=
-                        cosij*w[i]*w[j]*Foam::log(r&r);
-                }
-                else
-                {
-                    dIntFij +=
-                        cosij*w[i]*w[j]*Foam::log(r&r);
-                }
-
-            }
-        }
-
-        dIntFij *= (magSi/2) * (magSj/2);
-
-    }
-    return dIntFij;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::addNote
-    (
-        "Calculate view factors from face agglomeration array."
-        " The finalAgglom generated by faceAgglomerate utility."
-    );
-
-    #include "addRegionOption.H"
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createNamedMesh.H"
-
-    // Read view factor dictionary
-    IOdictionary viewFactorDict
-    (
-       IOobject
-       (
-            "viewFactorsDict",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-       )
-    );
-
-    const word viewFactorWall("viewFactorWall");
-
-    const bool writeViewFactors =
-        viewFactorDict.getOrDefault("writeViewFactorMatrix", false);
-
-    const bool dumpRays =
-        viewFactorDict.getOrDefault("dumpRays", false);
-
-    const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
-
-    const scalar tolGaussQuad =
-        viewFactorDict.getOrDefault<scalar>("tolGaussQuad", 0.01);
-
-    const scalar distTol =
-         viewFactorDict.getOrDefault<scalar>("distTol", 8);
-
-    const scalar alpha =
-         viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
-
-    const scalar intTol =
-         viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
-
-    // Read agglomeration map
-    labelListIOList finalAgglom
-    (
-        IOobject
-        (
-            "finalAgglom",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
-    // Create the coarse mesh  using agglomeration
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    if (debug)
-    {
-        Pout << "\nCreating single cell mesh..." << endl;
-    }
-
-    singleCellFvMesh coarseMesh
-    (
-        IOobject
-        (
-            "coarse:" + mesh.name(),
-            runTime.timeName(),
-            runTime,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        finalAgglom
-    );
-
-    if (debug)
-    {
-        Pout << "\nCreated single cell mesh..." << endl;
-    }
-
-
-    // Calculate total number of fine and coarse faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    label nCoarseFaces = 0;      //total number of coarse faces
-    label nFineFaces = 0;        //total number of fine faces
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
-
-    labelList viewFactorsPatches(patches.indices(viewFactorWall));
-    for (const label patchi : viewFactorsPatches)
-    {
-        nCoarseFaces += coarsePatches[patchi].size();
-        nFineFaces += patches[patchi].size();
-    }
-
-    // total number of coarse faces
-    label totalNCoarseFaces = nCoarseFaces;
-
-    reduce(totalNCoarseFaces, sumOp<label>());
-
-    if (Pstream::master())
-    {
-        Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
-    }
-
-    if (Pstream::master() && debug)
-    {
-        Pout << "\nView factor patches included in the calculation : "
-             << viewFactorsPatches << endl;
-    }
-
-    // Collect local Cf and Sf on coarse mesh
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    DynamicList<point> localCoarseCf(nCoarseFaces);
-    DynamicList<point> localCoarseSf(nCoarseFaces);
-    DynamicList<label> localAgg(nCoarseFaces);
-    labelHashSet includePatches;
-
-    for (const label patchID : viewFactorsPatches)
-    {
-        const polyPatch& pp = patches[patchID];
-        const labelList& agglom = finalAgglom[patchID];
-
-        includePatches.insert(patchID);
-
-        if (agglom.size() > 0)
-        {
-            label nAgglom = max(agglom)+1;
-            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
-            const labelList& coarsePatchFace =
-                coarseMesh.patchFaceMap()[patchID];
-
-            const pointField& coarseCf =
-                coarseMesh.Cf().boundaryField()[patchID];
-            const pointField& coarseSf =
-                coarseMesh.Sf().boundaryField()[patchID];
-
-            forAll(coarseCf, faceI)
-            {
-                point cf = coarseCf[faceI];
-
-                const label coarseFaceI = coarsePatchFace[faceI];
-                const labelList& fineFaces = coarseToFine[coarseFaceI];
-                const label agglomI =
-                    agglom[fineFaces[0]] + coarsePatches[patchID].start();
-
-                // Construct single face
-                uindirectPrimitivePatch upp
-                (
-                    UIndirectList<face>(pp, fineFaces),
-                    pp.points()
-                );
-
-                List<point> availablePoints
-                (
-                    upp.faceCentres().size()
-                  + upp.localPoints().size()
-                );
-
-                SubList<point>
-                (
-                    availablePoints,
-                    upp.faceCentres().size()
-                ) = upp.faceCentres();
-
-                SubList<point>
-                (
-                    availablePoints,
-                    upp.localPoints().size(),
-                    upp.faceCentres().size()
-                ) = upp.localPoints();
-
-                point cfo = cf;
-                scalar dist = GREAT;
-                forAll(availablePoints, iPoint)
-                {
-                    point cfFine = availablePoints[iPoint];
-                    if (mag(cfFine-cfo) < dist)
-                    {
-                        dist = mag(cfFine-cfo);
-                        cf = cfFine;
-                    }
-                }
-
-                point sf = coarseSf[faceI];
-                localCoarseCf.append(cf);
-                localCoarseSf.append(sf);
-                localAgg.append(agglomI);
-
-            }
-        }
-    }
-
-    // Distribute local coarse Cf and Sf for shooting rays
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    List<pointField> remoteCoarseCf(Pstream::nProcs());
-    List<pointField> remoteCoarseSf(Pstream::nProcs());
-    List<labelField> remoteCoarseAgg(Pstream::nProcs());
-
-    remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
-    remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
-    remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
-
-    Pstream::gatherList(remoteCoarseCf);
-    Pstream::scatterList(remoteCoarseCf);
-    Pstream::gatherList(remoteCoarseSf);
-    Pstream::scatterList(remoteCoarseSf);
-    Pstream::gatherList(remoteCoarseAgg);
-    Pstream::scatterList(remoteCoarseAgg);
-
-
-    globalIndex globalNumbering(nCoarseFaces);
-
-    // Set up searching engine for obstacles
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    #include "searchingEngineExt.H"
-
-    // Determine rays between coarse face centres
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
-
-    DynamicList<label> rayEndFace(rayStartFace.size());
-
-    // Return rayStartFace in local index and rayEndFace in global index
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-     #include "shootRaysExt.H"
-
-    // Calculate number of visible faces from local index
-    labelList nVisibleFaceFaces(nCoarseFaces, Zero);
-
-    forAll(rayStartFace, i)
-    {
-        nVisibleFaceFaces[rayStartFace[i]]++;
-    }
-
-    labelListList visibleFaceFaces(nCoarseFaces);
-
-    label nViewFactors = 0;
-    forAll(nVisibleFaceFaces, faceI)
-    {
-        visibleFaceFaces[faceI].setSize(nVisibleFaceFaces[faceI]);
-        nViewFactors += nVisibleFaceFaces[faceI];
-    }
-
-    // - 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
-
-    List<Map<label>> compactMap(Pstream::nProcs());
-
-    mapDistribute map(globalNumbering, rayEndFace, compactMap);
-
-    // visibleFaceFaces has:
-    //    (local face, local viewed face) = compact viewed face
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    nVisibleFaceFaces = 0;
-    forAll(rayStartFace, i)
-    {
-        label faceI = rayStartFace[i];
-        label compactI = rayEndFace[i];
-        visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
-    }
-
-    // Construct data in compact addressing
-    // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
-    // (2LI) need edges (li)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    pointField compactCoarseCf(map.constructSize(), Zero);
-    pointField compactCoarseSf(map.constructSize(), Zero);
-    List<List<point>> compactFineSf(map.constructSize());
-    List<List<point>> compactFineCf(map.constructSize());
-
-    DynamicList<List<point>> compactPoints(map.constructSize());
-
-    DynamicList<label> compactPatchId(map.constructSize());
-
-    // Insert my coarse local values
-    SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
-    SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
-
-    const faceList& faces = mesh.faces();
-
-    // Insert my fine local values
-    label compactI = 0;
-    forAll(viewFactorsPatches, i)
-    {
-        label patchID = viewFactorsPatches[i];
-
-        const labelList& agglom = finalAgglom[patchID];
-        if (agglom.size() > 0)
-        {
-            label nAgglom = max(agglom)+1;
-            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
-            const labelList& coarsePatchFace =
-                coarseMesh.patchFaceMap()[patchID];
-
-            const polyPatch& pp = patches[patchID];
-
-            forAll(coarseToFine, coarseI)
-            {
-                compactPatchId.append(patchID);
-                List<point>& fineCf = compactFineCf[compactI];
-                List<point>& fineSf = compactFineSf[compactI];
-
-                label startFace = pp.start();
-
-                const vectorField locPoints
-                (
-                    mesh.points(),
-                    faces[coarseI + startFace]
-                );
-
-                const label coarseFaceI = coarsePatchFace[coarseI];
-                const labelList& fineFaces = coarseToFine[coarseFaceI];
-
-                fineCf.setSize(fineFaces.size());
-                fineSf.setSize(fineFaces.size());
-
-                compactPoints.append(locPoints);
-
-                fineCf = UIndirectList<point>
-                (
-                    mesh.Cf().boundaryField()[patchID],
-                    coarseToFine[coarseFaceI]
-                );
-                fineSf = UIndirectList<point>
-                (
-                    mesh.Sf().boundaryField()[patchID],
-                    coarseToFine[coarseFaceI]
-                );
-
-                compactI++;
-            }
-        }
-    }
-
-    if (Pstream::master())
-    {
-        Info<< "map distribute..."  << endl;
-    }
-
-    // Do all swapping
-    map.distribute(compactCoarseSf);
-    map.distribute(compactCoarseCf);
-    map.distribute(compactFineCf);
-    map.distribute(compactFineSf);
-
-    map.distribute(compactPoints);
-
-    map.distribute(compactPatchId);
-
-    // Plot all rays between visible faces.
-    if (dumpRays)
-    {
-        writeRays
-        (
-            runTime.path()/"allVisibleFaces.obj",
-            compactCoarseCf,
-            remoteCoarseCf[Pstream::myProcNo()],
-            visibleFaceFaces
-        );
-    }
-
-
-    // Fill local view factor matrix
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalarListIOList F2LI
-    (
-        IOobject
-        (
-            "F",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        nCoarseFaces
-    );
-
-    label totalPatches = coarsePatches.size();
-    reduce(totalPatches, maxOp<label>());
-
-    // Matrix sum in j(Fij) for each i (if enclosure sum = 1)
-    scalarSquareMatrix sumViewFactorPatch
-    (
-        totalPatches,
-        0.0
-    );
-
-    scalarList patchArea(totalPatches, Zero);
-
-    if (Pstream::master())
-    {
-        Info<< "\nCalculating view factors..." << endl;
-    }
-
-    List<scalarList> GaussPoints(5);
-    GaussPoints[0].setSize(1);
-    GaussPoints[0] = 0;
-
-    GaussPoints[1].setSize(2);
-    GaussPoints[1][0] =  1/std::sqrt(3);
-    GaussPoints[1][1] = -1/std::sqrt(3);
-
-    GaussPoints[2].setSize(3);
-    GaussPoints[2][0] =  0;
-    GaussPoints[2][1] =  0.774597;
-    GaussPoints[2][2] = -0.774597;
-
-    GaussPoints[3].setSize(4);
-    GaussPoints[3][0] =  0.339981;
-    GaussPoints[3][1] = -0.339981;
-    GaussPoints[3][2] = 0.861136;
-    GaussPoints[3][3] = -0.861136;
-
-    GaussPoints[4].setSize(5);
-    GaussPoints[4][0] =  0;
-    GaussPoints[4][1] = 0.538469;
-    GaussPoints[4][2] = -0.538469;
-    GaussPoints[4][3] = 0.90618;
-    GaussPoints[4][4] = -0.90618;
-
-
-    List<scalarList> GaussWeights(5);
-    GaussWeights[0].setSize(1);
-    GaussWeights[0] = 2;
-
-    GaussWeights[1].setSize(2);
-    GaussWeights[1][0] =  1;
-    GaussWeights[1][1] =  1;
-
-    GaussWeights[2].setSize(3);
-    GaussWeights[2][0] =  0.888889;
-    GaussWeights[2][1] =  0.555556;
-    GaussWeights[2][2] =  0.555556;
-
-    GaussWeights[3].setSize(4);
-    GaussWeights[3][0] =  0.652145;
-    GaussWeights[3][1] =  0.652145;
-    GaussWeights[3][2] =  0.347855;
-    GaussWeights[3][3] =  0.347855;
-
-    GaussWeights[4].setSize(5);
-    GaussWeights[4][0] =  0.568889;
-    GaussWeights[4][1] =  0.478629;
-    GaussWeights[4][2] =  0.478629;
-    GaussWeights[4][3] =  0.236927;
-    GaussWeights[4][4] =  0.236927;
-
-    label maxQuadOrder = 5;
-
-    if (mesh.nSolutionD() == 3)
-    {
-        forAll(localCoarseSf, coarseFaceI)
-        {
-            const List<point>& localFineSf = compactFineSf[coarseFaceI];
-            const vector Ai = sum(localFineSf);
-            const List<point>& localFineCf = compactFineCf[coarseFaceI];
-            const label fromPatchId = compactPatchId[coarseFaceI];
-
-            const List<point>& lPoints = compactPoints[coarseFaceI];
-
-            patchArea[fromPatchId] += mag(Ai);
-
-            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
-
-            forAll(visCoarseFaces, visCoarseFaceI)
-            {
-                //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
-                F2LI[coarseFaceI].setSize(visCoarseFaces.size());
-                label compactJ = visCoarseFaces[visCoarseFaceI];
-                const List<point>& remoteFineSj = compactFineSf[compactJ];
-                const List<point>& remoteFineCj = compactFineCf[compactJ];
-
-                const List<point>& rPointsCj = compactPoints[compactJ];
-
-                const label toPatchId = compactPatchId[compactJ];
-
-                bool far(false);
-                // Relative distance
-                forAll(localFineSf, i)
-                {
-                    const scalar dAi =
-                        Foam::sqrt
-                        (
-                            mag(localFineSf[i])/constant::mathematical::pi
-                        );
-                    const vector& dCi = localFineCf[i];
-
-                    forAll(remoteFineSj, j)
-                    {
-                        const scalar dAj =
-                            Foam::sqrt
-                            (
-                                mag(remoteFineSj[j])/constant::mathematical::pi
-                            );
-                        const vector& dCj = remoteFineCj[j];
-
-                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
-
-                        if (dist > distTol)
-                        {
-                            far = true;
-                        }
-                    }
-                }
-
-                if (far)
-                {
-                    // 2AI method
-                    scalar F2AIij = 0;
-
-                    forAll(localFineSf, i)
-                    {
-                        const vector& dAi = localFineSf[i];
-                        const vector& dCi = localFineCf[i];
-
-                        forAll(remoteFineSj, j)
-                        {
-                            const vector& dAj = remoteFineSj[j];
-                            const vector& dCj = remoteFineCj[j];
-
-                            scalar dIntFij = calculateViewFactorFij2AI
-                            (
-                                dCi,
-                                dCj,
-                                dAi,
-                                dAj
-                            );
-
-                            F2AIij += dIntFij;
-                        }
-                    }
-                    F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
-                }
-                else
-                {
-                    // 2LI method
-                    label nLocal = lPoints.size();
-                    label nRemote = rPointsCj.size();
-
-                    // Using sub-divisions (quadrature)
-                    scalar oldEToeInt = 0;
-                    for (label gi=0; gi < maxQuadOrder; gi++)
-                    {
-                        scalar F2LIij = 0;
-                        for(label i=0; i<nLocal; i++)
-                        {
-                            vector si;
-                            vector ci;
-
-                            vector sj;
-                            vector cj;
-
-                            if (i == 0)
-                            {
-                                si = lPoints[i] - lPoints[nLocal-1];
-                                ci = (lPoints[i] + lPoints[nLocal-1])/2;
-                            }
-                            else
-                            {
-                                si = lPoints[i] - lPoints[i-1];
-                                ci = (lPoints[i] + lPoints[i-1])/2;
-                            }
-
-                            for(label j=0; j<nRemote; j++)
-                            {
-                                if (j == 0)
-                                {
-                                    sj = rPointsCj[j]-rPointsCj[nRemote-1];
-                                    cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
-                                }
-                                else
-                                {
-                                    sj = rPointsCj[j] - rPointsCj[j-1];
-                                    cj = (rPointsCj[j] + rPointsCj[j-1])/2;
-                                }
-
-
-                                scalar magSi = mag(si);
-                                scalar magSj = mag(sj);
-                                scalar cosij = (si & sj)/(magSi * magSj);
-
-                                vector di = si/magSi;
-                                vector dj = sj/magSj;
-
-                                label quadOrder = gi;
-                                const vector r(ci - cj);
-                                // Common edges use n = 0
-                                if (mag(r) < SMALL)
-                                {
-                                    quadOrder = 0;
-                                }
-
-                                scalar dIntFij =
-                                    GaussQuad
-                                    (
-                                        GaussWeights[gi],
-                                        GaussPoints[gi],
-                                        magSi,
-                                        magSj,
-                                        di,
-                                        dj,
-                                        ci,
-                                        cj,
-                                        cosij,
-                                        alpha,
-                                        quadOrder
-                                    );
-
-                                F2LIij += dIntFij;
-                            }
-                        }
-
-                        scalar err = (F2LIij-oldEToeInt)/F2LIij;
-
-                        if
-                        (
-                            (mag(err) < tolGaussQuad && gi > 0)
-                         || gi == maxQuadOrder-1
-                        )
-                        {
-                            //DebugVar(gi)
-                            F2LI[coarseFaceI][visCoarseFaceI] =
-                                F2LIij/mag(Ai)/4/constant::mathematical::pi;
-                            break;
-                        }
-                        else
-                        {
-                            oldEToeInt = F2LIij;
-                        }
-                    }
-                }
-
-                sumViewFactorPatch[fromPatchId][toPatchId] +=
-                    F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
-            }
-        }
-    }
-    else
-    {
-         FatalErrorInFunction
-            << " View factors are not available in 2D "
-            << exit(FatalError);
-    }
-
-
-
-    // Write view factors matrix in listlist form
-    F2LI.write();
-
-    reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
-    reduce(patchArea, sumOp<scalarList>());
-
-    if (Pstream::master() && debug)
-    {
-        forAll(viewFactorsPatches, i)
-        {
-            label patchI =  viewFactorsPatches[i];
-            for (label j=i; j<viewFactorsPatches.size(); j++)
-            {
-                label patchJ =  viewFactorsPatches[j];
-
-                Info << "F" << patchI << patchJ << ": "
-                     << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
-                     << endl;
-            }
-        }
-    }
-
-    if (writeViewFactors)
-    {
-        if (Pstream::master())
-        {
-            Info << "Writing view factor matrix..." << endl;
-        }
-
-        volScalarField viewFactorField
-        (
-            IOobject
-            (
-                "viewFactorField",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar(dimless, Zero)
-        );
-
-        label compactI = 0;
-
-        volScalarField::Boundary& vfbf = viewFactorField.boundaryFieldRef();
-        forAll(viewFactorsPatches, i)
-        {
-            label patchID = viewFactorsPatches[i];
-            const labelList& agglom = finalAgglom[patchID];
-            if (agglom.size() > 0)
-            {
-                label nAgglom = max(agglom)+1;
-                labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
-                const labelList& coarsePatchFace =
-                    coarseMesh.patchFaceMap()[patchID];
-
-                forAll(coarseToFine, coarseI)
-                {
-                    scalar FiSum = sum(F2LI[compactI]);
-
-                    const label coarseFaceID = coarsePatchFace[coarseI];
-                    const labelList& fineFaces = coarseToFine[coarseFaceID];
-                    forAll(fineFaces, fineId)
-                    {
-                        const label faceID = fineFaces[fineId];
-                        vfbf[patchID][faceID] = FiSum;
-                    }
-                    compactI++;
-                }
-            }
-        }
-        viewFactorField.write();
-    }
-
-
-    // Invert compactMap (from processor+localface to compact) to go
-    // from compact to processor+localface (expressed as a globalIndex)
-    // globalIndex globalCoarFaceNum(coarseMesh.nFaces());
-    labelList compactToGlobal(map.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()
-            );
-        }
-    }
-
-
-    labelListList globalFaceFaces(visibleFaceFaces.size());
-
-    // Create globalFaceFaces needed to insert view factors
-    // in F to the global matrix Fmatrix
-    forAll(globalFaceFaces, faceI)
-    {
-        globalFaceFaces[faceI] = renumber
-        (
-            compactToGlobal,
-            visibleFaceFaces[faceI]
-        );
-    }
-
-    labelListIOList IOglobalFaceFaces
-    (
-        IOobject
-        (
-            "globalFaceFaces",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        std::move(globalFaceFaces)
-    );
-    IOglobalFaceFaces.write();
-
-
-    IOmapDistribute IOmapDist
-    (
-        IOobject
-        (
-            "mapDist",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        std::move(map)
-    );
-
-    IOmapDist.write();
-
-    Info<< "End\n" << endl;
-    return 0;
-}
-
-
-// ************************************************************************* //

From f5598706c5e28eb6275800016ebc0d60e678801a Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 28 Jun 2022 16:04:14 -0700
Subject: [PATCH 04/10] ENH: Making agglomeration optional for viewFactor model

---
 .../radiationModels/viewFactor/viewFactor.C   | 101 ++++++++++++------
 .../radiationModels/viewFactor/viewFactor.H   |   4 +-
 2 files changed, 68 insertions(+), 37 deletions(-)

diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index cbfa9005be..4d8dfb3c87 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -57,7 +57,36 @@ const Foam::word Foam::radiation::viewFactor::viewFactorWalls
 
 void Foam::radiation::viewFactor::initialise()
 {
-    const polyBoundaryMesh& coarsePatches = coarseMesh_.boundaryMesh();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    if (!finalAgglom_.typeHeaderOk<labelListIOList>())
+    {
+        finalAgglom_.setSize(patches.size());
+        for (label patchi=0;  patchi < patches.size(); patchi++)
+        {
+            finalAgglom_[patchi] = identity(patches[patchi].size());
+        }
+    }
+
+    coarseMesh_.reset
+    (
+        new singleCellFvMesh
+        (
+            IOobject
+            (
+                "coarse:" + mesh_.name(),
+                mesh_.polyMesh::instance(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            finalAgglom_
+        )
+    );
+
+    const polyBoundaryMesh& coarsePatches = coarseMesh_->boundaryMesh();
 
     selectedPatches_ = mesh_.boundaryMesh().indices(viewFactorWalls);
 
@@ -82,6 +111,8 @@ void Foam::radiation::viewFactor::initialise()
 
     useDirect_ = coeffs_.getOrDefault<bool>("useDirectSolver", true);
 
+
+
     map_.reset
     (
         new IOmapDistribute
@@ -512,26 +543,26 @@ Foam::radiation::viewFactor::viewFactor(const volScalarField& T)
             "finalAgglom",
             mesh_.facesInstance(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE,
             false
         )
     ),
     map_(),
-    coarseMesh_
-    (
-        IOobject
-        (
-            "coarse:" + mesh_.name(),
-            mesh_.polyMesh::instance(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh_,
-        finalAgglom_
-    ),
+    coarseMesh_(),
+//     (
+//         IOobject
+//         (
+//             "coarse:" + mesh_.name(),
+//             mesh_.polyMesh::instance(),
+//             mesh_.time(),
+//             IOobject::NO_READ,
+//             IOobject::NO_WRITE,
+//             false
+//         ),
+//         mesh_,
+//         finalAgglom_
+//     ),
     qr_
     (
         IOobject
@@ -575,26 +606,26 @@ Foam::radiation::viewFactor::viewFactor
             "finalAgglom",
             mesh_.facesInstance(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE,
             false
         )
     ),
     map_(),
-    coarseMesh_
-    (
-        IOobject
-        (
-            "coarse:" + mesh_.name(),
-            mesh_.polyMesh::instance(),
-            mesh_.time(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh_,
-        finalAgglom_
-    ),
+    coarseMesh_(),
+//     (
+//         IOobject
+//         (
+//             "coarse:" + mesh_.name(),
+//             mesh_.polyMesh::instance(),
+//             mesh_.time(),
+//             IOobject::NO_READ,
+//             IOobject::NO_WRITE,
+//             false
+//         ),
+//         mesh_,
+//         finalAgglom_
+//     ),
     qr_
     (
         IOobject
@@ -721,10 +752,10 @@ void Foam::radiation::viewFactor::calculate()
             const tmp<scalarField> tHoi = qrp.qro(bandI);
             const scalarField& Hoi = tHoi();
 
-            const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
+            const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
 
             const labelList& coarsePatchFace =
-                coarseMesh_.patchFaceMap()[patchID];
+                coarseMesh_->patchFaceMap()[patchID];
 
             scalarList T4ave(pp.size(), 0.0);
             scalarList Eave(pp.size(), 0.0);
@@ -1046,7 +1077,7 @@ void Foam::radiation::viewFactor::calculate()
     label globCoarseId = 0;
     for (const label patchID : selectedPatches_)
     {
-        const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
+        const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
 
         if (pp.size() > 0)
         {
@@ -1058,7 +1089,7 @@ void Foam::radiation::viewFactor::calculate()
             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 
             const labelList& coarsePatchFace =
-                coarseMesh_.patchFaceMap()[patchID];
+                coarseMesh_->patchFaceMap()[patchID];
 
             //scalar heatFlux = 0.0;
             forAll(coarseToFine, coarseI)
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H
index 662c6586be..ec685cc881 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -107,7 +107,7 @@ protected:
         autoPtr<IOmapDistribute> map_;
 
         //- Coarse mesh
-        singleCellFvMesh coarseMesh_;
+        autoPtr<singleCellFvMesh> coarseMesh_;
 
         //- Net radiative heat flux [W/m2]
         volScalarField qr_;

From 3b40ee8f6b877dd293f577b595ac8c403140150c Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 28 Jun 2022 16:04:59 -0700
Subject: [PATCH 05/10] ENH: Updating tutorials for view factors generation

---
 .../externalSolarLoad/Allrun                  | 10 ++--
 .../externalSolarLoad/Allrun-parallel         | 10 ++--
 .../constant/air/viewFactorsDict              |  3 +-
 .../multiRegionHeaterRadiation/Allrun         | 10 ++--
 .../Allrun-parallel                           | 10 ++--
 .../constant/bottomAir/viewFactorsDict        | 50 ------------------
 .../constant/topAir/viewFactorsDict           | 51 -------------------
 7 files changed, 21 insertions(+), 123 deletions(-)

diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun
index a851883e79..0010d50a53 100755
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun
@@ -9,11 +9,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 #-- Run on single processor
 
 # Agglomerate patch faces
-for region in air
-do
-    runApplication -s $region \
-        faceAgglomerate -region $region -dict constant/viewFactorsDict
-done
+#for region in air
+#do
+#    runApplication -s $region \
+#        faceAgglomerate -region $region -dict constant/viewFactorsDict
+#done
 
 # Generate view factors
 for region in air
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun-parallel b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun-parallel
index 53466d7bdf..840f6a9549 100755
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun-parallel
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/Allrun-parallel
@@ -12,11 +12,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 runApplication decomposePar -allRegions -constant
 
 # Agglomerate patch faces
-for region in air
-do
-    runParallel -s $region \
-        faceAgglomerate -region $region -dict constant/viewFactorsDict
-done
+#for region in air
+#do
+#    runParallel -s $region \
+#        faceAgglomerate -region $region -dict constant/viewFactorsDict
+#done
 
 # Generate view factors
 for region in air
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/viewFactorsDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/viewFactorsDict
index c990ed6136..4d1ce89a5c 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/viewFactorsDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalSolarLoad/constant/air/viewFactorsDict
@@ -15,11 +15,10 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 writeViewFactorMatrix     true;
-writeFacesAgglomeration   true;
 writePatchViewFactors     false;
 // dumpRays    true;
 
-maxDynListLength     200000;
+maxDynListLength          200000;
 
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun
index cc06292527..f7400fbe4f 100755
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun
@@ -9,11 +9,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 #-- Run on single processor
 
 # Agglomerate patch faces
-for region in bottomAir topAir
-do
-    runApplication -s "$region" \
-        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
-done
+#for region in bottomAir topAir
+#do
+#    runApplication -s "$region" \
+#        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
+#done
 
 # Generate view factors
 for region in bottomAir topAir
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun-parallel b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun-parallel
index 88d9626612..26ed334b6a 100755
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun-parallel
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/Allrun-parallel
@@ -12,11 +12,11 @@ cd "${0%/*}" || exit                                # Run from this directory
 runApplication decomposePar -allRegions
 
 # Agglomerate patch faces
-for region in bottomAir topAir
-do
-    runParallel -s "$region" \
-        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
-done
+#for region in bottomAir topAir
+#do
+#    runParallel -s "$region" \
+#        faceAgglomerate -region "$region" -dict constant/viewFactorsDict
+#done
 
 # Generate view factors
 for region in bottomAir topAir
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/bottomAir/viewFactorsDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/bottomAir/viewFactorsDict
index 863e9e4a9d..cc24d07533 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/bottomAir/viewFactorsDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/bottomAir/viewFactorsDict
@@ -15,56 +15,6 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 writeViewFactorMatrix     true;
-writeFacesAgglomeration   true;
-writePatchViewFactors     false;
-
-bottomAir_to_heater
-{
-    nFacesInCoarsestLevel     30;
-    featureAngle              10;
-}
-
-bottomAir_to_leftSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-bottomAir_to_rightSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-minX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-minY
-{
-    nFacesInCoarsestLevel     30;
-    featureAngle              10;
-}
-
-minZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-maxX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-maxZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
 
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/topAir/viewFactorsDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/topAir/viewFactorsDict
index ebd7dcafdc..c885191e52 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/topAir/viewFactorsDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation/constant/topAir/viewFactorsDict
@@ -15,56 +15,5 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 writeViewFactorMatrix     true;
-writeFacesAgglomeration   true;
-writePatchViewFactors     false;
-
-topAir_to_heater
-{
-    nFacesInCoarsestLevel     24;
-    featureAngle              10;
-}
-
-topAir_to_leftSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-topAir_to_rightSolid
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-minX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-maxY
-{
-    nFacesInCoarsestLevel     40;
-    featureAngle              10;
-}
-
-minZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
-maxX
-{
-    nFacesInCoarsestLevel     10;
-    featureAngle              10;
-}
-
-maxZ
-{
-    nFacesInCoarsestLevel     20;
-    featureAngle              10;
-}
-
 
 // ************************************************************************* //

From bfef08a89f83898667c146841be291f5d3be47f7 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Fri, 1 Jul 2022 14:41:21 -0700
Subject: [PATCH 06/10] ENH: Using globalIndex to create full triSurface

---
 .../viewFactorsGen/searchingEngine.H          |   2 +-
 .../viewFactorsGen/viewFactorsGen.C           | 127 ++++++++----------
 2 files changed, 55 insertions(+), 74 deletions(-)

diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index 9b89b848c9..667470d3f2 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -47,7 +47,7 @@ distributedTriSurfaceMesh surfacesMesh
         "triSurface",           // instance
         runTime,                // registry
         IOobject::NO_READ,
-        IOobject::AUTO_WRITE
+        IOobject::NO_WRITE
     ),
     localSurface,
     dict
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 368fc8df1c..232cef45fa 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -124,8 +124,6 @@ using namespace Foam;
 using namespace Foam::constant;
 using namespace Foam::constant::mathematical;
 
-//using namespace pbrt;
-
 triSurface triangulate
 (
     const polyBoundaryMesh& bMesh,
@@ -193,62 +191,45 @@ triSurface triangulate
         rawSurface.localPoints()
     );
 
-    // Combine the triSurfaces across all processors
-    if (Pstream::parRun())
+
+    globalIndex globalFaceIdx(surface.size(), globalIndex::gatherOnly());
+    globalIndex globalPointIdx
+    (
+        surface.points().size(), globalIndex::gatherOnly()
+    );
+
+    List<labelledTri> globalSurfaceTris(globalFaceIdx.gather(surface));
+    pointField globalSurfacePoints(globalPointIdx.gather(surface.points()));
+
+    //label offset = 0;
+    for (const label proci : globalPointIdx.allProcs())
     {
-        List<List<labelledTri>> surfaceProcTris(Pstream::nProcs());
-        List<pointField> surfaceProcPoints(Pstream::nProcs());
+        const label offset = globalPointIdx.localStart(proci);
 
-        surfaceProcTris[Pstream::myProcNo()] = surface;
-        surfaceProcPoints[Pstream::myProcNo()] = surface.points();
-
-        Pstream::gatherList(surfaceProcTris);
-        Pstream::scatterList(surfaceProcTris);
-        Pstream::gatherList(surfaceProcPoints);
-        Pstream::scatterList(surfaceProcPoints);
-
-        label nTris = 0;
-        forAll(surfaceProcTris, i)
+        if (offset)
         {
-            nTris += surfaceProcTris[i].size();
-        }
-
-        List<labelledTri> globalSurfaceTris(nTris);
-        label trii = 0;
-        label offset = 0;
-        forAll(surfaceProcTris, i)
-        {
-            forAll(surfaceProcTris[i], j)
+            for
+            (
+                labelledTri& tri
+             :  globalSurfaceTris.slice(globalFaceIdx.range(proci))
+            )
             {
-                globalSurfaceTris[trii] = surfaceProcTris[i][j];
-                globalSurfaceTris[trii][0] += offset;
-                globalSurfaceTris[trii][1] += offset;
-                globalSurfaceTris[trii][2] += offset;
-                trii++;
-            }
-            offset += surfaceProcPoints[i].size();
-        }
-
-        label nPoints = 0;
-        forAll(surfaceProcPoints, i)
-        {
-            nPoints += surfaceProcPoints[i].size();
-        }
-
-        pointField globalSurfacePoints(nPoints);
-
-        label pointi = 0;
-        forAll(surfaceProcPoints, i)
-        {
-            forAll(surfaceProcPoints[i], j)
-            {
-                globalSurfacePoints[pointi++] = surfaceProcPoints[i][j];
+                tri[0] += offset;
+                tri[1] += offset;
+                tri[2] += offset;
             }
         }
-
-        surface = triSurface(globalSurfaceTris, globalSurfacePoints);
     }
 
+    surface =
+        triSurface
+        (
+            std::move(globalSurfaceTris),
+            std::move(globalSurfacePoints)
+        );
+
+    Pstream::broadcast(surface);
+
     // Add patch names to surface
     surface.patches().setSize(newPatchI);
 
@@ -857,21 +838,21 @@ int main(int argc, char *argv[])
 
     GaussPoints[2].setSize(3);
     GaussPoints[2][0] =  0;
-    GaussPoints[2][1] =  0.774597;
-    GaussPoints[2][2] = -0.774597;
+    GaussPoints[2][1] =  std::sqrt(3/5);
+    GaussPoints[2][2] = -std::sqrt(3/5);
 
     GaussPoints[3].setSize(4);
-    GaussPoints[3][0] =  0.339981;
-    GaussPoints[3][1] = -0.339981;
-    GaussPoints[3][2] = 0.861136;
-    GaussPoints[3][3] = -0.861136;
+    GaussPoints[3][0] = std::sqrt(3.0/7.0 - (2.0/7.0)*std::sqrt(6/5));
+    GaussPoints[3][1] = -GaussPoints[3][0];
+    GaussPoints[3][2] = std::sqrt(3.0/7.0 + (2.0/7.0)*std::sqrt(6/5));
+    GaussPoints[3][3] = -GaussPoints[3][2];
 
     GaussPoints[4].setSize(5);
     GaussPoints[4][0] =  0;
-    GaussPoints[4][1] = 0.538469;
-    GaussPoints[4][2] = -0.538469;
-    GaussPoints[4][3] = 0.90618;
-    GaussPoints[4][4] = -0.90618;
+    GaussPoints[4][1] = (1.0/3.0)*std::sqrt(5.0 - 2.0*std::sqrt(10/7));
+    GaussPoints[4][2] = -GaussPoints[4][1];
+    GaussPoints[4][3] = (1.0/3.0)*std::sqrt(5.0 + 2.0*std::sqrt(10/7));
+    GaussPoints[4][4] = -GaussPoints[4][3];
 
 
     List<scalarList> GaussWeights(5);
@@ -883,24 +864,24 @@ int main(int argc, char *argv[])
     GaussWeights[1][1] =  1;
 
     GaussWeights[2].setSize(3);
-    GaussWeights[2][0] =  0.888889;
-    GaussWeights[2][1] =  0.555556;
-    GaussWeights[2][2] =  0.555556;
+    GaussWeights[2][0] =  8.0/9.0;
+    GaussWeights[2][1] =  5.0/9.0;
+    GaussWeights[2][2] =  5.0/9.0;
 
     GaussWeights[3].setSize(4);
-    GaussWeights[3][0] =  0.652145;
-    GaussWeights[3][1] =  0.652145;
-    GaussWeights[3][2] =  0.347855;
-    GaussWeights[3][3] =  0.347855;
+    GaussWeights[3][0] =  (18.0 + std::sqrt(30))/36.0;
+    GaussWeights[3][1] =  (18.0 + std::sqrt(30))/36.0;
+    GaussWeights[3][2] =  (18.0 - std::sqrt(30))/36.0;
+    GaussWeights[3][3] =  (18.0 - std::sqrt(30))/36.0;
 
     GaussWeights[4].setSize(5);
-    GaussWeights[4][0] =  0.568889;
-    GaussWeights[4][1] =  0.478629;
-    GaussWeights[4][2] =  0.478629;
-    GaussWeights[4][3] =  0.236927;
-    GaussWeights[4][4] =  0.236927;
+    GaussWeights[4][0] =  128.0/225.0;
+    GaussWeights[4][1] =  (322.0 + 13.0*std::sqrt(70))/900.0;
+    GaussWeights[4][2] =  (322.0 + 13.0*std::sqrt(70))/900.0;
+    GaussWeights[4][3] =  (322.0 - 13.0*std::sqrt(70))/900.0;
+    GaussWeights[4][4] =  (322.0 - 13.0*std::sqrt(70))/900.0;
 
-    label maxQuadOrder = 5;
+    const label maxQuadOrder = 5;
 
     if (mesh.nSolutionD() == 3)
     {

From c652af4b82f0da4557499878735527ffd5d81b95 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Tue, 5 Jul 2022 10:15:05 -0700
Subject: [PATCH 07/10] STY: Minor style changes

---
 .../preProcessing/viewFactorsGen/viewFactorsGen.C    | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 232cef45fa..d913f7bf09 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -838,20 +838,20 @@ int main(int argc, char *argv[])
 
     GaussPoints[2].setSize(3);
     GaussPoints[2][0] =  0;
-    GaussPoints[2][1] =  std::sqrt(3/5);
-    GaussPoints[2][2] = -std::sqrt(3/5);
+    GaussPoints[2][1] =  std::sqrt(3.0/5.0);
+    GaussPoints[2][2] = -std::sqrt(3.0/5.0);
 
     GaussPoints[3].setSize(4);
-    GaussPoints[3][0] = std::sqrt(3.0/7.0 - (2.0/7.0)*std::sqrt(6/5));
+    GaussPoints[3][0] = std::sqrt(3.0/7.0 - (2.0/7.0)*std::sqrt(6.0/5.0));
     GaussPoints[3][1] = -GaussPoints[3][0];
-    GaussPoints[3][2] = std::sqrt(3.0/7.0 + (2.0/7.0)*std::sqrt(6/5));
+    GaussPoints[3][2] = std::sqrt(3.0/7.0 + (2.0/7.0)*std::sqrt(6.0/5.0));
     GaussPoints[3][3] = -GaussPoints[3][2];
 
     GaussPoints[4].setSize(5);
     GaussPoints[4][0] =  0;
-    GaussPoints[4][1] = (1.0/3.0)*std::sqrt(5.0 - 2.0*std::sqrt(10/7));
+    GaussPoints[4][1] = (1.0/3.0)*std::sqrt(5.0 - 2.0*std::sqrt(10.0/7.0));
     GaussPoints[4][2] = -GaussPoints[4][1];
-    GaussPoints[4][3] = (1.0/3.0)*std::sqrt(5.0 + 2.0*std::sqrt(10/7));
+    GaussPoints[4][3] = (1.0/3.0)*std::sqrt(5.0 + 2.0*std::sqrt(10.0/7.0));
     GaussPoints[4][4] = -GaussPoints[4][3];
 
 

From f2f71f6847b430a3c5e97e8fa65e60ce33488ac2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 20 Jul 2022 14:55:12 +0100
Subject: [PATCH 08/10] ENH: viewFactor: compile without lib

---
 .../utilities/preProcessing/viewFactorsGen/Make/options        | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/options b/applications/utilities/preProcessing/viewFactorsGen/Make/options
index 298501f447..6ac3efad51 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@@ -3,6 +3,7 @@ include $(GENERAL_RULES)/cgal
 EXE_INC = \
     -Wno-old-style-cast \
     ${CGAL_INC} \
+    -DCGAL_HEADER_ONLY \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
@@ -10,7 +11,7 @@ EXE_INC = \
 
 
 EXE_LIBS = \
-    ${CGAL_LIBS} \
+    /* ${CGAL_LIBS} */ \
     -lfiniteVolume \
     -lsurfMesh \
     -lmeshTools \

From 994addc5430c9337f898281702934629649206c3 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Fri, 22 Jul 2022 09:36:31 -0700
Subject: [PATCH 09/10] ENH: Changing key entry name to GaussQuadTol

---
 .../preProcessing/viewFactorsGen/viewFactorsGen.C        | 9 ++++-----
 .../radiation/radiationModels/viewFactor/viewFactor.C    | 8 ++++----
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index d913f7bf09..264f141cd2 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -40,7 +40,7 @@ Description
 
     The input from viewFactorsDict are:
 
-        tolGaussQuad              0.1;      // GaussQuad  error
+        GaussQuadTol              0.1;      // GaussQuad  error
         distTol                   8;        // R/Average(rm)
         alpha                     0.22;     // Use for common edges for 2LI
 
@@ -80,7 +80,6 @@ Description
 
 #include "indirectPrimitivePatch.H"
 #include "DynamicField.H"
-//#include "unitConversion.H"
 
 #include "scalarMatrices.H"
 #include "labelListIOList.H"
@@ -442,8 +441,8 @@ int main(int argc, char *argv[])
 
     const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
 
-    const scalar tolGaussQuad =
-        viewFactorDict.getOrDefault<scalar>("tolGaussQuad", 0.01);
+    const scalar GaussQuadTol =
+        viewFactorDict.getOrDefault<scalar>("GaussQuadTol", 0.01);
 
     const scalar distTol =
          viewFactorDict.getOrDefault<scalar>("distTol", 8);
@@ -1050,7 +1049,7 @@ int main(int argc, char *argv[])
 
                         if
                         (
-                            (mag(err) < tolGaussQuad && gi > 0)
+                            (mag(err) < GaussQuadTol && gi > 0)
                          || gi == maxQuadOrder-1
                         )
                         {
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 4d8dfb3c87..3c811ef202 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -219,7 +219,7 @@ void Foam::radiation::viewFactor::initialise()
         labelList upper(rays_.size(), -1);
         labelList lower(rays_.size(), -1);
 
-        const edgeList& raysLst = rays_.sortedToc();
+        const edgeList raysLst(rays_.sortedToc());
         label rayI = 0;
         for (const auto& e : raysLst)
         {
@@ -441,8 +441,8 @@ void Foam::radiation::viewFactor::initialise()
             totalDelta /= myF.size();
             reduce(totalDelta, sumOp<scalar>());
             reduce(maxDelta, maxOp<scalar>());
-            Info << "Smoothng average delta : " << totalDelta << endl;
-            Info << "Smoothng maximum delta : " << maxDelta << nl << endl;
+            Info << "Smoothing average delta : " << totalDelta << endl;
+            Info << "Smoothing maximum delta : " << maxDelta << nl << endl;
         }
     }
 
@@ -853,7 +853,7 @@ void Foam::radiation::viewFactor::calculate()
             // Local matrix coefficients
             if (!constEmissivity_ || iterCounter_ == 0)
             {
-                const edgeList& raysLst = rays_.sortedToc();
+                const edgeList raysLst(rays_.sortedToc());
 
                 label rayI = 0;
                 for (const auto& e : raysLst)

From c08fc5ecd9901dfb1df1fddc28a2107925c34f37 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 4 Aug 2022 14:24:50 +0100
Subject: [PATCH 10/10] ENH: viewFactorGen: revert to v2206 without CGAL

---
 .../preProcessing/viewFactorsGen/Allwmake     |  19 +++
 .../preProcessing/viewFactorsGen/Make/options |   3 +-
 .../viewFactorsGen/searchingEngine.H          |  31 +---
 .../viewFactorsGen/searchingEngine_CGAL.H     |  84 ++++++++++
 .../preProcessing/viewFactorsGen/shootRays.H  | 147 +++++++++++++-----
 .../viewFactorsGen/shootRays_CGAL.H           | 131 ++++++++++++++++
 .../viewFactorsGen/viewFactorsGen.C           |  41 +++--
 7 files changed, 371 insertions(+), 85 deletions(-)
 create mode 100755 applications/utilities/preProcessing/viewFactorsGen/Allwmake
 create mode 100644 applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
 create mode 100644 applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H

diff --git a/applications/utilities/preProcessing/viewFactorsGen/Allwmake b/applications/utilities/preProcessing/viewFactorsGen/Allwmake
new file mode 100755
index 0000000000..f0d604d4e3
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGen/Allwmake
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments # (error catching)
+. ${WM_PROJECT_DIR:?}/wmake/scripts/have_cgal
+
+#------------------------------------------------------------------------------
+unset COMP_FLAGS LINK_FLAGS
+
+if have_cgal
+then
+    echo "    found CGAL -- enabling CGAL support."
+else
+    echo "    did not find CGAL -- disabling CGAL functionality"
+    export COMP_FLAGS="-DNO_CGAL"
+fi
+
+wmake $targetType
+
+#------------------------------------------------------------------------------
diff --git a/applications/utilities/preProcessing/viewFactorsGen/Make/options b/applications/utilities/preProcessing/viewFactorsGen/Make/options
index 6ac3efad51..08e64c8a7d 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@@ -1,7 +1,8 @@
-include $(GENERAL_RULES)/cgal
+include $(GENERAL_RULES)/cgal-header-only
 
 EXE_INC = \
     -Wno-old-style-cast \
+    $(COMP_FLAGS) \
     ${CGAL_INC} \
     -DCGAL_HEADER_ONLY \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index 667470d3f2..df38326ded 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -25,7 +25,7 @@ dict.add("mergeDistance", SMALL);
 
 labelList triSurfaceToAgglom(5*nFineFaces);
 
-triSurface localSurface = triangulate
+const triSurface localSurface = triangulate
 (
     patches,
     includePatches,
@@ -36,8 +36,6 @@ triSurface localSurface = triangulate
 );
 
 
-// CGAL surface
-
 distributedTriSurfaceMesh surfacesMesh
 (
     IOobject
@@ -53,32 +51,7 @@ distributedTriSurfaceMesh surfacesMesh
     dict
 );
 
+
 triSurfaceToAgglom.resize(surfacesMesh.size());
 
 surfacesMesh.setField(triSurfaceToAgglom);
-
-const pointField& pts = surfacesMesh.localPoints();
-
-std::list<Triangle> triangles;
-
-for (const auto& triLocal : surfacesMesh.localFaces())
-{
-    point p1l = pts[triLocal[0]];
-    point p2l = pts[triLocal[1]];
-    point p3l = pts[triLocal[2]];
-
-    Point p1(p1l[0], p1l[1], p1l[2]);
-    Point p2(p2l[0], p2l[1], p2l[2]);
-    Point p3(p3l[0], p3l[1], p3l[2]);
-
-    Triangle tri(p1, p2, p3);
-
-    if (tri.is_degenerate())
-    {
-        std::cout << tri << std::endl;
-    }
-    triangles.push_back(tri);
-}
-
-// constructs AABB tree
-Tree tree(triangles.begin(), triangles.end());
diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
new file mode 100644
index 0000000000..667470d3f2
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
@@ -0,0 +1,84 @@
+Random rndGen(653213);
+
+// Determine mesh bounding boxes:
+List<treeBoundBox> meshBb
+(
+    1,
+    treeBoundBox
+    (
+        boundBox(coarseMesh.points(), false)
+    ).extend(rndGen, 1e-3)
+);
+
+// Dummy bounds dictionary
+dictionary dict;
+dict.add("bounds", meshBb);
+dict.add
+(
+    "distributionType",
+    distributedTriSurfaceMesh::distributionTypeNames_
+    [
+        distributedTriSurfaceMesh::FROZEN
+    ]
+);
+dict.add("mergeDistance", SMALL);
+
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+triSurface localSurface = triangulate
+(
+    patches,
+    includePatches,
+    finalAgglom,
+    triSurfaceToAgglom,
+    globalNumbering,
+    coarsePatches
+);
+
+
+// CGAL surface
+
+distributedTriSurfaceMesh surfacesMesh
+(
+    IOobject
+    (
+        "wallSurface.stl",
+        runTime.constant(),     // directory
+        "triSurface",           // instance
+        runTime,                // registry
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    localSurface,
+    dict
+);
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
+surfacesMesh.setField(triSurfaceToAgglom);
+
+const pointField& pts = surfacesMesh.localPoints();
+
+std::list<Triangle> triangles;
+
+for (const auto& triLocal : surfacesMesh.localFaces())
+{
+    point p1l = pts[triLocal[0]];
+    point p2l = pts[triLocal[1]];
+    point p3l = pts[triLocal[2]];
+
+    Point p1(p1l[0], p1l[1], p1l[2]);
+    Point p2(p2l[0], p2l[1], p2l[2]);
+    Point p3(p3l[0], p3l[1], p3l[2]);
+
+    Triangle tri(p1, p2, p3);
+
+    if (tri.is_degenerate())
+    {
+        std::cout << tri << std::endl;
+    }
+    triangles.push_back(tri);
+}
+
+// constructs AABB tree
+Tree tree(triangles.begin(), triangles.end());
diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
index c738784b79..a5f808c392 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -1,23 +1,35 @@
+// All rays expressed as start face (local) index end end face (global)
+// Pre-size by assuming a certain percentage is visible.
+
+
+if (Pstream::master())
+{
+    Info<< "\nShooting rays using distributedTriSurfaceMesh (no CGAL support)"
+        << endl;
+}
+
+
+
 // Maximum length for dynamicList
 const label maxDynListLength
 (
-    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000)
 );
 
 for (const int proci : Pstream::allProcs())
 {
-    std::vector<Point> start;
-    start.reserve(coarseMesh.nFaces());
-
-    std::vector<Point> end(start.size());
-    end.reserve(start.size());
+    // Shoot rays from me to proci. Note that even if processor has
+    // 0 faces we still need to call findLine to keep calls synced.
 
+    DynamicField<point> start(coarseMesh.nFaces());
+    DynamicField<point> end(start.size());
     DynamicList<label> startIndex(start.size());
     DynamicList<label> endIndex(start.size());
 
     DynamicList<label> startAgg(start.size());
     DynamicList<label> endAgg(start.size());
 
+
     const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
     const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
     const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
@@ -46,31 +58,15 @@ for (const int proci : Pstream::allProcs())
 
                     const vector d(remFc - fc);
 
-
-                    const vector nd = d/mag(d);
-                    const vector nfA = fA/mag(fA);
-                    const vector nremA = remA/mag(remA);
-
-                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
+                    if (((d & fA) < 0.) && ((d & remA) > 0))
                     {
-                        vector direction(d[0], d[1], d[2]);
-                        vector s(fc[0], fc[1], fc[2]);
-                        vector rayEnd(s + (1-intTol)*direction);
-                        end.push_back(Point(rayEnd[0], rayEnd[1], rayEnd[2]));
-
-                        s += vector(intTol*d[0], intTol*d[1], intTol*d[2]);
-
-                        start.push_back(Point(s[0], s[1], s[2]));
+                        start.append(fc + 0.001*d);
                         startIndex.append(i);
-                        if (useAgglomeration)
-                        {
-                            startAgg.append(globalNumbering.toGlobal(proci, fAgg));
-                            endAgg.append(globalNumbering.toGlobal(proci, remAgg));
-                        }
-
+                        startAgg.append(globalNumbering.toGlobal(proci, fAgg));
+                        end.append(fc + 0.999*d);
                         label globalI = globalNumbering.toGlobal(proci, j);
                         endIndex.append(globalI);
-
+                        endAgg.append(globalNumbering.toGlobal(proci, remAgg));
                         if (startIndex.size() > maxDynListLength)
                         {
                             FatalErrorInFunction
@@ -89,33 +85,98 @@ for (const int proci : Pstream::allProcs())
             }
         }
 
-    } while (i < myFc.size());
+    } while (returnReduce(i < myFc.size(), orOp<bool>()));
 
-    label totalRays(startIndex.size());
-    reduce(totalRays, sumOp<scalar>());
+    List<pointIndexHit> hitInfo(startIndex.size());
+    surfacesMesh.findLine(start, end, hitInfo);
 
-    if (debug)
+    // Return hit global agglo index
+    labelList aggHitIndex;
+    surfacesMesh.getField(hitInfo, aggHitIndex);
+
+    DynamicList<label> dRayIs;
+
+    // Collect the rays which have no obstacle in between rayStartFace
+    // and rayEndFace. If the ray hit itself, it gets stored in dRayIs
+    forAll(hitInfo, rayI)
     {
-        Pout<< "Number of rays :" << totalRays << endl;
-    }
-
-    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
-    {
-        Segment ray(start[rayI], end[rayI]);
-        Segment_intersection intersects = tree.any_intersection(ray);
-
-        if (!intersects)
+        if (!hitInfo[rayI].hit())
         {
             rayStartFace.append(startIndex[rayI]);
             rayEndFace.append(endIndex[rayI]);
         }
+        else if (aggHitIndex[rayI] == startAgg[rayI])
+        {
+            dRayIs.append(rayI);
+        }
     }
+
     start.clear();
 
-    if (debug)
+
+    // Continue rays which hit themself. If they hit the target
+    // agglomeration are added to rayStartFace and rayEndFace
+
+    bool firstLoop = true;
+    DynamicField<point> startHitItself;
+    DynamicField<point> endHitItself;
+    label iter = 0;
+
+    do
     {
-        Pout << "hits : " <<  rayStartFace.size()<< endl;
-    }
+        labelField rayIs;
+        rayIs.transfer(dRayIs);
+        dRayIs.clear();
+        forAll(rayIs, rayI)
+        {
+            const label rayID = rayIs[rayI];
+            label hitIndex = -1;
+
+            if (firstLoop)
+            {
+                hitIndex = rayIs[rayI];
+            }
+            else
+            {
+                hitIndex = rayI;
+            }
+
+            if (hitInfo[hitIndex].hit())
+            {
+                if (aggHitIndex[hitIndex] == startAgg[rayID])
+                {
+                    const vector& endP = end[rayID];
+                    const vector& startP = hitInfo[hitIndex].hitPoint();
+                    const vector& d = endP - startP;
+
+                    startHitItself.append(startP + 0.01*d);
+                    endHitItself.append(startP + 1.01*d);
+
+                    dRayIs.append(rayID);
+                }
+                else if (aggHitIndex[hitIndex] == endAgg[rayID])
+                {
+                    rayStartFace.append(startIndex[rayID]);
+                    rayEndFace.append(endIndex[rayID]);
+                }
+
+            }
+        }
+
+        hitInfo.clear();
+        hitInfo.resize(dRayIs.size());
+
+        surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
+
+        surfacesMesh.getField(hitInfo, aggHitIndex);
+
+
+        endHitItself.clear();
+        startHitItself.clear();
+        firstLoop = false;
+        iter ++;
+
+    } while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
 
     startIndex.clear();
     end.clear();
diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
new file mode 100644
index 0000000000..8be6e00c57
--- /dev/null
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
@@ -0,0 +1,131 @@
+// Maximum length for dynamicList
+const label maxDynListLength
+(
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+);
+
+for (const int proci : Pstream::allProcs())
+{
+    std::vector<Point> start;
+    start.reserve(coarseMesh.nFaces());
+
+    std::vector<Point> end(start.size());
+    end.reserve(start.size());
+
+    DynamicList<label> startIndex(start.size());
+    DynamicList<label> endIndex(start.size());
+
+    DynamicList<label> startAgg(start.size());
+    DynamicList<label> endAgg(start.size());
+
+    const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
+    const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
+    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
+
+    const pointField& remoteArea = remoteCoarseSf[proci];
+    const pointField& remoteFc = remoteCoarseCf[proci];
+    const labelField& remoteAgg = remoteCoarseAgg[proci];
+
+    label i = 0;
+    label j = 0;
+    do
+    {
+        for (; i < myFc.size(); i++)
+        {
+            const point& fc = myFc[i];
+            const vector& fA = myArea[i];
+            const label& fAgg = myAgg[i];
+
+            for (; j < remoteFc.size(); j++)
+            {
+                if (proci != Pstream::myProcNo() || i != j)
+                {
+                    const point& remFc = remoteFc[j];
+                    const vector& remA = remoteArea[j];
+                    const label& remAgg = remoteAgg[j];
+
+                    const vector d(remFc - fc);
+
+
+                    const vector nd = d/mag(d);
+                    const vector nfA = fA/mag(fA);
+                    const vector nremA = remA/mag(remA);
+
+                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
+                    {
+                        vector direction(d[0], d[1], d[2]);
+                        vector s(fc[0], fc[1], fc[2]);
+                        vector rayEnd(s + (1-intTol)*direction);
+                        end.push_back(Point(rayEnd[0], rayEnd[1], rayEnd[2]));
+
+                        s += vector(intTol*d[0], intTol*d[1], intTol*d[2]);
+
+                        start.push_back(Point(s[0], s[1], s[2]));
+                        startIndex.append(i);
+                        if (useAgglomeration)
+                        {
+                            startAgg.append
+                            (
+                                globalNumbering.toGlobal(proci, fAgg)
+                            );
+                            endAgg.append
+                            (
+                                globalNumbering.toGlobal(proci, remAgg)
+                            );
+                        }
+
+                        label globalI = globalNumbering.toGlobal(proci, j);
+                        endIndex.append(globalI);
+
+                        if (startIndex.size() > maxDynListLength)
+                        {
+                            FatalErrorInFunction
+                                << "Dynamic list need from capacity."
+                                << "Actual size maxDynListLength : "
+                                <<  maxDynListLength
+                                << abort(FatalError);
+                        }
+                    }
+                }
+            }
+
+            if (j == remoteFc.size())
+            {
+                j = 0;
+            }
+        }
+
+    } while (i < myFc.size());
+
+    label totalRays(startIndex.size());
+    reduce(totalRays, sumOp<scalar>());
+
+    if (debug)
+    {
+        Pout<< "Number of rays :" << totalRays << endl;
+    }
+
+    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
+    {
+        Segment ray(start[rayI], end[rayI]);
+        Segment_intersection intersects = tree.any_intersection(ray);
+
+        if (!intersects)
+        {
+            rayStartFace.append(startIndex[rayI]);
+            rayEndFace.append(endIndex[rayI]);
+        }
+    }
+    start.clear();
+
+    if (debug)
+    {
+        Pout << "hits : " <<  rayStartFace.size()<< endl;
+    }
+
+    startIndex.clear();
+    end.clear();
+    endIndex.clear();
+    startAgg.clear();
+    endAgg.clear();
+}
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 264f141cd2..d34a5eb0ac 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -4,7 +4,7 @@
    \\    /   O peration     |
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
-// -------------------------------------------------------------------------------
+-------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
     Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
@@ -25,7 +25,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    viewFactorsGenExt
+    viewFactorsGen
 
 Group
     grpPreProcessingUtilities
@@ -33,7 +33,7 @@ Group
 Description
     This view factors generation application uses a combined approach of
     double area integral (2AI) and double linear integral (2LI). 2AI is used
-    when the two surfaces are 'far' apart and 2LI whenre they are 'close'.
+    when the two surfaces are 'far' apart and 2LI when they are 'close'.
     2LI is integrated along edges using Gaussian quadrature.
     The distance between faces is calculating a ratio between averaged areas
     and the distance between face centres.
@@ -57,8 +57,8 @@ Description
         dumpRays                dumps rays
 
 
-    In order to specify the participants patches in the VF calculation the
-    keywaord viewFactorWall should be added to the boundary file.
+    The participating patches in the VF calculation have to be in the
+    'viewFactorWall' patch group (in the polyMesh/boundary file), e.g.
 
     floor
     {
@@ -68,6 +68,9 @@ Description
         startFace       3100;
     }
 
+    Compile with -DNO_CGAL only if no CGAL present - CGAL AABB tree performs
+    better than the built-in octree.
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
@@ -181,15 +184,13 @@ triSurface triangulate
     //triSurfaceToAgglom.resize(localTriFaceI-1);
 
     triangles.shrink();
+    triSurface surface(triangles, mesh.points());
+    surface.compactPoints();
 
-    triSurface rawSurface(triangles, mesh.points());
 
-    triSurface surface
-    (
-        rawSurface.localFaces(),
-        rawSurface.localPoints()
-    );
+#ifndef NO_CGAL
 
+    // CGAL : every processor has whole surface
 
     globalIndex globalFaceIdx(surface.size(), globalIndex::gatherOnly());
     globalIndex globalPointIdx
@@ -228,6 +229,7 @@ triSurface triangulate
         );
 
     Pstream::broadcast(surface);
+#endif
 
     // Add patch names to surface
     surface.patches().setSize(newPatchI);
@@ -337,6 +339,7 @@ void insertMatrixElements
     }
 }
 
+
 scalar GaussQuad
 (
 
@@ -403,6 +406,8 @@ scalar GaussQuad
     }
     return dIntFij;
 }
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -643,7 +648,13 @@ int main(int argc, char *argv[])
 
     // Set up searching engine for obstacles
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#ifdef NO_CGAL
+    // Using octree
     #include "searchingEngine.H"
+#else
+    // Using CGAL aabbtree (faster, more robust)
+    #include "searchingEngine_CGAL.H"
+#endif
 
     // Determine rays between coarse face centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -653,7 +664,13 @@ int main(int argc, char *argv[])
 
     // Return rayStartFace in local index and rayEndFace in global index
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#ifdef NO_CGAL
+    // Using octree, distributedTriSurfaceMesh
     #include "shootRays.H"
+#else
+    // Using CGAL aabbtree (faster, more robust)
+    #include "shootRays_CGAL.H"
+#endif
 
     // Calculate number of visible faces from local index
     labelList nVisibleFaceFaces(nCoarseFaces, Zero);
@@ -1135,7 +1152,7 @@ int main(int argc, char *argv[])
 
                 forAll(coarseToFine, coarseI)
                 {
-                    scalar FiSum = sum(F2LI[compactI]);
+                    const scalar FiSum = sum(F2LI[compactI]);
 
                     const label coarseFaceID = coarsePatchFace[coarseI];
                     const labelList& fineFaces = coarseToFine[coarseFaceID];