function ErrPP = Darcy_DG_TriP2_ErrPP(EqnBC,TriMesh,PermK,...
  CofBDM1,GAUSSQUAD)
%% Darcy: DG: TriP2, Postprocessed velocity: Computing all kinds of error
% Assuming PermK is an elementwise scalar 
% James Liu, ColoState; 2012/07--2016/05 

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

%% All kinds of error 
ErrPP = zeros(1,2);
% ErrPP(1): L^2-norm of pp.numer.vel. vs exact vel. over all elements 
% ErrPP(2): ... on mesh skeleton 

%% 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];
% 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);

%% Computing errors in numerical velocity inside elements 
% Exact velocity computed at quadrature points 
ErrPP(1) = 0;
BasFxn = zeros(NumEms,6,2);  % all elts., 2-vec., 6 edge fxns. 
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),:);
  pg = EqnBC.fxnpg(qp);  K = EqnBC.fxnK(qp);  vel = -[K,K].*pg;
  % vel = -[PermK,PermK].*pg;
  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 + [CofBDM1(:,j),CofBDM1(:,j)] .* squeeze(BasFxn(:,j,:));
  end
  DiffV = vel - NumerVel;
  DiffV2 = DiffV(:,1).^2 + DiffV(:,2).^2;
  tmp = tmp + GAUSSQUAD.TRIG(k,4)*DiffV2;
end
ErrPP(1) = sqrt(sum(tmp.*area));

return;