ENH: Adding viewFactorsGenExt which uses pbrt for ray tracing

applications/utilities/preProcessing/viewFactorsGen/shootRays.H

@@ -49,6 +49,7 @@ for (const int proci : Pstream::allProcs())
 
                     const vector d(remFc - fc);
 
+
                     if (((d & fA) < 0.) && ((d & remA) > 0))
                     {
                         start.append(fc + 0.001*d);

applications/utilities/preProcessing/viewFactorsGenExt/Make/files

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

applications/utilities/preProcessing/viewFactorsGenExt/Make/options

@@ -0,0 +1,39 @@
+PBRT = /home/sergio/pub/pbrt
+PBRTSRC = $(PBRT)/pbrt-v3/src
+EXTSRC = $(PBRT)/pbrt-v3/src/ext
+GLOGSRC = $(EXTSRC)/glog/src
+
+PBRTBUILD = /home/sergio/pub/pbrt/pbrt-v3-build
+EXTBUILD = $(PBRTBUILD)/src/ext
+GLOGBUILD = $(EXTBUILD)/glog
+
+EXE_INC = \
+    -Wno-old-style-cast \
+    -I$(PBRTSRC)/core \
+    -I$(PBRTSRC)/accelerators \
+    -I$(PBRTSRC)/shapes \
+    -I$(PBRTSRC) \
+    -I$(EXTSRC) \
+    -I$(GLOGSRC) \
+    -I$(GLOGBUILD) \
+    -I./rays \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/distributed/lnInclude \
+
+EXE_LIBS = \
+    -L$(GLOGBUILD) -lglog \
+    -L$(PBRTBUILD) -lpbrt \
+    -L$(EXTBUILD)/openexr/OpenEXR/IlmImf -lIlmImf \
+    -L$(EXTBUILD)/openexr/IlmBase/Imath -lImath \
+    -L$(EXTBUILD)/openexr/IlmBase/Half -lHalf \
+    -L$(EXTBUILD)/openexr/IlmBase/Iex -lIex \
+    -L$(EXTBUILD)/openexr/IlmBase/IexMath -lIexMath \
+    -L$(EXTBUILD)/openexr/IlmBase/IlmThread -lIlmThread \
+    -L$(EXTBUILD)/ptex/src/ptex -lPtex \
+    -lfiniteVolume \
+    -lsurfMesh \
+    -lmeshTools \
+    -ldistributed \
+    -lradiationModels

applications/utilities/preProcessing/viewFactorsGenExt/searchingEngineExt.H

@@ -0,0 +1,121 @@
+Random rndGen(653213);
+
+// Determine mesh bounding boxes:
+List<treeBoundBox> meshBb
+(
+    1,
+    treeBoundBox
+    (
+        boundBox(coarseMesh.points(), false)
+    ).extend(rndGen, 1e-3)
+);
+
+// Dummy bounds dictionary
+dictionary dict;
+dict.add("bounds", meshBb);
+dict.add
+(
+    "distributionType",
+    distributedTriSurfaceMesh::distributionTypeNames_
+    [
+        distributedTriSurfaceMesh::FROZEN
+    ]
+);
+dict.add("mergeDistance", SMALL);
+
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+const triSurface localSurface = triangulate
+(
+    patches,
+    includePatches,
+    finalAgglom,
+    triSurfaceToAgglom,
+    globalNumbering,
+    coarsePatches
+);
+
+
+distributedTriSurfaceMesh surfacesMesh
+(
+    IOobject
+    (
+        "wallSurface.stl",
+        runTime.constant(),     // directory
+        "triSurface",           // instance
+        runTime,                // registry
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    localSurface,
+    dict
+);
+
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
+// pbrt surface
+std::vector<Point3f> vertices;
+
+const pointField& pts = surfacesMesh.localPoints();
+for (const auto& pt : pts)
+{
+    vertices.push_back(Point3f(pt[0], pt[1], pt[2]));
+}
+
+std::vector<int> indices;
+for (const auto& tri : surfacesMesh.localFaces())
+{
+    indices.push_back(tri[0]);
+    indices.push_back(tri[1]);
+    indices.push_back(tri[2]);
+}
+
+std::vector<int> faceIndices;
+faceIndices.reserve(surfacesMesh.localFaces().size());
+
+for (label faceI = 0; faceI < surfacesMesh.localFaces().size(); faceI++)
+{
+    faceIndices.push_back(faceI);
+}
+
+Transform o2w;
+Transform w2o;
+
+std::vector<std::shared_ptr<Shape>> surfacesMesh_pbrt
+(
+    CreateTriangleMesh
+    (
+        &o2w,
+        &w2o,
+        false,                  // bool reverseOrientation
+        surfacesMesh.size(),    // int nTriangles
+        &indices[0],            // int *vertexIndices
+        vertices.size(),
+        &vertices[0],
+        nullptr,
+        nullptr,
+        nullptr,
+        nullptr,
+        nullptr,
+        &faceIndices[0]
+    )
+);
+
+std::vector<std::shared_ptr<Primitive>> prims;
+
+prims.reserve(surfacesMesh_pbrt.size());
+for (auto s : surfacesMesh_pbrt)
+{
+    prims.push_back
+    (
+        std::make_shared<GeometricPrimitive>
+            (s, nullptr, nullptr, nullptr)
+    );
+}
+
+std::shared_ptr<Primitive> accel;
+
+ParamSet paramSet;
+
+accel = CreateBVHAccelerator(std::move(prims), paramSet);

applications/utilities/preProcessing/viewFactorsGenExt/shootRaysExt.H

@@ -0,0 +1,196 @@
+// All rays expressed as start face (local) index end end face (global)
+// Pre-size by assuming a certain percentage is visible.
+
+// Maximum length for dynamicList
+const label maxDynListLength
+(
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
+);
+
+for (const int proci : Pstream::allProcs())
+{
+
+    std::vector<pbrt::Point3f> start;
+    start.reserve(coarseMesh.nFaces());
+    std::vector<pbrt::Vector3f> dir;
+    dir.reserve(start.size());
+
+    std::vector<pbrt::Point3f> end(start.size());
+    end.reserve(start.size());
+
+    DynamicList<label> startIndex(start.size());
+    DynamicList<label> endIndex(start.size());
+
+    DynamicList<label> startAgg(start.size());
+    DynamicList<label> endAgg(start.size());
+
+    const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
+    const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
+    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
+
+    const pointField& remoteArea = remoteCoarseSf[proci];
+    const pointField& remoteFc = remoteCoarseCf[proci];
+    const labelField& remoteAgg = remoteCoarseAgg[proci];
+
+    pbrt::Vector3f smallV(SMALL, SMALL, SMALL);
+
+    label i = 0;
+    label j = 0;
+    do
+    {
+        for (; i < myFc.size(); i++)
+        {
+            const point& fc = myFc[i];
+            const vector& fA = myArea[i];
+            const label& fAgg = myAgg[i];
+
+            for (; j < remoteFc.size(); j++)
+            {
+                if (proci != Pstream::myProcNo() || i != j)
+                {
+                    const point& remFc = remoteFc[j];
+                    const vector& remA = remoteArea[j];
+                    const label& remAgg = remoteAgg[j];
+
+                    const vector d = remFc - fc;
+
+                    const vector nd = d/mag(d);
+                    const vector nfA = fA/mag(fA);
+                    const vector nremA = remA/mag(remA);
+
+                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
+                    {
+                        pbrt::Vector3f direction(d[0], d[1], d[2]);
+                        direction += smallV;
+                        pbrt::Point3f s(fc[0], fc[1], fc[2]);
+                        end.push_back(s + direction);
+
+                        s += pbrt::Point3f(0.01*d[0], 0.01*d[1], 0.01*d[2]);
+
+                        start.push_back(s);
+                        startIndex.append(i);
+                        startAgg.append(globalNumbering.toGlobal(proci, fAgg));
+
+                        dir.push_back(direction);
+
+                        label globalI = globalNumbering.toGlobal(proci, j);
+                        endIndex.append(globalI);
+                        endAgg.append(globalNumbering.toGlobal(proci, remAgg));
+                        if (startIndex.size() > maxDynListLength)
+                        {
+                            FatalErrorInFunction
+                                << "Dynamic list need from capacity."
+                                << "Actual size maxDynListLength : "
+                                <<  maxDynListLength
+                                << abort(FatalError);
+                        }
+                    }
+                }
+            }
+
+            if (j == remoteFc.size())
+            {
+                j = 0;
+            }
+        }
+
+    } while (i < myFc.size());
+
+    label totalRays(startIndex.size());
+    reduce(totalRays, sumOp<scalar>());
+
+    if (Pstream::master() && debug)
+    {
+        Info<< "Number of rays :" << totalRays << endl;
+    }
+
+    DynamicList<label> dRayIs;
+
+    std::vector<pbrt::Ray> raysInt;
+    raysInt.reserve(start.size());
+
+    std::vector<int> raysTriIndex;
+    raysTriIndex.reserve(start.size());
+
+    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
+    {
+        pbrt::Ray ray(start[rayI], dir[rayI]);
+
+        pbrt::SurfaceInteraction isect;
+        bool intersects = accel->Intersect(ray, &isect);
+        const Triangle* tri = dynamic_cast<const Triangle *>(isect.shape);
+
+        if (intersects)
+        {
+            const int index = tri->faceIndex;
+
+            const Vector3f delta(ray(ray.tMax) - end[rayI]);
+
+            if (delta.Length() < 1e-4)
+            {
+                rayStartFace.append(startIndex[rayI]);
+                rayEndFace.append(endIndex[rayI]);
+            }
+            else if (triSurfaceToAgglom[index] == startAgg[rayI])
+            {
+                dRayIs.append(rayI);
+                raysInt.push_back(ray);
+                raysTriIndex.push_back(index);
+            }
+        }
+    }
+    start.clear();
+
+    labelField rayIs;
+    rayIs.transfer(dRayIs);
+    dRayIs.clear();
+
+    boolList converged(rayIs.size(), false);
+    label nConverged = 0;
+    label iter = 0;
+    do
+    {
+        forAll(rayIs, rayI)
+        {
+            const label rayID = rayIs[rayI];
+            pbrt::Ray& rayHit = raysInt[rayI];
+            const int index =  raysTriIndex[rayI];
+
+            if (!converged[rayI])
+            {
+                if (triSurfaceToAgglom[index] == startAgg[rayID])
+                {
+                    rayHit.tMax *= 1.1;
+                    rayHit.o = rayHit(rayHit.tMax);
+
+                    pbrt::SurfaceInteraction isect;
+                    bool intersects = accel->Intersect(rayHit, &isect);
+                    if (intersects)
+                    {
+                        const Triangle* tri =
+                            dynamic_cast<const Triangle *>(isect.shape);
+                        const int index = tri->faceIndex;
+
+                        raysTriIndex[rayI] = index;
+                    }
+                }
+                else if (triSurfaceToAgglom[index] == endAgg[rayID])
+                {
+                    converged[rayI] = true;
+                    nConverged++;
+                    rayStartFace.append(startIndex[rayID]);
+                    rayEndFace.append(endIndex[rayID]);
+                }
+            }
+        }
+
+        iter ++;
+
+    } while (nConverged < converged.size() && iter < 10);
+
+    startIndex.clear();
+    end.clear();
+    endIndex.clear();
+    startAgg.clear();
+    endAgg.clear();
+}