src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.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 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_)
  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().allPoints().size(), 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         // Create list of mesh faces part of the new zone
 288         labelList subAddressing
 289         (
 290             subset
 291             (
 292                 baseMesh().allFaces().size(),
 293                 fz,
 294                 faceMap()
 295             )
 296         );
 297
 298         // Flipmap for all mesh faces
 299         boolList fullFlipStatus(baseMesh().allFaces().size(), false);
 300         forAll(fz, j)
 301         {
 302             fullFlipStatus[fz[j]] = fz.flipMap()[j];
 303         }
 304         // Extract sub part
 305         boolList subFlipStatus(subAddressing.size(), false);
 306         forAll(subAddressing, j)
 307         {
 308             subFlipStatus[j] = fullFlipStatus[faceMap()[subAddressing[j]]];
 309         }
 310
 311         fZonePtrs[i] = new faceZone
 312         (
 313             fz.name(),
 314             subAddressing,
 315             subFlipStatus,
 316             i,
 317             fvMeshSubsetPtr_->faceZones()
 318         );
 319     }
 320
 321
 322     const cellZoneMesh& cellZones = baseMesh().cellZones();
 323
 324     List<cellZone*> cZonePtrs(cellZones.size());
 325
 326     forAll(cellZones, i)
 327     {
 328         const cellZone& cz = cellZones[i];
 329
 330         cZonePtrs[i] = new cellZone
 331         (
 332             cz.name(),
 333             subset(baseMesh().nCells(), cz, cellMap()),
 334             i,
 335             fvMeshSubsetPtr_->cellZones()
 336         );
 337     }
 338
 339
 340     // Add the zones
 341     fvMeshSubsetPtr_->addZones(pZonePtrs, fZonePtrs, cZonePtrs);
 342 }
 343
 344
 345 void Foam::fvMeshSubset::makeFvDictionaries()
 346 {
 347     // Try accessing dictionaries
 348     objectRegistry obr
 349     (
 350         IOobject
 351         (
 352             this->name() + "Subset",
 353             baseMesh().time().timeName(),
 354             baseMesh().time(),
 355             IOobject::NO_READ,
 356             IOobject::NO_WRITE
 357         )
 358     );
 359
 360     IOdictionary fvSchemes
 361     (
 362         IOobject
 363         (
 364             "fvSchemes",
 365             obr.time().system(),
 366             obr,
 367             IOobject::NO_READ,
 368             IOobject::NO_WRITE
 369         )
 370     );
 371
 372     if (!fvSchemes.headerOk())
 373     {
 374         Info << "Cannot read " << fvSchemes.path()
 375             << ".  Copy from base" << endl;
 376
 377         IOdictionary fvSchemesBase
 378         (
 379             IOobject
 380             (
 381                 "fvSchemes",
 382                 baseMesh().time().system(),
 383                 baseMesh().time(),
 384                 IOobject::MUST_READ,
 385                 IOobject::NO_WRITE
 386             )
 387         );
 388
 389         fvSchemes = fvSchemesBase;
 390         fvSchemes.regIOobject::write();
 391     }
 392
 393     IOdictionary fvSolution
 394     (
 395         IOobject
 396         (
 397             "fvSolution",
 398             obr.time().system(),
 399             obr,
 400             IOobject::NO_READ,
 401             IOobject::NO_WRITE
 402         )
 403     );
 404
 405     if (!fvSolution.headerOk())
 406     {
 407         Info << "Cannot read " << fvSolution.path()
 408             << ".  Copy from base" << endl;
 409
 410         IOdictionary fvSolutionBase
 411         (
 412             IOobject
 413             (
 414                 "fvSolution",
 415                 baseMesh().time().system(),
 416                 baseMesh().time(),
 417                 IOobject::MUST_READ,
 418                 IOobject::NO_WRITE
 419             )
 420         );
 421
 422         fvSolution = fvSolutionBase;
 423         fvSolution.regIOobject::write();
 424     }
 425 }
 426
 427
 428 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 429
 430 // Construct from components
 431 Foam::fvMeshSubset::fvMeshSubset
 432 (
 433     const IOobject& io,
 434     const fvMesh& baseMesh
 435 )
 436 :
 437     regIOobject(io),
 438     baseMesh_(baseMesh),
 439     fvMeshSubsetPtr_(NULL),
 440     pointMeshSubsetPtr_(NULL),
 441     pointMap_(0),
 442     faceMap_(0),
 443     cellMap_(0),
 444     patchMap_(0)
 445 {}
 446
 447
 448 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 449
 450 Foam::fvMeshSubset::~fvMeshSubset()
 451 {
 452     deleteDemandDrivenData(fvMeshSubsetPtr_);
 453     deleteDemandDrivenData(pointMeshSubsetPtr_);
 454 }
 455
 456
 457 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 458
 459 void Foam::fvMeshSubset::setCellSubset
 460 (
 461     const labelHashSet& globalCellMap,
 462     const label patchID
 463 )
 464 {
 465     // Initial check on patches before doing anything time consuming.
 466     const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
 467     const cellList& oldCells = baseMesh().cells();
 468     const faceList& oldFaces = baseMesh().allFaces();
 469     const pointField& oldPoints = baseMesh().points();
 470     const labelList& oldOwner = baseMesh().faceOwner();
 471     const labelList& oldNeighbour = baseMesh().faceNeighbour();
 472
 473     label wantedPatchID = patchID;
 474
 475     if (wantedPatchID == -1)
 476     {
 477         // No explicit patch specified. Put in oldInternalFaces patch.
 478         // Check if patch with this name already exists.
 479         wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
 480     }
 481     else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
 482     {
 483         FatalErrorIn
 484         (
 485             "fvMeshSubset::setCellSubset(const labelHashSet&"
 486             ", const label patchID)"
 487         )   << "Non-existing patch index " << wantedPatchID << endl
 488             << "Should be between 0 and " << oldPatches.size()-1
 489             << abort(FatalError);
 490     }
 491
 492
 493     cellMap_ = globalCellMap.toc();
 494
 495     // Sort the cell map in the ascending order
 496     sort(cellMap_);
 497
 498     label nCellsInSet = cellMap_.size();
 499
 500     // Mark all used faces
 501
 502     Map<label> facesToSubset(primitiveMesh::facesPerCell_*nCellsInSet);
 503
 504     forAll (cellMap_, cellI)
 505     {
 506         // Mark all faces from the cell
 507         const labelList& curFaces = oldCells[cellMap_[cellI]];
 508
 509         forAll (curFaces, faceI)
 510         {
 511             if (!facesToSubset.found(curFaces[faceI]))
 512             {
 513                 facesToSubset.insert(curFaces[faceI], 1);
 514             }
 515             else
 516             {
 517                 facesToSubset[curFaces[faceI]]++;
 518             }
 519         }
 520     }
 521
 522     // Mark all used points and make a global-to-local face map
 523     Map<label> globalFaceMap(facesToSubset.size());
 524
 525     // Make a global-to-local point map
 526     Map<label> globalPointMap
 527     (
 528         primitiveMesh::pointsPerFace_*facesToSubset.size()
 529     );
 530
 531     // This is done in two goes, so that the boundary faces are last
 532     // in the list.  Because of this, I need to create the face map
 533     // along the way rather than just grab the table of contents.
 534     labelList facesToc = facesToSubset.toc();
 535     sort(facesToc);
 536     faceMap_.setSize(facesToc.size());
 537
 538     // 1. Get all faces that will be internal to the submesh.
 539     forAll (facesToc, faceI)
 540     {
 541         if (facesToSubset[facesToc[faceI]] == 2)
 542         {
 543             // Mark face and increment number of points in set
 544             faceMap_[globalFaceMap.size()] = facesToc[faceI];
 545             globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
 546
 547             // Mark all points from the face
 548             markPoints(oldFaces[facesToc[faceI]], globalPointMap);
 549         }
 550     }
 551
 552     // These are all the internal faces in the mesh.
 553     label nInternalFaces = globalFaceMap.size();
 554
 555
 556     // Where to insert old internal faces.
 557     label oldPatchStart = labelMax;
 558     if (wantedPatchID != -1)
 559     {
 560         oldPatchStart = oldPatches[wantedPatchID].start();
 561     }
 562
 563
 564     label faceI = 0;
 565
 566     // 2. Boundary faces up to where we want to insert old internal faces
 567     for (; faceI< facesToc.size(); faceI++)
 568     {
 569         if (facesToc[faceI] >= oldPatchStart)
 570         {
 571             break;
 572         }
 573         if
 574         (
 575             !baseMesh().isInternalFace(facesToc[faceI])
 576          && facesToSubset[facesToc[faceI]] == 1
 577         )
 578         {
 579             // Mark face and increment number of points in set
 580             faceMap_[globalFaceMap.size()] = facesToc[faceI];
 581             globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
 582
 583             // Mark all points from the face
 584             markPoints(oldFaces[facesToc[faceI]], globalPointMap);
 585         }
 586     }
 587
 588     // 3. old internal faces
 589     forAll(facesToc, intFaceI)
 590     {
 591         if
 592         (
 593             baseMesh().isInternalFace(facesToc[intFaceI])
 594          && facesToSubset[facesToc[intFaceI]] == 1
 595         )
 596         {
 597             // Mark face and increment number of points in set
 598             faceMap_[globalFaceMap.size()] = facesToc[intFaceI];
 599             globalFaceMap.insert(facesToc[intFaceI], globalFaceMap.size());
 600
 601             // Mark all points from the face
 602             markPoints(oldFaces[facesToc[intFaceI]], globalPointMap);
 603         }
 604     }
 605
 606     // 4. Remaining boundary faces
 607     for (; faceI< facesToc.size(); faceI++)
 608     {
 609         if
 610         (
 611             !baseMesh().isInternalFace(facesToc[faceI])
 612          && facesToSubset[facesToc[faceI]] == 1
 613         )
 614         {
 615             // Mark face and increment number of points in set
 616             faceMap_[globalFaceMap.size()] = facesToc[faceI];
 617             globalFaceMap.insert(facesToc[faceI], globalFaceMap.size());
 618
 619             // Mark all points from the face
 620             markPoints(oldFaces[facesToc[faceI]], globalPointMap);
 621         }
 622     }
 623
 624
 625
 626     // Grab the points map
 627     pointMap_ = globalPointMap.toc();
 628     sort(pointMap_);
 629
 630     forAll (pointMap_, pointI)
 631     {
 632         globalPointMap[pointMap_[pointI]] = pointI;
 633     }
 634
 635     if (fvMesh::debug)
 636     {
 637         Pout<< "Number of cells in new mesh: " << nCellsInSet << nl
 638             << "Number of faces in new mesh: " << globalFaceMap.size() << nl
 639             << "Number of points in new mesh: " << globalPointMap.size()
 640             << endl;
 641     }
 642
 643     // Make a new mesh
 644     pointField newPoints(globalPointMap.size());
 645
 646     label nNewPoints = 0;
 647
 648     forAll (pointMap_, pointI)
 649     {
 650         newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
 651         nNewPoints++;
 652     }
 653
 654     faceList newFaces(globalFaceMap.size());
 655
 656     label nNewFaces = 0;
 657
 658     // Make internal faces
 659     for (label faceI = 0; faceI < nInternalFaces; faceI++)
 660     {
 661         const face& oldF = oldFaces[faceMap_[faceI]];
 662
 663         face newF(oldF.size());
 664
 665         forAll (newF, i)
 666         {
 667             newF[i] = globalPointMap[oldF[i]];
 668         }
 669
 670         newFaces[nNewFaces] = newF;
 671         nNewFaces++;
 672     }
 673
 674     // Make boundary faces
 675
 676     label nbSize = oldPatches.size();
 677     label oldInternalPatchID  = -1;
 678
 679     if (wantedPatchID == -1)
 680     {
 681         // Create 'oldInternalFaces' patch at the end
 682         // and put all exposed internal faces in there.
 683         oldInternalPatchID = nbSize;
 684         nbSize++;
 685
 686     }
 687     else
 688     {
 689         oldInternalPatchID = wantedPatchID;
 690     }
 691
 692
 693     // Grad size and start of each patch on the fly.  Because of the
 694     // structure of the underlying mesh, the patches will appear in the
 695     // ascending order
 696     labelList boundaryPatchSizes(nbSize, 0);
 697
 698     // Assign boundary faces. Visited in order of faceMap_.
 699     for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
 700     {
 701         label oldFaceI = faceMap_[faceI];
 702
 703         face oldF = oldFaces[oldFaceI];
 704
 705         // Turn the faces as necessary to point outwards
 706         if (baseMesh().isInternalFace(oldFaceI))
 707         {
 708             // Internal face. Possibly turned the wrong way round
 709             if
 710             (
 711                 !globalCellMap.found(oldOwner[oldFaceI])
 712              && globalCellMap.found(oldNeighbour[oldFaceI])
 713             )
 714             {
 715                 oldF = oldFaces[oldFaceI].reverseFace();
 716             }
 717
 718             // Update count for patch
 719             boundaryPatchSizes[oldInternalPatchID]++;
 720         }
 721         else
 722         {
 723             // Boundary face. Increment the appropriate patch
 724             label patchOfFace = oldPatches.whichPatch(oldFaceI);
 725
 726             // Update count for patch
 727             boundaryPatchSizes[patchOfFace]++;
 728         }
 729
 730         face newF(oldF.size());
 731
 732         forAll (newF, i)
 733         {
 734             newF[i] = globalPointMap[oldF[i]];
 735         }
 736
 737         newFaces[nNewFaces] = newF;
 738         nNewFaces++;
 739     }
 740
 741
 742
 743     // Create cells
 744     cellList newCells(nCellsInSet);
 745
 746     label nNewCells = 0;
 747
 748     forAll (cellMap_, cellI)
 749     {
 750         const labelList& oldC = oldCells[cellMap_[cellI]];
 751
 752         labelList newC(oldC.size());
 753
 754         forAll (newC, i)
 755         {
 756             newC[i] = globalFaceMap[oldC[i]];
 757         }
 758
 759         newCells[nNewCells] = cell(newC);
 760         nNewCells++;
 761     }
 762
 763
 764     // Make a new mesh
 765     fvMeshSubsetPtr_ = new fvMesh
 766     (
 767         IOobject
 768         (
 769             this->name() + "Subset",
 770             baseMesh().time().timeName(),
 771             baseMesh().time(),
 772             IOobject::NO_READ,
 773             IOobject::NO_WRITE
 774         ),
 775         xferMove(newPoints),
 776         xferMove(newFaces),
 777         xferMove(newCells)
 778     );
 779
 780     // Clear point mesh
 781     deleteDemandDrivenData(pointMeshSubsetPtr_);
 782
 783
 784     // Add old patches
 785     List<polyPatch*> newBoundary(nbSize);
 786     patchMap_.setSize(nbSize);
 787     label nNewPatches = 0;
 788     label patchStart = nInternalFaces;
 789
 790     forAll(oldPatches, patchI)
 791     {
 792         if (boundaryPatchSizes[patchI] > 0)
 793         {
 794             // Patch still exists. Add it
 795             newBoundary[nNewPatches] = oldPatches[patchI].clone
 796             (
 797                 fvMeshSubsetPtr_->boundaryMesh(),
 798                 nNewPatches,
 799                 boundaryPatchSizes[patchI],
 800                 patchStart
 801             ).ptr();
 802
 803             patchStart += boundaryPatchSizes[patchI];
 804             patchMap_[nNewPatches] = patchI;
 805             nNewPatches++;
 806         }
 807     }
 808
 809     if (wantedPatchID == -1)
 810     {
 811         // Newly created patch so is at end. Check if any faces in it.
 812         if (boundaryPatchSizes[oldInternalPatchID] > 0)
 813         {
 814             newBoundary[nNewPatches] = new polyPatch
 815             (
 816                 "oldInternalFaces",
 817                 boundaryPatchSizes[oldInternalPatchID],
 818                 patchStart,
 819                 nNewPatches,
 820                 fvMeshSubsetPtr_->boundaryMesh()
 821             );
 822
 823             // The index for the first patch is -1 as it originates from
 824             // the internal faces
 825             patchMap_[nNewPatches] = -1;
 826             nNewPatches++;
 827         }
 828     }
 829
 830     // Reset the patch lists
 831     newBoundary.setSize(nNewPatches);
 832     patchMap_.setSize(nNewPatches);
 833
 834     // Add the fvPatches
 835     fvMeshSubsetPtr_->addFvPatches(newBoundary);
 836
 837     // Subset and add any zones
 838     subsetZones();
 839 }
 840
 841
 842 void Foam::fvMeshSubset::setLargeCellSubset
 843 (
 844     const labelList& region,
 845     const label currentRegion,
 846     const label patchID,
 847     const bool syncPar
 848 )
 849 {
 850     const cellList& oldCells = baseMesh().cells();
 851     const faceList& oldFaces = baseMesh().faces();
 852     const pointField& oldPoints = baseMesh().points();
 853     const labelList& oldOwner = baseMesh().faceOwner();
 854     const labelList& oldNeighbour = baseMesh().faceNeighbour();
 855     const polyBoundaryMesh& oldPatches = baseMesh().boundaryMesh();
 856     const label oldNInternalFaces = baseMesh().nInternalFaces();
 857
 858     // Initial checks
 859
 860     if (region.size() != oldCells.size())
 861     {
 862         FatalErrorIn
 863         (
 864             "fvMeshSubset::setCellSubset(const labelList&"
 865             ", const label, const label, const bool)"
 866         )   << "Size of region " << region.size()
 867             << " is not equal to number of cells in mesh " << oldCells.size()
 868             << abort(FatalError);
 869     }
 870
 871
 872     label wantedPatchID = patchID;
 873
 874     if (wantedPatchID == -1)
 875     {
 876         // No explicit patch specified. Put in oldInternalFaces patch.
 877         // Check if patch with this name already exists.
 878         wantedPatchID = oldPatches.findPatchID("oldInternalFaces");
 879     }
 880     else if (wantedPatchID < 0 || wantedPatchID >= oldPatches.size())
 881     {
 882         FatalErrorIn
 883         (
 884             "fvMeshSubset::setCellSubset(const labelList&"
 885             ", const label, const label, const bool)"
 886         )   << "Non-existing patch index " << wantedPatchID << endl
 887             << "Should be between 0 and " << oldPatches.size()-1
 888             << abort(FatalError);
 889     }
 890
 891
 892     // Get the cells for the current region.
 893     cellMap_.setSize(oldCells.size());
 894     label nCellsInSet = 0;
 895
 896     forAll (region, oldCellI)
 897     {
 898         if (region[oldCellI] == currentRegion)
 899         {
 900             cellMap_[nCellsInSet++] = oldCellI;
 901         }
 902     }
 903     cellMap_.setSize(nCellsInSet);
 904
 905
 906     // Mark all used faces. Count number of cells using them
 907     // 0: face not used anymore
 908     // 1: face used by one cell, face becomes/stays boundary face
 909     // 2: face still used and remains internal face
 910     // 3: face coupled and used by one cell only (so should become normal,
 911     //    non-coupled patch face)
 912     //
 913     // Note that this is not really nessecary - but means we can size things
 914     // correctly. Also makes handling coupled faces much easier.
 915
 916     labelList nCellsUsingFace(oldFaces.size(), 0);
 917
 918     label nFacesInSet = 0;
 919     forAll (oldFaces, oldFaceI)
 920     {
 921         bool faceUsed = false;
 922
 923         if (region[oldOwner[oldFaceI]] == currentRegion)
 924         {
 925             nCellsUsingFace[oldFaceI]++;
 926             faceUsed = true;
 927         }
 928
 929         if
 930         (
 931             baseMesh().isInternalFace(oldFaceI)
 932          && (region[oldNeighbour[oldFaceI]] == currentRegion)
 933         )
 934         {
 935             nCellsUsingFace[oldFaceI]++;
 936             faceUsed = true;
 937         }
 938
 939         if (faceUsed)
 940         {
 941             nFacesInSet++;
 942         }
 943     }
 944     faceMap_.setSize(nFacesInSet);
 945
 946     // Handle coupled faces. Modifies patch faces to be uncoupled to 3.
 947     doCoupledPatches(syncPar, nCellsUsingFace);
 948
 949
 950     // See which patch to use for exposed internal faces.
 951     label oldInternalPatchID = 0;
 952
 953     // Insert faces before which patch
 954     label nextPatchID = oldPatches.size();
 955
 956     // old to new patches
 957     labelList globalPatchMap(oldPatches.size());
 958
 959     // New patch size
 960     label nbSize = oldPatches.size();
 961
 962     if (wantedPatchID == -1)
 963     {
 964         // Create 'oldInternalFaces' patch at the end (or before
 965         // processorPatches)
 966         // and put all exposed internal faces in there.
 967
 968         forAll(oldPatches, patchI)
 969         {
 970             if (isA<processorPolyPatch>(oldPatches[patchI]))
 971             {
 972                 nextPatchID = patchI;
 973                 break;
 974             }
 975             oldInternalPatchID++;
 976         }
 977
 978         nbSize++;
 979
 980         // adapt old to new patches for inserted patch
 981         for (label oldPatchI = 0; oldPatchI < nextPatchID; oldPatchI++)
 982         {
 983             globalPatchMap[oldPatchI] = oldPatchI;
 984         }
 985         for
 986         (
 987             label oldPatchI = nextPatchID;
 988             oldPatchI < oldPatches.size();
 989             oldPatchI++
 990         )
 991         {
 992             globalPatchMap[oldPatchI] = oldPatchI+1;
 993         }
 994     }
 995     else
 996     {
 997         oldInternalPatchID = wantedPatchID;
 998         nextPatchID = wantedPatchID+1;
 999
1000         // old to new patches
1001         globalPatchMap = identity(oldPatches.size());
1002     }
1003
1004     labelList boundaryPatchSizes(nbSize, 0);
1005
1006
1007     // Make a global-to-local point map
1008     labelList globalPointMap(oldPoints.size(), -1);
1009
1010     labelList globalFaceMap(oldFaces.size(), -1);
1011     label faceI = 0;
1012
1013     // 1. Pick up all preserved internal faces.
1014     for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
1015     {
1016         if (nCellsUsingFace[oldFaceI] == 2)
1017         {
1018             globalFaceMap[oldFaceI] = faceI;
1019             faceMap_[faceI++] = oldFaceI;
1020
1021             // Mark all points from the face
1022             markPoints(oldFaces[oldFaceI], globalPointMap);
1023         }
1024     }
1025
1026     // These are all the internal faces in the mesh.
1027     label nInternalFaces = faceI;
1028
1029     // 2. Boundary faces up to where we want to insert old internal faces
1030     for
1031     (
1032         label oldPatchI = 0;
1033         oldPatchI < oldPatches.size()
1034      && oldPatchI < nextPatchID;
1035         oldPatchI++
1036     )
1037     {
1038         const polyPatch& oldPatch = oldPatches[oldPatchI];
1039
1040         label oldFaceI = oldPatch.start();
1041
1042         forAll (oldPatch, i)
1043         {
1044             if (nCellsUsingFace[oldFaceI] == 1)
1045             {
1046                 // Boundary face is kept.
1047
1048                 // Mark face and increment number of points in set
1049                 globalFaceMap[oldFaceI] = faceI;
1050                 faceMap_[faceI++] = oldFaceI;
1051
1052                 // Mark all points from the face
1053                 markPoints(oldFaces[oldFaceI], globalPointMap);
1054
1055                 // Increment number of patch faces
1056                 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1057             }
1058             oldFaceI++;
1059         }
1060     }
1061
1062     // 3a. old internal faces that have become exposed.
1063     for (label oldFaceI = 0; oldFaceI < oldNInternalFaces; oldFaceI++)
1064     {
1065         if (nCellsUsingFace[oldFaceI] == 1)
1066         {
1067             globalFaceMap[oldFaceI] = faceI;
1068             faceMap_[faceI++] = oldFaceI;
1069
1070             // Mark all points from the face
1071             markPoints(oldFaces[oldFaceI], globalPointMap);
1072
1073             // Increment number of patch faces
1074             boundaryPatchSizes[oldInternalPatchID]++;
1075         }
1076     }
1077
1078     // 3b. coupled patch faces that have become uncoupled.
1079     for
1080     (
1081         label oldFaceI = oldNInternalFaces;
1082         oldFaceI < oldFaces.size();
1083         oldFaceI++
1084     )
1085     {
1086         if (nCellsUsingFace[oldFaceI] == 3)
1087         {
1088             globalFaceMap[oldFaceI] = faceI;
1089             faceMap_[faceI++] = oldFaceI;
1090
1091             // Mark all points from the face
1092             markPoints(oldFaces[oldFaceI], globalPointMap);
1093
1094             // Increment number of patch faces
1095             boundaryPatchSizes[oldInternalPatchID]++;
1096         }
1097     }
1098
1099     // 4. Remaining boundary faces
1100     for
1101     (
1102         label oldPatchI = nextPatchID;
1103         oldPatchI < oldPatches.size();
1104         oldPatchI++
1105     )
1106     {
1107         const polyPatch& oldPatch = oldPatches[oldPatchI];
1108
1109         label oldFaceI = oldPatch.start();
1110
1111         forAll (oldPatch, i)
1112         {
1113             if (nCellsUsingFace[oldFaceI] == 1)
1114             {
1115                 // Boundary face is kept.
1116
1117                 // Mark face and increment number of points in set
1118                 globalFaceMap[oldFaceI] = faceI;
1119                 faceMap_[faceI++] = oldFaceI;
1120
1121                 // Mark all points from the face
1122                 markPoints(oldFaces[oldFaceI], globalPointMap);
1123
1124                 // Increment number of patch faces
1125                 boundaryPatchSizes[globalPatchMap[oldPatchI]]++;
1126             }
1127             oldFaceI++;
1128         }
1129     }
1130
1131     if (faceI != nFacesInSet)
1132     {
1133         FatalErrorIn
1134         (
1135             "fvMeshSubset::setCellSubset(const labelList&"
1136             ", const label, const label, const bool)"
1137         )   << "Problem" << abort(FatalError);
1138     }
1139
1140
1141     // Grab the points map
1142     label nPointsInSet = 0;
1143
1144     forAll(globalPointMap, pointI)
1145     {
1146         if (globalPointMap[pointI] != -1)
1147         {
1148             nPointsInSet++;
1149         }
1150     }
1151     pointMap_.setSize(nPointsInSet);
1152
1153     nPointsInSet = 0;
1154
1155     forAll(globalPointMap, pointI)
1156     {
1157         if (globalPointMap[pointI] != -1)
1158         {
1159             pointMap_[nPointsInSet] = pointI;
1160             globalPointMap[pointI] = nPointsInSet;
1161             nPointsInSet++;
1162         }
1163     }
1164
1165     Pout<< "Number of cells in new mesh : " << cellMap_.size() << endl;
1166     Pout<< "Number of faces in new mesh : " << faceMap_.size() << endl;
1167     Pout<< "Number of points in new mesh: " << pointMap_.size() << endl;
1168
1169     // Make a new mesh
1170     pointField newPoints(pointMap_.size());
1171
1172     label nNewPoints = 0;
1173
1174     forAll (pointMap_, pointI)
1175     {
1176         newPoints[nNewPoints] = oldPoints[pointMap_[pointI]];
1177         nNewPoints++;
1178     }
1179
1180     faceList newFaces(faceMap_.size());
1181
1182     label nNewFaces = 0;
1183
1184     // Make internal faces
1185     for (label faceI = 0; faceI < nInternalFaces; faceI++)
1186     {
1187         const face& oldF = oldFaces[faceMap_[faceI]];
1188
1189         face newF(oldF.size());
1190
1191         forAll (newF, i)
1192         {
1193             newF[i] = globalPointMap[oldF[i]];
1194         }
1195
1196         newFaces[nNewFaces] = newF;
1197         nNewFaces++;
1198     }
1199
1200
1201     // Make boundary faces. (different from internal since might need to be
1202     // flipped)
1203     for (label faceI = nInternalFaces; faceI < faceMap_.size(); faceI++)
1204     {
1205         label oldFaceI = faceMap_[faceI];
1206
1207         face oldF = oldFaces[oldFaceI];
1208
1209         // Turn the faces as necessary to point outwards
1210         if (baseMesh().isInternalFace(oldFaceI))
1211         {
1212             // Was internal face. Possibly turned the wrong way round
1213             if
1214             (
1215                 region[oldOwner[oldFaceI]] != currentRegion
1216              && region[oldNeighbour[oldFaceI]] == currentRegion
1217             )
1218             {
1219                 oldF = oldFaces[oldFaceI].reverseFace();
1220             }
1221         }
1222
1223         // Relabel vertices of the (possibly turned) face.
1224         face newF(oldF.size());
1225
1226         forAll (newF, i)
1227         {
1228             newF[i] = globalPointMap[oldF[i]];
1229         }
1230
1231         newFaces[nNewFaces] = newF;
1232         nNewFaces++;
1233     }
1234
1235
1236
1237     // Create cells
1238     cellList newCells(nCellsInSet);
1239
1240     label nNewCells = 0;
1241
1242     forAll (cellMap_, cellI)
1243     {
1244         const labelList& oldC = oldCells[cellMap_[cellI]];
1245
1246         labelList newC(oldC.size());
1247
1248         forAll (newC, i)
1249         {
1250             newC[i] = globalFaceMap[oldC[i]];
1251         }
1252
1253         newCells[nNewCells] = cell(newC);
1254         nNewCells++;
1255     }
1256
1257
1258     // Make a new mesh
1259     makeFvDictionaries();
1260
1261     fvMeshSubsetPtr_ = new fvMesh
1262     (
1263         IOobject
1264         (
1265             this->name() + "Subset",
1266             baseMesh().time().timeName(),
1267             baseMesh().time(),
1268             IOobject::NO_READ,
1269             IOobject::NO_WRITE
1270         ),
1271         xferMove(newPoints),
1272         xferMove(newFaces),
1273         xferMove(newCells),
1274         syncPar           // parallel synchronisation
1275     );
1276
1277     // Clear point mesh
1278     deleteDemandDrivenData(pointMeshSubsetPtr_);
1279
1280     // Add old patches
1281     List<polyPatch*> newBoundary(nbSize);
1282     patchMap_.setSize(nbSize);
1283     label nNewPatches = 0;
1284     label patchStart = nInternalFaces;
1285
1286
1287     // For parallel: only remove patch if none of the processors has it.
1288     // This only gets done for patches before the one being inserted
1289     // (so patches < nextPatchID)
1290
1291     // Get sum of patch sizes. Zero if patch can be deleted.
1292     labelList globalPatchSizes(boundaryPatchSizes);
1293     globalPatchSizes.setSize(nextPatchID);
1294
1295     if (syncPar && Pstream::parRun())
1296     {
1297         // Get patch names (up to nextPatchID)
1298         List<wordList> patchNames(Pstream::nProcs());
1299         patchNames[Pstream::myProcNo()] = oldPatches.names();
1300         patchNames[Pstream::myProcNo()].setSize(nextPatchID);
1301         Pstream::gatherList(patchNames);
1302         Pstream::scatterList(patchNames);
1303
1304         // Get patch sizes (up to nextPatchID).
1305         // Note that up to nextPatchID the globalPatchMap is an identity so
1306         // no need to index through that.
1307         Pstream::listCombineGather(globalPatchSizes, plusEqOp<label>());
1308         Pstream::listCombineScatter(globalPatchSizes);
1309
1310         // Now all processors have all the patchnames.
1311         // Decide: if all processors have the same patch names and size is zero
1312         // everywhere remove the patch.
1313         bool samePatches = true;
1314
1315         for (label procI = 1; procI < patchNames.size(); procI++)
1316         {
1317             if (patchNames[procI] != patchNames[0])
1318             {
1319                 samePatches = false;
1320                 break;
1321             }
1322         }
1323
1324         if (!samePatches)
1325         {
1326             // Patchnames not sync on all processors so disable removal of
1327             // zero sized patches.
1328             globalPatchSizes = labelMax;
1329         }
1330     }
1331
1332
1333     // Old patches
1334
1335     for
1336     (
1337         label oldPatchI = 0;
1338         oldPatchI < oldPatches.size()
1339      && oldPatchI < nextPatchID;
1340         oldPatchI++
1341     )
1342     {
1343         label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1344
1345         // Clone (even if 0 size)
1346         newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1347         (
1348             fvMeshSubsetPtr_->boundaryMesh(),
1349             nNewPatches,
1350             newSize,
1351             patchStart
1352         ).ptr();
1353
1354         patchStart += newSize;
1355         patchMap_[nNewPatches] = oldPatchI;    // compact patchMap
1356         nNewPatches++;
1357     }
1358
1359     // Inserted patch
1360
1361     if (wantedPatchID == -1)
1362     {
1363         label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
1364
1365         if (syncPar)
1366         {
1367             reduce(oldInternalSize, sumOp<label>());
1368         }
1369
1370         // Newly created patch so is at end. Check if any faces in it.
1371         if (oldInternalSize > 0)
1372         {
1373             newBoundary[nNewPatches] = new polyPatch
1374             (
1375                 "oldInternalFaces",
1376                 boundaryPatchSizes[oldInternalPatchID],
1377                 patchStart,
1378                 nNewPatches,
1379                 fvMeshSubsetPtr_->boundaryMesh()
1380             );
1381
1382             Pout<< "    oldInternalFaces : "
1383                 << boundaryPatchSizes[oldInternalPatchID] << endl;
1384
1385             // The index for the first patch is -1 as it originates from
1386             // the internal faces
1387             patchStart += boundaryPatchSizes[oldInternalPatchID];
1388             patchMap_[nNewPatches] = -1;
1389             nNewPatches++;
1390         }
1391     }
1392
1393     // Old patches
1394
1395     for
1396     (
1397         label oldPatchI = nextPatchID;
1398         oldPatchI < oldPatches.size();
1399         oldPatchI++
1400     )
1401     {
1402         label newSize = boundaryPatchSizes[globalPatchMap[oldPatchI]];
1403
1404         // Patch still exists. Add it
1405         newBoundary[nNewPatches] = oldPatches[oldPatchI].clone
1406         (
1407             fvMeshSubsetPtr_->boundaryMesh(),
1408             nNewPatches,
1409             newSize,
1410             patchStart
1411         ).ptr();
1412
1413         //Pout<< "    " << oldPatches[oldPatchI].name() << " : "
1414         //    << newSize << endl;
1415
1416         patchStart += newSize;
1417         patchMap_[nNewPatches] = oldPatchI;    // compact patchMap
1418         nNewPatches++;
1419     }
1420
1421
1422     // Reset the patch lists
1423     newBoundary.setSize(nNewPatches);
1424     patchMap_.setSize(nNewPatches);
1425
1426
1427     // Add the fvPatches
1428     fvMeshSubsetPtr_->addFvPatches(newBoundary, syncPar);
1429
1430     // Subset and add any zones
1431     subsetZones();
1432 }
1433
1434
1435 void Foam::fvMeshSubset::setLargeCellSubset
1436 (
1437     const labelHashSet& globalCellMap,
1438     const label patchID,
1439     const bool syncPar
1440 )
1441 {
1442     labelList region(baseMesh().nCells(), 0);
1443
1444     forAllConstIter (labelHashSet, globalCellMap, iter)
1445     {
1446         region[iter.key()] = 1;
1447     }
1448     setLargeCellSubset(region, 1, patchID, syncPar);
1449 }
1450
1451
1452 const fvMesh& Foam::fvMeshSubset::subMesh() const
1453 {
1454     checkCellSubset();
1455
1456     return *fvMeshSubsetPtr_;
1457 }
1458
1459
1460 fvMesh& Foam::fvMeshSubset::subMesh()
1461 {
1462     checkCellSubset();
1463
1464     return *fvMeshSubsetPtr_;
1465 }
1466
1467
1468 const pointMesh& Foam::fvMeshSubset::subPointMesh() const
1469 {
1470     if (!pointMeshSubsetPtr_)
1471     {
1472         pointMeshSubsetPtr_ = new pointMesh(subMesh());
1473     }
1474
1475     return *pointMeshSubsetPtr_;
1476 }
1477
1478
1479 pointMesh& Foam::fvMeshSubset::subPointMesh()
1480 {
1481     if (!pointMeshSubsetPtr_)
1482     {
1483         pointMeshSubsetPtr_ = new pointMesh(subMesh());
1484     }
1485
1486     return *pointMeshSubsetPtr_;
1487 }
1488
1489
1490 const labelList& Foam::fvMeshSubset::pointMap() const
1491 {
1492     checkCellSubset();
1493
1494     return pointMap_;
1495 }
1496
1497
1498 const labelList& Foam::fvMeshSubset::faceMap() const
1499 {
1500     checkCellSubset();
1501
1502     return faceMap_;
1503 }
1504
1505
1506 const labelList& Foam::fvMeshSubset::cellMap() const
1507 {
1508     checkCellSubset();
1509
1510     return cellMap_;
1511 }
1512
1513
1514 const labelList& Foam::fvMeshSubset::patchMap() const
1515 {
1516     checkCellSubset();
1517
1518     return patchMap_;
1519 }
1520
1521
1522 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1523
1524 } // End namespace Foam
1525
1526 // ************************************************************************* //