function [NumerPresEmAvg, NumerVelEmCntr] = ... 
  DarcyGen_WG_QuadriP0P0AC0_PresVelFlux( ...
    EqnBC, QuadriMesh, PermK, sln, CofQhg, GAUSSQUAD) 
% James Liu, ColoState; 2014/09--2020/09 

%% Mesh info 
NumEms = QuadriMesh.NumEms;

%% For numerical pressure 
NumerPresEmAvg = sln(1:NumEms);

%% For numerical velocity elementwise averages 
% Elementwise 5 coeffs. in WG bas.fxns. 
cof = zeros(NumEms,5);
cof(:,1) = sln(1:NumEms);
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  cof(:,1+j) = sln(NumEms+k);
end
% Elementwise cofs. in AC0 nmlz.Piola.bas.: Disc.wk.grad. of WG0 bas.fxns.
% CDWGB = zeros(NumEms,5,4) 
CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
CDWGS = zeros(NumEms,4);
for j=1:4
  CDWGS(:,j) = sum(cof.*squeeze(CDWGB(:,:,j)),2);
end
GM  = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat( QuadriMesh,GAUSSQUAD);
GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh,PermK,GAUSSQUAD);
% Projection coeffs. in AC0 nmlz.Piola.bas. for permeability * bas.fxn.
ProjCof = zeros(NumEms,4,4);  % k-th fxn. projected onto l-th fxn. 
for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
  EltGM = squeeze(GM(ie,:,:));
  for k=1:4
    EltRHS = -squeeze(GMK(ie,:,k))';
    ProjCof(ie,k,:) = EltGM\EltRHS;
  end
end
% Numerical velocity: Coeffs. in AC0 elementwise 
NumerVelCofAC0 = zeros(NumEms,4);
for l=1:4
  NumerVelCofAC0(:,l) = sum(CDWGS.*squeeze(ProjCof(:,:,l)),2);
end
NumerVelCofAC0 = NumerVelCofAC0 - CofQhg;

%% Velcoity 
NumerVelEmCntr = zeros(NumEms,2);
xhat = 0.5;  yhat = 0.5;
FV = Hdiv_QuadriAC0_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat);
NumerVelEmCntr(:,1) = NumerVelCofAC0(:,1) + NumerVelCofAC0(:,3).*FV(:,3,1) + NumerVelCofAC0(:,4).*FV(:,4,1);
NumerVelEmCntr(:,2) = NumerVelCofAC0(:,2) + NumerVelCofAC0(:,3).*FV(:,3,2) + NumerVelCofAC0(:,4).*FV(:,4,2);

return;