function [PresAvgEm,NumerVelEmCntr,NumerFlux,NmlFlxEg,CofRT0]...
  = Darcy_DG_TriP1_PostProc(BndDescMat,TriMesh,PermK,CofNS)
%% Darcy: DG: TriP1: Postprocessing to get (RT0,P0) approx. 
% Assuming PermK is an elementwise scalar constant 
% James Liu, ColoState; 2012/07--2016/06 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag=3 
g2e = TriMesh.edge2elem;
e2g = TriMesh.elem2edge;

%% Computing pressure averages on elements: to get P0 
PresAvgEm = TriMesh_ConvDP1toP0(TriMesh,CofNS);

%% Computing numerical velocity at element centers 
% NumerVelEmCntr = zeros(NumEms,2);
NumerVelEmCntr = -[PermK,PermK].*CofNS(:,2:3);

%% Computing average numerical velocity on all edges 
NumerVelEgc = zeros(NumEgs,2);
% JL20131102: TO BE REVISED FOR EFFICIENCY 
for ig=1:NumEgs
  if (TriMesh.BndryEdge(ig)==0) 
    NumerVelEgc(ig,:) = 0.5*(NumerVelEmCntr(g2e(ig,1),:)...
      + NumerVelEmCntr(g2e(ig,2),:));
  end
  if (TriMesh.BndryEdge(ig)>0) 
    NumerVelEgc(ig,:) = NumerVelEmCntr(g2e(ig,1),:);
  end
end

%% Computing the normal "flux" (just the coefficients) over all edges 
% NmlFlxEg = zeros(NumEgs,1);
NmlFlxEg = dot(NumerVelEgc,TriMesh.NmlEg,2);
% element-edge normal flux 
NumerFlux.EmEgc = zeros(NumEms,3);
for j=1:3
  NumerFlux.EmEgc(:,j) = NmlFlxEg(e2g(:,j)).*TriMesh.SignEmEg(:,j);
end
% JL20131117: TO BE REVISED FOR EFFICIENCY 
NumBndryPcs = size(BndDescMat,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)*TriMesh.LenEg(ig);
  end
end

%% Computing coeffs. in RT0 edge-based bas.fxns.(local)
CofRT0 = zeros(NumEms,3);
for j=1:3
  CofRT0(:,j) = NmlFlxEg(e2g(:,j)).*TriMesh.SignEmEg(:,j);
end

return;