src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Shrinking mesh (part of adding cell layers)
  26
  27 \*----------------------------------------------------------------------------*/
  28
  29 #include "autoLayerDriver.H"
  30 #include "fvMesh.H"
  31 #include "Time.H"
  32 #include "pointFields.H"
  33 #include "motionSmoother.H"
  34 #include "pointData.H"
  35 #include "PointEdgeWave.H"
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 // Calculate inverse sum of edge weights (currently always 1.0)
  40 void Foam::autoLayerDriver::sumWeights
  41 (
  42     const PackedBoolList& isMasterEdge,
  43     const labelList& meshEdges,
  44     const labelList& meshPoints,
  45     const edgeList& edges,
  46     scalarField& invSumWeight
  47 ) const
  48 {
  49     invSumWeight = 0;
  50
  51     forAll(edges, edgeI)
  52     {
  53         if (isMasterEdge.get(meshEdges[edgeI]) == 1)
  54         {
  55             const edge& e = edges[edgeI];
  56             //scalar eWeight = edgeWeights[edgeI];
  57             scalar eWeight = 1.0;
  58
  59             invSumWeight[e[0]] += eWeight;
  60             invSumWeight[e[1]] += eWeight;
  61         }
  62     }
  63
  64     syncTools::syncPointList
  65     (
  66         meshRefiner_.mesh(),
  67         meshPoints,
  68         invSumWeight,
  69         plusEqOp<scalar>(),
  70         scalar(0.0)         // null value
  71     );
  72
  73     forAll(invSumWeight, pointI)
  74     {
  75         scalar w = invSumWeight[pointI];
  76
  77         if (w > 0.0)
  78         {
  79             invSumWeight[pointI] = 1.0/w;
  80         }
  81     }
  82 }
  83
  84
  85 // Smooth field on moving patch
  86 void Foam::autoLayerDriver::smoothField
  87 (
  88     const motionSmoother& meshMover,
  89     const PackedBoolList& isMasterEdge,
  90     const labelList& meshEdges,
  91     const scalarField& fieldMin,
  92     const label nSmoothDisp,
  93     scalarField& field
  94 ) const
  95 {
  96     const indirectPrimitivePatch& pp = meshMover.patch();
  97     const edgeList& edges = pp.edges();
  98     const labelList& meshPoints = pp.meshPoints();
  99
 100     scalarField invSumWeight(pp.nPoints());
 101     sumWeights
 102     (
 103         isMasterEdge,
 104         meshEdges,
 105         meshPoints,
 106         edges,
 107         invSumWeight
 108     );
 109
 110     // Get smoothly varying patch field.
 111     Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
 112
 113     for (label iter = 0; iter < nSmoothDisp; iter++)
 114     {
 115         scalarField average(pp.nPoints());
 116         averageNeighbours
 117         (
 118             meshMover.mesh(),
 119             isMasterEdge,
 120             meshEdges,
 121             meshPoints,
 122             pp.edges(),
 123             invSumWeight,
 124             field,
 125             average
 126         );
 127
 128         // Transfer to field
 129         forAll(field, pointI)
 130         {
 131             //full smoothing neighbours + point value
 132             average[pointI] = 0.5*(field[pointI]+average[pointI]);
 133
 134             // perform monotonic smoothing
 135             if
 136             (
 137                 average[pointI] < field[pointI]
 138              && average[pointI] >= fieldMin[pointI]
 139             )
 140             {
 141                 field[pointI] = average[pointI];
 142             }
 143         }
 144
 145         // Do residual calculation every so often.
 146         if ((iter % 10) == 0)
 147         {
 148             Info<< "    Iteration " << iter << "   residual "
 149                 <<  gSum(mag(field-average))
 150                    /returnReduce(average.size(), sumOp<label>())
 151                 << endl;
 152         }
 153     }
 154 }
 155
 156
 157 // Smooth normals on moving patch.
 158 void Foam::autoLayerDriver::smoothPatchNormals
 159 (
 160     const motionSmoother& meshMover,
 161     const PackedBoolList& isMasterEdge,
 162     const labelList& meshEdges,
 163     const label nSmoothDisp,
 164     pointField& normals
 165 ) const
 166 {
 167     const indirectPrimitivePatch& pp = meshMover.patch();
 168     const edgeList& edges = pp.edges();
 169     const labelList& meshPoints = pp.meshPoints();
 170
 171     // Calculate inverse sum of weights
 172
 173     scalarField invSumWeight(pp.nPoints());
 174     sumWeights
 175     (
 176         isMasterEdge,
 177         meshEdges,
 178         meshPoints,
 179         edges,
 180         invSumWeight
 181     );
 182
 183     // Get smoothly varying internal normals field.
 184     Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
 185
 186     for (label iter = 0; iter < nSmoothDisp; iter++)
 187     {
 188         vectorField average(pp.nPoints());
 189         averageNeighbours
 190         (
 191             meshMover.mesh(),
 192             isMasterEdge,
 193             meshEdges,
 194             meshPoints,
 195             pp.edges(),
 196             invSumWeight,
 197             normals,
 198             average
 199         );
 200
 201         // Do residual calculation every so often.
 202         if ((iter % 10) == 0)
 203         {
 204             Info<< "    Iteration " << iter << "   residual "
 205                 <<  gSum(mag(normals-average))
 206                    /returnReduce(average.size(), sumOp<label>())
 207                 << endl;
 208         }
 209
 210         // Transfer to normals vector field
 211         forAll(average, pointI)
 212         {
 213             // full smoothing neighbours + point value
 214             average[pointI] = 0.5*(normals[pointI]+average[pointI]);
 215             normals[pointI] = average[pointI];
 216             normals[pointI] /= mag(normals[pointI]) + VSMALL;
 217         }
 218     }
 219 }
 220
 221
 222 // Smooth normals in interior.
 223 void Foam::autoLayerDriver::smoothNormals
 224 (
 225     const label nSmoothDisp,
 226     const PackedBoolList& isMasterEdge,
 227     const labelList& fixedPoints,
 228     pointVectorField& normals
 229 ) const
 230 {
 231     // Get smoothly varying internal normals field.
 232     Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
 233
 234     const fvMesh& mesh = meshRefiner_.mesh();
 235     const edgeList& edges = mesh.edges();
 236
 237     // Points that do not change.
 238     PackedBoolList isFixedPoint(mesh.nPoints());
 239
 240     // Internal points that are fixed
 241     forAll(fixedPoints, i)
 242     {
 243         label meshPointI = fixedPoints[i];
 244         isFixedPoint.set(meshPointI, 1);
 245     }
 246
 247     // Correspondence between local edges/points and mesh edges/points
 248     const labelList meshEdges(identity(mesh.nEdges()));
 249     const labelList meshPoints(identity(mesh.nPoints()));
 250
 251     // Calculate inverse sum of weights
 252
 253     scalarField invSumWeight(mesh.nPoints(), 0);
 254     sumWeights
 255     (
 256         isMasterEdge,
 257         meshEdges,
 258         meshPoints,
 259         edges,
 260         invSumWeight
 261     );
 262
 263     Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
 264
 265     for (label iter = 0; iter < nSmoothDisp; iter++)
 266     {
 267         vectorField average(mesh.nPoints());
 268         averageNeighbours
 269         (
 270             mesh,
 271             isMasterEdge,
 272             meshEdges,
 273             meshPoints,
 274             edges,
 275             invSumWeight,
 276             normals,
 277             average
 278         );
 279
 280         // Do residual calculation every so often.
 281         if ((iter % 10) == 0)
 282         {
 283             Info<< "    Iteration " << iter << "   residual "
 284                 <<  gSum(mag(normals-average))
 285                    /returnReduce(average.size(), sumOp<label>())
 286                 << endl;
 287         }
 288
 289
 290         // Transfer to normals vector field
 291         forAll(average, pointI)
 292         {
 293             if (isFixedPoint.get(pointI) == 0)
 294             {
 295                 //full smoothing neighbours + point value
 296                 average[pointI] = 0.5*(normals[pointI]+average[pointI]);
 297                 normals[pointI] = average[pointI];
 298                 normals[pointI] /= mag(normals[pointI]) + VSMALL;
 299             }
 300         }
 301     }
 302 }
 303
 304
 305 // Tries and find a medial axis point. Done by comparing vectors to nearest
 306 // wall point for both vertices of edge.
 307 bool Foam::autoLayerDriver::isMaxEdge
 308 (
 309     const List<pointData>& pointWallDist,
 310     const label edgeI,
 311     const scalar minCos
 312 ) const
 313 {
 314     const fvMesh& mesh = meshRefiner_.mesh();
 315     const pointField& points = mesh.points();
 316
 317     // Do not mark edges with one side on moving wall.
 318
 319     const edge& e = mesh.edges()[edgeI];
 320
 321     vector v0(points[e[0]] - pointWallDist[e[0]].origin());
 322     scalar magV0(mag(v0));
 323
 324     if (magV0 < SMALL)
 325     {
 326         return false;
 327     }
 328
 329     vector v1(points[e[1]] - pointWallDist[e[1]].origin());
 330     scalar magV1(mag(v1));
 331
 332     if (magV1 < SMALL)
 333     {
 334         return false;
 335     }
 336
 337     v0 /= magV0;
 338     v1 /= magV1;
 339
 340     // Test angle.
 341     if ((v0 & v1) < minCos)
 342     {
 343         return true;
 344     }
 345     else
 346     {
 347         return false;
 348     }
 349 }
 350
 351
 352 // Stop layer growth where mesh wraps around edge with a
 353 // large feature angle
 354 void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
 355 (
 356     const indirectPrimitivePatch& pp,
 357     const scalar minCos,
 358     List<extrudeMode>& extrudeStatus,
 359     pointField& patchDisp,
 360     labelList& patchNLayers,
 361     label& nPointCounter
 362 ) const
 363 {
 364     // Mark faces that have all points extruded
 365     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 366
 367     boolList extrudedFaces(pp.size(), true);
 368
 369     forAll(pp.localFaces(), faceI)
 370     {
 371         const face& f = pp.localFaces()[faceI];
 372
 373         forAll(f, fp)
 374         {
 375             if (extrudeStatus[f[fp]] == NOEXTRUDE)
 376             {
 377                 extrudedFaces[faceI] = false;
 378                 break;
 379             }
 380         }
 381     }
 382
 383
 384     // Detect situation where two featureedge-neighbouring faces are partly or
 385     // not extruded and the edge itself is extruded. In this case unmark the
 386     // edge for extrusion.
 387
 388     forAll(pp.edgeFaces(), edgeI)
 389     {
 390         const labelList& eFaces = pp.edgeFaces()[edgeI];
 391
 392         if (eFaces.size() == 2)
 393         {
 394             const edge& e = pp.edges()[edgeI];
 395             label v0 = e[0];
 396             label v1 = e[1];
 397
 398             if
 399             (
 400                 extrudeStatus[v0] != NOEXTRUDE
 401              || extrudeStatus[v1] != NOEXTRUDE
 402             )
 403             {
 404                 if (!extrudedFaces[eFaces[0]] || !extrudedFaces[eFaces[1]])
 405                 {
 406                     const vector& n0 = pp.faceNormals()[eFaces[0]];
 407                     const vector& n1 = pp.faceNormals()[eFaces[1]];
 408
 409                     if ((n0 & n1) < minCos)
 410                     {
 411                         if
 412                         (
 413                             unmarkExtrusion
 414                             (
 415                                 v0,
 416                                 patchDisp,
 417                                 patchNLayers,
 418                                 extrudeStatus
 419                             )
 420                         )
 421                         {
 422                             nPointCounter++;
 423                         }
 424                         if
 425                         (
 426                             unmarkExtrusion
 427                             (
 428                                 v1,
 429                                 patchDisp,
 430                                 patchNLayers,
 431                                 extrudeStatus
 432                             )
 433                         )
 434                         {
 435                             nPointCounter++;
 436                         }
 437                     }
 438                 }
 439             }
 440         }
 441     }
 442 }
 443
 444
 445 // Find isolated islands (points, edges and faces and layer terminations)
 446 // in the layer mesh and stop any layer growth at these points.
 447 void Foam::autoLayerDriver::findIsolatedRegions
 448 (
 449     const indirectPrimitivePatch& pp,
 450     const PackedBoolList& isMasterEdge,
 451     const labelList& meshEdges,
 452     const scalar minCosLayerTermination,
 453     scalarField& field,
 454     List<extrudeMode>& extrudeStatus,
 455     pointField& patchDisp,
 456     labelList& patchNLayers
 457 ) const
 458 {
 459     const fvMesh& mesh = meshRefiner_.mesh();
 460
 461     Info<< "shrinkMeshDistance : Removing isolated regions ..." << endl;
 462
 463     // Keep count of number of points unextruded
 464     label nPointCounter = 0;
 465
 466     while (true)
 467     {
 468         // Stop layer growth where mesh wraps around edge with a
 469         // large feature angle
 470         handleFeatureAngleLayerTerminations
 471         (
 472             pp,
 473             minCosLayerTermination,
 474
 475             extrudeStatus,
 476             patchDisp,
 477             patchNLayers,
 478             nPointCounter
 479         );
 480
 481
 482         // Do not extrude from point where all neighbouring
 483         // faces are not grown
 484         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 485         boolList extrudedFaces(pp.size(), true);
 486         forAll(pp.localFaces(), faceI)
 487         {
 488             const face& f = pp.localFaces()[faceI];
 489             forAll(f, fp)
 490             {
 491                 if (extrudeStatus[f[fp]] == NOEXTRUDE)
 492                 {
 493                     extrudedFaces[faceI] = false;
 494                     break;
 495                 }
 496             }
 497         }
 498
 499         const labelListList& pointFaces = pp.pointFaces();
 500
 501         boolList keptPoints(pp.nPoints(), false);
 502         forAll(keptPoints, patchPointI)
 503         {
 504             const labelList& pFaces = pointFaces[patchPointI];
 505
 506             forAll(pFaces, i)
 507             {
 508                 label faceI = pFaces[i];
 509                 if (extrudedFaces[faceI])
 510                 {
 511                     keptPoints[patchPointI] = true;
 512                     break;
 513                 }
 514             }
 515         }
 516
 517         syncTools::syncPointList
 518         (
 519             mesh,
 520             pp.meshPoints(),
 521             keptPoints,
 522             orEqOp<bool>(),
 523             false               // null value
 524         );
 525
 526         label nChanged = 0;
 527
 528         forAll(keptPoints, patchPointI)
 529         {
 530             if (!keptPoints[patchPointI])
 531             {
 532                 if
 533                 (
 534                     unmarkExtrusion
 535                     (
 536                         patchPointI,
 537                         patchDisp,
 538                         patchNLayers,
 539                         extrudeStatus
 540                     )
 541                 )
 542                 {
 543                    nPointCounter++;
 544                    nChanged++;
 545                 }
 546            }
 547        }
 548
 549
 550        if (returnReduce(nChanged, sumOp<label>()) == 0)
 551        {
 552            break;
 553        }
 554    }
 555
 556     const edgeList& edges = pp.edges();
 557
 558
 559     // Count number of mesh edges using a point
 560     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 561
 562     labelList isolatedPoint(pp.nPoints(),0);
 563
 564     forAll(edges, edgeI)
 565     {
 566         if (isMasterEdge.get(meshEdges[edgeI]) == 1)
 567         {
 568             const edge& e = edges[edgeI];
 569
 570             label v0 = e[0];
 571             label v1 = e[1];
 572
 573             if (extrudeStatus[v1] != NOEXTRUDE)
 574             {
 575                 isolatedPoint[v0] += 1;
 576             }
 577             if (extrudeStatus[v0] != NOEXTRUDE)
 578             {
 579                 isolatedPoint[v1] += 1;
 580             }
 581         }
 582     }
 583
 584     syncTools::syncPointList
 585     (
 586         mesh,
 587         pp.meshPoints(),
 588         isolatedPoint,
 589         plusEqOp<label>(),
 590         0        // null value
 591     );
 592
 593     // stop layer growth on isolated faces
 594     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 595     forAll(pp, faceI)
 596     {
 597         const face& f = pp.localFaces()[faceI];
 598         bool failed = false;
 599         forAll(f, fp)
 600         {
 601             if (isolatedPoint[f[fp]] > 2)
 602             {
 603                 failed = true;
 604                 break;
 605             }
 606         }
 607         bool allPointsExtruded = true;
 608         if (!failed)
 609         {
 610             forAll(f, fp)
 611             {
 612                 if (extrudeStatus[f[fp]] == NOEXTRUDE)
 613                 {
 614                     allPointsExtruded = false;
 615                     break;
 616                 }
 617             }
 618
 619             if (allPointsExtruded)
 620             {
 621                 forAll(f, fp)
 622                 {
 623                     if
 624                     (
 625                         unmarkExtrusion
 626                         (
 627                             f[fp],
 628                             patchDisp,
 629                             patchNLayers,
 630                             extrudeStatus
 631                         )
 632                     )
 633                     {
 634                         field[f[fp]] = 0.0;
 635                         nPointCounter++;
 636                     }
 637                 }
 638             }
 639         }
 640     }
 641
 642     reduce(nPointCounter, sumOp<label>());
 643     Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
 644 }
 645
 646
 647 // Calculates medial axis fields:
 648 // dispVec     : normal of nearest wall. Where this normal changes direction
 649 //               defines the medial axis
 650 // medialDist  : distance to medial axis
 651 // medialRatio : ratio of medial distance to wall distance.
 652 //               (1 at wall, 0 at medial axis)
 653 void Foam::autoLayerDriver::medialAxisSmoothingInfo
 654 (
 655     const motionSmoother& meshMover,
 656     const label nSmoothNormals,
 657     const label nSmoothSurfaceNormals,
 658     const scalar minMedianAxisAngleCos,
 659
 660     pointVectorField& dispVec,
 661     pointScalarField& medialRatio,
 662     pointScalarField& medialDist
 663 ) const
 664 {
 665
 666     Info<< "medialAxisSmoothingInfo :"
 667         << " Calculate distance to Medial Axis ..." << endl;
 668
 669     const polyMesh& mesh = meshMover.mesh();
 670     const pointField& points = mesh.points();
 671
 672     const indirectPrimitivePatch& pp = meshMover.patch();
 673     const vectorField& faceNormals = pp.faceNormals();
 674     const labelList& meshPoints = pp.meshPoints();
 675
 676     // Predetermine mesh edges
 677     // ~~~~~~~~~~~~~~~~~~~~~~~
 678
 679     // Precalulate master edge (only relevant for shared edges)
 680     PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
 681     // Precalculate meshEdge per pp edge
 682     labelList meshEdges(pp.nEdges());
 683
 684     forAll(meshEdges, patchEdgeI)
 685     {
 686         const edge& e = pp.edges()[patchEdgeI];
 687
 688         label v0 = pp.meshPoints()[e[0]];
 689         label v1 = pp.meshPoints()[e[1]];
 690         meshEdges[patchEdgeI] = meshTools::findEdge
 691         (
 692             mesh.edges(),
 693             mesh.pointEdges()[v0],
 694             v0,
 695             v1
 696         );
 697     }
 698
 699
 700     // Determine pointNormal
 701     // ~~~~~~~~~~~~~~~~~~~~~
 702
 703     pointField pointNormals(pp.nPoints(), vector::zero);
 704     {
 705         labelList nPointFaces(pp.nPoints(), 0);
 706
 707         forAll(faceNormals, faceI)
 708         {
 709             const face& f = pp.localFaces()[faceI];
 710
 711             forAll(f, fp)
 712             {
 713                 pointNormals[f[fp]] += faceNormals[faceI];
 714                 nPointFaces[f[fp]] ++;
 715             }
 716         }
 717
 718         syncTools::syncPointList
 719         (
 720             mesh,
 721             meshPoints,
 722             pointNormals,
 723             plusEqOp<vector>(),
 724             vector::zero        // null value
 725         );
 726
 727         syncTools::syncPointList
 728         (
 729             mesh,
 730             meshPoints,
 731             nPointFaces,
 732             plusEqOp<label>(),
 733             0                   // null value
 734         );
 735
 736         forAll(pointNormals, i)
 737         {
 738             pointNormals[i] /= nPointFaces[i];
 739         }
 740     }
 741
 742     // Smooth patch normal vectors
 743     smoothPatchNormals
 744     (
 745         meshMover,
 746         isMasterEdge,
 747         meshEdges,
 748         nSmoothSurfaceNormals,
 749         pointNormals
 750     );
 751
 752
 753     // Calculate distance to pp points
 754     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 755
 756     // Distance to wall
 757     List<pointData> pointWallDist(mesh.nPoints());
 758
 759     // Dummy additional info for PointEdgeWave
 760     int dummyTrackData = 0;
 761
 762
 763     // 1. Calculate distance to points where displacement is specified.
 764     {
 765         // Seed data.
 766         List<pointData> wallInfo(meshPoints.size());
 767
 768         forAll(meshPoints, patchPointI)
 769         {
 770             label pointI = meshPoints[patchPointI];
 771             wallInfo[patchPointI] = pointData
 772             (
 773                 points[pointI],
 774                 0.0,
 775                 pointI,                       // passive scalar
 776                 pointNormals[patchPointI]     // surface normals
 777             );
 778         }
 779
 780         // Do all calculations
 781         List<pointData> edgeWallDist(mesh.nEdges());
 782         PointEdgeWave<pointData> wallDistCalc
 783         (
 784             mesh,
 785             meshPoints,
 786             wallInfo,
 787             pointWallDist,
 788             edgeWallDist,
 789             mesh.globalData().nTotalPoints(),   // max iterations
 790             dummyTrackData
 791         );
 792     }
 793
 794     // 2. Find points with max distance and transport information back to
 795     //    wall.
 796     {
 797         List<pointData> pointMedialDist(mesh.nPoints());
 798         List<pointData> edgeMedialDist(mesh.nEdges());
 799
 800         // Seed point data.
 801         DynamicList<pointData> maxInfo(meshPoints.size());
 802         DynamicList<label> maxPoints(meshPoints.size());
 803
 804         // 1. Medial axis points
 805
 806         const edgeList& edges = mesh.edges();
 807
 808         forAll(edges, edgeI)
 809         {
 810             if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos))
 811             {
 812                 // Both end points of edge have very different nearest wall
 813                 // point. Mark both points as medial axis points.
 814                 const edge& e = edges[edgeI];
 815
 816                 forAll(e, ep)
 817                 {
 818                     label pointI = e[ep];
 819
 820                     if (!pointMedialDist[pointI].valid(dummyTrackData))
 821                     {
 822                         maxPoints.append(pointI);
 823                         maxInfo.append
 824                         (
 825                             pointData
 826                             (
 827                                 points[pointI],
 828                                 0.0,
 829                                 pointI,         // passive data
 830                                 vector::zero    // passive data
 831                             )
 832                         );
 833                         pointMedialDist[pointI] = maxInfo.last();
 834                     }
 835                 }
 836             }
 837         }
 838
 839
 840         // 2. Seed non-adapt patches
 841         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 842
 843         labelHashSet adaptPatches(meshMover.adaptPatchIDs());
 844
 845         forAll(patches, patchI)
 846         {
 847             const polyPatch& pp = patches[patchI];
 848
 849             if
 850             (
 851                 !pp.coupled()
 852              && !isA<emptyPolyPatch>(pp)
 853              && !adaptPatches.found(patchI)
 854             )
 855             {
 856                 Info<< "Inserting points on patch " << pp.name() << endl;
 857
 858                 const labelList& meshPoints = pp.meshPoints();
 859
 860                 forAll(meshPoints, i)
 861                 {
 862                     label pointI = meshPoints[i];
 863
 864                     if (!pointMedialDist[pointI].valid(dummyTrackData))
 865                     {
 866                         maxPoints.append(pointI);
 867                         maxInfo.append
 868                         (
 869                             pointData
 870                             (
 871                                 points[pointI],
 872                                 0.0,
 873                                 pointI,         // passive data
 874                                 vector::zero    // passive data
 875                             )
 876                         );
 877                         pointMedialDist[pointI] = maxInfo.last();
 878                     }
 879                 }
 880             }
 881         }
 882
 883         maxInfo.shrink();
 884         maxPoints.shrink();
 885
 886         // Do all calculations
 887         PointEdgeWave<pointData> medialDistCalc
 888         (
 889             mesh,
 890             maxPoints,
 891             maxInfo,
 892
 893             pointMedialDist,
 894             edgeMedialDist,
 895             mesh.globalData().nTotalPoints(),   // max iterations
 896             dummyTrackData
 897         );
 898
 899         // Extract medial axis distance as pointScalarField
 900         forAll(pointMedialDist, pointI)
 901         {
 902             medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
 903         }
 904     }
 905
 906     // Extract transported surface normals as pointVectorField
 907     forAll(dispVec, i)
 908     {
 909         dispVec[i] = pointWallDist[i].v();
 910     }
 911
 912     // Smooth normal vectors. Do not change normals on pp.meshPoints
 913     smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
 914
 915     // Calculate ratio point medial distance to point wall distance
 916     forAll(medialRatio, pointI)
 917     {
 918         scalar wDist2 = pointWallDist[pointI].distSqr();
 919         scalar mDist = medialDist[pointI];
 920
 921         if (wDist2 < sqr(SMALL) && mDist < SMALL)
 922         {
 923             medialRatio[pointI] = 0.0;
 924         }
 925         else
 926         {
 927             medialRatio[pointI] = mDist / (Foam::sqrt(wDist2) + mDist);
 928         }
 929     }
 930
 931     if (debug)
 932     {
 933         Info<< "medialAxisSmoothingInfo :"
 934             << " Writing:" << nl
 935             << "    " << dispVec.name()
 936             << " : normalised direction of nearest displacement" << nl
 937             << "    " << medialDist.name()
 938             << " : distance to medial axis" << nl
 939             << "    " << medialRatio.name()
 940             << " : ratio of medial distance to wall distance" << nl
 941             << endl;
 942         dispVec.write();
 943         medialDist.write();
 944         medialRatio.write();
 945     }
 946 }
 947
 948
 949 void Foam::autoLayerDriver::shrinkMeshMedialDistance
 950 (
 951     motionSmoother& meshMover,
 952     const dictionary& meshQualityDict,
 953     const List<labelPair>& baffles,
 954     const label nSmoothThickness,
 955     const scalar maxThicknessToMedialRatio,
 956     const label nAllowableErrors,
 957     const label nSnap,
 958     const scalar minCosLayerTermination,
 959
 960     const scalarField& layerThickness,
 961     const scalarField& minThickness,
 962
 963     const pointVectorField& dispVec,
 964     const pointScalarField& medialRatio,
 965     const pointScalarField& medialDist,
 966
 967     List<extrudeMode>& extrudeStatus,
 968     pointField& patchDisp,
 969     labelList& patchNLayers
 970 ) const
 971 {
 972     Info<< "shrinkMeshMedialDistance : Smoothing using Medial Axis ..." << endl;
 973
 974     const polyMesh& mesh = meshMover.mesh();
 975
 976     const indirectPrimitivePatch& pp = meshMover.patch();
 977     const labelList& meshPoints = pp.meshPoints();
 978
 979     // Precalulate master edge (only relevant for shared edges)
 980     PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
 981     // Precalculate meshEdge per pp edge
 982     labelList meshEdges(pp.nEdges());
 983
 984     forAll(meshEdges, patchEdgeI)
 985     {
 986         const edge& e = pp.edges()[patchEdgeI];
 987
 988         label v0 = pp.meshPoints()[e[0]];
 989         label v1 = pp.meshPoints()[e[1]];
 990         meshEdges[patchEdgeI] = meshTools::findEdge
 991         (
 992             mesh.edges(),
 993             mesh.pointEdges()[v0],
 994             v0,
 995             v1
 996         );
 997     }
 998
 999
1000     scalarField thickness(layerThickness.size());
1001
1002     thickness = mag(patchDisp);
1003
1004     forAll(thickness, patchPointI)
1005     {
1006         if (extrudeStatus[patchPointI] == NOEXTRUDE)
1007         {
1008             thickness[patchPointI] = 0.0;
1009         }
1010     }
1011
1012     label numThicknessRatioExclude = 0;
1013
1014     // reduce thickness where thickness/medial axis distance large
1015     forAll(meshPoints, patchPointI)
1016     {
1017         if (extrudeStatus[patchPointI] != NOEXTRUDE)
1018         {
1019             label pointI = meshPoints[patchPointI];
1020
1021             scalar mDist = medialDist[pointI];
1022
1023             scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1024
1025             if (thicknessRatio > maxThicknessToMedialRatio)
1026             {
1027                 // Truncate thickness.
1028                 thickness[patchPointI] =
1029                     0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1030
1031                 patchDisp[patchPointI] =
1032                     thickness[patchPointI]
1033                   * patchDisp[patchPointI]
1034                   / (mag(patchDisp[patchPointI]) + VSMALL);
1035                 numThicknessRatioExclude++;
1036             }
1037         }
1038     }
1039
1040     reduce(numThicknessRatioExclude, sumOp<label>());
1041     Info<< "shrinkMeshMedialDistance : Reduce layer thickness at "
1042         << numThicknessRatioExclude
1043         << " nodes where thickness to medial axis distance is large " << endl;
1044
1045     // find points where layer growth isolated to a lone point, edge or face
1046
1047     findIsolatedRegions
1048     (
1049         pp,
1050         isMasterEdge,
1051         meshEdges,
1052         minCosLayerTermination,
1053         thickness,
1054         extrudeStatus,
1055         patchDisp,
1056         patchNLayers
1057     );
1058
1059     // smooth layer thickness on moving patch
1060     smoothField
1061     (
1062         meshMover,
1063         isMasterEdge,
1064         meshEdges,
1065         minThickness,
1066         nSmoothThickness,
1067         thickness
1068     );
1069
1070     List<pointData> pointWallDist(mesh.nPoints());
1071
1072     const pointField& points = mesh.points();
1073     // 1. Calculate distance to points where displacement is specified.
1074     // This wave is used to transport layer thickness
1075     {
1076         // Distance to wall and medial axis on edges.
1077         List<pointData> edgeWallDist(mesh.nEdges());
1078         labelList wallPoints(meshPoints.size());
1079
1080         // Seed data.
1081         List<pointData> wallInfo(meshPoints.size());
1082
1083         forAll(meshPoints, patchPointI)
1084         {
1085             label pointI = meshPoints[patchPointI];
1086             wallPoints[patchPointI] = pointI;
1087             wallInfo[patchPointI] = pointData
1088             (
1089                 points[pointI],
1090                 0.0,
1091                 thickness[patchPointI],       // transport layer thickness
1092                 vector::zero                  // passive vector
1093             );
1094         }
1095
1096         // Do all calculations
1097         PointEdgeWave<pointData> wallDistCalc
1098         (
1099             mesh,
1100             wallPoints,
1101             wallInfo,
1102             pointWallDist,
1103             edgeWallDist,
1104             mesh.globalData().nTotalPoints()    // max iterations
1105         );
1106     }
1107
1108     // Calculate scaled displacement vector
1109     pointVectorField& displacement = meshMover.displacement();
1110
1111     forAll(displacement, pointI)
1112     {
1113         // 1) displacement on nearest wall point, scaled by medialRatio
1114         //    (wall distance / medial distance)
1115         // 2) pointWallDist[pointI].s() is layer thickness transported
1116         //    from closest wall point.
1117         // 3) shrink in opposite direction of addedPoints
1118         displacement[pointI] =
1119             -medialRatio[pointI]
1120             *pointWallDist[pointI].s()
1121             *dispVec[pointI];
1122     }
1123
1124     // Current faces to check. Gets modified in meshMover.scaleMesh
1125     labelList checkFaces(identity(mesh.nFaces()));
1126
1127     Info<< "shrinkMeshMedialDistance : Moving mesh ..." << endl;
1128     scalar oldErrorReduction = -1;
1129
1130     for (label iter = 0; iter < 2*nSnap ; iter++)
1131     {
1132         Info<< "Iteration " << iter << endl;
1133         if (iter == nSnap)
1134         {
1135             Info<< "Displacement scaling for error reduction set to 0." << endl;
1136             oldErrorReduction = meshMover.setErrorReduction(0.0);
1137         }
1138
1139         if
1140         (
1141             meshMover.scaleMesh
1142             (
1143                 checkFaces,
1144                 baffles,
1145                 meshMover.paramDict(),
1146                 meshQualityDict,
1147                 true,
1148                 nAllowableErrors
1149             )
1150         )
1151         {
1152             Info<< "shrinkMeshMedialDistance : Successfully moved mesh" << endl;
1153             break;
1154         }
1155     }
1156
1157     if (oldErrorReduction >= 0)
1158     {
1159         meshMover.setErrorReduction(oldErrorReduction);
1160     }
1161
1162     Info<< "shrinkMeshMedialDistance : Finished moving mesh ..." << endl;
1163 }
1164
1165
1166 // ************************************************************************* //