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     }
 209     edge0.clear();
 210     edge1.clear();
 211
 212
 213     // Point can be deleted only if all processors want to delete it
 214     syncTools::syncPointList
 215     (
 216         mesh_,
 217         pointCanBeDeleted,
 218         andEqOp<bool>(),
 219         true,               // null value
 220         false               // no separation
 221     );
 222
 223     return returnReduce(nDeleted, sumOp<label>());
 224 }
 225
 226
 227 void Foam::removePoints::setRefinement
 228 (
 229     const boolList& pointCanBeDeleted,
 230     polyTopoChange& meshMod
 231 )
 232 {
 233     // Count deleted points
 234     label nDeleted = 0;
 235     forAll(pointCanBeDeleted, pointI)
 236     {
 237         if (pointCanBeDeleted[pointI])
 238         {
 239             nDeleted++;
 240         }
 241     }
 242
 243     // Faces (in mesh face labels) affected by points removed. Will hopefully
 244     // be only a few.
 245     labelHashSet facesAffected(4*nDeleted);
 246
 247
 248     // Undo: from global mesh point to index in savedPoint_
 249     Map<label> pointToSaved;
 250
 251     // Size undo storage
 252     if (undoable_)
 253     {
 254         savedPoints_.setSize(nDeleted);
 255         pointToSaved.resize(2*nDeleted);
 256     }
 257
 258
 259     // Remove points
 260     // ~~~~~~~~~~~~~
 261
 262     nDeleted = 0;
 263
 264     forAll(pointCanBeDeleted, pointI)
 265     {
 266         if (pointCanBeDeleted[pointI])
 267         {
 268             if (undoable_)
 269             {
 270                 pointToSaved.insert(pointI, nDeleted);
 271                 savedPoints_[nDeleted++] = mesh_.points()[pointI];
 272             }
 273             meshMod.setAction(polyRemovePoint(pointI));
 274
 275             // Store faces affected
 276             const labelList& pFaces = mesh_.pointFaces()[pointI];
 277
 278             forAll(pFaces, i)
 279             {
 280                 facesAffected.insert(pFaces[i]);
 281             }
 282         }
 283     }
 284
 285
 286
 287     // Update faces
 288     // ~~~~~~~~~~~~
 289
 290
 291     if (undoable_)
 292     {
 293         savedFaceLabels_.setSize(facesAffected.size());
 294         savedFaces_.setSize(facesAffected.size());
 295     }
 296     label nSaved = 0;
 297
 298     forAllConstIter(labelHashSet, facesAffected, iter)
 299     {
 300         label faceI = iter.key();
 301
 302         const face& f = mesh_.faces()[faceI];
 303
 304         face newFace(f.size());
 305
 306         label newI = 0;
 307
 308         forAll(f, fp)
 309         {
 310             label pointI = f[fp];
 311
 312             if (!pointCanBeDeleted[pointI])
 313             {
 314                 newFace[newI++] = pointI;
 315             }
 316         }
 317         newFace.setSize(newI);
 318
 319         // Actually change the face to the new vertices
 320         modifyFace(faceI, newFace, meshMod);
 321
 322         // Save the face. Negative indices are into savedPoints_
 323         if (undoable_)
 324         {
 325             savedFaceLabels_[nSaved] = faceI;
 326
 327             face& savedFace = savedFaces_[nSaved++];
 328             savedFace.setSize(f.size());
 329
 330             forAll(f, fp)
 331             {
 332                 label pointI = f[fp];
 333
 334                 if (pointCanBeDeleted[pointI])
 335                 {
 336                     savedFace[fp] = -pointToSaved[pointI]-1;
 337                 }
 338                 else
 339                 {
 340                     savedFace[fp] = pointI;
 341                 }
 342             }
 343         }
 344     }
 345
 346     if (undoable_)
 347     {
 348         // DEBUG: Compare the stored faces with the current ones.
 349         if (debug)
 350         {
 351             forAll(savedFaceLabels_, saveI)
 352             {
 353                 // Points from the mesh
 354                 List<point> meshPoints
 355                 (
 356                     UIndirectList<point>
 357                     (
 358                         mesh_.points(),
 359                         mesh_.faces()[savedFaceLabels_[saveI]]  // mesh face
 360                     )
 361                 );
 362
 363                 // Points from the saved face
 364                 List<point> keptPoints
 365                 (
 366                     BiIndirectList<point>
 367                     (
 368                         mesh_.points(),
 369                         savedPoints_,
 370                         savedFaces_[saveI]  // saved face
 371                     )
 372                 );
 373
 374                 if (meshPoints != keptPoints)
 375                 {
 376                     FatalErrorIn("setRefinement")
 377                         << "faceI:" << savedFaceLabels_[saveI] << nl
 378                         << "meshPoints:" << meshPoints << nl
 379                         << "keptPoints:" << keptPoints << nl
 380                         << abort(FatalError);
 381                 }
 382             }
 383         }
 384     }
 385 }
 386
 387
 388 void Foam::removePoints::updateMesh(const mapPolyMesh& map)
 389 {
 390     if (undoable_)
 391     {
 392         forAll(savedFaceLabels_, localI)
 393         {
 394             if (savedFaceLabels_[localI] >= 0)
 395             {
 396                 label newFaceI = map.reverseFaceMap()[savedFaceLabels_[localI]];
 397
 398                 if (newFaceI == -1)
 399                 {
 400                     FatalErrorIn
 401                     (
 402                         "removePoints::updateMesh(const mapPolyMesh&)"
 403                     )   << "Old face " << savedFaceLabels_[localI]
 404                         << " seems to have dissapeared."
 405                         << abort(FatalError);
 406                 }
 407                 savedFaceLabels_[localI] = newFaceI;
 408             }
 409         }
 410
 411
 412         // Renumber mesh vertices (indices >=0). Leave saved vertices
 413         // (<0) intact.
 414         forAll(savedFaces_, i)
 415         {
 416             face& f = savedFaces_[i];
 417
 418             forAll(f, fp)
 419             {
 420                 label pointI = f[fp];
 421
 422                 if (pointI >= 0)
 423                 {
 424                     f[fp] = map.reversePointMap()[pointI];
 425
 426                     if (f[fp] == -1)
 427                     {
 428                         FatalErrorIn
 429                         (
 430                             "removePoints::updateMesh(const mapPolyMesh&)"
 431                         )   << "Old point " << pointI
 432                             << " seems to have dissapeared."
 433                             << abort(FatalError);
 434                     }
 435                 }
 436             }
 437         }
 438
 439
 440         // DEBUG: Compare the stored faces with the current ones.
 441         if (debug)
 442         {
 443             forAll(savedFaceLabels_, saveI)
 444             {
 445                 if (savedFaceLabels_[saveI] >= 0)
 446                 {
 447                     const face& f = mesh_.faces()[savedFaceLabels_[saveI]];
 448
 449                     // Get kept points of saved faces.
 450                     const face& savedFace = savedFaces_[saveI];
 451
 452                     face keptFace(savedFace.size());
 453                     label keptFp = 0;
 454
 455                     forAll(savedFace, fp)
 456                     {
 457                         label pointI = savedFace[fp];
 458
 459                         if (pointI >= 0)
 460                         {
 461                             keptFace[keptFp++] = pointI;
 462                         }
 463                     }
 464                     keptFace.setSize(keptFp);
 465
 466                     // Compare as faces (since f might have rotated and
 467                     // face::operator== takes care of this)
 468                     if (keptFace != f)
 469                     {
 470                         FatalErrorIn("setRefinement")
 471                             << "faceI:" << savedFaceLabels_[saveI] << nl
 472                             << "face:" << f << nl
 473                             << "keptFace:" << keptFace << nl
 474                             << "saved points:"
 475                             <<  BiIndirectList<point>
 476                                 (
 477                                     mesh_.points(),
 478                                     savedPoints_,
 479                                     savedFace
 480                                 )() << nl
 481                             << abort(FatalError);
 482                     }
 483                 }
 484             }
 485         }
 486     }
 487 }
 488
 489
 490 // Given list of faces to undo picks up the local indices of the faces
 491 // to restore. Additionally it also picks up all the faces that use
 492 // any of the deleted points.
 493 void Foam::removePoints::getUnrefimentSet
 494 (
 495     const labelList& undoFaces,
 496     labelList& localFaces,
 497     labelList& localPoints
 498 ) const
 499 {
 500     if (!undoable_)
 501     {
 502         FatalErrorIn
 503         (
 504             "removePoints::getUnrefimentSet(const labelList&"
 505             ", labelList&, labelList&) const"
 506         )   << "removePoints not constructed with"
 507             << " unrefinement capability."
 508             << abort(FatalError);
 509     }
 510
 511     if (debug)
 512     {
 513         // Check if synced.
 514         faceSet undoFacesSet(mesh_, "undoFacesSet", undoFaces);
 515         label sz = undoFacesSet.size();
 516
 517         undoFacesSet.sync(mesh_);
 518         if (sz != undoFacesSet.size())
 519         {
 520             FatalErrorIn
 521             (
 522                 "removePoints::getUnrefimentSet(const labelList&"
 523                 ", labelList&, labelList&) const"
 524             )   << "undoFaces not synchronised across coupled faces." << endl
 525                 << "Size before sync:" << sz
 526                 << "  after sync:" << undoFacesSet.size()
 527                 << abort(FatalError);
 528         }
 529     }
 530
 531
 532     // Problem: input undoFaces are synced. But e.g.
 533     // two faces, A (uncoupled) and B(coupled) share a deleted point. A gets
 534     // passed in to be restored. Now picking up the deleted point and the
 535     // original faces using it picks up face B. But not its coupled neighbour!
 536     // Problem is that we cannot easily synchronise the deleted points
 537     // (e.g. syncPointList) since they're not in the mesh anymore - only the
 538     // faces are. So we have to collect the points-to-restore as indices
 539     // in the faces (which is information we can synchronise)
 540
 541
 542
 543     // Mark points-to-restore
 544     labelHashSet localPointsSet(undoFaces.size());
 545
 546     {
 547         // Create set of faces to be restored
 548         labelHashSet undoFacesSet(undoFaces);
 549
 550         forAll(savedFaceLabels_, saveI)
 551         {
 552             if (savedFaceLabels_[saveI] < 0)
 553             {
 554                 FatalErrorIn
 555                 (
 556                     "removePoints::getUnrefimentSet(const labelList&"
 557                     ", labelList&, labelList&) const"
 558                 )   << "Illegal face label " << savedFaceLabels_[saveI]
 559                     << " at index " << saveI
 560                     << abort(FatalError);
 561             }
 562
 563             if (undoFacesSet.found(savedFaceLabels_[saveI]))
 564             {
 565                 const face& savedFace = savedFaces_[saveI];
 566
 567                 forAll(savedFace, fp)
 568                 {
 569                     if (savedFace[fp] < 0)
 570                     {
 571                         label savedPointI = -savedFace[fp]-1;
 572
 573                         if (savedPoints_[savedPointI] == wallPoint::greatPoint)
 574                         {
 575                             FatalErrorIn
 576                             (
 577                                 "removePoints::getUnrefimentSet"
 578                                 "(const labelList&, labelList&, labelList&)"
 579                                 " const"
 580                             )   << "Trying to restore point " << savedPointI
 581                                 << " from mesh face " << savedFaceLabels_[saveI]
 582                                 << " saved face:" << savedFace
 583                                 << " which has already been undone."
 584                                 << abort(FatalError);
 585                         }
 586
 587                         localPointsSet.insert(savedPointI);
 588                     }
 589                 }
 590             }
 591         }
 592
 593
 594         // Per boundary face, per index in face whether the point needs
 595         // restoring. Note that this is over all saved faces, not just over
 596         // the ones in undoFaces.
 597
 598         boolListList faceVertexRestore(mesh_.nFaces()-mesh_.nInternalFaces());
 599
 600         // Populate with my local points-to-restore.
 601         forAll(savedFaces_, saveI)
 602         {
 603             label bFaceI = savedFaceLabels_[saveI] - mesh_.nInternalFaces();
 604
 605             if (bFaceI >= 0)
 606             {
 607                 const face& savedFace = savedFaces_[saveI];
 608
 609                 boolList& fRestore = faceVertexRestore[bFaceI];
 610
 611                 fRestore.setSize(savedFace.size());
 612                 fRestore = false;
 613
 614                 forAll(savedFace, fp)
 615                 {
 616                     if (savedFace[fp] < 0)
 617                     {
 618                         label savedPointI = -savedFace[fp]-1;
 619
 620                         if (localPointsSet.found(savedPointI))
 621                         {
 622                             fRestore[fp] = true;
 623                         }
 624                     }
 625                 }
 626             }
 627         }
 628
 629         syncTools::syncBoundaryFaceList
 630         (
 631             mesh_,
 632             faceVertexRestore,
 633             faceEqOp<bool, orEqOp>(),   // special operator to handle faces
 634             false                       // no separation
 635         );
 636
 637         // So now if any of the points-to-restore is used by any coupled face
 638         // anywhere the corresponding index in faceVertexRestore will be set.
 639
 640         // Now combine the localPointSet and the (sychronised)
 641         // boundary-points-to-restore.
 642
 643         forAll(savedFaces_, saveI)
 644         {
 645             label bFaceI = savedFaceLabels_[saveI] - mesh_.nInternalFaces();
 646
 647             if (bFaceI >= 0)
 648             {
 649                 const boolList& fRestore = faceVertexRestore[bFaceI];
 650
 651                 const face& savedFace = savedFaces_[saveI];
 652
 653                 forAll(fRestore, fp)
 654                 {
 655                     // Does neighbour require point restored?
 656                     if (fRestore[fp])
 657                     {
 658                         if (savedFace[fp] >= 0)
 659                         {
 660                             FatalErrorIn
 661                             (
 662                                 "removePoints::getUnrefimentSet"
 663                                 "(const labelList&, labelList&, labelList&)"
 664                                 " const"
 665                             )   << "Problem: on coupled face:"
 666                                 << savedFaceLabels_[saveI]
 667                                 << " fc:"
 668                                 << mesh_.faceCentres()[savedFaceLabels_[saveI]]
 669                                 << endl
 670                                 << " my neighbour tries to restore the vertex"
 671                                 << " at index " << fp
 672                                 << " whereas my saved face:" << savedFace
 673                                 << " does not indicate a deleted vertex"
 674                                 << " at that position."
 675                                 << abort(FatalError);
 676                         }
 677
 678                         label savedPointI = -savedFace[fp]-1;
 679
 680                         localPointsSet.insert(savedPointI);
 681                     }
 682                 }
 683             }
 684         }
 685     }
 686
 687     localPoints = localPointsSet.toc();
 688
 689
 690     // Collect all saved faces using any localPointsSet
 691     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 692
 693     labelHashSet localFacesSet(2*undoFaces.size());
 694
 695     forAll(savedFaces_, saveI)
 696     {
 697         const face& savedFace = savedFaces_[saveI];
 698
 699         forAll(savedFace, fp)
 700         {
 701             if (savedFace[fp] < 0)
 702             {
 703                 label savedPointI = -savedFace[fp]-1;
 704
 705                 if (localPointsSet.found(savedPointI))
 706                 {
 707                     localFacesSet.insert(saveI);
 708                 }
 709             }
 710         }
 711     }
 712     localFaces = localFacesSet.toc();
 713
 714
 715     // Note that at this point the localFaces to restore can contain points
 716     // that are not going to be restored! The localFaces though will
 717     // be guaranteed to be all the ones affected by the restoration of the
 718     // localPoints.
 719 }
 720
 721
 722 void Foam::removePoints::setUnrefinement
 723 (
 724     const labelList& localFaces,
 725     const labelList& localPoints,
 726     polyTopoChange& meshMod
 727 )
 728 {
 729     if (!undoable_)
 730     {
 731         FatalErrorIn
 732         (
 733             "removePoints::setUnrefinement(const labelList&"
 734             ", labelList&, polyTopoChange&)"
 735         )   << "removePoints not constructed with"
 736             << " unrefinement capability."
 737             << abort(FatalError);
 738     }
 739
 740
 741     // Per savedPoint -1 or the restored point label
 742     labelList addedPoints(savedPoints_.size(), -1);
 743
 744     forAll(localPoints, i)
 745     {
 746         label localI = localPoints[i];
 747
 748         if (savedPoints_[localI] == wallPoint::greatPoint)
 749         {
 750             FatalErrorIn
 751             (
 752                 "removePoints::setUnrefinement(const labelList&"
 753                 ", labelList&, polyTopoChange&)"
 754             )   << "Saved point " << localI << " already restored!"
 755                 << abort(FatalError);
 756         }
 757
 758         addedPoints[localI] = meshMod.setAction
 759         (
 760             polyAddPoint
 761             (
 762                 savedPoints_[localI],   // point
 763                 -1,                     // master point
 764                 -1,                     // zone for point
 765                 true                    // supports a cell
 766             )
 767         );
 768
 769         // Mark the restored points so they are not restored again.
 770         savedPoints_[localI] = wallPoint::greatPoint;
 771     }
 772
 773     forAll(localFaces, i)
 774     {
 775         label saveI = localFaces[i];
 776
 777         // Modify indices into saved points (so < 0) to point to the
 778         // added points.
 779         face& savedFace = savedFaces_[saveI];
 780
 781         face newFace(savedFace.size());
 782         label newFp = 0;
 783
 784         bool hasSavedPoints = false;
 785
 786         forAll(savedFace, fp)
 787         {
 788             if (savedFace[fp] < 0)
 789             {
 790                 label addedPointI = addedPoints[-savedFace[fp]-1];
 791
 792                 if (addedPointI != -1)
 793                 {
 794                     savedFace[fp] = addedPointI;
 795                     newFace[newFp++] = addedPointI;
 796                 }
 797                 else
 798                 {
 799                     hasSavedPoints = true;
 800                 }
 801             }
 802             else
 803             {
 804                 newFace[newFp++] = savedFace[fp];
 805             }
 806         }
 807         newFace.setSize(newFp);
 808
 809         modifyFace(savedFaceLabels_[saveI], newFace, meshMod);
 810
 811         if (!hasSavedPoints)
 812         {
 813             // Face fully restored. Mark for compaction later on
 814             savedFaceLabels_[saveI] = -1;
 815             savedFaces_[saveI].clear();
 816         }
 817     }
 818
 819
 820     // Compact face labels
 821     label newSaveI = 0;
 822
 823     forAll(savedFaceLabels_, saveI)
 824     {
 825         if (savedFaceLabels_[saveI] != -1)
 826         {
 827             if (newSaveI != saveI)
 828             {
 829                 savedFaceLabels_[newSaveI] = savedFaceLabels_[saveI];
 830                 savedFaces_[newSaveI].transfer(savedFaces_[saveI]);
 831             }
 832             newSaveI++;
 833         }
 834     }
 835
 836     savedFaceLabels_.setSize(newSaveI);
 837     savedFaces_.setSize(newSaveI);
 838
 839
 840     // Check that all faces have been restored that use any restored points
 841     if (debug)
 842     {
 843         forAll(savedFaceLabels_, saveI)
 844         {
 845             const face& savedFace = savedFaces_[saveI];
 846
 847             forAll(savedFace, fp)
 848             {
 849                 if (savedFace[fp] < 0)
 850                 {
 851                     label addedPointI = addedPoints[-savedFace[fp]-1];
 852
 853                     if (addedPointI != -1)
 854                     {
 855                         FatalErrorIn("setUnrefinement")
 856                             << "Face:" << savedFaceLabels_[saveI]
 857                             << " savedVerts:" << savedFace
 858                             << " uses restored point:" << -savedFace[fp]-1
 859                             << " with new pointlabel:" << addedPointI
 860                             << abort(FatalError);
 861                     }
 862                 }
 863             }
 864         }
 865     }
 866 }
 867
 868
 869 // ************************************************************************* //