function [L2ErrPres,L2ErrVel,L2ErrFlux] = Darcy_WG_HybridPQ0RT0_Err(...
  EqnBC,HybridMesh,NumerPresEm,NumerVelCofRT0,GAUSSQUAD)
%% Darcy: Hybrid mesh: Computing errors 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = HybridMesh.NumEms;
NumEmsQ = HybridMesh.NumEmsQ;
NumEmsT = HybridMesh.NumEmsT;
LblQ = HybridMesh.labelemQ;
LblT = HybridMesh.labelemT;
xc = zeros(NumEms,1);
yc = zeros(NumEms,1);

%% Auxiliary quantities for quadrilateral elements
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2);
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
x4 = HybridMesh.node(k4,1);  y4 = HybridMesh.node(k4,2);
areaQ = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
            + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));
CofA = zeros(NumEmsQ,4);      CofB = zeros(NumEmsQ,4);
CofA(:,1) = x1;               CofB(:,1) = y1;
CofA(:,2) = x2-x1;            CofB(:,2) = y2-y1;
CofA(:,3) = x4-x1;            CofB(:,3) = y4-y1;
CofA(:,4) = (x1+x3)-(x2+x4);  CofB(:,4) = (y1+y3)-(y2+y4);
xc(LblQ) = 0.25*(x1+x2+x3+x4);
yc(LblQ) = 0.25*(y1+y2+y3+y4);

%% Auxiliary quantities for triangular elements
k1 = HybridMesh.elemT(:,1);  
k2 = HybridMesh.elemT(:,2);  
k3 = HybridMesh.elemT(:,3);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2); 
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
areaT = 0.5*((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1));
xc(LblT) = (x1+x2+x3)/3;
yc(LblT) = (y1+y2+y3)/3;

