/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    snappyHexMesh

Description
    Automatic split hex mesher. Refines and snaps to surface.

\*---------------------------------------------------------------------------*/

#include "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "autoRefineDriver.H"
#include "autoSnapDriver.H"
#include "autoLayerDriver.H"
#include "searchableSurfaces.H"
#include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "shellSurfaces.H"
#include "decompositionMethod.H"
#include "noDecomp.H"
#include "fvMeshDistribute.H"
#include "wallPolyPatch.H"
#include "refinementParameters.H"
#include "snapParameters.H"
#include "layerParameters.H"
#include "vtkSetWriter.H"
#include "faceSet.H"
#include "motionSmoother.H"
#include "polyTopoChange.H"
#include "cellModeller.H"
#include "uindirectPrimitivePatch.H"
#include "surfZoneIdentifierList.H"
#include "UnsortedMeshedSurface.H"
#include "MeshedSurface.H"
#include "globalIndex.H"
#include "IOmanip.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Convert size (as fraction of defaultCellSize) to refinement level
label sizeCoeffToRefinement
(
    const scalar level0Coeff,   // ratio of hex cell size v.s. defaultCellSize
    const scalar sizeCoeff
)
{
     return round(::log(level0Coeff/sizeCoeff)/::log(2));
}


autoPtr<refinementSurfaces> createRefinementSurfaces
(
    const searchableSurfaces& allGeometry,
    const dictionary& surfacesDict,
    const dictionary& shapeControlDict,
    const label gapLevelIncrement,
    const scalar level0Coeff
)
{
    autoPtr<refinementSurfaces> surfacePtr;

    // Count number of surfaces.
    label surfI = 0;
    forAll(allGeometry.names(), geomI)
    {
        const word& geomName = allGeometry.names()[geomI];

        if (surfacesDict.found(geomName))
        {
            surfI++;
        }
    }

    labelList surfaces(surfI);
    wordList names(surfI);
    wordList faceZoneNames(surfI);
    wordList cellZoneNames(surfI);
    List<refinementSurfaces::areaSelectionAlgo> zoneInside
    (
        surfI,
        refinementSurfaces::NONE
    );
    pointField zoneInsidePoints(surfI);
    List<refinementSurfaces::faceZoneType> faceType
    (
        surfI,
        refinementSurfaces::INTERNAL
    );
    labelList regionOffset(surfI);

    labelList globalMinLevel(surfI, 0);
    labelList globalMaxLevel(surfI, 0);
    labelList globalLevelIncr(surfI, 0);
    PtrList<dictionary> globalPatchInfo(surfI);
    List<Map<label> > regionMinLevel(surfI);
    List<Map<label> > regionMaxLevel(surfI);
    List<Map<label> > regionLevelIncr(surfI);
    List<Map<scalar> > regionAngle(surfI);
    List<Map<autoPtr<dictionary> > > regionPatchInfo(surfI);

    surfI = 0;
    forAll(allGeometry.names(), geomI)
    {
        const word& geomName = allGeometry.names()[geomI];


        // Definition of surfaces to conform to
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        if (surfacesDict.found(geomName))
        {
            names[surfI] = geomName;
            surfaces[surfI] = geomI;

            const dictionary& shapeDict = shapeControlDict.subDict(geomName);

            // Find the index in shapeControlDict
            // Invert surfaceCellSize to get the refinementLevel

            const word scsFuncName =
                shapeDict.lookup("surfaceCellSizeFunction");
            const dictionary& scsDict =
                shapeDict.subDict(scsFuncName + "Coeffs");

            const scalar surfaceCellSize =
                readScalar(scsDict.lookup("surfaceCellSizeCoeff"));

            const label refLevel = sizeCoeffToRefinement
            (
                level0Coeff,
                surfaceCellSize
            );

            globalMinLevel[surfI] = refLevel;
            globalMaxLevel[surfI] = refLevel;
            globalLevelIncr[surfI] = gapLevelIncrement;

            const dictionary& dict = surfacesDict.subDict(geomName);

            // Global zone names per surface
            if (dict.readIfPresent("faceZone", faceZoneNames[surfI]))
            {
                // Read optional entry to determine inside of faceZone

                word method;
                bool hasSide = dict.readIfPresent("cellZoneInside", method);
                if (hasSide)
                {
                    zoneInside[surfI] =
                        refinementSurfaces::areaSelectionAlgoNames[method];
                    if (zoneInside[surfI] == refinementSurfaces::INSIDEPOINT)
                    {
                        dict.lookup("insidePoint") >> zoneInsidePoints[surfI];
                    }

                }
                else
                {
                    // Check old syntax
                    bool inside;
                    if (dict.readIfPresent("zoneInside", inside))
                    {
                        hasSide = true;
                        zoneInside[surfI] =
                            (
                                inside
                              ? refinementSurfaces::INSIDE
                              : refinementSurfaces::OUTSIDE
                            );
                    }
                }

                // Read optional cellZone name

                if (dict.readIfPresent("cellZone", cellZoneNames[surfI]))
                {
                    if
                    (
                        (
                            zoneInside[surfI] == refinementSurfaces::INSIDE
                         || zoneInside[surfI] == refinementSurfaces::OUTSIDE
                        )
                    && !allGeometry[surfaces[surfI]].hasVolumeType()
                    )
                    {
                        IOWarningIn
                        (
                            "createRefinementSurfaces(..)",
                            dict
                        )   << "Illegal entry zoneInside "
                            << refinementSurfaces::areaSelectionAlgoNames
                               [
                                   zoneInside[surfI]
                               ]
                            << " for faceZone "
                            << faceZoneNames[surfI]
                            << " since surface " << names[surfI]
                            << " is not closed." << endl;
                    }
                }
                else if (hasSide)
                {
                    IOWarningIn
                    (
                        "createRefinementSurfaces(..)",
                        dict
                    )   << "Unused entry zoneInside for faceZone "
                        << faceZoneNames[surfI]
                        << " since no cellZone specified."
                        << endl;
                }

                // How to handle faces on faceZone
                word faceTypeMethod;
                if (dict.readIfPresent("faceType", faceTypeMethod))
                {
                    faceType[surfI] =
                        refinementSurfaces::faceZoneTypeNames[faceTypeMethod];
                }
            }


            // Global perpendicular angle
            if (dict.found("patchInfo"))
            {
                globalPatchInfo.set
                (
                    surfI,
                    dict.subDict("patchInfo").clone()
                );
            }


            // Per region override of patchInfo

            if (dict.found("regions"))
            {
                const dictionary& regionsDict = dict.subDict("regions");
                const wordList& regionNames =
                    allGeometry[surfaces[surfI]].regions();

                forAll(regionNames, regionI)
                {
                    if (regionsDict.found(regionNames[regionI]))
                    {
                        // Get the dictionary for region
                        const dictionary& regionDict = regionsDict.subDict
                        (
                            regionNames[regionI]
                        );

                        if (regionDict.found("patchInfo"))
                        {
                            regionPatchInfo[surfI].insert
                            (
                                regionI,
                                regionDict.subDict("patchInfo").clone()
                            );
                        }
                    }
                }
            }

            // Per region override of cellSize
            if (shapeDict.found("regions"))
            {
                const dictionary& shapeControlRegionsDict =
                    shapeDict.subDict("regions");
                const wordList& regionNames =
                    allGeometry[surfaces[surfI]].regions();

                forAll(regionNames, regionI)
                {
                    if (shapeControlRegionsDict.found(regionNames[regionI]))
                    {
                        const dictionary& shapeControlRegionDict =
                            shapeControlRegionsDict.subDict
                            (
                                regionNames[regionI]
                            );

                        const word scsFuncName =
                            shapeControlRegionDict.lookup
                            (
                                "surfaceCellSizeFunction"
                            );
                        const dictionary& scsDict =
                            shapeControlRegionDict.subDict
                            (
                                scsFuncName + "Coeffs"
                            );

                        const scalar surfaceCellSize =
                            readScalar
                            (
                                scsDict.lookup("surfaceCellSizeCoeff")
                            );

                        const label refLevel = sizeCoeffToRefinement
                        (
                            level0Coeff,
                            surfaceCellSize
                        );

                        regionMinLevel[surfI].insert(regionI, refLevel);
                        regionMaxLevel[surfI].insert(regionI, refLevel);
                        regionLevelIncr[surfI].insert(regionI, 0);
                    }
                }
            }
            surfI++;
        }
    }

    // Calculate local to global region offset
    label nRegions = 0;

    forAll(surfaces, surfI)
    {
        regionOffset[surfI] = nRegions;
        nRegions += allGeometry[surfaces[surfI]].regions().size();
    }

    // Rework surface specific information into information per global region
    labelList minLevel(nRegions, 0);
    labelList maxLevel(nRegions, 0);
    labelList gapLevel(nRegions, -1);
    PtrList<dictionary> patchInfo(nRegions);

    forAll(globalMinLevel, surfI)
    {
        label nRegions = allGeometry[surfaces[surfI]].regions().size();

        // Initialise to global (i.e. per surface)
        for (label i = 0; i < nRegions; i++)
        {
            label globalRegionI = regionOffset[surfI] + i;
            minLevel[globalRegionI] = globalMinLevel[surfI];
            maxLevel[globalRegionI] = globalMaxLevel[surfI];
            gapLevel[globalRegionI] =
                maxLevel[globalRegionI]
              + globalLevelIncr[surfI];

            if (globalPatchInfo.set(surfI))
            {
                patchInfo.set
                (
                    globalRegionI,
                    globalPatchInfo[surfI].clone()
                );
            }
        }

        // Overwrite with region specific information
        forAllConstIter(Map<label>, regionMinLevel[surfI], iter)
        {
            label globalRegionI = regionOffset[surfI] + iter.key();

            minLevel[globalRegionI] = iter();
            maxLevel[globalRegionI] = regionMaxLevel[surfI][iter.key()];
            gapLevel[globalRegionI] =
                maxLevel[globalRegionI]
              + regionLevelIncr[surfI][iter.key()];
        }

        const Map<autoPtr<dictionary> >& localInfo = regionPatchInfo[surfI];
        forAllConstIter(Map<autoPtr<dictionary> >, localInfo, iter)
        {
            label globalRegionI = regionOffset[surfI] + iter.key();
            patchInfo.set(globalRegionI, iter()().clone());
        }
    }

    surfacePtr.set
    (
        new refinementSurfaces
        (
            allGeometry,
            surfaces,
            names,
            faceZoneNames,
            cellZoneNames,
            zoneInside,
            zoneInsidePoints,
            faceType,
            regionOffset,
            minLevel,
            maxLevel,
            gapLevel,
            scalarField(nRegions, -GREAT),  //perpendicularAngle,
            patchInfo
        )
    );


    const refinementSurfaces& rf = surfacePtr();

    Info<< setw(20) << "Region"
        << setw(10) << "Min Level"
        << setw(10) << "Max Level"
        << setw(10) << "Gap Level" << nl
        << setw(20) << "------"
        << setw(10) << "---------"
        << setw(10) << "---------"
        << setw(10) << "---------" << endl;

    forAll(rf.surfaces(), surfI)
    {
        label geomI = rf.surfaces()[surfI];

        Info<< rf.names()[surfI] << ':' << nl;

        const wordList& regionNames = allGeometry.regionNames()[geomI];

        forAll(regionNames, regionI)
        {
            label globalI = rf.globalRegion(surfI, regionI);

            Info<< setw(20) << regionNames[regionI]
                << setw(10) << rf.minLevel()[globalI]
                << setw(10) << rf.maxLevel()[globalI]
                << setw(10) << rf.gapLevel()[globalI] << endl;
        }
    }


    return surfacePtr;
}


