function [NumerPresEm,NumerVelCofRT0,NumerFlux,FluxDscp,LMCR] = ...
  Darcy_WG_TriP0P0RT0_PresVelFlux(TriMesh,PermK,sln)
%% Darcy: WG(P0,P0;RT0) tri.mesh: Computing numerical pressure, velocity, flux 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;

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

%% Preparation 
% Computing coeffs. in RT0 nmlz.bas. for disc.wk.grad. of numer.sln. 
% NOTE: CDWGS has size NumEms*3 
CDWGS = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradSln(TriMesh,sln);
% Computing projection coeffs. of (RT0 bas.fxns. multiplied by -PermKt) 
% NOTE: ProjCof has size NumEms*3*3 
ProjCof = Hdiv_TriRT0_NmlzBas_ProjMatK(TriMesh,-PermK);  % NOTE: Negative!

%% Computing numerical velocity: A structure  NumerVel .CofRT0 .EmCntr 
NumerVelCofRT0 = zeros(NumEms,3);
for j=1:3
  NumerVelCofRT0(:,j) = sum(CDWGS(:,1:3).*squeeze(ProjCof(:,1:3,j)),2);
end
NumerVelEmCntr = NumerVelCofRT0(:,1:2);

%% Computing all elts. RT0 3 nmlz.bas.fxns. outer normal fluxes on 3 edges 
% NOTE: NBNF has size NumEms*3*3 
NBNF = Hdiv_TriRT0_NmlzBas_NmlFlux(TriMesh);

%% Computing elementwise 3 edge outer-normal-fluxes (length incorporated) 
NumerFlux = zeros(NumEms,3);
for j=1:3
  NumerFlux(:,j) = sum(NumerVelCofRT0.*squeeze(NBNF(:,:,j)),2);
end

%% Computing (normal) flux discrepancies across all edges 
% This quantity is 0 by the theory of WG finite element methods 
FluxDscp = zeros(NumEgs,1);
% for j=1:3
%   k = TriMesh.elem2edge(:,j);
%   tmp = FluxDscp(k);
%   FluxDscp(k) = tmp + NumerFlux(:,j);
% end
% % Correction for boundary edges 
% k = find(TriMesh.BndryEdge>0);
% FluxDscp(k) = 0;

%% Computing local-mass-conservation residual 
% This quantity is 0 by the theory of WG finite element methods 
LMCR = zeros(NumEms,1);

return;