%% JL20210207: BEING REVISED By James 
function [L2ErrDisp,L2ErrStrs,L2ErrDiv] ... 
  = LinElas_WG_QuadriP02P02AC02P0_Err(EqnBC,QuadriMesh,sln,GAUSSQUAD)
%% Lin.Elas.: WG(P0^2,P0^2;AC0^2,P0)Quadri: Computing errors 
% James Liu, ColoState; 2017/02--2021/02 

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

%% Mesh info 
NumEms = QuadriMesh.NumEms;
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% Numer.sln.: Converted to elementwise coefficients in basis functions 
cof = zeros(NumEms,10);
cof(:,1) = sln(2*(0:NumEms-1)+1);
cof(:,2) = sln(2*(0:NumEms-1)+2);
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  cof(:,2+2*j-1) = sln(2*NumEms+2*k-1);
  cof(:,2+2*j  ) = sln(2*NumEms+2*k  );
end

%% Extracting numer.disp. in element interior P0^2, numer,div. P0 
NumerDispP02 = cof(:,1:2);
CDWDB = WG_QuadriP02P02P0_CofNmlzBas_DiscWkDivBasFxn(QuadriMesh);
NumerDivP0 = sum(cof.*CDWDB,2); 

%% Errors:  
%  L2ErrDisp: L^2-norm of error in displacement (for element interior) 
%  L2ErrStrs: L^2-norm of error in stress 
%  L2ErrDiv:  L^2-norm of error in divergence (dilation)

%% Computing L2-norm of error in displacement (element interiors) 
ErrDispEm = 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];
  uval = EqnBC.fxnu(qp);
  uerr = uval - NumerDispP02;
  err2 = uerr(:,1).^2 + uerr(:,2).^2;
  ErrDispEm = ErrDispEm + GAUSSQUAD.RECT(k,3) * jac .* err2;
end
L2ErrDisp = sqrt(sum(ErrDispEm));

%% JL20210207: TO BE REVISED FOR "NumerStrs" 
%% Computing L2-norm of error in stress 
ErrStrsEm = 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];
  % JL20210207: TO BE FINISHED 
  NumerStrs = zeros(NumEms,2,2);
  sval = EqnBC.fxnsigma(qp);
  serr = sval - NumerStrs;
  err2 = serr(:,1,1).^2 + serr(:,1,2).^2 + serr(:,2,1).^2 + serr(:,2,2).^2;
  ErrStrsEm = ErrStrsEm + GAUSSQUAD.RECT(k,3) * jac .* err2;
end
L2ErrStrs = sqrt(sum(ErrStrsEm));

%% JL20210206: TO BE REVISED FOR "NumerDiv" 
%% Computing L2-norm of error in dilation (divergence of displacement)
ErrDivEm = 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];
  dval = EqnBC.fxndivu(qp);
  derr = dval - NumerDivP0;
  err2 = derr.^2;
  ErrDivEm = ErrDivEm + GAUSSQUAD.RECT(k,3) * jac.*err2;
end
L2ErrDiv = sqrt(sum(ErrDivEm));

return;