ENH: viewFactorsGen:

Modification of shootRays algorithm. Rays that hit themselfs are continued
 	 until they hit(or not) the agglomeration target. If they do the surfaces can see
	 each other.
     faceAgglomerate:
	 The agglomeration now stops when the global nCoarse is reached (not per-processor)

applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,7 @@ Description
 #include "pairPatchAgglomeration.H"
 #include "labelListIOList.H"
 #include "syncTools.H"
+#include "globalIndex.H"
 
 using namespace Foam;
 
@@ -53,7 +54,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createNamedMesh.H"
 
-    const word dictName("faceAgglomerateDict");
+    const word dictName("viewFactorsDict");
 
     #include "setConstantMeshDictionaryIO.H"
 
@@ -79,16 +80,16 @@ int main(int argc, char *argv[])
         boundary.size()
     );
 
-    forAll(boundary, patchId)
+    label nCoarseFaces = 0;
+
+    forAllConstIter(dictionary, agglomDict, iter)
     {
-        const polyPatch& pp = boundary[patchId];
+        labelList patchIds = boundary.findIndices(iter().keyword());
-
+        forAll(patchIds, i)
-        label patchI = pp.index();
+        {
-        finalAgglom[patchI].setSize(pp.size(), 0);
+            label patchI =  patchIds[i];
-
+            const polyPatch& pp = boundary[patchI];
             if (!pp.coupled())
-        {
-            if (agglomDict.found(pp.name()))
             {
                 Info << "\nAgglomerating patch : " << pp.name() << endl;
                 pairPatchAgglomeration agglomObject
@@ -99,12 +100,11 @@ int main(int argc, char *argv[])
                 agglomObject.agglomerate();
                 finalAgglom[patchI] =
                     agglomObject.restrictTopBottomAddressing();
-            }
+
-            else
+                if (finalAgglom[patchI].size())
                 {
-                FatalErrorIn(args.executable())
+                    nCoarseFaces += max(finalAgglom[patchI] + 1);
-                    << "Patch " << pp.name() << " not found in dictionary: "
+                }
-                    << agglomDict.name() << exit(FatalError);
             }
         }
     }
@@ -144,6 +144,7 @@ int main(int argc, char *argv[])
 
     if (writeAgglom)
     {
+        globalIndex index(nCoarseFaces);
         volScalarField facesAgglomeration
         (
             IOobject
@@ -158,14 +159,26 @@ int main(int argc, char *argv[])
             dimensionedScalar("facesAgglomeration", dimless, 0)
         );
 
+        label coarsePatchIndex = 0;
         forAll(boundary, patchId)
+        {
+            const polyPatch& pp = boundary[patchId];
+            if (pp.size() > 0)
             {
                 fvPatchScalarField& bFacesAgglomeration =
                     facesAgglomeration.boundaryField()[patchId];
 
                 forAll(bFacesAgglomeration, j)
                 {
-                bFacesAgglomeration[j] = finalAgglom[patchId][j];
+                    bFacesAgglomeration[j] =
+                        index.toGlobal
+                        (
+                            Pstream::myProcNo(),
+                            finalAgglom[patchId][j] + coarsePatchIndex
+                        );
+                }
+
+                coarsePatchIndex += max(finalAgglom[patchId]) + 1;
             }
         }
 

applications/utilities/preProcessing/faceAgglomerate/viewFactorsDict

@@ -17,6 +17,12 @@ FoamFile
 // Write agglomeration as a volScalarField with calculated boundary values
 writeFacesAgglomeration   true;
 
+//Debug option
+debug                     0;
+
+//Dump connectivity rays
+dumpRays                  false;
+
 // Per patch (wildcard possible) the coarsening level
 bottomAir_to_heater
 {

applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H

@@ -32,6 +32,19 @@ forAll(patches, patchI)
     }
 }
 
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+const triSurface localSurface = triangulate
+(
+    patches,
+    includePatches,
+    finalAgglom,
+    triSurfaceToAgglom,
+    globalNumbering,
+    coarsePatches
+);
+
+
 distributedTriSurfaceMesh surfacesMesh
 (
     IOobject
@@ -43,12 +56,13 @@ distributedTriSurfaceMesh surfacesMesh
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
-    triSurfaceTools::triangulate
+    localSurface,
-    (
-        patches,
-        includePatches
-    ),
     dict
 );
 
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
 //surfacesMesh.searchableSurface::write();
+
+surfacesMesh.setField(triSurfaceToAgglom);

applications/utilities/preProcessing/viewFactorsGen/shootRays.H

@@ -2,7 +2,7 @@
 // Pre-size by assuming a certain percentage is visible.
 
 // Maximum lenght for dynamicList
