function [L2ErrPres] = Darcy_MFEM_RectRT0Q0_Err(EqnBC,RectMesh,sln,GAUSSQUAD)
%% Darcy: MFEM(RT[0],Q0) on a rectangular mesh: 
% Computing errors in pressure, etc. 
% JL20170211: TO BE FINISHED/REVISED 
% James Liu, ColoState; 2012/07--2017/02 

%% 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);
area = (x2-x1).*(y2-y1);

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

%% All kinds of errors 
err = zeros(1,6);
% err(1): L^2-norm of numerical vs exact velocity over all elements 
% err(2): L^\infty-norm of numerical vs exact pressure over all elements 
% err(5): L^2-norm of numerical vs exact pressure over all elements 
% err(6): L^\infty-norm of numerical vs exact pressure over all elements 

%% Computing L^2-norm of numerical vs exact velocity over all elements 
% Exact velocity computed at quadrature points 
% err(1) = 0;

%% Computing L^2-norm of error in pressure 
% Applying a Gaussian quadrature on rectangular elements 
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEms,2);
tmp = zeros(NumEms,1);
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;
  DiffP = EqnBC.fxnp(qp) - NumerPresEm;
  tmp = tmp + GAUSSQUAD.RECT(k,3)*(DiffP.^2);
end
err(5) = sqrt(sum(tmp.*area));
L2ErrPres = err(5);

return;