function [L2ErrPres,L2ErrVel] = Darcy_MFEM_TriBDM1P0_Err(...
  EqnBC,TriMesh,sln,GAUSSQUAD)
%% Darcy: MFEM(BDM1,P0): Computing errors 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag>2
area = TriMesh.area;
% LenEg = TriMesh.LenEg;
% e2g = TriMesh.elem2edge;

%% Auxiliary quantities 
k1 = TriMesh.elem(:,1);  k2 = TriMesh.elem(:,2);  k3 = TriMesh.elem(:,3);
x1 = TriMesh.node(k1,1);  y1 = TriMesh.node(k1,2); 
x2 = TriMesh.node(k2,1);  y2 = TriMesh.node(k2,2);
x3 = TriMesh.node(k3,1);  y3 = TriMesh.node(k3,2);
eg1 = [x3,y3] - [x2,y2];
eg2 = [x1,y1] - [x3,y3];
eg3 = [x2,y2] - [x1,y1];
lg1 = sqrt(eg1(:,1).^2+eg1(:,2).^2);
lg2 = sqrt(eg2(:,1).^2+eg2(:,2).^2);
lg3 = sqrt(eg3(:,1).^2+eg3(:,2).^2);
tan = zeros(NumEms,3,2);
tan(:,1,:) = eg1./[lg1,lg1];
tan(:,2,:) = eg2./[lg2,lg2];
tan(:,3,:) = eg3./[lg3,lg3];
% nml = zeros(NumEms,3,2);
% nml(:,1,:) = [tan(:,1,2),-tan(:,1,1)];
% nml(:,2,:) = [tan(:,2,2),-tan(:,2,1)];
% nml(:,3,:) = [tan(:,3,2),-tan(:,3,1)];
% cosine theta 
ct1 = (lg2.^2+lg3.^2-lg1.^2)./(2*lg2.*lg3);
ct2 = (lg3.^2+lg1.^2-lg2.^2)./(2*lg3.*lg1);
ct3 = (lg1.^2+lg2.^2-lg3.^2)./(2*lg1.*lg2);
% sine theta 
st1 = sqrt(1-ct1.^2);
st2 = sqrt(1-ct2.^2);
st3 = sqrt(1-ct3.^2);

%% Retrieving numerical pressure 
NumerPresEm = sln(2*NumEgs+(1:NumEms));

%% Retrieving info on numerical velocity 
CofBDM1 = zeros(NumEms,6);
for j=1:3
  CofBDM1(:,2*j-1) = sln(2*TriMesh.elem2edge(:,j)-1) .* TriMesh.SignEmEg(:,j);
  CofBDM1(:,2*j)   = sln(2*TriMesh.elem2edge(:,j))   .* TriMesh.SignEmEg(:,j);
end
% Adjusting: Swapping & changing sign 
for ie=1:NumEms
  for j=1:3
    if (TriMesh.SignEmEg(ie,j)==-1) 
      tmp = CofBDM1(ie,2*j-1);
      CofBDM1(ie,2*j-1) = CofBDM1(ie,2*j);
      CofBDM1(ie,2*j) = tmp;
    end
  end
end

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

%% Computing L^2-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,6,2);  % All elementts, 6 edge fxns., 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);
  BaryCrd = TriMesh_BaryCrd(TriMesh,qp);
  BasFxn(:,1,:) =  repmat(BaryCrd(:,2)./st2,1,2).* squeeze(tan(:,3,:));
  BasFxn(:,2,:) = -repmat(BaryCrd(:,3)./st3,1,2).* squeeze(tan(:,2,:));
  BasFxn(:,3,:) =  repmat(BaryCrd(:,3)./st3,1,2).* squeeze(tan(:,1,:));
  BasFxn(:,4,:) = -repmat(BaryCrd(:,1)./st1,1,2).* squeeze(tan(:,3,:));
  BasFxn(:,5,:) =  repmat(BaryCrd(:,1)./st1,1,2).* squeeze(tan(:,2,:));
  BasFxn(:,6,:) = -repmat(BaryCrd(:,2)./st2,1,2).* squeeze(tan(:,1,:));
  NumerVel = zeros(NumEms,2);
  for j=1:6
    NumerVel = NumerVel + repmat(CofBDM1(:,j),1,2) .* squeeze(BasFxn(:,j,:));
  end
  DiffV = vel - NumerVel;
  DiffV2 = DiffV(:,1).^2 + DiffV(:,2).^2;
  tmp = tmp + GAUSSQUAD.TRIG(k,4)*DiffV2;
end
err(1) = sqrt(sum(tmp.*area));
L2ErrVel = err(1);

% % Error in pressure in element interior in L^\infty-norm 
% err(4) = 0;
% NumQuadPts = size(GAUSSQUAD.TRIG,1);  % Need a high order GAUSSQUAD.TRIG 
% 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)-NumerPres.elem);
%   err(4) = max(err(4),norm(DiffP,inf));
% end

return;