applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.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 Application
  25     snappyHexMesh
  26
  27 Description
  28     Automatic split hex mesher. Refines and snaps to surface.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "argList.H"
  33 #include "Time.H"
  34 #include "fvMesh.H"
  35 #include "autoRefineDriver.H"
  36 #include "autoSnapDriver.H"
  37 #include "autoLayerDriver.H"
  38 #include "searchableSurfaces.H"
  39 #include "refinementSurfaces.H"
  40 #include "refinementFeatures.H"
  41 #include "shellSurfaces.H"
  42 #include "decompositionMethod.H"
  43 #include "fvMeshDistribute.H"
  44 #include "wallPolyPatch.H"
  45 #include "refinementParameters.H"
  46 #include "snapParameters.H"
  47 #include "layerParameters.H"
  48
  49
  50 using namespace Foam;
  51
  52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  53
  54 // Check writing tolerance before doing any serious work
  55 scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
  56 {
  57     const boundBox& meshBb = mesh.bounds();
  58     scalar mergeDist = mergeTol * meshBb.mag();
  59
  60     Info<< nl
  61         << "Overall mesh bounding box  : " << meshBb << nl
  62         << "Relative tolerance         : " << mergeTol << nl
  63         << "Absolute matching distance : " << mergeDist << nl
  64         << endl;
  65
  66     // check writing tolerance
  67     if (mesh.time().writeFormat() == IOstream::ASCII)
  68     {
  69         const scalar writeTol = std::pow
  70         (
  71             scalar(10.0),
  72             -scalar(IOstream::defaultPrecision())
  73         );
  74
  75         if (mergeTol < writeTol)
  76         {
  77             FatalErrorIn("getMergeDistance(const polyMesh&, const dictionary&)")
  78                 << "Your current settings specify ASCII writing with "
  79                 << IOstream::defaultPrecision() << " digits precision." << nl
  80                 << "Your merging tolerance (" << mergeTol
  81                 << ") is finer than this." << nl
  82                 << "Change to binary writeFormat, "
  83                 << "or increase the writePrecision" << endl
  84                 << "or adjust the merge tolerance (mergeTol)."
  85                 << exit(FatalError);
  86         }
  87     }
  88
  89     return mergeDist;
  90 }
  91
  92
  93 // Write mesh and additional information
  94 void writeMesh
  95 (
  96     const string& msg,
  97     const meshRefinement& meshRefiner,
  98     const label debug
  99 )
 100 {
 101     const fvMesh& mesh = meshRefiner.mesh();
 102
 103     meshRefiner.printMeshInfo(debug, msg);
 104     Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
 105
 106     meshRefiner.write(meshRefinement::MESH|meshRefinement::SCALARLEVELS, "");
 107     if (debug & meshRefinement::OBJINTERSECTIONS)
 108     {
 109         meshRefiner.write
 110         (
 111             meshRefinement::OBJINTERSECTIONS,
 112             mesh.time().path()/meshRefiner.timeName()
 113         );
 114     }
 115     Info<< "Wrote mesh in = "
 116         << mesh.time().cpuTimeIncrement() << " s." << endl;
 117 }
 118
 119
 120
 121 int main(int argc, char *argv[])
 122 {
 123 #   include "addOverwriteOption.H"
 124
 125 #   include "setRootCase.H"
 126 #   include "createTime.H"
 127     runTime.functionObjects().off();
 128 #   include "createMesh.H"
 129
 130     Info<< "Read mesh in = "
 131         << runTime.cpuTimeIncrement() << " s" << endl;
 132
 133     const bool overwrite = args.optionFound("overwrite");
 134
 135     // Check patches and faceZones are synchronised
 136     mesh.boundaryMesh().checkParallelSync(true);
 137     meshRefinement::checkCoupledFaceZones(mesh);
 138
 139
 140     // Read decomposePar dictionary
 141     IOdictionary decomposeDict
 142     (
 143         IOobject
 144         (
 145             "decomposeParDict",
 146             runTime.system(),
 147             mesh,
 148             IOobject::MUST_READ_IF_MODIFIED,
 149             IOobject::NO_WRITE
 150         )
 151     );
 152
 153     // Read meshing dictionary
 154     IOdictionary meshDict
 155     (
 156        IOobject
 157        (
 158             "snappyHexMeshDict",
 159             runTime.system(),
 160             mesh,
 161             IOobject::MUST_READ_IF_MODIFIED,
 162             IOobject::NO_WRITE
 163        )
 164     );
 165
 166     // all surface geometry
 167     const dictionary& geometryDict = meshDict.subDict("geometry");
 168
 169     // refinement parameters
 170     const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");
 171
 172     // mesh motion and mesh quality parameters
 173     const dictionary& motionDict = meshDict.subDict("meshQualityControls");
 174
 175     // snap-to-surface parameters
 176     const dictionary& snapDict = meshDict.subDict("snapControls");
 177
 178     // layer addition parameters
 179     const dictionary& layerDict = meshDict.subDict("addLayersControls");
 180
 181     // absolute merge distance
 182     const scalar mergeDist = getMergeDistance
 183     (
 184         mesh,
 185         readScalar(meshDict.lookup("mergeTolerance"))
 186     );
 187
 188
 189     // Debug
 190     // ~~~~~
 191
 192     const label debug = meshDict.lookupOrDefault<label>("debug", 0);
 193     if (debug > 0)
 194     {
 195         meshRefinement::debug   = debug;
 196         autoRefineDriver::debug = debug;
 197         autoSnapDriver::debug   = debug;
 198         autoLayerDriver::debug  = debug;
 199     }
 200
 201
 202     // Read geometry
 203     // ~~~~~~~~~~~~~
 204
 205     searchableSurfaces allGeometry
 206     (
 207         IOobject
 208         (
 209             "abc",                      // dummy name
 210             mesh.time().constant(),     // instance
 211             //mesh.time().findInstance("triSurface", word::null),// instance
 212             "triSurface",               // local
 213             mesh.time(),                // registry
 214             IOobject::MUST_READ,
 215             IOobject::NO_WRITE
 216         ),
 217         geometryDict
 218     );
 219
 220
 221     // Read refinement surfaces
 222     // ~~~~~~~~~~~~~~~~~~~~~~~~
 223
 224     Info<< "Reading refinement surfaces." << endl;
 225     refinementSurfaces surfaces
 226     (
 227         allGeometry,
 228         refineDict.subDict("refinementSurfaces")
 229     );
 230     Info<< "Read refinement surfaces in = "
 231         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 232
 233
 234     // Read refinement shells
 235     // ~~~~~~~~~~~~~~~~~~~~~~
 236
 237     Info<< "Reading refinement shells." << endl;
 238     shellSurfaces shells
 239     (
 240         allGeometry,
 241         refineDict.subDict("refinementRegions")
 242     );
 243     Info<< "Read refinement shells in = "
 244         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 245
 246
 247     Info<< "Setting refinement level of surface to be consistent"
 248         << " with shells." << endl;
 249     surfaces.setMinLevelFields(shells);
 250     Info<< "Checked shell refinement in = "
 251         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 252
 253
 254     // Read feature meshes
 255     // ~~~~~~~~~~~~~~~~~~~
 256
 257     Info<< "Reading features." << endl;
 258     refinementFeatures features
 259     (
 260         mesh,
 261         refineDict.lookup("features")
 262     );
 263     Info<< "Read features in = "
 264         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 265
 266
 267
 268     // Refinement engine
 269     // ~~~~~~~~~~~~~~~~~
 270
 271     Info<< nl
 272         << "Determining initial surface intersections" << nl
 273         << "-----------------------------------------" << nl
 274         << endl;
 275
 276     // Main refinement engine
 277     meshRefinement meshRefiner
 278     (
 279         mesh,
 280         mergeDist,          // tolerance used in sorting coordinates
 281         overwrite,          // overwrite mesh files?
 282         surfaces,           // for surface intersection refinement
 283         features,           // for feature edges/point based refinement
 284         shells              // for volume (inside/outside) refinement
 285     );
 286     Info<< "Calculated surface intersections in = "
 287         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 288
 289     // Some stats
 290     meshRefiner.printMeshInfo(debug, "Initial mesh");
 291
 292     meshRefiner.write
 293     (
 294         debug & meshRefinement::OBJINTERSECTIONS,
 295         mesh.time().path()/meshRefiner.timeName()
 296     );
 297
 298
 299     // Add all the surface regions as patches
 300     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 301
 302     labelList globalToPatch;
 303     {
 304         Info<< nl
 305             << "Adding patches for surface regions" << nl
 306             << "----------------------------------" << nl
 307             << endl;
 308
 309         // From global region number to mesh patch.
 310         globalToPatch.setSize(surfaces.nRegions(), -1);
 311
 312         Info<< "Patch\tType\tRegion" << nl
 313             << "-----\t----\t------"
 314             << endl;
 315
 316         const labelList& surfaceGeometry = surfaces.surfaces();
 317         const PtrList<dictionary>& surfacePatchInfo = surfaces.patchInfo();
 318
 319         forAll(surfaceGeometry, surfI)
 320         {
 321             label geomI = surfaceGeometry[surfI];
 322
 323             const wordList& regNames = allGeometry.regionNames()[geomI];
 324
 325             Info<< surfaces.names()[surfI] << ':' << nl << nl;
 326
 327             forAll(regNames, i)
 328             {
 329                 label globalRegionI = surfaces.globalRegion(surfI, i);
 330
 331                 label patchI;
 332
 333                 if (surfacePatchInfo.set(globalRegionI))
 334                 {
 335                     patchI = meshRefiner.addMeshedPatch
 336                     (
 337                         regNames[i],
 338                         surfacePatchInfo[globalRegionI]
 339                     );
 340                 }
 341                 else
 342                 {
 343                     dictionary patchInfo;
 344                     patchInfo.set("type", wallPolyPatch::typeName);
 345
 346                     patchI = meshRefiner.addMeshedPatch
 347                     (
 348                         regNames[i],
 349                         patchInfo
 350                     );
 351                 }
 352
 353                 Info<< patchI << '\t' << mesh.boundaryMesh()[patchI].type()
 354                     << '\t' << regNames[i] << nl;
 355
 356                 globalToPatch[globalRegionI] = patchI;
 357             }
 358
 359             Info<< nl;
 360         }
 361         Info<< "Added patches in = "
 362             << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 363     }
 364
 365
 366     // Parallel
 367     // ~~~~~~~~
 368
 369     // Decomposition
 370     autoPtr<decompositionMethod> decomposerPtr
 371     (
 372         decompositionMethod::New
 373         (
 374             decomposeDict
 375         )
 376     );
 377     decompositionMethod& decomposer = decomposerPtr();
 378
 379     if (Pstream::parRun() && !decomposer.parallelAware())
 380     {
 381         FatalErrorIn(args.executable())
 382             << "You have selected decomposition method "
 383             << decomposer.typeName
 384             << " which is not parallel aware." << endl
 385             << "Please select one that is (hierarchical, ptscotch)"
 386             << exit(FatalError);
 387     }
 388
 389     // Mesh distribution engine (uses tolerance to reconstruct meshes)
 390     fvMeshDistribute distributor(mesh, mergeDist);
 391
 392
 393
 394
 395
 396     // Now do the real work -refinement -snapping -layers
 397     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 398
 399     const Switch wantRefine(meshDict.lookup("castellatedMesh"));
 400     const Switch wantSnap(meshDict.lookup("snap"));
 401     const Switch wantLayers(meshDict.lookup("addLayers"));
 402
 403     if (wantRefine)
 404     {
 405         cpuTime timer;
 406
 407         autoRefineDriver refineDriver
 408         (
 409             meshRefiner,
 410             decomposer,
 411             distributor,
 412             globalToPatch
 413         );
 414
 415         // Refinement parameters
 416         refinementParameters refineParams(refineDict);
 417
 418         if (!overwrite)
 419         {
 420             const_cast<Time&>(mesh.time())++;
 421         }
 422
 423         refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
 424
 425         writeMesh
 426         (
 427             "Refined mesh",
 428             meshRefiner,
 429             debug
 430         );
 431
 432         Info<< "Mesh refined in = "
 433             << timer.cpuTimeIncrement() << " s." << endl;
 434     }
 435
 436     if (wantSnap)
 437     {
 438         cpuTime timer;
 439
 440         autoSnapDriver snapDriver
 441         (
 442             meshRefiner,
 443             globalToPatch
 444         );
 445
 446         // Snap parameters
 447         snapParameters snapParams(snapDict);
 448         // Temporary hack to get access to resolveFeatureAngle
 449         scalar curvature;
 450         {
 451             refinementParameters refineParams(refineDict);
 452             curvature = refineParams.curvature();
 453         }
 454
 455         if (!overwrite)
 456         {
 457             const_cast<Time&>(mesh.time())++;
 458         }
 459
 460         snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);
 461
 462         writeMesh
 463         (
 464             "Snapped mesh",
 465             meshRefiner,
 466             debug
 467         );
 468
 469         Info<< "Mesh snapped in = "
 470             << timer.cpuTimeIncrement() << " s." << endl;
 471     }
 472
 473     if (wantLayers)
 474     {
 475         cpuTime timer;
 476
 477         autoLayerDriver layerDriver(meshRefiner, globalToPatch);
 478
 479         // Layer addition parameters
 480         layerParameters layerParams(layerDict, mesh.boundaryMesh());
 481
 482         //!!! Temporary hack to get access to maxLocalCells
 483         bool preBalance;
 484         {
 485             refinementParameters refineParams(refineDict);
 486
 487             preBalance = returnReduce
 488             (
 489                 (mesh.nCells() >= refineParams.maxLocalCells()),
 490                 orOp<bool>()
 491             );
 492         }
 493
 494
 495         if (!overwrite)
 496         {
 497             const_cast<Time&>(mesh.time())++;
 498         }
 499
 500         layerDriver.doLayers
 501         (
 502             layerDict,
 503             motionDict,
 504             layerParams,
 505             preBalance,
 506             decomposer,
 507             distributor
 508         );
 509
 510         writeMesh
 511         (
 512             "Layer mesh",
 513             meshRefiner,
 514             debug
 515         );
 516
 517         Info<< "Layers added in = "
 518             << timer.cpuTimeIncrement() << " s." << endl;
 519     }
 520
 521
 522     Info<< "Finished meshing in = "
 523         << runTime.elapsedCpuTime() << " s." << endl;
 524
 525     Info<< "End\n" << endl;
 526
 527     return 0;
 528 }
 529
 530
 531 // ************************************************************************* //