%% JL20210116: BEING REVISED By James 
function [L2ErrDisp_GS, L2ErrPres_GS] ... 
  = LinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0_ErrGS( ...
    EqnBIC,QuadriMesh,tyme, NumerDisp_GS,NumerPres_GS, ... 
    old_NumerDisp_GS, old_NumerPres_GS, GAUSSQUAD,tnew) 
%% Computing errors in Gauss-Seidel iterations 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2017/07--2021/01 

% %% Equation info 
% lambda = EqnBIC.lambda;
%     mu = EqnBIC.mu;
%  alpha = EqnBIC.alpha;
%     c0 = EqnBIC.c0;

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

%% Tyme info 
NT = tyme.NT;  MT = tyme.MT;  Deltat = tyme.Deltat;

%% Computing erros at the nth. GS step 
L2ErrDisp_GS = 0;
L2ErrPres_GS = 0;
  
%% Elementwise cofficients for basis functions for displacement 
CofDispEm = zeros(NumEms,12);
for j=1:4  % Treating 4 vertices and 4 edges in one loop 
  k = QuadriMesh.elem(:,j);
  CofDispEm(:,2*j-1) = NumerDisp_GS(2*k-1);
  CofDispEm(:,2*j  ) = NumerDisp_GS(2*k);
  k = QuadriMesh.elem2edge(:,j);
  CofDispEm(:,8+j) = NumerDisp_GS(2*NumNds+k);    
end
% Adjusting signs for the edge bubble functions 
for j=1:4
  CofDispEm(:,8+j) = CofDispEm(:,8+j) .* QuadriMesh.EmEgSign(:,j);
end

old_CofDispEm = zeros(NumEms,12);
for j=1:4  % Treating 4 vertices and 4 edges in one loop 
  k = QuadriMesh.elem(:,j);
  old_CofDispEm(:,2*j-1) = old_NumerDisp_GS(2*k-1);
  old_CofDispEm(:,2*j  ) = old_NumerDisp_GS(2*k);
  k = QuadriMesh.elem2edge(:,j);
  old_CofDispEm(:,8+j) = old_NumerDisp_GS(2*NumNds+k);    
end
% Adjusting signs for the edge bubble functions 
for j=1:4
  old_CofDispEm(:,8+j) = old_CofDispEm(:,8+j) .* QuadriMesh.EmEgSign(:,j);
end

%% Elementwise constant numerical pressure for interior 
NumerPresEm = NumerPres_GS(1:NumEms);
old_NumerPresEm = old_NumerPres_GS(1:NumEms);

%% Computing errors on the entire mesh at this GS step
L2ErrDispEm = zeros(NumEms,1);
L2ErrPresEm = zeros(NumEms,1);  
MatJ = zeros(NumEms,2,2);
MatP = zeros(NumEms,2,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);
  J11 = CofA(:,2) + CofA(:,4)*yhat;  J12 = CofA(:,3) + CofA(:,4)*xhat;
  J21 = CofB(:,2) + CofB(:,4)*yhat;  J22 = CofB(:,3) + CofB(:,4)*xhat;
  MatJ(:,1,1) = J11;  MatJ(:,1,2) = J12;
  MatJ(:,2,1) = J21;  MatJ(:,2,2) = J22;
  jac = (CofB(:,3)+CofB(:,4)*xhat) .* (CofA(:,2)+CofA(:,4)*yhat)...
      - (CofA(:,3)+CofA(:,4)*xhat) .* (CofB(:,2)+CofB(:,4)*yhat);
  for i=1:2
    for j=1:2
      MatP(:,i,j) = MatJ(:,i,j)./jac;
    end
  end
  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];    
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  NumerDispEm = zeros(NumEms,2);
  for j=1:12
      NumerDispEm = NumerDispEm + CofDispEm(:,j).*BFVal(:,j);
  end
  old_NumerDispEm = zeros(NumEms,2);
  for j=1:12
    old_NumerDispEm = old_NumerDispEm + old_CofDispEm(:,j).*BFVal(:,j);
  end
  uerr = old_NumerDispEm - NumerDispEm;
  perr = old_NumerPresEm - NumerPresEm; 
  L2ErrDispEm = L2ErrDispEm ...
      + (uerr(:,1).^2+uerr(:,2).^2) .* jac * GAUSSQUAD.RECT(k,3);
  L2ErrPresEm = L2ErrPresEm + (perr.^2) .* jac * GAUSSQUAD.RECT(k,3);
end

%%
L2ErrDisp_GS = sqrt(sum(L2ErrDispEm));
L2ErrPres_GS = sqrt(sum(L2ErrPresEm));

return;