% JL20210116: BEING REVISED By James 
function [L2ErrDspl_all_sqrt, L2ErrPres_all_sqrt, ... 
  H1ErrDspl_all, L2ErrVel_all_sqrt, L2ErrStrs_all_sqrt] ... 
  = LinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0_Err( ...
    EqnBIC,QuadriMesh,tyme,NumerDisp,NumerPres,NumerVelCof,GAUSSQUAD)
%% Computing errors 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2012/07--2021/01 

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

%% 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 each time step  
L2ErrDspl = zeros(MT,1);
L2ErrPres = zeros(MT,1);
H1ErrDspl = zeros(MT,1); 
L2ErrVel = zeros(MT,1);
L2ErrStrs = zeros(MT,1);

  
L2ErrDspl_all = 0;
H1ErrDspl_all = 0;
L2ErrPres_all = 0;
L2ErrVel_all = 0;
L2ErrStrs_all = 0;

for n=1:MT
  n; 
  tnew = tyme.tm(n+1);
  
  %% Elementwise cofficients for basis functions for displacement 
  CofDsplEm = zeros(NumEms,12);
  for j=1:4  % Treating 4 vertices and 4 edges in one loop 
    k = QuadriMesh.elem(:,j);
    CofDsplEm(:,2*j-1) = NumerDisp(2*k-1,n+1);
    CofDsplEm(:,2*j  ) = NumerDisp(2*k,  n+1);
    k = QuadriMesh.elem2edge(:,j);
    CofDsplEm(:,8+j) = NumerDisp(2*NumNds+k,n+1);    
  end
  % Adjusting signs for the edge bubble functions 
  for j=1:4
    CofDsplEm(:,8+j) = CofDsplEm(:,8+j) .* QuadriMesh.EmEgSign(:,j);
  end
  
%% 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);
  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);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  for j=1:12
    BasFxnDivAvg(:,j) = BasFxnDivAvg(:,j) + BFDiv(:,j) .* jac * GAUSSQUAD.RECT(k,3);
  end
end
BasFxnDivAvg = BasFxnDivAvg./QuadriMesh.area;

  %% Elementwise constant numerical pressure for interior 
  NumerPresEm = NumerPres(1:NumEms,n+1);

  %% Computing errors on the entire mesh 
  L2ErrDsplEm = zeros(NumEms,1);
  H1SemiErrDsplEm = zeros(NumEms,1);
  L2ErrPresEm = zeros(NumEms,1);  
  L2ErrVelEm = zeros(NumEms,1);
  L2ErrStrsEm = 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);
    NumerDsplEm = zeros(NumEms,2);
    NumerDsplGradEm = zeros(NumEms,2,2);
%     for j=1:12
%       NumerDsplEm = NumerDsplEm + CofDsplEm(:,j).*BFVal(:,j);
%       NumerDsplGradEm = NumerDsplGradEm + CofDsplEm(:,j).*BFGrad(:,j,:,:);
%     end
    NumerDsplEm(:,1) = sum(CofDsplEm.*squeeze(BFVal(:,:,1)),2);
    NumerDsplEm(:,2) = sum(CofDsplEm.*squeeze(BFVal(:,:,2)),2);
    NumerDsplGradEm(:,1,1) = sum(CofDsplEm.*squeeze(BFGrad(:,:,1,1)),2);
    NumerDsplGradEm(:,1,2) = sum(CofDsplEm.*squeeze(BFGrad(:,:,1,2)),2);
    NumerDsplGradEm(:,2,1) = sum(CofDsplEm.*squeeze(BFGrad(:,:,2,1)),2);
    NumerDsplGradEm(:,2,2) = sum(CofDsplEm.*squeeze(BFGrad(:,:,2,2)),2);
    uval = EqnBIC.fxnu(qp,tnew);
    uGrad = EqnBIC.fxnuGrad(qp,tnew);
    pval = EqnBIC.fxnp(qp,tnew);
    vval = EqnBIC.fxnv(qp,tnew);
    uerr = uval - NumerDsplEm;
    uGraderr = uGrad - NumerDsplGradEm;
    perr = pval - NumerPresEm;
