applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "faMeshDecomposition.H"
  27 #include "foamTime.H"
  28 #include "dictionary.H"
  29 #include "labelIOList.H"
  30 #include "processorFaPatch.H"
  31 #include "faMesh.H"
  32 #include "OSspecific.H"
  33 #include "Map.H"
  34 #include "globalMeshData.H"
  35
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 void faMeshDecomposition::distributeFaces()
  40 {
  41     Info<< "\nCalculating distribution of faces" << endl;
  42
  43     cpuTime decompositionTime;
  44
  45     for (label procI = 0; procI < nProcs(); procI++)
  46     {
  47         Time processorDb
  48         (
  49             Time::controlDictName,
  50             time().rootPath(),
  51             time().caseName()/fileName(word("processor") + Foam::name(procI))
  52         );
  53
  54         fvMesh procMesh
  55         (
  56             IOobject
  57             (
  58                 fvMesh::defaultRegion,
  59                 processorDb.timeName(),
  60                 processorDb
  61             )
  62         );
  63
  64         labelHashSet faceProcAddressingHash
  65         (
  66             labelIOList
  67             (
  68                 IOobject
  69                 (
  70                     "faceProcAddressing",
  71                     "constant",
  72                     procMesh.meshSubDir,
  73                     procMesh,
  74                     IOobject::MUST_READ,
  75                     IOobject::NO_WRITE
  76                 )
  77             )
  78         );
  79
  80         forAll (faceLabels(), faceI)
  81         {
  82             if (faceProcAddressingHash.found(faceLabels()[faceI] + 1))
  83             {
  84                 faceToProc_[faceI] = procI;
  85             }
  86         }
  87     }
  88
  89     Info<< "\nFinished decomposition in "
  90         << decompositionTime.elapsedCpuTime()
  91         << " s" << endl;
  92 }
  93
  94 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  95
  96 // from components
  97 faMeshDecomposition::faMeshDecomposition(const fvMesh& mesh)
  98 :
  99     faMesh(mesh),
 100     decompositionDict_
 101     (
 102         IOobject
 103         (
 104             "decomposeParDict",
 105             time().system(),
 106             mesh,
 107             IOobject::MUST_READ,
 108             IOobject::NO_WRITE
 109         )
 110     ),
 111     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
 112     distributed_(false),
 113     faceToProc_(nFaces()),
 114     procFaceLabels_(nProcs_),
 115     procMeshEdgesMap_(nProcs_),
 116     procNInternalEdges_(nProcs_, 0),
 117     procPatchEdgeLabels_(nProcs_),
 118     procPatchPointAddressing_(nProcs_),
 119     procPatchEdgeAddressing_(nProcs_),
 120     procEdgeAddressing_(nProcs_),
 121     procFaceAddressing_(nProcs_),
 122     procBoundaryAddressing_(nProcs_),
 123     procPatchSize_(nProcs_),
 124     procPatchStartIndex_(nProcs_),
 125     procNeighbourProcessors_(nProcs_),
 126     procProcessorPatchSize_(nProcs_),
 127     procProcessorPatchStartIndex_(nProcs_),
 128     globallySharedPoints_(0),
 129     cyclicParallel_(false)
 130 {
 131     if (decompositionDict_.found("distributed"))
 132     {
 133         distributed_ = Switch(decompositionDict_.lookup("distributed"));
 134     }
 135 }
 136
 137
 138 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 139
 140 faMeshDecomposition::~faMeshDecomposition()
 141 {}
 142
 143
 144 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 145
 146 void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
 147 {
 148     // Decide which cell goes to which processor
 149     distributeFaces();
 150
 151     Info<< "\nDistributing faces to processors" << endl;
 152
 153     // Memory management
 154     {
 155         List<SLList<label> > procFaceList(nProcs());
 156
 157         forAll (faceToProc_, faceI)
 158         {
 159             if (faceToProc_[faceI] >= nProcs())
 160             {
 161                 FatalErrorIn("Finite area mesh decomposition")
 162                     << "Impossible processor label " << faceToProc_[faceI]
 163                     << "for face " << faceI
 164                     << abort(FatalError);
 165             }
 166             else
 167             {
 168                 procFaceList[faceToProc_[faceI]].append(faceI);
 169             }
 170         }
 171
 172         // Convert linked lists into normal lists
 173         forAll (procFaceList, procI)
 174         {
 175             procFaceAddressing_[procI] = procFaceList[procI];
 176         }
 177     }
 178
 179
 180     // Find processor mesh faceLabels and ...
 181
 182     for (label procI = 0; procI < nProcs(); procI++)
 183     {
 184         Time processorDb
 185         (
 186             Time::controlDictName,
 187             time().rootPath(),
 188             time().caseName()/fileName(word("processor") + Foam::name(procI))
 189         );
 190
 191         fvMesh procFvMesh
 192         (
 193             IOobject
 194             (
 195                 fvMesh::defaultRegion,
 196                 processorDb.timeName(),
 197                 processorDb
 198             )
 199         );
 200
 201         labelIOList fvPointProcAddressing
 202         (
 203             IOobject
 204             (
 205                 "pointProcAddressing",
 206                 "constant",
 207                 procFvMesh.meshSubDir,
 208                 procFvMesh,
 209                 IOobject::MUST_READ,
 210                 IOobject::NO_WRITE
 211             )
 212         );
 213
 214         HashTable<label, label, Hash<label> > fvFaceProcAddressingHash;
 215
 216         {
 217             labelIOList fvFaceProcAddressing
 218             (
 219                 IOobject
 220                 (
 221                     "faceProcAddressing",
 222                     "constant",
 223                     procFvMesh.meshSubDir,
 224                     procFvMesh,
 225                     IOobject::MUST_READ,
 226                     IOobject::NO_WRITE
 227                 )
 228             );
 229
 230             forAll(fvFaceProcAddressing, faceI)
 231             {
 232                  fvFaceProcAddressingHash.insert
 233                  (
 234                      fvFaceProcAddressing[faceI], faceI
 235                  );
 236             }
 237         };
 238
 239         const labelList& curProcFaceAddressing = procFaceAddressing_[procI];
 240
 241         labelList& curFaceLabels = procFaceLabels_[procI];
 242
 243         curFaceLabels = labelList(curProcFaceAddressing.size(), -1);
 244
 245         forAll(curProcFaceAddressing, faceI)
 246         {
 247             curFaceLabels[faceI] =
 248                 fvFaceProcAddressingHash.find
 249                 (
 250                     faceLabels()[curProcFaceAddressing[faceI]] + 1
 251                 )();
 252         }
 253
 254         // create processor finite area mesh
 255         faMesh procMesh
 256         (
 257             procFvMesh,
 258             procFaceLabels_[procI]
 259         );
 260
 261         const indirectPrimitivePatch& patch = this->patch();
 262         const Map<label>& map = patch.meshPointMap();
 263
 264         HashTable<label, edge, Hash<edge> > edgesHash;
 265
 266         label edgeI = -1;
 267
 268         label nIntEdges = patch.nInternalEdges();
 269
 270         for (label curEdge = 0; curEdge < nIntEdges; curEdge++)
 271         {
 272             edgesHash.insert(patch.edges()[curEdge], ++edgeI);
 273         }
 274
 275         forAll (boundary(), patchI)
 276         {
 277             // Include emptyFaPatch
 278
 279             label size = boundary()[patchI].labelList::size();
 280
 281             for(int eI=0; eI<size; eI++)
 282             {
 283                 edgesHash.insert(patch.edges()[boundary()[patchI][eI]], ++edgeI);
 284             }
 285         }
 286
 287
 288         const indirectPrimitivePatch& procPatch = procMesh.patch();
 289         const vectorField& procPoints = procPatch.localPoints();
 290         const labelList& procMeshPoints = procPatch.meshPoints();
 291         const edgeList& procEdges = procPatch.edges();
 292
 293         labelList& curPatchPointAddressing = procPatchPointAddressing_[procI];
 294         curPatchPointAddressing.setSize(procPoints.size(), -1);
 295
 296         forAll(procPoints, pointI)
 297         {
 298             curPatchPointAddressing[pointI] =
 299                 map[fvPointProcAddressing[procMeshPoints[pointI]]];
 300         }
 301
 302         labelList& curPatchEdgeAddressing = procPatchEdgeAddressing_[procI];
 303         curPatchEdgeAddressing.setSize(procEdges.size(), -1);
 304
 305         forAll(procEdges, edgeI)
 306         {
 307             edge curGlobalEdge = procEdges[edgeI];
 308             curGlobalEdge[0] = curPatchPointAddressing[curGlobalEdge[0]];
 309             curGlobalEdge[1] = curPatchPointAddressing[curGlobalEdge[1]];
 310
 311             curPatchEdgeAddressing[edgeI] = edgesHash.find(curGlobalEdge)();
 312         }
 313
 314         Map<label>& curMap = procMeshEdgesMap_[procI];
 315
 316         curMap = Map<label>(2*procEdges.size());
 317
 318         forAll(curPatchEdgeAddressing, edgeI)
 319         {
 320             curMap.insert(curPatchEdgeAddressing[edgeI], edgeI);
 321         }
 322
 323         procNInternalEdges_[procI] = procPatch.nInternalEdges();
 324     }
 325
 326
 327     Info << "\nDistributing edges to processors" << endl;
 328
 329     // Loop through all internal edges and decide which processor they
 330     // belong to. First visit all internal edges.
 331
 332     // set references to the original mesh
 333     const faBoundaryMesh& patches = boundary();
 334     const edgeList& edges = this->edges();
 335     const labelList& owner = edgeOwner();
 336     const labelList& neighbour = edgeNeighbour();
 337
 338     // Memory management
 339     {
 340         List<SLList<label> > procEdgeList(nProcs());
 341
 342         forAll(procEdgeList, procI)
 343         {
 344             for(label i=0; i<procNInternalEdges_[procI]; i++)
 345             {
 346                 procEdgeList[procI].append(procPatchEdgeAddressing_[procI][i]);
 347             }
 348         }
 349
 350
 351         // Detect inter-processor boundaries
 352
 353         // Neighbour processor for each subdomain
 354         List<SLList<label> > interProcBoundaries(nProcs());
 355
 356         // Edge labels belonging to each inter-processor boundary
 357         List<SLList<SLList<label> > > interProcBEdges(nProcs());
 358
 359         List<SLList<label> > procPatchIndex(nProcs());
 360
 361         forAll (neighbour, edgeI)
 362         {
 363             if (faceToProc_[owner[edgeI]] != faceToProc_[neighbour[edgeI]])
 364             {
 365                 // inter - processor patch edge found. Go through the list of
 366                 // inside boundaries for the owner processor and try to find
 367                 // this inter-processor patch.
 368
 369                 label ownerProc = faceToProc_[owner[edgeI]];
 370                 label neighbourProc = faceToProc_[neighbour[edgeI]];
 371
 372                 SLList<label>::iterator curInterProcBdrsOwnIter =
 373                     interProcBoundaries[ownerProc].begin();
 374
 375                 SLList<SLList<label> >::iterator curInterProcBEdgesOwnIter =
 376                     interProcBEdges[ownerProc].begin();
 377
 378                 bool interProcBouFound = false;
 379
 380                 // WARNING: Synchronous SLList iterators
 381
 382                 for
 383                 (
 384                     ;
 385                     curInterProcBdrsOwnIter
 386                  != interProcBoundaries[ownerProc].end()
 387                  && curInterProcBEdgesOwnIter
 388                  != interProcBEdges[ownerProc].end();
 389                     ++curInterProcBdrsOwnIter, ++curInterProcBEdgesOwnIter
 390                 )
 391                 {
 392                     if (curInterProcBdrsOwnIter() == neighbourProc)
 393                     {
 394                         // the inter - processor boundary exists. Add the face
 395                         interProcBouFound = true;
 396
 397                         curInterProcBEdgesOwnIter().append(edgeI);
 398
 399                         SLList<label>::iterator curInterProcBdrsNeiIter =
 400                             interProcBoundaries[neighbourProc].begin();
 401
 402                         SLList<SLList<label> >::iterator
 403                             curInterProcBEdgesNeiIter =
 404                             interProcBEdges[neighbourProc].begin();
 405
 406                         bool neighbourFound = false;
 407
 408                         // WARNING: Synchronous SLList iterators
 409
 410                         for
 411                         (
 412                             ;
 413                             curInterProcBdrsNeiIter !=
 414                             interProcBoundaries[neighbourProc].end()
 415                          && curInterProcBEdgesNeiIter !=
 416                             interProcBEdges[neighbourProc].end();
 417                             ++curInterProcBdrsNeiIter,
 418                             ++curInterProcBEdgesNeiIter
 419                         )
 420                         {
 421                             if (curInterProcBdrsNeiIter() == ownerProc)
 422                             {
 423                                 // boundary found. Add the face
 424                                 neighbourFound = true;
 425
 426                                 curInterProcBEdgesNeiIter().append(edgeI);
 427                             }
 428
 429                             if (neighbourFound) break;
 430                         }
 431
 432                         if (interProcBouFound && !neighbourFound)
 433                         {
 434                             FatalErrorIn
 435                                 ("faDomainDecomposition::decomposeMesh()")
 436                                 << "Inconsistency in inter - "
 437                                 << "processor boundary lists for processors "
 438                                 << ownerProc << " and " << neighbourProc
 439                                 << abort(FatalError);
 440                         }
 441                     }
 442
 443                     if (interProcBouFound) break;
 444                 }
 445
 446                 if (!interProcBouFound)
 447                 {
 448                     // inter - processor boundaries do not exist and need to
 449                     // be created
 450
 451                     // set the new addressing information
 452
 453                     // owner
 454                     interProcBoundaries[ownerProc].append(neighbourProc);
 455                     interProcBEdges[ownerProc].append(SLList<label>(edgeI));
 456
 457                     // neighbour
 458                     interProcBoundaries[neighbourProc].append(ownerProc);
 459                     interProcBEdges[neighbourProc].append
 460                         (
 461                             SLList<label>(edgeI)
 462                         );
 463                 }
 464             }
 465         }
 466
 467
 468         // Loop through patches. For cyclic boundaries detect inter-processor
 469         // edges; for all other, add edges to the edge list and remember start
 470         // and size of all patches.
 471
 472         // for all processors, set the size of start index and patch size
 473         // lists to the number of patches in the mesh
 474         forAll (procPatchSize_, procI)
 475         {
 476             procPatchSize_[procI].setSize(patches.size());
 477             procPatchStartIndex_[procI].setSize(patches.size());
 478         }
 479
 480         forAll (patches, patchI)
 481         {
 482             // Reset size and start index for all processors
 483             forAll (procPatchSize_, procI)
 484             {
 485                 procPatchSize_[procI][patchI] = 0;
 486                 procPatchStartIndex_[procI][patchI] =
 487                     procEdgeList[procI].size();
 488             }
 489
 490             const label patchStart = patches[patchI].start();
 491
 492 //             if (typeid(patches[patchI]) != typeid(cyclicFaPatch))
 493             if (true)
 494             {
 495                 // Normal patch. Add edges to processor where the face
 496                 // next to the edge lives
 497
 498                 const labelListList& eF = patch().edgeFaces();
 499
 500                 label size = patches[patchI].labelList::size();
 501
 502                 labelList patchEdgeFaces(size, -1);
 503
 504                 for(int eI=0; eI<size; eI++)
 505                 {
 506                     patchEdgeFaces[eI] = eF[patches[patchI][eI]][0];
 507                 }
 508
 509                 forAll (patchEdgeFaces, edgeI)
 510                 {
 511                     const label curProc = faceToProc_[patchEdgeFaces[edgeI]];
 512
 513                     // add the face
 514                     procEdgeList[curProc].append(patchStart + edgeI);
 515
 516                     // increment the number of edges for this patch
 517                     procPatchSize_[curProc][patchI]++;
 518                 }
 519             }
 520             else
 521             {
 522                 // Cyclic patch special treatment
 523
 524                 const faPatch& cPatch = patches[patchI];
 525
 526                 const label cycOffset = cPatch.size()/2;
 527
 528                 // Set reference to faceCells for both patches
 529                 const labelList::subList firstEdgeFaces
 530                 (
 531                     cPatch.edgeFaces(),
 532                     cycOffset
 533                 );
 534
 535                 const labelList::subList secondEdgeFaces
 536                 (
 537                     cPatch.edgeFaces(),
 538                     cycOffset,
 539                     cycOffset
 540                 );
 541
 542                 forAll (firstEdgeFaces, edgeI)
 543                 {
 544                     if
 545                     (
 546                         faceToProc_[firstEdgeFaces[edgeI]]
 547                      != faceToProc_[secondEdgeFaces[edgeI]]
 548                     )
 549                     {
 550                         // This edge becomes an inter-processor boundary edge
 551                         // inter - processor patch edge found. Go through
 552                         // the list of inside boundaries for the owner
 553                         // processor and try to find this inter-processor
 554                         // patch.
 555
 556                         cyclicParallel_ = true;
 557
 558                         label ownerProc = faceToProc_[firstEdgeFaces[edgeI]];
 559                         label neighbourProc =
 560                             faceToProc_[secondEdgeFaces[edgeI]];
 561
 562                         SLList<label>::iterator curInterProcBdrsOwnIter =
 563                             interProcBoundaries[ownerProc].begin();
 564
 565                         SLList<SLList<label> >::iterator
 566                             curInterProcBEdgesOwnIter =
 567                             interProcBEdges[ownerProc].begin();
 568
 569                         bool interProcBouFound = false;
 570
 571                         // WARNING: Synchronous SLList iterators
 572
 573                         for
 574                         (
 575                             ;
 576                             curInterProcBdrsOwnIter !=
 577                             interProcBoundaries[ownerProc].end()
 578                          && curInterProcBEdgesOwnIter !=
 579                             interProcBEdges[ownerProc].end();
 580                             ++curInterProcBdrsOwnIter,
 581                             ++curInterProcBEdgesOwnIter
 582                         )
 583                         {
 584                             if (curInterProcBdrsOwnIter() == neighbourProc)
 585                             {
 586                                 // the inter - processor boundary exists.
 587                                 // Add the face
 588                                 interProcBouFound = true;
 589
 590                                 curInterProcBEdgesOwnIter().append
 591                                     (patchStart + edgeI);
 592
 593                                 SLList<label>::iterator curInterProcBdrsNeiIter
 594                                    = interProcBoundaries[neighbourProc].begin();
 595
 596                                 SLList<SLList<label> >::iterator
 597                                     curInterProcBEdgesNeiIter =
 598                                     interProcBEdges[neighbourProc].begin();
 599
 600                                 bool neighbourFound = false;
 601
 602                                 // WARNING: Synchronous SLList iterators
 603
 604                                 for
 605                                 (
 606                                     ;
 607                                     curInterProcBdrsNeiIter
 608                                    != interProcBoundaries[neighbourProc].end()
 609                                  && curInterProcBEdgesNeiIter
 610                                    != interProcBEdges[neighbourProc].end();
 611                                     ++curInterProcBdrsNeiIter,
 612                                     ++curInterProcBEdgesNeiIter
 613                                 )
 614                                 {
 615                                     if (curInterProcBdrsNeiIter() == ownerProc)
 616                                     {
 617                                         // boundary found. Add the face
 618                                         neighbourFound = true;
 619
 620                                         curInterProcBEdgesNeiIter()
 621                                            .append
 622                                             (
 623                                                 patchStart
 624                                               + cycOffset
 625                                               + edgeI
 626                                             );
 627                                     }
 628
 629                                     if (neighbourFound) break;
 630                                 }
 631
 632                                 if (interProcBouFound && !neighbourFound)
 633                                 {
 634                                     FatalErrorIn
 635                                     (
 636                                         "faDomainDecomposition::decomposeMesh()"
 637                                     )   << "Inconsistency in inter-processor "
 638                                         << "boundary lists for processors "
 639                                         << ownerProc << " and " << neighbourProc
 640                                         << " in cyclic boundary matching"
 641                                         << abort(FatalError);
 642                                 }
 643                             }
 644
 645                             if (interProcBouFound) break;
 646                         }
 647
 648                         if (!interProcBouFound)
 649                         {
 650                             // inter - processor boundaries do not exist
 651                             // and need to be created
 652
 653                             // set the new addressing information
 654
 655                             // owner
 656                             interProcBoundaries[ownerProc]
 657                                 .append(neighbourProc);
 658                             interProcBEdges[ownerProc]
 659                                 .append(SLList<label>(patchStart + edgeI));
 660
 661                             // neighbour
 662                             interProcBoundaries[neighbourProc]
 663                                .append(ownerProc);
 664                             interProcBEdges[neighbourProc]
 665                                .append
 666                                 (
 667                                     SLList<label>
 668                                     (
 669                                         patchStart
 670                                       + cycOffset
 671                                       + edgeI
 672                                     )
 673                                 );
 674                         }
 675                     }
 676                     else
 677                     {
 678                         // This cyclic edge remains on the processor
 679                         label ownerProc = faceToProc_[firstEdgeFaces[edgeI]];
 680
 681                         // add the edge
 682                         procEdgeList[ownerProc].append(patchStart + edgeI);
 683
 684                         // increment the number of edges for this patch
 685                         procPatchSize_[ownerProc][patchI]++;
 686
 687                         // Note: I cannot add the other side of the cyclic
 688                         // boundary here because this would violate the order.
 689                         // They will be added in a separate loop below
 690                     }
 691                 }
 692
 693                 // Ordering in cyclic boundaries is important.
 694                 // Add the other half of cyclic edges for cyclic boundaries
 695                 // that remain on the processor
 696                 forAll (secondEdgeFaces, edgeI)
 697                 {
 698                     if
 699                     (
 700                         faceToProc_[firstEdgeFaces[edgeI]]
 701                      == faceToProc_[secondEdgeFaces[edgeI]]
 702                     )
 703                     {
 704                         // This cyclic edge remains on the processor
 705                         label ownerProc = faceToProc_[firstEdgeFaces[edgeI]];
 706
 707                         // add the second edge
 708                         procEdgeList[ownerProc].append
 709                             (patchStart + cycOffset + edgeI);
 710
 711                         // increment the number of edges for this patch
 712                         procPatchSize_[ownerProc][patchI]++;
 713                     }
 714                 }
 715             }
 716         }
 717
 718         // Convert linked lists into normal lists
 719         // Add inter-processor boundaries and remember start indices
 720         forAll (procEdgeList, procI)
 721         {
 722             // Get internal and regular boundary processor faces
 723             SLList<label>& curProcEdges = procEdgeList[procI];
 724
 725             // Get reference to processor edge addressing
 726             labelList& curProcEdgeAddressing = procEdgeAddressing_[procI];
 727
 728             labelList& curProcNeighbourProcessors =
 729                 procNeighbourProcessors_[procI];
 730
 731             labelList& curProcProcessorPatchSize =
 732                 procProcessorPatchSize_[procI];
 733
 734             labelList& curProcProcessorPatchStartIndex =
 735                 procProcessorPatchStartIndex_[procI];
 736
 737             // calculate the size
 738             label nEdgesOnProcessor = curProcEdges.size();
 739
 740             for
 741             (
 742                 SLList<SLList<label> >::iterator curInterProcBEdgesIter =
 743                     interProcBEdges[procI].begin();
 744                 curInterProcBEdgesIter != interProcBEdges[procI].end();
 745                 ++curInterProcBEdgesIter
 746             )
 747             {
 748                 nEdgesOnProcessor += curInterProcBEdgesIter().size();
 749             }
 750
 751             curProcEdgeAddressing.setSize(nEdgesOnProcessor);
 752
 753             // Fill in the list. Calculate turning index.
 754             // Turning index will be -1 only for some edges on processor
 755             // boundaries, i.e. the ones where the current processor ID
 756             // is in the face which is a edge neighbour.
 757             // Turning index is stored as the sign of the edge addressing list
 758
 759             label nEdges = 0;
 760
 761             // Add internal and boundary edges
 762             // Remember to increment the index by one such that the
 763             // turning index works properly.
 764             for
 765             (
 766                 SLList<label>::iterator curProcEdgeIter = curProcEdges.begin();
 767                 curProcEdgeIter != curProcEdges.end();
 768                 ++curProcEdgeIter
 769             )
 770             {
 771                 curProcEdgeAddressing[nEdges] = curProcEdgeIter();
 772 //                 curProcEdgeAddressing[nEdges] = curProcEdgeIter() + 1;
 773                 nEdges++;
 774             }
 775
 776             // Add inter-processor boundary edges. At the beginning of each
 777             // patch, grab the patch start index and size
 778
 779             curProcNeighbourProcessors.setSize
 780             (
 781                 interProcBoundaries[procI].size()
 782             );
 783
 784             curProcProcessorPatchSize.setSize
 785             (
 786                 interProcBoundaries[procI].size()
 787             );
 788
 789             curProcProcessorPatchStartIndex.setSize
 790             (
 791                 interProcBoundaries[procI].size()
 792             );
 793
 794             label nProcPatches = 0;
 795
 796             SLList<label>::iterator curInterProcBdrsIter =
 797                 interProcBoundaries[procI].begin();
 798
 799             SLList<SLList<label> >::iterator curInterProcBEdgesIter =
 800                 interProcBEdges[procI].begin();
 801
 802             for
 803             (
 804                 ;
 805                 curInterProcBdrsIter != interProcBoundaries[procI].end()
 806              && curInterProcBEdgesIter != interProcBEdges[procI].end();
 807                 ++curInterProcBdrsIter, ++curInterProcBEdgesIter
 808             )
 809             {
 810                 curProcNeighbourProcessors[nProcPatches] =
 811                     curInterProcBdrsIter();
 812
 813                 // Get start index for processor patch
 814                 curProcProcessorPatchStartIndex[nProcPatches] = nEdges;
 815
 816                 label& curSize =
 817                     curProcProcessorPatchSize[nProcPatches];
 818
 819                 curSize = 0;
 820
 821                 // add faces for this processor boundary
 822
 823                 for
 824                 (
 825                     SLList<label>::iterator curEdgesIter =
 826                         curInterProcBEdgesIter().begin();
 827                     curEdgesIter != curInterProcBEdgesIter().end();
 828                     ++curEdgesIter
 829                 )
 830                 {
 831                     // add the edges
 832
 833                     // Remember to increment the index by one such that the
 834                     // turning index works properly.
 835                     if (faceToProc_[owner[curEdgesIter()]] == procI)
 836                     {
 837                         curProcEdgeAddressing[nEdges] = curEdgesIter();
 838 //                         curProcEdgeAddressing[nEdges] = curEdgesIter() + 1;
 839                     }
 840                     else
 841                     {
 842                         // turning edge
 843                         curProcEdgeAddressing[nEdges] = curEdgesIter();
 844 //                         curProcEdgeAddressing[nEdges] = -(curEdgesIter() + 1);
 845                     }
 846
 847                     // increment the size
 848                     curSize++;
 849
 850                     nEdges++;
 851                 }
 852
 853                 nProcPatches++;
 854             }
 855         }
 856     }
 857
 858     Info << "\nCalculating processor boundary addressing" << endl;
 859     // For every patch of processor boundary, find the index of the original
 860     // patch. Mis-alignment is caused by the fact that patches with zero size
 861     // are omitted. For processor patches, set index to -1.
 862     // At the same time, filter the procPatchSize_ and procPatchStartIndex_
 863     // lists to exclude zero-size patches
 864     forAll (procPatchSize_, procI)
 865     {
 866         // Make a local copy of old lists
 867         const labelList oldPatchSizes = procPatchSize_[procI];
 868
 869         const labelList oldPatchStarts = procPatchStartIndex_[procI];
 870
 871         labelList& curPatchSizes = procPatchSize_[procI];
 872
 873         labelList& curPatchStarts = procPatchStartIndex_[procI];
 874
 875         const labelList& curProcessorPatchSizes =
 876             procProcessorPatchSize_[procI];
 877
 878         labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
 879
 880         curBoundaryAddressing.setSize
 881         (
 882             oldPatchSizes.size()
 883           + curProcessorPatchSizes.size()
 884         );
 885
 886         label nPatches = 0;
 887
 888         forAll (oldPatchSizes, patchI)
 889         {
 890             if (!filterEmptyPatches || oldPatchSizes[patchI] > 0)
 891             {
 892                 curBoundaryAddressing[nPatches] = patchI;
 893
 894                 curPatchSizes[nPatches] = oldPatchSizes[patchI];
 895
 896                 curPatchStarts[nPatches] = oldPatchStarts[patchI];
 897
 898                 nPatches++;
 899             }
 900         }
 901
 902         // reset to the size of live patches
 903         curPatchSizes.setSize(nPatches);
 904         curPatchStarts.setSize(nPatches);
 905
 906         forAll (curProcessorPatchSizes, procPatchI)
 907         {
 908             curBoundaryAddressing[nPatches] = -1;
 909
 910             nPatches++;
 911         }
 912
 913         curBoundaryAddressing.setSize(nPatches);
 914     }
 915
 916
 917     // Gather data about globally shared points
 918
 919     labelList globallySharedPoints_(0);
 920
 921     // Memory management
 922     {
 923         labelList pointsUsage(nPoints(), 0);
 924
 925         // Globally shared points are the ones used by more than 2 processors
 926         // Size the list approximately and gather the points
 927         labelHashSet gSharedPoints
 928         (
 929             min(100, nPoints()/1000)
 930         );
 931
 932         // Loop through all the processors and mark up points used by
 933         // processor boundaries.  When a point is used twice, it is a
 934         // globally shared point
 935
 936         for (label procI = 0; procI < nProcs(); procI++)
 937         {
 938             // Get list of edge labels
 939             const labelList& curEdgeLabels = procEdgeAddressing_[procI];
 940
 941             // Get start of processor faces
 942             const labelList& curProcessorPatchStarts =
 943                 procProcessorPatchStartIndex_[procI];
 944
 945             const labelList& curProcessorPatchSizes =
 946                 procProcessorPatchSize_[procI];
 947
 948             // Reset the lookup list
 949             pointsUsage = 0;
 950
 951             forAll (curProcessorPatchStarts, patchI)
 952             {
 953                 const label curStart = curProcessorPatchStarts[patchI];
 954                 const label curEnd = curStart + curProcessorPatchSizes[patchI];
 955
 956                 for
 957                 (
 958                     label edgeI = curStart;
 959                     edgeI < curEnd;
 960                     edgeI++
 961                 )
 962                 {
 963                     // Mark the original edge as used
 964                     // Remember to decrement the index by one (turning index)
 965                     const label curE = curEdgeLabels[edgeI];
 966
 967                     const edge& e = edges[curE];
 968
 969                     forAll (e, pointI)
 970                     {
 971                         if (pointsUsage[e[pointI]] == 0)
 972                         {
 973                             // Point not previously used
 974                             pointsUsage[e[pointI]] = patchI + 1;
 975                         }
 976                         else if (pointsUsage[e[pointI]] != patchI + 1)
 977                         {
 978                             // Point used by some other patch = global point!
 979                             gSharedPoints.insert(e[pointI]);
 980                         }
 981                     }
 982                 }
 983             }
 984         }
 985
 986         // Grab the result from the hash list
 987         globallySharedPoints_ = gSharedPoints.toc();
 988         sort(globallySharedPoints_);
 989     }
 990
 991
 992     // Edge label for faPatches
 993
 994     for (label procI = 0; procI < nProcs(); procI++)
 995     {
 996         fileName processorCasePath
 997         (
 998             time().caseName()/fileName(word("processor")
 999           + Foam::name(procI))
1000         );
1001
1002         // create a database
1003         Time processorDb
1004         (
1005             Time::controlDictName,
1006             time().rootPath(),
1007             processorCasePath
1008         );
1009
1010
1011         // read finite volume mesh
1012         fvMesh procFvMesh
1013         (
1014             IOobject
1015             (
1016                 fvMesh::defaultRegion,
1017                 processorDb.timeName(),
1018                 processorDb
1019             )
1020         );
1021
1022         // create finite area mesh
1023         faMesh procMesh
1024         (
1025             procFvMesh,
1026             procFaceLabels_[procI]
1027         );
1028
1029
1030         const labelList& curEdgeAddressing = procEdgeAddressing_[procI];
1031
1032         const labelList& curPatchStartIndex = procPatchStartIndex_[procI];
1033         const labelList& curPatchSize = procPatchSize_[procI];
1034
1035         const labelList& curProcessorPatchStartIndex =
1036             procProcessorPatchStartIndex_[procI];
1037
1038         const labelList& curProcessorPatchSize =
1039             procProcessorPatchSize_[procI];
1040
1041         labelListList& curPatchEdgeLabels = procPatchEdgeLabels_[procI];
1042         curPatchEdgeLabels =
1043             labelListList
1044             (
1045                 curPatchSize.size()
1046               + curProcessorPatchSize.size()
1047             );
1048
1049         forAll(curPatchSize, patchI)
1050         {
1051             labelList& curEdgeLabels = curPatchEdgeLabels[patchI];
1052             curEdgeLabels.setSize(curPatchSize[patchI], -1);
1053
1054             label edgeI = 0;
1055
1056             for
1057             (
1058                 int i=curPatchStartIndex[patchI];
1059                 i<(curPatchStartIndex[patchI]+curPatchSize[patchI]);
1060                 i++
1061             )
1062             {
1063                 curEdgeLabels[edgeI] =
1064                     procMeshEdgesMap_[procI][curEdgeAddressing[i]];
1065                 edgeI++;
1066             }
1067         }
1068
1069         forAll(curProcessorPatchSize, patchI)
1070         {
1071             labelList& curEdgeLabels  =
1072                 curPatchEdgeLabels[curPatchSize.size() + patchI];
1073             curEdgeLabels.setSize(curProcessorPatchSize[patchI], -1);
1074
1075             label edgeI = 0;
1076
1077             for
1078             (
1079                 int i=curProcessorPatchStartIndex[patchI];
1080                 i<(curProcessorPatchStartIndex[patchI]
1081                 +curProcessorPatchSize[patchI]);
1082                 i++
1083             )
1084             {
1085                 curEdgeLabels[edgeI] =
1086                     procMeshEdgesMap_[procI][curEdgeAddressing[i]];
1087                 edgeI++;
1088             }
1089         }
1090     }
1091 }
1092
1093
1094 bool faMeshDecomposition::writeDecomposition()
1095 {
1096     Info<< "\nConstructing processor FA meshes" << endl;
1097
1098
1099     // Make a lookup map for globally shared points
1100     Map<label> sharedPointLookup(2*globallySharedPoints_.size());
1101
1102     forAll (globallySharedPoints_, pointi)
1103     {
1104         sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
1105     }
1106
1107     label totProcEdges = 0;
1108     label maxProcPatches = 0;
1109     label maxProcEdges = 0;
1110
1111     // Write out the meshes
1112     for (label procI = 0; procI < nProcs(); procI++)
1113     {
1114         // Create processor mesh without a boundary
1115
1116         fileName processorCasePath
1117         (
1118             time().caseName()/fileName(word("processor") + Foam::name(procI))
1119         );
1120
1121         // create a database
1122         Time processorDb
1123         (
1124             Time::controlDictName,
1125             time().rootPath(),
1126             processorCasePath
1127         );
1128
1129         // read finite volume mesh
1130         fvMesh procFvMesh
1131         (
1132             IOobject
1133             (
1134                 fvMesh::defaultRegion,
1135                 processorDb.timeName(),
1136                 processorDb
1137             )
1138         );
1139
1140         labelIOList fvBoundaryProcAddressing
1141         (
1142             IOobject
1143             (
1144                 "boundaryProcAddressing",
1145                 "constant",
1146                 procFvMesh.meshSubDir,
1147                 procFvMesh,
1148                 IOobject::MUST_READ,
1149                 IOobject::NO_WRITE
1150             )
1151         );
1152
1153
1154         // create finite area mesh
1155         faMesh procMesh
1156         (
1157             procFvMesh,
1158             procFaceLabels_[procI]
1159         );
1160
1161         // Create processor boundary patches
1162         const labelList& curBoundaryAddressing =
1163             procBoundaryAddressing_[procI];
1164
1165         const labelList& curPatchSizes = procPatchSize_[procI];
1166
1167         const labelList& curNeighbourProcessors =
1168             procNeighbourProcessors_[procI];
1169
1170         const labelList& curProcessorPatchSizes =
1171             procProcessorPatchSize_[procI];
1172
1173         const labelListList& curPatchEdgeLabels =
1174             procPatchEdgeLabels_[procI];
1175
1176         const faPatchList& meshPatches = boundary();
1177
1178         List<faPatch*> procPatches
1179         (
1180             curPatchSizes.size()
1181           + curProcessorPatchSizes.size(),
1182             reinterpret_cast<faPatch*>(NULL)
1183         );
1184
1185         label nPatches = 0;
1186
1187         forAll (curPatchSizes, patchi)
1188         {
1189             const labelList& curEdgeLabels = curPatchEdgeLabels[nPatches];
1190
1191             label ngbPolyPatchIndex =
1192                 findIndex
1193                 (
1194                     fvBoundaryProcAddressing,
1195                     meshPatches[curBoundaryAddressing[patchi]]
1196                    .ngbPolyPatchIndex()
1197                 );
1198
1199             procPatches[nPatches] =
1200                 meshPatches[curBoundaryAddressing[patchi]].clone
1201                 (
1202                     procMesh.boundary(),
1203                     curEdgeLabels,
1204                     nPatches,
1205                     ngbPolyPatchIndex
1206                 ).ptr();
1207
1208             nPatches++;
1209         }
1210
1211         forAll (curProcessorPatchSizes, procPatchI)
1212         {
1213             const labelList& curEdgeLabels = curPatchEdgeLabels[nPatches];
1214
1215             procPatches[nPatches] =
1216                 new processorFaPatch
1217                 (
1218                     word("procBoundary") + Foam::name(procI)
1219                   + word("to")
1220                   + Foam::name(curNeighbourProcessors[procPatchI]),
1221                     curEdgeLabels,
1222                     nPatches,
1223                     procMesh.boundary(),
1224                     -1,
1225                     procI,
1226                     curNeighbourProcessors[procPatchI]
1227                 );
1228
1229             nPatches++;
1230         }
1231
1232         // Add boundary patches
1233         procMesh.addFaPatches(procPatches);
1234
1235         // Set the precision of the points data to 10
1236         IOstream::defaultPrecision(10);
1237
1238         procMesh.write();
1239
1240         Info<< endl
1241             << "Processor " << procI << nl
1242             << "    Number of faces = " << procMesh.nFaces()
1243             << endl;
1244
1245         label nBoundaryEdges = 0;
1246         label nProcPatches = 0;
1247         label nProcEdges = 0;
1248
1249         forAll (procMesh.boundary(), patchi)
1250         {
1251             if
1252             (
1253                 procMesh.boundary()[patchi].type()
1254              == processorFaPatch::typeName
1255             )
1256             {
1257                 const processorFaPatch& ppp =
1258                 refCast<const processorFaPatch>
1259                 (
1260                     procMesh.boundary()[patchi]
1261                 );
1262
1263                 Info<< "    Number of edges shared with processor "
1264                     << ppp.neighbProcNo() << " = " << ppp.size() << endl;
1265
1266                 nProcPatches++;
1267                 nProcEdges += ppp.size();
1268             }
1269             else
1270             {
1271                 nBoundaryEdges += procMesh.boundary()[patchi].size();
1272             }
1273         }
1274
1275         Info<< "    Number of processor patches = " << nProcPatches << nl
1276             << "    Number of processor edges = " << nProcEdges << nl
1277             << "    Number of boundary edges = " << nBoundaryEdges << endl;
1278
1279         totProcEdges += nProcEdges;
1280         maxProcPatches = max(maxProcPatches, nProcPatches);
1281         maxProcEdges = max(maxProcEdges, nProcEdges);
1282
1283         // create and write the addressing information
1284         labelIOList pointProcAddressing
1285         (
1286             IOobject
1287             (
1288                 "pointProcAddressing",
1289                 "constant",
1290                 procMesh.meshSubDir,
1291                 procFvMesh,
1292                 IOobject::NO_READ,
1293                 IOobject::NO_WRITE
1294             ),
1295             procPatchPointAddressing_[procI]
1296         );
1297         pointProcAddressing.write();
1298
1299         labelIOList edgeProcAddressing
1300         (
1301             IOobject
1302             (
1303                 "edgeProcAddressing",
1304                 "constant",
1305                 procMesh.meshSubDir,
1306                 procFvMesh,
1307                 IOobject::NO_READ,
1308                 IOobject::NO_WRITE
1309             ),
1310             procEdgeAddressing_[procI]
1311         );
1312         edgeProcAddressing.write();
1313
1314         labelIOList faceProcAddressing
1315         (
1316             IOobject
1317             (
1318                 "faceProcAddressing",
1319                 "constant",
1320                 procMesh.meshSubDir,
1321                 procFvMesh,
1322                 IOobject::NO_READ,
1323                 IOobject::NO_WRITE
1324             ),
1325             procFaceAddressing_[procI]
1326         );
1327         faceProcAddressing.write();
1328
1329         labelIOList boundaryProcAddressing
1330         (
1331             IOobject
1332             (
1333                 "boundaryProcAddressing",
1334                 "constant",
1335                 procMesh.meshSubDir,
1336                 procFvMesh,
1337                 IOobject::NO_READ,
1338                 IOobject::NO_WRITE
1339             ),
1340             procBoundaryAddressing_[procI]
1341         );
1342         boundaryProcAddressing.write();
1343     }
1344
1345     Info<< nl
1346         << "Number of processor edges = " << totProcEdges/2 << nl
1347         << "Max number of processor patches = " << maxProcPatches << nl
1348         << "Max number of faces between processors = " << maxProcEdges
1349         << endl;
1350
1351     return true;
1352 }