function [L2ErrPres,L2ErrVel,L2ErrFlux] = Darcy_WG_QuadriQ0Q0RT0_Err(...
  EqnBC,QuadriMesh,NumerPresEm,NumerVelCofRT0,GAUSSQUAD)
%% Darcy: WG(Q0,Q0,RT[0])Quadri: Computing errors 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
% Assuming TriMesh.flag>=2
% area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% All kinds of error 
%     L2ErrPres: L^2-norm of error in pressure (for element interior) 
%      L2ErrVel: L^2-norm of error in velocity 
%     L2ErrFlux: L^2-norm of error in flux 
% DiscL2ErrPres: Discrete L^2-norm of error in pressure (for element interior) 

%% 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
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  pval = EqnBC.fxnp(qp);
  perr = pval - NumerPresEm;
  ErrPresEm = ErrPresEm + GAUSSQUAD.RECT(k,3)*jac.*(perr.^2);
end
L2ErrPres = sqrt(sum(ErrPresEm));

%% Computing L2ErrVel: L^2-norm of the error in velocity 
% Using NumerVelCofRT0 
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
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp(:,1) = CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat;
  qp(:,2) = CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat;
  X = qp(:,1) - QuadriMesh.EmCntr(:,1);
  Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
  NV = [ NumerVelCofRT0(:,1) + NumerVelCofRT0(:,3).*X ,...
         NumerVelCofRT0(:,2) + NumerVelCofRT0(:,4).*Y ];
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat) ...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat) ;
  vval = EqnBC.fxnu(qp);
  verr = vval - NV;
  ErrVelEm = ErrVelEm + GAUSSQUAD.RECT(k,3)*jac.*(verr(:,1).^2+verr(:,2).^2);
end
L2ErrVel = sqrt(sum(ErrVelEm));

%% Computing L2ErrFlux: L2-norm of error in flux 
ErrFluxEg = zeros(NumEms,4);
for j=1:3
  x1 = QuadriMesh.node(QuadriMesh.elem(:,j),1);
  y1 = QuadriMesh.node(QuadriMesh.elem(:,j),2);
  x2 = QuadriMesh.node(QuadriMesh.elem(:,j+1),1);
  y2 = QuadriMesh.node(QuadriMesh.elem(:,j+1),2);
  PseudoNml = [y2-y1,x1-x2];  % Normal vector carrying with edge length 
  LenEg = sqrt(PseudoNml(:,1).^2+PseudoNml(:,2).^2);
  nml = PseudoNml./[LenEg,LenEg];  % Now unit normal vector   
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    X = qp(:,1) - QuadriMesh.EmCntr(:,1);
    Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
    NV = [ NumerVelCofRT0(:,1) + NumerVelCofRT0(:,3).*X ,...
           NumerVelCofRT0(:,2) + NumerVelCofRT0(:,4).*Y ];
    vval = EqnBC.fxnu(qp);
    verr = vval - NV;
    ferr = verr(:,1).*nml(:,1) + verr(:,2).*nml(:,2);
    ErrFluxEg(:,j) = ErrFluxEg(:,j) + GAUSSQUAD.LINE(k,3)*(ferr.^2);
  end
  ErrFluxEg(:,j) = ErrFluxEg(:,j) .* QuadriMesh.area;
end
for j=4:4
  x1 = QuadriMesh.node(QuadriMesh.elem(:,4),1);
  y1 = QuadriMesh.node(QuadriMesh.elem(:,4),2);
  x2 = QuadriMesh.node(QuadriMesh.elem(:,1),1);
  y2 = QuadriMesh.node(QuadriMesh.elem(:,1),2);
  PseudoNml = [y2-y1,x1-x2];  % Normal vector carrying with edge length 
  LenEg = sqrt(PseudoNml(:,1).^2+PseudoNml(:,2).^2);
  nml = PseudoNml./[LenEg,LenEg];  % Now unit normal vector   
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    X = qp(:,1) - QuadriMesh.EmCntr(:,1);
    Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
    NV = [ NumerVelCofRT0(:,1) + NumerVelCofRT0(:,3).*X ,...
           NumerVelCofRT0(:,2) + NumerVelCofRT0(:,4).*Y ];
    vval = EqnBC.fxnu(qp);
    verr = vval - NV;
    ferr = verr(:,1).*nml(:,1) + verr(:,2).*nml(:,2);
    ErrFluxEg(:,j) = ErrFluxEg(:,j) + GAUSSQUAD.LINE(k,3)*(ferr.^2);
  end
  ErrFluxEg(:,j) = ErrFluxEg(:,j) .* QuadriMesh.area;
end
L2ErrFlux = sqrt(sum(sum(ErrFluxEg,2)));

%% Computing DiscL2ErrPres: Discrete L^2-norm of error in pressure 
pval = EqnBC.fxnp(QuadriMesh.EmCntr);
perr = pval - NumerPresEm;
DiscL2ErrPres = sqrt(sum((perr.^2).*QuadriMesh.area));

return;