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] 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;
-}
-
-
-// ************************************************************************* //