function [L2ErrPres,L2ErrVel,L2ErrFlux] = Darcy_MFEM_QuadriRT0Q0_Err(...
  EqnBC,QuadriMesh,sln,GAUSSQUAD)
%% Darcy: MFEM(RT[0],Q0) on a quadrilateral mesh: Computing errors 
% Cf. Arnold,Boffi,Falk, SINUM, 2005 
% 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);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
% Coefficients for the bilinear mapping 
CofA = zeros(NumEms,4);       CofB = zeros(NumEms,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);

%% Computing the L2-norm of the error in pressure  
NumerPresEm = sln(NumEgs+(1:NumEms));
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
ErrPresEm = sqrt(ErrPresEm);
L2ErrPres = norm(ErrPresEm);

%% Computing the L2-norm of the error in velocity 
CofPiolaBasEm = zeros(NumEms,4);
for j=1:4
  CofPiolaBasEm(:,j) = sln(QuadriMesh.elem2edge(:,j));
end
Phihat = zeros(NumEms,4,2);
WO = ones(NumEms,1);
MatP = zeros(NumEms,2,2);
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;
  vel = EqnBC.fxnu(qp);
  Phihat(:,4,1) = (1-xhat)*WO;  Phihat(:,2,1) = xhat*WO;
  Phihat(:,1,2) = (1-yhat)*WO;  Phihat(:,3,2) = yhat*WO;
  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;
  MatP(:,1,1) = J11./jac;  MatP(:,1,2) = J12./jac;
  MatP(:,2,1) = J21./jac;  MatP(:,2,2) = J22./jac;
  NumerVelEm = zeros(NumEms,2);
  for j=1:4
    NumerVelEm(:,1) = NumerVelEm(:,1) + CofPiolaBasEm(:,j).*...
      squeeze(MatP(:,1,1).*Phihat(:,j,1) + MatP(:,1,2).*Phihat(:,j,2));
    NumerVelEm(:,2) = NumerVelEm(:,2) + CofPiolaBasEm(:,j).*...
      squeeze(MatP(:,2,1).*Phihat(:,j,1) + MatP(:,2,2).*Phihat(:,j,2));
  end
  Verr = vel - NumerVelEm;
  ErrVelEm = ErrVelEm + GAUSSQUAD.RECT(k,3)*jac.*(Verr(:,1).^2+Verr(:,2).^2);
end
ErrVelEm = sqrt(ErrVelEm);
L2ErrVel = norm(ErrVelEm);

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

return;