function [NumerPresEm,NumerVelEmCntr,FluxDscp,LMCR] ...
  = Darcy_CG_TriP1_PresVelFlux(EqnBC,TriMesh,PermK,sln,GAUSSQUAD)
%% Darcy: CG: TriP1: Computing numer.pres., vel.(em.cntr), flux dis., LMCR
% Computing: 
%   -- numer.pres.(emws.const.), 
%   -- vel.(em.cntr), 
%   -- flux discrepancy, 
%   -- LMCR (local-mass-conservation-residual) (For CG TriP1: LMCR = \int_E f)
% Assuming PermK is an elementwise constant 2x2 SPD matrix.  
% James Liu, ColoState; 2012/07--2018/12 

%% Retrieving mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;

%% Computing elementwise numerical pressure (const.) 
% by averaging the nodla values 
NumerPresEm = (1/3)*sum(sln(TriMesh.elem(:,1:3)),2);

%% Computing numerical velocity at element centers 
cof = sln(TriMesh.elem(:,1:3));
% NOTE: CP1BasGrad has size NumEms*3*2 
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
NumerPresGrad = zeros(NumEms,2);
NumerPresGrad(:,1) = sum(squeeze(CP1BasGrad(:,:,1)).*cof, 2);
NumerPresGrad(:,2) = sum(squeeze(CP1BasGrad(:,:,2)).*cof, 2);
NumerVelEmCntr = zeros(NumEms,2);
NumerVelEmCntr(:,1) = -PermK(:,1,1) .* NumerPresGrad(:,1) ... 
                      -PermK(:,1,2) .* NumerPresGrad(:,2) ;
NumerVelEmCntr(:,2) = -PermK(:,2,1) .* NumerPresGrad(:,1) ... 
                      -PermK(:,2,2) .* NumerPresGrad(:,2) ;

%% Computing flux discrepancy across interior edges 
FluxDscp = zeros(NumEgs,1);
vel1 = NumerVelEmCntr(TriMesh.edge2elem(:,1),:);
vel2 = zeros(NumEgs,2);
k = find(TriMesh.edge2elem(:,2)>0); 
vel2(k,:) = NumerVelEmCntr(TriMesh.edge2elem(k,2),:);
FluxDscp(k) = dot(vel1(k,:)-vel2(k,:), TriMesh.NmlEg(k,:), 2);

% for ig=1:NumEgs
%   if TriMesh.edge2elem(ig,2)>0 
%     vel2(ig,:) = NumerVelEmCntr(TriMesh.edge2elem(ig,2),:);
%   end
% end
% for ig=1:NumEgs
%   if TriMesh.edge2elem(ig,2)<=0 
%     FluxDscp(ig) = 0;
%   end
% end

%% Computing elementwise local-mass-conservation-residual (LMCR)
LMCR = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(:,3),:);
  fval = EqnBC.fxnf(qp);
  LMCR = LMCR + GAUSSQUAD.TRIG(k,4)*fval;
end
LMCR = LMCR .* TriMesh.area;

return;