function [L2ErrPres,L2ErrVel,L2ErrFlux,L2ErrDiv] ... 
  = Darcy_WG_QuadriP0P0AC0_Err(fxp, fxnu, ...
    QuadriMesh,NumerPresCofP0,NumerVelCofAC0,GAUSSQUAD)
% James Liu, ColoState; 2014/09--2021/01 

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

%% For numerical (interior) pressure and velocity 
EmErrPres = zeros(NumEms,1);
EmErrVel  = zeros(NumEms,1);
NumerPresEm = NumerPresCofP0(1:NumEms);
MatJ = zeros(NumEms,2,2);
MatP = zeros(NumEms,2,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);
  J11 = CofA(:,2) + CofA(:,4)*yhat;  
  J12 = CofA(:,3) + CofA(:,4)*xhat;
  J21 = CofB(:,2) + CofB(:,4)*yhat;  
  J22 = CofB(:,3) + CofB(:,4)*xhat;
  jac = J11.*J22 - J12.*J21;
  MatJ(:,1,1) = J11;  MatJ(:,1,2) = J12;
  MatJ(:,2,1) = J21;  MatJ(:,2,2) = J22;
  for i=1:2
    for j=1:2
      MatP(:,i,j) = MatJ(:,i,j)./jac;
    end
  end
  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];
  X = qp(:,1) - QuadriMesh.EmCntr(:,1);
  Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
  U = MatP(:,1,1) .* xhat - MatP(:,1,2) .* yhat;
  V = MatP(:,2,1) .* xhat - MatP(:,2,2) .* yhat;
  % For pressure 
  pval = fxnp(qp);
  perr = pval - NumerPresEm;
  EmErrPres = EmErrPres + (perr.^2) .* jac * GAUSSQUAD.RECT(k,3);
  % For velocity 
  uval = fxnu(qp);
  NumerVel = [NumerVelCofAC0(:,1) + NumerVelCofAC0(:,3).*X + NumerVelCofAC0(:,4).*U, ... 
              NumerVelCofAC0(:,2) + NumerVelCofAC0(:,3).*Y + NumerVelCofAC0(:,4).*V];
  uerr = uval - NumerVel;
  EmErrVel = EmErrVel + (uerr(:,1).^2+uerr(:,2).^2) .* jac * GAUSSQUAD.RECT(k,3);
end
L2ErrPres = sqrt(sum(EmErrPres));
L2ErrVel  = sqrt(sum(EmErrVel));

%% For numerical flux 
L2ErrFlux = 0;

%% For numerical divergence  
L2ErrDiv = 0;

return;