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