void extractSurface
(
    const polyMesh& mesh,
    const Time& runTime,
    const labelHashSet& includePatches,
    const fileName& outFileName
)
{
    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();

    // Collect sizes. Hash on names to handle local-only patches (e.g.
    //  processor patches)
    HashTable<label> patchSize(1000);
    label nFaces = 0;
    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const polyPatch& pp = bMesh[iter.key()];
        patchSize.insert(pp.name(), pp.size());
        nFaces += pp.size();
    }
    Pstream::mapCombineGather(patchSize, plusEqOp<label>());


    // Allocate zone/patch for all patches
    HashTable<label> compactZoneID(1000);
    forAllConstIter(HashTable<label>, patchSize, iter)
    {
        label sz = compactZoneID.size();
        compactZoneID.insert(iter.key(), sz);
    }
    Pstream::mapCombineScatter(compactZoneID);


    // Rework HashTable into labelList just for speed of conversion
    labelList patchToCompactZone(bMesh.size(), -1);
    forAllConstIter(HashTable<label>, compactZoneID, iter)
    {
        label patchI = bMesh.findPatchID(iter.key());
        if (patchI != -1)
        {
            patchToCompactZone[patchI] = iter();
        }
    }

    // Collect faces on zones
    DynamicList<label> faceLabels(nFaces);
    DynamicList<label> compactZones(nFaces);
    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const polyPatch& pp = bMesh[iter.key()];
        forAll(pp, i)
        {
            faceLabels.append(pp.start()+i);
            compactZones.append(patchToCompactZone[pp.index()]);
        }
    }

    // Addressing engine for all faces
    uindirectPrimitivePatch allBoundary
    (
        UIndirectList<face>(mesh.faces(), faceLabels),
        mesh.points()
    );


    // Find correspondence to master points
    labelList pointToGlobal;
    labelList uniqueMeshPoints;
    autoPtr<globalIndex> globalNumbers = mesh.globalData().mergePoints
    (
        allBoundary.meshPoints(),
        allBoundary.meshPointMap(),
        pointToGlobal,
        uniqueMeshPoints
    );

    // Gather all unique points on master
    List<pointField> gatheredPoints(Pstream::nProcs());
    gatheredPoints[Pstream::myProcNo()] = pointField
    (
        mesh.points(),
        uniqueMeshPoints
    );
    Pstream::gatherList(gatheredPoints);

    // Gather all faces
    List<faceList> gatheredFaces(Pstream::nProcs());
    gatheredFaces[Pstream::myProcNo()] = allBoundary.localFaces();
    forAll(gatheredFaces[Pstream::myProcNo()], i)
    {
        inplaceRenumber(pointToGlobal, gatheredFaces[Pstream::myProcNo()][i]);
    }
    Pstream::gatherList(gatheredFaces);

    // Gather all ZoneIDs
    List<labelList> gatheredZones(Pstream::nProcs());
    gatheredZones[Pstream::myProcNo()] = compactZones.xfer();
    Pstream::gatherList(gatheredZones);

    // On master combine all points, faces, zones
    if (Pstream::master())
    {
        pointField allPoints = ListListOps::combine<pointField>
        (
            gatheredPoints,
            accessOp<pointField>()
        );
        gatheredPoints.clear();

        faceList allFaces = ListListOps::combine<faceList>
        (
            gatheredFaces,
            accessOp<faceList>()
        );
        gatheredFaces.clear();

        labelList allZones = ListListOps::combine<labelList>
        (
            gatheredZones,
            accessOp<labelList>()
        );
        gatheredZones.clear();


        // Zones
        surfZoneIdentifierList surfZones(compactZoneID.size());
        forAllConstIter(HashTable<label>, compactZoneID, iter)
        {
            surfZones[iter()] = surfZoneIdentifier(iter.key(), iter());
            Info<< "surfZone " << iter()  <<  " : " << surfZones[iter()].name()
                << endl;
        }

        UnsortedMeshedSurface<face> unsortedFace
        (
            xferMove(allPoints),
            xferMove(allFaces),
            xferMove(allZones),
            xferMove(surfZones)
        );


        MeshedSurface<face> sortedFace(unsortedFace);

        fileName globalCasePath
        (
            runTime.processorCase()
          ? runTime.path()/".."/outFileName
          : runTime.path()/outFileName
        );

        Info<< "Writing merged surface to " << globalCasePath << endl;

        sortedFace.write(globalCasePath);
    }
}


