void Foam::fv::fixedTemperatureConstraint::constrain
(
    fvMatrix<scalar>& eqn,
    const label
)
{
    const basicThermo& thermo =
        mesh_.lookupObject<basicThermo>(basicThermo::dictName);

    switch (mode_)
    {
        case tmUniform:
        {
            const scalar t = mesh_.time().value();
            scalarField Tuni(cells_.size(), Tuniform_->value(t));
            eqn.setValues(cells_, thermo.he(thermo.p(), Tuni, cells_));

            break;
        }
        case tmLookup:
        {
            const volScalarField& T =
                mesh().lookupObject<volScalarField>(TName_);

            scalarField Tlkp(T, cells_);
            eqn.setValues(cells_, thermo.he(thermo.p(), Tlkp, cells_));

            break;
        }
        default:
        {
            // error handling done by NamedEnum
        }
    }
}

void Foam::fv::matrixChangeBeforeFvOption::setResidual(
    const fvMatrix<T>& matrix,
    autoPtr<GeometricField<T,fvPatchField,volMesh> >& val
) {
    Info << this->name() << " setting residual for "
        << matrix.psi().name() << " to " << this->fieldName() << endl;

    if(!val.valid()) {
        val.set(
            new GeometricField<T,fvPatchField,volMesh>(
                IOobject(
                    this->fieldName(),
                    matrix.psi().mesh().time().timeName(),
                    matrix.psi().mesh(),
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                matrix.psi().mesh(),
                dimensioned<T>(
                    "no",
                    matrix.dimensions()/dimVolume,
                    pTraits<T>::zero
                )
            )
        );
    }
    val()=this->calcResiduum(matrix);
}

void Foam::fv::explicitPorositySource::addSup
(
    fvMatrix<vector>& eqn,
    const label fieldi
)
{
    fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
    porosityPtr_->addResistance(porosityEqn);
    eqn -= porosityEqn;
}

void Foam::porousZone::modifyDdt(fvMatrix<Type>& m) const
{
    if (porosity_ < 1)
    {
        const labelList& cells = mesh_.cellZones()[cellZoneID_];

        forAll(cells, i)
        {
            m.diag()[cells[i]]   *= porosity_;
            m.source()[cells[i]] *= porosity_;
        }
    }

void Foam::fv::interRegionExplicitPorositySource::addSup
(
    fvMatrix<vector>& eqn,
    const label fieldI
)
{
    initialise();

    const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName_);

    const volVectorField& U = eqn.psi();

    volVectorField UNbr
    (
        IOobject
        (
            name_ + ":UNbr",
            nbrMesh.time().timeName(),
            nbrMesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        nbrMesh,
        dimensionedVector("zero", U.dimensions(), vector::zero)
    );

    // map local velocity onto neighbour region
    meshInterp().mapSrcToTgt
    (
        U.internalField(),
        plusEqOp<vector>(),
        UNbr.internalField()
    );

    fvMatrix<vector> nbrEqn(UNbr, eqn.dimensions());

    porosityPtr_->addResistance(nbrEqn);

    // convert source from neighbour to local region
    fvMatrix<vector> porosityEqn(U, eqn.dimensions());
    scalarField& Udiag = porosityEqn.diag();
    vectorField& Usource = porosityEqn.source();

    Udiag.setSize(eqn.diag().size(), 0.0);
    Usource.setSize(eqn.source().size(), vector::zero);

    meshInterp().mapTgtToSrc(nbrEqn.diag(), plusEqOp<scalar>(), Udiag);
    meshInterp().mapTgtToSrc(nbrEqn.source(), plusEqOp<vector>(), Usource);

    eqn -= porosityEqn;
}

void Foam::jjc2014Zone::modifyDdt(fvMatrix<Type>& m) const
{
    if (porosity_ < 1)
    {
        const labelList& cells = mesh_.cellZones()[cellZoneID_];

        scalar coeff( (1.0 + addedMassCoeff_ * ( 1 - porosity_ ) / porosity_) / porosity_ );

        forAll(cells, i)
        {
            m.diag()[cells[i]]   *= coeff;
            m.source()[cells[i]] *= coeff;
        }
    }

void Foam::fv::tabulatedAccelerationSource::addSup
(
    const RhoFieldType& rho,
    fvMatrix<vector>& eqn,
    const label fieldi
)
{
    Vector<vector> acceleration(motion_.acceleration());

    // If gravitational force is present combine with the linear acceleration
    if (mesh_.foundObject<uniformDimensionedVectorField>("g"))
    {
        uniformDimensionedVectorField& g =
            mesh_.lookupObjectRef<uniformDimensionedVectorField>("g");

        const uniformDimensionedScalarField& hRef =
            mesh_.lookupObject<uniformDimensionedScalarField>("hRef");

        g = g0_ - dimensionedVector("a", dimAcceleration, acceleration.x());

        dimensionedScalar ghRef
        (
            mag(g.value()) > SMALL
          ? g & (cmptMag(g.value())/mag(g.value()))*hRef
          : dimensionedScalar("ghRef", g.dimensions()*dimLength, 0)
        );

        mesh_.lookupObjectRef<volScalarField>("gh") = (g & mesh_.C()) - ghRef;

        mesh_.lookupObjectRef<surfaceScalarField>("ghf") =
            (g & mesh_.Cf()) - ghRef;
    }
    // ... otherwise include explicitly in the momentum equation
    else
    {
        eqn -= rho*dimensionedVector("a", dimAcceleration, acceleration.x());
    }

    dimensionedVector Omega
    (
        "Omega",
        dimensionSet(0, 0, -1, 0, 0),
        acceleration.y()
    );

    dimensionedVector dOmegaDT
    (
        "dOmegaDT",
        dimensionSet(0, 0, -2, 0, 0),
        acceleration.z()
    );

    eqn -=
    (
        rho*(2*Omega ^ eqn.psi())         // Coriolis force
      + rho*(Omega ^ (Omega ^ mesh_.C())) // Centrifugal force
      + rho*(dOmegaDT ^ mesh_.C())        // Angular tabulatedAcceleration force
    );
}

Foam::SolverPerformance<Type> Foam::solve
(
    fvMatrix<Type>& fvm,
    const word& name
)
{
    return fvm.solve(name);
}

void Foam::fixedInternalValueFvPatchField<Type>::manipulateMatrix
(
    fvMatrix<Type>& matrix
)
{
    // Apply the refValue as a constraint in the matrix
    matrix.setValues(this->patch().faceCells(), refValue_);
}

void Foam::actuationDiskSource::addSu(fvMatrix<vector>& UEqn)
{
    bool compressible = false;
    if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
    {
        compressible = true;
    }

    const scalarField& cellsV = this->mesh().V();
    vectorField& Usource = UEqn.source();
    const vectorField& U = UEqn.psi();

    if (V() > VSMALL)
    {
        if (compressible)
        {
            addActuationDiskAxialInertialResistance
            (
                Usource,
                cells_,
                cellsV,
                this->mesh().lookupObject<volScalarField>("rho"),
                U
            );
        }
        else
        {
            addActuationDiskAxialInertialResistance
            (
                Usource,
                cells_,
                cellsV,
                geometricOneField(),
                U
            );
        }
    }
}

void forceEquation<T>::operator()(fvMatrix<T> &eq)
{
    typedef GeometricField<T,fvPatchField,volMesh> resultField;

    clearVariables();

    DynamicList<label> cellIDs;

    if(!getMask(cellIDs,eq.psi().name())) {
        return;
    }

    Field<T> values(cellIDs.size());

    parse(valueExpression_);
    const resultField &calculated=getResult<resultField>();

    forAll(cellIDs,i) {
        values[i]=calculated[cellIDs[i]];
    }

    eq.setValues(cellIDs,values);
}

void Foam::fv::radiation::addSup
(
    const volScalarField& rho,
    fvMatrix<scalar>& eqn,
    const label fieldi
)
{
    const basicThermo& thermo =
        mesh_.lookupObject<basicThermo>(basicThermo::dictName);

    radiation_->correct();

    eqn += radiation_->Sh(thermo, eqn.psi());
}

void Foam::fv::explicitPorositySource::addSup
(
    fvMatrix<vector>& eqn,
    const label fieldI
)
{
    fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());

    if (eqn.dimensions() == dimForce)
    {
        const volScalarField& rho =
            mesh_.lookupObject<volScalarField>(rhoName_);
        const volScalarField& mu = mesh_.lookupObject<volScalarField>(muName_);

        porosityPtr_->addResistance(porosityEqn, rho, mu);
    }
    else
    {
        porosityPtr_->addResistance(porosityEqn);
    }

    eqn -= porosityEqn;
}

void SwakSetValue<T>::setValue
(
    fvMatrix<T>& eqn,
    const label fieldI
)
{
    this->driver().clearVariables();
    this->driver().parse(this->expressions_[fieldI]);
    if(
        !this->driver().
        FieldValueExpressionDriver::resultIsTyp<typename SwakSetValue<T>::resultField>()
    ) {
        FatalErrorIn("SwakSetValue<"+word(pTraits<T>::typeName)+">::setValue()")
            << "Result of " << this->expressions_[fieldI] << " is not a "
                << pTraits<T>::typeName
                << endl
                << exit(FatalError);
    }

    typename SwakSetValue<T>::resultField result(
        this->driver().
        FieldValueExpressionDriver::getResult<typename SwakSetValue<T>::resultField>()
    );

    DynamicList<label> cellIDs;

    if(useMaskExpression_) {
        if(!getMask(cellIDs,eqn.psi().name())) {
            return;
        }
    } else {
        cellIDs=this->cells_;
    }

    List<T> values(cellIDs.size());

    //    UIndirectList<Type>(values, cells_) = injectionRate_[fieldI];
    forAll(cellIDs,i)
    {
	label cellI=cellIDs[i];

        values[i]=result[cellI];
    }

void fvBlockMatrix<Type>::insertDiagSource
(
    const direction dir,
    fvMatrix<matrixType>& matrix
)
{
    matrix.completeAssembly();

    // Save a copy for different components
    scalarField& diag = matrix.diag();
    scalarField saveDiag(diag);

    // Add source boundary contribution
    Field<matrixType>& source = matrix.source();
    matrix.addBoundarySource(source, false);

    const direction nCmpts = pTraits<matrixType>::nComponents;
    direction localDir = dir;

    // Get reference to this source field of block system
    Field<Type>& b = this->source();

    if
    (
        // This is needed if the matrixType is <vector>, then you need to grab
        // coeffs as linear. Consider doing a matrixType check also.
        // VV, 17/March/2014
        this->diag().activeType() != blockCoeffBase::SQUARE
    )
    {
        typename CoeffField<Type>::linearTypeField& blockDiag =
            this->diag().asLinear();

        for (direction cmptI = 0; cmptI < nCmpts; cmptI++)
        {
            matrix.addBoundaryDiag(diag, cmptI);
            scalarField sourceCmpt(source.component(cmptI));

//            FieldField<Field, scalar> bouCoeffsCmpt
//            (
//                matrix.boundaryCoeffs().component(cmptI)
//            );

            // Possible problem for coupled non-aligned boundaries.
            // VV, 14/May/2014.
//            matrix.correctImplicitBoundarySource
//            (
//                bouCoeffsCmpt,
//                sourceCmpt,
//                cmptI
//            );

            forAll (diag, cellI)
            {
                blockDiag[cellI](localDir) = diag[cellI];
                b[cellI](localDir) += sourceCmpt[cellI];
            }

            localDir++;

            // Reset diagonal
            diag = saveDiag;
        }