function [L2ErrPres] = Darcy_MFEM_TriRT0P0_Err(EqnBC,TriMesh,sln,GAUSSQUAD)
%% Darcy: MFEM(RT0,P0) on a triangular mesh: 
% Computing errors in pressure, etc.  
% JL20170211: TO BE FINISHED/REVISED FOR L2ErrVel 
% James Liu, ColoState; 2012/07--2018/12 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag=3
area = TriMesh.area;
LenEg = TriMesh.LenEg;
e2g = TriMesh.elem2edge;
sn = TriMesh.SignEmEg;

%% "Computing" numerical pressure 
NumerPresEm = sln(NumEgs+(1:NumEms));

%% Computing numerical velocity: coeff. in RT0 edge-bas.bas.fxns.(local) 
CofRT0 = zeros(NumEms,3);
for j=1:3
  CofRT0(:,j) = sln(e2g(:,j)).*sn(:,j);
end

%% All kinds of error 
err = zeros(1,3);
% err(1): L^2-norm of numer.vel. vs exact vel. over all elements 
% err(2): 
% err(3): L^2-norm of numer.pres. vs exact pres. over all elements  

%% Computing L2-norm of error in pressure 
tmp = 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),:);
  DiffP = abs(EqnBC.fxnp(qp)-NumerPresEm);
  tmp = tmp + GAUSSQUAD.TRIG(k,4)*(DiffP.^2);
end
err(3) = sqrt(sum(tmp.*area));
L2ErrPres = err(3);

%% Computing L2-norm of error in velocity 
err(1) = 0; 
BasFxn = zeros(NumEms,3,2);  % All elts., 3 edges, 2-vec. 
tmp = 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),:);
  vel = EqnBC.fxnu(qp);
  for j=1:3
    coeff = LenEg(e2g(:,j))./(2*area);
    BasFxn(:,j,:) = [coeff,coeff].*(qp-TriMesh.node(TriMesh.elem(:,j),:));
  end
  NumerVel = [CofRT0(:,1),CofRT0(:,1)] .* squeeze(BasFxn(:,1,:))...
           + [CofRT0(:,2),CofRT0(:,2)] .* squeeze(BasFxn(:,2,:))...
           + [CofRT0(:,3),CofRT0(:,3)] .* squeeze(BasFxn(:,3,:));
  DiffV = vel - NumerVel;
  DiffV2 = DiffV(:,1).^2 + DiffV(:,2).^2;
  tmp = tmp + GAUSSQUAD.TRIG(k,4)*DiffV2;
end
err(1) = sqrt(sum(tmp.*area));

return;