function [L2ErrDisp,L2ErrStrs,L2ErrDiv,L2RHS] = ...
  LinElas_CG_QuadriEQ1_Err(EqnBC,QuadriMesh,sln,GAUSSQUAD) 
%% LinElas.: CG: Enriched Lagrangian Q1 on quadri.mesh (w/ red.integration)
% Computing errors 
% James Liu, ColoState; 2012/07--2020/08 

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

%% Mesh info 
NumNds = QuadriMesh.NumNds;
NumEms = QuadriMesh.NumEms;
% NumEgs = QuadriMesh.NumEgs;
area = QuadriMesh.area;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Auxiliary quantities 
% Elementwise cofficients in 12 BR1 basis functions for numerical solution 
cof = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  cof(:,2*j-1) = sln(2*k-1);
  cof(:,2*j  ) = sln(2*k  );
  k = QuadriMesh.elem2edge(:,j);
  cof(:,8+j) = sln(2*NumNds+k) .* QuadriMesh.EmEgSign(:,j);
end

%% Computing L2-norm of error in displacement (element interiors) 
ErrEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  NumerDspl = zeros(NumEms,2);
  NumerDspl(:,1) = sum(cof.*squeeze(BFVal(:,:,1)),2);
  NumerDspl(:,2) = sum(cof.*squeeze(BFVal(:,:,2)),2);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  uval = EqnBC.fxnu(qp);
  uerr = uval - NumerDspl;
  ErrEm = ErrEm + (uerr(:,1).^2+uerr(:,2).^2) .*jac * GAUSSQUAD.RECT(k,3);
end
ErrEm = sqrt(ErrEm);
L2ErrDisp = norm(ErrEm,2);

%% Auxiliary: Computing basis function div average (BasFxnDivAvg) 
BasFxnDivAvg = zeros(NumEms,12);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  for j=1:12
    BasFxnDivAvg(:,j) = BasFxnDivAvg(:,j) ... 
      + BFDiv(:,j) .*jac./area * GAUSSQUAD.RECT(k,3);
  end
end

%% Computing L2-norm of error in stress
ErrEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  StrnPhi = zeros(NumEms,12,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,12,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:12
    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
  serr = EqnBC.fxnsigma(qp) - NumerStrs;
  ErrEm = ErrEm + (serr(:,1,1).^2 + serr(:,1,2).^2 ...
                +  serr(:,2,1).^2 + serr(:,2,2).^2)...
                .* jac*GAUSSQUAD.RECT(k,3);
end
L2ErrStrs = sqrt(sum(ErrEm));

%% Computing L2-norm of error in dilation (divergence of displacement) 
ErrEm = zeros(NumEms,1);
NumerSlnDivAvg = sum(cof.*BasFxnDivAvg,2);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  divu = EqnBC.fxndivu(qp);
  derr = divu - NumerSlnDivAvg;
  ErrEm = ErrEm + (derr.^2) .* jac * GAUSSQUAD.RECT(k,3);
end
ErrEm = sqrt(ErrEm);
L2ErrDiv = norm(ErrEm,2);

%% Computing L2-norm of elsticity RHS 
ErrEm = zeros(NumEms,1);
NumQuadPts = size(GAUSSQUAD.RECT,1);
GAUSSQUADRECT1 = ones(NumQuadPts,2) - GAUSSQUAD.RECT(:,1:2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  qp = [CofA(:,1) + CofA(:,2)*xhat + CofA(:,3)*yhat + CofA(:,4)*xhat*yhat,...
        CofB(:,1) + CofB(:,2)*xhat + CofB(:,3)*yhat + CofB(:,4)*xhat*yhat];
  ferr = EqnBC.fxnf(qp);
  ErrEm = ErrEm + (ferr(:,1).^2+ferr(:,2).^2) .*jac * GAUSSQUAD.RECT(k,3);
end
ErrEm = sqrt(ErrEm);
L2RHS = norm(ErrEm,2);

return;