function [NumerPresEm, NumerVelEmCntr, NumerFlux, LMCR, FluxDscp] = ...
  Darcy_WG_QuadriP1P0P02_PresVelFlux(QuadriMesh, PermKt, sln)
%% Darcy: WG(P1,P0;P0^2)Quadri: 
% Computing (numerical) pressure, velocity, flux, etc. 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% JL20170212: TO BE FINISHED/REVISED 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
NumEgs = QuadriMesh.NumEgs;
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x = zeros(NumEms,5);  y = zeros(NumEms,5);
x(:,1) = QuadriMesh.node(k1,1);  y(:,1) = QuadriMesh.node(k1,2);
x(:,2) = QuadriMesh.node(k2,1);  y(:,2) = QuadriMesh.node(k2,2);
x(:,3) = QuadriMesh.node(k3,1);  y(:,3) = QuadriMesh.node(k3,2);
x(:,4) = QuadriMesh.node(k4,1);  y(:,4) = QuadriMesh.node(k4,2);
x(:,5) = x(:,1);  y(:,5) = y(:,1);
xc = 0.25*sum(x(:,1:4),2);
yc = 0.25*sum(y(:,1:4),2);
diag13 = sqrt((x(:,1)-x(:,3)).^2+(y(:,1)-y(:,3)).^2);
diag24 = sqrt((x(:,2)-x(:,4)).^2+(y(:,2)-y(:,4)).^2);
diam = max(diag13,diag24);
dm1 = 1./diam;

%% "Computing" the numerical pressure: Elementwise constants 
NumerPresEm = sln(3*(0:NumEms-1)+1);

%% Auxiliary quantities 
% Elementwise 7 coefficients in WG(P1,P0) basis functions 
CofWGB = zeros(NumEms,7);
CofWGB(:,1) = sln(3*(0:NumEms-1)+1);
CofWGB(:,2) = sln(3*(0:NumEms-1)+2);
CofWGB(:,3) = sln(3*(0:NumEms-1)+3);
for j=1:4
  CofWGB(:,3+j) = sln(3*NumEms+QuadriMesh.elem2edge(:,j));
end

%% Cofs. in P0^2 natural basis for disc.wk.grad. of WG(P1,P0)Q 7 bas.fxns.
% NOTE: CDWGB has size NumEms*7*2 
CDWGB = WG_QuadriP1P0P02_CofNatBas_DiscWkGradBasFxn(QuadriMesh);

%% Computing elementwise discrete weak gradient of the numerical solution 
DWGsln = zeros(NumEms,2);
DWGsln(:,1) = sum(CofWGB.*squeeze(CDWGB(:,:,1)),2);
DWGsln(:,2) = sum(CofWGB.*squeeze(CDWGB(:,:,2)),2);

%% Computing the numerical velocity at element centers 
NumerVelEmCntr = zeros(NumEms,2);
for i=1:2
  NumerVelEmCntr(:,i) = - ( squeeze(PermKt(:,i,1)).*DWGsln(:,1)...
                          + squeeze(PermKt(:,i,2)).*DWGsln(:,2));
end

%% "Computing" local-mass-conservation residuals for all elements 
LMCR = zeros(NumEms,1);

%% Computing elementwise 4 edge (total) outward normal fluxes 
NumerFlux = zeros(NumEms,4);
WO = ones(NumEms,1);
for j=1:4
  Xm = 0.5*(x(:,j)+x(:,j+1)) - xc;
  Ym = 0.5*(y(:,j)+y(:,j+1)) - yc;
  NmlEg = [y(:,j+1)-y(:,j), x(:,j)-x(:,j+1)];  % carrying length info 
  LenEg = sqrt((x(:,j)-x(:,j+1)).^2+(y(:,j)-y(:,j+1)).^2);
  UnitNmlEg = NmlEg./[LenEg,LenEg];
  NumerFlux(:,j) = sum(NumerVelEmCntr.*UnitNmlEg,2)...
                 + dm1 .* (sum(CofWGB(:,1:3).*[WO,Xm,Ym],2)-CofWGB(:,3+j));
end

%% Computing (normal) flux discrepancies across all edges 
FluxDscp = zeros(NumEgs,1);
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  tmp = FluxDscp(k);
  FluxDscp(k) = tmp - NumerFlux(:,j);
end
% Correction for boundary edges 
k = find(QuadriMesh.BndryEdge>0);
FluxDscp(k) = 0;

return;