function [L2ErrPres,L2ErrVel,L2ErrFlux] = Darcy_WG_RectQ0Q0RT0_Err(...
  EqnBC,RectMesh,NumerPresEm,NumerVelCofRT0,GAUSSQUAD)
%% Darcy: WG(Q0,Q0;RT[0]) on a rectangular mesh: Computing errors 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = RectMesh.NumEms;
% Assuming TriMesh.flag>=2
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);

%% Computing L2ErrPres: L^2-norm of error in pressure (for element interior) 
ErrPresEm = 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;
  pval = EqnBC.fxnp(qp);
  perr = pval - NumerPresEm;
  ErrPresEm = ErrPresEm + GAUSSQUAD.RECT(k,3)*(perr.^2);
end
L2ErrPres = sqrt(sum(ErrPresEm.*RectMesh.area));

%% Computing L2-norm of error in velocity 
ErrVelEm = 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;
  X = qp(:,1) - RectMesh.EmCntr(:,1);
  Y = qp(:,2) - RectMesh.EmCntr(:,2);
  NV = [ NumerVelCofRT0(:,1) + NumerVelCofRT0(:,3).*X ,...
         NumerVelCofRT0(:,2) + NumerVelCofRT0(:,4).*Y ];  
  uval = EqnBC.fxnu(qp);
  verr = uval - NV;
  ErrVelEm = ErrVelEm + GAUSSQUAD.RECT(k,3)*(verr(:,1).^2 + verr(:,2).^2);
end
L2ErrVel = sqrt(sum(ErrVelEm.*RectMesh.area));

%% Computing L2-norm of error in flux  
L2ErrFlux = 0;

return;