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