applications/utilities/parallelProcessing/reconstructParMesh/processorMeshesRebuild.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 "processorMeshesReconstructor.H"
  27 #include "processorPolyPatch.H"
  28 #include "globalMeshData.H"
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 bool Foam::processorMeshesReconstructor::readMapping()
  33 {
  34     // Check for mapping
  35     Info<< "Check for mesh mapping data for instance "
  36         << meshes_[0].facesInstance() << ".  ";
  37
  38     bool readOk = true;
  39
  40     forAll (meshes_, procI)
  41     {
  42         const fvMesh& procMesh = meshes_[procI];
  43
  44         IOobject pointProcAddressingHeader
  45         (
  46             "pointProcAddressing",
  47             procMesh.facesInstance(),
  48             procMesh.meshSubDir,
  49             procMesh,
  50             IOobject::MUST_READ
  51         );
  52
  53         IOobject faceProcAddressingHeader
  54         (
  55             "faceProcAddressing",
  56             procMesh.facesInstance(),
  57             procMesh.meshSubDir,
  58             procMesh,
  59             IOobject::NO_READ,
  60             IOobject::NO_WRITE
  61         );
  62
  63         IOobject cellProcAddressingHeader
  64         (
  65             "cellProcAddressing",
  66             procMesh.facesInstance(),
  67             procMesh.meshSubDir,
  68             procMesh,
  69             IOobject::NO_READ,
  70             IOobject::NO_WRITE
  71         );
  72
  73         IOobject boundaryProcAddressingHeader
  74         (
  75             "boundaryProcAddressing",
  76             procMesh.facesInstance(),
  77             procMesh.meshSubDir,
  78             procMesh,
  79             IOobject::NO_READ,
  80             IOobject::NO_WRITE
  81         );
  82
  83         if
  84         (
  85             !pointProcAddressingHeader.headerOk()
  86          || !faceProcAddressingHeader.headerOk()
  87          || !cellProcAddressingHeader.headerOk()
  88          || !boundaryProcAddressingHeader.headerOk()
  89         )
  90         {
  91             readOk = false;
  92             break;
  93         }
  94     }
  95
  96     // All processors are fine: read mapping data
  97     if (readOk)
  98     {
  99         Info<< "Mapping data present.  Reading." << endl;
 100
 101         // Size the mapping arrays
 102         pointProcAddressing_.setSize(meshes_.size());
 103         faceProcAddressing_.setSize(meshes_.size());
 104         cellProcAddressing_.setSize(meshes_.size());
 105         boundaryProcAddressing_.setSize(meshes_.size());
 106
 107         forAll (meshes_, procI)
 108         {
 109             const fvMesh& procMesh = meshes_[procI];
 110
 111             pointProcAddressing_.set
 112             (
 113                 procI,
 114                 new labelIOList
 115                 (
 116                     IOobject
 117                     (
 118                         "pointProcAddressing",
 119                         procMesh.facesInstance(),
 120                         procMesh.meshSubDir,
 121                         procMesh,
 122                         IOobject::MUST_READ,
 123                         IOobject::NO_WRITE
 124                     )
 125                 )
 126             );
 127
 128             faceProcAddressing_.set
 129             (
 130                 procI,
 131                 new labelIOList
 132                 (
 133                     IOobject
 134                     (
 135                         "faceProcAddressing",
 136                         procMesh.facesInstance(),
 137                         procMesh.meshSubDir,
 138                         procMesh,
 139                         IOobject::MUST_READ,
 140                         IOobject::NO_WRITE
 141                     )
 142                 )
 143             );
 144
 145             cellProcAddressing_.set
 146             (
 147                 procI,
 148                 new labelIOList
 149                 (
 150                     IOobject
 151                     (
 152                         "cellProcAddressing",
 153                         procMesh.facesInstance(),
 154                         procMesh.meshSubDir,
 155                         procMesh,
 156                         IOobject::MUST_READ,
 157                         IOobject::NO_WRITE
 158                     )
 159                 )
 160             );
 161
 162             boundaryProcAddressing_.set
 163             (
 164                 procI,
 165                 new labelIOList
 166                 (
 167                     IOobject
 168                     (
 169                         "boundaryProcAddressing",
 170                         procMesh.facesInstance(),
 171                         procMesh.meshSubDir,
 172                         procMesh,
 173                         IOobject::MUST_READ,
 174                         IOobject::NO_WRITE
 175                     )
 176                 )
 177             );
 178         }
 179
 180         Info<< "Addressing from files: " << nl;
 181         forAll (meshes_, procI)
 182         {
 183             Info<< "Proc " << procI
 184                 << " point addr: " << pointProcAddressing_[procI].size()
 185                 << " face addr: " << faceProcAddressing_[procI].size()
 186                 << " cell addr: " << cellProcAddressing_[procI].size()
 187                 << " boundary addr: " << boundaryProcAddressing_[procI].size()
 188                 << endl;
 189         }
 190     }
 191     else
 192     {
 193         Info<< "No mapping data available." << endl;
 194     }
 195
 196     return readOk;
 197 }
 198
 199
 200 void Foam::processorMeshesReconstructor::writeAddressing()
 201 {
 202     forAll (pointProcAddressing_, procI)
 203     {
 204         pointProcAddressing_[procI].write();
 205         faceProcAddressing_[procI].write();
 206         cellProcAddressing_[procI].write();
 207         boundaryProcAddressing_[procI].write();
 208     }
 209 }
 210
 211
 212 const Foam::processorPolyPatch&
 213 Foam::processorMeshesReconstructor::neighbourProcPatch
 214 (
 215     const processorPolyPatch& procPatch
 216 ) const
 217 {
 218     const label masterProcID = procPatch.neighbProcNo();
 219
 220     const polyMesh& masterMesh = meshes_[masterProcID];
 221
 222     bool found = false;
 223
 224     // Find the processor patch that corresponds to current patch
 225     const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
 226
 227     forAll (masterPatches, masterPatchI)
 228     {
 229         if
 230         (
 231             isA<processorPolyPatch>
 232             (
 233                 masterPatches[masterPatchI]
 234             )
 235         )
 236         {
 237             const processorPolyPatch& masterProcPatch =
 238                 refCast<const processorPolyPatch>
 239                 (
 240                     masterPatches[masterPatchI]
 241                 );
 242
 243             // Check neighbour processor index
 244             if (masterProcPatch.neighbProcNo() == procPatch.myProcNo())
 245             {
 246                 // Found matching patch
 247                 return masterProcPatch;
 248             }
 249         }
 250     }
 251
 252     if (!found)
 253     {
 254         FatalErrorIn
 255         (
 256             "const processorPolyPatch&\n"
 257             "processorMeshesReconstructor::neighbourProcPatch\n"
 258             "(\n"
 259             "    const processorPolyPatch& procPatch\n"
 260             ") const"
 261         )   << "Cannot find processor patch pair ("
 262             << procPatch.myProcNo() << " "
 263             << procPatch.neighbProcNo() << ") for merging"
 264             << abort(FatalError);
 265     }
 266
 267     // Dummy return
 268     return procPatch;
 269 }
 270
 271
 272 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 273
 274 void Foam::processorMeshesReconstructor::reconstructPoints(fvMesh& mesh) const
 275 {
 276     // Create the new points
 277     vectorField newPoints(mesh.nPoints());
 278
 279     forAll (meshes_, procI)
 280     {
 281         // Reconstruct only live points.  HJ, 7/Mar/2011
 282         const vectorField& procPoints = meshes_[procI].points();
 283
 284         // Set the cell values in the reconstructed field
 285
 286         const labelList& pointProcAddressingI = pointProcAddressing()[procI];
 287
 288         if (pointProcAddressingI.size() != procPoints.size())
 289         {
 290             FatalErrorIn("processorMeshes")
 291                 << "problem :"
 292                 << " pointProcAddressingI:" << pointProcAddressingI.size()
 293                 << " procPoints:" << procPoints.size()
 294                 << abort(FatalError);
 295         }
 296
 297         // Only live points carry reconstruction data.  Reconsider
 298         // HJ, 6/Sep/2009
 299         for (label pointI = 0; pointI < meshes_[procI].nPoints(); pointI++)
 300         {
 301             newPoints[pointProcAddressingI[pointI]] = procPoints[pointI];
 302         }
 303     }
 304
 305     mesh.movePoints(newPoints);
 306     mesh.write();
 307 }
 308
 309
 310 Foam::autoPtr<Foam::fvMesh>
 311 Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 312 {
 313     // Check for read
 314     if (readMapping())
 315     {
 316         // Mapping data present and read.  Reconstructed mesh may be
 317         // present as well
 318         bool readMesh = false;
 319
 320         {
 321             IOobject points
 322             (
 323                 "points",
 324                 db.timeName(),
 325                 meshes_[0].meshSubDir,
 326                 db,
 327                 IOobject::MUST_READ
 328             );
 329
 330             if (points.headerOk())
 331             {
 332                 readMesh = true;
 333             }
 334         }
 335
 336         if (readMesh)
 337         {
 338             Info<< "Global mesh present for time " << db.timeName()
 339                 << ".  Reading mesh." << endl;
 340
 341             autoPtr<fvMesh> globalMeshPtr
 342             (
 343                 new fvMesh
 344                 (
 345                     IOobject
 346                     (
 347                         meshName_,
 348                         db.timeName(),
 349                         db,
 350                         IOobject::MUST_READ,
 351                         IOobject::NO_WRITE
 352                     )
 353                 )
 354             );
 355
 356             return globalMeshPtr;
 357         }
 358     }
 359
 360     // Prepare patch reconstruction from processor 0
 361     wordList reconPatchNames;
 362     wordList reconPatchTypes;
 363     label nReconPatches = 0;
 364
 365     // Grab names and types of patches, excluding processor patches
 366     // Order must be identical on all processors
 367     {
 368         const polyBoundaryMesh& proc0Boundary = meshes_[0].boundaryMesh();
 369
 370         reconPatchNames.setSize(proc0Boundary.size());
 371         reconPatchTypes.setSize(proc0Boundary.size());
 372         nReconPatches = 0;
 373
 374         forAll (proc0Boundary, patchI)
 375         {
 376             if (!isA<processorPolyPatch>(proc0Boundary[patchI]))
 377             {
 378                 // Add patch name to list of patches
 379                 reconPatchNames[nReconPatches] = proc0Boundary[patchI].name();
 380                 reconPatchTypes[nReconPatches] = proc0Boundary[patchI].type();
 381                 nReconPatches++;
 382             }
 383         }
 384
 385         reconPatchNames.setSize(nReconPatches);
 386         reconPatchTypes.setSize(nReconPatches);
 387     }
 388
 389     // Prepare point, face and patch reconstruction
 390     label nReconPoints = 0;
 391     label nReconFaces = 0;
 392     label nReconCells = 0;
 393     labelList reconPatchSizes(reconPatchTypes.size(), 0);
 394
 395     forAll (meshes_, procI)
 396     {
 397         // Count total number of points and faces
 398         nReconPoints += meshes_[procI].nPoints();
 399         nReconFaces += meshes_[procI].allFaces().size();
 400         nReconCells += meshes_[procI].nCells();
 401
 402         const polyBoundaryMesh& procPatches = meshes_[procI].boundaryMesh();
 403
 404         nReconPatches = 0;
 405
 406         forAll (procPatches, patchI)
 407         {
 408             if (!isA<processorPolyPatch>(procPatches[patchI]))
 409             {
 410                 // Check patch name and type
 411                 if
 412                 (
 413                     procPatches[patchI].name()
 414                  != reconPatchNames[nReconPatches]
 415                  || procPatches[patchI].type()
 416                  != reconPatchTypes[nReconPatches]
 417                 )
 418                 {
 419                     FatalErrorIn
 420                     (
 421                         "autoPtr<fvMesh> processorMeshesReconstructor::"
 422                         "reconstructMesh(const Time& db)"
 423                     )   << "Patch names, types or ordering does not match "
 424                         << "across processors"
 425                         << abort(FatalError);
 426                 }
 427
 428                 // Record number of faces in patch
 429                 reconPatchSizes[nReconPatches] += procPatches[patchI].size();
 430                 nReconPatches++;
 431             }
 432         }
 433     }
 434
 435     Info<< "Estimated max global mesh size (with duplicates): " << nl
 436         << "    nPoints = " << nReconPoints << nl
 437         << "    nFaces = " << nReconFaces << nl
 438         << "    nCells = " << nReconCells << nl
 439         << "    nPatches = " << nReconPatches << nl
 440         << "    nPatchFaces = " << reconPatchSizes << endl;
 441
 442     // Note: for easier debugging, set owner and neighbour to -1
 443     pointField reconPoints(nReconPoints);
 444     faceList reconFaces(nReconFaces);
 445     labelList cellOffset(meshes_.size(), 0);
 446     labelList reconOwner(nReconFaces, -1);
 447     labelList reconNeighbour(nReconFaces, -1);
 448     faceListList reconPatchFaces(reconPatchTypes.size());
 449     labelListList reconPatchOwner(reconPatchTypes.size());
 450
 451     forAll (reconPatchFaces, patchI)
 452     {
 453         reconPatchFaces[patchI].setSize(reconPatchSizes[patchI]);
 454
 455         reconPatchOwner[patchI].setSize(reconPatchSizes[patchI]);
 456         reconPatchOwner[patchI] = -1;
 457     }
 458
 459     // Size the mapping arrays
 460     pointProcAddressing_.setSize(meshes_.size());
 461     faceProcAddressing_.setSize(meshes_.size());
 462     cellProcAddressing_.setSize(meshes_.size());
 463     boundaryProcAddressing_.setSize(meshes_.size());
 464
 465     // Allocate addressing arrays on all meshes
 466     forAll (meshes_, procI)
 467     {
 468         const fvMesh& procMesh = meshes_[procI];
 469
 470         pointProcAddressing_.set
 471         (
 472             procI,
 473             new labelIOList
 474             (
 475                 IOobject
 476                 (
 477                     "pointProcAddressing",
 478                     procMesh.facesInstance(),
 479                     procMesh.meshSubDir,
 480                     procMesh,
 481                     IOobject::NO_READ,
 482                     IOobject::NO_WRITE
 483                 ),
 484                 labelList(procMesh.nPoints(), -1)
 485             )
 486         );
 487
 488         faceProcAddressing_.set
 489         (
 490             procI,
 491             new labelIOList
 492             (
 493                 IOobject
 494                 (
 495                     "faceProcAddressing",
 496                     procMesh.facesInstance(),
 497                     procMesh.meshSubDir,
 498                     procMesh,
 499                     IOobject::NO_READ,
 500                     IOobject::NO_WRITE
 501                 ),
 502                 labelList(procMesh.allFaces().size(), -1)
 503             )
 504         );
 505
 506         cellProcAddressing_.set
 507         (
 508             procI,
 509             new labelIOList
 510             (
 511                 IOobject
 512                 (
 513                     "cellProcAddressing",
 514                     procMesh.facesInstance(),
 515                     procMesh.meshSubDir,
 516                     procMesh,
 517                     IOobject::NO_READ,
 518                     IOobject::NO_WRITE
 519                 ),
 520                 labelList(procMesh.nCells(), -1)
 521             )
 522         );
 523
 524         boundaryProcAddressing_.set
 525         (
 526             procI,
 527             new labelIOList
 528             (
 529                 IOobject
 530                 (
 531                     "boundaryProcAddressing",
 532                     procMesh.facesInstance(),
 533                     procMesh.meshSubDir,
 534                     procMesh,
 535                     IOobject::NO_READ,
 536                     IOobject::NO_WRITE
 537                 ),
 538                 labelList(procMesh.boundaryMesh().size(), -1)
 539             )
 540         );
 541     }
 542
 543     // Reset the counters
 544     nReconPoints = 0;
 545     nReconFaces = 0;
 546
 547     reconPatchSizes = 0;
 548
 549     // Prepare handling for globally shared points
 550     labelList globalPointMapping;
 551
 552     // Dump mesh zero without checking
 553     {
 554         cellOffset[0] = 0;
 555
 556         const polyMesh& curMesh = meshes_[0];
 557         labelList& ppAddr = pointProcAddressing_[0];
 558         labelList& fpAddr = faceProcAddressing_[0];
 559         labelList& cpAddr = cellProcAddressing_[0];
 560         labelList& bpAddr = boundaryProcAddressing_[0];
 561
 562         // Prepare handling for global mesh data: set to -1
 563         globalPointMapping.setSize(curMesh.globalData().nGlobalPoints());
 564         globalPointMapping = -1;
 565
 566         // Dump all points into the global point list
 567         // Reconstruct only live points.  HJ, 7/Mar/2011
 568         const pointField& curPoints = curMesh.points();
 569         ppAddr.setSize(curPoints.size());
 570
 571         forAll (curPoints, pointI)
 572         {
 573             reconPoints[nReconPoints] = curPoints[pointI];
 574             ppAddr[pointI] = nReconPoints;
 575             nReconPoints++;
 576         }
 577
 578         // Collect globally shared point labels
 579         const labelList& curSpl = curMesh.globalData().sharedPointLabels();
 580
 581         forAll (curSpl, splI)
 582         {
 583             // From processor 0, mark points without checking
 584             globalPointMapping[curSpl[splI]] = ppAddr[curSpl[splI]];
 585         }
 586
 587         // Dump all internal faces into the list
 588         const faceList& curFaces = curMesh.allFaces();
 589         const labelList& curOwner = curMesh.faceOwner();
 590         const labelList& curNeighbour = curMesh.faceNeighbour();
 591         fpAddr.setSize(curFaces.size());
 592
 593         for (label faceI = 0; faceI < curMesh.nInternalFaces(); faceI++)
 594         {
 595             // Renumber face in new vertices
 596             face newFace = curFaces[faceI];
 597             inplaceRenumber(ppAddr, newFace);
 598
 599             reconFaces[nReconFaces] = newFace;
 600             reconOwner[nReconFaces] = curOwner[faceI];//+cellOffset[0];
 601             reconNeighbour[nReconFaces] = curNeighbour[faceI];//+cellOffset[0];
 602
 603             // Face-processor addressing uses offset of 1 and a turning index
 604             // If the label is negative, it means the global face points
 605             // in direction opposite to decomposed face.  HJ, 16/Feb/2011
 606             fpAddr[faceI] = nReconFaces + 1;
 607             nReconFaces++;
 608         }
 609
 610         // Go through all patches.  For regular patches
 611         // dump the faces into patch lists
 612         const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();
 613
 614         forAll (procPatches, patchI)
 615         {
 616             if (isA<processorPolyPatch>(procPatches[patchI]))
 617             {
 618                 // Processor patch: faces become internal faces
 619                 const processorPolyPatch& procPatch =
 620                     refCast<const processorPolyPatch>(procPatches[patchI]);
 621
 622                 for
 623                 (
 624                     label faceI = procPatch.start();
 625                     faceI < procPatch.start() + procPatch.size();
 626                     faceI++
 627                 )
 628                 {
 629                     // Renumber face in new vertices
 630                     face newFace = curFaces[faceI];
 631                     inplaceRenumber(ppAddr, newFace);
 632
 633                     reconFaces[nReconFaces] = newFace;
 634
 635                     fpAddr[faceI] = nReconFaces + 1;
 636                     reconOwner[nReconFaces] = curOwner[faceI];//+cellOffset[0];
 637
 638                     // For partially completed neighbour, set nbr to -2
 639                     // for easier debugging
 640                     reconNeighbour[nReconFaces] = -2;
 641
 642                     nReconFaces++;
 643                 }
 644             }
 645             else
 646             {
 647                 // Regular patch: dump faces into patch face list
 648                 faceList& curRpFaces = reconPatchFaces[patchI];
 649                 labelList& curRpfOwner = reconPatchOwner[patchI];
 650                 label& nRpf = reconPatchSizes[patchI];
 651
 652                 const polyPatch& curPatch = procPatches[patchI];
 653
 654                 for
 655                 (
 656                     label faceI = curPatch.start();
 657                     faceI < curPatch.start() + curPatch.size();
 658                     faceI++
 659                 )
 660                 {
 661                     // Renumber face in new vertices
 662                     face newFace = curFaces[faceI];
 663                     inplaceRenumber(ppAddr, newFace);
 664
 665                     curRpFaces[nRpf] = newFace;
 666                     curRpfOwner[nRpf] = curOwner[faceI];//+cellOffset[0];
 667
 668                     // Temporarily record position of face in the patch.
 669                     // Offset for nInternalFaces will be added in the end
 670                     // when the complete list of faces is assembled
 671                     // HJ, 16/Feb/2011
 672                     fpAddr[faceI] = nRpf + 1;
 673                     nRpf++;
 674
 675                 }
 676             }
 677         }
 678
 679         // Cell-processor addressing
 680         forAll (cpAddr, cellI)
 681         {
 682             cpAddr[cellI] = cellI; // + cellOffset[0];
 683         }
 684
 685         // Sort out boundary addressing: i for live patches, -1 for processor
 686         bpAddr = -1;
 687
 688         // Note: loop over mapped patches
 689         forAll (reconPatchSizes, patchI)
 690         {
 691             bpAddr[patchI] = patchI;
 692         }
 693     }
 694
 695
 696     // Dump all other meshes, merging the processor boundaries
 697
 698     for (label procI = 1; procI < meshes_.size(); procI++)
 699     {
 700         // Grab cell offset from previous offset and mesh size
 701         cellOffset[procI] =
 702             cellOffset[procI - 1] + meshes_[procI - 1].nCells();
 703
 704         const polyMesh& curMesh = meshes_[procI];
 705         const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();
 706
 707         labelList& ppAddr = pointProcAddressing_[procI];
 708         labelList& fpAddr = faceProcAddressing_[procI];
 709         labelList& cpAddr = cellProcAddressing_[procI];
 710         labelList& bpAddr = boundaryProcAddressing_[procI];
 711
 712         // Point mapping
 713
 714         // Reconstruct only live points.  HJ, 7/Mar/2011
 715         const pointField& curPoints = curMesh.points();
 716
 717         // Set ppAddr to -1, to use as point usage indicators
 718         ppAddr.setSize(curPoints.size());
 719         ppAddr = -1;
 720
 721         // Find points already added via processor patches and mark them
 722         // in ppAddr
 723
 724         // Collect point-processor addressing for points on processor patches
 725
 726         // Go through all patches.  For neighbour patches, access
 727         // owner addressing and dump into ppAddr
 728         forAll (procPatches, patchI)
 729         {
 730             if (isA<processorPolyPatch>(procPatches[patchI]))
 731             {
 732                 // Processor patch: faces become internal faces
 733                 const processorPolyPatch& procPatch =
 734                     refCast<const processorPolyPatch>(procPatches[patchI]);
 735
 736                 // If patch is a neighbour, its master has already inserted
 737                 // the points
 738                 if (procPatch.neighbour())
 739                 {
 740                     const label masterProcID = procPatch.neighbProcNo();
 741
 742                     // Get the neighbour side patch
 743                     const processorPolyPatch& masterProcPatch =
 744                         neighbourProcPatch(procPatch);
 745
 746                     // Find the addressing of the master side
 747                     const labelList& masterPpAddr =
 748                         pointProcAddressing_[masterProcID];
 749
 750                     // Assemble neighbour mesh point addressing in matching
 751                     // order by reversing processor patch faces
 752                     faceList reversedFaces(procPatch.size());
 753
 754                     forAll (reversedFaces, faceI)
 755                     {
 756                         reversedFaces[faceI] = procPatch[faceI].reverseFace();
 757                     }
 758
 759                     primitiveFacePatch reversedPatch
 760                     (
 761                         reversedFaces,
 762                         procPatch.points()
 763                     );
 764
 765                     // Insert addressing from master side into
 766                     // local point addressing.  Each face of reversed patch
 767                     // now matches the master face.
 768                     // Note: this is done by visiting faces, since meshPoints
 769                     // are ordered in increasing order.  HJ, 10/Mar/2011
 770
 771                     forAll (reversedFaces, faceI)
 772                     {
 773                         // Current reverse face
 774                         const face& curRF = reversedFaces[faceI];
 775
 776                         // Current master face
 777                         const face& curMF = masterProcPatch[faceI];
 778
 779                         forAll (curRF, pointI)
 780                         {
 781                             // Mapping is established
 782                             ppAddr[curRF[pointI]] =
 783                                 masterPpAddr[curMF[pointI]];
 784                         }
 785                     }
 786                 } // End of "is neighbour"
 787             } // End of "is processor"
 788         }
 789
 790         // Dump unmarked points into the global point list
 791         label nMergedPoints = 0;
 792
 793         forAll (curPoints, pointI)
 794         {
 795             if (ppAddr[pointI] == -1)
 796             {
 797                 // Unmerged point
 798                 reconPoints[nReconPoints] = curPoints[pointI];
 799                 ppAddr[pointI] = nReconPoints;
 800                 nReconPoints++;
 801             }
 802             else
 803             {
 804                 nMergedPoints++;
 805             }
 806         }
 807
 808         Info<< "Processor " << procI << " merged " << nMergedPoints
 809             << " points out of local " << curPoints.size()
 810             << " and total " << nReconPoints << endl;
 811
 812         // Dump all internal faces into the list
 813         const faceList& curFaces = curMesh.allFaces();
 814         const labelList& curOwner = curMesh.faceOwner();
 815         const labelList& curNeighbour = curMesh.faceNeighbour();
 816         fpAddr.setSize(curFaces.size());
 817
 818         // Collect globally shared point labels
 819         const labelList& curSpl = curMesh.globalData().sharedPointLabels();
 820
 821         forAll (curSpl, splI)
 822         {
 823             // From other processors, check if point is already marked
 824             // If not, mark it; otherwise compare (and correct?) with local
 825             // mark
 826             if (globalPointMapping[curSpl[splI]] < 0)
 827             {
 828                 globalPointMapping[curSpl[splI]] = ppAddr[curSpl[splI]];
 829             }
 830             else
 831             {
 832                 // Compare.  Is this needed - should always be OK.
 833                 if (globalPointMapping[curSpl[splI]] != ppAddr[curSpl[splI]])
 834                 {
 835                     WarningIn
 836                     (
 837                         "autoPtr<fvMesh> "
 838                         "processorMeshesReconstructor::"
 839                         "reconstructMesh(const Time& db)"
 840                     )   << "Loss of sync???"
 841                         << abort(FatalError);
 842                 }
 843             }
 844         }
 845
 846
 847         for (label faceI = 0; faceI < curMesh.nInternalFaces(); faceI++)
 848         {
 849             // Renumber face in new vertices
 850             face newFace = curFaces[faceI];
 851             inplaceRenumber(ppAddr, newFace);
 852
 853             reconFaces[nReconFaces] = newFace;
 854             reconOwner[nReconFaces] = curOwner[faceI] + cellOffset[procI];
 855             reconNeighbour[nReconFaces] = curNeighbour[faceI]
 856                 + cellOffset[procI];
 857             fpAddr[faceI] = nReconFaces + 1;
 858             nReconFaces++;
 859         }
 860
 861         // Go through all patches.  For regular patches
 862         // dump the faces into patch lists
 863         forAll (procPatches, patchI)
 864         {
 865             if (isA<processorPolyPatch>(procPatches[patchI]))
 866             {
 867                 // Processor patch: faces become internal faces
 868                 const processorPolyPatch& procPatch =
 869                     refCast<const processorPolyPatch>(procPatches[patchI]);
 870
 871                 // If patch is a master, drop the faces and fill the
 872                 // owner side addressing
 873                 if (procPatch.owner())
 874                 {
 875                     for
 876                     (
 877                         label faceI = procPatch.start();
 878                         faceI < procPatch.start() + procPatch.size();
 879                         faceI++
 880                     )
 881                     {
 882                         // Renumber face in new vertices
 883                         face newFace = curFaces[faceI];
 884                         inplaceRenumber(ppAddr, newFace);
 885
 886                         reconFaces[nReconFaces] = newFace;
 887                         fpAddr[faceI] = nReconFaces + 1;
 888                         reconOwner[nReconFaces] = curOwner[faceI]
 889                             + cellOffset[procI];
 890
 891                         // For partially completed neighbour, set nbr to -2
 892                         // for easier debugging
 893                         reconNeighbour[nReconFaces] = -2;
 894
 895                         nReconFaces++;
 896                     }
 897                 }
 898                 else
 899                 {
 900                     // Fill the addressing for the neighbour side
 901
 902                     const label masterProcID = procPatch.neighbProcNo();
 903
 904                     // Get local face-cell addressing: it will become neighbour
 905                     // addressing for the already inserted faces
 906                     const labelList procFaceCells = procPatch.faceCells();
 907
 908                     // Get the neighbour side patch
 909                     const processorPolyPatch& masterProcPatch =
 910                         neighbourProcPatch(procPatch);
 911
 912                     // Find the addressing of the master side
 913                     // and insert the neighbour with offset
 914
 915                     const labelList& masterFaceAddr =
 916                         faceProcAddressing_[masterProcID];
 917
 918                     for
 919                     (
 920                         label faceI = procPatch.start();
 921                         faceI < procPatch.start() + procPatch.size();
 922                         faceI++
 923                     )
 924                     {
 925                         label faceInPatch = faceI - procPatch.start();
 926
 927                         // Calculate master index
 928                         label masterIndex = masterProcPatch.start()
 929                             + faceInPatch;
 930
 931                         label masterFp = masterFaceAddr[masterIndex] - 1;
 932
 933                         // Record face-cells for the neighbour
 934                         fpAddr[faceI] = -masterFaceAddr[masterIndex];
 935
 936                         reconNeighbour[masterFp] =
 937                             procFaceCells[faceInPatch] + cellOffset[procI];
 938                     }
 939                 }
 940             }
 941             else
 942             {
 943                 // Regular patch: dump faces into patch face list
 944                 faceList& curRpFaces = reconPatchFaces[patchI];
 945                 labelList& curRpfOwner = reconPatchOwner[patchI];
 946                 label& nRpf = reconPatchSizes[patchI];
 947
 948                 const polyPatch& curPatch = procPatches[patchI];
 949
 950                 for
 951                 (
 952                     label faceI = curPatch.start();
 953                     faceI < curPatch.start() + curPatch.size();
 954                     faceI++
 955                 )
 956                 {
 957                     // Renumber face in new vertices
 958                     face newFace = curFaces[faceI];
 959                     inplaceRenumber(ppAddr, newFace);
 960
 961                     curRpFaces[nRpf] = newFace;
 962                     curRpfOwner[nRpf] = curOwner[faceI] + cellOffset[procI];
 963
 964                     // Temporarily record position of face in the patch.
 965                     // Offset for nInternalFaces will be added in the end
 966                     // when the complete list of faces is assembled
 967                     // HJ, 16/Feb/2011
 968                     fpAddr[faceI] = nRpf + 1;
 969                     nRpf++;
 970
 971                 }
 972             }
 973         }
 974
 975         // Cell-processor addressing
 976         forAll (cpAddr, cellI)
 977         {
 978             cpAddr[cellI] = cellI + cellOffset[procI];
 979         }
 980
 981         // Sort out boundary addressing: i for live patches, -1 for processor
 982         bpAddr = -1;
 983
 984         // Note: loop over mapped patches
 985         forAll (reconPatchSizes, patchI)
 986         {
 987             bpAddr[patchI] = patchI;
 988         }
 989     }
 990
 991     // Resize the lists
 992     reconPoints.setSize(nReconPoints);
 993
 994     // Resize the neighbour list to the size of internalFaces
 995     label nInternalFaces = nReconFaces;
 996
 997     reconNeighbour.setSize(nInternalFaces);
 998
 999     // Resize the patch lists
1000     forAll (reconPatchFaces, patchI)
1001     {
1002         reconPatchFaces[patchI].setSize(reconPatchSizes[patchI]);
1003         reconPatchOwner[patchI].setSize(reconPatchSizes[patchI]);
1004     }
1005
1006     // Complete the global list of faces
1007     labelList reconPatchStarts(reconPatchSizes, 0);
1008
1009     // Copy the faces into face list
1010     forAll (reconPatchFaces, patchI)
1011     {
1012         reconPatchStarts[patchI] = nReconFaces;
1013
1014         const faceList& curPatchFaces = reconPatchFaces[patchI];
1015         const labelList& curPatchOwner = reconPatchOwner[patchI];
1016
1017         forAll (curPatchFaces, fI)
1018         {
1019             reconFaces[nReconFaces] = curPatchFaces[fI];
1020             reconOwner[nReconFaces] = curPatchOwner[fI];
1021             nReconFaces++;
1022         }
1023     }
1024
1025     reconFaces.setSize(nReconFaces);
1026     reconOwner.setSize(nReconFaces);
1027
1028     // Mesh assembly completed
1029
1030     Info<< "Global mesh size (final): " << nl
1031         << "    nPoints = " << reconPoints.size() << nl
1032         << "    nFaces = " << reconFaces.size() << nl
1033         << "    nCells = " << nReconCells << nl
1034         << "    nPatches = " << reconPatchSizes.size() << nl
1035         << "    nPatchFaces = " << reconPatchSizes << endl;
1036
1037     // Renumber the face-processor addressing list for all pieces
1038     // now that the number of internal faces is known
1039     forAll (meshes_, procI)
1040     {
1041         // Get processor mesh and boundary
1042         const polyMesh& curMesh = meshes_[procI];
1043         const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();
1044
1045         // Get face-processor addressing for corrent prorcessor
1046         labelList& fpAddr = faceProcAddressing_[procI];
1047
1048         const labelList& bpAddr = boundaryProcAddressing_[procI];
1049
1050         forAll (procPatches, patchI)
1051         {
1052             if (!isA<processorPolyPatch>(procPatches[patchI]))
1053             {
1054                 // Get master processor patch
1055                 const label reconPatchID = bpAddr[patchI];
1056
1057                 // Skip processor patches: bpAddr = -1
1058                 if (reconPatchID > -1)
1059                 {
1060                     const label reconStart = reconPatchStarts[reconPatchID];
1061
1062                     const polyPatch& curPatch = procPatches[patchI];
1063
1064                     for
1065                     (
1066                         label faceI = curPatch.start();
1067                         faceI < curPatch.start() + curPatch.size();
1068                         faceI++
1069                     )
1070                     {
1071                         // Add patch start
1072                         fpAddr[faceI] += reconStart;
1073                     }
1074                 }
1075             }
1076         }
1077     }
1078
1079     autoPtr<fvMesh> globalMeshPtr
1080     (
1081         new fvMesh
1082         (
1083             IOobject
1084             (
1085                 meshName_,
1086                 db.timeName(),
1087                 db,
1088                 IOobject::NO_READ
1089             ),
1090             xferCopy(reconPoints),
1091             xferCopy(reconFaces),
1092             xferCopy(reconOwner),
1093             xferCopy(reconNeighbour)
1094         )
1095     );
1096     fvMesh& globalMesh = globalMeshPtr();
1097
1098     // Create patch list using mesh from processor 0
1099     List<polyPatch*> reconPatches(nReconPatches);
1100
1101     {
1102         const polyBoundaryMesh& procPatches = meshes_[0].boundaryMesh();
1103
1104         reconPatches.setSize(reconPatchSizes.size());
1105
1106         forAll (reconPatchSizes, patchI)
1107         {
1108             reconPatches[patchI] =
1109                 procPatches[patchI].clone
1110                 (
1111                     globalMesh.boundaryMesh(),
1112                     patchI,
1113                     reconPatchSizes[patchI],
1114                     reconPatchStarts[patchI]
1115                 ).ptr();
1116         }
1117     }
1118
1119     // Add both poly and fv boundary patches
1120     globalMesh.addFvPatches(reconPatches);
1121
1122     // TODO: point, face and cell zones
1123
1124     Info<< "Reconstructed addressing: " << nl;
1125     forAll (meshes_, procI)
1126     {
1127         Info<< "Proc " << procI
1128             << " point addr: " << pointProcAddressing_[procI].size()
1129             << " face addr: " << faceProcAddressing_[procI].size()
1130             << " cell addr: " << cellProcAddressing_[procI].size()
1131             << " boundary addr: " << boundaryProcAddressing_[procI].size()
1132             << endl;
1133     }
1134
1135     return globalMeshPtr;
1136 }
1137
1138
1139
1140 // ************************************************************************* //