function [L2ErrPres,L2ErrVel,L2ErrFlux] = Darcy_WG_QuadriP1P0P02_Err(...
  EqnBC,QuadriMesh,PermK,sln,GAUSSQUAD)
%% Darcy: WG(P1,P0;P0^2)Quadri: 
% Cf.: MuWangYe_JCAM_2015 
% Computing errors in pressure, velocity, flux 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% JL20170212: TO BE FINISHED/REVISED 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = QuadriMesh.NumEms;
% NumEgs = QuadriMesh.NumEgs;
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x = zeros(NumEms,5);  y = zeros(NumEms,5);
x(:,1) = QuadriMesh.node(k1,1);  y(:,1) = QuadriMesh.node(k1,2);
x(:,2) = QuadriMesh.node(k2,1);  y(:,2) = QuadriMesh.node(k2,2);
x(:,3) = QuadriMesh.node(k3,1);  y(:,3) = QuadriMesh.node(k3,2);
x(:,4) = QuadriMesh.node(k4,1);  y(:,4) = QuadriMesh.node(k4,2);
x(:,5) = x(:,1);  y(:,5) = y(:,1);
xc = 0.25*sum(x(:,1:4),2);
yc = 0.25*sum(y(:,1:4),2);
% diag13 = sqrt((x(:,1)-x(:,3)).^2+(y(:,1)-y(:,3)).^2);
% diag24 = sqrt((x(:,2)-x(:,4)).^2+(y(:,2)-y(:,4)).^2);
% diam = max(diag13,diag24);
% dm1 = 1./diam;
CofA = zeros(NumEms,4);                       CofB = zeros(NumEms,4);
CofA(:,1) = x(:,1);                           CofB(:,1) = y(:,1);
CofA(:,2) = x(:,2)-x(:,1);                    CofB(:,2) = y(:,2)-y(:,1);
CofA(:,3) = x(:,4)-x(:,1);                    CofB(:,3) = y(:,4)-y(:,1);
CofA(:,4) = (x(:,1)+x(:,3))-(x(:,2)+x(:,4));  CofB(:,4) = (y(:,1)+y(:,3))-(y(:,2)+y(:,4));

%% Computing L2-norm of error in pressure (element interior)
CofNumerPres = zeros(NumEms,3);
CofNumerPres(:,3) = sln(3*(1:NumEms));
CofNumerPres(:,2) = sln(3*(1:NumEms)-1);
CofNumerPres(:,1) = sln(3*(1:NumEms)-2);
ErrPresEm = zeros(NumEms,1);
qp = 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);
  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;
  Y = qp(:,2) - yc;
  NumerPres = CofNumerPres(:,1) + CofNumerPres(:,2).*X + CofNumerPres(:,3).*Y;
  pval = EqnBC.fxnp(qp);
  perr = pval - NumerPres;
  ErrPresEm = ErrPresEm + GAUSSQUAD.RECT(k,3)*jac.*perr.^2;
end
ErrPresEm = sqrt(ErrPresEm);
L2ErrPres = norm(ErrPresEm,2);

%% JL20170212: TO BE INVESGATED: How is a numerical velocity defined?
%% Computing L2-norm of error in velocity 
L2ErrVel = 0;

%% JL20170212: TO BE FINISHED/REVISED 
%% Computing L2-norm of error in flux 
L2ErrFlux = 0;

return;