-const label maxDynListLength = 10000;
+const label maxDynListLength = 100000;
 
 for (label procI = 0; procI < Pstream::nProcs(); procI++)
 {
@@ -14,12 +14,17 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
     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;
@@ -29,6 +34,7 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
         {
             const point& fc = myFc[i];
             const vector& fA = myArea[i];
+            const label& fAgg = myAgg[i];
 
             for (; j < remoteFc.size(); j++)//
             {
@@ -36,26 +42,32 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
                 {
                     const point& remFc = remoteFc[j];
                     const vector& remA = remoteArea[j];
+                    const label& remAgg = remoteAgg[j];
+
                     const vector& d = remFc - fc;
 
                     if (((d & fA) < 0.) && ((d & remA) > 0))
                     {
-                        start.append(fc + 0.0001*d);
+                        start.append(fc + 0.001*d);
                         startIndex.append(i);
-                        end.append(fc + 0.9999*d);
+                        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)
                         {
-                            break;
+                            FatalErrorIn
+                            (
+                                "shootRays"
+                            )   << "Dynamic list need from capacity."
+                                << "Actual size maxDynListLength : "
+                                <<  maxDynListLength
+                                << abort(FatalError);
                         }
                     }
                 }
             }
-            if (startIndex.size() > maxDynListLength)
-            {
-                break;
-            }
 
             if (j == remoteFc.size())
             {
@@ -63,9 +75,19 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
             }
         }
 
+    }while (returnReduce(i < myFc.size(), orOp<bool>()));
+
     List<pointIndexHit> hitInfo(startIndex.size());
     surfacesMesh.findLine(start, end, hitInfo);
 
+    // Return hit global agglo index
+    labelList aggHitIndex;
+    surfacesMesh.getField(hitInfo, aggHitIndex);
+
+    DynamicList<label> dRayIs;
+
+    // Collect the rays which has not abstacle in bettween in rayStartFace
+    // and rayEndFace. If the ray hit itself get stored in dRayIs
     forAll (hitInfo, rayI)
     {
         if (!hitInfo[rayI].hit())
@@ -73,13 +95,82 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
             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
+    {
+        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();
     endIndex.clear();
-
+    startAgg.clear();
-    }while (returnReduce(i < myFc.size(), orOp<bool>()));
+    endAgg.clear();
-
 }

applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C

@@ -68,6 +68,101 @@ Description
 
 using namespace Foam;
 
+
+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.nFaces() - mesh.nInternalFaces()
+    );
+
+    label newPatchI = 0;
+    label localTriFaceI = 0;
+
+    forAllConstIter(labelHashSet, includePatches, iter)
+    {
+        const label patchI = iter.key();
+        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();
+
+    // Create globally numbered tri surface
+    triSurface rawSurface(triangles, mesh.points());
+
+    // Create locally numbered tri surface
+    triSurface surface
+    (
+        rawSurface.localFaces(),
+        rawSurface.localPoints()
+    );
+
+    // Add patch names to surface
+    surface.patches().setSize(newPatchI);
+
+    newPatchI = 0;
+
+    forAllConstIter(labelHashSet, includePatches, iter)
+    {
+        const label patchI = iter.key();
+        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,
@@ -213,7 +308,7 @@ int main(int argc, char *argv[])
 
     if (debug)
     {
-        Info << "\nCreating single cell mesh..." << endl;
+        Pout << "\nCreating single cell mesh..." << endl;
     }
 
     singleCellFvMesh coarseMesh
@@ -230,6 +325,11 @@ int main(int argc, char *argv[])
         finalAgglom
     );
 
+    if (debug)
+    {
+        Pout << "\nCreated single cell mesh..." << endl;
+    }
+
 
     // Calculate total number of fine and coarse faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -283,6 +383,8 @@ int main(int argc, char *argv[])
     DynamicList<point> localCoarseCf(nCoarseFaces);
     DynamicList<point> localCoarseSf(nCoarseFaces);
 
+    DynamicList<label> localAgg(nCoarseFaces);
+
     forAll (viewFactorsPatches, i)
     {
         const label patchID = viewFactorsPatches[i];
@@ -296,11 +398,18 @@ int main(int argc, char *argv[])
         const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
         const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
 
+        labelHashSet includePatches;
+        includePatches.insert(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
             (
@@ -308,6 +417,7 @@ int main(int argc, char *argv[])
                 pp.points()
             );
 
+
             List<point> availablePoints
             (
                 upp.faceCentres().size()
@@ -342,6 +452,7 @@ int main(int argc, char *argv[])
             point sf = coarseSf[faceI];
             localCoarseCf.append(cf);
             localCoarseSf.append(sf);
+            localAgg.append(agglomI);
         }
     }
 
@@ -350,9 +461,12 @@ int main(int argc, char *argv[])
 
     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;
+
 
     // Collect remote Cf and Sf on fine mesh
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -370,8 +484,12 @@ int main(int argc, char *argv[])
     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 "searchingEngine.H"
@@ -383,8 +501,6 @@ int main(int argc, char *argv[])
 
     DynamicList<label> rayEndFace(rayStartFace.size());
 
-    globalIndex globalNumbering(nCoarseFaces);
-
 
     // Return rayStartFace in local index andrayEndFace in global index
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,8 +43,10 @@ bool Foam::pairPatchAgglomeration::continueAgglomerating
 )
 {
     // Check the need for further agglomeration on all processors
-    bool contAgg = nCoarseFaces >= nFacesInCoarsestLevel_;
+    label localnCoarseFaces = nCoarseFaces;
-    reduce(contAgg, andOp<bool>());
+//     reduce(localnCoarseFaces, sumOp<label>());
+    bool contAgg = localnCoarseFaces >= nFacesInCoarsestLevel_;
+    //reduce(contAgg, andOp<bool>());
     return contAgg;
 }
 
@@ -263,6 +265,11 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
         )   << "min(fineToCoarse) == -1" << exit(FatalError);
     }
 