// Check writing tolerance before doing any serious work
scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
{
    const boundBox& meshBb = mesh.bounds();
    scalar mergeDist = mergeTol * meshBb.mag();

    Info<< nl
        << "Overall mesh bounding box  : " << meshBb << nl
        << "Relative tolerance         : " << mergeTol << nl
        << "Absolute matching distance : " << mergeDist << nl
        << endl;

    // check writing tolerance
    if (mesh.time().writeFormat() == IOstream::ASCII)
    {
        const scalar writeTol = std::pow
        (
            scalar(10.0),
            -scalar(IOstream::defaultPrecision())
        );

        if (mergeTol < writeTol)
        {
            FatalErrorIn("getMergeDistance(const polyMesh&, const dictionary&)")
                << "Your current settings specify ASCII writing with "
                << IOstream::defaultPrecision() << " digits precision." << nl
                << "Your merging tolerance (" << mergeTol
                << ") is finer than this." << nl
                << "Change to binary writeFormat, "
                << "or increase the writePrecision" << endl
                << "or adjust the merge tolerance (mergeTol)."
                << exit(FatalError);
        }
    }

    return mergeDist;
}


// Write mesh and additional information
void writeMesh
(
    const string& msg,
    const meshRefinement& meshRefiner,
    const bool writeLevel,
    const label debug
)
{
    const fvMesh& mesh = meshRefiner.mesh();

    meshRefiner.printMeshInfo(debug, msg);
    Info<< "Writing mesh to time " << meshRefiner.timeName() << endl;

    label flag = meshRefinement::MESH;
    if (writeLevel)
    {
        flag |= meshRefinement::SCALARLEVELS;
    }
    if (debug & meshRefinement::OBJINTERSECTIONS)
    {
        flag |= meshRefinement::OBJINTERSECTIONS;
    }
    meshRefiner.write(flag, mesh.time().path()/meshRefiner.timeName());
    Info<< "Wrote mesh in = "
        << mesh.time().cpuTimeIncrement() << " s." << endl;
}


