function [L2ErrVel,L2ErrPres] = ...
  Stokes_CG_QuadriBR1Q0_Err(EqnBC,QuadriMesh,sln,GAUSSQUAD)
%% Stokes_CG_QuadriBR1Q0_Err: Computing errors  
% James Liu, ColoState; 2012/07--2018/02 

%% Equation info 
% mu = EqnBC.mu;

%% Mesh info 
NumNds = QuadriMesh.NumNds;
NumEgs = QuadriMesh.NumEgs;
NumEms = QuadriMesh.NumEms;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Auxiliary quantities 
% Elementwise cofficients in 9 BR1 basis functions for numerical solution 
cof = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = QuadriMesh.elem2edge(:,j);
  cof(:,8+j) = sln(2*NumNds+k) .* QuadriMesh.EmEgSign(:,j);
end

%% Computing L2-norm of error in velocity 
ErrVelEm = zeros(NumEms,1);
NumerVel = zeros(NumEms,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  u = EqnBC.fxnu(qp);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  NumerVel(:,1) = sum(cof.*squeeze(BFVal(:,:,1)),2);
  NumerVel(:,2) = sum(cof.*squeeze(BFVal(:,:,2)),2);
  uerr = u - NumerVel;
  ErrVelEm = ErrVelEm + (uerr(:,1).^2+uerr(:,2).^2)...
    .* jac * GAUSSQUAD.RECT(k,3);
end
ErrVelEm = sqrt(ErrVelEm);
L2ErrVel = norm(ErrVelEm,2);

%% Computing L2-norm of error in pressure  
ErrPresEm = zeros(NumEms,1);
NumerPres = sln(2*NumNds+NumEgs+(1:NumEms));
avg = sum(NumerPres)/NumEms;
NumerPres = NumerPres - avg;
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  p = EqnBC.fxnp(qp);
  perr = p - NumerPres;
  ErrPresEm = ErrPresEm + perr(:).^2 .* jac * GAUSSQUAD.RECT(k,3);
end
ErrPresEm = sqrt(ErrPresEm);
L2ErrPres = norm(ErrPresEm,2);

return;