FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

snappyHexMesh.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Application
00025     snappyHexMesh
00026 
00027 Description
00028     Automatic split hex mesher. Refines and snaps to surface.
00029 
00030 Usage
00031 
00032     - snappyHexMesh [OPTIONS]
00033 
00034     @param -overwrite
00035     Overwrite existing data.
00036 
00037     @param -case <dir>\n
00038     Case directory.
00039 
00040     @param -parallel \n
00041     Run in parallel.
00042 
00043     @param -help \n
00044     Display help message.
00045 
00046     @param -doc \n
00047     Display Doxygen API documentation page for this application.
00048 
00049     @param -srcDoc \n
00050     Display Doxygen source documentation page for this application.
00051 
00052 \*---------------------------------------------------------------------------*/
00053 
00054 #include <OpenFOAM/argList.H>
00055 #include <OpenFOAM/Time.H>
00056 #include <finiteVolume/fvMesh.H>
00057 #include <autoMesh/autoRefineDriver.H>
00058 #include <autoMesh/autoSnapDriver.H>
00059 #include <autoMesh/autoLayerDriver.H>
00060 #include <meshTools/searchableSurfaces.H>
00061 #include <autoMesh/refinementSurfaces.H>
00062 #include <autoMesh/shellSurfaces.H>
00063 #include <decompositionMethods/decompositionMethod.H>
00064 #include <dynamicMesh/fvMeshDistribute.H>
00065 #include <OpenFOAM/wallPolyPatch.H>
00066 #include <autoMesh/refinementParameters.H>
00067 #include <autoMesh/snapParameters.H>
00068 #include <autoMesh/layerParameters.H>
00069 
00070 
00071 using namespace Foam;
00072 
00073 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00074 
00075 // Check writing tolerance before doing any serious work
00076 scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
00077 {
00078     const boundBox& meshBb = mesh.bounds();
00079     scalar mergeDist = mergeTol * meshBb.mag();
00080     scalar writeTol = std::pow
00081     (
00082         scalar(10.0),
00083        -scalar(IOstream::defaultPrecision())
00084     );
00085 
00086     Info<< nl
00087         << "Overall mesh bounding box  : " << meshBb << nl
00088         << "Relative tolerance         : " << mergeTol << nl
00089         << "Absolute matching distance : " << mergeDist << nl
00090         << endl;
00091 
00092     if (mesh.time().writeFormat() == IOstream::ASCII && mergeTol < writeTol)
00093     {
00094         FatalErrorIn("getMergeDistance(const polyMesh&, const scalar)")
00095             << "Your current settings specify ASCII writing with "
00096             << IOstream::defaultPrecision() << " digits precision." << endl
00097             << "Your merging tolerance (" << mergeTol << ") is finer than this."
00098             << endl
00099             << "Please change your writeFormat to binary"
00100             << " or increase the writePrecision" << endl
00101             << "or adjust the merge tolerance (-mergeTol)."
00102             << exit(FatalError);
00103     }
00104 
00105     return mergeDist;
00106 }
00107 
00108 
00109 // Write mesh and additional information
00110 void writeMesh
00111 (
00112     const string& msg,
00113     const meshRefinement& meshRefiner,
00114     const label debug
00115 )
00116 {
00117     const fvMesh& mesh = meshRefiner.mesh();
00118 
00119     meshRefiner.printMeshInfo(debug, msg);
00120     Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;
00121 
00122     meshRefiner.write(meshRefinement::MESH|meshRefinement::SCALARLEVELS, "");
00123     if (debug & meshRefinement::OBJINTERSECTIONS)
00124     {
00125         meshRefiner.write
00126         (
00127             meshRefinement::OBJINTERSECTIONS,
00128             mesh.time().path()/meshRefiner.timeName()
00129         );
00130     }
00131     Info<< "Written mesh in = "
00132         << mesh.time().cpuTimeIncrement() << " s." << endl;
00133 }
00134 
00135 
00136 
00137 int main(int argc, char *argv[])
00138 {
00139     argList::validOptions.insert("overwrite", "");
00140 #   include <OpenFOAM/setRootCase.H>
00141 #   include <OpenFOAM/createTime.H>
00142     runTime.functionObjects().off();
00143 #   include <OpenFOAM/createMesh.H>
00144 
00145     Info<< "Read mesh in = "
00146         << runTime.cpuTimeIncrement() << " s" << endl;
00147 
00148     const bool overwrite = args.optionFound("overwrite");
00149 
00150 
00151     // Check patches and faceZones are synchronised
00152     mesh.boundaryMesh().checkParallelSync(true);
00153     meshRefinement::checkCoupledFaceZones(mesh);
00154 
00155 
00156     // Read decomposePar dictionary
00157     IOdictionary decomposeDict
00158     (
00159         IOobject
00160         (
00161             "decomposeParDict",
00162             runTime.system(),
00163             mesh,
00164             IOobject::MUST_READ,
00165             IOobject::NO_WRITE
00166         )
00167     );
00168 
00169     // Read meshing dictionary
00170     IOdictionary meshDict
00171     (
00172        IOobject
00173        (
00174             "snappyHexMeshDict",
00175             runTime.system(),
00176             mesh,
00177             IOobject::MUST_READ,
00178             IOobject::NO_WRITE
00179        )
00180     );
00181 
00182     // all surface geometry
00183     const dictionary& geometryDict = meshDict.subDict("geometry");
00184 
00185     // refinement parameters
00186     const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");
00187 
00188     // mesh motion and mesh quality parameters
00189     const dictionary& motionDict = meshDict.subDict("meshQualityControls");
00190 
00191     // snap-to-surface parameters
00192     const dictionary& snapDict = meshDict.subDict("snapControls");
00193 
00194     // layer addition parameters
00195     const dictionary& layerDict = meshDict.subDict("addLayersControls");
00196 
00197 
00198     const scalar mergeDist = getMergeDistance
00199     (
00200         mesh,
00201         readScalar(meshDict.lookup("mergeTolerance"))
00202     );
00203 
00204 
00205 
00206     // Debug
00207     // ~~~~~
00208 
00209     const label debug(readLabel(meshDict.lookup("debug")));
00210     if (debug > 0)
00211     {
00212         meshRefinement::debug = debug;
00213         autoRefineDriver::debug = debug;
00214         autoSnapDriver::debug = debug;
00215         autoLayerDriver::debug = debug;
00216     }
00217 
00218 
00219     // Read geometry
00220     // ~~~~~~~~~~~~~
00221 
00222     searchableSurfaces allGeometry
00223     (
00224         IOobject
00225         (
00226             "abc",                      // dummy name
00227             mesh.time().constant(),     // instance
00228             //mesh.time().findInstance("triSurface", word::null),// instance
00229             "triSurface",               // local
00230             mesh.time(),                // registry
00231             IOobject::MUST_READ,
00232             IOobject::NO_WRITE
00233         ),
00234         geometryDict
00235     );
00236 
00237 
00238     // Read refinement surfaces
00239     // ~~~~~~~~~~~~~~~~~~~~~~~~
00240 
00241     Info<< "Reading refinement surfaces." << endl;
00242     refinementSurfaces surfaces
00243     (
00244         allGeometry,
00245         refineDict.subDict("refinementSurfaces")
00246     );
00247     Info<< "Read refinement surfaces in = "
00248         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
00249 
00250 
00251     // Read refinement shells
00252     // ~~~~~~~~~~~~~~~~~~~~~~
00253 
00254     Info<< "Reading refinement shells." << endl;
00255     shellSurfaces shells
00256     (
00257         allGeometry,
00258         refineDict.subDict("refinementRegions")
00259     );
00260     Info<< "Read refinement shells in = "
00261         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
00262 
00263 
00264     Info<< "Setting refinement level of surface to be consistent"
00265         << " with shells." << endl;
00266     surfaces.setMinLevelFields(shells);
00267     Info<< "Checked shell refinement in = "
00268         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
00269 
00270 
00271     // Refinement engine
00272     // ~~~~~~~~~~~~~~~~~
00273 
00274     Info<< nl
00275         << "Determining initial surface intersections" << nl
00276         << "-----------------------------------------" << nl
00277         << endl;
00278 
00279     // Main refinement engine
00280     meshRefinement meshRefiner
00281     (
00282         mesh,
00283         mergeDist,          // tolerance used in sorting coordinates
00284         overwrite,          // overwrite mesh files?
00285         surfaces,           // for surface intersection refinement
00286         shells              // for volume (inside/outside) refinement
00287     );
00288     Info<< "Calculated surface intersections in = "
00289         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
00290 
00291     // Some stats
00292     meshRefiner.printMeshInfo(debug, "Initial mesh");
00293 
00294     meshRefiner.write
00295     (
00296         debug&meshRefinement::OBJINTERSECTIONS,
00297         mesh.time().path()/meshRefiner.timeName()
00298     );
00299 
00300 
00301     // Add all the surface regions as patches
00302     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00303 
00304     labelList globalToPatch;
00305     {
00306         Info<< nl
00307             << "Adding patches for surface regions" << nl
00308             << "----------------------------------" << nl
00309             << endl;
00310 
00311         // From global region number to mesh patch.
00312         globalToPatch.setSize(surfaces.nRegions(), -1);
00313 
00314         Info<< "Patch\tRegion" << nl
00315             << "-----\t------"
00316             << endl;
00317 
00318         const labelList& surfaceGeometry = surfaces.surfaces();
00319         forAll(surfaceGeometry, surfI)
00320         {
00321             label geomI = surfaceGeometry[surfI];
00322 
00323             const wordList& regNames = allGeometry.regionNames()[geomI];
00324 
00325             Info<< surfaces.names()[surfI] << ':' << nl << nl;
00326 
00327             forAll(regNames, i)
00328             {
00329                 label patchI = meshRefiner.addMeshedPatch
00330                 (
00331                     regNames[i],
00332                     wallPolyPatch::typeName
00333                 );
00334 
00335                 Info<< patchI << '\t' << regNames[i] << nl;
00336 
00337                 globalToPatch[surfaces.globalRegion(surfI, i)] = patchI;
00338             }
00339 
00340             Info<< nl;
00341         }
00342         Info<< "Added patches in = "
00343             << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
00344     }
00345 
00346 
00347     // Parallel
00348     // ~~~~~~~~
00349 
00350     // Decomposition
00351     autoPtr<decompositionMethod> decomposerPtr
00352     (
00353         decompositionMethod::New
00354         (
00355             decomposeDict,
00356             mesh
00357         )
00358     );
00359     decompositionMethod& decomposer = decomposerPtr();
00360 
00361     if (Pstream::parRun() && !decomposer.parallelAware())
00362     {
00363         FatalErrorIn(args.executable())
00364             << "You have selected decomposition method "
00365             << decomposer.typeName
00366             << " which is not parallel aware." << endl
00367             << "Please select one that is (hierarchical, parMetis)"
00368             << exit(FatalError);
00369     }
00370 
00371     // Mesh distribution engine (uses tolerance to reconstruct meshes)
00372     fvMeshDistribute distributor(mesh, mergeDist);
00373 
00374 
00375 
00376 
00377 
00378     // Now do the real work -refinement -snapping -layers
00379     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00380 
00381     Switch wantRefine(meshDict.lookup("castellatedMesh"));
00382     Switch wantSnap(meshDict.lookup("snap"));
00383     Switch wantLayers(meshDict.lookup("addLayers"));
00384 
00385     if (wantRefine)
00386     {
00387         cpuTime timer;
00388 
00389         autoRefineDriver refineDriver
00390         (
00391             meshRefiner,
00392             decomposer,
00393             distributor,
00394             globalToPatch
00395         );
00396 
00397         // Refinement parameters
00398         refinementParameters refineParams(refineDict);
00399 
00400         if (!overwrite)
00401         {
00402             const_cast<Time&>(mesh.time())++;
00403         }
00404 
00405         refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
00406 
00407         writeMesh
00408         (
00409             "Refined mesh",
00410             meshRefiner,
00411             debug
00412         );
00413 
00414         Info<< "Mesh refined in = "
00415             << timer.cpuTimeIncrement() << " s." << endl;
00416     }
00417 
00418     if (wantSnap)
00419     {
00420         cpuTime timer;
00421 
00422         autoSnapDriver snapDriver
00423         (
00424             meshRefiner,
00425             globalToPatch
00426         );
00427 
00428         // Snap parameters
00429         snapParameters snapParams(snapDict);
00430 
00431         if (!overwrite)
00432         {
00433             const_cast<Time&>(mesh.time())++;
00434         }
00435 
00436         snapDriver.doSnap(snapDict, motionDict, snapParams);
00437 
00438         writeMesh
00439         (
00440             "Snapped mesh",
00441             meshRefiner,
00442             debug
00443         );
00444 
00445         Info<< "Mesh snapped in = "
00446             << timer.cpuTimeIncrement() << " s." << endl;
00447     }
00448 
00449     if (wantLayers)
00450     {
00451         cpuTime timer;
00452 
00453         autoLayerDriver layerDriver(meshRefiner);
00454 
00455         // Layer addition parameters
00456         layerParameters layerParams(layerDict, mesh.boundaryMesh());
00457 
00459         bool preBalance;
00460         {
00461             refinementParameters refineParams(refineDict);
00462 
00463             preBalance = returnReduce
00464             (
00465                 (mesh.nCells() >= refineParams.maxLocalCells()),
00466                 orOp<bool>()
00467             );
00468         }
00469 
00470 
00471         if (!overwrite)
00472         {
00473             const_cast<Time&>(mesh.time())++;
00474         }
00475 
00476         layerDriver.doLayers
00477         (
00478             layerDict,
00479             motionDict,
00480             layerParams,
00481             preBalance,
00482             decomposer,
00483             distributor
00484         );
00485 
00486         writeMesh
00487         (
00488             "Layer mesh",
00489             meshRefiner,
00490             debug
00491         );
00492 
00493         Info<< "Layers added in = "
00494             << timer.cpuTimeIncrement() << " s." << endl;
00495     }
00496 
00497 
00498     Info<< "Finished meshing in = "
00499         << runTime.elapsedCpuTime() << " s." << endl;
00500 
00501     Info<< "End\n" << endl;
00502 
00503     return(0);
00504 }
00505 
00506 
00507 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines