src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.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 "BiIndirectList.H"
  28 #include "removePoints.H"
  29 #include "PstreamReduceOps.H"
  30 #include "polyMesh.H"
  31 #include "polyTopoChange.H"
  32 #include "polyRemovePoint.H"
  33 #include "polyAddPoint.H"
  34 #include "polyModifyFace.H"
  35 #include "syncTools.H"
  36 #include "wallPoint.H"  // only to use 'greatPoint'
  37 #include "faceSet.H"
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 defineTypeNameAndDebug(removePoints, 0);
  45
  46 }
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 // Change the vertices of the face whilst keeping everything else the same.
  51 void Foam::removePoints::modifyFace
  52 (
  53     const label faceI,
  54     const face& newFace,
  55     polyTopoChange& meshMod
  56 ) const
  57 {
  58     // Get other face data.
  59     label patchI = -1;
  60     label owner = mesh_.faceOwner()[faceI];
  61     label neighbour = -1;
  62
  63     if (mesh_.isInternalFace(faceI))
  64     {
  65         neighbour = mesh_.faceNeighbour()[faceI];
  66     }
  67     else
  68     {
  69         patchI = mesh_.boundaryMesh().whichPatch(faceI);
  70     }
  71
  72     label zoneID = mesh_.faceZones().whichZone(faceI);
  73
  74     bool zoneFlip = false;
  75
  76     if (zoneID >= 0)
  77     {
  78         const faceZone& fZone = mesh_.faceZones()[zoneID];
  79
  80         zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
  81     }
  82
  83     meshMod.setAction
  84     (
  85         polyModifyFace
  86         (
  87             newFace,        // modified face
  88             faceI,          // label of face being modified
  89             owner,          // owner
  90             neighbour,      // neighbour
  91             false,          // face flip
  92             patchI,         // patch for face
  93             false,          // remove from zone
  94             zoneID,         // zone for face
  95             zoneFlip        // face flip in zone
  96         )
  97     );
  98 }
  99
 100
 101 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 102
 103 // Construct from mesh
 104 Foam::removePoints::removePoints
 105 (
 106     const polyMesh& mesh,
 107     const bool undoable
 108 )
 109 :
 110     mesh_(mesh),
 111     undoable_(undoable)
 112 {}
 113
 114
 115 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 116
 117 Foam::label Foam::removePoints::countPointUsage
 118 (
 119     const scalar minCos,
 120     boolList& pointCanBeDeleted
 121 ) const
 122 {
 123     // Containers to store two edges per point:
 124     // -1   : not filled
 125     // >= 0 : edge label
 126     // -2   : more than two edges for point
 127     labelList edge0(mesh_.nPoints(), -1);
 128     labelList edge1(mesh_.nPoints(), -1);
 129
 130     const edgeList& edges = mesh_.edges();
 131
 132     forAll(edges, edgeI)
 133     {
 134         const edge& e = edges[edgeI];
 135
 136         forAll(e, eI)
 137         {
 138             label pointI = e[eI];
 139
 140             if (edge0[pointI] == -2)
 141             {
 142                 // Already too many edges
 143             }
 144             else if (edge0[pointI] == -1)
 145             {
 146                 // Store first edge using point
 147                 edge0[pointI] = edgeI;
 148             }
 149             else
 150             {
 151                 // Already one edge using point. Check second container.
 152
 153                 if (edge1[pointI] == -1)
 154                 {
 155                     // Store second edge using point
 156                     edge1[pointI] = edgeI;
 157                 }
 158                 else
 159                 {
 160                     // Third edge using point. Mark.
 161                     edge0[pointI] = -2;
 162                     edge1[pointI] = -2;
 163                 }
 164             }
 165         }
 166     }
 167
 168
 169     // Check the ones used by only 2 edges that these are sufficiently in line.
 170     const pointField& points = mesh_.points();
 171
 172     pointCanBeDeleted.setSize(mesh_.nPoints());
 173     pointCanBeDeleted = false;
 174     label nDeleted = 0;
 175
 176     forAll(edge0, pointI)
 177     {
 178         if (edge0[pointI] >= 0 && edge1[pointI] >= 0)
 179         {
 180             // Point used by two edges exactly
 181
 182             const edge& e0 = edges[edge0[pointI]];
 183             const edge& e1 = edges[edge1[pointI]];
 184
 185             label common = e0.commonVertex(e1);
 186             label vLeft = e0.otherVertex(common);
 187             label vRight = e1.otherVertex(common);
 188
 189             vector e0Vec = points[common] - points[vLeft];
 190             e0Vec /= mag(e0Vec) + VSMALL;
 191
 192             vector e1Vec = points[vRight] - points[common];
 193             e1Vec /= mag(e1Vec) + VSMALL;
 194
 195             if ((e0Vec & e1Vec) > minCos)
 196             {
 197                 pointCanBeDeleted[pointI] = true;
 198                 nDeleted++;
 199             }
 200         }
 201         else if (edge0[pointI] == -1)
 202         {
 203             // point not used at all
 204             pointCanBeDeleted[pointI] = true;
 205             nDeleted++;
 206         }
 207     }
 208     edge0.clear();
 209     edge1.clear();
 210
 211
 212     // Protect any points on faces that would collapse down to nothing
 213     // No particular intelligence so might protect too many points
 214     forAll(mesh_.faces(), faceI)
 215     {
 216         const face& f = mesh_.faces()[faceI];
 217
 218         label nCollapse = 0;
 219         forAll(f, fp)
 220         {
 221             if (pointCanBeDeleted[f[fp]])
 222             {
 223                 nCollapse++;
 224             }
 225         }
 226
 227         if ((f.size() - nCollapse) < 3)
 228         {
 229             // Just unmark enough points
 230             forAll(f, fp)
 231             {
 232                 if (pointCanBeDeleted[f[fp]])
 233                 {
 234                     pointCanBeDeleted[f[fp]] = false;
 235                     --nCollapse;
 236                     if (nCollapse == 0)
 237                     {
 238                         break;
 239                     }
 240                 }
 241             }
 242         }
 243     }
 244
 245
 246     // Point can be deleted only if all processors want to delete it
 247     syncTools::syncPointList
 248     (
 249         mesh_,
 250         pointCanBeDeleted,
 251         andEqOp<bool>(),
 252         true,               // null value
 253         false               // no separation
 254     );
 255
 256     return returnReduce(nDeleted, sumOp<label>());
 257 }
 258
 259
 260 void Foam::removePoints::setRefinement
 261 (
 262     const boolList& pointCanBeDeleted,
 263     polyTopoChange& meshMod
 264 )
 265 {
 266     // Count deleted points
 267     label nDeleted = 0;
 268     forAll(pointCanBeDeleted, pointI)
 269     {
 270         if (pointCanBeDeleted[pointI])
 271         {
 272             nDeleted++;
 273         }
 274     }
 275
 276     // Faces (in mesh face labels) affected by points removed. Will hopefully
 277     // be only a few.
 278     labelHashSet facesAffected(4*nDeleted);
 279
 280
 281     // Undo: from global mesh point to index in savedPoint_
 282     Map<label> pointToSaved;
 283
 284     // Size undo storage
 285     if (undoable_)
 286     {
 287         savedPoints_.setSize(nDeleted);
 288         pointToSaved.resize(2*nDeleted);
 289     }
 290
 291
 292     // Remove points
 293     // ~~~~~~~~~~~~~
 294
 295     nDeleted = 0;
 296
 297     forAll(pointCanBeDeleted, pointI)
 298     {
 299         if (pointCanBeDeleted[pointI])
 300         {
 301             if (undoable_)
 302             {
 303                 pointToSaved.insert(pointI, nDeleted);
 304                 savedPoints_[nDeleted++] = mesh_.points()[pointI];
 305             }
 306             meshMod.setAction(polyRemovePoint(pointI));
 307
 308             // Store faces affected
 309             const labelList& pFaces = mesh_.pointFaces()[pointI];
 310
 311             forAll(pFaces, i)
 312             {
 313                 facesAffected.insert(pFaces[i]);
 314             }
 315         }
 316     }
 317
 318
 319
 320     // Update faces
 321     // ~~~~~~~~~~~~
 322
 323
 324     if (undoable_)
 325     {
 326         savedFaceLabels_.setSize(facesAffected.size());
 327         savedFaces_.setSize(facesAffected.size());
 328     }
 329     label nSaved = 0;
 330
 331     forAllConstIter(labelHashSet, facesAffected, iter)
 332     {
 333         label faceI = iter.key();
 334
 335         const face& f = mesh_.faces()[faceI];
 336
 337         face newFace(f.size());
 338
 339         label newI = 0;
 340
 341         forAll(f, fp)
 342         {
 343             label pointI = f[fp];
 344
 345             if (!pointCanBeDeleted[pointI])
 346             {
 347                 newFace[newI++] = pointI;
 348             }
 349         }
 350         newFace.setSize(newI);
 351
 352         // Actually change the face to the new vertices
 353         modifyFace(faceI, newFace, meshMod);
 354
 355         // Save the face. Negative indices are into savedPoints_
 356         if (undoable_)
 357         {
 358             savedFaceLabels_[nSaved] = faceI;
 359
 360             face& savedFace = savedFaces_[nSaved++];
 361             savedFace.setSize(f.size());
 362
 363             forAll(f, fp)
 364             {
 365                 label pointI = f[fp];
 366
 367                 if (pointCanBeDeleted[pointI])
 368                 {
 369                     savedFace[fp] = -pointToSaved[pointI]-1;
 370                 }
 371                 else
 372                 {
 373                     savedFace[fp] = pointI;
 374                 }
 375             }
 376         }
 377     }
 378
 379     if (undoable_)
 380     {
 381         // DEBUG: Compare the stored faces with the current ones.
 382         if (debug)
 383         {
 384             forAll(savedFaceLabels_, saveI)
 385             {
 386                 // Points from the mesh
 387                 List<point> meshPoints
 388                 (
 389                     UIndirectList<point>
 390                     (
 391                         mesh_.points(),
 392                         mesh_.faces()[savedFaceLabels_[saveI]]  // mesh face
 393                     )
 394                 );
 395
 396                 // Points from the saved face
 397                 List<point> keptPoints
 398                 (
 399                     BiIndirectList<point>
 400                     (
 401                         mesh_.points(),
 402                         savedPoints_,
 403                         savedFaces_[saveI]  // saved face
 404                     )
 405                 );
 406
 407                 if (meshPoints != keptPoints)
 408                 {
 409                     FatalErrorIn("setRefinement")
 410                         << "faceI:" << savedFaceLabels_[saveI] << nl
 411                         << "meshPoints:" << meshPoints << nl
 412                         << "keptPoints:" << keptPoints << nl
 413                         << abort(FatalError);
 414                 }
 415             }
 416         }
 417     }
 418 }
 419
 420
 421 void Foam::removePoints::updateMesh(const mapPolyMesh& map)
 422 {
 423     if (undoable_)
 424     {
 425         forAll(savedFaceLabels_, localI)
 426         {
 427             if (savedFaceLabels_[localI] >= 0)
 428             {
 429                 label newFaceI = map.reverseFaceMap()[savedFaceLabels_[localI]];
 430
 431                 if (newFaceI == -1)
 432                 {
 433                     FatalErrorIn
 434                     (
 435                         "removePoints::updateMesh(const mapPolyMesh&)"
 436                     )   << "Old face " << savedFaceLabels_[localI]
 437                         << " seems to have dissapeared."
 438                         << abort(FatalError);
 439                 }
 440                 savedFaceLabels_[localI] = newFaceI;
 441             }
 442         }
 443
 444
 445         // Renumber mesh vertices (indices >=0). Leave saved vertices
 446         // (<0) intact.
 447         forAll(savedFaces_, i)
 448         {
 449             face& f = savedFaces_[i];
 450
 451             forAll(f, fp)
 452             {
 453                 label pointI = f[fp];
 454
 455                 if (pointI >= 0)
 456                 {
 457                     f[fp] = map.reversePointMap()[pointI];
 458
 459                     if (f[fp] == -1)
 460                     {
 461                         FatalErrorIn
 462                         (
 463                             "removePoints::updateMesh(const mapPolyMesh&)"
 464                         )   << "Old point " << pointI
 465                             << " seems to have dissapeared."
 466                             << abort(FatalError);
 467                     }
 468                 }
 469             }
 470         }
 471
 472
 473         // DEBUG: Compare the stored faces with the current ones.
 474         if (debug)
 475         {
 476             forAll(savedFaceLabels_, saveI)
 477             {
 478                 if (savedFaceLabels_[saveI] >= 0)
 479                 {
 480                     const face& f = mesh_.faces()[savedFaceLabels_[saveI]];
 481
 482                     // Get kept points of saved faces.
 483                     const face& savedFace = savedFaces_[saveI];
 484
 485                     face keptFace(savedFace.size());
 486                     label keptFp = 0;
 487
 488                     forAll(savedFace, fp)
 489                     {
 490                         label pointI = savedFace[fp];
 491
 492                         if (pointI >= 0)
 493                         {
 494                             keptFace[keptFp++] = pointI;
 495                         }
 496                     }
 497                     keptFace.setSize(keptFp);
 498
 499                     // Compare as faces (since f might have rotated and
 500                     // face::operator== takes care of this)
 501                     if (keptFace != f)
 502                     {
 503                         FatalErrorIn("setRefinement")
 504                             << "faceI:" << savedFaceLabels_[saveI] << nl
 505                             << "face:" << f << nl
 506                             << "keptFace:" << keptFace << nl
 507                             << "saved points:"
 508                             <<  BiIndirectList<point>
 509                                 (
 510                                     mesh_.points(),
 511                                     savedPoints_,
 512                                     savedFace
 513                                 )() << nl
 514                             << abort(FatalError);
 515                     }
 516                 }
 517             }
 518         }
 519     }
 520 }
 521
 522
 523 // Given list of faces to undo picks up the local indices of the faces
 524 // to restore. Additionally it also picks up all the faces that use
 525 // any of the deleted points.
 526 void Foam::removePoints::getUnrefimentSet
 527 (
 528     const labelList& undoFaces,
 529     labelList& localFaces,
 530     labelList& localPoints
 531 ) const
 532 {
 533     if (!undoable_)
 534     {
 535         FatalErrorIn
 536         (
 537             "removePoints::getUnrefimentSet(const labelList&"
 538             ", labelList&, labelList&) const"
 539         )   << "removePoints not constructed with"
 540             << " unrefinement capability."
 541             << abort(FatalError);
 542     }
 543
 544     if (debug)
 545     {
 546         // Check if synced.
 547         faceSet undoFacesSet(mesh_, "undoFacesSet", undoFaces);
 548         label sz = undoFacesSet.size();
 549
 550         undoFacesSet.sync(mesh_);
 551         if (sz != undoFacesSet.size())
 552         {
 553             FatalErrorIn
 554             (
 555                 "removePoints::getUnrefimentSet(const labelList&"
 556                 ", labelList&, labelList&) const"
 557             )   << "undoFaces not synchronised across coupled faces." << endl
 558                 << "Size before sync:" << sz
 559                 << "  after sync:" << undoFacesSet.size()
 560                 << abort(FatalError);
 561         }
 562     }
 563
 564
 565     // Problem: input undoFaces are synced. But e.g.
 566     // two faces, A (uncoupled) and B(coupled) share a deleted point. A gets
 567     // passed in to be restored. Now picking up the deleted point and the
 568     // original faces using it picks up face B. But not its coupled neighbour!
 569     // Problem is that we cannot easily synchronise the deleted points
 570     // (e.g. syncPointList) since they're not in the mesh anymore - only the
 571     // faces are. So we have to collect the points-to-restore as indices
 572     // in the faces (which is information we can synchronise)
 573
 574
 575
 576     // Mark points-to-restore
 577     labelHashSet localPointsSet(undoFaces.size());
 578
 579     {
 580         // Create set of faces to be restored
 581         labelHashSet undoFacesSet(undoFaces);
 582
 583         forAll(savedFaceLabels_, saveI)
 584         {
 585             if (savedFaceLabels_[saveI] < 0)
 586             {
 587                 FatalErrorIn
 588                 (
 589                     "removePoints::getUnrefimentSet(const labelList&"
 590                     ", labelList&, labelList&) const"
 591                 )   << "Illegal face label " << savedFaceLabels_[saveI]
 592                     << " at index " << saveI
 593                     << abort(FatalError);
 594             }
 595
 596             if (undoFacesSet.found(savedFaceLabels_[saveI]))
 597             {
 598                 const face& savedFace = savedFaces_[saveI];
 599
 600                 forAll(savedFace, fp)
 601                 {
 602                     if (savedFace[fp] < 0)
 603                     {
 604                         label savedPointI = -savedFace[fp]-1;
 605
 606                         if (savedPoints_[savedPointI] == wallPoint::greatPoint)
 607                         {
 608                             FatalErrorIn
 609                             (
 610                                 "removePoints::getUnrefimentSet"
 611                                 "(const labelList&, labelList&, labelList&)"
 612                                 " const"
 613                             )   << "Trying to restore point " << savedPointI
 614                                 << " from mesh face " << savedFaceLabels_[saveI]
 615                                 << " saved face:" << savedFace
 616                                 << " which has already been undone."
 617                                 << abort(FatalError);
 618                         }
 619
 620                         localPointsSet.insert(savedPointI);
 621                     }
 622                 }
 623             }
 624         }
 625
 626
 627         // Per boundary face, per index in face whether the point needs
 628         // restoring. Note that this is over all saved faces, not just over
 629         // the ones in undoFaces.
 630
 631         boolListList faceVertexRestore(mesh_.nFaces()-mesh_.nInternalFaces());
 632
 633         // Populate with my local points-to-restore.
 634         forAll(savedFaces_, saveI)
 635         {
 636             label bFaceI = savedFaceLabels_[saveI] - mesh_.nInternalFaces();
 637
 638             if (bFaceI >= 0)
 639             {
 640                 const face& savedFace = savedFaces_[saveI];
 641
 642                 boolList& fRestore = faceVertexRestore[bFaceI];
 643
 644                 fRestore.setSize(savedFace.size());
 645                 fRestore = false;
 646
 647                 forAll(savedFace, fp)
 648                 {
 649                     if (savedFace[fp] < 0)
 650                     {
 651                         label savedPointI = -savedFace[fp]-1;
 652
 653                         if (localPointsSet.found(savedPointI))
 654                         {
 655                             fRestore[fp] = true;
 656                         }
 657                     }
 658                 }
 659             }
 660         }
 661
 662         syncTools::syncBoundaryFaceList
 663         (
 664             mesh_,
 665             faceVertexRestore,
 666             faceEqOp<bool, orEqOp>(),   // special operator to handle faces
 667             false                       // no separation
 668         );
 669
 670         // So now if any of the points-to-restore is used by any coupled face
 671         // anywhere the corresponding index in faceVertexRestore will be set.
 672
 673         // Now combine the localPointSet and the (sychronised)
 674         // boundary-points-to-restore.
 675
 676         forAll(savedFaces_, saveI)
 677         {
 678             label bFaceI = savedFaceLabels_[saveI] - mesh_.nInternalFaces();
 679
 680             if (bFaceI >= 0)
 681             {
 682                 const boolList& fRestore = faceVertexRestore[bFaceI];
 683
 684                 const face& savedFace = savedFaces_[saveI];
 685
 686                 forAll(fRestore, fp)
 687                 {
 688                     // Does neighbour require point restored?
 689                     if (fRestore[fp])
 690                     {
 691                         if (savedFace[fp] >= 0)
 692                         {
 693                             FatalErrorIn
 694                             (
 695                                 "removePoints::getUnrefimentSet"
 696                                 "(const labelList&, labelList&, labelList&)"
 697                                 " const"
 698                             )   << "Problem: on coupled face:"
 699                                 << savedFaceLabels_[saveI]
 700                                 << " fc:"
 701                                 << mesh_.faceCentres()[savedFaceLabels_[saveI]]
 702                                 << endl
 703                                 << " my neighbour tries to restore the vertex"
 704                                 << " at index " << fp
 705                                 << " whereas my saved face:" << savedFace
 706                                 << " does not indicate a deleted vertex"
 707                                 << " at that position."
 708                                 << abort(FatalError);
 709                         }
 710
 711                         label savedPointI = -savedFace[fp]-1;
 712
 713                         localPointsSet.insert(savedPointI);
 714                     }
 715                 }
 716             }
 717         }
 718     }
 719
 720     localPoints = localPointsSet.toc();
 721
 722
 723     // Collect all saved faces using any localPointsSet
 724     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 725
 726     labelHashSet localFacesSet(2*undoFaces.size());
 727
 728     forAll(savedFaces_, saveI)
 729     {
 730         const face& savedFace = savedFaces_[saveI];
 731
 732         forAll(savedFace, fp)
 733         {
 734             if (savedFace[fp] < 0)
 735             {
 736                 label savedPointI = -savedFace[fp]-1;
 737
 738                 if (localPointsSet.found(savedPointI))
 739                 {
 740                     localFacesSet.insert(saveI);
 741                 }
 742             }
 743         }
 744     }
 745     localFaces = localFacesSet.toc();
 746
 747
 748     // Note that at this point the localFaces to restore can contain points
 749     // that are not going to be restored! The localFaces though will
 750     // be guaranteed to be all the ones affected by the restoration of the
 751     // localPoints.
 752 }
 753
 754
 755 void Foam::removePoints::setUnrefinement
 756 (
 757     const labelList& localFaces,
 758     const labelList& localPoints,
 759     polyTopoChange& meshMod
 760 )
 761 {
 762     if (!undoable_)
 763     {
 764         FatalErrorIn
 765         (
 766             "removePoints::setUnrefinement(const labelList&"
 767             ", labelList&, polyTopoChange&)"
 768         )   << "removePoints not constructed with"
 769             << " unrefinement capability."
 770             << abort(FatalError);
 771     }
 772
 773
 774     // Per savedPoint -1 or the restored point label
 775     labelList addedPoints(savedPoints_.size(), -1);
 776
 777     forAll(localPoints, i)
 778     {
 779         label localI = localPoints[i];
 780
 781         if (savedPoints_[localI] == wallPoint::greatPoint)
 782         {
 783             FatalErrorIn
 784             (
 785                 "removePoints::setUnrefinement(const labelList&"
 786                 ", labelList&, polyTopoChange&)"
 787             )   << "Saved point " << localI << " already restored!"
 788                 << abort(FatalError);
 789         }
 790
 791         addedPoints[localI] = meshMod.setAction
 792         (
 793             polyAddPoint
 794             (
 795                 savedPoints_[localI],   // point
 796                 -1,                     // master point
 797                 -1,                     // zone for point
 798                 true                    // supports a cell
 799             )
 800         );
 801
 802         // Mark the restored points so they are not restored again.
 803         savedPoints_[localI] = wallPoint::greatPoint;
 804     }
 805
 806     forAll(localFaces, i)
 807     {
 808         label saveI = localFaces[i];
 809
 810         // Modify indices into saved points (so < 0) to point to the
 811         // added points.
 812         face& savedFace = savedFaces_[saveI];
 813
 814         face newFace(savedFace.size());
 815         label newFp = 0;
 816
 817         bool hasSavedPoints = false;
 818
 819         forAll(savedFace, fp)
 820         {
 821             if (savedFace[fp] < 0)
 822             {
 823                 label addedPointI = addedPoints[-savedFace[fp]-1];
 824
 825                 if (addedPointI != -1)
 826                 {
 827                     savedFace[fp] = addedPointI;
 828                     newFace[newFp++] = addedPointI;
 829                 }
 830                 else
 831                 {
 832                     hasSavedPoints = true;
 833                 }
 834             }
 835             else
 836             {
 837                 newFace[newFp++] = savedFace[fp];
 838             }
 839         }
 840         newFace.setSize(newFp);
 841
 842         modifyFace(savedFaceLabels_[saveI], newFace, meshMod);
 843
 844         if (!hasSavedPoints)
 845         {
 846             // Face fully restored. Mark for compaction later on
 847             savedFaceLabels_[saveI] = -1;
 848             savedFaces_[saveI].clear();
 849         }
 850     }
 851
 852
 853     // Compact face labels
 854     label newSaveI = 0;
 855
 856     forAll(savedFaceLabels_, saveI)
 857     {
 858         if (savedFaceLabels_[saveI] != -1)
 859         {
 860             if (newSaveI != saveI)
 861             {
 862                 savedFaceLabels_[newSaveI] = savedFaceLabels_[saveI];
 863                 savedFaces_[newSaveI].transfer(savedFaces_[saveI]);
 864             }
 865             newSaveI++;
 866         }
 867     }
 868
 869     savedFaceLabels_.setSize(newSaveI);
 870     savedFaces_.setSize(newSaveI);
 871
 872
 873     // Check that all faces have been restored that use any restored points
 874     if (debug)
 875     {
 876         forAll(savedFaceLabels_, saveI)
 877         {
 878             const face& savedFace = savedFaces_[saveI];
 879
 880             forAll(savedFace, fp)
 881             {
 882                 if (savedFace[fp] < 0)
 883                 {
 884                     label addedPointI = addedPoints[-savedFace[fp]-1];
 885
 886                     if (addedPointI != -1)
 887                     {
 888                         FatalErrorIn("setUnrefinement")
 889                             << "Face:" << savedFaceLabels_[saveI]
 890                             << " savedVerts:" << savedFace
 891                             << " uses restored point:" << -savedFace[fp]-1
 892                             << " with new pointlabel:" << addedPointI
 893                             << abort(FatalError);
 894                     }
 895                 }
 896             }
 897         }
 898     }
 899 }
 900
 901
 902 // ************************************************************************* //