src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.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 snapping to the surface
  26
  27 \*----------------------------------------------------------------------------*/
  28
  29 #include "autoSnapDriver.H"
  30 #include "motionSmoother.H"
  31 #include "polyTopoChange.H"
  32 #include "OFstream.H"
  33 #include "syncTools.H"
  34 #include "fvMesh.H"
  35 #include "Time.H"
  36 #include "OFstream.H"
  37 #include "mapPolyMesh.H"
  38 #include "pointEdgePoint.H"
  39 #include "PointEdgeWave.H"
  40 #include "mergePoints.H"
  41 #include "snapParameters.H"
  42 #include "refinementSurfaces.H"
  43 #include "unitConversion.H"
  44
  45 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  46
  47 namespace Foam
  48 {
  49
  50 defineTypeNameAndDebug(autoSnapDriver, 0);
  51
  52 } // End namespace Foam
  53
  54
  55 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  56
  57 // Calculate geometrically collocated points, Requires PackedList to be
  58 // sized and initalised!
  59 Foam::label Foam::autoSnapDriver::getCollocatedPoints
  60 (
  61     const scalar tol,
  62     const pointField& points,
  63     PackedBoolList& isCollocatedPoint
  64 )
  65 {
  66     labelList pointMap;
  67     pointField newPoints;
  68     bool hasMerged = mergePoints
  69     (
  70         points,                         // points
  71         tol,                            // mergeTol
  72         false,                          // verbose
  73         pointMap,
  74         newPoints
  75     );
  76
  77     if (!returnReduce(hasMerged, orOp<bool>()))
  78     {
  79         return 0;
  80     }
  81
  82     // Determine which newPoints are referenced more than once
  83     label nCollocated = 0;
  84
  85     // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
  86     // twice)
  87     labelList firstOldPoint(newPoints.size(), -1);
  88     forAll(pointMap, oldPointI)
  89     {
  90         label newPointI = pointMap[oldPointI];
  91
  92         if (firstOldPoint[newPointI] == -1)
  93         {
  94             // First use of oldPointI. Store.
  95             firstOldPoint[newPointI] = oldPointI;
  96         }
  97         else if (firstOldPoint[newPointI] == -2)
  98         {
  99             // Third or more reference of oldPointI -> non-manifold
 100             isCollocatedPoint.set(oldPointI, 1u);
 101             nCollocated++;
 102         }
 103         else
 104         {
 105             // Second reference of oldPointI -> non-manifold
 106             isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
 107             nCollocated++;
 108
 109             isCollocatedPoint.set(oldPointI, 1u);
 110             nCollocated++;
 111
 112             // Mark with special value to save checking next time round
 113             firstOldPoint[newPointI] = -2;
 114         }
 115     }
 116     return returnReduce(nCollocated, sumOp<label>());
 117 }
 118
 119
 120 // Calculate displacement as average of patch points.
 121 Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
 122 (
 123     const motionSmoother& meshMover,
 124     const List<labelPair>& baffles
 125 ) const
 126 {
 127     const indirectPrimitivePatch& pp = meshMover.patch();
 128
 129     // Calculate geometrically non-manifold points on the patch to be moved.
 130     PackedBoolList nonManifoldPoint(pp.nPoints());
 131     label nNonManifoldPoints = getCollocatedPoints
 132     (
 133         SMALL,
 134         pp.localPoints(),
 135         nonManifoldPoint
 136     );
 137     Info<< "Found " << nNonManifoldPoints << " non-mainfold point(s)."
 138         << endl;
 139
 140
 141     // Average points
 142     // ~~~~~~~~~~~~~~
 143
 144     // We determine three points:
 145     // - average of (centres of) connected patch faces
 146     // - average of (centres of) connected internal mesh faces
 147     // - as fallback: centre of any connected cell
 148     // so we can do something moderately sensible for non/manifold points.
 149
 150     // Note: the averages are calculated properly parallel. This is
 151     // necessary to get the points shared by processors correct.
 152
 153
 154     const labelListList& pointFaces = pp.pointFaces();
 155     const labelList& meshPoints = pp.meshPoints();
 156     const pointField& points = pp.points();
 157     const polyMesh& mesh = meshMover.mesh();
 158
 159     // Get labels of faces to count (master of coupled faces and baffle pairs)
 160     PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
 161
 162     {
 163         forAll(baffles, i)
 164         {
 165             label f0 = baffles[i].first();
 166             label f1 = baffles[i].second();
 167
 168             if (isMasterFace.get(f0))
 169             {
 170                 // Make f1 a slave
 171                 isMasterFace.unset(f1);
 172             }
 173             else if (isMasterFace.get(f1))
 174             {
 175                 isMasterFace.unset(f0);
 176             }
 177             else
 178             {
 179                 FatalErrorIn("autoSnapDriver::smoothPatchDisplacement(..)")
 180                     << "Both sides of baffle consisting of faces " << f0
 181                     << " and " << f1 << " are already slave faces."
 182                     << abort(FatalError);
 183             }
 184         }
 185     }
 186
 187
 188     // Get average position of boundary face centres
 189     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 190
 191     vectorField avgBoundary(pointFaces.size(), vector::zero);
 192     labelList nBoundary(pointFaces.size(), 0);
 193
 194     forAll(pointFaces, patchPointI)
 195     {
 196         const labelList& pFaces = pointFaces[patchPointI];
 197
 198         forAll(pFaces, pfI)
 199         {
 200             label faceI = pFaces[pfI];
 201
 202             if (isMasterFace.get(pp.addressing()[faceI]))
 203             {
 204                 avgBoundary[patchPointI] += pp[faceI].centre(points);
 205                 nBoundary[patchPointI]++;
 206             }
 207         }
 208     }
 209
 210     syncTools::syncPointList
 211     (
 212         mesh,
 213         pp.meshPoints(),
 214         avgBoundary,
 215         plusEqOp<point>(),  // combine op
 216         vector::zero        // null value
 217     );
 218     syncTools::syncPointList
 219     (
 220         mesh,
 221         pp.meshPoints(),
 222         nBoundary,
 223         plusEqOp<label>(),  // combine op
 224         0                   // null value
 225     );
 226
 227     forAll(avgBoundary, i)
 228     {
 229         avgBoundary[i] /= nBoundary[i];
 230     }
 231
 232
 233     // Get average position of internal face centres
 234     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 235
 236     vectorField avgInternal;
 237     labelList nInternal;
 238     {
 239         vectorField globalSum(mesh.nPoints(), vector::zero);
 240         labelList globalNum(mesh.nPoints(), 0);
 241
 242         // Note: no use of pointFaces
 243         const faceList& faces = mesh.faces();
 244
 245         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 246         {
 247             const face& f = faces[faceI];
 248             const point& fc = mesh.faceCentres()[faceI];
 249
 250             forAll(f, fp)
 251             {
 252                 globalSum[f[fp]] += fc;
 253                 globalNum[f[fp]]++;
 254             }
 255         }
 256
 257         // Count coupled faces as internal ones (but only once)
 258         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 259
 260         forAll(patches, patchI)
 261         {
 262             if
 263             (
 264                 patches[patchI].coupled()
 265              && refCast<const coupledPolyPatch>(patches[patchI]).owner()
 266             )
 267             {
 268                 const coupledPolyPatch& pp =
 269                     refCast<const coupledPolyPatch>(patches[patchI]);
 270
 271                 const vectorField::subField faceCentres = pp.faceCentres();
 272
 273                 forAll(pp, i)
 274                 {
 275                     const face& f = pp[i];
 276                     const point& fc = faceCentres[i];
 277
 278                     forAll(f, fp)
 279                     {
 280                         globalSum[f[fp]] += fc;
 281                         globalNum[f[fp]]++;
 282                     }
 283                 }
 284             }
 285         }
 286
 287         syncTools::syncPointList
 288         (
 289             mesh,
 290             globalSum,
 291             plusEqOp<vector>(), // combine op
 292             vector::zero        // null value
 293         );
 294         syncTools::syncPointList
 295         (
 296             mesh,
 297             globalNum,
 298             plusEqOp<label>(),  // combine op
 299             0                   // null value
 300         );
 301
 302         avgInternal.setSize(meshPoints.size());
 303         nInternal.setSize(meshPoints.size());
 304
 305         forAll(avgInternal, patchPointI)
 306         {
 307             label meshPointI = meshPoints[patchPointI];
 308
 309             nInternal[patchPointI] = globalNum[meshPointI];
 310
 311             if (nInternal[patchPointI] == 0)
 312             {
 313                 avgInternal[patchPointI] = globalSum[meshPointI];
 314             }
 315             else
 316             {
 317                 avgInternal[patchPointI] =
 318                     globalSum[meshPointI]
 319                   / nInternal[patchPointI];
 320             }
 321         }
 322     }
 323
 324
 325     // Precalculate any cell using mesh point (replacement of pointCells()[])
 326     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 327
 328     labelList anyCell(mesh.nPoints(), -1);
 329     forAll(mesh.faceNeighbour(), faceI)
 330     {
 331         label own = mesh.faceOwner()[faceI];
 332         const face& f = mesh.faces()[faceI];
 333
 334         forAll(f, fp)
 335         {
 336             anyCell[f[fp]] = own;
 337         }
 338     }
 339     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
 340     {
 341         label own = mesh.faceOwner()[faceI];
 342
 343         const face& f = mesh.faces()[faceI];
 344
 345         forAll(f, fp)
 346         {
 347             anyCell[f[fp]] = own;
 348         }
 349     }
 350
 351
 352     // Displacement to calculate.
 353     pointField patchDisp(meshPoints.size(), vector::zero);
 354
 355     forAll(pointFaces, i)
 356     {
 357         label meshPointI = meshPoints[i];
 358         const point& currentPos = pp.points()[meshPointI];
 359
 360         // Now we have the two average points: avgBoundary and avgInternal
 361         // and how many boundary/internal faces connect to the point
 362         // (nBoundary, nInternal)
 363         // Do some blending between the two.
 364         // Note: the following section has some reasoning behind it but the
 365         // blending factors can be experimented with.
 366
 367         point newPos;
 368
 369         if (!nonManifoldPoint.get(i))
 370         {
 371             // Points that are manifold. Weight the internal and boundary
 372             // by their number of faces and blend with
 373             scalar internalBlend = 0.1;
 374             scalar blend = 0.1;
 375
 376             point avgPos =
 377                 (
 378                    internalBlend*nInternal[i]*avgInternal[i]
 379                   +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
 380                 )
 381               / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
 382
 383             newPos = (1-blend)*avgPos + blend*currentPos;
 384         }
 385         else if (nInternal[i] == 0)
 386         {
 387             // Non-manifold without internal faces. Use any connected cell
 388             // as internal point instead. Use precalculated any cell to avoid
 389             // e.g. pointCells()[meshPointI][0]
 390
 391             const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
 392
 393             scalar cellCBlend = 0.8;
 394             scalar blend = 0.1;
 395
 396             point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
 397
 398             newPos = (1-blend)*avgPos + blend*currentPos;
 399         }
 400         else
 401         {
 402             // Non-manifold point with internal faces connected to them
 403             scalar internalBlend = 0.9;
 404             scalar blend = 0.1;
 405
 406             point avgPos =
 407                 internalBlend*avgInternal[i]
 408               + (1-internalBlend)*avgBoundary[i];
 409
 410             newPos = (1-blend)*avgPos + blend*currentPos;
 411         }
 412
 413         patchDisp[i] = newPos - currentPos;
 414     }
 415
 416     return patchDisp;
 417 }
 418
 419
 420 Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
 421 (
 422     const pointMesh& pMesh,
 423     const indirectPrimitivePatch& pp
 424 )
 425 {
 426     const polyMesh& mesh = pMesh();
 427
 428     // Set initial changed points to all the patch points
 429     List<pointEdgePoint> wallInfo(pp.nPoints());
 430
 431     forAll(pp.localPoints(), ppI)
 432     {
 433         wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
 434     }
 435
 436     // Current info on points
 437     List<pointEdgePoint> allPointInfo(mesh.nPoints());
 438
 439     // Current info on edges
 440     List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
 441
 442     PointEdgeWave<pointEdgePoint> wallCalc
 443     (
 444         mesh,
 445         pp.meshPoints(),
 446         wallInfo,
 447
 448         allPointInfo,
 449         allEdgeInfo,
 450         mesh.globalData().nTotalPoints()  // max iterations
 451     );
 452
 453     // Copy edge values into scalarField
 454     tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
 455     scalarField& edgeDist = tedgeDist();
 456
 457     forAll(allEdgeInfo, edgeI)
 458     {
 459         edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
 460     }
 461
 462
 463     //{
 464     //    // For debugging: dump to file
 465     //    pointScalarField pointDist
 466     //    (
 467     //        IOobject
 468     //        (
 469     //            "pointDist",
 470     //            meshRefiner_.timeName(),
 471     //            mesh.DB(),
 472     //            IOobject::NO_READ,
 473     //            IOobject::AUTO_WRITE
 474     //        ),
 475     //        pMesh,
 476     //        dimensionedScalar("pointDist", dimless, 0.0)
 477     //    );
 478     //
 479     //    forAll(allEdgeInfo, edgeI)
 480     //    {
 481     //        scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
 482     //
 483     //        const edge& e = mesh.edges()[edgeI];
 484     //
 485     //        pointDist[e[0]] += d;
 486     //        pointDist[e[1]] += d;
 487     //    }
 488     //    forAll(pointDist, pointI)
 489     //    {
 490     //        pointDist[pointI] /= mesh.pointEdges()[pointI].size();
 491     //    }
 492     //    Info<< "Writing patch distance to " << pointDist.name()
 493     //        << " at time " << meshRefiner_.timeName() << endl;
 494     //
 495     //    pointDist.write();
 496     //}
 497
 498     return tedgeDist;
 499 }
 500
 501
 502 void Foam::autoSnapDriver::dumpMove
 503 (
 504     const fileName& fName,
 505     const pointField& meshPts,
 506     const pointField& surfPts
 507 )
 508 {
 509     // Dump direction of growth into file
 510     Pout<< nl << "Dumping move direction to " << fName << endl;
 511
 512     OFstream nearestStream(fName);
 513
 514     label vertI = 0;
 515
 516     forAll(meshPts, ptI)
 517     {
 518         meshTools::writeOBJ(nearestStream, meshPts[ptI]);
 519         vertI++;
 520
 521         meshTools::writeOBJ(nearestStream, surfPts[ptI]);
 522         vertI++;
 523
 524         nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
 525     }
 526 }
 527
 528
 529 // Check whether all displacement vectors point outwards of patch. Return true
 530 // if so.
 531 bool Foam::autoSnapDriver::outwardsDisplacement
 532 (
 533     const indirectPrimitivePatch& pp,
 534     const vectorField& patchDisp
 535 )
 536 {
 537     const vectorField& faceNormals = pp.faceNormals();
 538     const labelListList& pointFaces = pp.pointFaces();
 539
 540     forAll(pointFaces, pointI)
 541     {
 542         const labelList& pFaces = pointFaces[pointI];
 543
 544         vector disp(patchDisp[pointI]);
 545
 546         scalar magDisp = mag(disp);
 547
 548         if (magDisp > SMALL)
 549         {
 550             disp /= magDisp;
 551
 552             bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
 553
 554             if (!outwards)
 555             {
 556                 Warning<< "Displacement " << patchDisp[pointI]
 557                     << " at mesh point " << pp.meshPoints()[pointI]
 558                     << " coord " << pp.points()[pp.meshPoints()[pointI]]
 559                     << " points through the surrounding patch faces" << endl;
 560                 return false;
 561             }
 562         }
 563         else
 564         {
 565             //? Displacement small but in wrong direction. Would probably be ok.
 566         }
 567     }
 568     return true;
 569 }
 570
 571
 572 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 573
 574 Foam::autoSnapDriver::autoSnapDriver
 575 (
 576     meshRefinement& meshRefiner,
 577     const labelList& globalToPatch
 578 )
 579 :
 580     meshRefiner_(meshRefiner),
 581     globalToPatch_(globalToPatch)
 582 {}
 583
 584
 585 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 586
 587 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
 588 (
 589     const List<labelPair>& baffles
 590 )
 591 {
 592     labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
 593
 594     autoPtr<mapPolyMesh> map;
 595
 596     // No need to sync; all processors will have all same zonedSurfaces.
 597     label nBaffles = returnReduce(baffles.size(), sumOp<label>());
 598     if (zonedSurfaces.size() && nBaffles > 0)
 599     {
 600         // Merge any baffles
 601         Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
 602             << endl;
 603
 604         map = meshRefiner_.mergeBaffles(baffles);
 605
 606         Info<< "Converted baffles in = "
 607             << meshRefiner_.mesh().time().cpuTimeIncrement()
 608             << " s\n" << nl << endl;
 609     }
 610
 611     return map;
 612 }
 613
 614
 615 Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 616 (
 617     const snapParameters& snapParams,
 618     const indirectPrimitivePatch& pp
 619 ) const
 620 {
 621     const edgeList& edges = pp.edges();
 622     const labelListList& pointEdges = pp.pointEdges();
 623     const pointField& localPoints = pp.localPoints();
 624     const fvMesh& mesh = meshRefiner_.mesh();
 625
 626     scalarField maxEdgeLen(localPoints.size(), -GREAT);
 627
 628     forAll(pointEdges, pointI)
 629     {
 630         const labelList& pEdges = pointEdges[pointI];
 631
 632         forAll(pEdges, pEdgeI)
 633         {
 634             const edge& e = edges[pEdges[pEdgeI]];
 635
 636             scalar len = e.mag(localPoints);
 637
 638             maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
 639         }
 640     }
 641
 642     syncTools::syncPointList
 643     (
 644         mesh,
 645         pp.meshPoints(),
 646         maxEdgeLen,
 647         maxEqOp<scalar>(),  // combine op
 648         -GREAT              // null value
 649     );
 650
 651     return scalarField(snapParams.snapTol()*maxEdgeLen);
 652 }
 653
 654
 655 void Foam::autoSnapDriver::preSmoothPatch
 656 (
 657     const snapParameters& snapParams,
 658     const label nInitErrors,
 659     const List<labelPair>& baffles,
 660     motionSmoother& meshMover
 661 ) const
 662 {
 663     const fvMesh& mesh = meshRefiner_.mesh();
 664
 665     labelList checkFaces;
 666
 667     Info<< "Smoothing patch points ..." << endl;
 668     for
 669     (
 670         label smoothIter = 0;
 671         smoothIter < snapParams.nSmoothPatch();
 672         smoothIter++
 673     )
 674     {
 675         Info<< "Smoothing iteration " << smoothIter << endl;
 676         checkFaces.setSize(mesh.nFaces());
 677         forAll(checkFaces, faceI)
 678         {
 679             checkFaces[faceI] = faceI;
 680         }
 681
 682         pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
 683
 684         // The current mesh is the starting mesh to smooth from.
 685         meshMover.setDisplacement(patchDisp);
 686
 687         meshMover.correct();
 688
 689         scalar oldErrorReduction = -1;
 690
 691         for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
 692         {
 693             Info<< nl << "Scaling iteration " << snapIter << endl;
 694
 695             if (snapIter == snapParams.nSnap())
 696             {
 697                 Info<< "Displacement scaling for error reduction set to 0."
 698                     << endl;
 699                 oldErrorReduction = meshMover.setErrorReduction(0.0);
 700             }
 701
 702             // Try to adapt mesh to obtain displacement by smoothly
 703             // decreasing displacement at error locations.
 704             if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
 705             {
 706                 Info<< "Successfully moved mesh" << endl;
 707                 break;
 708             }
 709         }
 710
 711         if (oldErrorReduction >= 0)
 712         {
 713             meshMover.setErrorReduction(oldErrorReduction);
 714         }
 715         Info<< endl;
 716     }
 717
 718
 719     // The current mesh is the starting mesh to smooth from.
 720     meshMover.correct();
 721
 722     if (debug)
 723     {
 724         const_cast<Time&>(mesh.time())++;
 725         Pout<< "Writing patch smoothed mesh to time "
 726             << meshRefiner_.timeName() << '.' << endl;
 727         meshRefiner_.write
 728         (
 729             debug,
 730             mesh.time().path()/meshRefiner_.timeName()
 731         );
 732         Pout<< "Dumped mesh in = "
 733             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 734     }
 735
 736     Info<< "Patch points smoothed in = "
 737         << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 738 }
 739
 740
 741 // Get (pp-local) indices of points that are both on zone and on patched surface
 742 Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 743 (
 744     const indirectPrimitivePatch& pp,
 745     const word& zoneName
 746 ) const
 747 {
 748     const fvMesh& mesh = meshRefiner_.mesh();
 749
 750     label zoneI = mesh.faceZones().findZoneID(zoneName);
 751
 752     if (zoneI == -1)
 753     {
 754         FatalErrorIn
 755         (
 756             "autoSnapDriver::getZoneSurfacePoints"
 757             "(const indirectPrimitivePatch&, const word&)"
 758         )   << "Cannot find zone " << zoneName
 759             << exit(FatalError);
 760     }
 761
 762     const faceZone& fZone = mesh.faceZones()[zoneI];
 763
 764
 765     // Could use PrimitivePatch & localFaces to extract points but might just
 766     // as well do it ourselves.
 767
 768     boolList pointOnZone(pp.nPoints(), false);
 769
 770     forAll(fZone, i)
 771     {
 772         const face& f = mesh.faces()[fZone[i]];
 773
 774         forAll(f, fp)
 775         {
 776             label meshPointI = f[fp];
 777
 778             Map<label>::const_iterator iter =
 779                 pp.meshPointMap().find(meshPointI);
 780
 781             if (iter != pp.meshPointMap().end())
 782             {
 783                 label pointI = iter();
 784                 pointOnZone[pointI] = true;
 785             }
 786         }
 787     }
 788
 789     return findIndices(pointOnZone, true);
 790 }
 791
 792
 793 Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 794 (
 795     const scalarField& snapDist,
 796     motionSmoother& meshMover
 797 ) const
 798 {
 799     Info<< "Calculating patchDisplacement as distance to nearest surface"
 800         << " point ..." << endl;
 801
 802     const indirectPrimitivePatch& pp = meshMover.patch();
 803     const pointField& localPoints = pp.localPoints();
 804     const refinementSurfaces& surfaces = meshRefiner_.surfaces();
 805     const fvMesh& mesh = meshRefiner_.mesh();
 806
 807     // Displacement per patch point
 808     vectorField patchDisp(localPoints.size(), vector::zero);
 809
 810     if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
 811     {
 812         // Current surface snapped to
 813         labelList snapSurf(localPoints.size(), -1);
 814
 815         // Divide surfaces into zoned and unzoned
 816         labelList zonedSurfaces =
 817             meshRefiner_.surfaces().getNamedSurfaces();
 818         labelList unzonedSurfaces =
 819             meshRefiner_.surfaces().getUnnamedSurfaces();
 820
 821
 822         // 1. All points to non-interface surfaces
 823         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 824
 825         {
 826             List<pointIndexHit> hitInfo;
 827             labelList hitSurface;
 828             surfaces.findNearest
 829             (
 830                 unzonedSurfaces,
 831                 localPoints,
 832                 sqr(snapDist),        // sqr of attract distance
 833                 hitSurface,
 834                 hitInfo
 835             );
 836
 837             forAll(hitInfo, pointI)
 838             {
 839                 if (hitInfo[pointI].hit())
 840                 {
 841                     patchDisp[pointI] =
 842                         hitInfo[pointI].hitPoint()
 843                       - localPoints[pointI];
 844
 845                     snapSurf[pointI] = hitSurface[pointI];
 846                 }
 847             }
 848         }
 849
 850
 851
 852         // 2. All points on zones to their respective surface
 853         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 854
 855         // Surfaces with zone information
 856         const wordList& faceZoneNames = surfaces.faceZoneNames();
 857
 858         // Current best snap distance
 859         scalarField minSnapDist(snapDist);
 860
 861         forAll(zonedSurfaces, i)
 862         {
 863             label zoneSurfI = zonedSurfaces[i];
 864
 865             const labelList surfacesToTest(1, zoneSurfI);
 866
 867             // Get indices of points both on faceZone and on pp.
 868             labelList zonePointIndices
 869             (
 870                 getZoneSurfacePoints
 871                 (
 872                     pp,
 873                     faceZoneNames[zoneSurfI]
 874                 )
 875             );
 876
 877             // Find nearest for points both on faceZone and pp.
 878             List<pointIndexHit> hitInfo;
 879             labelList hitSurface;
 880             surfaces.findNearest
 881             (
 882                 labelList(1, zoneSurfI),
 883                 pointField(localPoints, zonePointIndices),
 884                 sqr(scalarField(minSnapDist, zonePointIndices)),
 885                 hitSurface,
 886                 hitInfo
 887             );
 888
 889             forAll(hitInfo, i)
 890             {
 891                 label pointI = zonePointIndices[i];
 892
 893                 if (hitInfo[i].hit())
 894                 {
 895                     patchDisp[pointI] =
 896                         hitInfo[i].hitPoint()
 897                       - localPoints[pointI];
 898
 899                     minSnapDist[pointI] = min
 900                     (
 901                         minSnapDist[pointI],
 902                         mag(patchDisp[pointI])
 903                     );
 904
 905                     snapSurf[pointI] = zoneSurfI;
 906                 }
 907             }
 908         }
 909
 910         // Check if all points are being snapped
 911         forAll(snapSurf, pointI)
 912         {
 913             if (snapSurf[pointI] == -1)
 914             {
 915                 WarningIn("autoSnapDriver::calcNearestSurface(..)")
 916                     << "For point:" << pointI
 917                     << " coordinate:" << localPoints[pointI]
 918                     << " did not find any surface within:"
 919                     << minSnapDist[pointI]
 920                     << " meter." << endl;
 921             }
 922         }
 923
 924         {
 925             scalarField magDisp(mag(patchDisp));
 926
 927             Info<< "Wanted displacement : average:"
 928                 << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
 929                 << " min:" << gMin(magDisp)
 930                 << " max:" << gMax(magDisp) << endl;
 931         }
 932     }
 933
 934     Info<< "Calculated surface displacement in = "
 935         << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 936
 937
 938     // Limit amount of movement.
 939     forAll(patchDisp, patchPointI)
 940     {
 941         scalar magDisp = mag(patchDisp[patchPointI]);
 942
 943         if (magDisp > snapDist[patchPointI])
 944         {
 945             patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
 946
 947             Pout<< "Limiting displacement for " << patchPointI
 948                 << " from " << magDisp << " to " << snapDist[patchPointI]
 949                 << endl;
 950         }
 951     }
 952
 953     // Points on zones in one domain but only present as point on other
 954     // will not do condition 2 on all. Sync explicitly.
 955     syncTools::syncPointList
 956     (
 957         mesh,
 958         pp.meshPoints(),
 959         patchDisp,
 960         minMagSqrEqOp<point>(),         // combine op
 961         vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
 962     );
 963
 964     return patchDisp;
 965 }
 966
 967
 968 void Foam::autoSnapDriver::smoothDisplacement
 969 (
 970     const snapParameters& snapParams,
 971     motionSmoother& meshMover
 972 ) const
 973 {
 974     const fvMesh& mesh = meshRefiner_.mesh();
 975     const indirectPrimitivePatch& pp = meshMover.patch();
 976
 977     Info<< "Smoothing displacement ..." << endl;
 978
 979     // Set edge diffusivity as inverse of distance to patch
 980     scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
 981     //scalarField edgeGamma(mesh.nEdges(), 1.0);
 982     //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
 983
 984     // Get displacement field
 985     pointVectorField& disp = meshMover.displacement();
 986
 987     for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
 988     {
 989         if ((iter % 10) == 0)
 990         {
 991             Info<< "Iteration " << iter << endl;
 992         }
 993         pointVectorField oldDisp(disp);
 994         meshMover.smooth(oldDisp, edgeGamma, disp);
 995     }
 996     Info<< "Displacement smoothed in = "
 997         << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 998
 999     if (debug)
1000     {
1001         const_cast<Time&>(mesh.time())++;
1002         Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
1003             << endl;
1004
1005         // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
1006         // but this will also delete all pointMesh but not pointFields which
1007         // gives an illegal situation.
1008
1009         mesh.write();
1010
1011         Pout<< "Writing displacement field ..." << endl;
1012         disp.write();
1013         tmp<pointScalarField> magDisp(mag(disp));
1014         magDisp().write();
1015
1016         Pout<< "Writing actual patch displacement ..." << endl;
1017         vectorField actualPatchDisp(disp, pp.meshPoints());
1018         dumpMove
1019         (
1020             mesh.time().path()/"actualPatchDisplacement.obj",
1021             pp.localPoints(),
1022             pp.localPoints() + actualPatchDisp
1023         );
1024     }
1025 }
1026
1027
1028 bool Foam::autoSnapDriver::scaleMesh
1029 (
1030     const snapParameters& snapParams,
1031     const label nInitErrors,
1032     const List<labelPair>& baffles,
1033     motionSmoother& meshMover
1034 )
1035 {
1036     const fvMesh& mesh = meshRefiner_.mesh();
1037
1038     // Relax displacement until correct mesh
1039     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1040     labelList checkFaces(identity(mesh.nFaces()));
1041
1042     scalar oldErrorReduction = -1;
1043
1044     bool meshOk = false;
1045
1046     Info<< "Moving mesh ..." << endl;
1047     for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
1048     {
1049         Info<< nl << "Iteration " << iter << endl;
1050
1051         if (iter == snapParams.nSnap())
1052         {
1053             Info<< "Displacement scaling for error reduction set to 0." << endl;
1054             oldErrorReduction = meshMover.setErrorReduction(0.0);
1055         }
1056
1057         meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
1058
1059         if (meshOk)
1060         {
1061             Info<< "Successfully moved mesh" << endl;
1062             break;
1063         }
1064         if (debug)
1065         {
1066             const_cast<Time&>(mesh.time())++;
1067             Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName()
1068                 << endl;
1069             mesh.write();
1070
1071             Pout<< "Writing displacement field ..." << endl;
1072             meshMover.displacement().write();
1073             tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
1074             magDisp().write();
1075         }
1076     }
1077
1078     if (oldErrorReduction >= 0)
1079     {
1080         meshMover.setErrorReduction(oldErrorReduction);
1081     }
1082     Info<< "Moved mesh in = "
1083         << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1084
1085     return meshOk;
1086 }
1087
1088
1089 // After snapping: correct patching according to nearest surface.
1090 // Code is very similar to calcNearestSurface.
1091 // - calculate face-wise snap distance as max of point-wise
1092 // - calculate face-wise nearest surface point
1093 // - repatch face according to patch for surface point.
1094 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
1095 (
1096     const snapParameters& snapParams,
1097     const labelList& adaptPatchIDs
1098 )
1099 {
1100     const fvMesh& mesh = meshRefiner_.mesh();
1101     const refinementSurfaces& surfaces = meshRefiner_.surfaces();
1102
1103     Info<< "Repatching faces according to nearest surface ..." << endl;
1104
1105     // Get the labels of added patches.
1106     autoPtr<indirectPrimitivePatch> ppPtr
1107     (
1108         meshRefinement::makePatch
1109         (
1110             mesh,
1111             adaptPatchIDs
1112         )
1113     );
1114     indirectPrimitivePatch& pp = ppPtr();
1115
1116     // Divide surfaces into zoned and unzoned
1117     labelList zonedSurfaces = surfaces.getNamedSurfaces();
1118     labelList unzonedSurfaces = surfaces.getUnnamedSurfaces();
1119
1120
1121     // Faces that do not move
1122     PackedBoolList isZonedFace(mesh.nFaces());
1123     {
1124         // 1. All faces on zoned surfaces
1125         const wordList& faceZoneNames = surfaces.faceZoneNames();
1126         const faceZoneMesh& fZones = mesh.faceZones();
1127
1128         forAll(zonedSurfaces, i)
1129         {
1130             const label zoneSurfI = zonedSurfaces[i];
1131             const faceZone& fZone = fZones[faceZoneNames[zoneSurfI]];
1132
1133             forAll(fZone, i)
1134             {
1135                 isZonedFace.set(fZone[i], 1);
1136             }
1137         }
1138     }
1139
1140
1141     // Determine per pp face which patch it should be in
1142     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1143
1144     // Patch that face should be in
1145     labelList closestPatch(pp.size(), -1);
1146     {
1147         // face snap distance as max of point snap distance
1148         scalarField faceSnapDist(pp.size(), -GREAT);
1149         {
1150             // Distance to attract to nearest feature on surface
1151             const scalarField snapDist(calcSnapDistance(snapParams, pp));
1152
1153             const faceList& localFaces = pp.localFaces();
1154
1155             forAll(localFaces, faceI)
1156             {
1157                 const face& f = localFaces[faceI];
1158
1159                 forAll(f, fp)
1160                 {
1161                     faceSnapDist[faceI] = max
1162                     (
1163                         faceSnapDist[faceI],
1164                         snapDist[f[fp]]
1165                     );
1166                 }
1167             }
1168         }
1169
1170         pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
1171
1172         // Get nearest surface and region
1173         labelList hitSurface;
1174         labelList hitRegion;
1175         surfaces.findNearestRegion
1176         (
1177             unzonedSurfaces,
1178             localFaceCentres,
1179             sqr(faceSnapDist),    // sqr of attract distance
1180             hitSurface,
1181             hitRegion
1182         );
1183
1184         // Get patch
1185         forAll(pp, i)
1186         {
1187             label faceI = pp.addressing()[i];
1188
1189             if (hitSurface[i] != -1 && !isZonedFace.get(faceI))
1190             {
1191                 closestPatch[i] = globalToPatch_
1192                 [
1193                     surfaces.globalRegion
1194                     (
1195                         hitSurface[i],
1196                         hitRegion[i]
1197                     )
1198                 ];
1199             }
1200         }
1201     }
1202
1203
1204     // Change those faces for which there is a different closest patch
1205     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1206
1207     labelList ownPatch(mesh.nFaces(), -1);
1208     labelList neiPatch(mesh.nFaces(), -1);
1209
1210     const polyBoundaryMesh& patches = mesh.boundaryMesh();
1211
1212     forAll(patches, patchI)
1213     {
1214         const polyPatch& pp = patches[patchI];
1215
1216         forAll(pp, i)
1217         {
1218             ownPatch[pp.start()+i] = patchI;
1219             neiPatch[pp.start()+i] = patchI;
1220         }
1221     }
1222
1223     label nChanged = 0;
1224     forAll(closestPatch, i)
1225     {
1226         label faceI = pp.addressing()[i];
1227
1228         if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
1229         {
1230             ownPatch[faceI] = closestPatch[i];
1231             neiPatch[faceI] = closestPatch[i];
1232             nChanged++;
1233         }
1234     }
1235
1236     Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
1237         << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
1238         << endl;
1239
1240     return meshRefiner_.createBaffles(ownPatch, neiPatch);
1241 }
1242
1243
1244 void Foam::autoSnapDriver::doSnap
1245 (
1246     const dictionary& snapDict,
1247     const dictionary& motionDict,
1248     const scalar featureCos,
1249     const snapParameters& snapParams
1250 )
1251 {
1252     fvMesh& mesh = meshRefiner_.mesh();
1253
1254     Info<< nl
1255         << "Morphing phase" << nl
1256         << "--------------" << nl
1257         << endl;
1258
1259     // Get the labels of added patches.
1260     labelList adaptPatchIDs(meshRefiner_.meshedPatches());
1261
1262     // Create baffles (pairs of faces that share the same points)
1263     // Baffles stored as owner and neighbour face that have been created.
1264     List<labelPair> baffles;
1265     meshRefiner_.createZoneBaffles(globalToPatch_, baffles);
1266
1267
1268     bool doFeatures = false;
1269     label nFeatIter = 1;
1270     if (snapParams.nFeatureSnap() > 0)
1271     {
1272         doFeatures = true;
1273         nFeatIter = snapParams.nFeatureSnap();
1274
1275         Info<< "Snapping to features in " << nFeatIter
1276             << " iterations ..." << endl;
1277     }
1278
1279
1280     bool meshOk = false;
1281
1282     {
1283         autoPtr<indirectPrimitivePatch> ppPtr
1284         (
1285             meshRefinement::makePatch
1286             (
1287                 mesh,
1288                 adaptPatchIDs
1289             )
1290         );
1291
1292         // Distance to attract to nearest feature on surface
1293         const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
1294
1295
1296         // Construct iterative mesh mover.
1297         Info<< "Constructing mesh displacer ..." << endl;
1298         Info<< "Using mesh parameters " << motionDict << nl << endl;
1299
1300         const pointMesh& pMesh = pointMesh::New(mesh);
1301
1302         motionSmoother meshMover
1303         (
1304             mesh,
1305             ppPtr(),
1306             adaptPatchIDs,
1307             meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
1308             motionDict
1309         );
1310
1311
1312         // Check initial mesh
1313         Info<< "Checking initial mesh ..." << endl;
1314         labelHashSet wrongFaces(mesh.nFaces()/100);
1315         motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
1316         const label nInitErrors = returnReduce
1317         (
1318             wrongFaces.size(),
1319             sumOp<label>()
1320         );
1321
1322         Info<< "Detected " << nInitErrors << " illegal faces"
1323             << " (concave, zero area or negative cell pyramid volume)"
1324             << endl;
1325
1326
1327         Info<< "Checked initial mesh in = "
1328             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1329
1330         // Pre-smooth patch vertices (so before determining nearest)
1331         preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
1332
1333
1334         for (label iter = 0; iter < nFeatIter; iter++)
1335         {
1336             Info<< nl
1337                 << "Morph iteration " << iter << nl
1338                 << "-----------------" << endl;
1339
1340             // Calculate displacement at every patch point. Insert into
1341             // meshMover.
1342             vectorField disp = calcNearestSurface(snapDist, meshMover);
1343
1344             // Override displacement with feature edge attempt
1345             if (doFeatures)
1346             {
1347                 disp = calcNearestSurfaceFeature
1348                 (
1349                     iter,
1350                     featureCos,
1351                     scalar(iter+1)/nFeatIter,
1352                     snapDist,
1353                     disp,
1354                     meshMover
1355                 );
1356             }
1357
1358             // Check for displacement being outwards.
1359             outwardsDisplacement(ppPtr(), disp);
1360
1361             // Set initial distribution of displacement field (on patches)
1362             // from patchDisp and make displacement consistent with b.c.
1363             // on displacement pointVectorField.
1364             meshMover.setDisplacement(disp);
1365
1366
1367             if (debug&meshRefinement::OBJINTERSECTIONS)
1368             {
1369                 dumpMove
1370                 (
1371                     mesh.time().path()
1372                   / "patchDisplacement_" + name(iter) + ".obj",
1373                     ppPtr().localPoints(),
1374                     ppPtr().localPoints() + disp
1375                 );
1376             }
1377
1378             // Get smoothly varying internal displacement field.
1379             smoothDisplacement(snapParams, meshMover);
1380
1381             // Apply internal displacement to mesh.
1382             meshOk = scaleMesh
1383             (
1384                 snapParams,
1385                 nInitErrors,
1386                 baffles,
1387                 meshMover
1388             );
1389
1390             if (!meshOk)
1391             {
1392                 Info<< "Did not succesfully snap mesh. Giving up."
1393                     << nl << endl;
1394
1395                 // Use current mesh as base mesh
1396                 meshMover.correct();
1397
1398                 break;
1399             }
1400
1401             if (debug)
1402             {
1403                 const_cast<Time&>(mesh.time())++;
1404                 Pout<< "Writing scaled mesh to time "
1405                     << meshRefiner_.timeName() << endl;
1406                 meshRefiner_.write
1407                 (
1408                     debug,
1409                     mesh.time().path()/meshRefiner_.timeName()
1410                 );
1411                 Pout<< "Writing displacement field ..." << endl;
1412                 meshMover.displacement().write();
1413                 tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
1414                 magDisp().write();
1415             }
1416
1417             // Use current mesh as base mesh
1418             meshMover.correct();
1419         }
1420     }
1421
1422     // Merge any introduced baffles.
1423     mergeZoneBaffles(baffles);
1424
1425     // Repatch faces according to nearest.
1426     repatchToSurface(snapParams, adaptPatchIDs);
1427
1428     // Repatching might have caused faces to be on same patch and hence
1429     // mergeable so try again to merge coplanar faces
1430     label nChanged = meshRefiner_.mergePatchFacesUndo
1431     (
1432         featureCos,  // minCos
1433         featureCos,  // concaveCos
1434         meshRefiner_.meshedPatches(),
1435         motionDict
1436     );
1437
1438     nChanged += meshRefiner_.mergeEdgesUndo
1439     (
1440         featureCos,
1441         motionDict
1442     );
1443
1444     if (nChanged > 0 && debug)
1445     {
1446         const_cast<Time&>(mesh.time())++;
1447         Pout<< "Writing patchFace merged mesh to time "
1448             << meshRefiner_.timeName() << endl;
1449         meshRefiner_.write
1450         (
1451             debug,
1452             meshRefiner_.timeName()
1453         );
1454     }
1455 }
1456
1457
1458 // ************************************************************************* //