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