src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Post-processing mesh subset tool.  Given the original mesh and the
  27     list of selected cells, it creates the mesh consisting only of the
  28     desired cells, with the mapping list for points, faces, and cells.
  29
  30     MJ 23/03/05 on coupled faces change the patch of the face to the
  31     oldInternalFaces patch.
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "fvMeshSubset.H"
  36 #include "boolList.H"
  37 #include "Pstream.H"
  38 #include "emptyPolyPatch.H"
  39 #include "demandDrivenData.H"
  40 #include "cyclicPolyPatch.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 bool Foam::fvMeshSubset::checkCellSubset() const
  50 {
  51     if (fvMeshSubsetPtr_)
  52     {
  53         FatalErrorIn("bool fvMeshSubset::checkCellSubset() const")
  54             << "Mesh subset not set.  Please set the cell map using "
  55             << "void setCellSubset(const labelHashSet& cellsToSubset)" << endl
  56             << "before attempting to access subset data"
  57             << abort(FatalError);
  58
  59         return false;
  60     }
  61     else
  62     {
  63         return true;
  64     }
  65 }
  66
  67
  68 void Foam::fvMeshSubset::markPoints
  69 (
  70     const labelList& curPoints,
  71     Map<label>& pointMap
  72 )
  73 {
  74     forAll (curPoints, pointI)
  75     {
  76         // Note: insert will only insert if not yet there.
  77         pointMap.insert(curPoints[pointI], 0);
  78     }
  79 }
  80
  81
  82 void Foam::fvMeshSubset::markPoints
  83 (
  84     const labelList& curPoints,
  85     labelList& pointMap
  86 )
  87 {
  88     forAll (curPoints, pointI)
  89     {
  90         pointMap[curPoints[pointI]] = 0;
  91     }
  92 }
  93
  94
  95 // Synchronize nCellsUsingFace on both sides of coupled patches. Marks
  96 // faces that become 'uncoupled' with 3.
  97 void Foam::fvMeshSubset::doCoupledPatches
  98 (
  99     const bool syncPar,
 100     labelList& nCellsUsingFace
 101 ) const
 102 {
 103     const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
 104
 105     label nUncoupled = 0;
 106
 107     if (syncPar && Pstream::parRun())
 108     {
 109         // Send face usage across processor patches
 110         forAll (oldPatches, oldPatchI)
 111         {
 112             const polyPatch& pp = oldPatches[oldPatchI];
 113
 114             if (isA<processorPolyPatch>(pp))
 115             {
 116                 const processorPolyPatch& procPatch =
 117                     refCast<const processorPolyPatch>(pp);
 118
 119                 OPstream toNeighbour
 120                 (
 121                     Pstream::blocking,
 122                     procPatch.neighbProcNo()
 123                 );
 124
 125                 toNeighbour
 126                     << SubList<label>(nCellsUsingFace, pp.size(), pp.start());
 127             }
 128         }
 129
 130
 131         // Receive face usage count and check for faces that become uncoupled.
 132         forAll (oldPatches, oldPatchI)
 133         {
 134             const polyPatch& pp = oldPatches[oldPatchI];
 135
 136             if (isA<processorPolyPatch>(pp))
 137             {
 138                 const processorPolyPatch& procPatch =
 139                     refCast<const processorPolyPatch>(pp);
 140
 141                 IPstream fromNeighbour
 142                 (
 143                     Pstream::blocking,
 144                     procPatch.neighbProcNo()
 145                 );
 146
 147                 labelList nbrCellsUsingFace(fromNeighbour);
 148
 149                 // Combine with this side.
 150
 151                 forAll (pp, i)
 152                 {
 153                     if
 154                     (
 155                         nCellsUsingFace[pp.start()+i] == 1
 156                      && nbrCellsUsingFace[i] == 0
 157                     )
 158                     {
 159                         // Face's neighbour is no longer there. Mark face off
 160                         // as coupled
 161                         nCellsUsingFace[pp.start()+i] = 3;
 162                         nUncoupled++;
 163                     }
 164                 }
 165             }
 166         }
 167     }
 168
 169     // Do same for cyclics.
 170     forAll (oldPatches, oldPatchI)
 171     {
 172         const polyPatch& pp = oldPatches[oldPatchI];
 173
 174         if (isA<cyclicPolyPatch>(pp))
 175         {
 176             const cyclicPolyPatch& cycPatch =
 177                 refCast<const cyclicPolyPatch>(pp);
 178
 179             forAll (cycPatch, i)
 180             {
 181                 label thisFaceI = cycPatch.start() + i;
 182                 label otherFaceI = cycPatch.transformGlobalFace(thisFaceI);
 183
 184                 if
 185                 (
 186                     nCellsUsingFace[thisFaceI] == 1
 187                  && nCellsUsingFace[otherFaceI] == 0
 188                 )
 189                 {
 190                     nCellsUsingFace[thisFaceI] = 3;
 191                     nUncoupled++;
 192                 }
 193             }
 194         }
 195     }
 196
 197     if (syncPar)
 198     {
 199         reduce(nUncoupled, sumOp<label>());
 200     }
 201
 202     if (nUncoupled > 0)
 203     {
 204         Info<< "Uncoupled " << nUncoupled << " faces on coupled patches. "
 205             << "(processorPolyPatch, cyclicPolyPatch)" << nl
 206             << "You might need to run couplePatches to restore the patch face"
 207             << " ordering." << nl << endl;
 208     }
 209 }
 210
 211
 212 labelList Foam::fvMeshSubset::subset
 213 (
 214     const label nElems,
 215     const labelList& selectedElements,
 216     const labelList& subsetMap
 217 )
 218 {
 219     // Mark selected elements.
 220     boolList selected(nElems, false);
 221     forAll (selectedElements, i)
 222     {
 223         selected[selectedElements[i]] = true;
 224     }
 225
 226     // Count subset of selected elements
 227     label n = 0;
 228     forAll (subsetMap, i)
 229     {
 230         if (selected[subsetMap[i]])
 231         {
 232             n++;
 233         }
 234     }
 235
 236     // Collect selected elements
 237     labelList subsettedElements(n);
 238     n = 0;
 239
 240     forAll(subsetMap, i)
 241     {
 242         if (selected[subsetMap[i]])
 243         {
 244             subsettedElements[n++] = i;
 245         }
 246     }
 247
 248     return subsettedElements;
 249 }
 250
 251
 252 void Foam::fvMeshSubset::subsetZones()
 253 {
 254     // Keep all zones, even if zero size.
 255
 256     const pointZoneMesh& pointZones = baseMesh().pointZones();
 257
 258     // PointZones
 259     List<pointZone*> pZonePtrs(pointZones.size());
 260
 261     forAll(pointZones, i)
 262     {
 263         const pointZone& pz = pointZones[i];
 264
 265         pZonePtrs[i] = new pointZone
 266         (
 267             pz.name(),
 268             subset(baseMesh().allPoints().size(), pz, pointMap()),
 269             i,
 270             fvMeshSubsetPtr_->pointZones()
 271         );
 272     }
 273
 274
 275     // FaceZones
 276
 277     const faceZoneMesh& faceZones = baseMesh().faceZones();
 278
 279
 280     // Do we need to remove zones where the side we're interested in
 281     // no longer exists? Guess not.
 282     List<faceZone*> fZonePtrs(faceZones.size());
 283
 284     forAll(faceZones, i)
 285     {
 286         const faceZone& fz = faceZones[i];
 287
 288         // Create list of mesh faces part of the new zone
 289         labelList subAddressing
 290         (
 291             subset
 292             (
 293                 baseMesh().allFaces().size(),
 294                 fz,
 295                 faceMap()
 296             )
 297         );
 298
 299         // Flipmap for all mesh faces
 300         boolList fullFlipStatus(baseMesh().allFaces().size(), false);
 301         forAll(fz, j)
 302         {
 303             fullFlipStatus[fz[j]] = fz.flipMap()[j];
 304         }
 305         // Extract sub part
 306         boolList subFlipStatus(subAddressing.size(), false);
 307         forAll(subAddressing, j)
 308         {
 309             subFlipStatus[j] = fullFlipStatus[faceMap()[subAddressing[j]]];
 310         }
 311
 312         fZonePtrs[i] = new faceZone
 313         (
 314             fz.name(),
 315             subAddressing,
 316             subFlipStatus,
 317             i,
 318             fvMeshSubsetPtr_->faceZones()
 319         );
 320     }
 321
 322
 323     const cellZoneMesh& cellZones = baseMesh().cellZones();
 324
 325     List<cellZone*> cZonePtrs(cellZones.size());
 326
 327     forAll(cellZones, i)
 328     {
 329         const cellZone& cz = cellZones[i];
 330
 331         cZonePtrs[i] = new cellZone
 332         (
 333             cz.name(),
 334             subset(baseMesh().nCells(), cz, cellMap()),
 335             i,
 336             fvMeshSubsetPtr_->cellZones()
 337         );
 338     }
 339
 340
 341     // Add the zones
 342     fvMeshSubsetPtr_->addZones(pZonePtrs, fZonePtrs, cZonePtrs);
 343 }
 344
 345
 346 void Foam::fvMeshSubset::makeFvDictionaries()
 347 {
 348     // Try accessing dictionaries
 349     objectRegistry obr
 350     (
 351         IOobject
 352         (
 353             this->name() + "Subset",
 354             baseMesh().time().timeName(),
 355             baseMesh().time(),
 356             IOobject::NO_READ,
 357             IOobject::NO_WRITE
 358         )
 359     );
 360
 361     IOdictionary fvSchemes
 362     (
 363         IOobject
 364         (
 365             "fvSchemes",
 366             obr.time().system(),
 367             obr,
 368             IOobject::NO_READ,
 369             IOobject::NO_WRITE
 370         )
 371     );
 372
 373     if (!fvSchemes.headerOk())
 374     {
 375         Info << "Cannot read " << fvSchemes.path()
 376             << ".  Copy from base" << endl;
 377
 378         IOdictionary fvSchemesBase
 379         (
 380             IOobject
 381             (
 382                 "fvSchemes",
 383                 baseMesh().time().system(),
 384                 baseMesh().time(),
 385                 IOobject::MUST_READ,
 386                 IOobject::NO_WRITE
 387             )
 388         );
 389
 390         fvSchemes = fvSchemesBase;
 391         fvSchemes.regIOobject::write();
 392     }
 393
 394     IOdictionary fvSolution
 395     (
 396         IOobject
 397         (
 398             "fvSolution",
 399             obr.time().system(),
 400             obr,
 401             IOobject::NO_READ,
 402             IOobject::NO_WRITE
 403         )
 404     );
 405
 406     if (!fvSolution.headerOk())
 407     {
 408         Info << "Cannot read " << fvSolution.path()
 409             << ".  Copy from base" << endl;
 410
 411         IOdictionary fvSolutionBase
 412         (
 413             IOobject
 414             (
 415                 "fvSolution",
 416                 baseMesh().time().system(),
 417                 baseMesh().time(),
 418                 IOobject::MUST_READ,
 419                 IOobject::NO_WRITE
 420             )
 421         );
 422
 423         fvSolution = fvSolutionBase;
 424         fvSolution.regIOobject::write();
 425     }
 426 }
 427
 428
 429 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 430
 431 // Construct from components
 432 Foam::fvMeshSubset::fvMeshSubset
 433 (
 434     const IOobject& io,
 435     const fvMesh& baseMesh
 436 )
 437 :
 438     regIOobject(io),
 439     baseMesh_(baseMesh),
 440     fvMeshSubsetPtr_(NULL),
 441     pointMeshSubsetPtr_(NULL),
 442     pointMap_(0),
 443     faceMap_(0),
 444     cellMap_(0),
 445     patchMap_(0)
 446 {}
 447
 448
 449 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 450
 451 Foam::fvMeshSubset::~fvMeshSubset()
 452 {
 453     deleteDemandDrivenData(fvMeshSubsetPtr_);
 454     deleteDemandDrivenData(pointMeshSubsetPtr_);
 455 }
 456
 457
 458 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 459
 460 void Foam::fvMeshSubset::setCellSubset
 461 (
 462     const labelHashSet& globalCellMap,
 463     const label patchID
 464 )
 465 {
 466     // Initial check on patches before doing anything time consuming.
 467     const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
 468     const cellList& oldCells = baseMesh().cells();
 469     const faceList& oldFaces = baseMesh().allFaces();
 470     const pointField& oldPoints = baseMesh().points();
 471     const labelList& oldOwner = baseMesh().faceOwner();
 472     const labelList& oldNeighbour = baseMesh().faceNeighbour();
 473
 474     label wantedPatchID = patchID;
 475
 476     if (wantedPatchID == -1)
 477     {
 478         // No explicit patch specified. Put in oldInternalFaces patch.
 479         // Check if patch with this name already exists.
 480         wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
 481     }
 482     else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
 483     {
 484         FatalErrorIn
 485         (
 486             "fvMeshSubset::setCellSubset(const labelHashSet&"
 487             ", const label patchID)"
 488         )   << "Non-existing patch index " << wantedPatchID << endl
 489             << "Should be between 0 and " << oldPatches.size()-1
 490             << abort(FatalError);
 491     }
 492
 493
 494     cellMap_ = globalCellMap.toc();
 495
 496     // Sort the cell map in the ascending order
 497     sort(cellMap_);
 498
 499     label nCellsInSet = cellMap_.size();
 500
 501     // Mark all used faces
 502
 503     Map<label> facesToSubset(primitiveMesh::facesPerCell_*nCellsInSet);
 504
 505     forAll (cellMap_, cellI)
 506     {
 507         // Mark all faces from the cell
 508         const labelList& curFaces = oldCells[cellMap_[cellI]];
 509
 510         forAll (curFaces, faceI)
 511         {
 512             if (!facesToSubset.found(curFaces[faceI]))
 513             {
 514                 facesToSubset.insert(curFaces[faceI], 1);
 515             }
 516             else
 517             {
 518                 facesToSubset[curFaces[faceI]]++;
 519             }
 520         }
 521     }
 522
 523     // Mark all used points and make a global-to-local face map
 524     Map<label> globalFaceMap(facesToSubset.size());
 525
 526     // Make a global-to-local point map
 527     Map<label> globalPointMap
 528     (
 529         primitiveMesh::pointsPerFace_*facesToSubset.size()
 530     );
 531
 532     // This is done in two goes, so that the boundary faces are last
 533     // in the list.  Because of this, I need to create the face map
 534     // along the way rather than just grab the table of contents.
 535     labelList facesToc = facesToSubset.toc();
 536     sort(facesToc);
 537     faceMap_.setSize(facesToc.size());
 538
 539     // 1. Get all faces that will be internal to the submesh.
 540     forAll (facesToc, faceI)
 541     {
 542         if (facesToSubset[facesToc[faceI]] == 2)
 543         {
 544             // Mark face and increment number of points in set
 545             faceMap_[globalFaceMap.size()] = facesToc[faceI];
 546             globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
 547
 548             // Mark all points from the face
 549             markPoints(oldFaces[facesToc[faceI]], globalPointMap);
 550         }
 551     }
 552
 553     // These are all the internal faces in the mesh.
 554     label nInternalFaces = globalFaceMap.size();
 555
 556
 557     // Where to insert old internal faces.
 558     label oldPatchStart = labelMax;
 559     if (wantedPatchID != -1)
 560     {
 561         oldPatchStart = oldPatches[wantedPatchID].start();
 562     }
 563
 564
 565     label faceI = 0;
 566
 567     // 2. Boundary faces up to where we want to insert old internal faces
 568     for (; faceI< facesToc.size(); faceI++)
 569     {
 570         if (facesToc[faceI] >= oldPatchStart)
 571         {
 572             break;
 573         }
 574         if
 575         (
 576             !baseMesh().isInternalFace(facesToc[faceI])
 577          && facesToSubset[facesToc[faceI]] == 1
 578         )
 579         {
 580             // Mark face and increment number of points in set
 581             faceMap_[globalFaceMap.size()] = facesToc[faceI];
 582             globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
 583
 584             // Mark all points from the face
 585             markPoints(oldFaces[facesToc[faceI]], globalPointMap);
 586         }
 587     }
 588
 589     // 3. old internal faces
 590     forAll(facesToc, intFaceI)
 591     {
 592         if
 593         (
 594             baseMesh().isInternalFace(facesToc[intFaceI])
 595          && facesToSubset[facesToc[intFaceI]] == 1
 596         )
 597         {
 598             // Mark face and increment number of points in set
 599             faceMap_[globalFaceMap.size()] = facesToc[intFaceI];
 600             globalFaceMap.insert(facesToc[intFaceI], globalFaceMap.size());
 601
 602             // Mark all points from the face
 603             markPoints(oldFaces[facesToc[intFaceI]], globalPointMap);
 604         }
 605     }
 606
 607     // 4. Remaining boundary faces
 608     for (; faceI< facesToc.size(); faceI++)
 609     {
 610         if
 611         (
 612             !baseMesh().isInternalFace(facesToc[faceI])
 613          && facesToSubset[facesToc[faceI]] == 1
 614         )
 615         {
 616             // Mark face and increment number of points in set
 617             faceMap_[globalFaceMap.size()] = facesToc[faceI];
 618             globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
 619
 620             // Mark all points from the face
 621             markPoints(oldFaces[facesToc[faceI]], globalPointMap);
 622         }
 623     }
 624
 625
 626
 627     // Grab the points map
 628     pointMap_ = globalPointMap.toc();
 629     sort(pointMap_);
 630
 631     forAll (pointMap_, pointI)
 632     {
 633         globalPointMap[pointMap_[pointI]] = pointI;
 634     }
 635
 636     if (fvMesh::debug)
 637     {
 638         Pout<< "Number of cells in new mesh: " << nCellsInSet << nl
 639             << "Number of faces in new mesh: " << globalFaceMap.size() << nl
 640             << "Number of points in new mesh: " << globalPointMap.size()
 641             << endl;
 642     }
 643
 644     // Make a new mesh
 645     pointField newPoints(globalPointMap.size());
 646
 647     label nNewPoints = 0;
 648
 649     forAll (pointMap_, pointI)
 650     {
 651         newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
 652         nNewPoints++;
 653     }
 654
 655     faceList newFaces(globalFaceMap.size());
 656
 657     label nNewFaces = 0;
 658
 659     // Make internal faces
 660     for (label faceI = 0; faceI < nInternalFaces; faceI++)
 661     {
 662         const face& oldF = oldFaces[faceMap_[faceI]];
 663
 664         face newF(oldF.size());
 665
 666         forAll (newF, i)
 667         {
 668             newF[i] = globalPointMap[oldF[i]];
 669         }
 670
 671         newFaces[nNewFaces] = newF;
 672         nNewFaces++;
 673     }
 674
 675     // Make boundary faces
 676
 677     label nbSize = oldPatches.size();
 678     label oldInternalPatchID  = -1;
 679
 680     if (wantedPatchID == -1)
 681     {
 682         // Create 'oldInternalFaces' patch at the end
 683         // and put all exposed internal faces in there.
 684         oldInternalPatchID = nbSize;
 685         nbSize++;
 686
 687     }
 688     else
 689     {
 690         oldInternalPatchID = wantedPatchID;
 691     }
 692
 693
 694     // Grad size and start of each patch on the fly.  Because of the
 695     // structure of the underlying mesh, the patches will appear in the
 696     // ascending order
 697     labelList boundaryPatchSizes(nbSize, 0);
 698
 699     // Assign boundary faces. Visited in order of faceMap_.
 700     for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
 701     {
 702         label oldFaceI = faceMap_[faceI];
 703
 704         face oldF = oldFaces[oldFaceI];
 705
 706         // Turn the faces as necessary to point outwards
 707         if (baseMesh().isInternalFace(oldFaceI))
 708         {
 709             // Internal face. Possibly turned the wrong way round
 710             if
 711             (
 712                 !globalCellMap.found(oldOwner[oldFaceI])
 713              && globalCellMap.found(oldNeighbour[oldFaceI])
 714             )
 715             {
 716                 oldF = oldFaces[oldFaceI].reverseFace();
 717             }
 718
 719             // Update count for patch
 720             boundaryPatchSizes[oldInternalPatchID]++;
 721         }
 722         else
 723         {
 724             // Boundary face. Increment the appropriate patch
 725             label patchOfFace = oldPatches.whichPatch(oldFaceI);
 726
 727             // Update count for patch
 728             boundaryPatchSizes[patchOfFace]++;
 729         }
 730
 731         face newF(oldF.size());
 732
 733         forAll (newF, i)
 734         {
 735             newF[i] = globalPointMap[oldF[i]];
 736         }
 737
 738         newFaces[nNewFaces] = newF;
 739         nNewFaces++;
 740     }
 741
 742
 743
 744     // Create cells
 745     cellList newCells(nCellsInSet);
 746
 747     label nNewCells = 0;
 748
 749     forAll (cellMap_, cellI)
 750     {
 751         const labelList& oldC = oldCells[cellMap_[cellI]];
 752
 753         labelList newC(oldC.size());
 754
 755         forAll (newC, i)
 756         {
 757             newC[i] = globalFaceMap[oldC[i]];
 758         }
 759
 760         newCells[nNewCells] = cell(newC);
 761         nNewCells++;
 762     }
 763
 764
 765     // Make a new mesh
 766     fvMeshSubsetPtr_ = new fvMesh
 767     (
 768         IOobject
 769         (
 770             this->name() + "Subset",
 771             baseMesh().time().timeName(),
 772             baseMesh().time(),
 773             IOobject::NO_READ,
 774             IOobject::NO_WRITE
 775         ),
 776         xferMove(newPoints),
 777         xferMove(newFaces),
 778         xferMove(newCells)
 779     );
 780
 781     // Clear point mesh
 782     deleteDemandDrivenData(pointMeshSubsetPtr_);
 783
 784
 785     // Add old patches
 786     List<polyPatch*> newBoundary(nbSize);
 787     patchMap_.setSize(nbSize);
 788     label nNewPatches = 0;
 789     label patchStart = nInternalFaces;
 790
 791     forAll(oldPatches, patchI)
 792     {
 793         if (boundaryPatchSizes[patchI] > 0)
 794         {
 795             // Patch still exists. Add it
 796             newBoundary[nNewPatches] = oldPatches[patchI].clone
 797             (
 798                 fvMeshSubsetPtr_->boundaryMesh(),
 799                 nNewPatches,
 800                 boundaryPatchSizes[patchI],
 801                 patchStart
 802             ).ptr();
 803
 804             patchStart += boundaryPatchSizes[patchI];
 805             patchMap_[nNewPatches] = patchI;
 806             nNewPatches++;
 807         }
 808     }
 809
 810     if (wantedPatchID == -1)
 811     {
 812         // Newly created patch so is at end. Check if any faces in it.
 813         if (boundaryPatchSizes[oldInternalPatchID] > 0)
 814         {
 815             newBoundary[nNewPatches] = new polyPatch
 816             (
 817                 "oldInternalFaces",
 818                 boundaryPatchSizes[oldInternalPatchID],
 819                 patchStart,
 820                 nNewPatches,
 821                 fvMeshSubsetPtr_->boundaryMesh()
 822             );
 823
 824             // The index for the first patch is -1 as it originates from
 825             // the internal faces
 826             patchMap_[nNewPatches] = -1;
 827             nNewPatches++;
 828         }
 829     }
 830
 831     // Reset the patch lists
 832     newBoundary.setSize(nNewPatches);
 833     patchMap_.setSize(nNewPatches);
 834
 835     // Add the fvPatches
 836     fvMeshSubsetPtr_->addFvPatches(newBoundary);
 837
 838     // Subset and add any zones
 839     subsetZones();
 840 }
 841
 842
 843 void Foam::fvMeshSubset::setLargeCellSubset
 844 (
 845     const labelList& region,
 846     const label currentRegion,
 847     const label patchID,
 848     const bool syncPar
 849 )
 850 {
 851     const cellList& oldCells = baseMesh().cells();
 852     const faceList& oldFaces = baseMesh().faces();
 853     const pointField& oldPoints = baseMesh().points();
 854     const labelList& oldOwner = baseMesh().faceOwner();
 855     const labelList& oldNeighbour = baseMesh().faceNeighbour();
 856     const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
 857     const label oldNInternalFaces = baseMesh().nInternalFaces();
 858
 859     // Initial checks
 860
 861     if (region.size() != oldCells.size())
 862     {
 863         FatalErrorIn
 864         (
 865             "fvMeshSubset::setCellSubset(const labelList&"
 866             ", const label, const label, const bool)"
 867         )   << "Size of region " << region.size()
 868             << " is not equal to number of cells in mesh " << oldCells.size()
 869             << abort(FatalError);
 870     }
 871
 872
 873     label wantedPatchID = patchID;
 874
 875     if (wantedPatchID == -1)
 876     {
 877         // No explicit patch specified. Put in oldInternalFaces patch.
 878         // Check if patch with this name already exists.
 879         wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
 880     }
 881     else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
 882     {
 883         FatalErrorIn
 884         (
 885             "fvMeshSubset::setCellSubset(const labelList&"
 886             ", const label, const label, const bool)"
 887         )   << "Non-existing patch index " << wantedPatchID << endl
 888             << "Should be between 0 and " << oldPatches.size()-1
 889             << abort(FatalError);
 890     }
 891
 892
 893     // Get the cells for the current region.
 894     cellMap_.setSize(oldCells.size());
 895     label nCellsInSet = 0;
 896
 897     forAll (region, oldCellI)
 898     {
 899         if (region[oldCellI] == currentRegion)
 900         {
 901             cellMap_[nCellsInSet++] = oldCellI;
 902         }
 903     }
 904     cellMap_.setSize(nCellsInSet);
 905
 906
 907     // Mark all used faces. Count number of cells using them
 908     // 0: face not used anymore
 909     // 1: face used by one cell, face becomes/stays boundary face
 910     // 2: face still used and remains internal face
 911     // 3: face coupled and used by one cell only (so should become normal,
 912     //    non-coupled patch face)
 913     //
 914     // Note that this is not really nessecary - but means we can size things
 915     // correctly. Also makes handling coupled faces much easier.
 916
 917     labelList nCellsUsingFace(oldFaces.size(), 0);
 918
 919     label nFacesInSet = 0;
 920     forAll (oldFaces, oldFaceI)
 921     {
 922         bool faceUsed = false;
 923
 924         if (region[oldOwner[oldFaceI]] == currentRegion)
 925         {
 926             nCellsUsingFace[oldFaceI]++;
 927             faceUsed = true;
 928         }
 929
 930         if
 931         (
 932             baseMesh().isInternalFace(oldFaceI)
 933          && (region[oldNeighbour[oldFaceI]] == currentRegion)
 934         )
 935         {
 936             nCellsUsingFace[oldFaceI]++;
 937             faceUsed = true;
 938         }
 939
 940         if (faceUsed)
 941         {
 942             nFacesInSet++;
 943         }
 944     }
 945     faceMap_.setSize(nFacesInSet);
 946
 947     // Handle coupled faces. Modifies patch faces to be uncoupled to 3.
 948     doCoupledPatches(syncPar, nCellsUsingFace);
 949
 950
 951     // See which patch to use for exposed internal faces.
 952     label oldInternalPatchID = 0;
 953
 954     // Insert faces before which patch
 955     label nextPatchID = oldPatches.size();
 956
 957     // old to new patches
 958     labelList globalPatchMap(oldPatches.size());
 959
 960     // New patch size
 961     label nbSize = oldPatches.size();
 962
 963     if (wantedPatchID == -1)
 964     {
 965         // Create 'oldInternalFaces' patch at the end (or before
 966         // processorPatches)
 967         // and put all exposed internal faces in there.
 968
 969         forAll(oldPatches, patchI)
 970         {
 971             if (isA<processorPolyPatch>(oldPatches[patchI]))
 972             {
 973                 nextPatchID = patchI;
 974                 break;
 975             }
 976             oldInternalPatchID++;
 977         }
 978
 979         nbSize++;
 980
 981         // adapt old to new patches for inserted patch
 982         for (label oldPatchI = 0; oldPatchI < nextPatchID; oldPatchI++)
 983         {
 984             globalPatchMap[oldPatchI] = oldPatchI;
 985         }
 986         for
 987         (
 988             label oldPatchI = nextPatchID;
 989             oldPatchI < oldPatches.size();
 990             oldPatchI++
 991         )
 992         {
 993             globalPatchMap[oldPatchI] = oldPatchI+1;
 994         }
 995     }
 996     else
 997     {
 998         oldInternalPatchID = wantedPatchID;
 999         nextPatchID = wantedPatchID+1;
1000
1001         // old to new patches
1002         globalPatchMap = identity(oldPatches.size());
1003     }
1004
1005     labelList boundaryPatchSizes(nbSize, 0);
1006
1007
1008     // Make a global-to-local point map
1009     labelList globalPointMap(oldPoints.size(), -1);
1010
1011     labelList globalFaceMap(oldFaces.size(), -1);
1012     label faceI = 0;
1013
1014     // 1. Pick up all preserved internal faces.
1015     for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
1016     {
1017         if (nCellsUsingFace[oldFaceI] == 2)
1018         {
1019             globalFaceMap[oldFaceI] = faceI;
1020             faceMap_[faceI++] = oldFaceI;
1021
1022             // Mark all points from the face
1023             markPoints(oldFaces[oldFaceI], globalPointMap);
1024         }
1025     }
1026
1027     // These are all the internal faces in the mesh.
1028     label nInternalFaces = faceI;
1029
1030     // 2. Boundary faces up to where we want to insert old internal faces
1031     for
1032     (
1033         label oldPatchI = 0;
1034         oldPatchI < oldPatches.size()
1035      && oldPatchI < nextPatchID;
1036         oldPatchI++
1037     )
1038     {
1039         const polyPatch& oldPatch = oldPatches[oldPatchI];
1040
1041         label oldFaceI = oldPatch.start();
1042
1043         forAll (oldPatch, i)
1044         {
1045             if (nCellsUsingFace[oldFaceI] == 1)
1046             {
1047                 // Boundary face is kept.
1048
1049                 // Mark face and increment number of points in set
1050                 globalFaceMap[oldFaceI] = faceI;
1051                 faceMap_[faceI++] = oldFaceI;
1052
1053                 // Mark all points from the face
1054                 markPoints(oldFaces[oldFaceI], globalPointMap);
1055
1056                 // Increment number of patch faces
1057                 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1058             }
1059             oldFaceI++;
1060         }
1061     }
1062
1063     // 3a. old internal faces that have become exposed.
1064     for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
1065     {
1066         if (nCellsUsingFace[oldFaceI] == 1)
1067         {
1068             globalFaceMap[oldFaceI] = faceI;
1069             faceMap_[faceI++] = oldFaceI;
1070
1071             // Mark all points from the face
1072             markPoints(oldFaces[oldFaceI], globalPointMap);
1073
1074             // Increment number of patch faces
1075             boundaryPatchSizes[oldInternalPatchID]++;
1076         }
1077     }
1078
1079     // 3b. coupled patch faces that have become uncoupled.
1080     for
1081     (
1082         label oldFaceI = oldNInternalFaces;
1083         oldFaceI < oldFaces.size();
1084         oldFaceI++
1085     )
1086     {
1087         if (nCellsUsingFace[oldFaceI] == 3)
1088         {
1089             globalFaceMap[oldFaceI] = faceI;
1090             faceMap_[faceI++] = oldFaceI;
1091
1092             // Mark all points from the face
1093             markPoints(oldFaces[oldFaceI], globalPointMap);
1094
1095             // Increment number of patch faces
1096             boundaryPatchSizes[oldInternalPatchID]++;
1097         }
1098     }
1099
1100     // 4. Remaining boundary faces
1101     for
1102     (
1103         label oldPatchI = nextPatchID;
1104         oldPatchI < oldPatches.size();
1105         oldPatchI++
1106     )
1107     {
1108         const polyPatch& oldPatch = oldPatches[oldPatchI];
1109
1110         label oldFaceI = oldPatch.start();
1111
1112         forAll (oldPatch, i)
1113         {
1114             if (nCellsUsingFace[oldFaceI] == 1)
1115             {
1116                 // Boundary face is kept.
1117
1118                 // Mark face and increment number of points in set
1119                 globalFaceMap[oldFaceI] = faceI;
1120                 faceMap_[faceI++] = oldFaceI;
1121
1122                 // Mark all points from the face
1123                 markPoints(oldFaces[oldFaceI], globalPointMap);
1124
1125                 // Increment number of patch faces
1126                 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1127             }
1128             oldFaceI++;
1129         }
1130     }
1131
1132     if (faceI != nFacesInSet)
1133     {
1134         FatalErrorIn
1135         (
1136             "fvMeshSubset::setCellSubset(const labelList&"
1137             ", const label, const label, const bool)"
1138         )   << "Problem" << abort(FatalError);
1139     }
1140
1141
1142     // Grab the points map
1143     label nPointsInSet = 0;
1144
1145     forAll(globalPointMap, pointI)
1146     {
1147         if (globalPointMap[pointI] != -1)
1148         {
1149             nPointsInSet++;
1150         }
1151     }
1152     pointMap_.setSize(nPointsInSet);
1153
1154     nPointsInSet = 0;
1155
1156     forAll(globalPointMap, pointI)
1157     {
1158         if (globalPointMap[pointI] != -1)
1159         {
1160             pointMap_[nPointsInSet] = pointI;
1161             globalPointMap[pointI] = nPointsInSet;
1162             nPointsInSet++;
1163         }
1164     }
1165
1166     Pout<< "Number of cells in new mesh : " << cellMap_.size() << endl;
1167     Pout<< "Number of faces in new mesh : " << faceMap_.size() << endl;
1168     Pout<< "Number of points in new mesh: " << pointMap_.size() << endl;
1169
1170     // Make a new mesh
1171     pointField newPoints(pointMap_.size());
1172
1173     label nNewPoints = 0;
1174
1175     forAll (pointMap_, pointI)
1176     {
1177         newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
1178         nNewPoints++;
1179     }
1180
1181     faceList newFaces(faceMap_.size());
1182
1183     label nNewFaces = 0;
1184
1185     // Make internal faces
1186     for (label faceI = 0; faceI < nInternalFaces; faceI++)
1187     {
1188         const face& oldF = oldFaces[faceMap_[faceI]];
1189
1190         face newF(oldF.size());
1191
1192         forAll (newF, i)
1193         {
1194             newF[i] = globalPointMap[oldF[i]];
1195         }
1196
1197         newFaces[nNewFaces] = newF;
1198         nNewFaces++;
1199     }
1200
1201
1202     // Make boundary faces. (different from internal since might need to be
1203     // flipped)
1204     for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
1205     {
1206         label oldFaceI = faceMap_[faceI];
1207
1208         face oldF = oldFaces[oldFaceI];
1209
1210         // Turn the faces as necessary to point outwards
1211         if (baseMesh().isInternalFace(oldFaceI))
1212         {
1213             // Was internal face. Possibly turned the wrong way round
1214             if
1215             (
1216                 region[oldOwner[oldFaceI]] != currentRegion
1217              && region[oldNeighbour[oldFaceI]] == currentRegion
1218             )
1219             {
1220                 oldF = oldFaces[oldFaceI].reverseFace();
1221             }
1222         }
1223
1224         // Relabel vertices of the (possibly turned) face.
1225         face newF(oldF.size());
1226
1227         forAll (newF, i)
1228         {
1229             newF[i] = globalPointMap[oldF[i]];
1230         }
1231
1232         newFaces[nNewFaces] = newF;
1233         nNewFaces++;
1234     }
1235
1236
1237
1238     // Create cells
1239     cellList newCells(nCellsInSet);
1240
1241     label nNewCells = 0;
1242
1243     forAll (cellMap_, cellI)
1244     {
1245         const labelList& oldC = oldCells[cellMap_[cellI]];
1246
1247         labelList newC(oldC.size());
1248
1249         forAll (newC, i)
1250         {
1251             newC[i] = globalFaceMap[oldC[i]];
1252         }
1253
1254         newCells[nNewCells] = cell(newC);
1255         nNewCells++;
1256     }
1257
1258
1259     // Make a new mesh
1260     makeFvDictionaries();
1261
1262     fvMeshSubsetPtr_ = new fvMesh
1263     (
1264         IOobject
1265         (
1266             this->name() + "Subset",
1267             baseMesh().time().timeName(),
1268             baseMesh().time(),
1269             IOobject::NO_READ,
1270             IOobject::NO_WRITE
1271         ),
1272         xferMove(newPoints),
1273         xferMove(newFaces),
1274         xferMove(newCells),
1275         syncPar           // parallel synchronisation
1276     );
1277
1278     // Clear point mesh
1279     deleteDemandDrivenData(pointMeshSubsetPtr_);
1280
1281     // Add old patches
1282     List<polyPatch*> newBoundary(nbSize);
1283     patchMap_.setSize(nbSize);
1284     label nNewPatches = 0;
1285     label patchStart = nInternalFaces;
1286
1287
1288     // For parallel: only remove patch if none of the processors has it.
1289     // This only gets done for patches before the one being inserted
1290     // (so patches < nextPatchID)
1291
1292     // Get sum of patch sizes. Zero if patch can be deleted.
1293     labelList globalPatchSizes(boundaryPatchSizes);
1294     globalPatchSizes.setSize(nextPatchID);
1295
1296     if (syncPar && Pstream::parRun())
1297     {
1298         // Get patch names (up to nextPatchID)
1299         List<wordList> patchNames(Pstream::nProcs());
1300         patchNames[Pstream::myProcNo()] = oldPatches.names();
1301         patchNames[Pstream::myProcNo()].setSize(nextPatchID);
1302         Pstream::gatherList(patchNames);
1303         Pstream::scatterList(patchNames);
1304
1305         // Get patch sizes (up to nextPatchID).
1306         // Note that up to nextPatchID the globalPatchMap is an identity so
1307         // no need to index through that.
1308         Pstream::listCombineGather(globalPatchSizes, plusEqOp<label>());
1309         Pstream::listCombineScatter(globalPatchSizes);
1310
1311         // Now all processors have all the patchnames.
1312         // Decide: if all processors have the same patch names and size is zero
1313         // everywhere remove the patch.
1314         bool samePatches = true;
1315
1316         for (label procI = 1; procI < patchNames.size(); procI++)
1317         {
1318             if (patchNames[procI] != patchNames[0])
1319             {
1320                 samePatches = false;
1321                 break;
1322             }
1323         }
1324
1325         if (!samePatches)
1326         {
1327             // Patchnames not sync on all processors so disable removal of
1328             // zero sized patches.
1329             globalPatchSizes = labelMax;
1330         }
1331     }
1332
1333
1334     // Old patches
1335
1336     for
1337     (
1338         label oldPatchI = 0;
1339         oldPatchI < oldPatches.size()
1340      && oldPatchI < nextPatchID;
1341         oldPatchI++
1342     )
1343     {
1344         label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1345
1346         // Clone (even if 0 size)
1347         newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1348         (
1349             fvMeshSubsetPtr_->boundaryMesh(),
1350             nNewPatches,
1351             newSize,
1352             patchStart
1353         ).ptr();
1354
1355         patchStart += newSize;
1356         patchMap_[nNewPatches] = oldPatchI;    // compact patchMap
1357         nNewPatches++;
1358     }
1359
1360     // Inserted patch
1361
1362     if (wantedPatchID == -1)
1363     {
1364         label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1365
1366         if (syncPar)
1367         {
1368             reduce(oldInternalSize, sumOp<label>());
1369         }
1370
1371         // Newly created patch so is at end. Check if any faces in it.
1372         if (oldInternalSize > 0)
1373         {
1374             newBoundary[nNewPatches] = new polyPatch
1375             (
1376                 "oldInternalFaces",
1377                 boundaryPatchSizes[oldInternalPatchID],
1378                 patchStart,
1379                 nNewPatches,
1380                 fvMeshSubsetPtr_->boundaryMesh()
1381             );
1382
1383             Pout<< "    oldInternalFaces : "
1384                 << boundaryPatchSizes[oldInternalPatchID] << endl;
1385
1386             // The index for the first patch is -1 as it originates from
1387             // the internal faces
1388             patchStart += boundaryPatchSizes[oldInternalPatchID];
1389             patchMap_[nNewPatches] = -1;
1390             nNewPatches++;
1391         }
1392     }
1393
1394     // Old patches
1395
1396     for
1397     (
1398         label oldPatchI = nextPatchID;
1399         oldPatchI < oldPatches.size();
1400         oldPatchI++
1401     )
1402     {
1403         label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1404
1405         // Patch still exists. Add it
1406         newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1407         (
1408             fvMeshSubsetPtr_->boundaryMesh(),
1409             nNewPatches,
1410             newSize,
1411             patchStart
1412         ).ptr();
1413
1414         //Pout<< "    " << oldPatches[oldPatchI].name() << " : "
1415         //    << newSize << endl;
1416
1417         patchStart += newSize;
1418         patchMap_[nNewPatches] = oldPatchI;    // compact patchMap
1419         nNewPatches++;
1420     }
1421
1422
1423     // Reset the patch lists
1424     newBoundary.setSize(nNewPatches);
1425     patchMap_.setSize(nNewPatches);
1426
1427
1428     // Add the fvPatches
1429     fvMeshSubsetPtr_->addFvPatches(newBoundary, syncPar);
1430
1431     // Subset and add any zones
1432     subsetZones();
1433 }
1434
1435
1436 void Foam::fvMeshSubset::setLargeCellSubset
1437 (
1438     const labelHashSet& globalCellMap,
1439     const label patchID,
1440     const bool syncPar
1441 )
1442 {
1443     labelList region(baseMesh().nCells(), 0);
1444
1445     forAllConstIter (labelHashSet, globalCellMap, iter)
1446     {
1447         region[iter.key()] = 1;
1448     }
1449     setLargeCellSubset(region, 1, patchID, syncPar);
1450 }
1451
1452
1453 const fvMesh& Foam::fvMeshSubset::subMesh() const
1454 {
1455     checkCellSubset();
1456
1457     return *fvMeshSubsetPtr_;
1458 }
1459
1460
1461 fvMesh& Foam::fvMeshSubset::subMesh()
1462 {
1463     checkCellSubset();
1464
1465     return *fvMeshSubsetPtr_;
1466 }
1467
1468
1469 const pointMesh& Foam::fvMeshSubset::subPointMesh() const
1470 {
1471     if (!pointMeshSubsetPtr_)
1472     {
1473         pointMeshSubsetPtr_ = new pointMesh(subMesh());
1474     }
1475
1476     return *pointMeshSubsetPtr_;
1477 }
1478
1479
1480 pointMesh& Foam::fvMeshSubset::subPointMesh()
1481 {
1482     if (!pointMeshSubsetPtr_)
1483     {
1484         pointMeshSubsetPtr_ = new pointMesh(subMesh());
1485     }
1486
1487     return *pointMeshSubsetPtr_;
1488 }
1489
1490
1491 const labelList& Foam::fvMeshSubset::pointMap() const
1492 {
1493     checkCellSubset();
1494
1495     return pointMap_;
1496 }
1497
1498
1499 const labelList& Foam::fvMeshSubset::faceMap() const
1500 {
1501     checkCellSubset();
1502
1503     return faceMap_;
1504 }
1505
1506
1507 const labelList& Foam::fvMeshSubset::cellMap() const
1508 {
1509     checkCellSubset();
1510
1511     return cellMap_;
1512 }
1513
1514
1515 const labelList& Foam::fvMeshSubset::patchMap() const
1516 {
1517     checkCellSubset();
1518
1519     return patchMap_;
1520 }
1521
1522
1523 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1524
1525 } // End namespace Foam
1526
1527 // ************************************************************************* //