%% JL20180113: TO BE REVISED 
function [L2ErrDspl,L2ErrStrs,L2ErrDiv] = LinElas_CG_TriP12_Err( ...
  EqnBC,TriMesh,sln,GAUSSQUAD)
%% Computing L2-norm of errors in displacement, stress, dilation 
% James Liu, ColoState; 2012/07--2018/01 

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

%% Mesh info 
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

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

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

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

return;