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-2011 OpenCFD Ltd.
   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 polyPatch& pp = patches[patchI];
 251         const fvPatchScalarField& QrpI = Qrb[patchI];
 252
 253         if ((isA<fixedValueFvPatchScalarField>(QrpI)) && (pp.size() > 0))
 254         {
 255             viewFactorsPatches[count] = QrpI.patch().index();
 256             nCoarseFaces += coarsePatches[patchI].size();
 257             nFineFaces += patches[patchI].size();
 258             count ++;
 259         }
 260     }
 261
 262     viewFactorsPatches.resize(count--);
 263
 264     // total number of coarse faces
 265     label totalNCoarseFaces = nCoarseFaces;
 266
 267     reduce(totalNCoarseFaces, sumOp<label>());
 268
 269     if (Pstream::master())
 270     {
 271         Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
 272     }
 273
 274     if (Pstream::master() && debug)
 275     {
 276         Pout << "\nView factor patches included in the calculation : "
 277              << viewFactorsPatches << endl;
 278     }
 279
 280     // Collect local Cf and Sf on coarse mesh
 281     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 282
 283     DynamicList<point> localCoarseCf(nCoarseFaces);
 284     DynamicList<point> localCoarseSf(nCoarseFaces);
 285
 286     forAll (viewFactorsPatches, i)
 287     {
 288         const label patchID = viewFactorsPatches[i];
 289
 290         const polyPatch& pp = patches[patchID];
 291         const labelList& agglom = finalAgglom[patchID];
 292         label nAgglom = max(agglom)+1;
 293         labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 294         const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchID];
 295
 296         const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
 297         const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
 298
 299         forAll(coarseCf, faceI)
 300         {
 301             point cf = coarseCf[faceI];
 302             const label coarseFaceI = coarsePatchFace[faceI];
 303             const labelList& fineFaces = coarseToFine[coarseFaceI];
 304             // Construct single face
 305             uindirectPrimitivePatch upp
 306             (
 307                 UIndirectList<face>(pp, fineFaces),
 308                 pp.points()
 309             );
 310
 311             List<point> availablePoints
 312             (
 313                 upp.faceCentres().size()
 314               + upp.localPoints().size()
 315             );
 316
 317             SubList<point>
 318             (
 319                 availablePoints,
 320                 upp.faceCentres().size()
 321             ).assign(upp.faceCentres());
 322
 323             SubList<point>
 324             (
 325                 availablePoints,
 326                 upp.localPoints().size(),
 327                 upp.faceCentres().size()
 328             ).assign(upp.localPoints());
 329
 330             point cfo = cf;
 331             scalar dist = GREAT;
 332             forAll(availablePoints, iPoint)
 333             {
 334                 point cfFine = availablePoints[iPoint];
 335                 if(mag(cfFine-cfo) < dist)
 336                 {
 337                     dist = mag(cfFine-cfo);
 338                     cf = cfFine;
 339                 }
 340             }
 341
 342             point sf = coarseSf[faceI];
 343             localCoarseCf.append(cf);
 344             localCoarseSf.append(sf);
 345         }
 346     }
 347
 348     // Collect remote Cf and Sf on coarse mesh
 349     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 350
 351     List<pointField> remoteCoarseCf(Pstream::nProcs());
 352     List<pointField> remoteCoarseSf(Pstream::nProcs());
 353
 354     remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
 355     remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
 356
 357     // Collect remote Cf and Sf on fine mesh
 358     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 359
 360     List<pointField> remoteFineCf(Pstream::nProcs());
 361     List<pointField> remoteFineSf(Pstream::nProcs());
 362
 363     remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
 364     remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
 365
 366     // Distribute local coarse Cf and Sf for shooting rays
 367     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 368
 369     Pstream::gatherList(remoteCoarseCf);
 370     Pstream::scatterList(remoteCoarseCf);
 371     Pstream::gatherList(remoteCoarseSf);
 372     Pstream::scatterList(remoteCoarseSf);
 373
 374
 375     // Set up searching engine for obstacles
 376     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 377     #include "searchingEngine.H"
 378
 379
 380     // Determine rays between coarse face centres
 381     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 382     DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
 383
 384     DynamicList<label> rayEndFace(rayStartFace.size());
 385
 386     globalIndex globalNumbering(nCoarseFaces);
 387
 388
 389     // Return rayStartFace in local index andrayEndFace in global index
 390     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 391
 392     #include "shootRays.H"
 393
 394     // Calculate number of visible faces from local index
 395     labelList nVisibleFaceFaces(nCoarseFaces, 0);
 396
 397     forAll(rayStartFace, i)
 398     {
 399         nVisibleFaceFaces[rayStartFace[i]]++;
 400     }
 401
 402     labelListList visibleFaceFaces(nCoarseFaces);
 403
 404     label nViewFactors = 0;
 405     forAll(nVisibleFaceFaces, faceI)
 406     {
 407         visibleFaceFaces[faceI].setSize(nVisibleFaceFaces[faceI]);
 408         nViewFactors += nVisibleFaceFaces[faceI];
 409     }
 410
 411
 412     // - Construct compact numbering
 413     // - return map from remote to compact indices
 414     //   (per processor (!= myProcNo) a map from remote index to compact index)
 415     // - construct distribute map
 416     // - renumber rayEndFace into compact addressing
 417
 418     List<Map<label> > compactMap(Pstream::nProcs());
 419
 420     mapDistribute map(globalNumbering, rayEndFace, compactMap);
 421
 422     labelListIOList IOsubMap
 423     (
 424         IOobject
 425         (
 426             "subMap",
 427             mesh.facesInstance(),
 428             mesh,
 429             IOobject::NO_READ,
 430             IOobject::NO_WRITE,
 431             false
 432         ),
 433         map.subMap()
 434     );
 435     IOsubMap.write();
 436
 437
 438     labelListIOList IOconstructMap
 439     (
 440         IOobject
 441         (
 442             "constructMap",
 443             mesh.facesInstance(),
 444             mesh,
 445             IOobject::NO_READ,
 446             IOobject::NO_WRITE,
 447             false
 448         ),
 449         map.constructMap()
 450     );
 451     IOconstructMap.write();
 452
 453
 454     IOList<label> consMapDim
 455     (
 456         IOobject
 457         (
 458             "constructMapDim",
 459             mesh.facesInstance(),
 460             mesh,
 461             IOobject::NO_READ,
 462             IOobject::NO_WRITE,
 463             false
 464         ),
 465         List<label>(1, map.constructSize())
 466     );
 467     consMapDim.write();
 468
 469
 470     // visibleFaceFaces has:
 471     //    (local face, local viewed face) = compact viewed face
 472     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 473
 474     nVisibleFaceFaces = 0;
 475     forAll(rayStartFace, i)
 476     {
 477         label faceI = rayStartFace[i];
 478         label compactI = rayEndFace[i];
 479         visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
 480     }
 481
 482
 483     // Construct data in compact addressing
 484     // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
 485     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 486
 487     pointField compactCoarseCf(map.constructSize(), pTraits<vector>::zero);
 488     pointField compactCoarseSf(map.constructSize(), pTraits<vector>::zero);
 489     List<List<point> > compactFineSf(map.constructSize());
 490     List<List<point> > compactFineCf(map.constructSize());
 491
 492     DynamicList<label> compactPatchId(map.constructSize());
 493
 494     // Insert my coarse local values
 495     SubList<point>(compactCoarseSf, nCoarseFaces).assign(localCoarseSf);
 496     SubList<point>(compactCoarseCf, nCoarseFaces).assign(localCoarseCf);
 497
 498     // Insert my fine local values
 499     label compactI = 0;
 500     forAll(viewFactorsPatches, i)
 501     {
 502         label patchID = viewFactorsPatches[i];
 503         const labelList& agglom = finalAgglom[patchID];
 504         label nAgglom = max(agglom)+1;
 505         labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 506         const labelList& coarsePatchFace = coarseMesh.patchFaceMap()[patchID];
 507
 508         forAll(coarseToFine, coarseI)
 509         {
 510             compactPatchId.append(patchID);
 511             List<point>& fineCf = compactFineCf[compactI];
 512             List<point>& fineSf = compactFineSf[compactI++];
 513
 514             const label coarseFaceI = coarsePatchFace[coarseI];
 515             const labelList& fineFaces = coarseToFine[coarseFaceI];
 516
 517             fineCf.setSize(fineFaces.size());
 518             fineSf.setSize(fineFaces.size());
 519
 520             fineCf = UIndirectList<point>
 521             (
 522                 mesh.Cf().boundaryField()[patchID],
 523                 coarseToFine[coarseFaceI]
 524             );
 525             fineSf = UIndirectList<point>
 526             (
 527                 mesh.Sf().boundaryField()[patchID],
 528                 coarseToFine[coarseFaceI]
 529             );
 530         }
 531     }
 532
 533     // Do all swapping
 534     map.distribute(compactCoarseSf);
 535     map.distribute(compactCoarseCf);
 536     map.distribute(compactFineCf);
 537     map.distribute(compactFineSf);
 538
 539     map.distribute(compactPatchId);
 540
 541
 542     // Plot all rays between visible faces.
 543     if (dumpRays)
 544     {
 545         writeRays
 546         (
 547             runTime.path()/"allVisibleFaces.obj",
 548             compactCoarseCf,
 549             remoteCoarseCf[Pstream::myProcNo()],
 550             visibleFaceFaces
 551         );
 552     }
 553
 554
 555     // Fill local view factor matrix
 556     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 557
 558     scalarListIOList F
 559     (
 560         IOobject
 561         (
 562             "F",
 563             mesh.facesInstance(),
 564             mesh,
 565             IOobject::NO_READ,
 566             IOobject::NO_WRITE,
 567             false
 568         ),
 569         nCoarseFaces
 570     );
 571
 572     label totalPatches = coarsePatches.size();
 573     reduce(totalPatches, maxOp<label>());
 574
 575     // Matrix sum in j(Fij) for each i (if enclosure sum = 1
 576     scalarSquareMatrix sumViewFactorPatch
 577     (
 578         totalPatches,
 579         totalPatches,
 580         0.0
 581     );
 582
 583     scalarList patchArea(totalPatches, 0.0);
 584
 585     if (Pstream::master())
 586     {
 587         Info<< "\nCalculating view factors..." << endl;
 588     }
 589
 590     if (mesh.nSolutionD() == 3)
 591     {
 592         forAll (localCoarseSf, coarseFaceI)
 593         {
 594             const List<point>& localFineSf = compactFineSf[coarseFaceI];
 595             const vector Ai = sum(localFineSf);
 596             const List<point>& localFineCf = compactFineCf[coarseFaceI];
 597             const label fromPatchId = compactPatchId[coarseFaceI];
 598             patchArea[fromPatchId] += mag(Ai);
 599
 600             const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 601
 602             forAll(visCoarseFaces, visCoarseFaceI)
 603             {
 604                 F[coarseFaceI].setSize(visCoarseFaces.size());
 605                 label compactJ = visCoarseFaces[visCoarseFaceI];
 606                 const List<point>& remoteFineSj = compactFineSf[compactJ];
 607                 const List<point>& remoteFineCj = compactFineCf[compactJ];
 608
 609                 const label toPatchId = compactPatchId[compactJ];
 610
 611                 scalar Fij = 0;
 612                 forAll (localFineSf, i)
 613                 {
 614                     const vector& dAi = localFineSf[i];
 615                     const vector& dCi = localFineCf[i];
 616
 617                     forAll (remoteFineSj, j)
 618                     {
 619                         const vector& dAj = remoteFineSj[j];
 620                         const vector& dCj = remoteFineCj[j];
 621
 622                         scalar dIntFij = calculateViewFactorFij
 623                         (
 624                             dCi,
 625                             dCj,
 626                             dAi,
 627                             dAj
 628                         );
 629
 630                         Fij += dIntFij;
 631                     }
 632                 }
 633                 F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
 634                 sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
 635             }
 636         }
 637     }
 638     else if (mesh.nSolutionD() == 2)
 639     {
 640         const boundBox& box = mesh.bounds();
 641         const Vector<label>& dirs = mesh.geometricD();
 642         vector emptyDir = vector::zero;
 643         forAll(dirs, i)
 644         {
 645             if (dirs[i] == -1)
 646             {
 647                 emptyDir[i] = 1.0;
 648             }
 649         }
 650
 651         scalar wideBy2 = (box.span() & emptyDir)*2.0;
 652
 653         forAll(localCoarseSf, coarseFaceI)
 654         {
 655             const vector& Ai = localCoarseSf[coarseFaceI];
 656             const vector& Ci = localCoarseCf[coarseFaceI];
 657             vector Ain = Ai/mag(Ai);
 658             vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
 659             vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
 660
 661             const label fromPatchId = compactPatchId[coarseFaceI];
 662             patchArea[fromPatchId] += mag(Ai);
 663
 664             const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 665             forAll (visCoarseFaces, visCoarseFaceI)
 666             {
 667                 F[coarseFaceI].setSize(visCoarseFaces.size());
 668                 label compactJ = visCoarseFaces[visCoarseFaceI];
 669                 const vector& Aj = compactCoarseSf[compactJ];
 670                 const vector& Cj = compactCoarseCf[compactJ];
 671
 672                 const label toPatchId = compactPatchId[compactJ];
 673
 674                 vector Ajn = Aj/mag(Aj);
 675                 vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
 676                 vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
 677
 678                 scalar d1 = mag(R1i - R2j);
 679                 scalar d2 = mag(R2i - R1j);
 680                 scalar s1 = mag(R1i - R1j);
 681                 scalar s2 = mag(R2i - R2j);
 682
 683                 scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
 684
 685                 F[coarseFaceI][visCoarseFaceI] = Fij;
 686                 sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
 687             }
 688         }
 689     }
 690
 691     if (Pstream::master())
 692     {
 693         Info << "Writing view factor matrix..." << endl;
 694     }
 695
 696     // Write view factors matrix in listlist form
 697     F.write();
 698
 699     reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
 700     reduce(patchArea, sumOp<scalarList>());
 701
 702
 703     if (Pstream::master() && debug)
 704     {
 705         forAll(viewFactorsPatches, i)
 706         {
 707             label patchI =  viewFactorsPatches[i];
 708             forAll(viewFactorsPatches, i)
 709             {
 710                 label patchJ =  viewFactorsPatches[i];
 711                 Info << "F" << patchI << patchJ << ": "
 712                      << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
 713                      << endl;
 714             }
 715         }
 716     }
 717
 718
 719     if (writeViewFactors)
 720     {
 721         volScalarField viewFactorField
 722         (
 723             IOobject
 724             (
 725                 "viewFactorField",
 726                 mesh.time().timeName(),
 727                 mesh,
 728                 IOobject::NO_READ,
 729                 IOobject::NO_WRITE
 730             ),
 731             mesh,
 732             dimensionedScalar("viewFactorField", dimless, 0)
 733         );
 734
 735         label compactI = 0;
 736         forAll(viewFactorsPatches, i)
 737         {
 738             label patchID = viewFactorsPatches[i];
 739             const labelList& agglom = finalAgglom[patchID];
 740             label nAgglom = max(agglom)+1;
 741             labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
 742             const labelList& coarsePatchFace =
 743                 coarseMesh.patchFaceMap()[patchID];
 744
 745             forAll(coarseToFine, coarseI)
 746             {
 747                 const scalar Fij = sum(F[compactI]);
 748                 const label coarseFaceID = coarsePatchFace[coarseI];
 749                 const labelList& fineFaces = coarseToFine[coarseFaceID];
 750                 forAll (fineFaces, fineId)
 751                 {
 752                     const label faceID = fineFaces[fineId];
 753                     viewFactorField.boundaryField()[patchID][faceID] = Fij;
 754                 }
 755                 compactI++;
 756             }
 757         }
 758         viewFactorField.write();
 759     }
 760
 761
 762     // Invert compactMap (from processor+localface to compact) to go
 763     // from compact to processor+localface (expressed as a globalIndex)
 764     // globalIndex globalCoarFaceNum(coarseMesh.nFaces());
 765     labelList compactToGlobal(map.constructSize());
 766
 767     // Local indices first (note: are not in compactMap)
 768     for (label i = 0; i < globalNumbering.localSize(); i++)
 769     {
 770         compactToGlobal[i] = globalNumbering.toGlobal(i);
 771     }
 772
 773
 774     forAll(compactMap, procI)
 775     {
 776         const Map<label>& localToCompactMap = compactMap[procI];
 777
 778         forAllConstIter(Map<label>, localToCompactMap, iter)
 779         {
 780             compactToGlobal[iter()] = globalNumbering.toGlobal
 781             (
 782                 procI,
 783                 iter.key()
 784             );
 785         }
 786     }
 787
 788
 789     if (Pstream::master())
 790     {
 791         scalarSquareMatrix Fmatrix(totalNCoarseFaces, totalNCoarseFaces, 0.0);
 792
 793         labelListList globalFaceFaces(visibleFaceFaces.size());
 794
 795         // Create globalFaceFaces needed to insert view factors
 796         // in F to the global matrix Fmatrix
 797         forAll(globalFaceFaces, faceI)
 798         {
 799             globalFaceFaces[faceI] = renumber
 800             (
 801                 compactToGlobal,
 802                 visibleFaceFaces[faceI]
 803             );
 804         }
 805
 806         labelListIOList IOglobalFaceFaces
 807         (
 808             IOobject
 809             (
 810                 "globalFaceFaces",
 811                 mesh.facesInstance(),
 812                 mesh,
 813                 IOobject::NO_READ,
 814                 IOobject::NO_WRITE,
 815                 false
 816             ),
 817             globalFaceFaces
 818         );
 819         IOglobalFaceFaces.write();
 820     }
 821     else
 822     {
 823         labelListList globalFaceFaces(visibleFaceFaces.size());
 824         forAll(globalFaceFaces, faceI)
 825         {
 826             globalFaceFaces[faceI] = renumber
 827             (
 828                 compactToGlobal,
 829                 visibleFaceFaces[faceI]
 830             );
 831         }
 832
 833         labelListIOList IOglobalFaceFaces
 834         (
 835             IOobject
 836             (
 837                 "globalFaceFaces",
 838                 mesh.facesInstance(),
 839                 mesh,
 840                 IOobject::NO_READ,
 841                 IOobject::NO_WRITE,
 842                 false
 843             ),
 844             globalFaceFaces
 845         );
 846
 847         IOglobalFaceFaces.write();
 848     }
 849
 850     Info<< "End\n" << endl;
 851     return 0;
 852 }
 853
 854
 855 // ************************************************************************* //