function [L2ErrPres,L2ErrVel] = Darcy_WG_TriP0P0RT0_Err(...
  EqnBC,TriMesh,NumerPresEm,NumerVelCofRT0,GAUSSQUAD)
%% Darcy: WG(P0,P0;RT0) on a triangular mesh: Computing errors 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;

%% Computing L2-norm of error in pressure 
ErrPresEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(:,1),:) ...
     + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(:,2),:) ...
     + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(:,3),:) ;
  perr = EqnBC.fxnp(qp) - NumerPresEm;
  ErrPresEm = ErrPresEm + GAUSSQUAD.TRIG(k,4)*(perr.^2);
end
L2ErrPres = sqrt(sum(ErrPresEm.*TriMesh.area));

%% Computing L2-norm of error in velocity 
ErrVelEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1) * TriMesh.node(TriMesh.elem(:,1),:) ...
     + GAUSSQUAD.TRIG(k,2) * TriMesh.node(TriMesh.elem(:,2),:) ...
     + GAUSSQUAD.TRIG(k,3) * TriMesh.node(TriMesh.elem(:,3),:) ;
  X = qp(:,1) - TriMesh.EmCntr(:,1);
  Y = qp(:,2) - TriMesh.EmCntr(:,2);
  NV = [ NumerVelCofRT0(:,1) + NumerVelCofRT0(:,3).*X ,...
         NumerVelCofRT0(:,2) + NumerVelCofRT0(:,3).*Y ];  
  uval = EqnBC.fxnu(qp);
  verr = uval - NV;
  ErrVelEm = ErrVelEm + GAUSSQUAD.TRIG(k,4)*(verr(:,1).^2 + verr(:,2).^2);
end
L2ErrVel = sqrt(sum(ErrVelEm.*TriMesh.area));

return;