function [NumerPresEm,NumerVelEmCntr,NumerFlux,FluxDscp,LMCR] = ...
  Darcy_MFEM_RectRT0Q0_PresVelFlux(RectMesh,sln,BndryDescMat)
%% Darcy: MFEM(RT[0],Q0) on a rectangular mesh: 
% Computing numerical pressure, velocity, flux, etc. 
% JL20170211: TO BE FINISHED/REVISED 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;
% Assuming RectMesh.flag=3
LenEg = RectMesh.LenEg;
e2g = RectMesh.elem2edge;
sn = RectMesh.SignEmEg;

%% "Computing" numerical pressure: Elementwise constants (cell averages)
NumerPresEm = sln(NumEgs+(1:NumEms));

%% Computing numerical velocity at element centers 
cof = zeros(NumEms,4);
for j=1:4
  cof(:,j) = sln(RectMesh.elem2edge(:,j));
end
NumerVelEmCntr = zeros(NumEms,2);
NumerVelEmCntr(:,1) = -0.5*(cof(:,2)+cof(:,4));
NumerVelEmCntr(:,2) = -0.5*(cof(:,1)+cof(:,3));
% % JL20170211: TO BE REVISED FOR EFFICIENCY 
% for ie=1:NumEms
%   LblEg = RectMesh.elem2edge(ie,1:4);
%   coeff = sln(LblEg);
%   % JL20131026: Somewhere need reverse sign 
%   NumerVelEmCntr(ie,:) = -( coeff(1)*[0,0.5] + coeff(3)*[0,0.5] ...
%                           + coeff(2)*[0.5,0] + coeff(4)*[0.5,0] );
% end

%% Computing numerical flux: Elementwise 4 edges bulk normal fluxes (outward)
NumerFlux = zeros(NumEms,4);
coeff = sln(1:NumEgs);
for j=1:4
  NumerFlux(:,j) = coeff(e2g(:,j)).*sn(:,j).*LenEg(e2g(:,j));
end

% %% JL20140629: TO BE REVISED FOR EFFICIENCY 
% NumBndryPcs = size(BndryDescMat,1);
% NumerFlux.BndryFlux = zeros(NumBndryPcs,1);
% for ig=1:NumEgs
%   k = RectMesh.BndryEdge(ig);
%   if k>0 
%     ie1 = RectMesh.edge2elem(ig,1);
%     lbl = RectMesh.WhichEdge(ig,1);
%     NumerFlux.BndryFlux(k) = NumerFlux.BndryFlux(k)...
%       + NumerFlux.EmEgcNml(ie1,lbl)*LenEg(ig);
%   end
% end

%% Computing flux discrepancy aross edges: 0 theoretically  
FluxDscp = zeros(RectMesh.NumEgs,1);

%% Computing elementwise local mass-conservation residuals: 0 theoretically 
LMCR = zeros(RectMesh.NumEms,1);

return;