%     
    X = qp(:,1) - QuadriMesh.EmCntr(:,1);
    Y = qp(:,2) - QuadriMesh.EmCntr(:,2);
    U = MatP(:,1,1) .* xhat - MatP(:,1,2) .* yhat;
    V = MatP(:,2,1) .* xhat - MatP(:,2,2) .* yhat;
    NV(:,1) = NumerVelCof(:,1,n) + NumerVelCof(:,3,n).*X + NumerVelCof(:,4,n).*U;
    NV(:,2) = NumerVelCof(:,2,n) + NumerVelCof(:,3,n).*Y + NumerVelCof(:,4,n).*V;
    verr = vval - NV;  
    
% The following is for stress tensor errors.
  % 
  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) + CofDsplEm(:,i).*StrsPhi(:,i,1,1);
    NumerStrs(:,2,2) = NumerStrs(:,2,2) + CofDsplEm(:,i).*StrsPhi(:,i,2,2);
    NumerStrs(:,1,2) = NumerStrs(:,1,2) + CofDsplEm(:,i).*StrsPhi(:,i,1,2);  
    NumerStrs(:,2,1) = NumerStrs(:,2,1) + CofDsplEm(:,i).*StrsPhi(:,i,2,1);  
  end
  serr = EqnBIC.fxnsigma(qp,tnew) - NumerStrs;
  
    L2ErrStrsEm = L2ErrStrsEm + (serr(:,1,1).^2 + serr(:,1,2).^2 ...
                 + serr(:,2,1).^2 + serr(:,2,2).^2) .*jac * GAUSSQUAD.RECT(k,3);
    
    L2ErrDsplEm = L2ErrDsplEm ...
      + (uerr(:,1).^2+uerr(:,2).^2).*jac * GAUSSQUAD.RECT(k,3);
    H1SemiErrDsplEm = H1SemiErrDsplEm ...
      + (uGraderr(:,1,1).^2+uGraderr(:,1,2).^2+uGraderr(:,2,1).^2+uGraderr(:,2,2).^2) ...
      .*jac * GAUSSQUAD.RECT(k,3);
    L2ErrPresEm = L2ErrPresEm + (perr.^2).*jac * GAUSSQUAD.RECT(k,3);
    L2ErrVelEm = L2ErrVelEm + (verr(:,1).^2+verr(:,2).^2).*jac * GAUSSQUAD.RECT(k,3);
  end
  H1ErrDsplEm = L2ErrDsplEm + H1SemiErrDsplEm;

  L2ErrDspl(n) = sqrt(sum(L2ErrDsplEm));
  L2ErrPres(n) = sqrt(sum(L2ErrPresEm));
  H1ErrDspl(n) = sqrt(sum(H1ErrDsplEm));
  L2ErrVel(n) = sqrt(sum(L2ErrVelEm));
  L2ErrStrs(n) = sqrt(sum(L2ErrStrsEm));

  L2ErrDspl_all = L2ErrDspl_all + Deltat*L2ErrDspl(n)^2;
  L2ErrPres_all = L2ErrPres_all + Deltat*L2ErrPres(n)^2; 
  L2ErrVel_all = L2ErrVel_all + Deltat*L2ErrVel(n)^2;
  L2ErrStrs_all = L2ErrStrs_all + Deltat*L2ErrStrs(n)^2;
  
  if H1ErrDspl(n) > H1ErrDspl_all
      H1ErrDspl_all = H1ErrDspl(n);
%   else
%       H1ErrDspl_all = H1ErrDspl_all;
  end
  
end

% L2ErrDspl_all = 0;
% H1ErrDspl_all = H1ErrDspl(1);
% L2ErrPres_all = 0;
% L2ErrVel_all = 0;
% for n=1:MT
%     L2ErrDspl_all = L2ErrDspl_all + Deltat*L2ErrDspl(n)*L2ErrDspl(n);
%     L2ErrPres_all = L2ErrPres_all + Deltat*L2ErrPres(n)*L2ErrPres(n); 
%     if H1ErrDspl_all < H1ErrDspl(n)
%         H1ErrDspl_all = H1ErrDspl(n);
%     end
%     L2ErrVel_all = L2ErrVel_all + Deltat*L2ErrVel(n)*L2ErrVel(n);
% end

L2ErrDspl_all_sqrt = sqrt(L2ErrDspl_all);
L2ErrPres_all_sqrt = sqrt(L2ErrPres_all);
L2ErrVel_all_sqrt = sqrt(L2ErrVel_all);
L2ErrStrs_all_sqrt = sqrt(L2ErrStrs_all);


return;