src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*----------------------------------------------------------------------------*/
  25
  26 #include "meshRefinement.H"
  27 #include "combineFaces.H"
  28 #include "polyTopoChange.H"
  29 #include "removePoints.H"
  30 #include "faceSet.H"
  31 #include "Time.H"
  32 #include "motionSmoother.H"
  33 #include "syncTools.H"
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36
  37 //// Merge faces that are in-line.
  38 //Foam::label Foam::meshRefinement::mergePatchFaces
  39 //(
  40 //    const scalar minCos,
  41 //    const scalar concaveCos,
  42 //    const labelList& patchIDs
  43 //)
  44 //{
  45 //    // Patch face merging engine
  46 //    combineFaces faceCombiner(mesh_);
  47 //
  48 //    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  49 //
  50 //    // Pick up all candidate cells on boundary
  51 //    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
  52 //
  53 //    forAll(patchIDs, i)
  54 //    {
  55 //        label patchI = patchIDs[i];
  56 //
  57 //        const polyPatch& patch = patches[patchI];
  58 //
  59 //        if (!patch.coupled())
  60 //        {
  61 //            forAll(patch, i)
  62 //            {
  63 //                boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
  64 //            }
  65 //        }
  66 //    }
  67 //
  68 //    // Get all sets of faces that can be merged
  69 //    labelListList mergeSets
  70 //    (
  71 //        faceCombiner.getMergeSets
  72 //        (
  73 //            minCos,
  74 //            concaveCos,
  75 //            boundaryCells
  76 //        )
  77 //    );
  78 //
  79 //    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
  80 //
  81 //    Info<< "mergePatchFaces : Merging " << nFaceSets
  82 //        << " sets of faces." << endl;
  83 //
  84 //    if (nFaceSets > 0)
  85 //    {
  86 //        // Topology changes container
  87 //        polyTopoChange meshMod(mesh_);
  88 //
  89 //        // Merge all faces of a set into the first face of the set. Remove
  90 //        // unused points.
  91 //        faceCombiner.setRefinement(mergeSets, meshMod);
  92 //
  93 //        // Change the mesh (no inflation)
  94 //        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
  95 //
  96 //        // Update fields
  97 //        mesh_.updateMesh(map);
  98 //
  99 //        // Move mesh (since morphing does not do this)
 100 //        if (map().hasMotionPoints())
 101 //        {
 102 //            mesh_.movePoints(map().preMotionPoints());
 103 //        }
 104 //        else
 105 //        {
 106 //            // Delete mesh volumes. No other way to do this?
 107 //            mesh_.clearOut();
 108 //        }
 109 //
 110 //
 111 //        // Reset the instance for if in overwrite mode
 112 //        mesh_.setInstance(timeName());
 113 //
 114 //        faceCombiner.updateMesh(map);
 115 //
 116 //        // Get the kept faces that need to be recalculated.
 117 //        // Merging two boundary faces might shift the cell centre
 118 //        // (unless the faces are absolutely planar)
 119 //        labelHashSet retestFaces(6*mergeSets.size());
 120 //
 121 //        forAll(mergeSets, setI)
 122 //        {
 123 //            label oldMasterI = mergeSets[setI][0];
 124 //
 125 //            label faceI = map().reverseFaceMap()[oldMasterI];
 126 //
 127 //            // faceI is always uncoupled boundary face
 128 //            const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
 129 //
 130 //            forAll(cFaces, i)
 131 //            {
 132 //                retestFaces.insert(cFaces[i]);
 133 //            }
 134 //        }
 135 //        updateMesh(map, retestFaces.toc());
 136 //    }
 137 //
 138 //
 139 //    return nFaceSets;
 140 //}
 141 //
 142 //
 143 //// Remove points not used by any face or points used by only two faces where
 144 //// the edges are in line
 145 //Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
 146 //(
 147 //    const scalar minCos
 148 //)
 149 //{
 150 //    // Point removal analysis engine
 151 //    removePoints pointRemover(mesh_);
 152 //
 153 //    // Count usage of points
 154 //    boolList pointCanBeDeleted;
 155 //    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
 156 //
 157 //    Info<< "Removing " << nRemove
 158 //        << " straight edge points." << endl;
 159 //
 160 //    autoPtr<mapPolyMesh> map;
 161 //
 162 //    if (nRemove > 0)
 163 //    {
 164 //        // Save my local faces that will change. These changed faces might
 165 //        // cause a shift in the cell centre which needs to be retested.
 166 //        // Have to do this before changing mesh since point will be removed.
 167 //        labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
 168 //
 169 //        {
 170 //            const faceList& faces = mesh_.faces();
 171 //
 172 //            forAll(faces, faceI)
 173 //            {
 174 //                const face& f = faces[faceI];
 175 //
 176 //                forAll(f, fp)
 177 //                {
 178 //                    if (pointCanBeDeleted[f[fp]])
 179 //                    {
 180 //                        retestOldFaces.insert(faceI);
 181 //                        break;
 182 //                    }
 183 //                }
 184 //            }
 185 //        }
 186 //
 187 //        // Topology changes container
 188 //        polyTopoChange meshMod(mesh_);
 189 //
 190 //        pointRemover.setRefinement(pointCanBeDeleted, meshMod);
 191 //
 192 //        // Change the mesh (no inflation)
 193 //        map = meshMod.changeMesh(mesh_, false, true);
 194 //
 195 //        // Update fields
 196 //        mesh_.updateMesh(map);
 197 //
 198 //        // Move mesh (since morphing does not do this)
 199 //        if (map().hasMotionPoints())
 200 //        {
 201 //            mesh_.movePoints(map().preMotionPoints());
 202 //        }
 203 //        else
 204 //        {
 205 //            // Delete mesh volumes. No other way to do this?
 206 //            mesh_.clearOut();
 207 //        }
 208 //
 209 //        // Reset the instance for if in overwrite mode
 210 //        mesh_.setInstance(timeName());
 211 //
 212 //        pointRemover.updateMesh(map);
 213 //
 214 //        // Get the kept faces that need to be recalculated.
 215 //        labelHashSet retestFaces(6*retestOldFaces.size());
 216 //
 217 //        const cellList& cells = mesh_.cells();
 218 //
 219 //        forAllConstIter(labelHashSet, retestOldFaces, iter)
 220 //        {
 221 //            label faceI = map().reverseFaceMap()[iter.key()];
 222 //
 223 //            const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
 224 //
 225 //            forAll(ownFaces, i)
 226 //            {
 227 //                retestFaces.insert(ownFaces[i]);
 228 //            }
 229 //
 230 //            if (mesh_.isInternalFace(faceI))
 231 //            {
 232 //                const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
 233 //
 234 //                forAll(neiFaces, i)
 235 //                {
 236 //                    retestFaces.insert(neiFaces[i]);
 237 //                }
 238 //            }
 239 //        }
 240 //        updateMesh(map, retestFaces.toc());
 241 //    }
 242 //
 243 //    return map;
 244 //}
 245
 246
 247 Foam::label Foam::meshRefinement::mergePatchFacesUndo
 248 (
 249     const scalar minCos,
 250     const scalar concaveCos,
 251     const labelList& patchIDs,
 252     const dictionary& motionDict
 253 )
 254 {
 255     // Patch face merging engine
 256     combineFaces faceCombiner(mesh_, true);
 257
 258     // Pick up all candidate cells on boundary
 259     labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
 260
 261     {
 262         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 263
 264         forAll(patchIDs, i)
 265         {
 266             label patchI = patchIDs[i];
 267
 268             const polyPatch& patch = patches[patchI];
 269
 270             if (!patch.coupled())
 271             {
 272                 forAll(patch, i)
 273                 {
 274                     boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
 275                 }
 276             }
 277         }
 278     }
 279
 280     // Get all sets of faces that can be merged. Since only faces on the same
 281     // patch get merged there is no risk of e.g. patchID faces getting merged
 282     // with original patches (or even processor patches). There is a risk
 283     // though of original patchfaces getting merged if they make a small
 284     // angle. Would be pretty weird starting mesh though.
 285     labelListList allFaceSets
 286     (
 287         faceCombiner.getMergeSets
 288         (
 289             minCos,
 290             concaveCos,
 291             boundaryCells
 292         )
 293     );
 294
 295     label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
 296
 297     Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
 298
 299     if (nFaceSets > 0)
 300     {
 301         if (debug)
 302         {
 303             faceSet allSets(mesh_, "allFaceSets", allFaceSets.size());
 304             forAll(allFaceSets, setI)
 305             {
 306                 forAll(allFaceSets[setI], i)
 307                 {
 308                     allSets.insert(allFaceSets[setI][i]);
 309                 }
 310             }
 311             Pout<< "Writing all faces to be merged to set "
 312                 << allSets.objectPath() << endl;
 313             allSets.instance() = timeName();
 314             allSets.write();
 315
 316             const_cast<Time&>(mesh_.time())++;
 317         }
 318
 319
 320         // Topology changes container
 321         polyTopoChange meshMod(mesh_);
 322
 323         // Merge all faces of a set into the first face of the set.
 324         faceCombiner.setRefinement(allFaceSets, meshMod);
 325
 326         // Experimental: store data for all the points that have been deleted
 327         storeData
 328         (
 329             faceCombiner.savedPointLabels(),    // points to store
 330             labelList(0),                       // faces to store
 331             labelList(0)                        // cells to store
 332         );
 333
 334         // Change the mesh (no inflation)
 335         autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 336
 337         // Update fields
 338         mesh_.updateMesh(map);
 339
 340         // Move mesh (since morphing does not do this)
 341         if (map().hasMotionPoints())
 342         {
 343             mesh_.movePoints(map().preMotionPoints());
 344         }
 345         else
 346         {
 347             // Delete mesh volumes.
 348             mesh_.clearOut();
 349         }
 350
 351         // Reset the instance for if in overwrite mode
 352         mesh_.setInstance(timeName());
 353
 354         faceCombiner.updateMesh(map);
 355
 356         // Get the kept faces that need to be recalculated.
 357         // Merging two boundary faces might shift the cell centre
 358         // (unless the faces are absolutely planar)
 359         labelHashSet retestFaces(2*allFaceSets.size());
 360
 361         forAll(allFaceSets, setI)
 362         {
 363             label oldMasterI = allFaceSets[setI][0];
 364             retestFaces.insert(map().reverseFaceMap()[oldMasterI]);
 365         }
 366         updateMesh(map, growFaceCellFace(retestFaces));
 367
 368         if (debug)
 369         {
 370             // Check sync
 371             Pout<< "Checking sync after initial merging " << nFaceSets
 372                 << " faces." << endl;
 373             checkData();
 374
 375             Pout<< "Writing initial merged-faces mesh to time "
 376                 << timeName() << nl << endl;
 377             write();
 378         }
 379
 380         for (label iteration = 0; iteration < 100; iteration++)
 381         {
 382             Info<< nl
 383                 << "Undo iteration " << iteration << nl
 384                 << "----------------" << endl;
 385
 386
 387             // Check mesh for errors
 388             // ~~~~~~~~~~~~~~~~~~~~~
 389
 390             faceSet errorFaces
 391             (
 392                 mesh_,
 393                 "errorFaces",
 394                 mesh_.nFaces()-mesh_.nInternalFaces()
 395             );
 396             bool hasErrors = motionSmoother::checkMesh
 397             (
 398                 false,  // report
 399                 mesh_,
 400                 motionDict,
 401                 errorFaces
 402             );
 403
 404             //if (checkEdgeConnectivity)
 405             //{
 406             //    Info<< "Checking edge-face connectivity (duplicate faces"
 407             //        << " or non-consecutive shared vertices)" << endl;
 408             //
 409             //    label nOldSize = errorFaces.size();
 410             //
 411             //    hasErrors =
 412             //        mesh_.checkFaceFaces
 413             //        (
 414             //            false,
 415             //            &errorFaces
 416             //        )
 417             //     || hasErrors;
 418             //
 419             //    Info<< "Detected additional "
 420             //        <<  returnReduce
 421             //            (
 422             //                errorFaces.size() - nOldSize,
 423             //                sumOp<label>()
 424             //            )
 425             //        << " faces with illegal face-face connectivity"
 426             //        << endl;
 427             //}
 428
 429             if (!hasErrors)
 430             {
 431                 break;
 432             }
 433
 434
 435             if (debug)
 436             {
 437                 errorFaces.instance() = timeName();
 438                 Pout<< "Writing all faces in error to faceSet "
 439                     << errorFaces.objectPath() << nl << endl;
 440                 errorFaces.write();
 441             }
 442
 443
 444             // Check any master cells for using any of the error faces
 445             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 446
 447             DynamicList<label> mastersToRestore(allFaceSets.size());
 448
 449             forAll(allFaceSets, setI)
 450             {
 451                 label masterFaceI = faceCombiner.masterFace()[setI];
 452
 453                 if (masterFaceI != -1)
 454                 {
 455                     label masterCellII = mesh_.faceOwner()[masterFaceI];
 456
 457                     const cell& cFaces = mesh_.cells()[masterCellII];
 458
 459                     forAll(cFaces, i)
 460                     {
 461                         if (errorFaces.found(cFaces[i]))
 462                         {
 463                             mastersToRestore.append(masterFaceI);
 464                             break;
 465                         }
 466                     }
 467                 }
 468             }
 469             mastersToRestore.shrink();
 470
 471             label nRestore = returnReduce
 472             (
 473                 mastersToRestore.size(),
 474                 sumOp<label>()
 475             );
 476
 477             Info<< "Masters that need to be restored:"
 478                 << nRestore << endl;
 479
 480             if (debug)
 481             {
 482                 faceSet restoreSet(mesh_, "mastersToRestore", mastersToRestore);
 483                 restoreSet.instance() = timeName();
 484                 Pout<< "Writing all " << mastersToRestore.size()
 485                     << " masterfaces to be restored to set "
 486                     << restoreSet.objectPath() << endl;
 487                 restoreSet.write();
 488             }
 489
 490
 491             if (nRestore == 0)
 492             {
 493                 break;
 494             }
 495
 496
 497             // Undo
 498             // ~~~~
 499
 500             if (debug)
 501             {
 502                 const_cast<Time&>(mesh_.time())++;
 503             }
 504
 505             // Topology changes container
 506             polyTopoChange meshMod(mesh_);
 507
 508             // Merge all faces of a set into the first face of the set.
 509             // Experimental:mark all points/faces/cells that have been restored.
 510             Map<label> restoredPoints(0);
 511             Map<label> restoredFaces(0);
 512             Map<label> restoredCells(0);
 513
 514             faceCombiner.setUnrefinement
 515             (
 516                 mastersToRestore,
 517                 meshMod,
 518                 restoredPoints,
 519                 restoredFaces,
 520                 restoredCells
 521             );
 522
 523             // Change the mesh (no inflation)
 524             autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 525
 526             // Update fields
 527             mesh_.updateMesh(map);
 528
 529             // Move mesh (since morphing does not do this)
 530             if (map().hasMotionPoints())
 531             {
 532                 mesh_.movePoints(map().preMotionPoints());
 533             }
 534             else
 535             {
 536                 // Delete mesh volumes.
 537                 mesh_.clearOut();
 538             }
 539
 540
 541             // Reset the instance for if in overwrite mode
 542             mesh_.setInstance(timeName());
 543
 544             faceCombiner.updateMesh(map);
 545
 546             // Renumber restore maps
 547             inplaceMapKey(map().reversePointMap(), restoredPoints);
 548             inplaceMapKey(map().reverseFaceMap(), restoredFaces);
 549             inplaceMapKey(map().reverseCellMap(), restoredCells);
 550
 551
 552             // Get the kept faces that need to be recalculated.
 553             // Merging two boundary faces might shift the cell centre
 554             // (unless the faces are absolutely planar)
 555             labelHashSet retestFaces(2*restoredFaces.size());
 556
 557             forAllConstIter(Map<label>, restoredFaces, iter)
 558             {
 559                 retestFaces.insert(iter.key());
 560             }
 561
 562             // Experimental:restore all points/face/cells in maps
 563             updateMesh
 564             (
 565                 map,
 566                 growFaceCellFace(retestFaces),
 567                 restoredPoints,
 568                 restoredFaces,
 569                 restoredCells
 570             );
 571
 572             if (debug)
 573             {
 574                 // Check sync
 575                 Pout<< "Checking sync after restoring " << retestFaces.size()
 576                     << " faces." << endl;
 577                 checkData();
 578
 579                 Pout<< "Writing merged-faces mesh to time "
 580                     << timeName() << nl << endl;
 581                 write();
 582             }
 583
 584             Info<< endl;
 585         }
 586     }
 587     else
 588     {
 589         Info<< "No faces merged ..." << endl;
 590     }
 591
 592     return nFaceSets;
 593 }
 594
 595
 596 // Remove points. pointCanBeDeleted is parallel synchronised.
 597 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemovePoints
 598 (
 599     removePoints& pointRemover,
 600     const boolList& pointCanBeDeleted
 601 )
 602 {
 603     // Topology changes container
 604     polyTopoChange meshMod(mesh_);
 605
 606     pointRemover.setRefinement(pointCanBeDeleted, meshMod);
 607
 608     // Change the mesh (no inflation)
 609     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 610
 611     // Update fields
 612     mesh_.updateMesh(map);
 613
 614     // Move mesh (since morphing does not do this)
 615     if (map().hasMotionPoints())
 616     {
 617         mesh_.movePoints(map().preMotionPoints());
 618     }
 619     else
 620     {
 621         // Delete mesh volumes.
 622         mesh_.clearOut();
 623     }
 624
 625     // Reset the instance for if in overwrite mode
 626     mesh_.setInstance(timeName());
 627
 628     pointRemover.updateMesh(map);
 629
 630
 631     // Retest all affected faces and all the cells using them
 632     labelHashSet retestFaces(pointRemover.savedFaceLabels().size());
 633     forAll(pointRemover.savedFaceLabels(), i)
 634     {
 635         label faceI = pointRemover.savedFaceLabels()[i];
 636         if (faceI >= 0)
 637         {
 638             retestFaces.insert(faceI);
 639         }
 640     }
 641     updateMesh(map, growFaceCellFace(retestFaces));
 642
 643     if (debug)
 644     {
 645         // Check sync
 646         Pout<< "Checking sync after removing points." << endl;
 647         checkData();
 648     }
 649
 650     return map;
 651 }
 652
 653
 654 // Restore faces (which contain removed points)
 655 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRestorePoints
 656 (
 657     removePoints& pointRemover,
 658     const labelList& facesToRestore
 659 )
 660 {
 661     // Topology changes container
 662     polyTopoChange meshMod(mesh_);
 663
 664     // Determine sets of points and faces to restore
 665     labelList localFaces, localPoints;
 666     pointRemover.getUnrefimentSet
 667     (
 668         facesToRestore,
 669         localFaces,
 670         localPoints
 671     );
 672
 673     // Undo the changes on the faces that are in error.
 674     pointRemover.setUnrefinement
 675     (
 676         localFaces,
 677         localPoints,
 678         meshMod
 679     );
 680
 681     // Change the mesh (no inflation)
 682     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 683
 684     // Update fields
 685     mesh_.updateMesh(map);
 686
 687     // Move mesh (since morphing does not do this)
 688     if (map().hasMotionPoints())
 689     {
 690         mesh_.movePoints(map().preMotionPoints());
 691     }
 692     else
 693     {
 694         // Delete mesh volumes.
 695         mesh_.clearOut();
 696     }
 697
 698     // Reset the instance for if in overwrite mode
 699     mesh_.setInstance(timeName());
 700
 701     pointRemover.updateMesh(map);
 702
 703     labelHashSet retestFaces(2*facesToRestore.size());
 704     forAll(facesToRestore, i)
 705     {
 706         label faceI = map().reverseFaceMap()[facesToRestore[i]];
 707         if (faceI >= 0)
 708         {
 709             retestFaces.insert(faceI);
 710         }
 711     }
 712     updateMesh(map, growFaceCellFace(retestFaces));
 713
 714     if (debug)
 715     {
 716         // Check sync
 717         Pout<< "Checking sync after restoring points on "
 718             << facesToRestore.size() << " faces." << endl;
 719         checkData();
 720     }
 721
 722     return map;
 723 }
 724
 725
 726 // Collect all faces that are both in candidateFaces and in the set.
 727 // If coupled face also collects the coupled face.
 728 Foam::labelList Foam::meshRefinement::collectFaces
 729 (
 730     const labelList& candidateFaces,
 731     const labelHashSet& set
 732 ) const
 733 {
 734     // Has face been selected?
 735     boolList selected(mesh_.nFaces(), false);
 736
 737     forAll(candidateFaces, i)
 738     {
 739         label faceI = candidateFaces[i];
 740
 741         if (set.found(faceI))
 742         {
 743             selected[faceI] = true;
 744         }
 745     }
 746     syncTools::syncFaceList
 747     (
 748         mesh_,
 749         selected,
 750         orEqOp<bool>()      // combine operator
 751     );
 752
 753     labelList selectedFaces(findIndices(selected, true));
 754
 755     return selectedFaces;
 756 }
 757
 758
 759 // Pick up faces of cells of faces in set.
 760 Foam::labelList Foam::meshRefinement::growFaceCellFace
 761 (
 762     const labelHashSet& set
 763 ) const
 764 {
 765     boolList selected(mesh_.nFaces(), false);
 766
 767     forAllConstIter(faceSet, set, iter)
 768     {
 769         label faceI = iter.key();
 770
 771         label own = mesh_.faceOwner()[faceI];
 772
 773         const cell& ownFaces = mesh_.cells()[own];
 774         forAll(ownFaces, ownFaceI)
 775         {
 776             selected[ownFaces[ownFaceI]] = true;
 777         }
 778
 779         if (mesh_.isInternalFace(faceI))
 780         {
 781             label nbr = mesh_.faceNeighbour()[faceI];
 782
 783             const cell& nbrFaces = mesh_.cells()[nbr];
 784             forAll(nbrFaces, nbrFaceI)
 785             {
 786                 selected[nbrFaces[nbrFaceI]] = true;
 787             }
 788         }
 789     }
 790     syncTools::syncFaceList
 791     (
 792         mesh_,
 793         selected,
 794         orEqOp<bool>()      // combine operator
 795     );
 796     return findIndices(selected, true);
 797 }
 798
 799
 800 // Remove points not used by any face or points used by only two faces where
 801 // the edges are in line
 802 Foam::label Foam::meshRefinement::mergeEdgesUndo
 803 (
 804     const scalar minCos,
 805     const dictionary& motionDict
 806 )
 807 {
 808     Info<< nl
 809         << "Merging all points on surface that" << nl
 810         << "- are used by only two boundary faces and" << nl
 811         << "- make an angle with a cosine of more than " << minCos
 812         << "." << nl << endl;
 813
 814     // Point removal analysis engine with undo
 815     removePoints pointRemover(mesh_, true);
 816
 817     // Count usage of points
 818     boolList pointCanBeDeleted;
 819     label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
 820
 821     if (nRemove > 0)
 822     {
 823         Info<< "Removing " << nRemove
 824             << " straight edge points ..." << nl << endl;
 825
 826         // Remove points
 827         // ~~~~~~~~~~~~~
 828
 829         doRemovePoints(pointRemover, pointCanBeDeleted);
 830
 831
 832         for (label iteration = 0; iteration < 100; iteration++)
 833         {
 834             Info<< nl
 835                 << "Undo iteration " << iteration << nl
 836                 << "----------------" << endl;
 837
 838
 839             // Check mesh for errors
 840             // ~~~~~~~~~~~~~~~~~~~~~
 841
 842             faceSet errorFaces
 843             (
 844                 mesh_,
 845                 "errorFaces",
 846                 mesh_.nFaces()-mesh_.nInternalFaces()
 847             );
 848             bool hasErrors = motionSmoother::checkMesh
 849             (
 850                 false,  // report
 851                 mesh_,
 852                 motionDict,
 853                 errorFaces
 854             );
 855             //if (checkEdgeConnectivity)
 856             //{
 857             //    Info<< "Checking edge-face connectivity (duplicate faces"
 858             //        << " or non-consecutive shared vertices)" << endl;
 859             //
 860             //    label nOldSize = errorFaces.size();
 861             //
 862             //    hasErrors =
 863             //        mesh_.checkFaceFaces
 864             //        (
 865             //            false,
 866             //            &errorFaces
 867             //        )
 868             //     || hasErrors;
 869             //
 870             //    Info<< "Detected additional "
 871             //        << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
 872             //        << " faces with illegal face-face connectivity"
 873             //        << endl;
 874             //}
 875
 876             if (!hasErrors)
 877             {
 878                 break;
 879             }
 880
 881             if (debug)
 882             {
 883                 errorFaces.instance() = timeName();
 884                 Pout<< "**Writing all faces in error to faceSet "
 885                     << errorFaces.objectPath() << nl << endl;
 886                 errorFaces.write();
 887             }
 888
 889             labelList masterErrorFaces
 890             (
 891                 collectFaces
 892                 (
 893                     pointRemover.savedFaceLabels(),
 894                     errorFaces
 895                 )
 896             );
 897
 898             label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
 899
 900             Info<< "Detected " << n
 901                 << " error faces on boundaries that have been merged."
 902                 << " These will be restored to their original faces." << nl
 903                 << endl;
 904
 905             if (n == 0)
 906             {
 907                 if (hasErrors)
 908                 {
 909                     Info<< "Detected "
 910                         << returnReduce(errorFaces.size(), sumOp<label>())
 911                         << " error faces in mesh."
 912                         << " Restoring neighbours of faces in error." << nl
 913                         << endl;
 914
 915                     labelList expandedErrorFaces
 916                     (
 917                         growFaceCellFace
 918                         (
 919                             errorFaces
 920                         )
 921                     );
 922
 923                     doRestorePoints(pointRemover, expandedErrorFaces);
 924                 }
 925
 926                 break;
 927             }
 928
 929             doRestorePoints(pointRemover, masterErrorFaces);
 930         }
 931
 932         if (debug)
 933         {
 934             const_cast<Time&>(mesh_.time())++;
 935             Pout<< "Writing merged-edges mesh to time "
 936                 << timeName() << nl << endl;
 937             write();
 938         }
 939     }
 940     else
 941     {
 942         Info<< "No straight edges simplified and no points removed ..." << endl;
 943     }
 944
 945     return nRemove;
 946 }
 947
 948
 949 // ************************************************************************* //