//- (optionally destructively) construct from components
Foam::mapDistributePolyMesh::mapDistributePolyMesh
(
    const polyMesh& mesh,
    const label nOldPoints,
    const label nOldFaces,
    const label nOldCells,
    labelList& oldPatchStarts,
    labelList& oldPatchNMeshPoints,

    labelListList& subPointMap,
    labelListList& subFaceMap,
    labelListList& subCellMap,
    labelListList& subPatchMap,
    labelListList& constructPointMap,
    labelListList& constructFaceMap,
    labelListList& constructCellMap,
    labelListList& constructPatchMap,
    const bool reUse                // clone or reuse
)
:
    mesh_(mesh),
    nOldPoints_(nOldPoints),
    nOldFaces_(nOldFaces),
    nOldCells_(nOldCells),
    oldPatchSizes_(oldPatchStarts.size()),
    oldPatchStarts_(oldPatchStarts, reUse),
    oldPatchNMeshPoints_(oldPatchNMeshPoints, reUse),

    pointMap_(mesh.nPoints(), subPointMap, constructPointMap, reUse),
    faceMap_(mesh.nFaces(), subFaceMap, constructFaceMap, reUse),
    cellMap_(mesh.nCells(), subCellMap, constructCellMap, reUse),
    patchMap_(mesh.boundaryMesh().size(), subPatchMap, constructPatchMap, reUse)
{
    calcPatchSizes();
}

Foam::labelList Foam::structuredRenumber::renumber
(
    const polyMesh& mesh,
    const pointField& points
) const
{
    if (points.size() != mesh.nCells())
    {
        FatalErrorInFunction
            << "Number of points " << points.size()
            << " should equal the number of cells " << mesh.nCells()
            << exit(FatalError);
    }

    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
    const labelHashSet patchIDs(pbm.patchSet(patches_));

    label nFaces = 0;
    forAllConstIter(labelHashSet, patchIDs, iter)
    {
        nFaces += pbm[iter.key()].size();
    }


    // Extract a submesh.
    labelHashSet patchCells(2*nFaces);
    forAllConstIter(labelHashSet, patchIDs, iter)
    {
        const labelUList& fc = pbm[iter.key()].faceCells();
        forAll(fc, i)
        {
            patchCells.insert(fc[i]);
        }
    }

Foam::labelList Foam::metisDecomp::decompose
(
    const polyMesh& mesh,
    const pointField& points,
    const scalarField& pointWeights
)
{
    if (points.size() != mesh.nCells())
    {
        FatalErrorIn
        (
            "metisDecomp::decompose(const pointField&,const scalarField&)"
        )   << "Can use this decomposition method only for the whole mesh"
            << endl
            << "and supply one coordinate (cellCentre) for every cell." << endl
            << "The number of coordinates " << points.size() << endl
            << "The number of cells in the mesh " << mesh.nCells()
            << exit(FatalError);
    }

    CompactListList<label> cellCells;
    calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);

    // Decompose using default weights
    labelList decomp;
    decompose(cellCells.m(), cellCells.offsets(), pointWeights, decomp);

    return decomp;
}

Foam::label Foam::checkMeshQuality
(
    const polyMesh& mesh,
    const dictionary& dict
)
{
    label noFailedChecks = 0;

    {
        faceSet faces(mesh, "meshQualityFaces", mesh.nFaces()/100+1);
        motionSmoother::checkMesh(false, mesh, dict, faces);

        label nFaces = returnReduce(faces.size(), sumOp<label>());

        if (nFaces > 0)
        {
            noFailedChecks++;

            Info<< "  <<Writing " << nFaces
                << " faces in error to set " << faces.name() << endl;
            faces.instance() = mesh.pointsInstance();
            faces.write();
        }
    }

    return noFailedChecks;
}

// Write set to VTK readable files
void writeVTK
(
    const polyMesh& mesh,
    const topoSet& currentSet,
    const fileName& vtkName
)
{
    if (isA<faceSet>(currentSet))
    {
        // Faces of set with OpenFOAM faceID as value

        faceList setFaces(currentSet.size());
        labelList faceValues(currentSet.size());
        label setFaceI = 0;

        forAllConstIter(topoSet, currentSet, iter)
        {
            setFaces[setFaceI] = mesh.faces()[iter.key()];
            faceValues[setFaceI] = iter.key();
            setFaceI++;
        }

        primitiveFacePatch fp(setFaces, mesh.points());

        writePatch
        (
            true,
            currentSet.name(),
            fp,
            "faceID",
            faceValues,
            mesh.time().path()/vtkName
        );
    }

// 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();
    scalar writeTol = std::pow
    (
        scalar(10.0),
       -scalar(IOstream::defaultPrecision())
    );

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

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

    return mergeDist;
}

