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) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Application
  26     snappyHexMesh
  27
  28 Description
  29     Automatic split hex mesher. Refines and snaps to surface.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "argList.H"
  34 #include "Time.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         autoRefineDriver refineDriver
 367         (
 368             meshRefiner,
 369             decomposer,
 370             distributor,
 371             globalToPatch
 372         );
 373
 374         // Refinement parameters
 375         refinementParameters refineParams(refineDict);
 376
 377         if (!overwrite)
 378         {
 379             const_cast<Time&>(mesh.time())++;
 380         }
 381
 382         refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
 383
 384         writeMesh
 385         (
 386             "Refined mesh",
 387             meshRefiner,
 388             debug
 389         );
 390     }
 391
 392     if (wantSnap)
 393     {
 394         autoSnapDriver snapDriver
 395         (
 396             meshRefiner,
 397             globalToPatch
 398         );
 399
 400         // Snap parameters
 401         snapParameters snapParams(snapDict);
 402
 403         if (!overwrite)
 404         {
 405             const_cast<Time&>(mesh.time())++;
 406         }
 407
 408         snapDriver.doSnap(snapDict, motionDict, snapParams);
 409
 410         writeMesh
 411         (
 412             "Snapped mesh",
 413             meshRefiner,
 414             debug
 415         );
 416     }
 417
 418     if (wantLayers)
 419     {
 420         autoLayerDriver layerDriver(meshRefiner);
 421
 422         // Layer addition parameters
 423         layerParameters layerParams(layerDict, mesh.boundaryMesh());
 424
 425         if (!overwrite)
 426         {
 427             const_cast<Time&>(mesh.time())++;
 428         }
 429
 430         layerDriver.doLayers
 431         (
 432             layerDict,
 433             motionDict,
 434             layerParams,
 435             decomposer,
 436             distributor
 437         );
 438
 439         writeMesh
 440         (
 441             "Layer mesh",
 442             meshRefiner,
 443             debug
 444         );
 445     }
 446
 447
 448     Info<< "Finished meshing in = "
 449         << runTime.elapsedCpuTime() << " s." << endl;
 450
 451     Info<< "End\n" << endl;
 452
 453     return(0);
 454 }
 455
 456
 457 // ************************************************************************* //