function [NumerVelEmAvg,NumerDivEmAvg,NumerPresEmAvg,NumerVelNd] = ... 
  Stokes_CG_QuadriBR1Q0_Presn(QuadriMesh,sln,GAUSSQUAD)
%% Stokes_CG_QuadriBR1Q0_Prsn: Computing quantities for presentation 
% James Liu, ColoState; 2012/07--2020/08 

%% Equation info 
% nu = EqnBC.nu;

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

%% Auxiliary quantities 
% Elementwise cofficients in 9 BR1 basis functions for numerical solution 
cof = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = QuadriMesh.elem2edge(:,j);
  cof(:,8+j) = sln(2*NumNds+k) .* QuadriMesh.EmEgSign(:,j);
end

%% Computing numerical velocity elementwise averages 
NumerVelEmAvg = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  NumerVelEmAvg(:,1) = NumerVelEmAvg(:,1) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area .* sum(cof.*squeeze(BFVal(:,:,1)),2);
  NumerVelEmAvg(:,2) = NumerVelEmAvg(:,2) + GAUSSQUAD.RECT(k,3) ...
    * jac ./area .* sum(cof.*squeeze(BFVal(:,:,2)),2);
end

%% Computing numerical divergence elementwise averages 
% Using Gaussian quadrature 
NumerDivEmAvg = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  NumerDivEmAvg = NumerDivEmAvg + GAUSSQUAD.RECT(k,3) ... 
    * jac ./ area .* sum(cof.*BFDiv,2);
end

%% "Computing" numerical pressure elementwise averages 
NumerPresEmAvg = sln(2*NumNds+NumEgs+(1:NumEms));
% % Further treatment for satisfying L0 condition 
% avg = sum(NumerPresEmAvg)/NumEms;
% NumerPresEmAvg = NumerPresEmAvg - avg;

%% "Computing"/extracting numerical velocity at all nodes 
NumerVelNd = zeros(NumNds,2);
NumerVelNd(:,1) = sln(2*(1:NumNds)-1);
NumerVelNd(:,2) = sln(2*(1:NumNds)  );

return;