Foam::Kmesh::Kmesh(const fvMesh& mesh)
:
    vectorField(mesh.V().size()),
    nn_(vector::dim)
{
    boundBox box = mesh.bounds();
    l_ = box.span();

    vector cornerCellCentre = ::Foam::max(mesh.C().internalField());
    vector cellL = 2*(box.max() - cornerCellCentre);

    vector rdeltaByL;
    label nTot = 1;

    forAll(nn_, i)
    {
        nn_[i] = label(l_[i]/cellL[i] + 0.5);
        nTot *= nn_[i];

        if (nn_[i] > 1)
        {
            l_[i] -= cellL[i];
        }

        rdeltaByL[i] = nn_[i]/(l_[i]*l_[i]);
    }

void replaceBoundaryType
(
    const fvMesh& mesh,
    const word& fieldName,
    const word& boundaryType,
    const string& boundaryValue
)
{
    IOobject header
    (
        fieldName,
        mesh.time().timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    );

    if (!header.headerOk())
    {
        return;
    }

    Info<< "Updating boundary types for field " << header.name() << endl;

    const word oldTypeName = IOdictionary::typeName;
    const_cast<word&>(IOdictionary::typeName) = word::null;

    IOdictionary dict(header);

    const_cast<word&>(IOdictionary::typeName) = oldTypeName;
    const_cast<word&>(dict.type()) = dict.headerClassName();

    // Make a backup of the old file
    if (mvBak(dict.objectPath(), "old"))
    {
        Info<< "    Backup original file to "
            << (dict.objectPath() + ".old") << endl;
    }

    // Loop through boundary patches and update
    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
    dictionary& boundaryDict = dict.subDict("boundaryField");
    forAll(bMesh, patchI)
    {
        if (isA<wallPolyPatch>(bMesh[patchI]))
        {
            word patchName = bMesh[patchI].name();
            dictionary& oldPatch = boundaryDict.subDict(patchName);

            dictionary newPatch(dictionary::null);
            newPatch.add("type", boundaryType);
            newPatch.add("value", ("uniform " + boundaryValue).c_str());

            oldPatch = newPatch;
        }
    }

    Info<< "    writing updated " << dict.name() << nl << endl;
    dict.regIOobject::write();
}

Foam::anisotropicFilter::anisotropicFilter
(
    const fvMesh& mesh,
    const dictionary& bd
)
:
    LESfilter(mesh),
    widthCoeff_(readScalar(bd.subDict(type() + "Coeffs").lookup("widthCoeff"))),
    coeff_
    (
        IOobject
        (
            "anisotropicFilterCoeff",
            mesh.time().timeName(),
            mesh
        ),
        mesh,
        dimensionedVector("zero", dimLength*dimLength, Zero),
        calculatedFvPatchScalarField::typeName
    )
{
    for (direction d=0; d<vector::nComponents; d++)
    {
        coeff_.internalField().replace
        (
            d,
            (1/widthCoeff_)*
            sqr
            (
                2.0*mesh.V()
               /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
            )
        );
    }
}

Foam::psiThermo::psiThermo(const fvMesh& mesh, const word& phaseName)
:
    fluidThermo(mesh, phaseName),

    psi_
    (
        IOobject
        (
            phasePropertyName("thermo:psi"),
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionSet(0, -2, 2, 0, 0)
    ),

    mu_
    (
        IOobject
        (
            phasePropertyName("thermo:mu"),
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionSet(1, -1, -1, 0, 0)
    )
{}

// Construct from components
    Foam::contactPatchPair::contactPatchPair
(
    const fvMesh& m,
    const label master,
    const label slave,
    const scalar tolerance
)
:
    mesh_(m),
    masterPatchIndex_(master),
    slavePatchIndex_(slave),
    masterInterpolate_(m.boundaryMesh()[masterPatchIndex_]),
    slaveInterpolate_(m.boundaryMesh()[slavePatchIndex_]),
    patchToPatchInterpolate_
    (
        m.boundaryMesh()[masterPatchIndex_],
        m.boundaryMesh()[slavePatchIndex_],
        intersection::FULL_RAY,
        intersection::CONTACT_SPHERE
    ),
    tol_(tolerance),
    touchFraction_
    (
        mesh_.boundaryMesh()[slavePatchIndex_].size(),
        pTraits<scalar>::zero
    ),
    slaveDisplacement_
    (
        mesh_.boundaryMesh()[slavePatchIndex_].size(),
        pTraits<vector>::zero
    )
{}

Foam::baseAtmosphere::baseAtmosphere
(
    const wordList& partNames,
    const fvMesh& mesh,
    const dictionary dict
)
:
    PtrList<fluidSpecie>(partNames.size())
{
    for(label ip = 0; ip < size(); ip++)
    {
        set
        (
            ip,
            new fluidSpecie
            (
                partNames[ip],
                IOobject(partNames[ip]+"VapourRho", mesh.time().timeName(), mesh,
                         IOobject::MUST_READ, IOobject::AUTO_WRITE),
                IOobject(partNames[ip]+"LiquidFrac", mesh.time().timeName(), mesh,
                         IOobject::MUST_READ, IOobject::AUTO_WRITE),
                mesh,
                dict.subDict(partNames[ip])
            )
        );
    }
}

Foam::basicChemistryModel::basicChemistryModel
(
    const fvMesh& mesh,
    const word& phaseName
)
    :
    IOdictionary
    (
       IOobject
       (
           IOobject::groupName("chemistryProperties", phaseName),
           mesh.time().constant(),
           mesh,
           IOobject::MUST_READ_IF_MODIFIED,
           IOobject::NO_WRITE
       )
    ),
    mesh_(mesh),
    chemistry_(lookup("chemistry")),
    deltaTChemIni_(readScalar(lookup("initialChemicalTimeStep"))),
    deltaTChem_
    (
       IOobject
       (
           IOobject::groupName("deltaTChem", phaseName),
           mesh.time().constant(),
           mesh,
           IOobject::NO_READ,
           IOobject::NO_WRITE
       ),
       mesh,
       dimensionedScalar("deltaTChem0", dimTime, deltaTChemIni_)
    )
{}

Foam::anisotropicFilter::anisotropicFilter
(
    const fvMesh& mesh,
    scalar widthCoeff
)
:
    LESfilter(mesh),
    widthCoeff_(widthCoeff),
    coeff_
    (
        IOobject
        (
            "anisotropicFilterCoeff",
            mesh.time().timeName(),
            mesh
        ),
        mesh,
        dimensionedVector("zero", dimLength*dimLength, vector::zero),
        calculatedFvPatchVectorField::typeName
    )
{
    for (direction d=0; d<vector::nComponents; d++)
    {
        coeff_.internalField().replace
        (
            d,
            (2.0/widthCoeff_)*mesh.V()
           /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
        );
    }
}

hhuCombustionThermo::hhuCombustionThermo(const fvMesh& mesh)
:
    hCombustionThermo(mesh),

    Tu_
    (
        IOobject
        (
            "Tu",
            mesh.time().timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    ),

    hu_
    (
        IOobject
        (
            "hu",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionSet(0, 2, -2, 0, 0),
        huBoundaryTypes()
    )
{}

Foam::scalar Foam::compressibleCourantNo(
    const fvMesh& mesh,
    const Time& runTime,
    const volScalarField& rho,
    const surfaceScalarField& phi
)
{
    scalar CoNum = 0.0;
    scalar meanCoNum = 0.0;

    //- Can have fluid domains with 0 cells so do not test.
    //if (mesh.nInternalFaces())
    {
        surfaceScalarField SfUfbyDelta =
            mesh.surfaceInterpolation::deltaCoeffs()
          * mag(phi)
          / fvc::interpolate(rho);

        CoNum = max(SfUfbyDelta/mesh.magSf())
            .value()*runTime.deltaT().value();

        meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
            .value()*runTime.deltaT().value();
    }

    Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
        << " max: " << CoNum << endl;

    return CoNum;
}

vectorField fieldOperations::
getWallPointMotion
(
  const fvMesh& mesh,
  const volScalarField& C, 
  const label movingPatchID
)
{
  // interpolate the concentration from cells to wall faces
  coupledPatchInterpolation patchInterpolator
  (
    mesh.boundaryMesh()[movingPatchID], mesh 
  );

  // concentration and normals on the faces
  scalarField pointCface = -C.boundaryField()[movingPatchID].snGrad();
  vectorField pointNface = mesh.boundaryMesh()[movingPatchID].faceNormals();
  
  scalarField motionC = patchInterpolator.faceToPointInterpolate(pointCface);
  vectorField motionN = patchInterpolator.faceToPointInterpolate(pointNface);
  
  // normalize point normals to 1
  forAll(motionN, ii) motionN[ii]/=mag(motionN[ii]);
  
  return motionC*motionN;
}

// Inplace add mesh1 to mesh0
Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::fvMeshAdder::add
(
    fvMesh& mesh0,
    const fvMesh& mesh1,
    const faceCoupleInfo& coupleInfo,
    const bool validBoundary
)
{
    mesh0.clearOut();

    // Resulting merged mesh (polyMesh only!)
    autoPtr<mapAddedPolyMesh> mapPtr
    (
        polyMeshAdder::add
        (
            mesh0,
            mesh1,
            coupleInfo,
            validBoundary
        )
    );

    // Adjust the fvMesh part.
    const polyBoundaryMesh& patches = mesh0.boundaryMesh();

    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh0.boundary());
    fvPatches.setSize(patches.size());
    forAll(patches, patchI)
    {
        fvPatches.set(patchI, fvPatch::New(patches[patchI], fvPatches));
    }

void Foam::patchProbes::findElements(const fvMesh& mesh)
{
    (void)mesh.tetBasePtIs();

    const polyBoundaryMesh& bm = mesh.boundaryMesh();

    label patchi = bm.findPatchID(patchName_);

    if (patchi == -1)
    {
        FatalErrorInFunction
            << " Unknown patch name "
            << patchName_ << endl
            << exit(FatalError);
    }

     // All the info for nearest. Construct to miss
    List<mappedPatchBase::nearInfo> nearest(this->size());

    const polyPatch& pp = bm[patchi];

    if (pp.size() > 0)
    {
        labelList bndFaces(pp.size());
        forAll(bndFaces, i)
        {
            bndFaces[i] =  pp.start() + i;
        }

Foam::twoPhaseMixture::twoPhaseMixture
(
    const fvMesh& mesh,
    const dictionary& dict
)
:
    phase1Name_(wordList(dict.lookup("phases"))[0]),
    phase2Name_(wordList(dict.lookup("phases"))[1]),

    alpha1_
    (
        IOobject
        (
            IOobject::groupName("alpha", phase1Name_),
            mesh.time().timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    ),

    alpha2_
    (
        IOobject
        (
            IOobject::groupName("alpha", phase2Name_),
            mesh.time().timeName(),
            mesh
        ),
        1.0 - alpha1_
    )
{}

void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{

    Info<< "\tdomain volume = " << gSum(mesh.V()) << endl;

    const cellZoneMesh& cellZones = mesh.cellZones();
    scalar cellZoneVolTotal(0);
    if (cellZones.size())
    {
        Info  << "cellZones:" << endl;
        forAll(cellZones, i)
        {
            const cellZone& zone = mesh.cellZones()[i];
            const cellZoneMesh& zoneMesh = zone.zoneMesh();
            const labelList& cellsZone = zoneMesh[i];
            scalar cellZoneVol(0);
            //float cellZoneVol=0.0;
            forAll(cellsZone, cI)
            { 
                cellZoneVol += mesh.V()[cellsZone[cI]];
            }
            Info  << '\t' << zone.name() << " volume = " << cellZoneVol << endl;
            cellZoneVolTotal += cellZoneVol;
        }
    }

// Construct from components
cfdemCloudIB::cfdemCloudIB
(
    const fvMesh& mesh
)
:
    cfdemCloud(mesh),
    angularVelocities_(NULL),
    pRefCell_(readLabel(mesh.solutionDict().subDict("PISO").lookup("pRefCell"))),
    pRefValue_(readScalar(mesh.solutionDict().subDict("PISO").lookup("pRefValue")))
{}

Foam::AveragingMethods::Basic<Type>::Basic
(
    const IOobject& io,
    const dictionary& dict,
    const fvMesh& mesh
)
:
    AveragingMethod<Type>(io, dict, mesh, labelList(1, mesh.nCells())),
    data_(FieldField<Field, Type>::operator[](0)),
    dataGrad_(mesh.nCells())
{}

void Foam::quadraticFitSnGradData::findFaceDirs
(
    vector& idir,        // value changed in return
    vector& jdir,        // value changed in return
    vector& kdir,        // value changed in return
    const fvMesh& mesh,
    const label faci
)
{
    idir = mesh.Sf()[faci];
    idir /= mag(idir);

    #ifndef SPHERICAL_GEOMETRY
        if (mesh.nGeometricD() <= 2) // find the normal direcion
        {
            if (mesh.geometricD()[0] == -1)
            {
                kdir = vector(1, 0, 0);
            }
            else if (mesh.geometricD()[1] == -1)
            {
                kdir = vector(0, 1, 0);
            }
            else
            {
                kdir = vector(0, 0, 1);
            }
        }
        else // 3D so find a direction in the plane of the face
        {
            const face& f = mesh.faces()[faci];
            kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
        }
    #else
        // Spherical geometry so kdir is the radial direction
        kdir = mesh.Cf()[faci];
    #endif

    if (mesh.nGeometricD() == 3)
    {
        // Remove the idir component from kdir and normalise
        kdir -= (idir & kdir)*idir;

        scalar magk = mag(kdir);

        if (magk < SMALL)
        {
            FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
                << exit(FatalError);
        }
        else
        {
            kdir /= magk;
        }
    }

    jdir = kdir ^ idir;
}

void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{
    IOobject phiHeader
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ
    );

    IOobject faceIbMaskHeader
    (
        "faceIbMask",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ
    );

    if (phiHeader.headerOk() && faceIbMaskHeader.headerOk())
    {
        Info<< "    Reading phi" << endl;
        surfaceScalarField phi(phiHeader, mesh);

        Info<< "    Reading faceIbMask" << endl;
        surfaceScalarField faceIbMask(faceIbMaskHeader, mesh);

        volScalarField contErr = fvc::div(faceIbMask*phi);

        scalar sumLocalContErr = runTime.deltaT().value()*
            mag(contErr)().weightedAverage(mesh.V()).value();

        scalar globalContErr = runTime.deltaT().value()*
            contErr.weightedAverage(mesh.V()).value();

        Info<< "IB time step continuity errors : sum local = "
            << sumLocalContErr
            << ", global = " << globalContErr
            << endl;

        volScalarField magContErr
        (
            "magContErr",
            mag(contErr)
        );
        magContErr.write();

    }
    else
    {
        Info<< "    No phi or faceIbMask" << endl;
    }

}

void Foam::fvMeshTools::addPatchFields
(
    fvMesh& mesh,
    const dictionary& patchFieldDict,
    const word& defaultPatchFieldType,
    const typename GeoField::value_type& defaultPatchValue
)
{
    HashTable<GeoField*> flds
    (
        mesh.objectRegistry::lookupClass<GeoField>()
    );

    forAllIter(typename HashTable<GeoField*>, flds, iter)
    {
        GeoField& fld = *iter();

        typename GeoField::GeometricBoundaryField& bfld =
            fld.boundaryField();

        label sz = bfld.size();
        bfld.setSize(sz+1);

        if (patchFieldDict.found(fld.name()))
        {
            bfld.set
            (
                sz,
                GeoField::PatchFieldType::New
                (
                    mesh.boundary()[sz],
                    fld.dimensionedInternalField(),
                    patchFieldDict.subDict(fld.name())
                )
            );
        }
        else
        {
            bfld.set
            (
                sz,
                GeoField::PatchFieldType::New
                (
                    defaultPatchFieldType,
                    mesh.boundary()[sz],
                    fld.dimensionedInternalField()
                )
            );
            bfld[sz] == defaultPatchValue;
        }
    }
}