src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "removeFaces.H"
  28 #include "polyMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "meshTools.H"
  31 #include "polyModifyFace.H"
  32 #include "polyRemoveFace.H"
  33 #include "polyRemoveCell.H"
  34 #include "polyRemovePoint.H"
  35 #include "syncTools.H"
  36 #include "OFstream.H"
  37 #include "indirectPrimitivePatch.H"
  38 #include "Time.H"
  39 #include "faceSet.H"
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45
  46 defineTypeNameAndDebug(removeFaces, 0);
  47
  48 }
  49
  50 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  51
  52 // Changes region of connected set of cells. Can be recursive since hopefully
  53 // only small area of faces removed in one go.
  54 void Foam::removeFaces::changeCellRegion
  55 (
  56     const label cellI,
  57     const label oldRegion,
  58     const label newRegion,
  59     labelList& cellRegion
  60 ) const
  61 {
  62     if (cellRegion[cellI] == oldRegion)
  63     {
  64         cellRegion[cellI] = newRegion;
  65
  66         // Step to neighbouring cells
  67
  68         const labelList& cCells = mesh_.cellCells()[cellI];
  69
  70         forAll(cCells, i)
  71         {
  72             changeCellRegion(cCells[i], oldRegion, newRegion, cellRegion);
  73         }
  74     }
  75 }
  76
  77
  78 // Changes region of connected set of faces. Returns number of changed faces.
  79 Foam::label Foam::removeFaces::changeFaceRegion
  80 (
  81     const labelList& cellRegion,
  82     const boolList& removedFace,
  83     const labelList& nFacesPerEdge,
  84     const label faceI,
  85     const label newRegion,
  86     const labelList& fEdges,
  87     labelList& faceRegion
  88 ) const
  89 {
  90     label nChanged = 0;
  91
  92     if (faceRegion[faceI] == -1 && !removedFace[faceI])
  93     {
  94         faceRegion[faceI] = newRegion;
  95
  96         nChanged = 1;
  97
  98         // Storage for on-the-fly addressing
  99         DynamicList<label> fe;
 100         DynamicList<label> ef;
 101
 102         // Step to neighbouring faces across edges that will get removed
 103         forAll(fEdges, i)
 104         {
 105             label edgeI = fEdges[i];
 106
 107             if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
 108             {
 109                 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 110
 111                 forAll(eFaces, j)
 112                 {
 113                     label nbrFaceI = eFaces[j];
 114
 115                     const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
 116
 117                     nChanged += changeFaceRegion
 118                     (
 119                         cellRegion,
 120                         removedFace,
 121                         nFacesPerEdge,
 122                         nbrFaceI,
 123                         newRegion,
 124                         fEdges1,
 125                         faceRegion
 126                     );
 127                 }
 128             }
 129         }
 130     }
 131     return nChanged;
 132 }
 133
 134
 135 // Mark all faces affected in any way by
 136 // - removal of cells
 137 // - removal of faces
 138 // - removal of edges
 139 // - removal of points
 140 Foam::boolList Foam::removeFaces::getFacesAffected
 141 (
 142     const labelList& cellRegion,
 143     const labelList& cellRegionMaster,
 144     const labelList& facesToRemove,
 145     const labelHashSet& edgesToRemove,
 146     const labelHashSet& pointsToRemove
 147 ) const
 148 {
 149     boolList affectedFace(mesh_.nFaces(), false);
 150
 151     // Mark faces affected by removal of cells
 152     forAll(cellRegion, cellI)
 153     {
 154         label region = cellRegion[cellI];
 155
 156         if (region != -1 && (cellI != cellRegionMaster[region]))
 157         {
 158             const labelList& cFaces = mesh_.cells()[cellI];
 159
 160             forAll(cFaces, cFaceI)
 161             {
 162                 affectedFace[cFaces[cFaceI]] = true;
 163             }
 164         }
 165     }
 166
 167     // Mark faces affected by removal of face.
 168     forAll(facesToRemove, i)
 169     {
 170          affectedFace[facesToRemove[i]] = true;
 171     }
 172
 173     //  Mark faces affected by removal of edges
 174     forAllConstIter(labelHashSet, edgesToRemove, iter)
 175     {
 176         const labelList& eFaces = mesh_.edgeFaces(iter.key());
 177
 178         forAll(eFaces, eFaceI)
 179         {
 180             affectedFace[eFaces[eFaceI]] = true;
 181         }
 182     }
 183
 184     // Mark faces affected by removal of points
 185     forAllConstIter(labelHashSet, pointsToRemove, iter)
 186     {
 187         label pointI = iter.key();
 188
 189         const labelList& pFaces = mesh_.pointFaces()[pointI];
 190
 191         forAll(pFaces, pFaceI)
 192         {
 193             affectedFace[pFaces[pFaceI]] = true;
 194         }
 195     }
 196     return affectedFace;
 197 }
 198
 199
 200 void Foam::removeFaces::writeOBJ
 201 (
 202     const indirectPrimitivePatch& fp,
 203     const fileName& fName
 204 )
 205 {
 206     OFstream str(fName);
 207     Pout<< "removeFaces::writeOBJ : Writing faces to file "
 208         << str.name() << endl;
 209
 210     const pointField& localPoints = fp.localPoints();
 211
 212     forAll(localPoints, i)
 213     {
 214         meshTools::writeOBJ(str, localPoints[i]);
 215     }
 216
 217     const faceList& localFaces = fp.localFaces();
 218
 219     forAll(localFaces, i)
 220     {
 221         const face& f = localFaces[i];
 222
 223         str<< 'f';
 224
 225         forAll(f, fp)
 226         {
 227             str<< ' ' << f[fp]+1;
 228         }
 229         str<< nl;
 230     }
 231 }
 232
 233
 234 // Inserts commands to merge faceLabels into one face.
 235 void Foam::removeFaces::mergeFaces
 236 (
 237     const labelList& cellRegion,
 238     const labelList& cellRegionMaster,
 239     const labelHashSet& pointsToRemove,
 240     const labelList& faceLabels,
 241     polyTopoChange& meshMod
 242 ) const
 243 {
 244     // Construct addressing engine from faceLabels (in order of faceLabels as
 245     // well)
 246     indirectPrimitivePatch fp
 247     (
 248         IndirectList<face>
 249         (
 250             mesh_.faces(),
 251             faceLabels
 252         ),
 253         mesh_.points()
 254     );
 255
 256     // Get outside vertices (in local vertex numbering)
 257
 258     if (fp.edgeLoops().size() != 1)
 259     {
 260         writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
 261         FatalErrorIn("removeFaces::mergeFaces")
 262             << "Cannot merge faces " << faceLabels
 263             << " into single face since outside vertices " << fp.edgeLoops()
 264             << " do not form single loop but form " << fp.edgeLoops().size()
 265             << " loops instead." << abort(FatalError);
 266     }
 267
 268     const labelList& edgeLoop = fp.edgeLoops()[0];
 269
 270     // Get outside vertices in order of one of the faces in faceLabels.
 271     // (this becomes the master face)
 272     // Find the first face that uses edgeLoop[0] and edgeLoop[1] as consecutive
 273     // vertices.
 274
 275     label masterIndex = -1;
 276     bool reverseLoop = false;
 277
 278     const labelList& pFaces = fp.pointFaces()[edgeLoop[0]];
 279
 280     // Find face among pFaces which uses edgeLoop[1]
 281     forAll(pFaces, i)
 282     {
 283         label faceI = pFaces[i];
 284
 285         const face& f = fp.localFaces()[faceI];
 286
 287         label index1 = findIndex(f, edgeLoop[1]);
 288
 289         if (index1 != -1)
 290         {
 291             // Check whether consecutive to edgeLoop[0]
 292             label index0 = findIndex(f, edgeLoop[0]);
 293
 294             if (index0 != -1)
 295             {
 296                 if (index1 == f.fcIndex(index0))
 297                 {
 298                     masterIndex = faceI;
 299                     reverseLoop = false;
 300                     break;
 301                 }
 302                 else if (index1 == f.rcIndex(index0))
 303                 {
 304                     masterIndex = faceI;
 305                     reverseLoop = true;
 306                     break;
 307                 }
 308             }
 309         }
 310     }
 311
 312     if (masterIndex == -1)
 313     {
 314         writeOBJ(fp, mesh_.time().path()/"facesToBeMerged.obj");
 315         FatalErrorIn("removeFaces::mergeFaces")
 316             << "Problem" << abort(FatalError);
 317     }
 318
 319
 320     // Modify the master face.
 321     // ~~~~~~~~~~~~~~~~~~~~~~~
 322
 323     // Modify first face.
 324     label faceI = faceLabels[masterIndex];
 325
 326     label own = mesh_.faceOwner()[faceI];
 327
 328     if (cellRegion[own] != -1)
 329     {
 330         own = cellRegionMaster[cellRegion[own]];
 331     }
 332
 333     label patchID, zoneID, zoneFlip;
 334
 335     getFaceInfo(faceI, patchID, zoneID, zoneFlip);
 336
 337     label nei = -1;
 338
 339     if (mesh_.isInternalFace(faceI))
 340     {
 341         nei = mesh_.faceNeighbour()[faceI];
 342
 343         if (cellRegion[nei] != -1)
 344         {
 345             nei = cellRegionMaster[cellRegion[nei]];
 346         }
 347     }
 348
 349
 350     DynamicList<label> faceVerts(edgeLoop.size());
 351
 352     forAll(edgeLoop, i)
 353     {
 354         label pointI = fp.meshPoints()[edgeLoop[i]];
 355
 356         if (pointsToRemove.found(pointI))
 357         {
 358             //Pout<< "**Removing point " << pointI << " from "
 359             //    << edgeLoop << endl;
 360         }
 361         else
 362         {
 363             faceVerts.append(pointI);
 364         }
 365     }
 366
 367     face mergedFace;
 368     mergedFace.transfer(faceVerts);
 369
 370     if (reverseLoop)
 371     {
 372         reverse(mergedFace);
 373     }
 374
 375     //{
 376     //    Pout<< "Modifying masterface " << faceI
 377     //        << " from faces:" << faceLabels
 378     //        << " old verts:" << UIndirectList<face>(mesh_.faces(), faceLabels)
 379     //        << " for new verts:"
 380     //        << mergedFace
 381     //        << " possibly new owner " << own
 382     //        << " or new nei " << nei
 383     //        << endl;
 384     //}
 385
 386     modFace
 387     (
 388         mergedFace,         // modified face
 389         faceI,              // label of face being modified
 390         own,                // owner
 391         nei,                // neighbour
 392         false,              // face flip
 393         patchID,            // patch for face
 394         false,              // remove from zone
 395         zoneID,             // zone for face
 396         zoneFlip,           // face flip in zone
 397
 398         meshMod
 399     );
 400
 401
 402     // Remove all but master face.
 403     forAll(faceLabels, patchFaceI)
 404     {
 405         if (patchFaceI != masterIndex)
 406         {
 407             //Pout<< "Removing face " << faceLabels[patchFaceI] << endl;
 408
 409             meshMod.setAction(polyRemoveFace(faceLabels[patchFaceI], faceI));
 410         }
 411     }
 412 }
 413
 414
 415 // Get patch, zone info for faceI
 416 void Foam::removeFaces::getFaceInfo
 417 (
 418     const label faceI,
 419
 420     label& patchID,
 421     label& zoneID,
 422     label& zoneFlip
 423 ) const
 424 {
 425     patchID = -1;
 426
 427     if (!mesh_.isInternalFace(faceI))
 428     {
 429         patchID = mesh_.boundaryMesh().whichPatch(faceI);
 430     }
 431
 432     zoneID = mesh_.faceZones().whichZone(faceI);
 433
 434     zoneFlip = false;
 435
 436     if (zoneID >= 0)
 437     {
 438         const faceZone& fZone = mesh_.faceZones()[zoneID];
 439
 440         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 441     }
 442 }
 443
 444
 445 // Return face with all pointsToRemove removed.
 446 Foam::face Foam::removeFaces::filterFace
 447 (
 448     const labelHashSet& pointsToRemove,
 449     const label faceI
 450 ) const
 451 {
 452     const face& f = mesh_.faces()[faceI];
 453
 454     labelList newFace(f.size(), -1);
 455
 456     label newFp = 0;
 457
 458     forAll(f, fp)
 459     {
 460         label vertI = f[fp];
 461
 462         if (!pointsToRemove.found(vertI))
 463         {
 464             newFace[newFp++] = vertI;
 465         }
 466     }
 467
 468     newFace.setSize(newFp);
 469
 470     return face(newFace);
 471 }
 472
 473
 474 // Wrapper for meshMod.modifyFace. Reverses face if own>nei.
 475 void Foam::removeFaces::modFace
 476 (
 477     const face& f,
 478     const label masterFaceID,
 479     const label own,
 480     const label nei,
 481     const bool flipFaceFlux,
 482     const label newPatchID,
 483     const bool removeFromZone,
 484     const label zoneID,
 485     const bool zoneFlip,
 486
 487     polyTopoChange& meshMod
 488 ) const
 489 {
 490     if ((nei == -1) || (own < nei))
 491     {
 492 //        if (debug)
 493 //        {
 494 //            Pout<< "ModifyFace (unreversed) :"
 495 //                << "  faceI:" << masterFaceID
 496 //                << "  f:" << f
 497 //                << "  own:" << own
 498 //                << "  nei:" << nei
 499 //                << "  flipFaceFlux:" << flipFaceFlux
 500 //                << "  newPatchID:" << newPatchID
 501 //                << "  removeFromZone:" << removeFromZone
 502 //                << "  zoneID:" << zoneID
 503 //                << "  zoneFlip:" << zoneFlip
 504 //                << endl;
 505 //        }
 506
 507         meshMod.setAction
 508         (
 509             polyModifyFace
 510             (
 511                 f,              // modified face
 512                 masterFaceID,   // label of face being modified
 513                 own,            // owner
 514                 nei,            // neighbour
 515                 flipFaceFlux,   // face flip
 516                 newPatchID,     // patch for face
 517                 removeFromZone, // remove from zone
 518                 zoneID,         // zone for face
 519                 zoneFlip        // face flip in zone
 520             )
 521         );
 522     }
 523     else
 524     {
 525 //        if (debug)
 526 //        {
 527 //            Pout<< "ModifyFace (!reversed) :"
 528 //                << "  faceI:" << masterFaceID
 529 //                << "  f:" << f.reverseFace()
 530 //                << "  own:" << nei
 531 //                << "  nei:" << own
 532 //                << "  flipFaceFlux:" << flipFaceFlux
 533 //                << "  newPatchID:" << newPatchID
 534 //                << "  removeFromZone:" << removeFromZone
 535 //                << "  zoneID:" << zoneID
 536 //                << "  zoneFlip:" << zoneFlip
 537 //                << endl;
 538 //        }
 539
 540         meshMod.setAction
 541         (
 542             polyModifyFace
 543             (
 544                 f.reverseFace(),// modified face
 545                 masterFaceID,   // label of face being modified
 546                 nei,            // owner
 547                 own,            // neighbour
 548                 flipFaceFlux,   // face flip
 549                 newPatchID,     // patch for face
 550                 removeFromZone, // remove from zone
 551                 zoneID,         // zone for face
 552                 zoneFlip        // face flip in zone
 553             )
 554         );
 555     }
 556 }
 557
 558
 559 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 560
 561 // Construct from mesh
 562 Foam::removeFaces::removeFaces
 563 (
 564     const polyMesh& mesh,
 565     const scalar minCos
 566 )
 567 :
 568     mesh_(mesh),
 569     minCos_(minCos)
 570 {}
 571
 572
 573 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 574
 575 // Removing face connects cells. This function works out a consistent set of
 576 // cell regions.
 577 // - returns faces to remove. Can be extended with additional faces
 578 //   (if owner would become neighbour)
 579 // - sets cellRegion to -1 or to region number
 580 // - regionMaster contains for every region number a master cell.
 581 Foam::label Foam::removeFaces::compatibleRemoves
 582 (
 583     const labelList& facesToRemove,
 584     labelList& cellRegion,
 585     labelList& regionMaster,
 586     labelList& newFacesToRemove
 587 ) const
 588 {
 589     const labelList& faceOwner = mesh_.faceOwner();
 590     const labelList& faceNeighbour = mesh_.faceNeighbour();
 591
 592     cellRegion.setSize(mesh_.nCells());
 593     cellRegion = -1;
 594
 595     regionMaster.setSize(mesh_.nCells());
 596     regionMaster = -1;
 597
 598     label nRegions = 0;
 599
 600     forAll(facesToRemove, i)
 601     {
 602         label faceI = facesToRemove[i];
 603
 604         if (!mesh_.isInternalFace(faceI))
 605         {
 606             FatalErrorIn
 607             (
 608                 "removeFaces::compatibleRemoves(const labelList&"
 609                 ", labelList&, labelList&, labelList&)"
 610             )   << "Not internal face:" << faceI << abort(FatalError);
 611         }
 612
 613
 614         label own = faceOwner[faceI];
 615         label nei = faceNeighbour[faceI];
 616
 617         label region0 = cellRegion[own];
 618         label region1 = cellRegion[nei];
 619
 620         if (region0 == -1)
 621         {
 622             if (region1 == -1)
 623             {
 624                 // Create new region
 625                 cellRegion[own] = nRegions;
 626                 cellRegion[nei] = nRegions;
 627
 628                 // Make owner (lowest numbered!) the master of the region
 629                 regionMaster[nRegions] = own;
 630                 nRegions++;
 631             }
 632             else
 633             {
 634                 // Add owner to neighbour region
 635                 cellRegion[own] = region1;
 636                 // See if owner becomes the master of the region
 637                 regionMaster[region1] = min(own, regionMaster[region1]);
 638             }
 639         }
 640         else
 641         {
 642             if (region1 == -1)
 643             {
 644                 // Add neighbour to owner region
 645                 cellRegion[nei] = region0;
 646                 // nei is higher numbered than own so guaranteed not lower
 647                 // than master of region0.
 648             }
 649             else if (region0 != region1)
 650             {
 651                 // Both have regions. Keep lowest numbered region and master.
 652                 label freedRegion = -1;
 653                 label keptRegion = -1;
 654
 655                 if (region0 < region1)
 656                 {
 657                     changeCellRegion
 658                     (
 659                         nei,
 660                         region1,    // old region
 661                         region0,    // new region
 662                         cellRegion
 663                     );
 664
 665                     keptRegion = region0;
 666                     freedRegion = region1;
 667                 }
 668                 else if (region1 < region0)
 669                 {
 670                     changeCellRegion
 671                     (
 672                         own,
 673                         region0,    // old region
 674                         region1,    // new region
 675                         cellRegion
 676                     );
 677
 678                     keptRegion = region1;
 679                     freedRegion = region0;
 680                 }
 681
 682                 label master0 = regionMaster[region0];
 683                 label master1 = regionMaster[region1];
 684
 685                 regionMaster[freedRegion] = -1;
 686                 regionMaster[keptRegion] = min(master0, master1);
 687             }
 688         }
 689     }
 690
 691     regionMaster.setSize(nRegions);
 692
 693
 694     // Various checks
 695     // - master is lowest numbered in any region
 696     // - regions have more than 1 cell
 697     {
 698         labelList nCells(regionMaster.size(), 0);
 699
 700         forAll(cellRegion, cellI)
 701         {
 702             label r = cellRegion[cellI];
 703
 704             if (r != -1)
 705             {
 706                 nCells[r]++;
 707
 708                 if (cellI < regionMaster[r])
 709                 {
 710                     FatalErrorIn
 711                     (
 712                         "removeFaces::compatibleRemoves(const labelList&"
 713                         ", labelList&, labelList&, labelList&)"
 714                     )   << "Not lowest numbered : cell:" << cellI
 715                         << " region:" << r
 716                         << " regionmaster:" << regionMaster[r]
 717                         << abort(FatalError);
 718                 }
 719             }
 720         }
 721
 722         forAll(nCells, region)
 723         {
 724             if (nCells[region] == 1)
 725             {
 726                 FatalErrorIn
 727                 (
 728                     "removeFaces::compatibleRemoves(const labelList&"
 729                     ", labelList&, labelList&, labelList&)"
 730                 )   << "Region " << region
 731                     << " has only " << nCells[region] << " cells in it"
 732                     << abort(FatalError);
 733             }
 734         }
 735     }
 736
 737
 738     // Count number of used regions
 739     label nUsedRegions = 0;
 740
 741     forAll(regionMaster, i)
 742     {
 743         if (regionMaster[i] != -1)
 744         {
 745             nUsedRegions++;
 746         }
 747     }
 748
 749     // Recreate facesToRemove to be consistent with the cellRegions.
 750     DynamicList<label> allFacesToRemove(facesToRemove.size());
 751
 752     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 753     {
 754         label own = faceOwner[faceI];
 755         label nei = faceNeighbour[faceI];
 756
 757         if (cellRegion[own] != -1 && cellRegion[own] == cellRegion[nei])
 758         {
 759             // Both will become the same cell so add face to list of faces
 760             // to be removed.
 761             allFacesToRemove.append(faceI);
 762         }
 763     }
 764
 765     newFacesToRemove.transfer(allFacesToRemove);
 766
 767     return nUsedRegions;
 768 }
 769
 770
 771 void Foam::removeFaces::setRefinement
 772 (
 773     const labelList& faceLabels,
 774     const labelList& cellRegion,
 775     const labelList& cellRegionMaster,
 776     polyTopoChange& meshMod
 777 ) const
 778 {
 779     if (debug)
 780     {
 781         faceSet facesToRemove(mesh_, "facesToRemove", faceLabels);
 782         Pout<< "Writing faces to remove to faceSet " << facesToRemove.name()
 783             << endl;
 784         facesToRemove.write();
 785     }
 786
 787     // Make map of all faces to be removed
 788     boolList removedFace(mesh_.nFaces(), false);
 789
 790     forAll(faceLabels, i)
 791     {
 792         label faceI = faceLabels[i];
 793
 794         if (!mesh_.isInternalFace(faceI))
 795         {
 796             FatalErrorIn
 797             (
 798                 "removeFaces::setRefinement(const labelList&"
 799                 ", const labelList&, const labelList&, polyTopoChange&)"
 800             )   << "Face to remove is not internal face:" << faceI
 801                 << abort(FatalError);
 802         }
 803
 804         removedFace[faceI] = true;
 805     }
 806
 807
 808     // Edges to be removed
 809     // ~~~~~~~~~~~~~~~~~~~
 810
 811
 812     // Edges to remove
 813     labelHashSet edgesToRemove(faceLabels.size());
 814
 815     // Per face the region it is in. -1 for removed faces, -2 for regions
 816     // consisting of single face only.
 817     labelList faceRegion(mesh_.nFaces(), -1);
 818
 819     // Number of connected face regions
 820     label nRegions = 0;
 821
 822     // Storage for on-the-fly addressing
 823     DynamicList<label> fe;
 824     DynamicList<label> ef;
 825
 826
 827     {
 828         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 829
 830         // Usage of edges by non-removed faces.
 831         // See below about initialization.
 832         labelList nFacesPerEdge(mesh_.nEdges(), -1);
 833
 834         // Count usage of edges by non-removed faces.
 835         forAll(faceLabels, i)
 836         {
 837             label faceI = faceLabels[i];
 838
 839             const labelList& fEdges = mesh_.faceEdges(faceI, fe);
 840
 841             forAll(fEdges, i)
 842             {
 843                 label edgeI = fEdges[i];
 844
 845                 if (nFacesPerEdge[edgeI] == -1)
 846                 {
 847                     nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
 848                 }
 849                 else
 850                 {
 851                     nFacesPerEdge[edgeI]--;
 852                 }
 853             }
 854         }
 855
 856         // Count usage for edges not on faces-to-be-removed.
 857         // Note that this only needs to be done for possibly coupled edges
 858         // so we could choose to loop only over boundary faces and use faceEdges
 859         // of those.
 860
 861         forAll(mesh_.edges(), edgeI)
 862         {
 863             if (nFacesPerEdge[edgeI] == -1)
 864             {
 865                 // Edge not yet handled in loop above so is not used by any
 866                 // face to be removed.
 867
 868                 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 869
 870                 if (eFaces.size() > 2)
 871                 {
 872                     nFacesPerEdge[edgeI] = eFaces.size();
 873                 }
 874                 else if (eFaces.size() == 2)
 875                 {
 876                     // nFacesPerEdge already -1 so do nothing.
 877                 }
 878                 else
 879                 {
 880                     const edge& e = mesh_.edges()[edgeI];
 881
 882                     FatalErrorIn("removeFaces::setRefinement")
 883                         << "Problem : edge has too few face neighbours:"
 884                         << eFaces << endl
 885                         << "edge:" << edgeI
 886                         << " vertices:" << e
 887                         << " coords:" << mesh_.points()[e[0]]
 888                         << mesh_.points()[e[1]]
 889                         << abort(FatalError);
 890                 }
 891             }
 892         }
 893
 894
 895
 896         if (debug)
 897         {
 898             OFstream str(mesh_.time().path()/"edgesWithTwoFaces.obj");
 899             Pout<< "Dumping edgesWithTwoFaces to " << str.name() << endl;
 900             label vertI = 0;
 901
 902             forAll(nFacesPerEdge, edgeI)
 903             {
 904                 if (nFacesPerEdge[edgeI] == 2)
 905                 {
 906                     // Edge will get removed.
 907                     const edge& e = mesh_.edges()[edgeI];
 908
 909                     meshTools::writeOBJ(str, mesh_.points()[e[0]]);
 910                     vertI++;
 911                     meshTools::writeOBJ(str, mesh_.points()[e[1]]);
 912                     vertI++;
 913                     str<< "l " << vertI-1 << ' ' << vertI << nl;
 914                 }
 915             }
 916         }
 917
 918
 919         // Now all unaffected edges will have labelMax, all affected edges the
 920         // number of unremoved faces.
 921
 922         // Filter for edges inbetween two remaining boundary faces that
 923         // make too big an angle.
 924         forAll(nFacesPerEdge, edgeI)
 925         {
 926             if (nFacesPerEdge[edgeI] == 2)
 927             {
 928                 // See if they are two boundary faces
 929                 label f0 = -1;
 930                 label f1 = -1;
 931
 932                 const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 933
 934                 forAll(eFaces, i)
 935                 {
 936                     label faceI = eFaces[i];
 937
 938                     if (!removedFace[faceI] && !mesh_.isInternalFace(faceI))
 939                     {
 940                         if (f0 == -1)
 941                         {
 942                             f0 = faceI;
 943                         }
 944                         else
 945                         {
 946                             f1 = faceI;
 947                             break;
 948                         }
 949                     }
 950                 }
 951
 952                 if (f0 != -1 && f1 != -1)
 953                 {
 954                     // Edge has two boundary faces remaining.
 955                     // See if should be merged.
 956
 957                     label patch0 = patches.whichPatch(f0);
 958                     label patch1 = patches.whichPatch(f1);
 959
 960                     if (patch0 != patch1)
 961                     {
 962                         // Different patches. Do not merge edge.
 963                         WarningIn("removeFaces::setRefinement")
 964                             << "not merging faces " << f0 << " and "
 965                             << f1 << " across patch boundary edge " << edgeI
 966                             << endl;
 967
 968                         // Mark so it gets preserved
 969                         nFacesPerEdge[edgeI] = 3;
 970                     }
 971                     else if (minCos_ < 1 && minCos_ > -1)
 972                     {
 973                         const polyPatch& pp0 = patches[patch0];
 974                         const vectorField& n0 = pp0.faceNormals();
 975
 976                         if
 977                         (
 978                             mag
 979                             (
 980                                 n0[f0 - pp0.start()]
 981                               & n0[f1 - pp0.start()]
 982                             )
 983                             < minCos_
 984                         )
 985                         {
 986                             WarningIn("removeFaces::setRefinement")
 987                                 << "not merging faces " << f0 << " and "
 988                                 << f1 << " across edge " << edgeI
 989                                 << endl;
 990
 991                             // Angle between two remaining faces too large.
 992                             // Mark so it gets preserved
 993                             nFacesPerEdge[edgeI] = 3;
 994                         }
 995                     }
 996                 }
 997                 else if (f0 != -1 || f1 != -1)
 998                 {
 999                     const edge& e = mesh_.edges()[edgeI];
1000
1001                     // Only found one boundary face. Problem.
1002                     FatalErrorIn("removeFaces::setRefinement")
1003                         << "Problem : edge would have one boundary face"
1004                         << " and one internal face using it." << endl
1005                         << "Your remove pattern is probably incorrect." << endl
1006                         << "edge:" << edgeI
1007                         << " nFaces:" << nFacesPerEdge[edgeI]
1008                         << " vertices:" << e
1009                         << " coords:" << mesh_.points()[e[0]]
1010                         << mesh_.points()[e[1]]
1011                         << " face0:" << f0
1012                         << " face1:" << f1
1013                         << abort(FatalError);
1014                 }
1015             }
1016         }
1017
1018
1019
1020         // Check locally (before synchronizing) for strangeness
1021         forAll(nFacesPerEdge, edgeI)
1022         {
1023             if (nFacesPerEdge[edgeI] == 1)
1024             {
1025                 const edge& e = mesh_.edges()[edgeI];
1026
1027                 FatalErrorIn("removeFaces::setRefinement")
1028                     << "Problem : edge would get 1 face using it only"
1029                     << " edge:" << edgeI
1030                     << " nFaces:" << nFacesPerEdge[edgeI]
1031                     << " vertices:" << e
1032                     << " coords:" << mesh_.points()[e[0]]
1033                     << ' ' << mesh_.points()[e[1]]
1034                     << abort(FatalError);
1035             }
1036             // Could check here for boundary edge with <=1 faces remaining.
1037         }
1038
1039
1040         // Synchronize edge usage. This is to make sure that both sides remove
1041         // (or not remove) an edge on the boundary at the same time.
1042         //
1043         // Coupled edges (edge0, edge1 are opposite each other)
1044         // a. edge not on face to be removed, edge has >= 3 faces
1045         // b.  ,,                             edge has 2 faces
1046         // c. edge has >= 3 remaining faces
1047         // d. edge has 2 remaining faces (assume angle>minCos already handled)
1048         //
1049         // - a + a: do not remove edge
1050         // - a + b: do not remove edge
1051         // - a + c: do not remove edge
1052         // - a + d: do not remove edge
1053         //
1054         // - b + b: do not remove edge
1055         // - b + c: do not remove edge
1056         // - b + d: remove edge
1057         //
1058         // - c + c: do not remove edge
1059         // - c + d: do not remove edge
1060         // - d + d: remove edge
1061         //
1062         //
1063         // So code situation a. with >= 3
1064         //                   b. with -1
1065         //                   c. with >=3
1066         //                   d. with 2
1067         // then do max and check result.
1068         //
1069         // a+a : max(3,3) = 3. do not remove
1070         // a+b : max(3,-1) = 3. do not remove
1071         // a+c : max(3,3) = 3. do not remove
1072         // a+d : max(3,2) = 3. do not remove
1073         //
1074         // b+b : max(-1,-1) = -1. do not remove
1075         // b+c : max(-1,3) = 3. do not remove
1076         // b+d : max(-1,2) = 2. remove
1077         //
1078         // c+c : max(3,3) = 3. do not remove
1079         // c+d : max(3,2) = 3. do not remove
1080         //
1081         // d+d : max(2,2) = 2. remove
1082
1083         syncTools::syncEdgeList
1084         (
1085             mesh_,
1086             nFacesPerEdge,
1087             maxEqOp<label>(),
1088             labelMin,               // guaranteed to be overridden by maxEqOp
1089             false                   // no separation
1090         );
1091
1092         // Convert to labelHashSet
1093         forAll(nFacesPerEdge, edgeI)
1094         {
1095             if (nFacesPerEdge[edgeI] == 0)
1096             {
1097                 // 0: edge not used anymore.
1098                 edgesToRemove.insert(edgeI);
1099             }
1100             else if (nFacesPerEdge[edgeI] == 1)
1101             {
1102                 // 1: illegal. Tested above.
1103             }
1104             else if (nFacesPerEdge[edgeI] == 2)
1105             {
1106                 // 2: merge faces.
1107                 edgesToRemove.insert(edgeI);
1108             }
1109         }
1110
1111         if (debug)
1112         {
1113             OFstream str(mesh_.time().path()/"edgesToRemove.obj");
1114             Pout<< "Dumping edgesToRemove to " << str.name() << endl;
1115             label vertI = 0;
1116
1117             forAllConstIter(labelHashSet, edgesToRemove, iter)
1118             {
1119                 // Edge will get removed.
1120                 const edge& e = mesh_.edges()[iter.key()];
1121
1122                 meshTools::writeOBJ(str, mesh_.points()[e[0]]);
1123                 vertI++;
1124                 meshTools::writeOBJ(str, mesh_.points()[e[1]]);
1125                 vertI++;
1126                 str<< "l " << vertI-1 << ' ' << vertI << nl;
1127             }
1128         }
1129
1130
1131         // Walk to fill faceRegion with faces that will be connected across
1132         // edges that will be removed.
1133
1134         label startFaceI = 0;
1135
1136         while (true)
1137         {
1138             // Find unset region.
1139             for (; startFaceI < mesh_.nFaces(); startFaceI++)
1140             {
1141                 if (faceRegion[startFaceI] == -1 && !removedFace[startFaceI])
1142                 {
1143                     break;
1144                 }
1145             }
1146
1147             if (startFaceI == mesh_.nFaces())
1148             {
1149                 break;
1150             }
1151
1152             // Start walking face-edge-face, crossing edges that will get
1153             // removed. Every thus connected region will get single region
1154             // number.
1155             label nRegion = changeFaceRegion
1156             (
1157                 cellRegion,
1158                 removedFace,
1159                 nFacesPerEdge,
1160                 startFaceI,
1161                 nRegions,
1162                 mesh_.faceEdges(startFaceI, fe),
1163                 faceRegion
1164             );
1165
1166             if (nRegion < 1)
1167             {
1168                 FatalErrorIn("setRefinement") << "Problem" << abort(FatalError);
1169             }
1170             else if (nRegion == 1)
1171             {
1172                 // Reset face to be single region.
1173                 faceRegion[startFaceI] = -2;
1174             }
1175             else
1176             {
1177                 nRegions++;
1178             }
1179         }
1180
1181
1182         // Check we're deciding the same on both sides. Since the regioning
1183         // is done based on nFacesPerEdge (which is synced) this should
1184         // indeed be the case.
1185
1186         labelList nbrFaceRegion(faceRegion);
1187
1188         syncTools::swapFaceList
1189         (
1190             mesh_,
1191             nbrFaceRegion,
1192             false                   // no separation
1193         );
1194
1195         labelList toNbrRegion(nRegions, -1);
1196
1197         for
1198         (
1199             label faceI = mesh_.nInternalFaces();
1200             faceI < mesh_.nFaces();
1201             faceI++
1202         )
1203         {
1204             // Get the neighbouring region.
1205             label nbrRegion = nbrFaceRegion[faceI];
1206             label myRegion = faceRegion[faceI];
1207
1208             if (myRegion <= -1 || nbrRegion <= -1)
1209             {
1210                 if (nbrRegion != myRegion)
1211                 {
1212                     FatalErrorIn("removeFaces::setRefinement")
1213                         << "Inconsistent face region across coupled patches."
1214                         << endl
1215                         << "This side has for faceI:" << faceI
1216                         << " region:" << myRegion << endl
1217                         << "The other side has region:" << nbrRegion
1218                         << endl
1219                         << "(region -1 means face is to be deleted)"
1220                         << abort(FatalError);
1221                 }
1222             }
1223             else if (toNbrRegion[myRegion] == -1)
1224             {
1225                 // First visit of region. Store correspondence.
1226                 toNbrRegion[myRegion] = nbrRegion;
1227             }
1228             else
1229             {
1230                 // Second visit of this region.
1231                 if (toNbrRegion[myRegion] != nbrRegion)
1232                 {
1233                     FatalErrorIn("removeFaces::setRefinement")
1234                         << "Inconsistent face region across coupled patches."
1235                         << endl
1236                         << "This side has for faceI:" << faceI
1237                         << " region:" << myRegion
1238                         << " with coupled neighbouring regions:"
1239                         << toNbrRegion[myRegion] << " and "
1240                         << nbrRegion
1241                         << abort(FatalError);
1242                 }
1243             }
1244         }
1245     }
1246
1247     //if (debug)
1248     //{
1249     //    labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1250     //
1251     //    forAll(regionToFaces, regionI)
1252     //    {
1253     //        Pout<< "    " << regionI << " faces:" << regionToFaces[regionI]
1254     //            << endl;
1255     //    }
1256     //}
1257
1258
1259     // Points to be removed
1260     // ~~~~~~~~~~~~~~~~~~~~
1261
1262     labelHashSet pointsToRemove(4*faceLabels.size());
1263
1264
1265     // Per point count the number of unremoved edges. Store the ones that
1266     // are only used by 2 unremoved edges.
1267     {
1268         // Usage of points by non-removed edges.
1269         labelList nEdgesPerPoint(mesh_.nPoints());
1270
1271         const labelListList& pointEdges = mesh_.pointEdges();
1272
1273         forAll(pointEdges, pointI)
1274         {
1275             nEdgesPerPoint[pointI] = pointEdges[pointI].size();
1276         }
1277
1278         forAllConstIter(labelHashSet, edgesToRemove, iter)
1279         {
1280             // Edge will get removed.
1281             const edge& e = mesh_.edges()[iter.key()];
1282
1283             forAll(e, i)
1284             {
1285                 nEdgesPerPoint[e[i]]--;
1286             }
1287         }
1288
1289         // Check locally (before synchronizing) for strangeness
1290         forAll(nEdgesPerPoint, pointI)
1291         {
1292             if (nEdgesPerPoint[pointI] == 1)
1293             {
1294                 FatalErrorIn("removeFaces::setRefinement")
1295                     << "Problem : point would get 1 edge using it only."
1296                     << " pointI:" << pointI
1297                     << " coord:" << mesh_.points()[pointI]
1298                     << abort(FatalError);
1299             }
1300         }
1301
1302         // Synchronize point usage. This is to make sure that both sides remove
1303         // (or not remove) a point on the boundary at the same time.
1304         syncTools::syncPointList
1305         (
1306             mesh_,
1307             nEdgesPerPoint,
1308             maxEqOp<label>(),
1309             labelMin,
1310             false                   // no separation
1311         );
1312
1313         forAll(nEdgesPerPoint, pointI)
1314         {
1315             if (nEdgesPerPoint[pointI] == 0)
1316             {
1317                 pointsToRemove.insert(pointI);
1318             }
1319             else if (nEdgesPerPoint[pointI] == 1)
1320             {
1321                 // Already checked before
1322             }
1323             else if (nEdgesPerPoint[pointI] == 2)
1324             {
1325                 // Remove point and merge edges.
1326                 pointsToRemove.insert(pointI);
1327             }
1328         }
1329     }
1330
1331
1332     if (debug)
1333     {
1334         OFstream str(mesh_.time().path()/"pointsToRemove.obj");
1335         Pout<< "Dumping pointsToRemove to " << str.name() << endl;
1336
1337         forAllConstIter(labelHashSet, pointsToRemove, iter)
1338         {
1339             meshTools::writeOBJ(str, mesh_.points()[iter.key()]);
1340         }
1341     }
1342
1343
1344     // All faces affected in any way
1345     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1346
1347     // Get all faces affected in any way by removal of points/edges/faces/cells
1348     boolList affectedFace
1349     (
1350         getFacesAffected
1351         (
1352             cellRegion,
1353             cellRegionMaster,
1354             faceLabels,
1355             edgesToRemove,
1356             pointsToRemove
1357         )
1358     );
1359
1360     //
1361     // Now we know
1362     // - faceLabels         : faces to remove (sync since no boundary faces)
1363     // - cellRegion/Master  : cells to remove (sync since cells)
1364     // - pointsToRemove     : points to remove (sync)
1365     // - faceRegion         : connected face region of faces to be merged (sync)
1366     // - affectedFace       : faces with points removed and/or owner/neighbour
1367     //                        changed (non sync)
1368
1369
1370     // Start modifying mesh and keep track of faces changed.
1371
1372
1373     // Do all removals
1374     // ~~~~~~~~~~~~~~~
1375
1376     // Remove split faces.
1377     forAll(faceLabels, labelI)
1378     {
1379         label faceI = faceLabels[labelI];
1380
1381         // Remove face if not yet uptodate (which is never; but want to be
1382         // consistent with rest of face removals/modifications)
1383         if (affectedFace[faceI])
1384         {
1385             affectedFace[faceI] = false;
1386
1387             meshMod.setAction(polyRemoveFace(faceI, -1));
1388         }
1389     }
1390
1391
1392     // Remove points.
1393     forAllConstIter(labelHashSet, pointsToRemove, iter)
1394     {
1395         label pointI = iter.key();
1396
1397         meshMod.setAction(polyRemovePoint(pointI, -1));
1398     }
1399
1400
1401     // Remove cells.
1402     forAll(cellRegion, cellI)
1403     {
1404         label region = cellRegion[cellI];
1405
1406         if (region != -1 && (cellI != cellRegionMaster[region]))
1407         {
1408             meshMod.setAction(polyRemoveCell(cellI, cellRegionMaster[region]));
1409         }
1410     }
1411
1412
1413
1414     // Merge faces across edges to be merged
1415     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1416
1417     // Invert faceRegion so we get region to faces.
1418     {
1419         labelListList regionToFaces(invertOneToMany(nRegions, faceRegion));
1420
1421         forAll(regionToFaces, regionI)
1422         {
1423             const labelList& rFaces = regionToFaces[regionI];
1424
1425             if (rFaces.size() <= 1)
1426             {
1427                 FatalErrorIn("setRefinement")
1428                     << "Region:" << regionI
1429                     << " contains only faces " << rFaces
1430                     << abort(FatalError);
1431             }
1432
1433             // rFaces[0] is master, rest gets removed.
1434             mergeFaces
1435             (
1436                 cellRegion,
1437                 cellRegionMaster,
1438                 pointsToRemove,
1439                 rFaces,
1440                 meshMod
1441             );
1442
1443             forAll(rFaces, i)
1444             {
1445                 affectedFace[rFaces[i]] = false;
1446             }
1447         }
1448     }
1449
1450
1451     // Remaining affected faces
1452     // ~~~~~~~~~~~~~~~~~~~~~~~~
1453
1454     // Check if any remaining faces have not been updated for new slave/master
1455     // or points removed.
1456     forAll(affectedFace, faceI)
1457     {
1458         if (affectedFace[faceI])
1459         {
1460             affectedFace[faceI] = false;
1461
1462             face f(filterFace(pointsToRemove, faceI));
1463
1464             label own = mesh_.faceOwner()[faceI];
1465
1466             if (cellRegion[own] != -1)
1467             {
1468                 own = cellRegionMaster[cellRegion[own]];
1469             }
1470
1471             label patchID, zoneID, zoneFlip;
1472
1473             getFaceInfo(faceI, patchID, zoneID, zoneFlip);
1474
1475             label nei = -1;
1476
1477             if (mesh_.isInternalFace(faceI))
1478             {
1479                 nei = mesh_.faceNeighbour()[faceI];
1480
1481                 if (cellRegion[nei] != -1)
1482                 {
1483                     nei = cellRegionMaster[cellRegion[nei]];
1484                 }
1485             }
1486
1487 //            if (debug)
1488 //            {
1489 //                Pout<< "Modifying " << faceI
1490 //                    << " old verts:" << mesh_.faces()[faceI]
1491 //                    << " for new verts:" << f
1492 //                    << " or for new owner " << own << " or for new nei "
1493 //                    << nei
1494 //                    << endl;
1495 //            }
1496
1497             modFace
1498             (
1499                 f,                  // modified face
1500                 faceI,              // label of face being modified
1501                 own,                // owner
1502                 nei,                // neighbour
1503                 false,              // face flip
1504                 patchID,            // patch for face
1505                 false,              // remove from zone
1506                 zoneID,             // zone for face
1507                 zoneFlip,           // face flip in zone
1508
1509                 meshMod
1510             );
1511         }
1512     }
1513 }
1514
1515
1516 // ************************************************************************* //