function [NumerPresEm,NumerVelCofRT0,NumerFlux,FluxDscp,LMCR] = ...
  Darcy_WG_RectQ0Q0RT0_PresVelFlux(RectMesh,PermK,sln)
%% Darcy: WG(Q0,Q0;RT[0]): Computing numer. pressure, velocity, flux, etc. 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;
% Assuming RectMesh.flag=3

%% "Computing" numerical pressure 
NumerPresEm = sln(1:NumEms);

%% Comput. emws.coeffs. in RT[0] nmlz.bas. for disc.wk.grad. of numer.sln. 
% NOTE: CofDWGsln has size NumEms*4 
CofDWGsln = WG_RectQ0Q0RT0_CofRT0NmlzBas_DiscWkGradSln(RectMesh,sln);

%% Computing proj.coeffs.(PermKt * RT[0]nmlz.bas.fxn.) in RT[0]nmlz.bas.fxn.
% NOTE: ProjCof has size NumEms*4*4 
ProjCof = Hdiv_RectRT0_NmlzBas_ProjMatK(RectMesh, -PermK);  % Note: Negative!

%% Computing numerical velocity coefficients in RT[0] nmlz.bas. 
NumerVelCofRT0 = zeros(NumEms,4);
for j=1:4  % The j-th RT[0] bas.fxn.
  NumerVelCofRT0(:,j) = sum(CofDWGsln.*squeeze(ProjCof(:,:,j)),2); 
end

%% Computing all elts. RT[0] 4 nmlz.bas.fxns. outer normal fluxes on 4 edges 
% NOTE: NBNF has size NumEms*4*4 
NBNF = Hdiv_RectRT0_NmlzBas_NmlFlux(RectMesh);

%% Computing: For all elements, outer normal fluxes on 4 edges 
NumerFlux = zeros(NumEms,4);
for j=1:4  % The j-th edge 
  NumerFlux(:,j) = sum(NumerVelCofRT0.*squeeze(NBNF(:,:,j)),2);
end

%% Computing (normal) flux discrepancies across all edges 
FluxDscp = zeros(NumEgs,1);
for j=1:4
  k = RectMesh.elem2edge(:,j);
  tmp = FluxDscp(k);
  FluxDscp(k) = tmp - NumerFlux(:,j);
end
% Correction for boundary edges 
k = find(RectMesh.BndryEdge>0);
FluxDscp(k) = 0;

%% "Computing" local-mass-conservation residuals for all elements 
% This quantity is 0 according to the theory of WG finiet element methods 
LMCR = zeros(NumEms,1);

return;