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];