function [NumerPresEm,NumerVelEmCntr,NumerFlux,LMCR,FluxDscp] = ...
  Darcy_MFEM_QuadriRT0Q0_PresVelFlux(QuadriMesh,sln)
%% Darcy: MFEM(RT[0],Q0) on a quadrilateral mesh: 
% Computing numerical pressure, velocity, flux, etc. 
% 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);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
% Coefficients for the bilinear mapping 
CofA = zeros(NumEms,4);       CofB = zeros(NumEms,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);

%% "Computing" (retrieving) numerical pressure 
NumerPresEm = sln(NumEgs+(1:NumEms));

%% Computing numerical velocity at element centers 
CofPiolaBasEm = zeros(NumEms,4);
for j=1:4
  CofPiolaBasEm(:,j) = sln(QuadriMesh.elem2edge(:,j));
end
xhat = 0.5;  yhat = 0.5;  % Center of the reference element [0,1]^2 
Phihat = zeros(NumEms,4,2);
WO = ones(NumEms,1);
Phihat(:,4,1) = (1-xhat)*WO;  Phihat(:,2,1) = xhat*WO;
Phihat(:,1,2) = (1-yhat)*WO;  Phihat(:,3,2) = yhat*WO;
J11 = CofA(:,2) + CofA(:,4)*yhat;  J12 = CofA(:,3) + CofA(:,4)*xhat;
J21 = CofB(:,2) + CofB(:,4)*yhat;  J22 = CofB(:,3) + CofB(:,4)*xhat;
jac = J11.*J22 - J12.*J21;
% Piola matrix = JacobianMatrix/JacobianDet 
MatP = zeros(NumEms,2,2);
MatP(:,1,1) = J11./jac;  MatP(:,1,2) = J12./jac;
MatP(:,2,1) = J21./jac;  MatP(:,2,2) = J22./jac;
NumerVelEmCntr = zeros(NumEms,2);
for j=1:4
  NumerVelEmCntr(:,1) = NumerVelEmCntr(:,1) + CofPiolaBasEm(:,j).*...
    squeeze(MatP(:,1,1).*Phihat(:,j,1) + MatP(:,1,2).*Phihat(:,j,2));
  NumerVelEmCntr(:,2) = NumerVelEmCntr(:,2) + CofPiolaBasEm(:,j).*...
    squeeze(MatP(:,2,1).*Phihat(:,j,1) + MatP(:,2,2).*Phihat(:,j,2));
end

%% JL20170211: TO BE FINISHED/REVISED 
%% Computing numerical flux 
NumerFlux = zeros(NumEms,4);

%% "Computing" Local-Mass-Conservation-Residual: 0 theoretically  
LMCR = zeros(NumEms,1);

%% "Computing" flux discrepancy across edges 0 theoretically  
FluxDscp = zeros(NumEgs,1);

return;