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

%% Mesh info 
% Basic mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;
% Assuming TriMesh.flag=3
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 

%% Computing errors in numerical velocity inside elements 
% Exact velocity computed at quadrature points 
ErrPP(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 = -repmat(EqnBC.fxnK(qp),1,2).*EqnBC.fxnpg(qp);
  % pg = EqnBC.fxnpg(qp);  
  % K = EqnBC.fxnK(qp);  vel = -[K,K].*pg;
  % vel = -[PermK,PermK].*pg;
  for j=1:3
    coeff = LenEg(e2g(:,j))./(2*area);  % coefficients in local bas. fxns.
    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
ErrPP(1) = sqrt(sum(tmp.*area));

% %% Computing error in numerical velocity on the skeleton 
% % Exact velocity computed at quadrature points 
% ErrPP(2) = 0;
% x1 = TriMesh.node(TriMesh.edge(:,1),1);
% y1 = TriMesh.node(TriMesh.edge(:,1),2);
% x2 = TriMesh.node(TriMesh.edge(:,2),1);
% y2 = TriMesh.node(TriMesh.edge(:,2),2);
% tmp = zeros(NumEgs,1);
% NumQuadPts = size(GAUSSQUAD.LINE,1);
% for k=1:NumQuadPts
%   qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
%   vel = -repmat(EqnBC.fxnK(qp),1,2).*EqnBC.fxnpg(qp);
%   pg = EqnBC.fxnpg(qp);
%   K = EqnBC.fxnK(qp);  % Using averages across edges?
%   vel = -[K,K].*pg;
%   vn = dot(vel,TriMesh.NmlEg,2);
%   tmp = tmp + GAUSSQUAD.LINE(k,3)*((vn-NmlFlxEg).^2);
% end
% ErrPP(2) = sqrt(sum(tmp.*LenEg));

return;