label addPointZone(const polyMesh& mesh, const word& name)
{
    label zoneID = mesh.pointZones().findZoneID(name);

    if (zoneID != -1)
    {
        Info<< "Reusing existing pointZone "
            << mesh.pointZones()[zoneID].name()
            << " at index " << zoneID << endl;
    }
    else
    {
        pointZoneMesh& pointZones = const_cast<polyMesh&>(mesh).pointZones();
        zoneID = pointZones.size();
        Info<< "Adding pointZone " << name << " at index " << zoneID << endl;

        pointZones.setSize(zoneID+1);
        pointZones.set
        (
            zoneID,
            new pointZone
            (
                name,
                labelList(0),
                zoneID,
                pointZones
            )
        );
    }
    return zoneID;
}

    bool MeshDistFromPatch::updateCell
    (
        const polyMesh& mesh,
        const label thisCellI,
        const label neighbourFaceI,
        const MeshDistFromPatch& neighbourInfo,
        const scalar tol
#ifdef FOAM_FACECELLWAVE_HAS_TRACKINGDATA
        ,TrackingData &td
#endif
    )
        {
            const scalar d=mag(
                mesh.cellCentres()[thisCellI]
                -
                mesh.faceCentres()[neighbourFaceI]
            );

            if(!valid(TRACKDATA)) {
                dist_=d+neighbourInfo.dist();
                return true;
            } else {
                const scalar nd=d+neighbourInfo.dist();
                if(nd<dist_) {
                    dist_=nd;
                    return true;
                } else {
                    return false;
                }
            }
        }

label addFaceZone(const polyMesh& mesh, const word& name)
{
    label zoneID = mesh.faceZones().findZoneID(name);

    if (zoneID != -1)
    {
        Info<< "Reusing existing faceZone " << mesh.faceZones()[zoneID].name()
            << " at index " << zoneID << endl;
    }
    else
    {
        faceZoneMesh& faceZones = const_cast<polyMesh&>(mesh).faceZones();
        zoneID = faceZones.size();
        Info<< "Adding faceZone " << name << " at index " << zoneID << endl;

        faceZones.setSize(zoneID+1);
        faceZones.set
        (
            zoneID,
            new faceZone
            (
                name,
                labelList(0),
                boolList(),
                zoneID,
                faceZones
            )
        );
    }
    return zoneID;
}

Foam::labelList Foam::CuthillMcKeeRenumber::renumber
(
    const polyMesh& mesh,
    const pointField& points
) const
{
    CompactListList<label> cellCells;
    decompositionMethod::calcCellCells
    (
        mesh,
        identity(mesh.nCells()),
        mesh.nCells(),
        false,                      // local only
        cellCells
    );

    labelList orderedToOld = bandCompression(cellCells());

    if (reverse_)
    {
        reverse(orderedToOld);
    }

    return orderedToOld;
}

