function [L2L2ErrDisp,L2L2ErrPres,L2L2ErrVel] ... 
  = LinPoroElas_MS_QuadriIE_WGP02P02AC02P0_WGP0P0AC0_Err( ...
    EqnBIC,QuadriMesh,tyme,NumerDisp,NumerPres,GAUSSQUAD) 
% Computing errors 
% James Liu, ColoState; Zhuoran Wang, SYSU; 2020/09--2021/04 

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

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

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

%% Temporary 
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBIC.fxnK, QuadriMesh,GAUSSQUAD);

%% As going through time-marching 
L2L2ErrDisp = 0;  
L2L2ErrPres = 0;
L2L2ErrVel  = 0;
for n=1:MT
  tnew = tyme.tm(n+1);

  % For step-n: 
  % Extracting elementwise constant displacement (interiors) 
  NumerDispEmP02 = zeros(NumEms,2);
  NumerDispEmP02(:,1) = NumerDisp(2*(1:NumEms)-1, n+1);
  NumerDispEmP02(:,2) = NumerDisp(2*(1:NumEms),   n+1);
  % Extracting elementwise constant pressure (interiors) 
  NumerPresEmP0 = NumerPres(1:NumEms, n+1);
  % Computing numerical Darcy velocity 
  [NumerPresCofP0P0, NumerVelCofAC0] ... 
    = Darcy_WG_QuadriP0P0AC0_PresVelFlux(... 
      QuadriMesh, PermK, NumerPres(:,n+1), GAUSSQUAD);

  % Now computing spatial errors in displacement & pressure simultaneously 
  ErrEmDisp2 = zeros(NumEms,1);
  ErrEmPres2 = zeros(NumEms,1);
  ErrEmVel2 = 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];
    % For displacement 
    uval = EqnBIC.fxnu(qp,tnew);
    uerr = uval - NumerDispEmP02;
    err2 = uerr(:,1).^2 + uerr(:,2).^2;
    ErrEmDisp2 = ErrEmDisp2  + GAUSSQUAD.RECT(k,3) * jac.*err2;
    % For pressure 
    pval = EqnBIC.fxnp(qp, tnew);
    perr = pval - NumerPresEmP0;
    err2 = perr.^2;
    ErrEmPres2 = ErrEmPres2  + GAUSSQUAD.RECT(k,3) * jac.*err2;
    % For velocity 
    FV = Hdiv_QuadriAC0_NmlzPiolaBas_FxnVal(QuadriMesh,xhat,yhat);
    NumerVel = zeros(NumEms,2);
    NumerVel(:,1) = sum(NumerVelCofAC0.*squeeze(FV(:,:,1)), 2);
    NumerVel(:,2) = sum(NumerVelCofAC0.*squeeze(FV(:,:,2)), 2);
    vval = EqnBIC.fxnv(qp, tnew);
    verr = vval - NumerVel;
    err2 = verr(:,1).^2 + verr(:,2).^2;
    ErrEmVel2 = ErrEmVel2  + GAUSSQUAD.RECT(k,3) * jac.*err2;    
  end
  L2L2ErrDisp = L2L2ErrDisp + Deltat * sum(ErrEmDisp2);
  L2L2ErrPres = L2L2ErrPres + Deltat * sum(ErrEmPres2);
  L2L2ErrVel  = L2L2ErrVel  + Deltat * sum(ErrEmVel2 );
end
L2L2ErrDisp = sqrt(L2L2ErrDisp);
L2L2ErrPres = sqrt(L2L2ErrPres);
L2L2ErrVel  = sqrt(L2L2ErrVel );

return;