src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     All to do with adding cell layers
  27
  28 \*----------------------------------------------------------------------------*/
  29
  30 #include "autoLayerDriver.H"
  31 #include "fvMesh.H"
  32 #include "Time.H"
  33 #include "meshRefinement.H"
  34 #include "removePoints.H"
  35 #include "pointFields.H"
  36 #include "motionSmoother.H"
  37 #include "mathematicalConstants.H"
  38 #include "pointSet.H"
  39 #include "faceSet.H"
  40 #include "cellSet.H"
  41 #include "polyTopoChange.H"
  42 #include "mapPolyMesh.H"
  43 #include "addPatchCellLayer.H"
  44 #include "mapDistributePolyMesh.H"
  45 #include "OFstream.H"
  46 #include "layerParameters.H"
  47 #include "combineFaces.H"
  48
  49 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  50
  51 namespace Foam
  52 {
  53
  54 defineTypeNameAndDebug(autoLayerDriver, 0);
  55
  56 } // End namespace Foam
  57
  58
  59 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  60
  61 Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
  62 (
  63     const scalar minCos,
  64     const scalar concaveCos,
  65     const dictionary& motionDict
  66 )
  67 {
  68     fvMesh& mesh = meshRefiner_.mesh();
  69
  70     // Patch face merging engine
  71     combineFaces faceCombiner(mesh, true);
  72
  73     // Pick up all candidate cells on boundary
  74     labelHashSet boundaryCells(mesh.nFaces()-mesh.nInternalFaces());
  75
  76     {
  77         labelList patchIDs(meshRefinement::addedPatches(globalToPatch_));
  78
  79         const polyBoundaryMesh& patches = mesh.boundaryMesh();
  80
  81         forAll(patchIDs, i)
  82         {
  83             label patchI = patchIDs[i];
  84
  85             const polyPatch& patch = patches[patchI];
  86
  87             if (!patch.coupled())
  88             {
  89                 forAll(patch, i)
  90                 {
  91                     boundaryCells.insert(mesh.faceOwner()[patch.start()+i]);
  92                 }
  93             }
  94         }
  95     }
  96
  97     // Get all sets of faces that can be merged
  98     labelListList allFaceSets
  99     (
 100         faceCombiner.getMergeSets
 101         (
 102             minCos,
 103             concaveCos,
 104             boundaryCells
 105         )
 106     );
 107
 108     label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
 109
 110     Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
 111
 112     if (nFaceSets > 0)
 113     {
 114         if (debug)
 115         {
 116             faceSet allSets(mesh, "allFaceSets", allFaceSets.size());
 117             forAll(allFaceSets, setI)
 118             {
 119                 forAll(allFaceSets[setI], i)
 120                 {
 121                     allSets.insert(allFaceSets[setI][i]);
 122                 }
 123             }
 124             Pout<< "Writing all faces to be merged to set "
 125                 << allSets.objectPath() << endl;
 126             allSets.write();
 127         }
 128
 129
 130         // Topology changes container
 131         polyTopoChange meshMod(mesh);
 132
 133         // Merge all faces of a set into the first face of the set.
 134         faceCombiner.setRefinement(allFaceSets, meshMod);
 135
 136         // Experimental: store data for all the points that have been deleted
 137         meshRefiner_.storeData
 138         (
 139             faceCombiner.savedPointLabels(),    // points to store
 140             labelList(0),                       // faces to store
 141             labelList(0)                        // cells to store
 142         );
 143
 144         // Change the mesh (no inflation)
 145         autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
 146
 147         // Update fields
 148         mesh.updateMesh(map);
 149
 150         // Move mesh (since morphing does not do this)
 151         if (map().hasMotionPoints())
 152         {
 153             mesh.movePoints(map().preMotionPoints());
 154         }
 155
 156         faceCombiner.updateMesh(map);
 157
 158         meshRefiner_.updateMesh(map, labelList(0));
 159
 160
 161
 162         for (label iteration = 0; iteration < 100; iteration++)
 163         {
 164             Info<< nl
 165                 << "Undo iteration " << iteration << nl
 166                 << "----------------" << endl;
 167
 168
 169             // Check mesh for errors
 170             // ~~~~~~~~~~~~~~~~~~~~~
 171
 172             faceSet errorFaces
 173             (
 174                 mesh,
 175                 "errorFaces",
 176                 mesh.nFaces()-mesh.nInternalFaces()
 177             );
 178             bool hasErrors = motionSmoother::checkMesh
 179             (
 180                 false,  // report
 181                 mesh,
 182                 motionDict,
 183                 errorFaces
 184             );
 185             //if (checkEdgeConnectivity)
 186             //{
 187             //    Info<< "Checking edge-face connectivity (duplicate faces"
 188             //        << " or non-consecutive shared vertices)" << endl;
 189             //
 190             //    label nOldSize = errorFaces.size();
 191             //
 192             //    hasErrors =
 193             //        mesh.checkFaceFaces
 194             //        (
 195             //            false,
 196             //            &errorFaces
 197             //        )
 198             //     || hasErrors;
 199             //
 200             //    Info<< "Detected additional "
 201             //        << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
 202             //        << " faces with illegal face-face connectivity"
 203             //        << endl;
 204             //}
 205
 206             if (!hasErrors)
 207             {
 208                 break;
 209             }
 210
 211
 212             if (debug)
 213             {
 214                 Pout<< "Writing all faces in error to faceSet "
 215                     << errorFaces.objectPath() << nl << endl;
 216                 errorFaces.write();
 217             }
 218
 219
 220             // Check any master cells for using any of the error faces
 221             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 222
 223             DynamicList<label> mastersToRestore(allFaceSets.size());
 224
 225             forAll(allFaceSets, setI)
 226             {
 227                 label masterFaceI = faceCombiner.masterFace()[setI];
 228
 229                 if (masterFaceI != -1)
 230                 {
 231                     label masterCellII = mesh.faceOwner()[masterFaceI];
 232
 233                     const cell& cFaces = mesh.cells()[masterCellII];
 234
 235                     forAll(cFaces, i)
 236                     {
 237                         if (errorFaces.found(cFaces[i]))
 238                         {
 239                             mastersToRestore.append(masterFaceI);
 240                             break;
 241                         }
 242                     }
 243                 }
 244             }
 245             mastersToRestore.shrink();
 246
 247             label nRestore = returnReduce
 248             (
 249                 mastersToRestore.size(),
 250                 sumOp<label>()
 251             );
 252
 253             Info<< "Masters that need to be restored:"
 254                 << nRestore << endl;
 255
 256             if (debug)
 257             {
 258                 faceSet restoreSet(mesh, "mastersToRestore", mastersToRestore);
 259                 Pout<< "Writing all " << mastersToRestore.size()
 260                     << " masterfaces to be restored to set "
 261                     << restoreSet.objectPath() << endl;
 262                 restoreSet.write();
 263             }
 264
 265
 266             if (nRestore == 0)
 267             {
 268                 break;
 269             }
 270
 271
 272             // Undo
 273             // ~~~~
 274
 275             // Topology changes container
 276             polyTopoChange meshMod(mesh);
 277
 278             // Merge all faces of a set into the first face of the set.
 279             // Experimental:mark all points/faces/cells that have been restored.
 280             Map<label> restoredPoints(0);
 281             Map<label> restoredFaces(0);
 282             Map<label> restoredCells(0);
 283
 284             faceCombiner.setUnrefinement
 285             (
 286                 mastersToRestore,
 287                 meshMod,
 288                 restoredPoints,
 289                 restoredFaces,
 290                 restoredCells
 291             );
 292
 293             // Change the mesh (no inflation)
 294             autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
 295
 296             // Update fields
 297             mesh.updateMesh(map);
 298
 299             // Move mesh (since morphing does not do this)
 300             if (map().hasMotionPoints())
 301             {
 302                 mesh.movePoints(map().preMotionPoints());
 303             }
 304
 305             faceCombiner.updateMesh(map);
 306
 307             // Renumber restore maps
 308             inplaceMapKey(map().reversePointMap(), restoredPoints);
 309             inplaceMapKey(map().reverseFaceMap(), restoredFaces);
 310             inplaceMapKey(map().reverseCellMap(), restoredCells);
 311
 312             // Experimental:restore all points/face/cells in maps
 313             meshRefiner_.updateMesh
 314             (
 315                 map,
 316                 labelList(0),       // changedFaces
 317                 restoredPoints,
 318                 restoredFaces,
 319                 restoredCells
 320             );
 321
 322             Info<< endl;
 323         }
 324
 325         if (debug)
 326         {
 327             Pout<< "Writing merged-faces mesh to time "
 328                 << mesh.time().timeName() << nl << endl;
 329             mesh.write();
 330         }
 331     }
 332     else
 333     {
 334         Info<< "No faces merged ..." << endl;
 335     }
 336
 337     return nFaceSets;
 338 }
 339
 340
 341 // Remove points. pointCanBeDeleted is parallel synchronised.
 342 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
 343 (
 344     removePoints& pointRemover,
 345     const boolList& pointCanBeDeleted
 346 )
 347 {
 348     fvMesh& mesh = meshRefiner_.mesh();
 349
 350     // Topology changes container
 351     polyTopoChange meshMod(mesh);
 352
 353     pointRemover.setRefinement(pointCanBeDeleted, meshMod);
 354
 355     // Change the mesh (no inflation)
 356     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
 357
 358     // Update fields
 359     mesh.updateMesh(map);
 360
 361     // Move mesh (since morphing does not do this)
 362     if (map().hasMotionPoints())
 363     {
 364         mesh.movePoints(map().preMotionPoints());
 365     }
 366
 367     pointRemover.updateMesh(map);
 368     meshRefiner_.updateMesh(map, labelList(0));
 369
 370     return map;
 371 }
 372
 373
 374 // Restore faces (which contain removed points)
 375 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
 376 (
 377     removePoints& pointRemover,
 378     const labelList& facesToRestore
 379 )
 380 {
 381     fvMesh& mesh = meshRefiner_.mesh();
 382
 383     // Topology changes container
 384     polyTopoChange meshMod(mesh);
 385
 386     // Determine sets of points and faces to restore
 387     labelList localFaces, localPoints;
 388     pointRemover.getUnrefimentSet
 389     (
 390         facesToRestore,
 391         localFaces,
 392         localPoints
 393     );
 394
 395     // Undo the changes on the faces that are in error.
 396     pointRemover.setUnrefinement
 397     (
 398         localFaces,
 399         localPoints,
 400         meshMod
 401     );
 402
 403     // Change the mesh (no inflation)
 404     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
 405
 406     // Update fields
 407     mesh.updateMesh(map);
 408
 409     // Move mesh (since morphing does not do this)
 410     if (map().hasMotionPoints())
 411     {
 412         mesh.movePoints(map().preMotionPoints());
 413     }
 414
 415     pointRemover.updateMesh(map);
 416     meshRefiner_.updateMesh(map, labelList(0));
 417
 418     return map;
 419 }
 420
 421
 422 // Collect all faces that are both in candidateFaces and in the set.
 423 // If coupled face also collects the coupled face.
 424 Foam::labelList Foam::autoLayerDriver::collectFaces
 425 (
 426     const labelList& candidateFaces,
 427     const labelHashSet& set
 428 ) const
 429 {
 430     const fvMesh& mesh = meshRefiner_.mesh();
 431
 432     // Has face been selected?
 433     boolList selected(mesh.nFaces(), false);
 434
 435     forAll(candidateFaces, i)
 436     {
 437         label faceI = candidateFaces[i];
 438
 439         if (set.found(faceI))
 440         {
 441             selected[faceI] = true;
 442         }
 443     }
 444     syncTools::syncFaceList
 445     (
 446         mesh,
 447         selected,
 448         orEqOp<bool>(),     // combine operator
 449         false               // separation
 450     );
 451
 452     labelList selectedFaces(findIndices(selected, true));
 453
 454     return selectedFaces;
 455 }
 456
 457
 458 // Pick up faces of cells of faces in set.
 459 Foam::labelList Foam::autoLayerDriver::growFaceCellFace
 460 (
 461     const labelHashSet& set
 462 ) const
 463 {
 464     const fvMesh& mesh = meshRefiner_.mesh();
 465
 466     boolList selected(mesh.nFaces(), false);
 467
 468     forAllConstIter(faceSet, set, iter)
 469     {
 470         label faceI = iter.key();
 471
 472         label own = mesh.faceOwner()[faceI];
 473
 474         const cell& ownFaces = mesh.cells()[own];
 475         forAll(ownFaces, ownFaceI)
 476         {
 477             selected[ownFaces[ownFaceI]] = true;
 478         }
 479
 480         if (mesh.isInternalFace(faceI))
 481         {
 482             label nbr = mesh.faceNeighbour()[faceI];
 483
 484             const cell& nbrFaces = mesh.cells()[nbr];
 485             forAll(nbrFaces, nbrFaceI)
 486             {
 487                 selected[nbrFaces[nbrFaceI]] = true;
 488             }
 489         }
 490     }
 491     syncTools::syncFaceList
 492     (
 493         mesh,
 494         selected,
 495         orEqOp<bool>(),     // combine operator
 496         false               // separation
 497     );
 498     return findIndices(selected, true);
 499 }
 500
 501
 502 // Remove points not used by any face or points used by only two faces where
 503 // the edges are in line
 504 Foam::label Foam::autoLayerDriver::mergeEdgesUndo
 505 (
 506     const scalar minCos,
 507     const dictionary& motionDict
 508 )
 509 {
 510     fvMesh& mesh = meshRefiner_.mesh();
 511
 512     Info<< nl
 513         << "Merging all points on surface that" << nl
 514         << "- are used by only two boundary faces and" << nl
 515         << "- make an angle with a cosine of more than " << minCos
 516         << "." << nl << endl;
 517
 518     // Point removal analysis engine with undo
 519     removePoints pointRemover(mesh, true);
 520
 521     // Count usage of points
 522     boolList pointCanBeDeleted;
 523     label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
 524
 525     if (nRemove > 0)
 526     {
 527         Info<< "Removing " << nRemove
 528             << " straight edge points ..." << nl << endl;
 529
 530         // Remove points
 531         // ~~~~~~~~~~~~~
 532
 533         doRemovePoints(pointRemover, pointCanBeDeleted);
 534
 535
 536         for (label iteration = 0; iteration < 100; iteration++)
 537         {
 538             Info<< nl
 539                 << "Undo iteration " << iteration << nl
 540                 << "----------------" << endl;
 541
 542
 543             // Check mesh for errors
 544             // ~~~~~~~~~~~~~~~~~~~~~
 545
 546             faceSet errorFaces
 547             (
 548                 mesh,
 549                 "errorFaces",
 550                 mesh.nFaces()-mesh.nInternalFaces()
 551             );
 552             bool hasErrors = motionSmoother::checkMesh
 553             (
 554                 false,  // report
 555                 mesh,
 556                 motionDict,
 557                 errorFaces
 558             );
 559             //if (checkEdgeConnectivity)
 560             //{
 561             //    Info<< "Checking edge-face connectivity (duplicate faces"
 562             //        << " or non-consecutive shared vertices)" << endl;
 563             //
 564             //    label nOldSize = errorFaces.size();
 565             //
 566             //    hasErrors =
 567             //        mesh.checkFaceFaces
 568             //        (
 569             //            false,
 570             //            &errorFaces
 571             //        )
 572             //     || hasErrors;
 573             //
 574             //    Info<< "Detected additional "
 575             //        << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
 576             //        << " faces with illegal face-face connectivity"
 577             //        << endl;
 578             //}
 579
 580             if (!hasErrors)
 581             {
 582                 break;
 583             }
 584
 585             if (debug)
 586             {
 587                 Pout<< "**Writing all faces in error to faceSet "
 588                     << errorFaces.objectPath() << nl << endl;
 589                 errorFaces.write();
 590             }
 591
 592             labelList masterErrorFaces
 593             (
 594                 collectFaces
 595                 (
 596                     pointRemover.savedFaceLabels(),
 597                     errorFaces
 598                 )
 599             );
 600
 601             label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
 602
 603             Info<< "Detected " << n
 604                 << " error faces on boundaries that have been merged."
 605                 << " These will be restored to their original faces." << nl
 606                 << endl;
 607
 608             if (n == 0)
 609             {
 610                 if (hasErrors)
 611                 {
 612                     Info<< "Detected "
 613                         << returnReduce(errorFaces.size(), sumOp<label>())
 614                         << " error faces in mesh."
 615                         << " Restoring neighbours of faces in error." << nl
 616                         << endl;
 617
 618                     labelList expandedErrorFaces
 619                     (
 620                         growFaceCellFace
 621                         (
 622                             errorFaces
 623                         )
 624                     );
 625
 626                     doRestorePoints(pointRemover, expandedErrorFaces);
 627                 }
 628
 629                 break;
 630             }
 631
 632             doRestorePoints(pointRemover, masterErrorFaces);
 633         }
 634
 635         if (debug)
 636         {
 637             Pout<< "Writing merged-edges mesh to time "
 638                 << mesh.time().timeName() << nl << endl;
 639             mesh.write();
 640         }
 641     }
 642     else
 643     {
 644         Info<< "No straight edges simplified and no points removed ..." << endl;
 645     }
 646
 647     return nRemove;
 648 }
 649
 650
 651 // For debugging: Dump displacement to .obj files
 652 void Foam::autoLayerDriver::dumpDisplacement
 653 (
 654     const fileName& prefix,
 655     const indirectPrimitivePatch& pp,
 656     const vectorField& patchDisp,
 657     const List<extrudeMode>& extrudeStatus
 658 )
 659 {
 660     OFstream dispStr(prefix + "_disp.obj");
 661     Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
 662
 663     label vertI = 0;
 664
 665     forAll(patchDisp, patchPointI)
 666     {
 667         const point& pt = pp.localPoints()[patchPointI];
 668
 669         meshTools::writeOBJ(dispStr, pt); vertI++;
 670         meshTools::writeOBJ(dispStr, pt + patchDisp[patchPointI]); vertI++;
 671
 672         dispStr << "l " << vertI-1 << ' ' << vertI << nl;
 673     }
 674
 675
 676     OFstream illStr(prefix + "_illegal.obj");
 677     Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
 678
 679     vertI = 0;
 680
 681     forAll(patchDisp, patchPointI)
 682     {
 683         if (extrudeStatus[patchPointI] != EXTRUDE)
 684         {
 685             const point& pt = pp.localPoints()[patchPointI];
 686
 687             meshTools::writeOBJ(illStr, pt); vertI++;
 688             meshTools::writeOBJ(illStr, pt + patchDisp[patchPointI]); vertI++;
 689
 690             illStr << "l " << vertI-1 << ' ' << vertI << nl;
 691         }
 692     }
 693 }
 694
 695
 696 // Check that primitivePatch is not multiply connected. Collect non-manifold
 697 // points in pointSet.
 698 void Foam::autoLayerDriver::checkManifold
 699 (
 700     const indirectPrimitivePatch& fp,
 701     pointSet& nonManifoldPoints
 702 )
 703 {
 704     // Check for non-manifold points (surface pinched at point)
 705     fp.checkPointManifold(false, &nonManifoldPoints);
 706
 707     // Check for edge-faces (surface pinched at edge)
 708     const labelListList& edgeFaces = fp.edgeFaces();
 709
 710     forAll(edgeFaces, edgeI)
 711     {
 712         const labelList& eFaces = edgeFaces[edgeI];
 713
 714         if (eFaces.size() > 2)
 715         {
 716             const edge& e = fp.edges()[edgeI];
 717
 718             nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
 719             nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
 720         }
 721     }
 722 }
 723
 724
 725 void Foam::autoLayerDriver::checkMeshManifold() const
 726 {
 727     const fvMesh& mesh = meshRefiner_.mesh();
 728
 729     Info<< nl << "Checking mesh manifoldness ..." << endl;
 730
 731     // Get all outside faces
 732     labelList outsideFaces(mesh.nFaces() - mesh.nInternalFaces());
 733
 734     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 735     {
 736          outsideFaces[faceI - mesh.nInternalFaces()] = faceI;
 737     }
 738
 739     pointSet nonManifoldPoints
 740     (
 741         mesh,
 742         "nonManifoldPoints",
 743         mesh.nPoints() / 100
 744     );
 745
 746     // Build primitivePatch out of faces and check it for problems.
 747     checkManifold
 748     (
 749         indirectPrimitivePatch
 750         (
 751             IndirectList<face>(mesh.faces(), outsideFaces),
 752             mesh.points()
 753         ),
 754         nonManifoldPoints
 755     );
 756
 757     label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
 758
 759     if (nNonManif > 0)
 760     {
 761         Info<< "Outside of mesh is multiply connected across edges or"
 762             << " points." << nl
 763             << "This is not a fatal error but might cause some unexpected"
 764             << " behaviour." << nl
 765             << "Writing " << nNonManif
 766             << " points where this happens to pointSet "
 767             << nonManifoldPoints.name() << endl;
 768
 769         nonManifoldPoints.write();
 770     }
 771     Info<< endl;
 772 }
 773
 774
 775
 776 // Unset extrusion on point. Returns true if anything unset.
 777 bool Foam::autoLayerDriver::unmarkExtrusion
 778 (
 779     const label patchPointI,
 780     pointField& patchDisp,
 781     labelList& patchNLayers,
 782     List<extrudeMode>& extrudeStatus
 783 )
 784 {
 785     if (extrudeStatus[patchPointI] == EXTRUDE)
 786     {
 787         extrudeStatus[patchPointI] = NOEXTRUDE;
 788         patchNLayers[patchPointI] = 0;
 789         patchDisp[patchPointI] = vector::zero;
 790         return true;
 791     }
 792     else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
 793     {
 794         extrudeStatus[patchPointI] = NOEXTRUDE;
 795         patchNLayers[patchPointI] = 0;
 796         patchDisp[patchPointI] = vector::zero;
 797         return true;
 798     }
 799     else
 800     {
 801         return false;
 802     }
 803 }
 804
 805
 806 // Unset extrusion on face. Returns true if anything unset.
 807 bool Foam::autoLayerDriver::unmarkExtrusion
 808 (
 809     const face& localFace,
 810     pointField& patchDisp,
 811     labelList& patchNLayers,
 812     List<extrudeMode>& extrudeStatus
 813 )
 814 {
 815     bool unextruded = false;
 816
 817     forAll(localFace, fp)
 818     {
 819         if
 820         (
 821             unmarkExtrusion
 822             (
 823                 localFace[fp],
 824                 patchDisp,
 825                 patchNLayers,
 826                 extrudeStatus
 827             )
 828         )
 829         {
 830             unextruded = true;
 831         }
 832     }
 833     return unextruded;
 834 }
 835
 836
 837 // No extrusion at non-manifold points.
 838 void Foam::autoLayerDriver::handleNonManifolds
 839 (
 840     const indirectPrimitivePatch& pp,
 841     const labelList& meshEdges,
 842     pointField& patchDisp,
 843     labelList& patchNLayers,
 844     List<extrudeMode>& extrudeStatus
 845 ) const
 846 {
 847     const fvMesh& mesh = meshRefiner_.mesh();
 848
 849     Info<< nl << "Handling non-manifold points ..." << endl;
 850
 851     // Detect non-manifold points
 852     Info<< nl << "Checking patch manifoldness ..." << endl;
 853
 854     pointSet nonManifoldPoints(mesh, "nonManifoldPoints", pp.nPoints());
 855
 856     // 1. Local check
 857     checkManifold(pp, nonManifoldPoints);
 858
 859     label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
 860
 861     Info<< "Outside of local patch is multiply connected across edges or"
 862         << " points at " << nNonManif << " points." << endl;
 863
 864
 865     const labelList& meshPoints = pp.meshPoints();
 866
 867
 868     // 2. Parallel check
 869     // For now disable extrusion at any shared edges.
 870     const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels());
 871
 872     forAll(pp.edges(), edgeI)
 873     {
 874         if (sharedEdgeSet.found(meshEdges[edgeI]))
 875         {
 876             const edge& e = mesh.edges()[meshEdges[edgeI]];
 877
 878             Pout<< "Disabling extrusion at edge "
 879                 << mesh.points()[e[0]] << mesh.points()[e[1]]
 880                 << " since it is non-manifold across coupled patches."
 881                 << endl;
 882
 883             nonManifoldPoints.insert(e[0]);
 884             nonManifoldPoints.insert(e[1]);
 885         }
 886     }
 887
 888     // 3b. extrusion can produce multiple faces between 2 cells
 889     // across processor boundary
 890     // This occurs when a coupled face shares more than 1 edge with a
 891     // non-coupled boundary face.
 892     // This is now correctly handled by addPatchCellLayer in that it
 893     // extrudes a single face from the stringed up edges.
 894
 895
 896     nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
 897
 898     if (nNonManif > 0)
 899     {
 900         Info<< "Outside of patches is multiply connected across edges or"
 901             << " points at " << nNonManif << " points." << nl
 902             << "Writing " << nNonManif
 903             << " points where this happens to pointSet "
 904             << nonManifoldPoints.name() << nl
 905             << "and setting layer thickness to zero on these points."
 906             << endl;
 907
 908         nonManifoldPoints.write();
 909
 910         forAll(meshPoints, patchPointI)
 911         {
 912             if (nonManifoldPoints.found(meshPoints[patchPointI]))
 913             {
 914                 unmarkExtrusion
 915                 (
 916                     patchPointI,
 917                     patchDisp,
 918                     patchNLayers,
 919                     extrudeStatus
 920                 );
 921             }
 922         }
 923     }
 924
 925     Info<< "Set displacement to zero for all " << nNonManif
 926         << " non-manifold points" << endl;
 927 }
 928
 929
 930 // Parallel feature edge detection. Assumes non-manifold edges already handled.
 931 void Foam::autoLayerDriver::handleFeatureAngle
 932 (
 933     const indirectPrimitivePatch& pp,
 934     const labelList& meshEdges,
 935     const scalar minCos,
 936     pointField& patchDisp,
 937     labelList& patchNLayers,
 938     List<extrudeMode>& extrudeStatus
 939 ) const
 940 {
 941     const fvMesh& mesh = meshRefiner_.mesh();
 942
 943     Info<< nl << "Handling feature edges ..." << endl;
 944
 945     if (minCos < 1-SMALL)
 946     {
 947         // Normal component of normals of connected faces.
 948         vectorField edgeNormal(mesh.nEdges(), wallPoint::greatPoint);
 949
 950         const labelListList& edgeFaces = pp.edgeFaces();
 951
 952         forAll(edgeFaces, edgeI)
 953         {
 954             const labelList& eFaces = pp.edgeFaces()[edgeI];
 955
 956             label meshEdgeI = meshEdges[edgeI];
 957
 958             forAll(eFaces, i)
 959             {
 960                 nomalsCombine()
 961                 (
 962                     edgeNormal[meshEdgeI],
 963                     pp.faceNormals()[eFaces[i]]
 964                 );
 965             }
 966         }
 967
 968         syncTools::syncEdgeList
 969         (
 970             mesh,
 971             edgeNormal,
 972             nomalsCombine(),
 973             wallPoint::greatPoint,  // null value
 974             false                   // no separation
 975         );
 976
 977         label vertI = 0;
 978         autoPtr<OFstream> str;
 979         if (debug)
 980         {
 981             str.reset(new OFstream(mesh.time().path()/"featureEdges.obj"));
 982             Info<< "Writing feature edges to " << str().name() << endl;
 983         }
 984
 985         label nFeats = 0;
 986
 987         // Now on coupled edges the edgeNormal will have been truncated and
 988         // only be still be the old value where two faces have the same normal
 989         forAll(edgeFaces, edgeI)
 990         {
 991             const labelList& eFaces = pp.edgeFaces()[edgeI];
 992
 993             label meshEdgeI = meshEdges[edgeI];
 994
 995             const vector& n = edgeNormal[meshEdgeI];
 996
 997             if (n != wallPoint::greatPoint)
 998             {
 999                 scalar cos = n & pp.faceNormals()[eFaces[0]];
1000
1001                 if (cos < minCos)
1002                 {
1003                     const edge& e = pp.edges()[edgeI];
1004
1005                     unmarkExtrusion
1006                     (
1007                         e[0],
1008                         patchDisp,
1009                         patchNLayers,
1010                         extrudeStatus
1011                     );
1012                     unmarkExtrusion
1013                     (
1014                         e[1],
1015                         patchDisp,
1016                         patchNLayers,
1017                         extrudeStatus
1018                     );
1019
1020                     nFeats++;
1021
1022                     if (str.valid())
1023                     {
1024                         meshTools::writeOBJ(str(), pp.localPoints()[e[0]]);
1025                         vertI++;
1026                         meshTools::writeOBJ(str(), pp.localPoints()[e[1]]);
1027                         vertI++;
1028                         str()<< "l " << vertI-1 << ' ' << vertI << nl;
1029                     }
1030                 }
1031             }
1032         }
1033
1034         Info<< "Set displacement to zero for points on "
1035             << returnReduce(nFeats, sumOp<label>())
1036             << " feature edges" << endl;
1037     }
1038 }
1039
1040
1041 // No extrusion on cells with warped faces. Calculates the thickness of the
1042 // layer and compares it to the space the warped face takes up. Disables
1043 // extrusion if layer thickness is more than faceRatio of the thickness of
1044 // the face.
1045 void Foam::autoLayerDriver::handleWarpedFaces
1046 (
1047     const indirectPrimitivePatch& pp,
1048     const scalar faceRatio,
1049     const scalar edge0Len,
1050     const labelList& cellLevel,
1051     pointField& patchDisp,
1052     labelList& patchNLayers,
1053     List<extrudeMode>& extrudeStatus
1054 ) const
1055 {
1056     const fvMesh& mesh = meshRefiner_.mesh();
1057
1058     Info<< nl << "Handling cells with warped patch faces ..." << nl;
1059
1060     const pointField& points = mesh.points();
1061
1062     label nWarpedFaces = 0;
1063
1064     forAll(pp, i)
1065     {
1066         const face& f = pp[i];
1067
1068         if (f.size() > 3)
1069         {
1070             label faceI = pp.addressing()[i];
1071
1072             label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
1073             scalar edgeLen = edge0Len/(1<<ownLevel);
1074
1075             // Normal distance to face centre plane
1076             const point& fc = mesh.faceCentres()[faceI];
1077             const vector& fn = pp.faceNormals()[i];
1078
1079             scalarField vProj(f.size());
1080
1081             forAll(f, fp)
1082             {
1083                 vector n = points[f[fp]] - fc;
1084                 vProj[fp] = (n & fn);
1085             }
1086
1087             // Get normal 'span' of face
1088             scalar minVal = min(vProj);
1089             scalar maxVal = max(vProj);
1090
1091             if ((maxVal - minVal) > faceRatio * edgeLen)
1092             {
1093                 if
1094                 (
1095                     unmarkExtrusion
1096                     (
1097                         pp.localFaces()[i],
1098                         patchDisp,
1099                         patchNLayers,
1100                         extrudeStatus
1101                     )
1102                 )
1103                 {
1104                     nWarpedFaces++;
1105                 }
1106             }
1107         }
1108     }
1109
1110     Info<< "Set displacement to zero on "
1111         << returnReduce(nWarpedFaces, sumOp<label>())
1112         << " warped faces since layer would be > " << faceRatio
1113         << " of the size of the bounding box." << endl;
1114 }
1115
1116
1117
1118 //// No extrusion on cells with multiple patch faces. There ususally is a reason
1119 //// why combinePatchFaces hasn't succeeded.
1120 //void Foam::autoLayerDriver::handleMultiplePatchFaces
1121 //(
1122 //    const indirectPrimitivePatch& pp,
1123 //    pointField& patchDisp,
1124 //    labelList& patchNLayers,
1125 //    List<extrudeMode>& extrudeStatus
1126 //) const
1127 //{
1128 //    const fvMesh& mesh = meshRefiner_.mesh();
1129 //
1130 //    Info<< nl << "Handling cells with multiple patch faces ..." << nl;
1131 //
1132 //    const labelListList& pointFaces = pp.pointFaces();
1133 //
1134 //    // Cells that should not get an extrusion layer
1135 //    cellSet multiPatchCells(mesh, "multiPatchCells", pp.size());
1136 //
1137 //    // Detect points that use multiple faces on same cell.
1138 //    forAll(pointFaces, patchPointI)
1139 //    {
1140 //        const labelList& pFaces = pointFaces[patchPointI];
1141 //
1142 //        labelHashSet pointCells(pFaces.size());
1143 //
1144 //        forAll(pFaces, i)
1145 //        {
1146 //            label cellI = mesh.faceOwner()[pp.addressing()[pFaces[i]]];
1147 //
1148 //            if (!pointCells.insert(cellI))
1149 //            {
1150 //                // Second or more occurrence of cell so cell has two or more
1151 //                // pp faces connected to this point.
1152 //                multiPatchCells.insert(cellI);
1153 //            }
1154 //        }
1155 //    }
1156 //
1157 //    label nMultiPatchCells = returnReduce
1158 //    (
1159 //        multiPatchCells.size(),
1160 //        sumOp<label>()
1161 //    );
1162 //
1163 //    Info<< "Detected " << nMultiPatchCells
1164 //        << " cells with multiple (connected) patch faces." << endl;
1165 //
1166 //    label nChanged = 0;
1167 //
1168 //    if (nMultiPatchCells > 0)
1169 //    {
1170 //        Info<< "Writing " << nMultiPatchCells
1171 //            << " cells with multiple (connected) patch faces to cellSet "
1172 //            << multiPatchCells.objectPath() << endl;
1173 //        multiPatchCells.write();
1174 //
1175 //
1176 //        // Go through all points and remove extrusion on any cell in
1177 //        // multiPatchCells
1178 //        // (has to be done in separate loop since having one point on
1179 //        // multipatches has to reset extrusion on all points of cell)
1180 //
1181 //        forAll(pointFaces, patchPointI)
1182 //        {
1183 //            if (extrudeStatus[patchPointI] != NOEXTRUDE)
1184 //            {
1185 //                const labelList& pFaces = pointFaces[patchPointI];
1186 //
1187 //                forAll(pFaces, i)
1188 //                {
1189 //                    label cellI =
1190 //                        mesh.faceOwner()[pp.addressing()[pFaces[i]]];
1191 //
1192 //                    if (multiPatchCells.found(cellI))
1193 //                    {
1194 //                        if
1195 //                        (
1196 //                            unmarkExtrusion
1197 //                            (
1198 //                                patchPointI,
1199 //                                patchDisp,
1200 //                                patchNLayers,
1201 //                                extrudeStatus
1202 //                            )
1203 //                        )
1204 //                        {
1205 //                            nChanged++;
1206 //                        }
1207 //                    }
1208 //                }
1209 //            }
1210 //        }
1211 //
1212 //        reduce(nChanged, sumOp<label>());
1213 //    }
1214 //
1215 //    Info<< "Prevented extrusion on " << nChanged
1216 //        << " points due to multiple patch faces." << nl << endl;
1217 //}
1218
1219
1220 // No extrusion on faces with differing number of layers for points
1221 void Foam::autoLayerDriver::setNumLayers
1222 (
1223     const labelList& patchToNLayers,
1224     const labelList& patchIDs,
1225     const indirectPrimitivePatch& pp,
1226     pointField& patchDisp,
1227     labelList& patchNLayers,
1228     List<extrudeMode>& extrudeStatus
1229 ) const
1230 {
1231     const fvMesh& mesh = meshRefiner_.mesh();
1232
1233     Info<< nl << "Handling points with inconsistent layer specification ..."
1234         << endl;
1235
1236     // Get for every point (really only nessecary on patch external points)
1237     // the max and min of any patch faces using it.
1238     labelList maxLayers(patchNLayers.size(), labelMin);
1239     labelList minLayers(patchNLayers.size(), labelMax);
1240
1241     forAll(patchIDs, i)
1242     {
1243         label patchI = patchIDs[i];
1244
1245         const labelList& meshPoints = mesh.boundaryMesh()[patchI].meshPoints();
1246
1247         label wantedLayers = patchToNLayers[patchI];
1248
1249         forAll(meshPoints, patchPointI)
1250         {
1251             label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
1252
1253             maxLayers[ppPointI] = max(wantedLayers, maxLayers[ppPointI]);
1254             minLayers[ppPointI] = min(wantedLayers, minLayers[ppPointI]);
1255         }
1256     }
1257
1258     syncTools::syncPointList
1259     (
1260         mesh,
1261         pp.meshPoints(),
1262         maxLayers,
1263         maxEqOp<label>(),
1264         labelMin,           // null value
1265         false               // no separation
1266     );
1267     syncTools::syncPointList
1268     (
1269         mesh,
1270         pp.meshPoints(),
1271         minLayers,
1272         minEqOp<label>(),
1273         labelMax,           // null value
1274         false               // no separation
1275     );
1276
1277     // Unmark any point with different min and max
1278     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1279
1280     //label nConflicts = 0;
1281
1282     forAll(maxLayers, i)
1283     {
1284         if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
1285         {
1286             FatalErrorIn("setNumLayers(..)")
1287                 << "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
1288                 << " maxLayers:" << maxLayers
1289                 << " minLayers:" << minLayers
1290                 << abort(FatalError);
1291         }
1292         else if (maxLayers[i] == minLayers[i])
1293         {
1294             // Ok setting.
1295             patchNLayers[i] = maxLayers[i];
1296         }
1297         else
1298         {
1299             // Inconsistent num layers between patch faces using point
1300             //if
1301             //(
1302             //    unmarkExtrusion
1303             //    (
1304             //        i,
1305             //        patchDisp,
1306             //        patchNLayers,
1307             //        extrudeStatus
1308             //    )
1309             //)
1310             //{
1311             //    nConflicts++;
1312             //}
1313             patchNLayers[i] = maxLayers[i];
1314         }
1315     }
1316
1317     //reduce(nConflicts, sumOp<label>());
1318     //
1319     //Info<< "Set displacement to zero for " << nConflicts
1320     //    << " points due to points being on multiple regions"
1321     //    << " with inconsistent nLayers specification." << endl;
1322 }
1323
1324
1325 // Grow no-extrusion layer
1326 void Foam::autoLayerDriver::growNoExtrusion
1327 (
1328     const indirectPrimitivePatch& pp,
1329     pointField& patchDisp,
1330     labelList& patchNLayers,
1331     List<extrudeMode>& extrudeStatus
1332 )
1333 {
1334     Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
1335
1336     List<extrudeMode> grownExtrudeStatus(extrudeStatus);
1337
1338     const faceList& localFaces = pp.localFaces();
1339
1340     label nGrown = 0;
1341
1342     forAll(localFaces, faceI)
1343     {
1344         const face& f = localFaces[faceI];
1345
1346         bool hasSqueeze = false;
1347         forAll(f, fp)
1348         {
1349             if (extrudeStatus[f[fp]] == NOEXTRUDE)
1350             {
1351                 hasSqueeze = true;
1352                 break;
1353             }
1354         }
1355
1356         if (hasSqueeze)
1357         {
1358             // Squeeze all points of face
1359             forAll(f, fp)
1360             {
1361                 if
1362                 (
1363                     extrudeStatus[f[fp]] == NOEXTRUDE
1364                  && grownExtrudeStatus[f[fp]] != NOEXTRUDE
1365                 )
1366                 {
1367                     grownExtrudeStatus[f[fp]] = NOEXTRUDE;
1368                     nGrown++;
1369                 }
1370             }
1371         }
1372     }
1373
1374     extrudeStatus.transfer(grownExtrudeStatus);
1375
1376     forAll(extrudeStatus, patchPointI)
1377     {
1378         if (extrudeStatus[patchPointI] == NOEXTRUDE)
1379         {
1380             patchDisp[patchPointI] = vector::zero;
1381             patchNLayers[patchPointI] = 0;
1382         }
1383     }
1384
1385     reduce(nGrown, sumOp<label>());
1386
1387     Info<< "Set displacement to zero for an additional " << nGrown
1388         << " points." << endl;
1389 }
1390
1391
1392 void Foam::autoLayerDriver::calculateLayerThickness
1393 (
1394     const indirectPrimitivePatch& pp,
1395     const labelList& patchIDs,
1396     const scalarField& patchExpansionRatio,
1397     const scalarField& patchFinalLayerRatio,
1398     const scalarField& patchRelMinThickness,
1399     const labelList& cellLevel,
1400     const labelList& patchNLayers,
1401     const scalar edge0Len,
1402
1403     scalarField& thickness,
1404     scalarField& minThickness,
1405     scalarField& expansionRatio
1406 ) const
1407 {
1408     const fvMesh& mesh = meshRefiner_.mesh();
1409     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1410
1411     if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
1412     {
1413         FatalErrorIn("calculateLayerThickness(..)")
1414             << "Thickness should be factor of local undistorted cell size."
1415             << " Valid values are [0..2]." << nl
1416             << " minThickness:" << patchRelMinThickness
1417             << exit(FatalError);
1418     }
1419
1420
1421     // Per point the max cell level of connected cells
1422     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1423
1424     labelList maxPointLevel(pp.nPoints(), labelMin);
1425
1426     {
1427         forAll(pp, i)
1428         {
1429             label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
1430
1431             const face& f = pp.localFaces()[i];
1432
1433             forAll(f, fp)
1434             {
1435                 maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
1436             }
1437         }
1438
1439         syncTools::syncPointList
1440         (
1441             mesh,
1442             pp.meshPoints(),
1443             maxPointLevel,
1444             maxEqOp<label>(),
1445             labelMin,           // null value
1446             false               // no separation
1447         );
1448     }
1449
1450
1451     // Rework patch-wise layer parameters into minimum per point
1452     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1453
1454     expansionRatio.setSize(pp.nPoints());
1455     expansionRatio = GREAT;
1456     scalarField finalLayerRatio(pp.nPoints(), GREAT);
1457     scalarField relMinThickness(pp.nPoints(), GREAT);
1458
1459     {
1460         forAll(patchIDs, i)
1461         {
1462             label patchI = patchIDs[i];
1463
1464             const labelList& meshPoints = patches[patchI].meshPoints();
1465
1466             forAll(meshPoints, patchPointI)
1467             {
1468                 label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
1469
1470                 expansionRatio[ppPointI] = min
1471                 (
1472                     expansionRatio[ppPointI],
1473                     patchExpansionRatio[patchI]
1474                 );
1475                 finalLayerRatio[ppPointI] = min
1476                 (
1477                     finalLayerRatio[ppPointI],
1478                     patchFinalLayerRatio[patchI]
1479                 );
1480                 relMinThickness[ppPointI] = min
1481                 (
1482                     relMinThickness[ppPointI],
1483                     patchRelMinThickness[patchI]
1484                 );
1485             }
1486         }
1487
1488         syncTools::syncPointList
1489         (
1490             mesh,
1491             pp.meshPoints(),
1492             expansionRatio,
1493             minEqOp<scalar>(),
1494             GREAT,              // null value
1495             false               // no separation
1496         );
1497         syncTools::syncPointList
1498         (
1499             mesh,
1500             pp.meshPoints(),
1501             finalLayerRatio,
1502             minEqOp<scalar>(),
1503             GREAT,              // null value
1504             false               // no separation
1505         );
1506         syncTools::syncPointList
1507         (
1508             mesh,
1509             pp.meshPoints(),
1510             relMinThickness,
1511             minEqOp<scalar>(),
1512             GREAT,              // null value
1513             false               // no separation
1514         );
1515     }
1516
1517
1518
1519     // Per mesh point the expansion parameters
1520     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1521
1522     thickness.setSize(pp.nPoints());
1523     minThickness.setSize(pp.nPoints());
1524
1525     forAll(maxPointLevel, pointI)
1526     {
1527         // Find undistorted edge size for this level.
1528         scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
1529
1530         // Calculate layer thickness based on expansion ratio
1531         // and final layer height
1532         if (expansionRatio[pointI] == 1)
1533         {
1534             thickness[pointI] =
1535                 finalLayerRatio[pointI]
1536               * patchNLayers[pointI]
1537               * edgeLen;
1538             minThickness[pointI] = relMinThickness[pointI]*edgeLen;
1539         }
1540         else
1541         {
1542             scalar invExpansion = 1.0 / expansionRatio[pointI];
1543             label nLay = patchNLayers[pointI];
1544             thickness[pointI] =
1545                 finalLayerRatio[pointI]
1546               * edgeLen
1547               * (1.0 - pow(invExpansion, nLay))
1548               / (1.0 - invExpansion);
1549             minThickness[pointI] = relMinThickness[pointI]*edgeLen;
1550         }
1551     }
1552
1553     Info<< "calculateLayerThickness : min:" << gMin(thickness)
1554         << " max:" << gMax(thickness) << endl;
1555 }
1556
1557
1558 // Synchronize displacement among coupled patches.
1559 void Foam::autoLayerDriver::syncPatchDisplacement
1560 (
1561     const motionSmoother& meshMover,
1562     const scalarField& minThickness,
1563     pointField& patchDisp,
1564     labelList& patchNLayers,
1565     List<extrudeMode>& extrudeStatus
1566 ) const
1567 {
1568     const fvMesh& mesh = meshRefiner_.mesh();
1569     const labelList& meshPoints = meshMover.patch().meshPoints();
1570
1571     label nChangedTotal = 0;
1572
1573     while (true)
1574     {
1575         label nChanged = 0;
1576
1577         // Sync displacement (by taking min)
1578         syncTools::syncPointList
1579         (
1580             mesh,
1581             meshPoints,
1582             patchDisp,
1583             minEqOp<vector>(),
1584             wallPoint::greatPoint,          // null value
1585             false                           // no separation
1586         );
1587
1588         // Unmark if displacement too small
1589         forAll(patchDisp, i)
1590         {
1591             if (mag(patchDisp[i]) < minThickness[i])
1592             {
1593                 if
1594                 (
1595                     unmarkExtrusion
1596                     (
1597                         i,
1598                         patchDisp,
1599                         patchNLayers,
1600                         extrudeStatus
1601                     )
1602                 )
1603                 {
1604                     nChanged++;
1605                 }
1606             }
1607         }
1608
1609         labelList syncPatchNLayers(patchNLayers);
1610
1611         syncTools::syncPointList
1612         (
1613             mesh,
1614             meshPoints,
1615             syncPatchNLayers,
1616             minEqOp<label>(),
1617             labelMax,           // null value
1618             false               // no separation
1619         );
1620
1621         // Reset if differs
1622         forAll(syncPatchNLayers, i)
1623         {
1624             if (syncPatchNLayers[i] != patchNLayers[i])
1625             {
1626                 if
1627                 (
1628                     unmarkExtrusion
1629                     (
1630                         i,
1631                         patchDisp,
1632                         patchNLayers,
1633                         extrudeStatus
1634                     )
1635                 )
1636                 {
1637                     nChanged++;
1638                 }
1639             }
1640         }
1641
1642         syncTools::syncPointList
1643         (
1644             mesh,
1645             meshPoints,
1646             syncPatchNLayers,
1647             maxEqOp<label>(),
1648             labelMin,           // null value
1649             false               // no separation
1650         );
1651
1652         // Reset if differs
1653         forAll(syncPatchNLayers, i)
1654         {
1655             if (syncPatchNLayers[i] != patchNLayers[i])
1656             {
1657                 if
1658                 (
1659                     unmarkExtrusion
1660                     (
1661                         i,
1662                         patchDisp,
1663                         patchNLayers,
1664                         extrudeStatus
1665                     )
1666                 )
1667                 {
1668                     nChanged++;
1669                 }
1670             }
1671         }
1672         nChangedTotal += nChanged;
1673
1674         if (!returnReduce(nChanged, sumOp<label>()))
1675         {
1676             break;
1677         }
1678     }
1679
1680     Info<< "Prevented extrusion on "
1681         << returnReduce(nChangedTotal, sumOp<label>())
1682         << " coupled patch points during syncPatchDisplacement." << endl;
1683 }
1684
1685
1686 // Calculate displacement vector for all patch points. Uses pointNormal.
1687 // Checks that displaced patch point would be visible from all centres
1688 // of the faces using it.
1689 // extrudeStatus is both input and output and gives the status of each
1690 // patch point.
1691 void Foam::autoLayerDriver::getPatchDisplacement
1692 (
1693     const motionSmoother& meshMover,
1694     const scalarField& thickness,
1695     const scalarField& minThickness,
1696     pointField& patchDisp,
1697     labelList& patchNLayers,
1698     List<extrudeMode>& extrudeStatus
1699 ) const
1700 {
1701     Info<< nl << "Determining displacement for added points"
1702         << " according to pointNormal ..." << endl;
1703
1704     const fvMesh& mesh = meshRefiner_.mesh();
1705     const indirectPrimitivePatch& pp = meshMover.patch();
1706     const vectorField& faceNormals = pp.faceNormals();
1707     const labelListList& pointFaces = pp.pointFaces();
1708     const pointField& localPoints = pp.localPoints();
1709     const labelList& meshPoints = pp.meshPoints();
1710
1711     // Determine pointNormal
1712     // ~~~~~~~~~~~~~~~~~~~~~
1713
1714     pointField pointNormals(pp.nPoints(), vector::zero);
1715     {
1716         labelList nPointFaces(pp.nPoints(), 0);
1717
1718         forAll(faceNormals, faceI)
1719         {
1720             const face& f = pp.localFaces()[faceI];
1721
1722             forAll(f, fp)
1723             {
1724                 pointNormals[f[fp]] += faceNormals[faceI];
1725                 nPointFaces[f[fp]] ++;
1726             }
1727         }
1728
1729         syncTools::syncPointList
1730         (
1731             mesh,
1732             meshPoints,
1733             pointNormals,
1734             plusEqOp<vector>(),
1735             vector::zero,       // null value
1736             false               // no separation
1737         );
1738
1739         syncTools::syncPointList
1740         (
1741             mesh,
1742             meshPoints,
1743             nPointFaces,
1744             plusEqOp<label>(),
1745             0,                  // null value
1746             false               // no separation
1747         );
1748
1749         forAll(pointNormals, i)
1750         {
1751             pointNormals[i] /= nPointFaces[i];
1752         }
1753     }
1754
1755
1756     // Determine local length scale on patch
1757     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1758
1759     // Start off from same thickness everywhere (except where no extrusion)
1760     patchDisp = thickness*pointNormals;
1761
1762     // Check if no extrude possible.
1763     forAll(pointNormals, patchPointI)
1764     {
1765         label meshPointI = pp.meshPoints()[patchPointI];
1766
1767         if (extrudeStatus[patchPointI] == NOEXTRUDE)
1768         {
1769             // Do not use unmarkExtrusion; forcibly set to zero extrusion.
1770             patchNLayers[patchPointI] = 0;
1771             patchDisp[patchPointI] = vector::zero;
1772         }
1773         else
1774         {
1775             // Get normal
1776             const vector& n = pointNormals[patchPointI];
1777
1778             if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
1779             {
1780                 Pout<< "No valid normal for point " << meshPointI
1781                     << ' ' << pp.points()[meshPointI]
1782                     << "; setting displacement to " << patchDisp[patchPointI]
1783                     << endl;
1784
1785                 extrudeStatus[patchPointI] = EXTRUDEREMOVE;
1786             }
1787         }
1788     }
1789
1790     // At illegal points make displacement average of new neighbour positions
1791     forAll(extrudeStatus, patchPointI)
1792     {
1793         if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
1794         {
1795             point avg(vector::zero);
1796             label nPoints = 0;
1797
1798             const labelList& pEdges = pp.pointEdges()[patchPointI];
1799
1800             forAll(pEdges, i)
1801             {
1802                 label edgeI = pEdges[i];
1803
1804                 label otherPointI = pp.edges()[edgeI].otherVertex(patchPointI);
1805
1806                 if (extrudeStatus[otherPointI] != NOEXTRUDE)
1807                 {
1808                     avg += localPoints[otherPointI] + patchDisp[otherPointI];
1809                     nPoints++;
1810                 }
1811             }
1812
1813             if (nPoints > 0)
1814             {
1815                 Pout<< "Displacement at illegal point "
1816                     << localPoints[patchPointI]
1817                     << " set to " << (avg / nPoints - localPoints[patchPointI])
1818                     << endl;
1819
1820                 patchDisp[patchPointI] =
1821                     avg / nPoints
1822                   - localPoints[patchPointI];
1823             }
1824         }
1825     }
1826
1827     // Make sure displacement is equal on both sides of coupled patches.
1828     syncPatchDisplacement
1829     (
1830         meshMover,
1831         minThickness,
1832         patchDisp,
1833         patchNLayers,
1834         extrudeStatus
1835     );
1836
1837     Info<< endl;
1838 }
1839
1840
1841 // Truncates displacement
1842 // - for all patchFaces in the faceset displacement gets set to zero
1843 // - all displacement < minThickness gets set to zero
1844 Foam::label Foam::autoLayerDriver::truncateDisplacement
1845 (
1846     const motionSmoother& meshMover,
1847     const scalarField& minThickness,
1848     const faceSet& illegalPatchFaces,
1849     pointField& patchDisp,
1850     labelList& patchNLayers,
1851     List<extrudeMode>& extrudeStatus
1852 ) const
1853 {
1854     const polyMesh& mesh = meshMover.mesh();
1855     const indirectPrimitivePatch& pp = meshMover.patch();
1856
1857     label nChanged = 0;
1858
1859     const Map<label>& meshPointMap = pp.meshPointMap();
1860
1861     forAllConstIter(faceSet, illegalPatchFaces, iter)
1862     {
1863         label faceI = iter.key();
1864
1865         if (mesh.isInternalFace(faceI))
1866         {
1867             FatalErrorIn("truncateDisplacement(..)")
1868                 << "Faceset " << illegalPatchFaces.name()
1869                 << " contains internal face " << faceI << nl
1870                 << "It should only contain patch faces" << abort(FatalError);
1871         }
1872
1873         const face& f = mesh.faces()[faceI];
1874
1875
1876         forAll(f, fp)
1877         {
1878             if (meshPointMap.found(f[fp]))
1879             {
1880                 label patchPointI = meshPointMap[f[fp]];
1881
1882                 if (extrudeStatus[patchPointI] != NOEXTRUDE)
1883                 {
1884                     unmarkExtrusion
1885                     (
1886                         patchPointI,
1887                         patchDisp,
1888                         patchNLayers,
1889                         extrudeStatus
1890                     );
1891                     nChanged++;
1892                 }
1893             }
1894         }
1895     }
1896
1897     forAll(patchDisp, patchPointI)
1898     {
1899         if (mag(patchDisp[patchPointI]) < minThickness[patchPointI])
1900         {
1901             if
1902             (
1903                 unmarkExtrusion
1904                 (
1905                     patchPointI,
1906                     patchDisp,
1907                     patchNLayers,
1908                     extrudeStatus
1909                 )
1910             )
1911             {
1912                 nChanged++;
1913             }
1914         }
1915         else if (extrudeStatus[patchPointI] == NOEXTRUDE)
1916         {
1917             // Make sure displacement is 0. Should already be so but ...
1918             patchDisp[patchPointI] = vector::zero;
1919             patchNLayers[patchPointI] = 0;
1920         }
1921     }
1922
1923
1924     const faceList& localFaces = pp.localFaces();
1925
1926     while (true)
1927     {
1928         syncPatchDisplacement
1929         (
1930             meshMover,
1931             minThickness,
1932             patchDisp,
1933             patchNLayers,
1934             extrudeStatus
1935         );
1936
1937         // Make sure that a face doesn't have two non-consecutive areas
1938         // not extruded (e.g. quad where vertex 0 and 2 are not extruded
1939         // but 1 and 3 are) since this gives topological errors.
1940
1941         label nPinched = 0;
1942
1943         forAll(localFaces, i)
1944         {
1945             const face& localF = localFaces[i];
1946
1947             // Count number of transitions from unsnapped to snapped.
1948             label nTrans = 0;
1949
1950             extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
1951
1952             forAll(localF, fp)
1953             {
1954                 extrudeMode fpMode = extrudeStatus[localF[fp]];
1955
1956                 if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
1957                 {
1958                     nTrans++;
1959                 }
1960                 prevMode = fpMode;
1961             }
1962
1963             if (nTrans > 1)
1964             {
1965                 // Multiple pinches. Reset whole face as unextruded.
1966                 if
1967                 (
1968                     unmarkExtrusion
1969                     (
1970                         localF,
1971                         patchDisp,
1972                         patchNLayers,
1973                         extrudeStatus
1974                     )
1975                 )
1976                 {
1977                     nPinched++;
1978                     nChanged++;
1979                 }
1980             }
1981         }
1982
1983         reduce(nPinched, sumOp<label>());
1984
1985         Info<< "truncateDisplacement : Unextruded " << nPinched
1986             << " faces due to non-consecutive vertices being extruded." << endl;
1987
1988
1989         // Make sure that a face has consistent number of layers for all
1990         // its vertices.
1991
1992         label nDiffering = 0;
1993
1994         //forAll(localFaces, i)
1995         //{
1996         //    const face& localF = localFaces[i];
1997         //
1998         //    label numLayers = -1;
1999         //
2000         //    forAll(localF, fp)
2001         //    {
2002         //        if (patchNLayers[localF[fp]] > 0)
2003         //        {
2004         //            if (numLayers == -1)
2005         //            {
2006         //                numLayers = patchNLayers[localF[fp]];
2007         //            }
2008         //            else if (numLayers != patchNLayers[localF[fp]])
2009         //            {
2010         //                // Differing number of layers
2011         //                if
2012         //                (
2013         //                    unmarkExtrusion
2014         //                    (
2015         //                        localF,
2016         //                        patchDisp,
2017         //                        patchNLayers,
2018         //                        extrudeStatus
2019         //                    )
2020         //                )
2021         //                {
2022         //                    nDiffering++;
2023         //                    nChanged++;
2024         //                }
2025         //                break;
2026         //            }
2027         //        }
2028         //    }
2029         //}
2030         //
2031         //reduce(nDiffering, sumOp<label>());
2032         //
2033         //Info<< "truncateDisplacement : Unextruded " << nDiffering
2034         //    << " faces due to having differing number of layers." << endl;
2035
2036         if (nPinched+nDiffering == 0)
2037         {
2038             break;
2039         }
2040     }
2041
2042     return nChanged;
2043 }
2044
2045
2046 // Setup layer information (at points and faces) to modify mesh topology in
2047 // regions where layer mesh terminates.
2048 void Foam::autoLayerDriver::setupLayerInfoTruncation
2049 (
2050     const motionSmoother& meshMover,
2051     const labelList& patchNLayers,
2052     const List<extrudeMode>& extrudeStatus,
2053     const label nBufferCellsNoExtrude,
2054     labelList& nPatchPointLayers,
2055     labelList& nPatchFaceLayers
2056 ) const
2057 {
2058     Info<< nl << "Setting up information for layer truncation ..." << endl;
2059
2060     const indirectPrimitivePatch& pp = meshMover.patch();
2061     const polyMesh& mesh = meshMover.mesh();
2062
2063     if (nBufferCellsNoExtrude < 0)
2064     {
2065         Info<< nl << "Performing no layer truncation."
2066             << " nBufferCellsNoExtrude set to less than 0  ..." << endl;
2067
2068         // Face layers if any point get extruded
2069         forAll(pp.localFaces(), patchFaceI)
2070         {
2071             const face& f = pp.localFaces()[patchFaceI];
2072
2073             forAll(f, fp)
2074             {
2075                 if (patchNLayers[f[fp]] > 0)
2076                 {
2077                     nPatchFaceLayers[patchFaceI] = patchNLayers[f[fp]];
2078                     break;
2079                 }
2080             }
2081         }
2082         nPatchPointLayers = patchNLayers;
2083     }
2084     else
2085     {
2086         // Determine max point layers per face.
2087         labelList maxLevel(pp.size(), 0);
2088
2089         forAll(pp.localFaces(), patchFaceI)
2090         {
2091             const face& f = pp.localFaces()[patchFaceI];
2092
2093             // find patch faces where layer terminates (i.e contains extrude
2094             // and noextrude points).
2095
2096             bool noExtrude = false;
2097             label mLevel = 0;
2098
2099             forAll(f, fp)
2100             {
2101                 if (extrudeStatus[f[fp]] == NOEXTRUDE)
2102                 {
2103                     noExtrude = true;
2104                 }
2105                 mLevel = max(mLevel, patchNLayers[f[fp]]);
2106             }
2107
2108             if (mLevel > 0)
2109             {
2110                 // So one of the points is extruded. Check if all are extruded
2111                 // or is a mix.
2112
2113                 if (noExtrude)
2114                 {
2115                     nPatchFaceLayers[patchFaceI] = 1;
2116                     maxLevel[patchFaceI] = mLevel;
2117                 }
2118                 else
2119                 {
2120                     maxLevel[patchFaceI] = mLevel;
2121                 }
2122             }
2123         }
2124
2125         // We have the seed faces (faces with nPatchFaceLayers != maxLevel)
2126         // Now do a meshwave across the patch where we pick up neighbours
2127         // of seed faces.
2128         // Note: quite inefficient. Could probably be coded better.
2129
2130         const labelListList& pointFaces = pp.pointFaces();
2131
2132         label nLevels = gMax(patchNLayers);
2133
2134         // flag neighbouring patch faces with number of layers to grow
2135         for (label ilevel = 1; ilevel < nLevels; ilevel++)
2136         {
2137             label nBuffer;
2138
2139             if (ilevel == 1)
2140             {
2141                 nBuffer = nBufferCellsNoExtrude - 1;
2142             }
2143             else
2144             {
2145                 nBuffer = nBufferCellsNoExtrude;
2146             }
2147
2148             for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
2149             {
2150                 labelList tempCounter(nPatchFaceLayers);
2151
2152                 boolList foundNeighbour(pp.nPoints(), false);
2153
2154                 forAll(pp.meshPoints(), patchPointI)
2155                 {
2156                     forAll(pointFaces[patchPointI], pointFaceI)
2157                     {
2158                         label faceI = pointFaces[patchPointI][pointFaceI];
2159
2160                         if
2161                         (
2162                             nPatchFaceLayers[faceI] != -1
2163                          && maxLevel[faceI] > 0
2164                         )
2165                         {
2166                             foundNeighbour[patchPointI] = true;
2167                             break;
2168                         }
2169                     }
2170                 }
2171
2172                 syncTools::syncPointList
2173                 (
2174                     mesh,
2175                     pp.meshPoints(),
2176                     foundNeighbour,
2177                     orEqOp<bool>(),
2178                     false,              // null value
2179                     false               // no separation
2180                 );
2181
2182                 forAll(pp.meshPoints(), patchPointI)
2183                 {
2184                     if (foundNeighbour[patchPointI])
2185                     {
2186                         forAll(pointFaces[patchPointI], pointFaceI)
2187                         {
2188                             label faceI = pointFaces[patchPointI][pointFaceI];
2189                             if
2190                             (
2191                                 nPatchFaceLayers[faceI] == -1
2192                              && maxLevel[faceI] > 0
2193                              && ilevel < maxLevel[faceI]
2194                             )
2195                             {
2196                                 tempCounter[faceI] = ilevel;
2197                             }
2198                         }
2199                     }
2200                 }
2201                 nPatchFaceLayers = tempCounter;
2202             }
2203         }
2204
2205         forAll(pp.localFaces(), patchFaceI)
2206         {
2207             if (nPatchFaceLayers[patchFaceI] == -1)
2208             {
2209                 nPatchFaceLayers[patchFaceI] = maxLevel[patchFaceI];
2210             }
2211         }
2212
2213         forAll(pp.meshPoints(), patchPointI)
2214         {
2215             if (extrudeStatus[patchPointI] != NOEXTRUDE)
2216             {
2217                 forAll(pointFaces[patchPointI], pointFaceI)
2218                 {
2219                     label face = pointFaces[patchPointI][pointFaceI];
2220                     nPatchPointLayers[patchPointI] = max
2221                     (
2222                         nPatchPointLayers[patchPointI],
2223                         nPatchFaceLayers[face]
2224                     );
2225                 }
2226             }
2227             else
2228             {
2229                 nPatchPointLayers[patchPointI] = 0;
2230             }
2231         }
2232         syncTools::syncPointList
2233         (
2234             mesh,
2235             pp.meshPoints(),
2236             nPatchPointLayers,
2237             maxEqOp<label>(),
2238             0,                  // null value
2239             false               // no separation
2240         );
2241     }
2242 }
2243
2244
2245 // Does any of the cells use a face from faces?
2246 bool Foam::autoLayerDriver::cellsUseFace
2247 (
2248     const polyMesh& mesh,
2249     const labelList& cellLabels,
2250     const labelHashSet& faces
2251 )
2252 {
2253     forAll(cellLabels, i)
2254     {
2255         const cell& cFaces = mesh.cells()[cellLabels[i]];
2256
2257         forAll(cFaces, cFaceI)
2258         {
2259             if (faces.found(cFaces[cFaceI]))
2260             {
2261                 return true;
2262             }
2263         }
2264     }
2265     return false;
2266 }
2267
2268
2269 // Checks the newly added cells and locally unmarks points so they
2270 // will not get extruded next time round. Returns global number of unmarked
2271 // points (0 if all was fine)
2272 Foam::label Foam::autoLayerDriver::checkAndUnmark
2273 (
2274     const addPatchCellLayer& addLayer,
2275     const dictionary& motionDict,
2276     const indirectPrimitivePatch& pp,
2277     const fvMesh& newMesh,
2278
2279     pointField& patchDisp,
2280     labelList& patchNLayers,
2281     List<extrudeMode>& extrudeStatus
2282 )
2283 {
2284     // Check the resulting mesh for errors
2285     Info<< nl << "Checking mesh with layer ..." << endl;
2286     faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
2287     motionSmoother::checkMesh(false, newMesh, motionDict, wrongFaces);
2288     Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
2289         << " illegal faces"
2290         << " (concave, zero area or negative cell pyramid volume)"
2291         << endl;
2292
2293     // Undo local extrusion if
2294     // - any of the added cells in error
2295
2296     label nChanged = 0;
2297
2298     // Get all cells in the layer.
2299     labelListList addedCells
2300     (
2301         addPatchCellLayer::addedCells
2302         (
2303             newMesh,
2304             addLayer.layerFaces()
2305         )
2306     );
2307
2308     // Check if any of the faces in error uses any face of an added cell
2309     forAll(addedCells, oldPatchFaceI)
2310     {
2311         // Get the cells (in newMesh labels) per old patch face (in mesh
2312         // labels)
2313         const labelList& fCells = addedCells[oldPatchFaceI];
2314
2315         if (cellsUseFace(newMesh, fCells, wrongFaces))
2316         {
2317             // Unmark points on old mesh
2318             if
2319             (
2320                 unmarkExtrusion
2321                 (
2322                     pp.localFaces()[oldPatchFaceI],
2323                     patchDisp,
2324                     patchNLayers,
2325                     extrudeStatus
2326                 )
2327             )
2328             {
2329                 nChanged++;
2330             }
2331         }
2332     }
2333
2334     return returnReduce(nChanged, sumOp<label>());
2335 }
2336
2337
2338 //- Count global number of extruded faces
2339 Foam::label Foam::autoLayerDriver::countExtrusion
2340 (
2341     const indirectPrimitivePatch& pp,
2342     const List<extrudeMode>& extrudeStatus
2343 )
2344 {
2345     // Count number of extruded patch faces
2346     label nExtruded = 0;
2347     {
2348         const faceList& localFaces = pp.localFaces();
2349
2350         forAll(localFaces, i)
2351         {
2352             const face& localFace = localFaces[i];
2353
2354             forAll(localFace, fp)
2355             {
2356                 if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
2357                 {
2358                     nExtruded++;
2359                     break;
2360                 }
2361             }
2362         }
2363     }
2364
2365     return returnReduce(nExtruded, sumOp<label>());
2366 }
2367
2368
2369 // Collect layer faces and layer cells into bools for ease of handling
2370 void Foam::autoLayerDriver::getLayerCellsFaces
2371 (
2372     const polyMesh& mesh,
2373     const addPatchCellLayer& addLayer,
2374     boolList& flaggedCells,
2375     boolList& flaggedFaces
2376 )
2377 {
2378     flaggedCells.setSize(mesh.nCells());
2379     flaggedCells = false;
2380     flaggedFaces.setSize(mesh.nFaces());
2381     flaggedFaces = false;
2382
2383     // Mark all faces in the layer
2384     const labelListList& layerFaces = addLayer.layerFaces();
2385
2386     // Mark all cells in the layer.
2387     labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
2388
2389     forAll(addedCells, oldPatchFaceI)
2390     {
2391         const labelList& added = addedCells[oldPatchFaceI];
2392
2393         forAll(added, i)
2394         {
2395             flaggedCells[added[i]] = true;
2396         }
2397     }
2398
2399     forAll(layerFaces, oldPatchFaceI)
2400     {
2401         const labelList& layer = layerFaces[oldPatchFaceI];
2402
2403         if (layer.size() > 0)
2404         {
2405             for (label i = 1; i < layer.size()-1; i++)
2406             {
2407                 flaggedFaces[layer[i]] = true;
2408             }
2409         }
2410     }
2411 }
2412
2413
2414 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
2415
2416 Foam::autoLayerDriver::autoLayerDriver
2417 (
2418     meshRefinement& meshRefiner,
2419     const labelList& globalToPatch
2420 )
2421 :
2422     meshRefiner_(meshRefiner),
2423     globalToPatch_(globalToPatch)
2424 {}
2425
2426
2427 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
2428
2429 void Foam::autoLayerDriver::mergePatchFacesUndo
2430 (
2431     const layerParameters& layerParams,
2432     const dictionary& motionDict
2433 )
2434 {
2435     scalar minCos = Foam::cos
2436     (
2437         layerParams.featureAngle()
2438       * mathematicalConstant::pi/180.0
2439     );
2440
2441     scalar concaveCos = Foam::cos
2442     (
2443         layerParams.concaveAngle()
2444       * mathematicalConstant::pi/180.0
2445     );
2446
2447     Info<< nl
2448         << "Merging all faces of a cell" << nl
2449         << "---------------------------" << nl
2450         << "    - which are on the same patch" << nl
2451         << "    - which make an angle < " << layerParams.featureAngle()
2452         << " degrees"
2453         << nl
2454         << "      (cos:" << minCos << ')' << nl
2455         << "    - as long as the resulting face doesn't become concave"
2456         << " by more than "
2457         << layerParams.concaveAngle()
2458         << " degrees (0=straight, 180=fully concave)" << nl
2459         << endl;
2460
2461     label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
2462
2463     nChanged += mergeEdgesUndo(minCos, motionDict);
2464 }
2465
2466
2467 void Foam::autoLayerDriver::addLayers
2468 (
2469     const layerParameters& layerParams,
2470     const dictionary& motionDict,
2471     const label nAllowableErrors,
2472     motionSmoother& meshMover,
2473     decompositionMethod& decomposer,
2474     fvMeshDistribute& distributor
2475 )
2476 {
2477     fvMesh& mesh = meshRefiner_.mesh();
2478     const indirectPrimitivePatch& pp = meshMover.patch();
2479     const labelList& meshPoints = pp.meshPoints();
2480
2481     // Precalculate mesh edge labels for patch edges
2482     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2483
2484     labelList meshEdges(pp.nEdges());
2485     forAll(pp.edges(), edgeI)
2486     {
2487         const edge& ppEdge = pp.edges()[edgeI];
2488         label v0 = meshPoints[ppEdge[0]];
2489         label v1 = meshPoints[ppEdge[1]];
2490         meshEdges[edgeI] = meshTools::findEdge
2491         (
2492             mesh.edges(),
2493             mesh.pointEdges()[v0],
2494             v0,
2495             v1
2496         );
2497     }
2498
2499     // Displacement for all pp.localPoints.
2500     vectorField patchDisp(pp.nPoints(), vector::one);
2501
2502     // Number of layers for all pp.localPoints. Note: only valid if
2503     // extrudeStatus = EXTRUDE.
2504     labelList patchNLayers(pp.nPoints(), 0);
2505
2506     // Whether to add edge for all pp.localPoints.
2507     List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
2508
2509
2510     // Get number of layer per point from number of layers per patch
2511     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2512
2513     setNumLayers
2514     (
2515         layerParams.numLayers(),    // per patch the num layers
2516         meshMover.adaptPatchIDs(),  // patches that are being moved
2517         pp,                         // indirectpatch for all faces moving
2518
2519         patchDisp,
2520         patchNLayers,
2521         extrudeStatus
2522     );
2523
2524     // Disable extrusion on non-manifold points
2525     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2526
2527     handleNonManifolds
2528     (
2529         pp,
2530         meshEdges,
2531
2532         patchDisp,
2533         patchNLayers,
2534         extrudeStatus
2535     );
2536
2537     // Disable extrusion on feature angles
2538     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2539
2540     handleFeatureAngle
2541     (
2542         pp,
2543         meshEdges,
2544         layerParams.featureAngle()*mathematicalConstant::pi/180.0,
2545
2546         patchDisp,
2547         patchNLayers,
2548         extrudeStatus
2549     );
2550
2551     // Disable extrusion on warped faces
2552     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2553
2554     // Undistorted edge length
2555     const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
2556     const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
2557
2558     handleWarpedFaces
2559     (
2560         pp,
2561         layerParams.maxFaceThicknessRatio(),
2562         edge0Len,
2563         cellLevel,
2564
2565         patchDisp,
2566         patchNLayers,
2567         extrudeStatus
2568     );
2569
2570     //// Disable extrusion on cells with multiple patch faces
2571     //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2572     //
2573     //handleMultiplePatchFaces
2574     //(
2575     //    pp,
2576     //
2577     //    patchDisp,
2578     //    patchNLayers,
2579     //    extrudeStatus
2580     //);
2581
2582
2583     // Grow out region of non-extrusion
2584     for (label i = 0; i < layerParams.nGrow(); i++)
2585     {
2586         growNoExtrusion
2587         (
2588             pp,
2589             patchDisp,
2590             patchNLayers,
2591             extrudeStatus
2592         );
2593     }
2594
2595     // Determine (wanted) point-wise layer thickness and expansion ratio
2596     scalarField thickness(pp.nPoints());
2597     scalarField minThickness(pp.nPoints());
2598     scalarField expansionRatio(pp.nPoints());
2599     calculateLayerThickness
2600     (
2601         pp,
2602         meshMover.adaptPatchIDs(),
2603         layerParams.expansionRatio(),
2604         layerParams.finalLayerRatio(),
2605         layerParams.minThickness(),
2606         cellLevel,
2607         patchNLayers,
2608         edge0Len,
2609
2610         thickness,
2611         minThickness,
2612         expansionRatio
2613     );
2614
2615     // Calculate wall to medial axis distance for smoothing displacement
2616     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2617
2618     pointScalarField pointMedialDist
2619     (
2620         IOobject
2621         (
2622             "pointMedialDist",
2623             mesh.time().timeName(),
2624             mesh,
2625             IOobject::NO_READ,
2626             IOobject::NO_WRITE,
2627             false
2628         ),
2629         meshMover.pMesh(),
2630         dimensionedScalar("pointMedialDist", dimless, 0.0)
2631     );
2632
2633     pointVectorField dispVec
2634     (
2635         IOobject
2636         (
2637             "dispVec",
2638             mesh.time().timeName(),
2639             mesh,
2640             IOobject::NO_READ,
2641             IOobject::NO_WRITE,
2642             false
2643         ),
2644         meshMover.pMesh(),
2645         dimensionedVector("dispVec", dimless, vector::zero)
2646     );
2647
2648     pointScalarField medialRatio
2649     (
2650         IOobject
2651         (
2652             "medialRatio",
2653             mesh.time().timeName(),
2654             mesh,
2655             IOobject::NO_READ,
2656             IOobject::NO_WRITE,
2657             false
2658         ),
2659         meshMover.pMesh(),
2660         dimensionedScalar("medialRatio", dimless, 0.0)
2661     );
2662
2663     // Setup information for medial axis smoothing. Calculates medial axis
2664     // and a smoothed displacement direction.
2665     // - pointMedialDist : distance to medial axis
2666     // - dispVec : normalised direction of nearest displacement
2667     // - medialRatio : ratio of medial distance to wall distance.
2668     //   (1 at wall, 0 at medial axis)
2669     medialAxisSmoothingInfo
2670     (
2671         meshMover,
2672         layerParams.nSmoothNormals(),
2673         layerParams.nSmoothSurfaceNormals(),
2674         layerParams.minMedianAxisAngleCos(),
2675
2676         dispVec,
2677         medialRatio,
2678         pointMedialDist
2679     );
2680
2681
2682
2683     // Saved old points
2684     pointField oldPoints(mesh.points());
2685
2686     // Last set of topology changes. (changing mesh clears out polyTopoChange)
2687     polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
2688
2689     boolList flaggedCells;
2690     boolList flaggedFaces;
2691
2692     while (true)
2693     {
2694         // Make sure displacement is equal on both sides of coupled patches.
2695         syncPatchDisplacement
2696         (
2697             meshMover,
2698             minThickness,
2699             patchDisp,
2700             patchNLayers,
2701             extrudeStatus
2702         );
2703
2704         // Displacement acc. to pointnormals
2705         getPatchDisplacement
2706         (
2707             meshMover,
2708             thickness,
2709             minThickness,
2710             patchDisp,
2711             patchNLayers,
2712             extrudeStatus
2713         );
2714
2715         // Shrink mesh by displacement value first.
2716         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2717
2718         {
2719             pointField oldPatchPos(pp.localPoints());
2720
2721             //// Laplacian displacement shrinking.
2722             //shrinkMeshDistance
2723             //(
2724             //    debug,
2725             //    meshMover,
2726             //    -patchDisp,     // Shrink in opposite direction of addedPoints
2727             //    layerParams.nSmoothDisp(),
2728             //    layerParams.nSnap()
2729             //);
2730
2731             // Medial axis based shrinking
2732             shrinkMeshMedialDistance
2733             (
2734                 meshMover,
2735
2736                 layerParams.nSmoothThickness(),
2737                 layerParams.maxThicknessToMedialRatio(),
2738                 nAllowableErrors,
2739                 layerParams.nSnap(),
2740                 layerParams.layerTerminationCos(),
2741
2742                 thickness,
2743                 minThickness,
2744
2745                 dispVec,
2746                 medialRatio,
2747                 pointMedialDist,
2748
2749                 extrudeStatus,
2750                 patchDisp,
2751                 patchNLayers
2752             );
2753
2754             // Update patchDisp (since not all might have been honoured)
2755             patchDisp = oldPatchPos - pp.localPoints();
2756         }
2757
2758         // Truncate displacements that are too small (this will do internal
2759         // ones, coupled ones have already been truncated by syncPatch)
2760         faceSet dummySet(mesh, "wrongPatchFaces", 0);
2761         truncateDisplacement
2762         (
2763             meshMover,
2764             minThickness,
2765             dummySet,
2766             patchDisp,
2767             patchNLayers,
2768             extrudeStatus
2769         );
2770
2771
2772         // Dump to .obj file for debugging.
2773         if (debug)
2774         {
2775             dumpDisplacement
2776             (
2777                 mesh.time().path()/"layer",
2778                 pp,
2779                 patchDisp,
2780                 extrudeStatus
2781             );
2782
2783             const_cast<Time&>(mesh.time())++;
2784             Info<< "Writing shrunk mesh to " << mesh.time().timeName() << endl;
2785             mesh.write();
2786         }
2787
2788
2789         // Mesh topo change engine
2790         polyTopoChange meshMod(mesh);
2791
2792         // Grow layer of cells on to patch. Handles zero sized displacement.
2793         addPatchCellLayer addLayer(mesh);
2794
2795         // Determine per point/per face number of layers to extrude. Also
2796         // handles the slow termination of layers when going switching layers
2797
2798         labelList nPatchPointLayers(pp.nPoints(),-1);
2799         labelList nPatchFaceLayers(pp.localFaces().size(),-1);
2800         setupLayerInfoTruncation
2801         (
2802             meshMover,
2803             patchNLayers,
2804             extrudeStatus,
2805             layerParams.nBufferCellsNoExtrude(),
2806             nPatchPointLayers,
2807             nPatchFaceLayers
2808         );
2809
2810         // Calculate displacement for first layer for addPatchLayer
2811         vectorField firstDisp(patchNLayers.size(), vector::zero);
2812
2813         forAll(patchNLayers, i)
2814         {
2815             if (patchNLayers[i] > 0)
2816             {
2817                 if (expansionRatio[i] == 1.0)
2818                 {
2819                     firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
2820                 }
2821                 else
2822                 {
2823                     label nLay = nPatchPointLayers[i];
2824                     scalar h =
2825                         pow(expansionRatio[i], nLay - 1)
2826                       * (mag(patchDisp[i])*(1.0 - expansionRatio[i]))
2827                       / (1.0 - pow(expansionRatio[i], nLay));
2828                     firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
2829                 }
2830             }
2831         }
2832
2833         scalarField invExpansionRatio = 1.0 / expansionRatio;
2834
2835         // Add topo regardless of whether extrudeStatus is extruderemove.
2836         // Not add layer if patchDisp is zero.
2837         addLayer.setRefinement
2838         (
2839             invExpansionRatio,
2840             pp,
2841             nPatchFaceLayers,   // layers per face
2842             nPatchPointLayers,  // layers per point
2843             firstDisp,          // thickness of first layer
2844             meshMod
2845         );
2846
2847         if (debug)
2848         {
2849             const_cast<Time&>(mesh.time())++;
2850         }
2851
2852         // Store mesh changes for if mesh is correct.
2853         savedMeshMod = meshMod;
2854
2855
2856         // With the stored topo changes we create a new mesh so we can
2857         // undo if neccesary.
2858
2859         autoPtr<fvMesh> newMeshPtr;
2860         autoPtr<mapPolyMesh> map = meshMod.makeMesh
2861         (
2862             newMeshPtr,
2863             IOobject
2864             (
2865                 //mesh.name()+"_layer",
2866                 mesh.name(),
2867                 static_cast<polyMesh&>(mesh).instance(),
2868                 mesh.db(),
2869                 static_cast<polyMesh&>(mesh).readOpt(),
2870                 static_cast<polyMesh&>(mesh).writeOpt()
2871             ),          // io params from original mesh but new name
2872             mesh,      // original mesh
2873             true        // parallel sync
2874         );
2875         fvMesh& newMesh = newMeshPtr();
2876
2877         //?neccesary? Update fields
2878         newMesh.updateMesh(map);
2879
2880         // Update numbering on addLayer:
2881         // - cell/point labels to be newMesh.
2882         // - patchFaces to remain in oldMesh order.
2883         addLayer.updateMesh
2884         (
2885             map,
2886             identity(pp.size()),
2887             identity(pp.nPoints())
2888         );
2889
2890         // Collect layer faces and cells for outside loop.
2891         getLayerCellsFaces
2892         (
2893             newMesh,
2894             addLayer,
2895             flaggedCells,
2896             flaggedFaces
2897         );
2898
2899
2900         if (debug)
2901         {
2902             Info<< "Writing layer mesh to " << mesh.time().timeName() << endl;
2903             newMesh.write();
2904             cellSet addedCellSet
2905             (
2906                 newMesh,
2907                 "addedCells",
2908                 findIndices(flaggedCells, true)
2909             );
2910             Info<< "Writing "
2911                 << returnReduce(addedCellSet.size(), sumOp<label>())
2912                 << " added cells to cellSet "
2913                 << addedCellSet.name() << endl;
2914             addedCellSet.write();
2915
2916             faceSet layerFacesSet
2917             (
2918                 newMesh,
2919                 "layerFaces",
2920                 findIndices(flaggedCells, true)
2921             );
2922             Info<< "Writing "
2923                 << returnReduce(layerFacesSet.size(), sumOp<label>())
2924                 << " faces inside added layer to faceSet "
2925                 << layerFacesSet.name() << endl;
2926             layerFacesSet.write();
2927         }
2928
2929         // Unset the extrusion at the pp.
2930         label nTotChanged = checkAndUnmark
2931         (
2932             addLayer,
2933             motionDict,
2934             pp,
2935             newMesh,
2936
2937             patchDisp,
2938             patchNLayers,
2939             extrudeStatus
2940         );
2941
2942         Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
2943             << " out of " << returnReduce(pp.size(), sumOp<label>())
2944             << " faces. Removed extrusion at " << nTotChanged << " faces."
2945             << endl;
2946
2947         if (nTotChanged == 0)
2948         {
2949             break;
2950         }
2951
2952         // Reset mesh points and start again
2953         meshMover.movePoints(oldPoints);
2954         meshMover.correct();
2955     }
2956
2957
2958     // At this point we have a (shrunk) mesh and a set of topology changes
2959     // which will make a valid mesh with layer. Apply these changes to the
2960     // current mesh.
2961
2962     // Apply the stored topo changes to the current mesh.
2963     autoPtr<mapPolyMesh> map = savedMeshMod.changeMesh(mesh, false);
2964
2965     // Update fields
2966     mesh.updateMesh(map);
2967
2968     // Move mesh (since morphing does not do this)
2969     if (map().hasMotionPoints())
2970     {
2971         mesh.movePoints(map().preMotionPoints());
2972     }
2973
2974     meshRefiner_.updateMesh(map, labelList(0));
2975
2976
2977     // Do final balancing
2978     // ~~~~~~~~~~~~~~~~~~
2979
2980     if (Pstream::parRun())
2981     {
2982         Info<< nl
2983             << "Doing final balancing" << nl
2984             << "---------------------" << nl
2985             << endl;
2986
2987         if (debug)
2988         {
2989             const_cast<Time&>(mesh.time())++;
2990         }
2991
2992         // Balance. No restriction on face zones and baffles.
2993         autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
2994         (
2995             false,
2996             false,
2997             decomposer,
2998             distributor
2999         );
3000
3001         // Re-distribute flag of layer faces and cells
3002         map().distributeCellData(flaggedCells);
3003         map().distributeFaceData(flaggedFaces);
3004     }
3005
3006
3007     // Write mesh
3008     // ~~~~~~~~~~
3009
3010     cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
3011     Info<< "Writing "
3012         << returnReduce(addedCellSet.size(), sumOp<label>())
3013         << " added cells to cellSet "
3014         << addedCellSet.name() << endl;
3015     addedCellSet.write();
3016
3017     faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
3018     Info<< "Writing "
3019         << returnReduce(layerFacesSet.size(), sumOp<label>())
3020         << " faces inside added layer to faceSet "
3021         << layerFacesSet.name() << endl;
3022     layerFacesSet.write();
3023 }
3024
3025
3026 void Foam::autoLayerDriver::doLayers
3027 (
3028     const dictionary& shrinkDict,
3029     const dictionary& motionDict,
3030     const layerParameters& layerParams,
3031     decompositionMethod& decomposer,
3032     fvMeshDistribute& distributor
3033 )
3034 {
3035     fvMesh& mesh = meshRefiner_.mesh();
3036
3037     Info<< nl
3038         << "Shrinking and layer addition phase" << nl
3039         << "----------------------------------" << nl
3040         << endl;
3041
3042     const_cast<Time&>(mesh.time())++;
3043
3044     Info<< "Using mesh parameters " << motionDict << nl << endl;
3045
3046     // Merge coplanar boundary faces
3047     mergePatchFacesUndo(layerParams, motionDict);
3048
3049     // Per patch the number of layers (0 if no layer)
3050     const labelList& numLayers = layerParams.numLayers();
3051
3052     // Patches that need to get a layer
3053     DynamicList<label> patchIDs(numLayers.size());
3054     label nFacesWithLayers = 0;
3055     forAll(numLayers, patchI)
3056     {
3057         if (numLayers[patchI] > 0)
3058         {
3059             patchIDs.append(patchI);
3060             nFacesWithLayers += mesh.boundaryMesh()[patchI].size();
3061         }
3062     }
3063     patchIDs.shrink();
3064
3065     if (returnReduce(nFacesWithLayers, sumOp<label>()) == 0)
3066     {
3067         Info<< nl << "No layers to generate ..." << endl;
3068     }
3069     else
3070     {
3071         autoPtr<indirectPrimitivePatch> ppPtr
3072         (
3073             meshRefinement::makePatch
3074             (
3075                 mesh,
3076                 patchIDs
3077             )
3078         );
3079         indirectPrimitivePatch& pp = ppPtr();
3080
3081         // Construct iterative mesh mover.
3082         Info<< "Constructing mesh displacer ..." << endl;
3083
3084         {
3085             pointMesh pMesh(mesh);
3086
3087             motionSmoother meshMover
3088             (
3089                 mesh,
3090                 pp,
3091                 patchIDs,
3092                 meshRefinement::makeDisplacementField(pMesh, patchIDs),
3093                 motionDict
3094             );
3095
3096             // Check that outside of mesh is not multiply connected.
3097             checkMeshManifold();
3098
3099             // Check initial mesh
3100             Info<< "Checking initial mesh ..." << endl;
3101             labelHashSet wrongFaces(mesh.nFaces()/100);
3102             motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
3103             const label nInitErrors = returnReduce
3104             (
3105                 wrongFaces.size(),
3106                 sumOp<label>()
3107             );
3108
3109             Info<< "Detected " << nInitErrors << " illegal faces"
3110                 << " (concave, zero area or negative cell pyramid volume)"
3111                 << endl;
3112
3113             // Do all topo changes
3114             addLayers
3115             (
3116                 layerParams,
3117                 motionDict,
3118                 nInitErrors,
3119                 meshMover,
3120                 decomposer,
3121                 distributor
3122             );
3123         }
3124     }
3125 }
3126
3127
3128 // ************************************************************************* //