function [AhErr] = LinElas_CG_TriEP1_ErrAh(EqnBC,TriMesh,sln,GAUSSQUAD)
%% Computing (Ah) energy-norm of the error 
% Solving lin.elasticity by 1st ord. enriched Lagrange elements on tri.mesh 
% Need displacement gradient 
% James Liu, ColoState; 2012/07--2020/11 

%% Equation info 
lambda = EqnBC.lambda;
    mu = EqnBC.mu;

%% Mesh info 
NumNds = TriMesh.NumNds;
NumEms = TriMesh.NumEms;
area = TriMesh.area;

%% Auxiliary quantities 
% Elementwise coeffs. in 9 BR1 basis functions for the numerical solution 
cof = zeros(NumEms,9);
for j=1:3  % Treating vertices and edges in one loop 
  k = TriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = TriMesh.elem2edge(:,j);
  cof(:,6+j) = sln(2*NumNds+k) .* TriMesh.EmEgSign(:,j);
end

%% Part I: Computing errors in strain squared 
ErrEmStrn2 = zeros(NumEms,1);
StrnPhi = zeros(NumEms,9,2,2);
ExactStrn = zeros(NumEms,2,2);
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),:);
  ug = EqnBC.fxnug(qp);
  ExactStrn(:,1,1) = ug(:,1,1);
  ExactStrn(:,2,2) = ug(:,2,2);
  ExactStrn(:,1,2) = 0.5 * (ug(:,1,2)+ug(:,2,1));
  ExactStrn(:,2,1) = ExactStrn(:,1,2);
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh, ...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  StrnPhi(:,:,1,1) = BFGrad(:,:,1,1);
  StrnPhi(:,:,2,2) = BFGrad(:,:,2,2);
  StrnPhi(:,:,1,2) = 0.5 * (BFGrad(:,:,1,2)+BFGrad(:,:,2,1));
  StrnPhi(:,:,2,1) = StrnPhi(:,:,1,2);
  NumerStrn = zeros(NumEms,2,2);
  for j=1:9
    NumerStrn(:,1,1) = NumerStrn(:,1,1) + cof(:,j).*StrnPhi(:,j,1,1);
    NumerStrn(:,1,2) = NumerStrn(:,1,2) + cof(:,j).*StrnPhi(:,j,1,2);
    NumerStrn(:,2,1) = NumerStrn(:,2,1) + cof(:,j).*StrnPhi(:,j,2,1);
    NumerStrn(:,2,2) = NumerStrn(:,2,2) + cof(:,j).*StrnPhi(:,j,2,2);
  end 
  serr = ExactStrn - NumerStrn;  
  ErrEmStrn2 = ErrEmStrn2 + GAUSSQUAD.TRIG(k,4) ...
    *( serr(:,1,1).^2 + serr(:,1,2).^2 + serr(:,2,1).^2 + serr(:,2,2).^2 );
end

%% Part II: Computing error average in dilation (div of displacement) 
ErrEmDivAvg = 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),:);
  divu = EqnBC.fxndivu(qp);
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh, ...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  NumerDiv = sum(cof.*BFDiv,2);
  derr = divu - NumerDiv;
  ErrEmDivAvg = ErrEmDivAvg + derr * GAUSSQUAD.TRIG(k,4);
end

%% Finally: Combined 
AhErr = sqrt(2*mu*sum(ErrEmStrn2.*area)+lambda*sum(ErrEmDivAvg.^2.*area));

return;