// Solve the Momentum equation
const unallocLabelList& owner = mesh.owner();
const unallocLabelList& neighbour = mesh.neighbour();

const scalarField& volume = mesh.V();

const surfaceVectorField& Sf = mesh.Sf();

tmp<surfaceScalarField> tweights = mesh.weights();
const surfaceScalarField& weights = tweights();

tmp<surfaceVectorField> tdelta = mesh.delta();
const surfaceVectorField& delta = tdelta();



volVectorField B = - beta*(T-To)*g;
surfaceVectorField B_f = linearInterpolate(B);
volVectorField B_reconstructed (
                                IOobject
                                (
                                 "Breconstructed",
                                 runTime.timeName(),
                                 mesh,
                                 IOobject::NO_READ,
                                 IOobject::AUTO_WRITE
                                 ),
                                mesh,
                                dimensionedVector("zero", B.dimensions(), vector::zero),
                                "zeroGradient"
                                
                                );
// assemble B_reconstruted
forAll(owner, facei)
{
    vector SfBfdf = (B_f[facei]&delta[facei])*Sf[facei];
    
    B_reconstructed[owner[facei]] +=  weights[facei]*SfBfdf;
    B_reconstructed[neighbour[facei]] += (1.0-weights[facei])*SfBfdf;
}


forAll(mesh.boundary(), patchi)
{
    const labelUList& pFaceCells = mesh.boundary()[patchi].faceCells();
    
    const vectorField& pSf = Sf.boundaryField()[patchi];
    const vectorField& pdeltaf = delta.boundaryField()[patchi];
    const vectorField& pB_f = B_f.boundaryField()[patchi];
    
    forAll(mesh.boundary()[patchi], facei)
    {
        vector SfBfdf =pSf[facei]*(pB_f[facei]&pdeltaf[facei]);
        B_reconstructed[pFaceCells[facei]] += SfBfdf;
    }
}


B_reconstructed.internalField() /= volume;
B_reconstructed.correctBoundaryConditions();


fvVectorMatrix UEqn
(
        fvm::ddt(U)
      + fvm::div(mdotf, U)
      - fvm::laplacian(nu, U)
      + fvm::SuSp(-fvc::div(mdotf),U)
 == - fvc::grad(p) + B_reconstructed

);

//- Implicit relaxation
UEqn.relax();

//- Solution of the linear solver adding the pressure gradient as source term
UEqn.solve();
