function [PresP0,CofBDM1,NumerFlux] = Darcy_DG_TriP2_PostProc(...
  BndDescMat,TriMesh,PermK,CofNS,GAUSSQUAD)
%% Darcy: DG: TriP2: Postprocessing 
% Assuming PermK is an elementwise scalar 
% James Liu, ColoState; 2012/07--2016/06 

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

%% Computing pressure averages over all elements 
PresP0 = zeros(NumEms,1);
% BasFxn = zeros(NumEms,6);
WO = ones(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),:);
  X = qp(:,1);  Y = qp(:,2);
  BasFxn = [WO, X, Y, X.^2, X.*Y, Y.^2];
  PresP0 = PresP0 + GAUSSQUAD.TRIG(k,4) * dot(CofNS,BasFxn,2);  % */area
end

%% Elementwise coeffs. for numerical velocity in nat.bas. (1,x,y) 
% First as coeffs. for grad. 
CofNumerVel = zeros(NumEms,3,2);
CofNumerVel(:,:,1) = [CofNS(:,2), 2*CofNS(:,4), CofNS(:,5)];
CofNumerVel(:,:,2) = [CofNS(:,3), CofNS(:,5), 2*CofNS(:,6)];
% Now as coeffs. for numer.vel.
CofNumerVel(:,:,1) = -[PermK,PermK,PermK].*CofNumerVel(:,:,1);
CofNumerVel(:,:,2) = -[PermK,PermK,PermK].*CofNumerVel(:,:,2);

%% JL20131109: TO BE REVISED FOR EFFICIENCY 
%% Edgewise averaged coeffs. for numerical velocity in nat.bas.(1,x,y) 
EgCofNV = zeros(NumEgs,3,2);
for ig=1:NumEgs
  if (TriMesh.BndryEdge(ig)==0) 
    EgCofNV(ig,:,:) = 0.5*(CofNumerVel(g2e(ig,1),:,:)...
      + CofNumerVel(g2e(ig,2),:,:));
  end
  if (TriMesh.BndryEdge(ig)>0) 
    EgCofNV(ig,:,:) = CofNumerVel(g2e(ig,1),:,:);
  end
end

%% Auxiliary quantities from the mesh 
xcrd = zeros(NumEms,3);  ycrd = zeros(NumEms,3);
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
xcrd(:,1) = TriMesh.node(k1,1);  ycrd(:,1) = TriMesh.node(k1,2); 
xcrd(:,2) = TriMesh.node(k2,1);  ycrd(:,2) = TriMesh.node(k2,2);
xcrd(:,3) = TriMesh.node(k3,1);  ycrd(:,3) = TriMesh.node(k3,2);
eg1 = [xcrd(:,3),ycrd(:,3)] - [xcrd(:,2),ycrd(:,2)];
eg2 = [xcrd(:,1),ycrd(:,1)] - [xcrd(:,3),ycrd(:,3)];
eg3 = [xcrd(:,2),ycrd(:,2)] - [xcrd(:,1),ycrd(:,1)];
lg1 = sqrt(eg1(:,1).^2+eg1(:,2).^2);
lg2 = sqrt(eg2(:,1).^2+eg2(:,2).^2);
lg3 = sqrt(eg3(:,1).^2+eg3(:,2).^2);
tan = zeros(NumEms,3,2);
tan(:,1,:) = eg1./[lg1,lg1];
tan(:,2,:) = eg2./[lg2,lg2];
tan(:,3,:) = eg3./[lg3,lg3];
nml = zeros(NumEms,3,2);
nml(:,1,:) = [tan(:,1,2),-tan(:,1,1)];
nml(:,2,:) = [tan(:,2,2),-tan(:,2,1)];
nml(:,3,:) = [tan(:,3,2),-tan(:,3,1)];

%% Projecting(intplt) numer.vel.: coeffs. in BDM1 eg.bas.bas.fxns.(local)
CofBDM1 = zeros(NumEms,6);  % all elements, 3 edges, 2-vec. 
NumerVel = zeros(NumEms,2);  % working array 
WO = ones(NumEms,1);
j1s = [2,3,1];
j2e = [3,1,2];
for j=1:3
  % Edge-j starting point 
  js = j1s(j);
  NumerVel(:,1) = dot(CofNumerVel(:,:,1),[WO,xcrd(:,js),ycrd(:,js)],2);
  NumerVel(:,2) = dot(CofNumerVel(:,:,2),[WO,xcrd(:,js),ycrd(:,js)],2);
  CofBDM1(:,2*j-1) = dot(NumerVel,squeeze(nml(:,j,:)),2);
  % Edge-j ending point 
  je = j2e(j);
  NumerVel(:,1) = dot(CofNumerVel(:,:,1),[WO,xcrd(:,je),ycrd(:,je)],2);
  NumerVel(:,2) = dot(CofNumerVel(:,:,2),[WO,xcrd(:,je),ycrd(:,je)],2);
  CofBDM1(:,2*j) = dot(NumerVel,squeeze(nml(:,j,:)),2);  
end

%% JL20131120: TO BE REVISED FOR EFFICIENCY 
% element-edge normal flux 
NumerFlux.EmEgc = zeros(NumEms,3);
for j=1:3
  NumerFlux.EmEgc(:,j) = 0.5*(CofBDM1(:,2*j-1)+CofBDM1(:,2*j));
end
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

return;