+    if (fineToCoarse.size() == 0)
+    {
+        return true;
+    }
+
     if (fineToCoarse.size() != patch.size())
     {
         FatalErrorIn
@@ -282,6 +289,7 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
     const label nCoarseI =  max(fineToCoarse)+1;
     List<face> patchFaces(nCoarseI);
 
+
     // Patch faces per agglomeration
     labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));
 
@@ -335,6 +343,7 @@ bool Foam::pairPatchAgglomeration::agglomeratePatch
             patch.points()
         )
     );
+
     return true;
 }
 
@@ -343,40 +352,34 @@ void Foam::pairPatchAgglomeration:: agglomerate()
 {
     label nPairLevels = 0;
     label nCreatedLevels = 1; //0 level is the base patch
+    label nCoarseFaces = 0;
+    label nCoarseFacesOld = 0;
 
     while (nCreatedLevels < maxLevels_)
     {
-        label nCoarseCells = -1;
-
         const bPatch& patch = patchLevels_[nCreatedLevels - 1];
         tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
 
         bool agglomOK = false;
-        while (!agglomOK && patch.size())
+        while (!agglomOK)
         {
             finalAgglomPtr = agglomerateOneLevel
             (
-                nCoarseCells,
+                nCoarseFaces,
                 patch
             );
 
+            if (nCoarseFaces > 0)
+            {
                 agglomOK = agglomeratePatch
                 (
                     patch,
                     finalAgglomPtr,
                     nCreatedLevels
                 );
-        }
 
-        nFaces_[nCreatedLevels] = nCoarseCells;
                 restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
                 mapBaseToTopAgglom(nCreatedLevels);
-
-        if (!continueAgglomerating(nCoarseCells))
-        {
-            break;
-        }
-
                 setEdgeWeights(nCreatedLevels);
 
                 if (nPairLevels % mergeLevels_)
@@ -390,12 +393,32 @@ void Foam::pairPatchAgglomeration:: agglomerate()
 
                 nPairLevels++;
             }
+            else
+            {
+                agglomOK = true;
+            }
+            reduce(nCoarseFaces, sumOp<label>());
+        }
+
+        nFaces_[nCreatedLevels] = nCoarseFaces;
+
+        if
+        (
+            !continueAgglomerating(nCoarseFaces)
+          || (nCoarseFacesOld ==  nCoarseFaces)
+        )
+        {
+            break;
+        }
+
+        nCoarseFacesOld = nCoarseFaces;
+    }
 }
 
 
 Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
 (
-    label& nCoarseCells,
+    label& nCoarseFaces,
     const bPatch& patch
 )
 {
@@ -406,7 +429,7 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
 
     const labelListList& faceFaces = patch.faceFaces();
 
-    nCoarseCells = 0;
+    nCoarseFaces = 0;
 
     forAll(faceFaces, facei)
     {
@@ -440,9 +463,9 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
             if (matchFaceNo >= 0)
             {
                 // Make a new group
-                coarseCellMap[matchFaceNo] = nCoarseCells;
+                coarseCellMap[matchFaceNo] = nCoarseFaces;
-                coarseCellMap[matchFaceNeibNo] = nCoarseCells;
+                coarseCellMap[matchFaceNeibNo] = nCoarseFaces;
-                nCoarseCells++;
+                nCoarseFaces++;
             }
             else
             {
@@ -475,8 +498,8 @@ Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
                 else
                 {
                     // if not create single-cell "clusters" for each
-                    coarseCellMap[facei] = nCoarseCells;
+                    coarseCellMap[facei] = nCoarseFaces;
-                    nCoarseCells ++;
+                    nCoarseFaces ++;
                 }
             }
         }