src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.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 "combineFaces.H"
  28 #include "polyMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "polyRemoveFace.H"
  31 #include "polyAddFace.H"
  32 #include "polyModifyFace.H"
  33 #include "polyRemovePoint.H"
  34 #include "polyAddPoint.H"
  35 #include "syncTools.H"
  36 #include "meshTools.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 defineTypeNameAndDebug(combineFaces, 0);
  44
  45 }
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // Test face for (almost) convexeness. Allows certain convexity before
  50 // complaining.
  51 // For every two consecutive edges calculate the normal. If it differs too
  52 // much from the average face normal complain.
  53 bool Foam::combineFaces::convexFace
  54 (
  55     const scalar minConcaveCos,
  56     const pointField& points,
  57     const face& f
  58 )
  59 {
  60     // Get outwards pointing normal of f.
  61     vector n = f.normal(points);
  62     n /= mag(n);
  63
  64     // Get edge from f[0] to f[size-1];
  65     vector ePrev(points[f[0]] - points[f[f.size()-1]]);
  66     scalar magEPrev = mag(ePrev);
  67     ePrev /= magEPrev + VSMALL;
  68
  69     forAll(f, fp0)
  70     {
  71         // Get vertex after fp
  72         label fp1 = (fp0 + 1) % f.size();
  73
  74         // Normalized vector between two consecutive points
  75         vector e10(points[f[fp1]] - points[f[fp0]]);
  76         scalar magE10 = mag(e10);
  77         e10 /= magE10 + VSMALL;
  78
  79         if (magEPrev > SMALL && magE10 > SMALL)
  80         {
  81             vector edgeNormal = ePrev ^ e10;
  82
  83             if ((edgeNormal & n) < 0)
  84             {
  85                 // Concave. Check angle.
  86                 if ((ePrev & e10) < minConcaveCos)
  87                 {
  88                     return false;
  89                 }
  90             }
  91         }
  92
  93         ePrev = e10;
  94         magEPrev = magE10;
  95     }
  96
  97     // Not a single internal angle is concave so face is convex.
  98     return true;
  99 }
 100
 101
 102 // Determines if set of faces is valid to collapse into single face.
 103 bool Foam::combineFaces::validFace
 104 (
 105     const scalar minConcaveCos,
 106     const indirectPrimitivePatch& bigFace
 107 )
 108 {
 109     // Get outside vertices (in local vertex numbering)
 110     const labelListList& edgeLoops = bigFace.edgeLoops();
 111
 112     if (edgeLoops.size() > 1)
 113     {
 114         return false;
 115     }
 116
 117     // Check for convexness
 118     face f(getOutsideFace(bigFace));
 119
 120     return convexFace(minConcaveCos, bigFace.points(), f);
 121 }
 122
 123
 124 void Foam::combineFaces::regioniseFaces
 125 (
 126     const scalar minCos,
 127     const label cellI,
 128     Map<label>& faceRegion
 129 ) const
 130 {
 131     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 132     const labelList& cEdges = mesh_.cellEdges()[cellI];
 133
 134     forAll(cEdges, i)
 135     {
 136         label edgeI = cEdges[i];
 137
 138         label f0, f1;
 139         meshTools::getEdgeFaces(mesh_, cellI, edgeI, f0, f1);
 140
 141         label p0 = patches.whichPatch(f0);
 142         label p1 = patches.whichPatch(f1);
 143
 144         // Face can be merged if
 145         // - same non-coupled patch
 146         // - small angle
 147         if (p0 != -1 && p0 == p1 && !patches[p0].coupled())
 148         {
 149             vector f0Normal = mesh_.faceAreas()[f0];
 150             f0Normal /= mag(f0Normal);
 151             vector f1Normal = mesh_.faceAreas()[f1];
 152             f1Normal /= mag(f1Normal);
 153
 154             if ((f0Normal&f1Normal) > minCos)
 155             {
 156                 Map<label>::const_iterator f0Fnd = faceRegion.find(f0);
 157
 158                 label region0 = -1;
 159                 if (f0Fnd != faceRegion.end())
 160                 {
 161                     region0 = f0Fnd();
 162                 }
 163
 164                 Map<label>::const_iterator f1Fnd = faceRegion.find(f1);
 165
 166                 label region1 = -1;
 167                 if (f1Fnd != faceRegion.end())
 168                 {
 169                     region1 = f1Fnd();
 170                 }
 171
 172                 if (region0 == -1)
 173                 {
 174                     if (region1 == -1)
 175                     {
 176                         label useRegion = faceRegion.size();
 177                         faceRegion.insert(f0, useRegion);
 178                         faceRegion.insert(f1, useRegion);
 179                     }
 180                     else
 181                     {
 182                         faceRegion.insert(f0, region1);
 183                     }
 184                 }
 185                 else
 186                 {
 187                     if (region1 == -1)
 188                     {
 189                         faceRegion.insert(f1, region0);
 190                     }
 191                     else if (region0 != region1)
 192                     {
 193                         // Merge the two regions
 194                         label useRegion = min(region0, region1);
 195                         label freeRegion = max(region0, region1);
 196
 197                         forAllIter(Map<label>, faceRegion, iter)
 198                         {
 199                             if (iter() == freeRegion)
 200                             {
 201                                 iter() = useRegion;
 202                             }
 203                         }
 204                     }
 205                 }
 206             }
 207         }
 208     }
 209 }
 210
 211
 212 bool Foam::combineFaces::faceNeighboursValid
 213 (
 214     const label cellI,
 215     const Map<label>& faceRegion
 216 ) const
 217 {
 218     if (faceRegion.size() <= 1)
 219     {
 220         return true;
 221     }
 222
 223     const labelListList& faceEdges = mesh_.faceEdges();
 224     const cell& cFaces = mesh_.cells()[cellI];
 225
 226     // Test for face collapsing to edge since too many neighbours merged.
 227     forAll(cFaces, cFaceI)
 228     {
 229         label faceI = cFaces[cFaceI];
 230
 231         if (!faceRegion.found(faceI))
 232         {
 233             const labelList& fEdges = faceEdges[faceI];
 234
 235             // Count number of remaining faces neighbouring faceI. This has
 236             // to be 3 or more.
 237
 238             // Unregioned neighbouring faces
 239             DynamicList<label> neighbourFaces(cFaces.size());
 240             // Regioned neighbouring faces
 241             labelHashSet neighbourRegions;
 242
 243             forAll(fEdges, i)
 244             {
 245                 label edgeI = fEdges[i];
 246                 label nbrI = meshTools::otherFace(mesh_, cellI, faceI, edgeI);
 247
 248                 Map<label>::const_iterator iter = faceRegion.find(nbrI);
 249
 250                 if (iter == faceRegion.end())
 251                 {
 252                     if (findIndex(neighbourFaces, nbrI) == -1)
 253                     {
 254                         neighbourFaces.append(nbrI);
 255                     }
 256                 }
 257                 else
 258                 {
 259                     neighbourRegions.insert(iter());
 260                 }
 261             }
 262
 263             if ((neighbourFaces.size()+neighbourRegions.size()) < 3)
 264             {
 265                 return false;
 266             }
 267         }
 268     }
 269     return true;
 270 }
 271
 272
 273 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 274
 275 // Construct from mesh
 276 Foam::combineFaces::combineFaces
 277 (
 278     const polyMesh& mesh,
 279     const bool undoable
 280 )
 281 :
 282     mesh_(mesh),
 283     undoable_(undoable),
 284     masterFace_(0),
 285     faceSetsVertices_(0),
 286     savedPointLabels_(0),
 287     savedPoints_(0)
 288 {}
 289
 290
 291 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 292
 293 Foam::labelListList Foam::combineFaces::getMergeSets
 294 (
 295     const scalar featureCos,
 296     const scalar minConcaveCos,
 297     const labelHashSet& boundaryCells
 298 ) const
 299 {
 300     // Lists of faces that can be merged.
 301     DynamicList<labelList> allFaceSets(boundaryCells.size() / 10);
 302
 303     // On all cells regionise the faces
 304     forAllConstIter(labelHashSet, boundaryCells, iter)
 305     {
 306         label cellI = iter.key();
 307
 308         const cell& cFaces = mesh_.cells()[cellI];
 309
 310         // Region per face
 311         Map<label> faceRegion(cFaces.size());
 312         regioniseFaces(featureCos, cellI, faceRegion);
 313
 314         // Now we have in faceRegion for every face the region with planar
 315         // face sharing the same region. We now check whether the resulting
 316         // sets cause a face
 317         // - to become a set of edges since too many faces are merged.
 318         // - to become convex
 319
 320         if (faceNeighboursValid(cellI, faceRegion))
 321         {
 322             // Create region-to-faces addressing
 323             Map<labelList> regionToFaces(faceRegion.size());
 324
 325             forAllConstIter(Map<label>, faceRegion, iter)
 326             {
 327                 label faceI = iter.key();
 328                 label region = iter();
 329
 330                 Map<labelList>::iterator regionFnd = regionToFaces.find(region);
 331
 332                 if (regionFnd != regionToFaces.end())
 333                 {
 334                     labelList& setFaces = regionFnd();
 335                     label sz = setFaces.size();
 336                     setFaces.setSize(sz+1);
 337                     setFaces[sz] = faceI;
 338                 }
 339                 else
 340                 {
 341                     regionToFaces.insert(region, labelList(1, faceI));
 342                 }
 343             }
 344
 345             // For every set check if it forms a valid convex face
 346
 347             forAllConstIter(Map<labelList>, regionToFaces, iter)
 348             {
 349                 // Make face out of setFaces
 350                 indirectPrimitivePatch bigFace
 351                 (
 352                     IndirectList<face>
 353                     (
 354                         mesh_.faces(),
 355                         iter()
 356                     ),
 357                     mesh_.points()
 358                 );
 359
 360                 // Only store if -only one outside loop -which forms convex face
 361                 if (validFace(minConcaveCos, bigFace))
 362                 {
 363                     allFaceSets.append(iter());
 364                 }
 365             }
 366         }
 367     }
 368
 369     return allFaceSets.shrink();
 370 }
 371
 372
 373 Foam::labelListList Foam::combineFaces::getMergeSets
 374 (
 375     const scalar featureCos,
 376     const scalar minConcaveCos
 377 ) const
 378 {
 379     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 380
 381     // Pick up all cells on boundary
 382     labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
 383
 384     forAll(patches, patchI)
 385     {
 386         const polyPatch& patch = patches[patchI];
 387
 388         if (!patch.coupled())
 389         {
 390             forAll(patch, i)
 391             {
 392                 boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
 393             }
 394         }
 395     }
 396
 397     return getMergeSets(featureCos, minConcaveCos, boundaryCells);
 398 }
 399
 400
 401 // Gets outside edgeloop as a face
 402 // - in same order as faces
 403 // - in mesh vertex labels
 404 Foam::face Foam::combineFaces::getOutsideFace
 405 (
 406     const indirectPrimitivePatch& fp
 407 )
 408 {
 409     if (fp.edgeLoops().size() != 1)
 410     {
 411         FatalErrorIn
 412         (
 413             "combineFaces::getOutsideFace(const indirectPrimitivePatch&)"
 414         )   << "Multiple outside loops:" << fp.edgeLoops()
 415             << abort(FatalError);
 416     }
 417
 418     // Get first boundary edge. Since guaranteed one edgeLoop when in here this
 419     // edge must be on it.
 420     label bEdgeI = fp.nInternalEdges();
 421
 422     const edge& e = fp.edges()[bEdgeI];
 423
 424     const labelList& eFaces = fp.edgeFaces()[bEdgeI];
 425
 426     if (eFaces.size() != 1)
 427     {
 428         FatalErrorIn
 429         (
 430             "combineFaces::getOutsideFace(const indirectPrimitivePatch&)"
 431         )   << "boundary edge:" << bEdgeI
 432             << " points:" << fp.meshPoints()[e[0]]
 433             << ' ' << fp.meshPoints()[e[1]]
 434             << " on indirectPrimitivePatch has " << eFaces.size()
 435             << " faces using it" << abort(FatalError);
 436     }
 437
 438
 439     // Outside loop
 440     const labelList& outsideLoop = fp.edgeLoops()[0];
 441
 442
 443     // Get order of edge e in outside loop
 444     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 445
 446     // edgeLoopConsistent : true  = edge in same order as outsideloop
 447     //                      false = opposite order
 448     bool edgeLoopConsistent = false;
 449
 450     {
 451         label index0 = findIndex(outsideLoop, e[0]);
 452         label index1 = findIndex(outsideLoop, e[1]);
 453
 454         if (index0 == -1 || index1 == -1)
 455         {
 456             FatalErrorIn
 457             (
 458                 "combineFaces::getOutsideFace"
 459                 "(const indirectPrimitivePatch&)"
 460             )   << "Cannot find boundary edge:" << e
 461                 << " points:" << fp.meshPoints()[e[0]]
 462                 << ' ' << fp.meshPoints()[e[1]]
 463                 << " in edgeLoop:" << outsideLoop << abort(FatalError);
 464         }
 465         else if (index1 == outsideLoop.fcIndex(index0))
 466         {
 467             edgeLoopConsistent = true;
 468         }
 469         else if (index0 == outsideLoop.fcIndex(index1))
 470         {
 471             edgeLoopConsistent = false;
 472         }
 473         else
 474         {
 475             FatalErrorIn
 476             (
 477                 "combineFaces::getOutsideFace"
 478                 "(const indirectPrimitivePatch&)"
 479             )   << "Cannot find boundary edge:" << e
 480                 << " points:" << fp.meshPoints()[e[0]]
 481                 << ' ' << fp.meshPoints()[e[1]]
 482                 << " on consecutive points in edgeLoop:"
 483                 << outsideLoop << abort(FatalError);
 484         }
 485     }
 486
 487
 488     // Get face in local vertices
 489     const face& localF = fp.localFaces()[eFaces[0]];
 490
 491     // Get order of edge in localF
 492     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 493
 494     // faceEdgeConsistent : true  = edge in same order as localF
 495     //                      false = opposite order
 496     bool faceEdgeConsistent = false;
 497
 498     {
 499         // Find edge in face.
 500         label index = findIndex(fp.faceEdges()[eFaces[0]], bEdgeI);
 501
 502         if (index == -1)
 503         {
 504             FatalErrorIn
 505             (
 506                 "combineFaces::getOutsideFace"
 507                 "(const indirectPrimitivePatch&)"
 508             )   << "Cannot find boundary edge:" << e
 509                 << " points:" << fp.meshPoints()[e[0]]
 510                 << ' ' << fp.meshPoints()[e[1]]
 511                 << " in face:" << eFaces[0]
 512                 << " edges:" << fp.faceEdges()[eFaces[0]]
 513                 << abort(FatalError);
 514         }
 515         else if (localF[index] == e[0] && localF.nextLabel(index) == e[1])
 516         {
 517             faceEdgeConsistent = true;
 518         }
 519         else if (localF[index] == e[1] && localF.nextLabel(index) == e[0])
 520         {
 521             faceEdgeConsistent = false;
 522         }
 523         else
 524         {
 525             FatalErrorIn
 526             (
 527                 "combineFaces::getOutsideFace"
 528                 "(const indirectPrimitivePatch&)"
 529             )   << "Cannot find boundary edge:" << e
 530                 << " points:" << fp.meshPoints()[e[0]]
 531                 << ' ' << fp.meshPoints()[e[1]]
 532                 << " in face:" << eFaces[0] << " verts:" << localF
 533                 << abort(FatalError);
 534         }
 535     }
 536
 537     // Get face in mesh points.
 538     face meshFace(renumber(fp.meshPoints(), outsideLoop));
 539
 540     if (faceEdgeConsistent != edgeLoopConsistent)
 541     {
 542         reverse(meshFace);
 543     }
 544     return meshFace;
 545 }
 546
 547
 548 void Foam::combineFaces::setRefinement
 549 (
 550     const labelListList& faceSets,
 551     polyTopoChange& meshMod
 552 )
 553 {
 554     if (undoable_)
 555     {
 556         masterFace_.setSize(faceSets.size());
 557         faceSetsVertices_.setSize(faceSets.size());
 558         forAll(faceSets, setI)
 559         {
 560             const labelList& setFaces = faceSets[setI];
 561
 562             masterFace_[setI] = setFaces[0];
 563             faceSetsVertices_[setI] = IndirectList<face>
 564             (
 565                 mesh_.faces(),
 566                 setFaces
 567             );
 568         }
 569     }
 570
 571     // Running count of number of faces using a point
 572     labelList nPointFaces(mesh_.nPoints(), 0);
 573
 574     const labelListList& pointFaces = mesh_.pointFaces();
 575
 576     forAll(pointFaces, pointI)
 577     {
 578         nPointFaces[pointI] = pointFaces[pointI].size();
 579     }
 580
 581     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 582
 583     forAll(faceSets, setI)
 584     {
 585         const labelList& setFaces = faceSets[setI];
 586
 587         // Check
 588         if (debug)
 589         {
 590             forAll(setFaces, i)
 591             {
 592                 label patchI = patches.whichPatch(setFaces[i]);
 593
 594                 if (patchI == -1 || patches[patchI].coupled())
 595                 {
 596                     FatalErrorIn
 597                     (
 598                         "combineFaces::setRefinement"
 599                         "(const bool, const labelListList&"
 600                         ", polyTopoChange&)"
 601                     )   << "Can only merge non-coupled boundary faces"
 602                         << " but found internal or coupled face:"
 603                         << setFaces[i] << " in set " << setI
 604                         << abort(FatalError);
 605                 }
 606             }
 607         }
 608
 609         // Make face out of setFaces
 610         indirectPrimitivePatch bigFace
 611         (
 612             IndirectList<face>
 613             (
 614                 mesh_.faces(),
 615                 setFaces
 616             ),
 617             mesh_.points()
 618         );
 619
 620         // Get outside vertices (in local vertex numbering)
 621         const labelListList& edgeLoops = bigFace.edgeLoops();
 622
 623         if (edgeLoops.size() != 1)
 624         {
 625             FatalErrorIn
 626             (
 627                 "combineFaces::setRefinement"
 628                 "(const bool, const labelListList&, polyTopoChange&)"
 629             )   << "Faces to-be-merged " << setFaces
 630                 << " do not form a single big face." << nl
 631                 << abort(FatalError);
 632         }
 633
 634
 635         // Delete all faces except master, whose face gets modified.
 636
 637         // Modify master face
 638         // ~~~~~~~~~~~~~~~~~~
 639
 640         label masterFaceI = setFaces[0];
 641
 642         // Get outside face in mesh vertex labels
 643         face outsideFace(getOutsideFace(bigFace));
 644
 645         label zoneID = mesh_.faceZones().whichZone(masterFaceI);
 646
 647         bool zoneFlip = false;
 648
 649         if (zoneID >= 0)
 650         {
 651             const faceZone& fZone = mesh_.faceZones()[zoneID];
 652
 653             zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
 654         }
 655
 656         label patchI = mesh_.boundaryMesh().whichPatch(masterFaceI);
 657
 658         meshMod.setAction
 659         (
 660             polyModifyFace
 661             (
 662                 outsideFace,                    // modified face
 663                 masterFaceI,                    // label of face being modified
 664                 mesh_.faceOwner()[masterFaceI], // owner
 665                 -1,                             // neighbour
 666                 false,                          // face flip
 667                 patchI,                         // patch for face
 668                 false,                          // remove from zone
 669                 zoneID,                         // zone for face
 670                 zoneFlip                        // face flip in zone
 671             )
 672         );
 673
 674
 675         // Delete all non-master faces
 676         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 677
 678         for (label i = 1; i < setFaces.size(); i++)
 679         {
 680             meshMod.setAction(polyRemoveFace(setFaces[i]));
 681         }
 682
 683
 684         // Mark unused points for deletion
 685         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 686
 687         // Remove count of all faces combined
 688         forAll(setFaces, i)
 689         {
 690             const face& f = mesh_.faces()[setFaces[i]];
 691
 692             forAll(f, fp)
 693             {
 694                 nPointFaces[f[fp]]--;
 695             }
 696         }
 697         // But keep count on new outside face
 698         forAll(outsideFace, fp)
 699         {
 700             nPointFaces[outsideFace[fp]]++;
 701         }
 702     }
 703
 704
 705     // If one or more processors use the point it needs to be kept.
 706     syncTools::syncPointList
 707     (
 708         mesh_,
 709         nPointFaces,
 710         plusEqOp<label>(),
 711         0,                  // null value
 712         false               // no separation
 713     );
 714
 715     // Remove all unused points. Store position if undoable.
 716     if (!undoable_)
 717     {
 718         forAll(nPointFaces, pointI)
 719         {
 720             if (nPointFaces[pointI] == 0)
 721             {
 722                 meshMod.setAction(polyRemovePoint(pointI));
 723             }
 724         }
 725     }
 726     else
 727     {
 728         // Count removed points
 729         label n = 0;
 730         forAll(nPointFaces, pointI)
 731         {
 732             if (nPointFaces[pointI] == 0)
 733             {
 734                 n++;
 735             }
 736         }
 737
 738         savedPointLabels_.setSize(n);
 739         savedPoints_.setSize(n);
 740         Map<label> meshToSaved(2*n);
 741
 742         // Remove points and store position
 743         n = 0;
 744         forAll(nPointFaces, pointI)
 745         {
 746             if (nPointFaces[pointI] == 0)
 747             {
 748                 meshMod.setAction(polyRemovePoint(pointI));
 749
 750                 savedPointLabels_[n] = pointI;
 751                 savedPoints_[n] = mesh_.points()[pointI];
 752
 753                 meshToSaved.insert(pointI, n);
 754                 n++;
 755             }
 756         }
 757
 758         // Update stored vertex labels. Negative indices index into local points
 759         forAll(faceSetsVertices_, setI)
 760         {
 761             faceList& setFaces = faceSetsVertices_[setI];
 762
 763             forAll(setFaces, i)
 764             {
 765                 face& f = setFaces[i];
 766
 767                 forAll(f, fp)
 768                 {
 769                     label pointI = f[fp];
 770
 771                     if (nPointFaces[pointI] == 0)
 772                     {
 773                         f[fp] = -meshToSaved[pointI]-1;
 774                     }
 775                 }
 776             }
 777         }
 778     }
 779 }
 780
 781
 782 void Foam::combineFaces::updateMesh(const mapPolyMesh& map)
 783 {
 784     if (undoable_)
 785     {
 786         // Master face just renumbering of point labels
 787         inplaceRenumber(map.reverseFaceMap(), masterFace_);
 788
 789         // Stored faces refer to backed-up vertices (not changed)
 790         // and normal vertices (need to be renumbered)
 791         forAll(faceSetsVertices_, setI)
 792         {
 793             faceList& faces = faceSetsVertices_[setI];
 794
 795             forAll(faces, i)
 796             {
 797                 // Note: faces[i] can have negative elements.
 798                 face& f = faces[i];
 799
 800                 forAll(f, fp)
 801                 {
 802                     label pointI = f[fp];
 803
 804                     if (pointI >= 0)
 805                     {
 806                         f[fp] = map.reversePointMap()[pointI];
 807
 808                         if (f[fp] < 0)
 809                         {
 810                             FatalErrorIn
 811                             (
 812                                 "combineFaces::updateMesh"
 813                                 "(const mapPolyMesh&)"
 814                             )   << "In set " << setI << " at position " << i
 815                                 << " with master face "
 816                                 << masterFace_[setI] << nl
 817                                 << "the points of the slave face " << faces[i]
 818                                 << " don't exist anymore."
 819                                 << abort(FatalError);
 820                         }
 821                     }
 822                 }
 823             }
 824         }
 825     }
 826 }
 827
 828
 829 // Note that faces on coupled patches are never combined (or at least
 830 // getMergeSets never picks them up. Thus the points on them will never get
 831 // deleted since still used by the coupled faces. So never a risk of undoing
 832 // things (faces/points) on coupled patches.
 833 void Foam::combineFaces::setUnrefinement
 834 (
 835     const labelList& masterFaces,
 836     polyTopoChange& meshMod,
 837     Map<label>& restoredPoints,
 838     Map<label>& restoredFaces,
 839     Map<label>& restoredCells
 840 )
 841 {
 842     if (!undoable_)
 843     {
 844         FatalErrorIn
 845         (
 846             "combineFaces::setUnrefinement"
 847             "(const labelList&, polyTopoChange&"
 848             ", Map<label>&, Map<label>&, Map<label>&)"
 849         )   << "Can only call setUnrefinement if constructed with"
 850             << " unrefinement capability." << exit(FatalError);
 851     }
 852
 853
 854     // Restored points
 855     labelList addedPoints(savedPoints_.size(), -1);
 856
 857     // Invert set-to-master-face
 858     Map<label> masterToSet(masterFace_.size());
 859
 860     forAll(masterFace_, setI)
 861     {
 862         if (masterFace_[setI] >= 0)
 863         {
 864             masterToSet.insert(masterFace_[setI], setI);
 865         }
 866     }
 867
 868     forAll(masterFaces, i)
 869     {
 870         label masterFaceI = masterFaces[i];
 871
 872         Map<label>::const_iterator iter = masterToSet.find(masterFaceI);
 873
 874         if (iter == masterToSet.end())
 875         {
 876             FatalErrorIn
 877             (
 878                 "combineFaces::setUnrefinement"
 879                 "(const labelList&, polyTopoChange&"
 880                 ", Map<label>&, Map<label>&, Map<label>&)"
 881             )   << "Master face " << masterFaceI
 882                 << " is not the master of one of the merge sets"
 883                 << " or has already been merged"
 884                 << abort(FatalError);
 885         }
 886
 887         label setI = iter();
 888
 889
 890         // Update faces of the merge setI for reintroduced vertices
 891         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 892
 893         faceList& faces = faceSetsVertices_[setI];
 894
 895         if (faces.size() == 0)
 896         {
 897             FatalErrorIn
 898             (
 899                 "combineFaces::setUnrefinement"
 900                 "(const labelList&, polyTopoChange&"
 901                 ", Map<label>&, Map<label>&, Map<label>&)"
 902             )   << "Set " << setI << " with master face " << masterFaceI
 903                 << " has already been merged." << abort(FatalError);
 904         }
 905
 906         forAll(faces, i)
 907         {
 908             face& f = faces[i];
 909
 910             forAll(f, fp)
 911             {
 912                 label pointI = f[fp];
 913
 914                 if (pointI < 0)
 915                 {
 916                     label localI = -pointI-1;
 917
 918                     if (addedPoints[localI] == -1)
 919                     {
 920                         // First occurrence of saved point. Reintroduce point
 921                         addedPoints[localI] = meshMod.setAction
 922                         (
 923                             polyAddPoint
 924                             (
 925                                 savedPoints_[localI],   // point
 926                                 -1,                     // master point
 927                                 -1,                     // zone for point
 928                                 true                    // supports a cell
 929                             )
 930                         );
 931                         restoredPoints.insert
 932                         (
 933                             addedPoints[localI],        // current point label
 934                             savedPointLabels_[localI]   // point label when it
 935                                                         // was stored
 936                         );
 937                     }
 938                     f[fp] = addedPoints[localI];
 939                 }
 940             }
 941         }
 942
 943
 944         // Restore
 945         // ~~~~~~~
 946
 947         label own = mesh_.faceOwner()[masterFaceI];
 948         label zoneID = mesh_.faceZones().whichZone(masterFaceI);
 949         bool zoneFlip = false;
 950         if (zoneID >= 0)
 951         {
 952             const faceZone& fZone = mesh_.faceZones()[zoneID];
 953             zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
 954         }
 955         label patchI = mesh_.boundaryMesh().whichPatch(masterFaceI);
 956
 957         if (mesh_.boundaryMesh()[patchI].coupled())
 958         {
 959             FatalErrorIn
 960             (
 961                 "combineFaces::setUnrefinement"
 962                 "(const labelList&, polyTopoChange&"
 963                 ", Map<label>&, Map<label>&, Map<label>&)"
 964             )   << "Master face " << masterFaceI << " is on coupled patch "
 965                 << mesh_.boundaryMesh()[patchI].name()
 966                 << abort(FatalError);
 967         }
 968
 969         //Pout<< "Restoring new master face " << masterFaceI
 970         //    << " to vertices " << faces[0] << endl;
 971
 972         // Modify the master face.
 973         meshMod.setAction
 974         (
 975             polyModifyFace
 976             (
 977                 faces[0],                       // original face
 978                 masterFaceI,                    // label of face
 979                 own,                            // owner
 980                 -1,                             // neighbour
 981                 false,                          // face flip
 982                 patchI,                         // patch for face
 983                 false,                          // remove from zone
 984                 zoneID,                         // zone for face
 985                 zoneFlip                        // face flip in zone
 986             )
 987         );
 988
 989         // Add the previously removed faces
 990         for (label i = 1; i < faces.size(); i++)
 991         {
 992             //Pout<< "Restoring removed face with vertices " << faces[i]
 993             //    << endl;
 994
 995             meshMod.setAction
 996             (
 997                 polyAddFace
 998                 (
 999                     faces[i],        // vertices
1000                     own,                    // owner,
1001                     -1,                     // neighbour,
1002                     -1,                     // masterPointID,
1003                     -1,                     // masterEdgeID,
1004                     masterFaceI,             // masterFaceID,
1005                     false,                  // flipFaceFlux,
1006                     patchI,                 // patchID,
1007                     zoneID,                 // zoneID,
1008                     zoneFlip                // zoneFlip
1009                 )
1010             );
1011         }
1012
1013         // Clear out restored set
1014         faceSetsVertices_[setI].clear();
1015         masterFace_[setI] = -1;
1016     }
1017 }
1018
1019
1020 // ************************************************************************* //