function [NumerPresEmAvg, NumerVelEmAvg, NumerVelNd] = ... 
  DarcyGen_WG_QuadriP0P0AC0_Presn( ...
    EqnBC, QuadriMesh, PermK, sln, CofQhg, GAUSSQUAD) 
%% Darcy flow: Computing quantities for presentation on a quadri.mesh 
% James Liu, ColoState; 2014/09--2020/08 

%% Mesh info 
NumNds = QuadriMesh.NumNds;
NumEms = QuadriMesh.NumEms;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% 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 = 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 
CofAC0 = zeros(NumEms,4);
for k=1:4
  CofAC0(:,k) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);
end
CofAC0 = CofAC0 - CofQhg;

%% 
NumerVelEmAvg = CofAC0(:,1:2);

%% Computing nodal velocity after auxiliary elementwise vertexwise velocity 
NumerVelNd = zeros(NumNds,2);
cnt = zeros(NumNds,1);
xhat = zeros(4,1);  xhat(1) = 0;  xhat(2) = 1;  xhat(3) = 1;  xhat(4) = 0; 
yhat = zeros(4,1);  yhat(1) = 0;  yhat(2) = 0;  yhat(3) = 1;  yhat(4) = 1;   
EltMatJ = zeros(2,2);
EltMatP = zeros(2,2);
for ie=1:NumEms
  for j=1:4
    k = QuadriMesh.elem(ie,j);
    cnt(k) = cnt(k) + 1;
    x = QuadriMesh.node(k,1);
    y = QuadriMesh.node(k,2);
    X = x - QuadriMesh.EmCntr(ie,1);
    Y = y - QuadriMesh.EmCntr(ie,2);
    J11 = CofA(ie,2) + CofA(ie,4)*yhat(j);  
    J12 = CofA(ie,3) + CofA(ie,4)*xhat(j);
    J21 = CofB(ie,2) + CofB(ie,4)*yhat(j);  
    J22 = CofB(ie,3) + CofB(ie,4)*xhat(j);
    jac = J11.*J22 - J12.*J21;
    EltMatJ(1,1) = J11;  EltMatJ(1,2) = J12;
    EltMatJ(2,1) = J21;  EltMatJ(2,2) = J22;
    for ii=1:2
      for jj=1:2
        EltMatP(ii,jj) = EltMatJ(ii,jj)./jac;
      end
    end
    U = EltMatP(1,1) .* xhat(j) - EltMatP(1,2) .* yhat(j);
    V = EltMatP(2,1) .* xhat(j) - EltMatP(2,2) .* yhat(j);
    v1tmp = CofAC0(ie,1) + CofAC0(ie,3)*X + CofAC0(ie,4)*U; 
    v2tmp = CofAC0(ie,2) + CofAC0(ie,3)*Y + CofAC0(ie,4)*V; 
    NumerVelNd(k,1) = NumerVelNd(k,1) + v1tmp;
    NumerVelNd(k,2) = NumerVelNd(k,2) + v2tmp;
  end
end
NumerVelNd(:,1) = NumerVelNd(:,1)./cnt;
NumerVelNd(:,2) = NumerVelNd(:,2)./cnt;

return;