function [P,V]=Pres(Grid,S,Fluid,q)

% Compute K*lambda(S)
[Mw,Mo]=RelPerm(S,Fluid);
Mt=Mw+Mo;
KM = reshape([Mt,Mt,Mt]',3,Grid.Nx,Grid.Ny,Grid.Nz).*Grid.K;

% Compute pressure and extract fluxes
[P,V]=TPFA(Grid,KM,q);
