% JL20210211: BEING REVISED By James 
function [L2L2ErrDisp, L2L2ErrPres] ...
  = LinPoroElas_MS_QuadriIE_CGEQ1_WGP0P0AC0_Err( ...
    EqnBIC, QuadriMesh, tyme, NumerDisp, NumerPres, GAUSSQUAD)
%% Computing errors when exact solutions are known 
% James Liu, ColoState; Zhuoran Wang, SYSU; 2017/07--2021/02 

% %% 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 spatial errors for time steps 
L2ErrDisp  = zeros(MT,1);
H1sErrDisp = zeros(MT,1);
L2ErrPres  = zeros(MT,1);
for n=1:MT
  tnew = tyme.tm(n+1);
  
  %% For numer.displacement: Elementwise cofficients of basis functions 
  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(2*k-1,n+1);
    CofDispEm(:,2*j  ) = NumerDisp(2*k,  n+1);
    k = QuadriMesh.elem2edge(:,j);
    CofDispEm(:,8+j) = NumerDisp(2*NumNds+k,n+1);    
  end
  % Adjusting signs for the edge bubble functions 
  for j=1:4
    CofDispEm(:,8+j) = CofDispEm(:,8+j) .* QuadriMesh.EmEgSign(:,j);
  end

  %% For numer.pressure: Elementwise constants 
  NumerPresEmP0 = NumerPres(1:NumEms,n+1);

  %% Computing spatial errors on the entire mesh 
  L2ErrDispEm = zeros(NumEms,1);
  L2ErrPresEm = zeros(NumEms,1);  
  H1sErrDispEm = 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];
    % 
    BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
    NumerDispQuadPt = zeros(NumEms,2);
    for j=1:12
      NumerDispQuadPt = NumerDispQuadPt + CofDispEm(:,j).*BFVal(:,j);
    end
    uval = EqnBIC.fxnu(qp,tnew);
    uerr = uval - NumerDispQuadPt;
    err2 = uerr(:,1).^2 + uerr(:,2).^2;
    L2ErrDispEm = L2ErrDispEm + GAUSSQUAD.RECT(k,3) * jac .* err2;
    % 
    [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
    NumerDispGradQuadPt = zeros(NumEms,2,2);
    for j=1:12
      for ii=1:2
        for jj=1:2 
          NumerDispGradQuadPt(:,ii,jj) = NumerDispGradQuadPt(:,ii,jj) ... 
            + CofDispEm(:,j).*BFGrad(:,j,ii,jj);
        end
      end
    end
    ug = EqnBIC.fxngradu(qp,tnew);
    ugerr = ug - NumerDispGradQuadPt;
    err2 = ugerr(:,1,1).^2 + ugerr(:,1,2).^2 + ugerr(:,2,1).^2 + ugerr(:,2,2).^2;
    H1sErrDispEm = H1sErrDispEm + GAUSSQUAD.RECT(k,3) * jac .* err2;
    % 
    pval = EqnBIC.fxnp(qp,tnew);
    perr = pval - NumerPresEmP0;
    err2 = perr.^2;
    L2ErrPresEm = L2ErrPresEm + GAUSSQUAD.RECT(k,3) * jac .* err2;
  end
  L2ErrDisp(n) = sqrt(sum(L2ErrDispEm));
  L2ErrPres(n) = sqrt(sum(L2ErrPresEm));
  H1sErrDisp(n) = sqrt(sum(H1sErrDispEm)); 
end

%% 
L2L2ErrDisp = sqrt(Deltat*sum(L2ErrDisp.^2));
% L0H1ErrDisp = max(sqrt(L2ErrDisp.^2+H1sErrDisp.^2));
L2L2ErrPres = sqrt(Deltat*sum(L2ErrPres.^2));

return;