int main(int argc, char *argv[])
{
#   include "addOverwriteOption.H"
    Foam::argList::addBoolOption
    (
        "checkGeometry",
        "check all surface geometry for quality"
    );
    Foam::argList::addOption
    (
        "surfaceSimplify",
        "boundBox",
        "simplify the surface using snappyHexMesh starting from a boundBox"
    );
    Foam::argList::addBoolOption
    (
        "writeLevel",
        "write pointLevel and cellLevel postprocessing files"
    );
    Foam::argList::addOption
    (
        "patches",
        "(patch0 .. patchN)",
        "only triangulate selected patches (wildcards supported)"
    );
    Foam::argList::addOption
    (
        "outFile",
        "fileName",
        "name of the file to save the simplified surface to"
    );
#   include "addDictOption.H"

#   include "setRootCase.H"
#   include "createTime.H"
    runTime.functionObjects().off();

    const bool overwrite = args.optionFound("overwrite");
    const bool checkGeometry = args.optionFound("checkGeometry");
    const bool surfaceSimplify = args.optionFound("surfaceSimplify");
    const bool writeLevel = args.optionFound("writeLevel");

    autoPtr<fvMesh> meshPtr;

    if (surfaceSimplify)
    {
        IOdictionary foamyHexMeshDict
        (
           IOobject
           (
                "foamyHexMeshDict",
                runTime.system(),
                runTime,
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE
           )
        );

        const dictionary& motionDict =
            foamyHexMeshDict.subDict("motionControl");

        const scalar defaultCellSize =
            readScalar(motionDict.lookup("defaultCellSize"));

        Info<< "Constructing single cell mesh from boundBox" << nl << endl;

        boundBox bb(args.optionRead<boundBox>("surfaceSimplify"));

        labelList owner(6, label(0));
        labelList neighbour(0);

        const cellModel& hexa = *(cellModeller::lookup("hex"));
        faceList faces = hexa.modelFaces();

        meshPtr.set
        (
            new fvMesh
            (
                IOobject
                (
                    fvMesh::defaultRegion,
                    runTime.timeName(),
                    runTime,
                    IOobject::NO_READ
                ),
                xferMove<Field<vector> >(bb.points()()),
                faces.xfer(),
                owner.xfer(),
                neighbour.xfer()
            )
        );

        List<polyPatch*> patches(1);

        patches[0] = new wallPolyPatch
        (
            "boundary",
            6,
            0,
            0,
            meshPtr().boundaryMesh(),
            emptyPolyPatch::typeName
        );

        meshPtr().addFvPatches(patches);

        const scalar initialCellSize = ::pow(meshPtr().V()[0], 1.0/3.0);
        const label initialRefLevels =
            ::log(initialCellSize/defaultCellSize)/::log(2);

        Info<< "Default cell size = " << defaultCellSize << endl;
        Info<< "Initial cell size = " << initialCellSize << endl;

        Info<< "Initial refinement levels = " << initialRefLevels << endl;

        Info<< "Mesh starting size = " << meshPtr().nCells() << endl;

        // meshCutter must be destroyed before writing the mesh otherwise it
        // writes the cellLevel/pointLevel files
        {
            hexRef8 meshCutter(meshPtr(), false);

            for (label refineI = 0; refineI < initialRefLevels; ++refineI)
            {
                // Mesh changing engine.
                polyTopoChange meshMod(meshPtr(), true);

                // Play refinement commands into mesh changer.
                meshCutter.setRefinement(identity(meshPtr().nCells()), meshMod);

                // Create mesh (no inflation), return map from old to new mesh.
                autoPtr<mapPolyMesh> map = meshMod.changeMesh(meshPtr(), false);

                // Update fields
                meshPtr().updateMesh(map);

                // Delete mesh volumes.
                meshPtr().clearOut();

                Info<< "Refinement Iteration " << refineI + 1
                    << ", Mesh size = " << meshPtr().nCells() << endl;
            }
        }

        Info<< "Mesh end size = " << meshPtr().nCells() << endl;

        meshPtr().write();
    }
    else
    {
        Foam::Info
            << "Create mesh for time = "
            << runTime.timeName() << Foam::nl << Foam::endl;

        meshPtr.set
        (
            new fvMesh
            (
                Foam::IOobject
                (
                    Foam::fvMesh::defaultRegion,
                    runTime.timeName(),
                    runTime,
                    Foam::IOobject::MUST_READ
                )
            )
        );
    }

    fvMesh& mesh = meshPtr();

    Info<< "Read mesh in = "
        << runTime.cpuTimeIncrement() << " s" << endl;

    // Check patches and faceZones are synchronised
    mesh.boundaryMesh().checkParallelSync(true);
    meshRefinement::checkCoupledFaceZones(mesh);


    // Read meshing dictionary
    const word dictName("snappyHexMeshDict");
    #include "setSystemMeshDictionaryIO.H"
    const IOdictionary meshDict(dictIO);


    // all surface geometry
    const dictionary& geometryDict = meshDict.subDict("geometry");

    // refinement parameters
    const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");

    // mesh motion and mesh quality parameters
    const dictionary& motionDict = meshDict.subDict("meshQualityControls");

    // snap-to-surface parameters
    const dictionary& snapDict = meshDict.subDict("snapControls");

    // layer addition parameters
    const dictionary& layerDict = meshDict.subDict("addLayersControls");

    // absolute merge distance
    const scalar mergeDist = getMergeDistance
    (
        mesh,
        readScalar(meshDict.lookup("mergeTolerance"))
    );



    // Read decomposePar dictionary
    dictionary decomposeDict;
    {
        if (Pstream::parRun())
        {
            decomposeDict = IOdictionary
            (
                IOobject
                (
                    "decomposeParDict",
                    runTime.system(),
                    mesh,
                    IOobject::MUST_READ_IF_MODIFIED,
                    IOobject::NO_WRITE
                )
            );
        }
        else
        {
            decomposeDict.add("method", "none");
            decomposeDict.add("numberOfSubdomains", 1);
        }
    }


    // Debug
    // ~~~~~

    const label debug = meshDict.lookupOrDefault<label>("debug", 0);
    if (debug > 0)
    {
        meshRefinement::debug   = debug;
        autoRefineDriver::debug = debug;
        autoSnapDriver::debug   = debug;
        autoLayerDriver::debug  = debug;
    }


    // Read geometry
    // ~~~~~~~~~~~~~

    searchableSurfaces allGeometry
    (
        IOobject
        (
            "abc",                      // dummy name
            mesh.time().constant(),     // instance
            //mesh.time().findInstance("triSurface", word::null),// instance
            "triSurface",               // local
            mesh.time(),                // registry
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
        geometryDict
    );


    // Read refinement surfaces
    // ~~~~~~~~~~~~~~~~~~~~~~~~

    autoPtr<refinementSurfaces> surfacesPtr;

    Info<< "Reading refinement surfaces." << endl;

    if (surfaceSimplify)
    {
        IOdictionary foamyHexMeshDict
        (
           IOobject
           (
                "foamyHexMeshDict",
                runTime.system(),
                runTime,
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE
           )
        );

        const dictionary& conformationDict =
            foamyHexMeshDict.subDict("surfaceConformation").subDict
            (
                "geometryToConformTo"
            );

        const dictionary& motionDict =
            foamyHexMeshDict.subDict("motionControl");

        const dictionary& shapeControlDict =
            motionDict.subDict("shapeControlFunctions");

        // Calculate current ratio of hex cells v.s. wanted cell size
        const scalar defaultCellSize =
            readScalar(motionDict.lookup("defaultCellSize"));

        const scalar initialCellSize = ::pow(meshPtr().V()[0], 1.0/3.0);

        //Info<< "Wanted cell size  = " << defaultCellSize << endl;
        //Info<< "Current cell size = " << initialCellSize << endl;
        //Info<< "Fraction          = " << initialCellSize/defaultCellSize
        //    << endl;

        surfacesPtr =
            createRefinementSurfaces
            (
                allGeometry,
                conformationDict,
                shapeControlDict,
                refineDict.lookupOrDefault("gapLevelIncrement", 0),
                initialCellSize/defaultCellSize
            );
    }
    else
    {
        surfacesPtr.set
        (
            new refinementSurfaces
            (
                allGeometry,
                refineDict.subDict("refinementSurfaces"),
                refineDict.lookupOrDefault("gapLevelIncrement", 0)
            )
        );

        Info<< "Read refinement surfaces in = "
            << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
    }

    refinementSurfaces& surfaces = surfacesPtr();


    // Checking only?

    if (checkGeometry)
    {
        // Extract patchInfo
        List<wordList> patchTypes(allGeometry.size());

        const PtrList<dictionary>& patchInfo = surfaces.patchInfo();
        const labelList& surfaceGeometry = surfaces.surfaces();
        forAll(surfaceGeometry, surfI)
        {
            label geomI = surfaceGeometry[surfI];
            const wordList& regNames = allGeometry.regionNames()[geomI];

            patchTypes[geomI].setSize(regNames.size());
            forAll(regNames, regionI)
            {
                label globalRegionI = surfaces.globalRegion(surfI, regionI);

                if (patchInfo.set(globalRegionI))
                {
                    patchTypes[geomI][regionI] =
                        word(patchInfo[globalRegionI].lookup("type"));
                }
                else
                {
                    patchTypes[geomI][regionI] = wallPolyPatch::typeName;
                }
            }
        }

        // Write some stats
        allGeometry.writeStats(patchTypes, Info);
        // Check topology
        allGeometry.checkTopology(true);
        // Check geometry
        allGeometry.checkGeometry
        (
            100.0,      // max size ratio
            1e-9,       // intersection tolerance
            autoPtr<writer<scalar> >(new vtkSetWriter<scalar>()),
            0.01,       // min triangle quality
            true
        );

        return 0;
    }



    // Read refinement shells
    // ~~~~~~~~~~~~~~~~~~~~~~

    Info<< "Reading refinement shells." << endl;
    shellSurfaces shells
    (
        allGeometry,
        refineDict.subDict("refinementRegions")
    );
    Info<< "Read refinement shells in = "
        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;


    Info<< "Setting refinement level of surface to be consistent"
        << " with shells." << endl;
    surfaces.setMinLevelFields(shells);
    Info<< "Checked shell refinement in = "
        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;


    // Read feature meshes
    // ~~~~~~~~~~~~~~~~~~~

    Info<< "Reading features." << endl;
    refinementFeatures features
    (
        mesh,
        refineDict.lookup("features")
    );
    Info<< "Read features in = "
        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;



    // Refinement engine
    // ~~~~~~~~~~~~~~~~~

    Info<< nl
        << "Determining initial surface intersections" << nl
        << "-----------------------------------------" << nl
        << endl;

    // Main refinement engine
    meshRefinement meshRefiner
    (
        mesh,
        mergeDist,          // tolerance used in sorting coordinates
        overwrite,          // overwrite mesh files?
        surfaces,           // for surface intersection refinement
        features,           // for feature edges/point based refinement
        shells              // for volume (inside/outside) refinement
    );
    Info<< "Calculated surface intersections in = "
        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;

    // Some stats
    meshRefiner.printMeshInfo(debug, "Initial mesh");

    meshRefiner.write
    (
        debug & meshRefinement::OBJINTERSECTIONS,
        mesh.time().path()/meshRefiner.timeName()
    );


    // Add all the surface regions as patches
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    //- Global surface region to patch (non faceZone surface) or patches
    //  (faceZone surfaces)
    labelList globalToMasterPatch;
    labelList globalToSlavePatch;
    {
        Info<< nl
            << "Adding patches for surface regions" << nl
            << "----------------------------------" << nl
            << endl;

        // From global region number to mesh patch.
        globalToMasterPatch.setSize(surfaces.nRegions(), -1);
        globalToSlavePatch.setSize(surfaces.nRegions(), -1);

        Info<< setf(ios_base::left)
            << setw(6) << "Patch"
            << setw(20) << "Type"
            << setw(30) << "Region" << nl
            << setw(6) << "-----"
            << setw(20) << "----"
            << setw(30) << "------" << endl;

        const labelList& surfaceGeometry = surfaces.surfaces();
        const PtrList<dictionary>& surfacePatchInfo = surfaces.patchInfo();

        forAll(surfaceGeometry, surfI)
        {
            label geomI = surfaceGeometry[surfI];

            const wordList& regNames = allGeometry.regionNames()[geomI];

            Info<< surfaces.names()[surfI] << ':' << nl << nl;

            if (surfaces.faceZoneNames()[surfI].empty())
            {
                // 'Normal' surface
                forAll(regNames, i)
                {
                    label globalRegionI = surfaces.globalRegion(surfI, i);

                    label patchI;

                    if (surfacePatchInfo.set(globalRegionI))
                    {
                        patchI = meshRefiner.addMeshedPatch
                        (
                            regNames[i],
                            surfacePatchInfo[globalRegionI]
                        );
                    }
                    else
                    {
                        dictionary patchInfo;
                        patchInfo.set("type", wallPolyPatch::typeName);

                        patchI = meshRefiner.addMeshedPatch
                        (
                            regNames[i],
                            patchInfo
                        );
                    }

                    Info<< setf(ios_base::left)
                        << setw(6) << patchI
                        << setw(20) << mesh.boundaryMesh()[patchI].type()
                        << setw(30) << regNames[i] << nl;

                    globalToMasterPatch[globalRegionI] = patchI;
                    globalToSlavePatch[globalRegionI] = patchI;
                }
            }
            else
            {
                // Zoned surface
                forAll(regNames, i)
                {
                    label globalRegionI = surfaces.globalRegion(surfI, i);

                    // Add master side patch
                    {
                        label patchI;

                        if (surfacePatchInfo.set(globalRegionI))
                        {
                            patchI = meshRefiner.addMeshedPatch
                            (
                                regNames[i],
                                surfacePatchInfo[globalRegionI]
                            );
                        }
                        else
                        {
                            dictionary patchInfo;
                            patchInfo.set("type", wallPolyPatch::typeName);

                            patchI = meshRefiner.addMeshedPatch
                            (
                                regNames[i],
                                patchInfo
                            );
                        }

                        Info<< setf(ios_base::left)
                            << setw(6) << patchI
                            << setw(20) << mesh.boundaryMesh()[patchI].type()
                            << setw(30) << regNames[i] << nl;

                        globalToMasterPatch[globalRegionI] = patchI;
                    }
                    // Add slave side patch
                    {
                        const word slaveName = regNames[i] + "_slave";
                        label patchI;

                        if (surfacePatchInfo.set(globalRegionI))
                        {
                            patchI = meshRefiner.addMeshedPatch
                            (
                                slaveName,
                                surfacePatchInfo[globalRegionI]
                            );
                        }
                        else
                        {
                            dictionary patchInfo;
                            patchInfo.set("type", wallPolyPatch::typeName);

                            patchI = meshRefiner.addMeshedPatch
                            (
                                slaveName,
                                patchInfo
                            );
                        }

                        Info<< setf(ios_base::left)
                            << setw(6) << patchI
                            << setw(20) << mesh.boundaryMesh()[patchI].type()
                            << setw(30) << slaveName << nl;

                        globalToSlavePatch[globalRegionI] = patchI;
                    }
                }
            }

            Info<< nl;
        }
        Info<< "Added patches in = "
            << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
    }


    // Parallel
    // ~~~~~~~~

    // Decomposition
    autoPtr<decompositionMethod> decomposerPtr
    (
        decompositionMethod::New
        (
            decomposeDict
        )
    );
    decompositionMethod& decomposer = decomposerPtr();

    if (Pstream::parRun() && !decomposer.parallelAware())
    {
        FatalErrorIn(args.executable())
            << "You have selected decomposition method "
            << decomposer.typeName
            << " which is not parallel aware." << endl
            << "Please select one that is (hierarchical, ptscotch)"
            << exit(FatalError);
    }

    // Mesh distribution engine (uses tolerance to reconstruct meshes)
    fvMeshDistribute distributor(mesh, mergeDist);





    // Now do the real work -refinement -snapping -layers
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    const Switch wantRefine(meshDict.lookup("castellatedMesh"));
    const Switch wantSnap(meshDict.lookup("snap"));
    const Switch wantLayers(meshDict.lookup("addLayers"));

    // Refinement parameters
    const refinementParameters refineParams(refineDict);

    // Snap parameters
    const snapParameters snapParams(snapDict);


    if (wantRefine)
    {
        cpuTime timer;

        autoRefineDriver refineDriver
        (
            meshRefiner,
            decomposer,
            distributor,
            globalToMasterPatch,
            globalToSlavePatch
        );


        if (!overwrite && !debug)
        {
            const_cast<Time&>(mesh.time())++;
        }

        refineDriver.doRefine
        (
            refineDict,
            refineParams,
            snapParams,
            wantSnap,
            motionDict
        );

        writeMesh
        (
            "Refined mesh",
            meshRefiner,
            writeLevel,
            debug
        );

        Info<< "Mesh refined in = "
            << timer.cpuTimeIncrement() << " s." << endl;
    }

    if (wantSnap)
    {
        cpuTime timer;

        autoSnapDriver snapDriver
        (
            meshRefiner,
            globalToMasterPatch,
            globalToSlavePatch
        );

        if (!overwrite && !debug)
        {
            const_cast<Time&>(mesh.time())++;
        }

        // Use the resolveFeatureAngle from the refinement parameters
        scalar curvature = refineParams.curvature();

        snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);

        writeMesh
        (
            "Snapped mesh",
            meshRefiner,
            writeLevel,
            debug
        );

        Info<< "Mesh snapped in = "
            << timer.cpuTimeIncrement() << " s." << endl;
    }

    if (wantLayers)
    {
        cpuTime timer;

        autoLayerDriver layerDriver
        (
            meshRefiner,
            globalToMasterPatch,
            globalToSlavePatch
        );

        // Layer addition parameters
        layerParameters layerParams(layerDict, mesh.boundaryMesh());

        // Use the maxLocalCells from the refinement parameters
        bool preBalance = returnReduce
        (
            (mesh.nCells() >= refineParams.maxLocalCells()),
            orOp<bool>()
        );


        if (!overwrite &&  !debug)
        {
            const_cast<Time&>(mesh.time())++;
        }

        layerDriver.doLayers
        (
            layerDict,
            motionDict,
            layerParams,
            preBalance,
            decomposer,
            distributor
        );

        writeMesh
        (
            "Layer mesh",
            meshRefiner,
            writeLevel,
            debug
        );

        Info<< "Layers added in = "
            << timer.cpuTimeIncrement() << " s." << endl;
    }


    {
        // Check final mesh
        Info<< "Checking final mesh ..." << endl;
        faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
        const label nErrors = returnReduce
        (
            wrongFaces.size(),
            sumOp<label>()
        );

        if (nErrors > 0)
        {
            Info<< "Finished meshing with " << nErrors << " illegal faces"
                << " (concave, zero area or negative cell pyramid volume)"
                << endl;
            wrongFaces.write();
        }
        else
        {
            Info<< "Finished meshing without any errors" << endl;
        }
    }


    if (surfaceSimplify)
    {
        const polyBoundaryMesh& bMesh = mesh.boundaryMesh();

        labelHashSet includePatches(bMesh.size());

        if (args.optionFound("patches"))
        {
            includePatches = bMesh.patchSet
            (
                wordReList(args.optionLookup("patches")())
            );
        }
        else
        {
            forAll(bMesh, patchI)
            {
                const polyPatch& patch = bMesh[patchI];

                if (!isA<processorPolyPatch>(patch))
                {
                    includePatches.insert(patchI);
                }
            }
        }

        fileName outFileName
        (
            args.optionLookupOrDefault<fileName>
            (
                "outFile",
                "constant/triSurface/simplifiedSurface.stl"
            )
        );

        extractSurface
        (
            mesh,
            runTime,
            includePatches,
            outFileName
        );
    }


    Info<< "Finished meshing in = "
        << runTime.elapsedCpuTime() << " s." << endl;

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //