% JL20151231: TO BE FINISHED/REVISED 
function err = Darcy_DG_TriP2_Err(EqnBC,TriMesh,PermK,CofNS,GAUSSQUAD)
%% Darcy: DG: TriP2: Computing all kinds of errors 
% Assuming PermK is an elementwise scalar 
% James Liu, ColoState; 2012/07--2016/05 

%% Mesh info 
NumEms = TriMesh.NumEms;
% Assuming TriMesh.flag>=2
area = TriMesh.area;

%% All kinds of error 
err = zeros(1,6);
% err(2): L^2-norm of numer.pres. vs exact pres. over all elements 
% err(3): L^\infty-norm of numer.pres. vs exact pres. over all elements 
% Need a high order GAUSSQUAD.TRIG for err(3) 

%% Computing error against the exact/true pressure 
err(2) = 0;  
err(3) = 0;
tmp = zeros(NumEms,1);
WO = ones(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);  Y = qp(:,2);
  BasFxn = [WO, X, Y, X.^2, X.*Y, Y.^2];
  NumerPres = dot(CofNS(:,1:6),BasFxn,2);
  DiffP = abs(EqnBC.fxnp(qp)-NumerPres);
  tmp = tmp + GAUSSQUAD.TRIG(k,4)*(DiffP.^2);
  err(3) = max(err(3),norm(DiffP,inf));
end
err(2) = sqrt(sum(tmp.*area));

return;