applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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 "objectRegistry.H"
  34 #include "foamTime.H"
  35 #include "fvMesh.H"
  36 #include "autoRefineDriver.H"
  37 #include "autoSnapDriver.H"
  38 #include "autoLayerDriver.H"
  39 #include "searchableSurfaces.H"
  40 #include "refinementSurfaces.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     scalar writeTol = std::pow
  60     (
  61         scalar(10.0),
  62        -scalar(IOstream::defaultPrecision())
  63     );
  64
  65     Info<< nl
  66         << "Overall mesh bounding box  : " << meshBb << nl
  67         << "Relative tolerance         : " << mergeTol << nl
  68         << "Absolute matching distance : " << mergeDist << nl
  69         << endl;
  70
  71     if (mesh.time().writeFormat() == IOstream::ASCII && mergeTol < writeTol)
  72     {
  73         FatalErrorIn("getMergeDistance(const polyMesh&, const scalar)")
  74             << "Your current settings specify ASCII writing with "
  75             << IOstream::defaultPrecision() << " digits precision." << endl
  76             << "Your merging tolerance (" << mergeTol << ") is finer than this."
  77             << endl
  78             << "Please change your writeFormat to binary"
  79             << " or increase the writePrecision" << endl
  80             << "or adjust the merge tolerance (-mergeTol)."
  81             << exit(FatalError);
  82     }
  83
  84     return mergeDist;
  85 }
  86
  87
  88 // Write mesh and additional information
  89 void writeMesh
  90 (
  91     const string& msg,
  92     const meshRefinement& meshRefiner,
  93     const label debug
  94 )
  95 {
  96     const fvMesh& mesh = meshRefiner.mesh();
  97
  98     meshRefiner.printMeshInfo(debug, msg);
  99     Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
 100
 101     meshRefiner.write(meshRefinement::MESH|meshRefinement::SCALARLEVELS, "");
 102     if (debug & meshRefinement::OBJINTERSECTIONS)
 103     {
 104         meshRefiner.write
 105         (
 106             meshRefinement::OBJINTERSECTIONS,
 107             mesh.time().path()/meshRefiner.timeName()
 108         );
 109     }
 110     Info<< "Written mesh in = "
 111         << mesh.time().cpuTimeIncrement() << " s." << endl;
 112 }
 113
 114
 115
 116 int main(int argc, char *argv[])
 117 {
 118     argList::validOptions.insert("overwrite", "");
 119 #   include "setRootCase.H"
 120 #   include "createTime.H"
 121     runTime.functionObjects().off();
 122 #   include "createMesh.H"
 123
 124     Info<< "Read mesh in = "
 125         << runTime.cpuTimeIncrement() << " s" << endl;
 126
 127     const bool overwrite = args.optionFound("overwrite");
 128
 129
 130     // Check patches and faceZones are synchronised
 131     mesh.boundaryMesh().checkParallelSync(true);
 132     meshRefinement::checkCoupledFaceZones(mesh);
 133
 134
 135     // Read decomposePar dictionary
 136     IOdictionary decomposeDict
 137     (
 138         IOobject
 139         (
 140             "decomposeParDict",
 141             runTime.system(),
 142             mesh,
 143             IOobject::MUST_READ,
 144             IOobject::NO_WRITE
 145         )
 146     );
 147
 148     // Read meshing dictionary
 149     IOdictionary meshDict
 150     (
 151        IOobject
 152        (
 153             "snappyHexMeshDict",
 154             runTime.system(),
 155             mesh,
 156             IOobject::MUST_READ,
 157             IOobject::NO_WRITE
 158        )
 159     );
 160
 161     // all surface geometry
 162     const dictionary& geometryDict = meshDict.subDict("geometry");
 163
 164     // refinement parameters
 165     const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");
 166
 167     // mesh motion and mesh quality parameters
 168     const dictionary& motionDict = meshDict.subDict("meshQualityControls");
 169
 170     // snap-to-surface parameters
 171     const dictionary& snapDict = meshDict.subDict("snapControls");
 172
 173     // layer addition parameters
 174     const dictionary& layerDict = meshDict.subDict("addLayersControls");
 175
 176
 177     const scalar mergeDist = getMergeDistance
 178     (
 179         mesh,
 180         readScalar(meshDict.lookup("mergeTolerance"))
 181     );
 182
 183
 184
 185     // Debug
 186     // ~~~~~
 187
 188     const label debug(readLabel(meshDict.lookup("debug")));
 189     if (debug > 0)
 190     {
 191         meshRefinement::debug = debug;
 192         autoRefineDriver::debug = debug;
 193         autoSnapDriver::debug = debug;
 194         autoLayerDriver::debug = debug;
 195     }
 196
 197
 198     // Read geometry
 199     // ~~~~~~~~~~~~~
 200
 201     searchableSurfaces allGeometry
 202     (
 203         IOobject
 204         (
 205             "abc",                      // dummy name
 206             mesh.time().constant(),     // instance
 207             //mesh.time().findInstance("triSurface", word::null),// instance
 208             "triSurface",               // local
 209             mesh.time(),                // registry
 210             IOobject::MUST_READ,
 211             IOobject::NO_WRITE
 212         ),
 213         geometryDict
 214     );
 215
 216
 217     // Read refinement surfaces
 218     // ~~~~~~~~~~~~~~~~~~~~~~~~
 219
 220     Info<< "Reading refinement surfaces." << endl;
 221     refinementSurfaces surfaces
 222     (
 223         allGeometry,
 224         refineDict.subDict("refinementSurfaces")
 225     );
 226     Info<< "Read refinement surfaces in = "
 227         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 228
 229
 230     // Read refinement shells
 231     // ~~~~~~~~~~~~~~~~~~~~~~
 232
 233     Info<< "Reading refinement shells." << endl;
 234     shellSurfaces shells
 235     (
 236         allGeometry,
 237         refineDict.subDict("refinementRegions")
 238     );
 239     Info<< "Read refinement shells in = "
 240         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 241
 242
 243     Info<< "Setting refinement level of surface to be consistent"
 244         << " with shells." << endl;
 245     surfaces.setMinLevelFields(shells);
 246     Info<< "Checked shell refinement in = "
 247         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 248
 249
 250     // Refinement engine
 251     // ~~~~~~~~~~~~~~~~~
 252
 253     Info<< nl
 254         << "Determining initial surface intersections" << nl
 255         << "-----------------------------------------" << nl
 256         << endl;
 257
 258     // Main refinement engine
 259     meshRefinement meshRefiner
 260     (
 261         mesh,
 262         mergeDist,          // tolerance used in sorting coordinates
 263         overwrite,          // overwrite mesh files?
 264         surfaces,           // for surface intersection refinement
 265         shells              // for volume (inside/outside) refinement
 266     );
 267     Info<< "Calculated surface intersections in = "
 268         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 269
 270     // Some stats
 271     meshRefiner.printMeshInfo(debug, "Initial mesh");
 272
 273     meshRefiner.write
 274     (
 275         debug&meshRefinement::OBJINTERSECTIONS,
 276         mesh.time().path()/meshRefiner.timeName()
 277     );
 278
 279
 280     // Add all the surface regions as patches
 281     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 282
 283     labelList globalToPatch;
 284     {
 285         Info<< nl
 286             << "Adding patches for surface regions" << nl
 287             << "----------------------------------" << nl
 288             << endl;
 289
 290         // From global region number to mesh patch.
 291         globalToPatch.setSize(surfaces.nRegions(), -1);
 292
 293         Info<< "Patch\tRegion" << nl
 294             << "-----\t------"
 295             << endl;
 296
 297         const labelList& surfaceGeometry = surfaces.surfaces();
 298         forAll(surfaceGeometry, surfI)
 299         {
 300             label geomI = surfaceGeometry[surfI];
 301
 302             const wordList& regNames = allGeometry.regionNames()[geomI];
 303
 304             Info<< surfaces.names()[surfI] << ':' << nl << nl;
 305
 306             forAll(regNames, i)
 307             {
 308                 label patchI = meshRefiner.addMeshedPatch
 309                 (
 310                     regNames[i],
 311                     wallPolyPatch::typeName
 312                 );
 313
 314                 Info<< patchI << '\t' << regNames[i] << nl;
 315
 316                 globalToPatch[surfaces.globalRegion(surfI, i)] = patchI;
 317             }
 318
 319             Info<< nl;
 320         }
 321         Info<< "Added patches in = "
 322             << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 323     }
 324
 325
 326     // Parallel
 327     // ~~~~~~~~
 328
 329     // Decomposition
 330     autoPtr<decompositionMethod> decomposerPtr
 331     (
 332         decompositionMethod::New
 333         (
 334             decomposeDict,
 335             mesh
 336         )
 337     );
 338     decompositionMethod& decomposer = decomposerPtr();
 339
 340     if (Pstream::parRun() && !decomposer.parallelAware())
 341     {
 342         FatalErrorIn(args.executable())
 343             << "You have selected decomposition method "
 344             << decomposer.typeName
 345             << " which is not parallel aware." << endl
 346             << "Please select one that is (hierarchical, parMetis)"
 347             << exit(FatalError);
 348     }
 349
 350     // Mesh distribution engine (uses tolerance to reconstruct meshes)
 351     fvMeshDistribute distributor(mesh, mergeDist);
 352
 353
 354
 355
 356
 357     // Now do the real work -refinement -snapping -layers
 358     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 359
 360     Switch wantRefine(meshDict.lookup("castellatedMesh"));
 361     Switch wantSnap(meshDict.lookup("snap"));
 362     Switch wantLayers(meshDict.lookup("addLayers"));
 363
 364     if (wantRefine)
 365     {
 366         cpuTime timer;
 367
 368         autoRefineDriver refineDriver
 369         (
 370             meshRefiner,
 371             decomposer,
 372             distributor,
 373             globalToPatch
 374         );
 375
 376         // Refinement parameters
 377         refinementParameters refineParams(refineDict);
 378
 379         if (!overwrite)
 380         {
 381             const_cast<Time&>(mesh.time())++;
 382         }
 383
 384         refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
 385
 386         writeMesh
 387         (
 388             "Refined mesh",
 389             meshRefiner,
 390             debug
 391         );
 392
 393         Info<< "Mesh refined in = "
 394             << timer.cpuTimeIncrement() << " s." << endl;
 395     }
 396
 397     if (wantSnap)
 398     {
 399         cpuTime timer;
 400
 401         autoSnapDriver snapDriver
 402         (
 403             meshRefiner,
 404             globalToPatch
 405         );
 406
 407         // Snap parameters
 408         snapParameters snapParams(snapDict);
 409
 410         if (!overwrite)
 411         {
 412             const_cast<Time&>(mesh.time())++;
 413         }
 414
 415         snapDriver.doSnap(snapDict, motionDict, snapParams);
 416
 417         writeMesh
 418         (
 419             "Snapped mesh",
 420             meshRefiner,
 421             debug
 422         );
 423
 424         Info<< "Mesh snapped in = "
 425             << timer.cpuTimeIncrement() << " s." << endl;
 426     }
 427
 428     if (wantLayers)
 429     {
 430         cpuTime timer;
 431
 432         autoLayerDriver layerDriver(meshRefiner);
 433
 434         // Layer addition parameters
 435         layerParameters layerParams(layerDict, mesh.boundaryMesh());
 436
 437         //!!! Temporary hack to get access to maxLocalCells
 438         bool preBalance;
 439         {
 440             refinementParameters refineParams(refineDict);
 441
 442             preBalance = returnReduce
 443             (
 444                 (mesh.nCells() >= refineParams.maxLocalCells()),
 445                 orOp<bool>()
 446             );
 447         }
 448
 449
 450         if (!overwrite)
 451         {
 452             const_cast<Time&>(mesh.time())++;
 453         }
 454
 455         layerDriver.doLayers
 456         (
 457             layerDict,
 458             motionDict,
 459             layerParams,
 460             preBalance,
 461             decomposer,
 462             distributor
 463         );
 464
 465         writeMesh
 466         (
 467             "Layer mesh",
 468             meshRefiner,
 469             debug
 470         );
 471
 472         Info<< "Layers added in = "
 473             << timer.cpuTimeIncrement() << " s." << endl;
 474     }
 475
 476
 477     Info<< "Finished meshing in = "
 478         << runTime.elapsedCpuTime() << " s." << endl;
 479
 480     Info<< "End\n" << endl;
 481
 482     return(0);
 483 }
 484
 485
 486 // ************************************************************************* //