function [L2ErrDisp, L2ErrDiv, L2ErrStrs, L2RHS] ... 
  = LinElas_CG_TriEP1_Err(EqnBC,TriMesh,sln,GAUSSQUAD)
%% Computing L2-norms of errors in displacement, dilation, stress 
% 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 cofficients in 9 BR1 basis functions for 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

%% Auxiliary: Computing basis function div average (BasFxnDivAvg) 
BasFxnDivAvg = zeros(NumEms,9);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh, ...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  for j=1:9
    BasFxnDivAvg(:,j) = BasFxnDivAvg(:,j) + BFDiv(:,j) * GAUSSQUAD.TRIG(k,4);
  end
end

%% Auxiliary: 
NumerSlnDivAvg = sum(cof.*BasFxnDivAvg,2);

%% Computing L2-norm of error in displacement 
ErrEm = 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),:);
  BFVal = CG_TriBR1_BasFxnVal(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  NumerDisp = zeros(NumEms,2);
  NumerDisp(:,1) = sum(cof.*squeeze(BFVal(:,:,1)),2);
  NumerDisp(:,2) = sum(cof.*squeeze(BFVal(:,:,2)),2);
  uval = EqnBC.fxnu(qp);
  uerr = uval - NumerDisp;
  ErrEm = ErrEm + (uerr(:,1).^2+uerr(:,2).^2) * GAUSSQUAD.TRIG(k,4);
end
ErrEm = sqrt(ErrEm.*area);
L2ErrDisp = norm(ErrEm,2);

%% Computing L2-norm of error in dilation (div of displacement) 
ErrEm = 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);
  derr = divu - NumerSlnDivAvg;
  ErrEm = ErrEm + (derr.^2) * GAUSSQUAD.TRIG(k,4);
end
ErrEm = sqrt(ErrEm.*area);
L2ErrDiv = norm(ErrEm,2);

%% Computing L2-norm of error in stress
ErrEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh, ...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  % 
  StrnPhi = zeros(NumEms,9,2,2);
  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);
  % 
  StrsPhi = zeros(NumEms,9,2,2);
  StrsPhi(:,:,1,1) = (2*mu)*StrnPhi(:,:,1,1) + lambda * BasFxnDivAvg;
  StrsPhi(:,:,2,2) = (2*mu)*StrnPhi(:,:,2,2) + lambda * BasFxnDivAvg;
  StrsPhi(:,:,1,2) = (2*mu)*StrnPhi(:,:,1,2);
  StrsPhi(:,:,2,1) = StrsPhi(:,:,1,2);
  % 
  NumerStrs = zeros(NumEms,2,2);
  for i=1:9
    NumerStrs(:,1,1) = NumerStrs(:,1,1) + cof(:,i).*StrsPhi(:,i,1,1);
    NumerStrs(:,2,2) = NumerStrs(:,2,2) + cof(:,i).*StrsPhi(:,i,2,2);
    NumerStrs(:,1,2) = NumerStrs(:,1,2) + cof(:,i).*StrsPhi(:,i,1,2);  
    NumerStrs(:,2,1) = NumerStrs(:,2,1) + cof(:,i).*StrsPhi(:,i,2,1);  
  end
  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),:);
  serr = EqnBC.fxnsigma(qp) - NumerStrs;
  ErrEm = ErrEm + (serr(:,1,1).^2 + serr(:,1,2).^2 ...
                 + serr(:,2,1).^2 + serr(:,2,2).^2) * GAUSSQUAD.TRIG(k,4);
end
ErrEm = sqrt(ErrEm.*area);
L2ErrStrs = norm(ErrEm,2);

%% Computing L2-norm of the elsticity RHS 
ErrEm = 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),:);
  ferr = EqnBC.fxnf(qp);
  ErrEm = ErrEm + (ferr(:,1).^2+ferr(:,2).^2) * GAUSSQUAD.TRIG(k,4);
end
ErrEm = sqrt(ErrEm.*area);
L2RHS = norm(ErrEm,2);

return;