function [NumerPresEm,NumerVelEmCntr,NumerFlux,FluxDscp,LMCR]...
  = Darcy_CG_RectQ1_PresVelFlux(EqnBC,RectMesh,PermK,sln,GAUSSQUAD)
%% Darcy: CG: RectQ1: Computing numer.pres., vel., flux dscp., LMCR
% Computing numerical pressure as elementwise constants,
%           numerical velocity at all element centers,
%           flux discrepancy across all interior edges,
%           local-mass-conservation-residual (LMCR) 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info
NumEms = RectMesh.NumEms;
NumEgs = RectMesh.NumEgs;
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);
Deltax = x2 - x1;
Deltay = y2 - y1;
area = Deltax .* Deltay;

%% Computing elementwise const numerical pressure 
% by averaging the nodal values
NumerPresEm = 0.25*sum(sln(RectMesh.elem),2);  % 0.25=1/4

%% Computing numerical velocity at element centers 
% See notes for derivation 
cof = sln(RectMesh.elem);
A1 = -cof(:,1) + cof(:,2) + cof(:,3) - cof(:,4);
A2 = -cof(:,1) - cof(:,2) + cof(:,3) + cof(:,4);
NumerVelEmCntr = zeros(NumEms,2);
NumerVelEmCntr(:,1) = -0.5*(PermK(:,1)./Deltax.*A1...
                           +PermK(:,2)./Deltay.*A2);
NumerVelEmCntr(:,2) = -0.5*(PermK(:,3)./Deltax.*A1...
                           +PermK(:,4)./Deltay.*A2);

%% JL2016013: TO BE FINISHED 
%% Computing elementwise outer-normal-fluxes on four(4) edges 
% Normal fluxes of pseudo basis: See notes for derivation (without -)
hxy = Deltax./Deltay;
hyx = Deltay./Deltax;
%
NmlFluxPsudBas = zeros(NumEms,4,4);
%
NmlFluxPsudBas(:,1,1) =  PermK(:,2,1) + 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,1,2) =                - 0.5*PermK(:,1,1).*hyx;
NmlFluxPsudBas(:,1,3) =                - 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,1,4) =  PermK(:,1,2) + 0.5*PermK(:,1,1).*hyx;
% 
NmlFluxPsudBas(:,2,1) = -PermK(:,2,1) + 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,2,2) = -PermK(:,1,2) + 0.5*PermK(:,1,1).*hyx;
NmlFluxPsudBas(:,2,3) =                - 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,2,4) =                - 0.5*PermK(:,1,1).*hyx;
% 
NmlFluxPsudBas(:,3,1) =                - 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,3,2) =  PermK(:,1,2) + 0.5*PermK(:,1,1).*hyx;
NmlFluxPsudBas(:,3,3) =  PermK(:,2,1) + 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,3,4) =                - 0.5*PermK(:,1,1).*hyx;
% 
NmlFluxPsudBas(:,4,1) =                - 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,4,2) =                - 0.5*PermK(:,1,1).*hyx;
NmlFluxPsudBas(:,4,3) = -PermK(:,2,1) + 0.5*PermK(:,2,2).*hxy;
NmlFluxPsudBas(:,4,4) = -PermK(:,1,2) + 0.5*PermK(:,1,1).*hyx;
% Now the outer-normal fluxes of the numerical solution 
NumerFlux = zeros(NumEms,4);
for j=1:4  % The j-th edge 
  NumerFlux(:,j) = -sum(cof.*squeeze(NmlFluxPsudBas(:,:,j)),2);  % Negative
end

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

%% JL20150722: TO BE REVISED FOR EFFECIENCY !!
%% Computing flux discrepancy across edges 
FluxDscp = zeros(NumEgs,1);
for ig=1:NumEgs
  if RectMesh.WhichEdge(ig,2)>0 
    % [ig, RectMesh.edge2elem(ig,1), RectMesh.edge2elem(ig,2)]
    FluxDscp(ig) = ...
      NumerFlux(RectMesh.edge2elem(ig,1),RectMesh.WhichEdge(ig,1)) + ...
      NumerFlux(RectMesh.edge2elem(ig,2),RectMesh.WhichEdge(ig,2));      
  end
end

%% Computing elementwise local-mass-conservation-residual (LMCR)
% First computing elementwise integrals of the source term 
IntgrlSrc = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
for k=1:NumQuadPts
  qp(:,1) = GAUSSQUAD.RECT(k,1)*x1 + GAUSSQUADRECT1(k,1)*x2;
  qp(:,2) = GAUSSQUAD.RECT(k,2)*y1 + GAUSSQUADRECT1(k,2)*y2;
  IntgrlSrc = IntgrlSrc + GAUSSQUAD.RECT(k,3) * EqnBC.fxnf(qp);
end
IntgrlSrc = IntgrlSrc .* area;
% Now comuting LMCR 
LMCR = IntgrlSrc - sum(NumerFlux,2);

return;