function [L2ErrDspl_all_sqrt, L2ErrPres_all_sqrt,H1ErrDspl_all, ... 
  L2ErrVel_all_sqrt,L2ErrStrs_all_sqrt] ...
  = LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_Err( ...
    EqnBIC,TriMesh,tyme,NumerDspl,NumerPres,NumerVelCof,GAUSSQUAD)
% function [L2ErrDspl_all_sqrt, L2ErrPres_all_sqrt,H1ErrDspl_all] = ...
%   LinPoroElas_MS_TriIE_CGBR1_WGP0P0RT0_Err( ...
%   EqnBIC,TriMesh,tyme,NumerDspl,NumerPres,GAUSSQUAD)
%% Computing errors 
% James Liu, ColoState; 2012/07--2021/04 

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

%% Mesh info 
NumNds = TriMesh.NumNds;
NumEms = TriMesh.NumEms;
% NumEgs = TriMesh.NumEgs;
area = TriMesh.area;

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

%% Auxiliary: Computing basis function div average (BasFxnDivAvg) 
BasFxnDivAvg = zeros(NumEms,9);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.TRIG(k,1)*TriMesh.node(TriMesh.elem(:,1),:)...
     + GAUSSQUAD.TRIG(k,2)*TriMesh.node(TriMesh.elem(:,2),:)...
     + GAUSSQUAD.TRIG(k,3)*TriMesh.node(TriMesh.elem(:,3),:);   
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,...
        GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  for j=1:9
    BasFxnDivAvg(:,j) = BasFxnDivAvg(:,j) + BFDiv(:,j) .* GAUSSQUAD.TRIG(k,4).*area;
  end
end
BasFxnDivAvg = BasFxnDivAvg./area;

%% Computing erros at each time step  
L2ErrDspl = zeros(MT,1);
L2ErrPres = zeros(MT,1);
H1ErrDspl = zeros(MT,1); 
L2ErrVel = 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,9);
  for j=1:3  % Treating 4 vertices and 4 edges in one loop 
    k = TriMesh.elem(:,j);
    CofDsplEm(:,2*j-1) = NumerDspl(2*k-1,n+1);
    CofDsplEm(:,2*j  ) = NumerDspl(2*k,  n+1);
    k = TriMesh.elem2edge(:,j);
    CofDsplEm(:,6+j) = NumerDspl(2*NumNds+k,n+1);    
  end
  % Adjusting signs for the edge bubble functions 
  for j=1:3
    CofDsplEm(:,6+j) = CofDsplEm(:,6+j) .* TriMesh.EmEgSign(:,j);
  end

  %% 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);
  NumQuadPts = size(GAUSSQUAD.TRIG,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.TRIG(k,1)*TriMesh.node(TriMesh.elem(:,1),:)...
       + GAUSSQUAD.TRIG(k,2)*TriMesh.node(TriMesh.elem(:,2),:)...
       + GAUSSQUAD.TRIG(k,3)*TriMesh.node(TriMesh.elem(:,3),:);   
    BFVal = CG_TriBR1_BasFxnVal(TriMesh,...
         GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
    [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,...
        GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
    NumerDsplEm = zeros(NumEms,2);
    NumerDsplGradEm = zeros(NumEms,2,2);
    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);
    StrnPhi = zeros(NumEms,9,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,9,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:9
      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
    uval = EqnBIC.fxnu(qp,tnew);
    uGrad = EqnBIC.fxngradu(qp,tnew);
    pval = EqnBIC.fxnp(qp,tnew);
    vval = EqnBIC.fxnv(qp,tnew);
    sval = EqnBIC.fxnsigma(qp,tnew);
    uerr = uval - NumerDsplEm;
    uGraderr = uGrad - NumerDsplGradEm;
    perr = pval - NumerPresEm;
    X = qp(:,1) - TriMesh.EmCntr(:,1);
    Y = qp(:,2) - TriMesh.EmCntr(:,2);
    NV = [NumerVelCof(:,1,n) + NumerVelCof(:,3,n).*X,...
          NumerVelCof(:,2,n) + NumerVelCof(:,3,n).*Y];
    verr = vval - NV;  
    serr = sval - NumerStrs;
    L2ErrDsplEm = L2ErrDsplEm ...
      + GAUSSQUAD.TRIG(k,4)* (uerr(:,1).^2+uerr(:,2).^2).*area;
    H1SemiErrDsplEm = H1SemiErrDsplEm ...
      + (uGraderr(:,1,1).^2+uGraderr(:,1,2).^2+uGraderr(:,2,1).^2+uGraderr(:,2,2).^2) ...
      .*GAUSSQUAD.TRIG(k,4).*area;
    L2ErrPresEm = L2ErrPresEm + (perr.^2).*GAUSSQUAD.TRIG(k,4).*area;
    L2ErrVelEm = L2ErrVelEm + (verr(:,1).^2+verr(:,2).^2).*GAUSSQUAD.TRIG(k,4).*area;
    L2ErrStrsEm = L2ErrStrsEm + (serr(:,1,1).^2 + serr(:,1,2).^2 ...
                +  serr(:,2,1).^2 + serr(:,2,2).^2)...
                .* GAUSSQUAD.TRIG(k,4).*area;
  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);
  end
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;