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