function [NumerPresEm,NumerVelCofBDM1] = Darcy_WG_TriP0P1BDM1s_PresVelFlux(...
  TriMesh,PermKs,sln,GAUSSQUAD)
%% Darcy: WG(P0,P1;BDM1): Computing numerical pressure, velocity, etc. 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag>=2
% area = TriMesh.area;
% LenEg = TriMesh.LenEg;

%% Retrieving solution info: Part 1 
NumerPresEm = sln(1:NumEms);
NumerPresEg = reshape(sln(NumEms+1:end),2,NumEgs)';

%% Generating and adjusting local coefficients for WG basis functions 
coeff = zeros(NumEms,7);
for j=1:3
  coeff(:,2*j-1) = NumerPresEg(TriMesh.elem2edge(:,j),1);
  coeff(:,2*j)   = NumerPresEg(TriMesh.elem2edge(:,j),2);  
end
% JL20170209: TO BE REVISED FOR EFFECIENCY 
for ie=1:NumEms
  for j=1:3
    if TriMesh.SignEmEg(ie,j)==-1 
      tmp = coeff(ie,2*j-1);
      coeff(ie,2*j-1) = coeff(ie,2*j);
      coeff(ie,2*j)   = tmp;
    end
  end
end
coeff(:,7) = NumerPresEm(:);

%% Computing coeff. of disc.wk.grad. in BDM1 eg.bas.fxns. of WG bas.fxns.
CDWGB = WG_TriP0P1BDM1_CofBDM1EgBas_DiscWkGradBasFxn(TriMesh);

%% Retrieving solution info: Part 2 
% "abuse" coeffs. for the discrete weak gradients 
NumerVelCofBDM1 = zeros(NumEms,6);
for j=1:6
  NumerVelCofBDM1(:,j) = sum(coeff(:,1:7).*squeeze(CDWGB(:,j,1:7)),2);
end
% now coeffs. for numerical velocity 
NumerVelCofBDM1 = -repmat(PermKs,1,6).*NumerVelCofBDM1;

%% JL20170209: GAUSSQUAD TO BE USED LATER ON FOR L2 PROJECTION 

return;