src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C

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