applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     viewFactorGenerator
  26
  27 Description
  28     View factors are calculated based on a face agglomeration array
  29     (finalAgglom generated by faceAgglomerate utility).
  30
  31     Each view factor between the agglomerated faces i and j (Fij) is calculated
  32     using a double integral of the sub-areas composing the agglomaration.
  33
  34     The patches involved in the view factor calculation are taken from the Qr
  35     volScalarField (radiative flux) when is greyDiffusiveRadiationViewFactor
  36     otherwise they are not included.
  37
  38 \*---------------------------------------------------------------------------*/
  39
  40
  41 #include "argList.H"
  42 #include "fvMesh.H"
  43 #include "Time.H"
  44 #include "volFields.H"
  45 #include "surfaceFields.H"
  46 #include "distributedTriSurfaceMesh.H"
  47 #include "triSurfaceTools.H"
  48 #include "mapDistribute.H"
  49
  50 #include "OFstream.H"
  51 #include "meshTools.H"
  52 #include "plane.H"
  53 #include "uindirectPrimitivePatch.H"
  54 #include "DynamicField.H"
  55 #include "IFstream.H"
  56 #include "unitConversion.H"
  57
  58 #include "mathematicalConstants.H"
  59 #include "scalarMatrices.H"
  60 #include "CompactListList.H"
  61 #include "labelIOList.H"
  62 #include "labelListIOList.H"
  63 #include "scalarListIOList.H"
  64
  65 #include "singleCellFvMesh.H"
  66 #include "IOdictionary.H"
  67 #include "fixedValueFvPatchFields.H"
  68
  69 using namespace Foam;
  70
  71 void writeRays
  72 (
  73     const fileName& fName,
  74     const pointField& compactCf,
  75     const pointField& myFc,
  76     const labelListList& visibleFaceFaces
  77 )
  78 {
  79     OFstream str(fName);
  80     label vertI = 0;
  81
  82     Pout<< "Dumping rays to " << str.name() << endl;
  83
  84     forAll(myFc, faceI)
  85     {
  86         const labelList visFaces = visibleFaceFaces[faceI];
  87         forAll(visFaces, faceRemote)
  88         {
  89             label compactI = visFaces[faceRemote];
  90             const point& remoteFc = compactCf[compactI];
  91
  92             meshTools::writeOBJ(str, myFc[faceI]);
  93             vertI++;
  94             meshTools::writeOBJ(str, remoteFc);
  95             vertI++;
  96             str << "l " << vertI-1 << ' ' << vertI << nl;
  97         }
  98     }
  99     string cmd("objToVTK " + fName + " " + fName.lessExt() + ".vtk");
 100     Pout<< "cmd:" << cmd << endl;
 101     system(cmd);
 102 }
 103
 104
 105 scalar calculateViewFactorFij
 106 (
 107     const vector& i,
 108     const vector& j,
 109     const vector& dAi,
 110     const vector& dAj
 111 )
 112 {
 113     vector r = i - j;
 114     scalar rMag = mag(r);
 115     scalar dAiMag = mag(dAi);
 116     scalar dAjMag = mag(dAj);
 117
 118     vector ni = dAi/dAiMag;
 119     vector nj = dAj/dAjMag;
 120     scalar cosThetaJ = mag(nj & r)/rMag;
 121     scalar cosThetaI = mag(ni & r)/rMag;
 122
 123     return
 124     (
 125         (cosThetaI*cosThetaJ*dAjMag*dAiMag)
 126        /(sqr(rMag)*constant::mathematical::pi)
 127     );
 128 }
 129
 130
 131 void insertMatrixElements
 132 (
 133     const globalIndex& globalNumbering,
 134     const label fromProcI,
 135     const labelListList& globalFaceFaces,
 136     const scalarListList& viewFactors,
 137     scalarSquareMatrix& matrix
 138 )
 139 {
 140     forAll(viewFactors, faceI)
 141     {
 142         const scalarList& vf = viewFactors[faceI];
 143         const labelList& globalFaces = globalFaceFaces[faceI];
 144
 145         label globalI = globalNumbering.toGlobal(fromProcI, faceI);
 146         forAll(globalFaces, i)
 147         {
 148             matrix[globalI][globalFaces[i]] = vf[i];
 149         }
 150     }
 151 }
 152
 153
 154 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 155
 156 int main(int argc, char *argv[])
 157 {
 158     #include "addRegionOption.H"
 159     #include "setRootCase.H"
 160     #include "createTime.H"
 161     #include "createNamedMesh.H"
 162
 163     // Read view factor dictionary
 164     IOdictionary viewFactorDict
 165     (
 166        IOobject
 167        (
 168             "viewFactorsDict",
 169             runTime.constant(),
 170             mesh,
 171             IOobject::MUST_READ_IF_MODIFIED,
 172             IOobject::NO_WRITE
 173        )
 174     );
 175
 176     const bool writeViewFactors =
 177         viewFactorDict.lookupOrDefault<bool>("writeViewFactorMatrix", false);
 178
 179     const bool dumpRays =
 180         viewFactorDict.lookupOrDefault<bool>("dumpRays", false);
 181
 182     const label debug = viewFactorDict.lookupOrDefault<label>("debug", 0);
 183
 184     volScalarField Qr
 185     (
 186         IOobject
 187         (
 188             "Qr",
 189             runTime.timeName(),
 190             mesh,
 191             IOobject::MUST_READ,
 192             IOobject::NO_WRITE
 193         ),
 194         mesh
 195     );
 196
 197     // Read agglomeration map
 198     labelListIOList finalAgglom
 199     (
 200         IOobject
 201         (
 202             "finalAgglom",
 203             mesh.facesInstance(),
 204             mesh,
 205             IOobject::MUST_READ,
 206             IOobject::NO_WRITE,
 207             false
 208         )
 209     );
 210
 211     // Create the coarse mesh  using agglomeration
 212     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 213
 214     if (debug)
 215     {
 216         Info << "\nCreating single cell mesh..." << endl;
 217     }
 218
 219     singleCellFvMesh coarseMesh
 220     (
 221         IOobject
 222         (
 223             mesh.name(),
 224             runTime.timeName(),
 225             runTime,
 226             IOobject::NO_READ,
 227             IOobject::NO_WRITE
 228         ),
 229         mesh,
 230         finalAgglom
 231     );
 232
 233
 234     // Calculate total number of fine and coarse faces
 235     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 236
 237     label nCoarseFaces = 0;      //total number of coarse faces
 238     label nFineFaces = 0;        //total number of fine faces
 239
 240     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 241     const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
 242
 243     labelList viewFactorsPatches(patches.size());
 244
 245     const volScalarField::GeometricBoundaryField& Qrb = Qr.boundaryField();
 246
 247     label count = 0;
 248     forAll(Qrb, patchI)
 249     {
 250         const fvPatchScalarField& QrpI = Qrb[patchI];
 251
 252         if ((isA<fixedValueFvPatchScalarField>(QrpI)))// && (pp.size() > 0))
 253         {
 254             viewFactorsPatches[count] = QrpI.patch().index();
 255             nCoarseFaces += coarsePatches[patchI].size();
 256             nFineFaces += patches[patchI].size();
 257             count ++;
 258         }
 259     }
 260
 261     viewFactorsPatches.resize(count--);
 262
 263     // total number of coarse faces
 264     label totalNCoarseFaces = nCoarseFaces;
 265
 266     reduce(totalNCoarseFaces, sumOp<label>());
 267
 268     if (Pstream::master())
 269     {
 270         Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
 271     }
 272
 273     if (Pstream::master() && debug)
 274     {
 275         Pout << "\nView factor patches included in the calculation : "
 276              << viewFactorsPatches << endl;
 277     }
 278
 279     // Collect local Cf and Sf on coarse mesh
 280     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 281
 282     DynamicList<point> localCoarseCf(nCoarseFaces);
 283     DynamicList<point> localCoarseSf(nCoarseFaces);
 284
 285     forAll (viewFactorsPatches, i)
 286     {
 287         const label patchID = viewFactorsPatches[i];
 288
 289         const polyPatch& pp = patches[patchID];
 290
 291         if (pp.size() > 0)
 292         {
 293             const labelList& agglom = finalAgglom[patchID];
 294             label nAgglom = max(agglom)+1;
 295             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 296             const labelList& coarsePatchFace =
 297                 coarseMesh.patchFaceMap()[patchID];
 298
 299             const pointField& coarseCf =
 300                 coarseMesh.Cf().boundaryField()[patchID];
 301
 302             const pointField& coarseSf =
 303                 coarseMesh.Sf().boundaryField()[patchID];
 304
 305             forAll(coarseCf, faceI)
 306             {
 307                 point cf = coarseCf[faceI];
 308                 const label coarseFaceI = coarsePatchFace[faceI];
 309                 const labelList& fineFaces = coarseToFine[coarseFaceI];
 310                 // Construct single face
 311                 uindirectPrimitivePatch upp
 312                 (
 313                     UIndirectList<face>(pp, fineFaces),
 314                     pp.points()
 315                 );
 316
 317                 List<point> availablePoints
 318                 (
 319                     upp.faceCentres().size()
 320                   + upp.localPoints().size()
 321                 );
 322
 323                 SubList<point>
 324                 (
 325                     availablePoints,
 326                     upp.faceCentres().size()
 327                 ).assign(upp.faceCentres());
 328
 329                 SubList<point>
 330                 (
 331                     availablePoints,
 332                     upp.localPoints().size(),
 333                     upp.faceCentres().size()
 334                 ).assign(upp.localPoints());
 335
 336                 point cfo = cf;
 337                 scalar dist = GREAT;
 338                 forAll(availablePoints, iPoint)
 339                 {
 340                     point cfFine = availablePoints[iPoint];
 341                     if(mag(cfFine-cfo) < dist)
 342                     {
 343                         dist = mag(cfFine-cfo);
 344                         cf = cfFine;
 345                     }
 346                 }
 347
 348                 point sf = coarseSf[faceI];
 349                 localCoarseCf.append(cf);
 350                 localCoarseSf.append(sf);
 351             }
 352         }
 353     }
 354
 355     // Collect remote Cf and Sf on coarse mesh
 356     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 357
 358     List<pointField> remoteCoarseCf(Pstream::nProcs());
 359     List<pointField> remoteCoarseSf(Pstream::nProcs());
 360
 361     remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
 362     remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
 363
 364     // Collect remote Cf and Sf on fine mesh
 365     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 366
 367     List<pointField> remoteFineCf(Pstream::nProcs());
 368     List<pointField> remoteFineSf(Pstream::nProcs());
 369
 370     remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
 371     remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
 372
 373     // Distribute local coarse Cf and Sf for shooting rays
 374     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 375
 376     Pstream::gatherList(remoteCoarseCf);
 377     Pstream::scatterList(remoteCoarseCf);
 378     Pstream::gatherList(remoteCoarseSf);
 379     Pstream::scatterList(remoteCoarseSf);
 380
 381
 382     // Set up searching engine for obstacles
 383     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 384     #include "searchingEngine.H"
 385
 386
 387     // Determine rays between coarse face centres
 388     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 389     DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
 390
 391     DynamicList<label> rayEndFace(rayStartFace.size());
 392
 393     globalIndex globalNumbering(nCoarseFaces);
 394
 395
 396     // Return rayStartFace in local index andrayEndFace in global index
 397     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 398     #include "shootRays.H"
 399
 400     // Calculate number of visible faces from local index
 401     labelList nVisibleFaceFaces(nCoarseFaces, 0);
 402
 403     forAll(rayStartFace, i)
 404     {
 405         nVisibleFaceFaces[rayStartFace[i]]++;
 406     }
 407
 408     labelListList visibleFaceFaces(nCoarseFaces);
 409
 410     label nViewFactors = 0;
 411     forAll(nVisibleFaceFaces, faceI)
 412     {
 413         visibleFaceFaces[faceI].setSize(nVisibleFaceFaces[faceI]);
 414         nViewFactors += nVisibleFaceFaces[faceI];
 415     }
 416
 417
 418     // - Construct compact numbering
 419     // - return map from remote to compact indices
 420     //   (per processor (!= myProcNo) a map from remote index to compact index)
 421     // - construct distribute map
 422     // - renumber rayEndFace into compact addressing
 423
 424     List<Map<label> > compactMap(Pstream::nProcs());
 425
 426     mapDistribute map(globalNumbering, rayEndFace, compactMap);
 427
 428     labelListIOList IOsubMap
 429     (
 430         IOobject
 431         (
 432             "subMap",
 433             mesh.facesInstance(),
 434             mesh,
 435             IOobject::NO_READ,
 436             IOobject::NO_WRITE,
 437             false
 438         ),
 439         map.subMap()
 440     );
 441     IOsubMap.write();
 442
 443
 444     labelListIOList IOconstructMap
 445     (
 446         IOobject
 447         (
 448             "constructMap",
 449             mesh.facesInstance(),
 450             mesh,
 451             IOobject::NO_READ,
 452             IOobject::NO_WRITE,
 453             false
 454         ),
 455         map.constructMap()
 456     );
 457     IOconstructMap.write();
 458
 459
 460     IOList<label> consMapDim
 461     (
 462         IOobject
 463         (
 464             "constructMapDim",
 465             mesh.facesInstance(),
 466             mesh,
 467             IOobject::NO_READ,
 468             IOobject::NO_WRITE,
 469             false
 470         ),
 471         List<label>(1, map.constructSize())
 472     );
 473     consMapDim.write();
 474
 475     // visibleFaceFaces has:
 476     //    (local face, local viewed face) = compact viewed face
 477     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 478
 479     nVisibleFaceFaces = 0;
 480     forAll(rayStartFace, i)
 481     {
 482         label faceI = rayStartFace[i];
 483         label compactI = rayEndFace[i];
 484         visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
 485     }
 486
 487
 488     // Construct data in compact addressing
 489     // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
 490     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 491
 492     pointField compactCoarseCf(map.constructSize(), pTraits<vector>::zero);
 493     pointField compactCoarseSf(map.constructSize(), pTraits<vector>::zero);
 494     List<List<point> > compactFineSf(map.constructSize());
 495     List<List<point> > compactFineCf(map.constructSize());
 496
 497     DynamicList<label> compactPatchId(map.constructSize());
 498
 499     // Insert my coarse local values
 500     SubList<point>(compactCoarseSf, nCoarseFaces).assign(localCoarseSf);
 501     SubList<point>(compactCoarseCf, nCoarseFaces).assign(localCoarseCf);
 502
 503     // Insert my fine local values
 504     label compactI = 0;
 505     forAll(viewFactorsPatches, i)
 506     {
 507         label patchID = viewFactorsPatches[i];
 508         const polyPatch& pp = patches[patchID];
 509
 510         if (pp.size() > 0)
 511         {
 512             const labelList& agglom = finalAgglom[patchID];
 513             label nAgglom = max(agglom)+1;
 514             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 515             const labelList& coarsePatchFace =
 516                 coarseMesh.patchFaceMap()[patchID];
 517
 518             forAll(coarseToFine, coarseI)
 519             {
 520                 compactPatchId.append(patchID);
 521                 List<point>& fineCf = compactFineCf[compactI];
 522                 List<point>& fineSf = compactFineSf[compactI++];
 523
 524                 const label coarseFaceI = coarsePatchFace[coarseI];
 525                 const labelList& fineFaces = coarseToFine[coarseFaceI];
 526
 527                 fineCf.setSize(fineFaces.size());
 528                 fineSf.setSize(fineFaces.size());
 529
 530                 fineCf = UIndirectList<point>
 531                 (
 532                     mesh.Cf().boundaryField()[patchID],
 533                     coarseToFine[coarseFaceI]
 534                 );
 535                 fineSf = UIndirectList<point>
 536                 (
 537                     mesh.Sf().boundaryField()[patchID],
 538                     coarseToFine[coarseFaceI]
 539                 );
 540             }
 541         }
 542     }
 543
 544     // Do all swapping
 545     map.distribute(compactCoarseSf);
 546     map.distribute(compactCoarseCf);
 547     map.distribute(compactFineCf);
 548     map.distribute(compactFineSf);
 549
 550     map.distribute(compactPatchId);
 551
 552
 553     // Plot all rays between visible faces.
 554     if (dumpRays)
 555     {
 556         writeRays
 557         (
 558             runTime.path()/"allVisibleFaces.obj",
 559             compactCoarseCf,
 560             remoteCoarseCf[Pstream::myProcNo()],
 561             visibleFaceFaces
 562         );
 563     }
 564
 565
 566     // Fill local view factor matrix
 567     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 568
 569     scalarListIOList F
 570     (
 571         IOobject
 572         (
 573             "F",
 574             mesh.facesInstance(),
 575             mesh,
 576             IOobject::NO_READ,
 577             IOobject::NO_WRITE,
 578             false
 579         ),
 580         nCoarseFaces
 581     );
 582
 583     label totalPatches = coarsePatches.size();
 584     reduce(totalPatches, maxOp<label>());
 585
 586     // Matrix sum in j(Fij) for each i (if enclosure sum = 1
 587     scalarSquareMatrix sumViewFactorPatch
 588     (
 589         totalPatches,
 590         totalPatches,
 591         0.0
 592     );
 593
 594     scalarList patchArea(totalPatches, 0.0);
 595
 596     if (Pstream::master())
 597     {
 598         Info<< "\nCalculating view factors..." << endl;
 599     }
 600
 601     if (mesh.nSolutionD() == 3)
 602     {
 603         forAll (localCoarseSf, coarseFaceI)
 604         {
 605             const List<point>& localFineSf = compactFineSf[coarseFaceI];
 606             const vector Ai = sum(localFineSf);
 607             const List<point>& localFineCf = compactFineCf[coarseFaceI];
 608             const label fromPatchId = compactPatchId[coarseFaceI];
 609             patchArea[fromPatchId] += mag(Ai);
 610
 611             const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 612
 613             forAll(visCoarseFaces, visCoarseFaceI)
 614             {
 615                 F[coarseFaceI].setSize(visCoarseFaces.size());
 616                 label compactJ = visCoarseFaces[visCoarseFaceI];
 617                 const List<point>& remoteFineSj = compactFineSf[compactJ];
 618                 const List<point>& remoteFineCj = compactFineCf[compactJ];
 619
 620                 const label toPatchId = compactPatchId[compactJ];
 621
 622                 scalar Fij = 0;
 623                 forAll (localFineSf, i)
 624                 {
 625                     const vector& dAi = localFineSf[i];
 626                     const vector& dCi = localFineCf[i];
 627
 628                     forAll (remoteFineSj, j)
 629                     {
 630                         const vector& dAj = remoteFineSj[j];
 631                         const vector& dCj = remoteFineCj[j];
 632
 633                         scalar dIntFij = calculateViewFactorFij
 634                         (
 635                             dCi,
 636                             dCj,
 637                             dAi,
 638                             dAj
 639                         );
 640
 641                         Fij += dIntFij;
 642                     }
 643                 }
 644                 F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
 645                 sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
 646             }
 647         }
 648     }
 649     else if (mesh.nSolutionD() == 2)
 650     {
 651         const boundBox& box = mesh.bounds();
 652         const Vector<label>& dirs = mesh.geometricD();
 653         vector emptyDir = vector::zero;
 654         forAll(dirs, i)
 655         {
 656             if (dirs[i] == -1)
 657             {
 658                 emptyDir[i] = 1.0;
 659             }
 660         }
 661
 662         scalar wideBy2 = (box.span() & emptyDir)*2.0;
 663
 664         forAll(localCoarseSf, coarseFaceI)
 665         {
 666             const vector& Ai = localCoarseSf[coarseFaceI];
 667             const vector& Ci = localCoarseCf[coarseFaceI];
 668             vector Ain = Ai/mag(Ai);
 669             vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
 670             vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
 671
 672             const label fromPatchId = compactPatchId[coarseFaceI];
 673             patchArea[fromPatchId] += mag(Ai);
 674
 675             const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 676             forAll (visCoarseFaces, visCoarseFaceI)
 677             {
 678                 F[coarseFaceI].setSize(visCoarseFaces.size());
 679                 label compactJ = visCoarseFaces[visCoarseFaceI];
 680                 const vector& Aj = compactCoarseSf[compactJ];
 681                 const vector& Cj = compactCoarseCf[compactJ];
 682
 683                 const label toPatchId = compactPatchId[compactJ];
 684
 685                 vector Ajn = Aj/mag(Aj);
 686                 vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
 687                 vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
 688
 689                 scalar d1 = mag(R1i - R2j);
 690                 scalar d2 = mag(R2i - R1j);
 691                 scalar s1 = mag(R1i - R1j);
 692                 scalar s2 = mag(R2i - R2j);
 693
 694                 scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
 695
 696                 F[coarseFaceI][visCoarseFaceI] = Fij;
 697                 sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
 698             }
 699         }
 700     }
 701
 702     if (Pstream::master())
 703     {
 704         Info << "Writing view factor matrix..." << endl;
 705     }
 706
 707     // Write view factors matrix in listlist form
 708     F.write();
 709
 710     reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
 711     reduce(patchArea, sumOp<scalarList>());
 712
 713
 714     if (Pstream::master() && debug)
 715     {
 716         forAll(viewFactorsPatches, i)
 717         {
 718             label patchI =  viewFactorsPatches[i];
 719             forAll(viewFactorsPatches, i)
 720             {
 721                 label patchJ =  viewFactorsPatches[i];
 722                 Info << "F" << patchI << patchJ << ": "
 723                      << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
 724                      << endl;
 725             }
 726         }
 727     }
 728
 729
 730     if (writeViewFactors)
 731     {
 732         volScalarField viewFactorField
 733         (
 734             IOobject
 735             (
 736                 "viewFactorField",
 737                 mesh.time().timeName(),
 738                 mesh,
 739                 IOobject::NO_READ,
 740                 IOobject::NO_WRITE
 741             ),
 742             mesh,
 743             dimensionedScalar("viewFactorField", dimless, 0)
 744         );
 745
 746         label compactI = 0;
 747         forAll(viewFactorsPatches, i)
 748         {
 749             label patchID = viewFactorsPatches[i];
 750             const polyPatch& pp = patches[patchID];
 751
 752             if (pp.size() > 0)
 753             {
 754                 const labelList& agglom = finalAgglom[patchID];
 755                 label nAgglom = max(agglom)+1;
 756                 labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 757                 const labelList& coarsePatchFace =
 758                     coarseMesh.patchFaceMap()[patchID];
 759
 760                 forAll(coarseToFine, coarseI)
 761                 {
 762                     const scalar Fij = sum(F[compactI]);
 763                     const label coarseFaceID = coarsePatchFace[coarseI];
 764                     const labelList& fineFaces = coarseToFine[coarseFaceID];
 765                     forAll (fineFaces, fineId)
 766                     {
 767                         const label faceID = fineFaces[fineId];
 768                         viewFactorField.boundaryField()[patchID][faceID] = Fij;
 769                     }
 770                     compactI++;
 771                 }
 772             }
 773         }
 774         viewFactorField.write();
 775     }
 776
 777
 778     // Invert compactMap (from processor+localface to compact) to go
 779     // from compact to processor+localface (expressed as a globalIndex)
 780     // globalIndex globalCoarFaceNum(coarseMesh.nFaces());
 781     labelList compactToGlobal(map.constructSize());
 782
 783     // Local indices first (note: are not in compactMap)
 784     for (label i = 0; i < globalNumbering.localSize(); i++)
 785     {
 786         compactToGlobal[i] = globalNumbering.toGlobal(i);
 787     }
 788
 789
 790     forAll(compactMap, procI)
 791     {
 792         const Map<label>& localToCompactMap = compactMap[procI];
 793
 794         forAllConstIter(Map<label>, localToCompactMap, iter)
 795         {
 796             compactToGlobal[iter()] = globalNumbering.toGlobal
 797             (
 798                 procI,
 799                 iter.key()
 800             );
 801         }
 802     }
 803
 804
 805     labelListList globalFaceFaces(visibleFaceFaces.size());
 806
 807     // Create globalFaceFaces needed to insert view factors
 808     // in F to the global matrix Fmatrix
 809     forAll(globalFaceFaces, faceI)
 810     {
 811         globalFaceFaces[faceI] = renumber
 812         (
 813             compactToGlobal,
 814             visibleFaceFaces[faceI]
 815         );
 816     }
 817
 818     labelListIOList IOglobalFaceFaces
 819     (
 820         IOobject
 821         (
 822             "globalFaceFaces",
 823             mesh.facesInstance(),
 824             mesh,
 825             IOobject::NO_READ,
 826             IOobject::NO_WRITE,
 827             false
 828         ),
 829         globalFaceFaces
 830     );
 831     IOglobalFaceFaces.write();
 832
 833     Info<< "End\n" << endl;
 834     return 0;
 835 }
 836
 837 // ************************************************************************* //