function err = Darcy_CG_RectQ1_Err(EqnBC,RectMesh,PermK,sln,GAUSSQUAD)
%% Darcy: CG: Q1 on a rectangular mesh: Computing all errors 
% James Liu, ColoState; 2012/07--2017/02 

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

%% All kinds of error 
err = zeros(1,6);
% err(5): L^2-norm of numer.pres. vs exact pres. over all elements 

%% Computing errors against the exact/true pressure 
cof = sln(RectMesh.elem);
BasFxnVal = zeros(NumEms,4);
tmp = 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)-x1)./Deltax;  Y = (qp(:,2)-y1)./Deltay;
  BasFxnVal(:,1) = (1-X).*(1-Y);
  BasFxnVal(:,2) =     X.*(1-Y);
  BasFxnVal(:,3) =     X.*Y;
  BasFxnVal(:,4) = (1-X).*Y;
  NumerPres = sum(cof.*BasFxnVal,2);
  DiffP = abs(EqnBC.fxnp(qp)-NumerPres);
  tmp = tmp + GAUSSQUAD.RECT(k,3)*(DiffP.^2);
end
err(5) = sqrt(sum(tmp.*area));

return;