Foam::sampledPatchInternalField::sampledPatchInternalField
(
    const word& name,
    const polyMesh& mesh,
    const dictionary& dict
)
:
    sampledPatch(name, mesh, dict),
    mappers_(patchIDs().size())
{
    const scalar distance = readScalar(dict.lookup("distance"));

    forAll(patchIDs(), i)
    {
        label patchI = patchIDs()[i];
        mappers_.set
        (
            i,
            new directMappedPatchBase
            (
                mesh.boundaryMesh()[patchI],
                mesh.name(),                        // sampleRegion
                directMappedPatchBase::NEARESTCELL, // sampleMode
                word::null,                         // samplePatch
                -distance                           // sample inside my domain
            )
        );
    }

Foam::points0MotionSolver::points0MotionSolver
(
    const polyMesh& mesh,
    const IOdictionary& dict,
    const word& type
)
:
    motionSolver(mesh, dict, type),
    points0_(pointIOField(points0IO(mesh)))
{
    if (points0_.size() != mesh.nPoints())
    {
        FatalErrorInFunction
            << "Number of points in mesh " << mesh.nPoints()
            << " differs from number of points " << points0_.size()
            << " read from file "
            <<  typeFilePath<pointIOField>
                (
                    IOobject
                    (
                        "points",
                        time().constant(),
                        polyMesh::meshSubDir,
                        mesh,
                        IOobject::MUST_READ,
                        IOobject::NO_WRITE,
                        false
                    )
                )
            << exit(FatalError);
    }
}

Foam::patchInjectionBase::patchInjectionBase
(
    const polyMesh& mesh,
    const word& patchName
)
:
    patchName_(patchName),
    patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
    patchArea_(0.0),
    patchNormal_(),
    cellOwners_(),
    triFace_(),
    triToFace_(),
    triCumulativeMagSf_(),
    sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
{
    if (patchId_ < 0)
    {
        FatalErrorInFunction
            << "Requested patch " << patchName_ << " not found" << nl
            << "Available patches are: " << mesh.boundaryMesh().names() << nl
            << exit(FatalError);
    }

    updateMesh(mesh);
}

// Construct from mesh.  No morphing data: the mesh has not changed
// HJ, 27/Nov/2009
Foam::mapPolyMesh::mapPolyMesh(const polyMesh& mesh)
:
    mesh_(mesh),
    morphing_(false),
    nOldPoints_(mesh.nPoints()),
    nOldFaces_(mesh.nFaces()),
    nOldCells_(mesh.nCells())
{}

Foam::sampledIsoSurface::sampledIsoSurface
(
    const word& name,
    const polyMesh& mesh,
    const dictionary& dict
)
:
    sampledSurface(name, mesh, dict),
    isoField_(dict.lookup("isoField")),
    isoVal_(readScalar(dict.lookup("isoValue"))),
    mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
    regularise_(dict.lookupOrDefault("regularise", true)),
    average_(dict.lookupOrDefault("average", false)),
    zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
    exposedPatchName_(word::null),
    surfPtr_(NULL),
    facesPtr_(NULL),
    prevTimeIndex_(-1),
    storedVolFieldPtr_(NULL),
    volFieldPtr_(NULL),
    storedPointFieldPtr_(NULL),
    pointFieldPtr_(NULL)
{
    if (!sampledSurface::interpolate())
    {
        FatalIOErrorIn
        (
            "sampledIsoSurface::sampledIsoSurface"
            "(const word&, const polyMesh&, const dictionary&)",
            dict
        )   << "Non-interpolated iso surface not supported since triangles"
            << " span across cells." << exit(FatalIOError);
    }

    if (zoneID_.index() != -1)
    {
        dict.lookup("exposedPatchName") >> exposedPatchName_;

        if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
        {
            FatalIOErrorIn
            (
                "sampledIsoSurface::sampledIsoSurface"
                "(const word&, const polyMesh&, const dictionary&)",
                dict
            )   << "Cannot find patch " << exposedPatchName_
                << " in which to put exposed faces." << endl
                << "Valid patches are " << mesh.boundaryMesh().names()
                << exit(FatalIOError);
        }

        if (debug && zoneID_.index() != -1)
        {
            Info<< "Restricting to cellZone " << zoneID_.name()
                << " with exposed internal faces into patch "
                << exposedPatchName_ << endl;
        }
    }

Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
const
{
    // Force calculation of tet-diag decomposition (for use in findCell)
    (void)mesh.tetBasePtIs();

    // Global calculation engine
    globalIndex globalCells(mesh.nCells());

    // Cell label per point
    labelList cellLabels(keepPoints_.size());

    forAll(keepPoints_, i)
    {
        const point& keepPoint = keepPoints_[i];

        label localCellI = mesh.findCell(keepPoint);

        label globalCellI = -1;

        if (localCellI != -1)
        {
            globalCellI = globalCells.toGlobal(localCellI);
        }

        reduce(globalCellI, maxOp<label>());

        if (globalCellI == -1)
        {
            FatalErrorInFunction
                << "Point " << keepPoint
                << " is not inside the mesh or on a face or edge." << nl
                << "Bounding box of the mesh:" << mesh.bounds()
                << exit(FatalError);
        }


        label procI = globalCells.whichProcID(globalCellI);
        label procCellI = globalCells.toLocal(procI, globalCellI);

        Info<< "Found point " << keepPoint << " in cell " << procCellI
            << " on processor " << procI << endl;


        if (globalCells.isLocal(globalCellI))
        {
            cellLabels[i] = localCellI;
        }
        else
        {
            cellLabels[i] = -1;
        }
    }
    return cellLabels;
}

bool Foam::meshStructure::isStructuredCell
(
    const polyMesh& mesh,
    const label layerI,
    const label cellI
) const
{
    const cell& cFaces = mesh.cells()[cellI];

    // Count number of side faces
    label nSide = 0;
    forAll(cFaces, i)
    {
        if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
        {
            nSide++;
        }
    }

    if (nSide != cFaces.size()-2)
    {
        return false;
    }

    // Check that side faces have correct point layers
    forAll(cFaces, i)
    {
        if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
        {
            const face& f = mesh.faces()[cFaces[i]];

            label nLayer = 0;
            label nLayerPlus1 = 0;
            forAll(f, fp)
            {
                label pointI = f[fp];
                if (pointLayer_[pointI] == layerI)
                {
                    nLayer++;
                }
                else if (pointLayer_[pointI] == layerI+1)
                {
                    nLayerPlus1++;
                }
            }

            if (f.size() != 4 || (nLayer+nLayerPlus1 != 4))
            {
                return false;
            }
        }
    }

Foam::labelList Foam::ptscotchDecomp::decompose
(
    const polyMesh& mesh,
    const labelList& agglom,
    const pointField& agglomPoints,
    const scalarField& pointWeights
)
{
    if (agglom.size() != mesh.nCells())
    {
        FatalErrorIn
        (
            "ptscotchDecomp::decompose(const labelList&, const pointField&)"
        )   << "Size of cell-to-coarse map " << agglom.size()
            << " differs from number of cells in mesh " << mesh.nCells()
            << exit(FatalError);
    }

//    // For running sequential ...
//    if (Pstream::nProcs() <= 1)
//    {
//        return scotchDecomp(decompositionDict_, mesh)
//            .decompose(agglom, agglomPoints, pointWeights);
//    }

    // Make Metis CSR (Compressed Storage Format) storage
    //   adjncy      : contains neighbours (= edges in graph)
    //   xadj(celli) : start of information in adjncy for celli
    CompactListList<label> cellCells;
    calcCellCells(mesh, agglom, agglomPoints.size(), cellCells);

    // Decompose using weights
    List<int> finalDecomp;
    decomposeZeroDomains
    (
        mesh.time().path()/mesh.name(),
        cellCells.m(),
        cellCells.offsets(),
        pointWeights,
        finalDecomp
    );

    // Rework back into decomposition for original mesh
    labelList fineDistribution(agglom.size());

    forAll(fineDistribution, i)
    {
        fineDistribution[i] = finalDecomp[agglom[i]];
    }

    return fineDistribution;
}

label findPatchID(const polyMesh& mesh, const word& name)
{
    label patchI = mesh.boundaryMesh().findPatchID(name);

    if (patchI == -1)
    {
        FatalErrorIn("findPatchID(const polyMesh&, const word&)")
            << "Cannot find patch " << name << endl
            << "Valid patches are " << mesh.boundaryMesh().names()
            << exit(FatalError);
    }
    return patchI;
}

void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
{
    // Create dummy system/fv*
    {
        IOobject io
        (
            "fvSchemes",
            mesh.time().system(),
            regionName,
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        );

        Info<< "Testing:" << io.objectPath() << endl;

        if (!io.headerOk())
        {
            Info<< "Writing dummy " << regionName/io.name() << endl;
            dictionary dummyDict;
            dictionary divDict;
            dummyDict.add("divSchemes", divDict);
            dictionary gradDict;
            dummyDict.add("gradSchemes", gradDict);
            dictionary laplDict;
            dummyDict.add("laplacianSchemes", laplDict);

            IOdictionary(io, dummyDict).regIOobject::write();
        }
    }
    {
        IOobject io
        (
            "fvSolution",
            mesh.time().system(),
            regionName,
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        );

        if (!io.headerOk())
        {
            Info<< "Writing dummy " << regionName/io.name() << endl;
            dictionary dummyDict;
            IOdictionary(io, dummyDict).regIOobject::write();
        }
    }
}