%% Computing L2-norm of error in pressure 
% For quadrilaterals 
NumerPresEmQ = NumerPresEm(LblQ);
ErrPresEmQtmp2 = zeros(NumEmsQ,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEmsQ,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 - NumerPresEmQ;
  ErrPresEmQtmp2 = ErrPresEmQtmp2 + GAUSSQUAD.RECT(k,3)*jac.*(perr.^2);
end
ErrPresEmQ2 = sum(ErrPresEmQtmp2);
% For triangles 
NumerPresEmT = NumerPresEm(LblT);
ErrPresEmTtmp2 = zeros(NumEmsT,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * HybridMesh.node(HybridMesh.elemT(:,1),:) ...
     + GAUSSQUAD.TRIG(k,2) * HybridMesh.node(HybridMesh.elemT(:,2),:) ...
     + GAUSSQUAD.TRIG(k,3) * HybridMesh.node(HybridMesh.elemT(:,3),:) ;
  perr = EqnBC.fxnp(qp) - NumerPresEmT;
  ErrPresEmTtmp2 = ErrPresEmTtmp2 + GAUSSQUAD.TRIG(k,4)*(perr.^2);
end
ErrPresEmT2 = sum(ErrPresEmTtmp2.*areaT);
% Now putting together
L2ErrPres = sqrt(ErrPresEmQ2+ErrPresEmT2);

%% Computing L2-norm of error in velocity 
% For quadrilaterals 
ErrVelEmQtmp2 = zeros(NumEmsQ,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
qp = zeros(NumEmsQ,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);
  X = qp(:,1) - xc(LblQ,1);
  Y = qp(:,2) - yc(LblQ);
  NV = [ NumerVelCofRT0.quadri(:,1) + NumerVelCofRT0.quadri(:,3).*X ,...
         NumerVelCofRT0.quadri(:,2) + NumerVelCofRT0.quadri(:,4).*Y ];
  uval = EqnBC.fxnu(qp);
  verr = uval - NV;
  ErrVelEmQtmp2 = ErrVelEmQtmp2 + GAUSSQUAD.RECT(k,3)*jac.*(verr(:,1).^2+verr(:,2).^2);
end
ErrVelEmQ2 = sum(ErrVelEmQtmp2);
% For triangles
ErrVelEmTtmp2 = zeros(NumEmsT,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1)*HybridMesh.node(HybridMesh.elemT(:,1),:)...
     + GAUSSQUAD.TRIG(k,2)*HybridMesh.node(HybridMesh.elemT(:,2),:)...
     + GAUSSQUAD.TRIG(k,3)*HybridMesh.node(HybridMesh.elemT(:,3),:);
  X = qp(:,1) - xc(LblT);
  Y = qp(:,2) - yc(LblT);
  NV = [ NumerVelCofRT0.tri(:,1) + NumerVelCofRT0.tri(:,3).*X ,...
         NumerVelCofRT0.tri(:,2) + NumerVelCofRT0.tri(:,3).*Y ];  
  uval = EqnBC.fxnu(qp);
  verr = uval - NV;
  ErrVelEmTtmp2 = ErrVelEmTtmp2 + GAUSSQUAD.TRIG(k,4)*(verr(:,1).^2 + verr(:,2).^2);
end
ErrVelEmT2 = sum(ErrVelEmTtmp2.*areaT);
% Now putting together
L2ErrVel = sqrt(ErrVelEmQ2+ErrVelEmT2);

%% Computing L2-norm of flux error 
% For quadrilaterals
ErrFluxEg.quadri = zeros(NumEmsQ,4);
for j=1:3
  x1 = HybridMesh.node(HybridMesh.elemQ(:,j),1);
  y1 = HybridMesh.node(HybridMesh.elemQ(:,j),2);
  x2 = HybridMesh.node(HybridMesh.elemQ(:,j+1),1);
  y2 = HybridMesh.node(HybridMesh.elemQ(:,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) - xc(LblQ);
    Y = qp(:,2) - yc(LblQ);
    NV = [ NumerVelCofRT0.quadri(:,1) + NumerVelCofRT0.quadri(:,3).*X ,...
           NumerVelCofRT0.quadri(:,2) + NumerVelCofRT0.quadri(:,4).*Y ];
    vval = EqnBC.fxnu(qp);
    verr = vval - NV;
    ferr = verr(:,1).*nml(:,1) + verr(:,2).*nml(:,2);
    ErrFluxEg.quadri(:,j) = ErrFluxEg.quadri(:,j) + GAUSSQUAD.LINE(k,3)*(ferr.^2);
  end
  ErrFluxEg.quadri(:,j) = ErrFluxEg.quadri(:,j) .* areaQ;
end
for j=4:4
  x1 = HybridMesh.node(HybridMesh.elemQ(:,4),1);
  y1 = HybridMesh.node(HybridMesh.elemQ(:,4),2);
  x2 = HybridMesh.node(HybridMesh.elemQ(:,1),1);
  y2 = HybridMesh.node(HybridMesh.elemQ(:,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) - xc(LblQ);
    Y = qp(:,2) - yc(LblQ);
    NV = [ NumerVelCofRT0.quadri(:,1) + NumerVelCofRT0.quadri(:,3).*X ,...
           NumerVelCofRT0.quadri(:,2) + NumerVelCofRT0.quadri(:,4).*Y ];
    vval = EqnBC.fxnu(qp);
    verr = vval - NV;
    ferr = verr(:,1).*nml(:,1) + verr(:,2).*nml(:,2);
    ErrFluxEg.quadri(:,j) = ErrFluxEg.quadri(:,j) + GAUSSQUAD.LINE(k,3)*(ferr.^2);
  end
  ErrFluxEg.quadri(:,j) = ErrFluxEg.quadri(:,j) .* areaQ;
end
% For triangles 
ErrFluxEg.tri = zeros(NumEmsT,3);
for j=1:2
  x1 = HybridMesh.node(HybridMesh.elemT(:,j),1);
  y1 = HybridMesh.node(HybridMesh.elemT(:,j),2);
  x2 = HybridMesh.node(HybridMesh.elemT(:,j+1),1);
  y2 = HybridMesh.node(HybridMesh.elemT(:,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) - xc(LblT);
    Y = qp(:,2) - yc(LblT);
    NV = [ NumerVelCofRT0.tri(:,1) + NumerVelCofRT0.tri(:,3).*X ,...
           NumerVelCofRT0.tri(:,2) + NumerVelCofRT0.tri(:,3).*Y ];
    vval = EqnBC.fxnu(qp);
    verr = vval - NV;
    ferr = verr(:,1).*nml(:,1) + verr(:,2).*nml(:,2);
    ErrFluxEg.tri(:,j) = ErrFluxEg.tri(:,j) + GAUSSQUAD.LINE(k,3)*(ferr.^2);
  end
  ErrFluxEg.tri(:,j) = ErrFluxEg.tri(:,j) .* areaT;
end
for j=3:3
  x1 = HybridMesh.node(HybridMesh.elemT(:,3),1);
  y1 = HybridMesh.node(HybridMesh.elemT(:,3),2);
  x2 = HybridMesh.node(HybridMesh.elemT(:,1),1);
  y2 = HybridMesh.node(HybridMesh.elemT(:,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) - xc(LblT);
    Y = qp(:,2) - yc(LblT);
    NV = [ NumerVelCofRT0.tri(:,1) + NumerVelCofRT0.tri(:,3).*X ,...
           NumerVelCofRT0.tri(:,2) + NumerVelCofRT0.tri(:,3).*Y ];
    vval = EqnBC.fxnu(qp);
    verr = vval - NV;
    ferr = verr(:,1).*nml(:,1) + verr(:,2).*nml(:,2);
    ErrFluxEg.tri(:,j) = ErrFluxEg.tri(:,j) + GAUSSQUAD.LINE(k,3)*(ferr.^2);
  end
  ErrFluxEg.tri(:,j) = ErrFluxEg.tri(:,j) .* areaT;
end
% Now putting together 
L2ErrFlux = sqrt(sum(sum(ErrFluxEg.quadri,2))+sum(sum(ErrFluxEg.tri,2)));

return