ENH: Adding new viewFactor utility using CGAL

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 \

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());

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
+    std::vector<Point> start;
-    // 0 faces we still need to call findLine to keep calls synced.
+    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));
+                        if (useAgglomeration)
-                        end.append(fc + 0.999*d);
+                        {
+                            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());
+    label totalRays(startIndex.size());
-    surfacesMesh.findLine(start, end, hitInfo);
+    reduce(totalRays, sumOp<scalar>());
 
-    // Return hit global agglo index
+    if (debug)
-    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 (!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();
 
-
+    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
     {
-        labelField rayIs;
+        Pout << "hits : " <<  rayStartFace.size()<< endl;
-        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();

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
+    This view factors generation application uses a combined approach of
-    (finalAgglom generated by faceAgglomerate utility).
+    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
+    The input from viewFactorsDict are:
-    using a double integral of the sub-areas composing the agglomeration.
 
-    The patches involved in the view factor calculation are taken from the
+        tolGaussQuad              0.1;      // GaussQuad  error
-    boundary file and should be part on the group viewFactorWall. ie.:
+        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 F2LI
-    scalarListIOList F
     (
         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
+                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
-                        (
-                            dCi,
-                            dCj,
-                            dAi,
-                            dAj
-                        );
 
-                        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();
+         FatalErrorInFunction
-        const Vector<label>& dirs = mesh.geometricD();
+            << " View factors are not available in 2D "
-        vector emptyDir = Zero;
+            << exit(FatalError);
-        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;
     }
 
     // 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++;
                 }

applications/utilities/preProcessing/viewFactorsGenExt/Make/files

@@ -1,3 +0,0 @@
-viewFactorsGenExt.C
-
-EXE = $(FOAM_APPBIN)/viewFactorsGenExt

applications/utilities/preProcessing/viewFactorsGenExt/Make/options

@@ -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

applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H

@@ -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);

applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H

@@ -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();
-}