function [NumerPresCofP0P0, NumerVelCofAC0] ... 
  = Darcy_WG_QuadriP0P0AC0_PresVelFlux(QuadriMesh,PermK,sln,GAUSSQUAD) 
% James Liu, ColoState; 2017/07--2021/04 

%% Mesh info 
NumEms = QuadriMesh.NumEms;

%% For numerical pressure 
NumerPresCofP0P0 = sln;

%% 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 = zeros(NumEms,4,4);
% GMK = zeros(NumEms,4,4);
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);  % j-th fxn. projected onto k-th fxn. 
for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
  EltGM = squeeze(GM(ie,:,:));
  for j=1:4
    EltRHS = squeeze(GMK(ie,:,j));
    ProjCof(ie,j,:) = EltGM \ EltRHS';
  end
end
% Numerical velocity: Coeffs. in AC0 elementwise 
NumerVelCofAC0 = zeros(NumEms,4);
for k=1:4
  NumerVelCofAC0(:,k) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);
end

return;