function [NumerPresEm,NumerVelEmCntr,NumerFlux,FluxDscp,LMCR] = ...
  Darcy_MFEM_TriRT0P0_PresVelFlux(BndryDescMat,TriMesh,sln)
% Darcy: MFEM(RT0,P0): Computing (numerical) pressure, velocity, flux, etc. 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
LenEg = TriMesh.LenEg;
e2g = TriMesh.elem2edge;
sn = TriMesh.SignEmEg;
area = TriMesh.area;

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

%% Computing numerical velocity at element centers 
CofRT0 = zeros(NumEms,3);
for j=1:3
  CofRT0(:,j) = sln(TriMesh.elem2edge(:,j)).*TriMesh.SignEmEg(:,j);
end
BasFxn = zeros(NumEms,3,2);  % All elts., 3 edges, 2-vec. 
cntr = (1/3)*(TriMesh.node(TriMesh.elem(:,1),:)...
            + TriMesh.node(TriMesh.elem(:,2),:)...
            + TriMesh.node(TriMesh.elem(:,3),:));
for j=1:3
  %% JL20160614: THERE IS AN ISSUE ABOUT NEGATIVE SIGN !!!
  coeff = -LenEg(TriMesh.elem2edge(:,j))./(2*area);
  BasFxn(:,j,:) = [coeff,coeff].*(cntr-TriMesh.node(TriMesh.elem(:,j),:));
end
% NumerVelEmCntr = zeros(NumEms,2);
NumerVelEmCntr = [CofRT0(:,1),CofRT0(:,1)] .* squeeze(BasFxn(:,1,:))...
               + [CofRT0(:,2),CofRT0(:,2)] .* squeeze(BasFxn(:,2,:))...
               + [CofRT0(:,3),CofRT0(:,3)] .* squeeze(BasFxn(:,3,:));
NumerVelEmCntr = -NumerVelEmCntr;
% % JL20131003: OLD-STYLE 
% for ie=1:NumEms
%   LblVrtx = TriMesh.elem(ie,1:3);
%   vrtx = TriMesh.node(LblVrtx,1:2);
%   cntr = (1.0/3)*sum(vrtx,1);
%   area = 0.5*abs(det([1,1,1;vrtx']));
%   LblEg = TriMesh.elem2edge(ie,1:3);
%   signum = ones(1,3);
%   signum(find(ie==TriMesh.edge2elem(LblEg,2))) = -1;
%   LenEg(1) = norm(vrtx(2,:)-vrtx(3,:),2);
%   LenEg(2) = norm(vrtx(3,:)-vrtx(1,:),2);
%   LenEg(3) = norm(vrtx(1,:)-vrtx(2,:),2);
%   coeff = sln(LblEg);
%   NumerVel.EmCntr(ie,:) = ...
%     coeff(1)*signum(1)*LenEg(1)/(2*area)*(cntr-vrtx(1,:))+...
%     coeff(2)*signum(2)*LenEg(2)/(2*area)*(cntr-vrtx(2,:))+...
%     coeff(3)*signum(3)*LenEg(3)/(2*area)*(cntr-vrtx(3,:));
% end

%% Computing numerical flux: For all elements, 3 edges 
NumerFlux = zeros(NumEms,3);
coeff = sln(1:NumEgs);
for j=1:3
  NumerFlux(:,j) = coeff(e2g(:,j)) .* sn(:,j) .* LenEg(e2g(:,j));
end

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

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

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

return;