%% JL20210116: BEING REVISED By James 
function [NumerDisp,NumerPres,NumerVelCof] ...
  = NlinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0_AsmSlv( ...
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,QuadriMesh,tyme,GAUSSQUAD)
%% Solving nlin. PDE system by operator-splitting (OS) / Gauss-Seidel (GS) 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2012/07--2021/01 

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

%% Mesh info 
NumNds = QuadriMesh.NumNds;
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;

%% Gauss-Seidel 
GS_step = 60000;

%% Sorting out boundary info 
% For Elas: Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdgeElas = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==1);
NeumannEdgeElas   = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==2);
SymmetryEdgeX     = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==4);
SymmetryEdgeY     = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==5);
NumDirichletEgsElas = size(DirichletEdgeElas,1);
NumNeumannEgsElas = size(NeumannEdgeElas,1);
NumSymmetryEgsX   = size(SymmetryEdgeX,1);
NumSymmetryEgsY   = size(SymmetryEdgeY,1);
% JL20160808: TO BE REVISED FOR EFFICIENCY !
DirichletNodeElasFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgsElas 
  k1 = QuadriMesh.edge(DirichletEdgeElas(ig),1);
  k2 = QuadriMesh.edge(DirichletEdgeElas(ig),2);
  DirichletNodeElasFlag(k1) = 1;
  DirichletNodeElasFlag(k2) = 1;
end
DirichletNodeElas = find(DirichletNodeElasFlag);
NumDirichletNdsElas = sum(DirichletNodeElasFlag);
% For Darcy: Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdgeDarcy = find(BndryCondTypeDarcy(QuadriMesh.BndryEdge+1)==1);
NeumannEdgeDarcy   = find(BndryCondTypeDarcy(QuadriMesh.BndryEdge+1)==2);
NumDirichletEgsDarcy = size(DirichletEdgeDarcy,1);
NumNeumannEgsDarcy = size(NeumannEdgeDarcy,1);

%% Setup for the discrete linear system 
DOFsDisp = 2*NumNds + NumEgs;
DOFsPres = NumEms + NumEgs;
NumerDisp = zeros(DOFsDisp,NT+1);
NumerPres = zeros(DOFsPres,NT+1);
NumerVelCof = zeros(NumEms,4,NT);

% Displacement and pressure GS solutions at nth. time step
% They are used to construct system at nth time step, so I suppose no need
% to store them, just use them and delete.
% Might need to adjust later.
old_NumerDspl_GS = zeros(DOFsDisp,1);
old_NumerPres_GS = zeros(DOFsPres,1);
NumerDisp_GS = zeros(DOFsDisp,GS_step+1);
NumerPres_GS = zeros(DOFsPres,GS_step+1);


%% Assembling GlbMatElas 
% Element matrix: strain-strain 
EltMatStrnStrn = zeros(NumEms,12,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);
  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);
  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);
  for i=1:12
    for j=1:12
      EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j) ...
        + ( StrnPhi(:,i,1,1).*StrnPhi(:,j,1,1) ...
          + StrnPhi(:,i,1,2).*StrnPhi(:,j,1,2) ...
          + StrnPhi(:,i,2,1).*StrnPhi(:,j,2,1) ...
          + StrnPhi(:,i,2,2).*StrnPhi(:,j,2,2) ) .* jac * GAUSSQUAD.RECT(k,3);
    end
  end
end
% Element matrix: div(avg)-div(avg) 
EltMatDivDiv = zeros(NumEms,12,12);
EltDivAvg = 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);
  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
    EltDivAvg(:,j) = EltDivAvg(:,j) ...
      + GAUSSQUAD.RECT(k,3) * BFDiv(:,j) .* jac ./ area;
  end
end
for i=1:12
  for j=1:12
    EltMatDivDiv(:,i,j) = EltDivAvg(:,i) .* EltDivAvg(:,j) .* area;
  end
end
% Altogether 
% EltMat = zeros(NumEms,12,12);
EltMat = (2*mu)*EltMatStrnStrn + lambda*EltMatDivDiv;
% Modification with signs for edge normal vectors 
for i=1:4
  for j=1:12
    EltMat(:,8+i,j) = EltMat(:,8+i,j) .* QuadriMesh.EmEgSign(:,i);
  end
end
for i=1:12
  for j=1:4
    EltMat(:,i,8+j) = EltMat(:,i,8+j) .* QuadriMesh.EmEgSign(:,j);
  end
end
% Positions 
GlbMatElas = sparse(DOFsDisp,DOFsDisp);
pos = zeros(NumEms,12);
for j=1:4  % Treating vertices and edges in one loop 
  k = QuadriMesh.elem(:,j);
  pos(:,2*j-1) = 2*k - 1;
  pos(:,2*j  ) = 2*k;
  k = QuadriMesh.elem2edge(:,j);
  pos(:,8+j) = 2*NumNds + k;
end
% Now assembling 
for i=1:12 
  II = pos(:,i);
  for j=1:12 
    JJ = pos(:,j);
    GlbMatElas = GlbMatElas + sparse(II,JJ, EltMat(:,i,j), DOFsDisp,DOFsDisp);
  end
end

%% Assembling GlbMatDarcy 
% GlbMatDarcyMass = sparse(DOFsPres,DOFsPres);
II = (1:NumEms);  JJ = II;
GlbMatDarcyMass = sparse(II,JJ, area, DOFsPres,DOFsPres);
GlbMatDarcy = c0 * GlbMatDarcyMass;

%% JL20171231: TO BE REVISED 
%% Assembling GlbMatDarcyElas 
% Coupling: P0 (test fxns.) vs. averages of BR1 (trial fxns.) 
EltMatP0DivBR1 = zeros(NumEms,12);
for j=1:12
  EltMatP0DivBR1(:,j) = EltDivAvg(:,j) .* area;
end
% Adjusting for signs of edge normal vectors 
for j=1:4
  EltMatP0DivBR1(:,8+j) = EltMatP0DivBR1(:,8+j) .* QuadriMesh.EmEgSign(:,j);
end
GlbMatDarcyElas = sparse(DOFsPres,DOFsDisp);
pos = zeros(NumEms,12);
for j=1:4  % For nodes and edges (treated simultaneously) 
 k = QuadriMesh.elem(:,j);
 pos(:,2*j-1) = 2*k - 1;
 pos(:,2*j  ) = 2*k;
 k = QuadriMesh.elem2edge(:,j);
 pos(:,8+j) = 2*NumNds + k;
end
II = (1:NumEms);
for j=1:12
  JJ = pos(:,j);
  GlbMatDarcyElas = GlbMatDarcyElas ... 
    + sparse(II,JJ, alpha*EltMatP0DivBR1(:,j), DOFsPres,DOFsDisp);
end

%% Producing GlbMatElasDarcy: Coupling Elas (test fxns.) w/ Darcy (trial fxns.)
GlbMatElasDarcy = -GlbMatDarcyElas';

%% JL20210116: TEMPORARILY 0 
%% Computing projections/interpolations 
NumerDisp(:,1) = 0;
NumerPres(:,1) = 0;

%% Now time-marching 
fprintf('Time-marching step: ');
StepInTimeMarching = 0;
for n=1:MT
  fprintf(repmat('\b', 1, StepInTimeMarching));
  StepInTimeMarching = fprintf('%d', n);
  tnew = tyme.tm(n+1);
  
  old_NumerDspl_GS(:,1) = NumerDisp(:,n);
  old_NumerPres_GS(:,1) = NumerPres(:,n);
  
  % Sampling permeability 
  PermK = Darcy_SmplnNlinPerm_QuadriMesh( ...
    EqnBIC,QuadriMesh,GAUSSQUAD,NumerDisp(:,n));
  GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh,PermK,GAUSSQUAD);
  CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
  EltMat = zeros(NumEms,5,5);
  % JL20171231: TO BE REVISED FOR EFFECIENCY 
  for ie=1:NumEms
    EltGM = squeeze(GMK(ie,:,:));
    for i=1:5
      for j=1:5
        EltMat(ie,i,j) = squeeze(CDWGB(ie,i,:))' *EltGM* squeeze(CDWGB(ie,j,:));
      end
    end
  end
  pos = zeros(NumEms,5);
  pos(:,1) = (1:NumEms);
  for j=1:4
    pos(:,1+j) = NumEms + QuadriMesh.elem2edge(:,j);
  end
  GlbMatDarcyStf = sparse(DOFsPres,DOFsPres);
  for i=1:5
    II = pos(:,i);
    for j=1:5
      JJ = pos(:,j);
      GlbMatDarcyStf = GlbMatDarcyStf ...
        + sparse(II,JJ, EltMat(:,i,j), DOFsPres,DOFsPres);
    end
  end
  GlbMatDarcy = Deltat * GlbMatDarcyStf;
 
  for m = 1:GS_step
    % JL20171221: TO BE REVISED 
    % Assembling GlbRHS: For solid 
    % Computing elementwise contributions of 12 basis functions 
    EltRHSElas = 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);
      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];
      fval = EqnBIC.fxnf(qp,tnew);
      BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
      for j=1:12
        EltRHSElas(:,j) = EltRHSElas(:,j) + GAUSSQUAD.RECT(k,3)...
          * jac .*( BFVal(:,j,1).*fval(:,1) + BFVal(:,j,2).*fval(:,2));
      end  
    end 
    % 
    for j=1:12  %  alpha(p,div) 
      EltRHSElas(:,j) = EltRHSElas(:,j) + ...
        alpha * area .* old_NumerPres_GS(1:NumEms,1) .* EltDivAvg(:,j);
    end  
    % Modification for 4 edge-based basis functions 
    for j=1:4
      EltRHSElas(:,8+j) = EltRHSElas(:,8+j) .* QuadriMesh.EmEgSign(:,j);
    end
  
    % Assembling contributions of elementwise 12 basis functions 
    GlbRhsDisp = zeros(DOFsDisp,1);
    for ie=1:NumEms
      for j=1:4  % Treating vertices and edges in one loop 
        k = QuadriMesh.elem(ie,j);
        GlbRhsDisp(2*k-1) = GlbRhsDisp(2*k-1) + EltRHSElas(ie,2*j-1);
        GlbRhsDisp(2*k  ) = GlbRhsDisp(2*k  ) + EltRHSElas(ie,2*j  );
        k = QuadriMesh.elem2edge(ie,j);
        GlbRhsDisp(2*NumNds+k) = GlbRhsDisp(2*NumNds+k) + EltRHSElas(ie,8+j);
      end
    end

    if NumNeumannEgsElas>0
      % For end-nodes 
      NeumannBC = zeros(NumNeumannEgsElas,2,2);
      j1 = QuadriMesh.edge(NeumannEdgeElas,1);  
      j2 = QuadriMesh.edge(NeumannEdgeElas,2);
      x1 = QuadriMesh.node(j1,1);  y1 = QuadriMesh.node(j1,2);
      x2 = QuadriMesh.node(j2,1);  y2 = QuadriMesh.node(j2,2);
      LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
      NumQuadPts = size(GAUSSQUAD.LINE,1);
      for k=1:NumQuadPts
        qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
        tNval = EqnBIC.fxntN(qp,tnew);
        NeumannBC(:,1,1) = NeumannBC(:,1,1) + GAUSSQUAD.LINE(k,3)...
          * GAUSSQUAD.LINE(k,1) * tNval(:,1) .* LenEgNeumann;
        NeumannBC(:,1,2) = NeumannBC(:,1,2) + GAUSSQUAD.LINE(k,3)...
          * GAUSSQUAD.LINE(k,1) * tNval(:,2) .* LenEgNeumann;
        NeumannBC(:,2,1) = NeumannBC(:,2,1) + GAUSSQUAD.LINE(k,3)...
          * GAUSSQUAD.LINE(k,2) * tNval(:,1) .* LenEgNeumann;
        NeumannBC(:,2,2) = NeumannBC(:,2,2) + GAUSSQUAD.LINE(k,3)...
          * GAUSSQUAD.LINE(k,2) * tNval(:,2) .* LenEgNeumann;
      end
      % JL20171231: TO BE REVISED 
      for ig=1:NumNeumannEgsElas 
        j1 = QuadriMesh.edge(NeumannEdgeElas(ig),1);
        j2 = QuadriMesh.edge(NeumannEdgeElas(ig),2);
        GlbRhsDisp(2*j1-1) = GlbRhsDisp(2*j1-1) + NeumannBC(ig,1,1);
        GlbRhsDisp(2*j1  ) = GlbRhsDisp(2*j1  ) + NeumannBC(ig,1,2);
        GlbRhsDisp(2*j2-1) = GlbRhsDisp(2*j2-1) + NeumannBC(ig,2,1);
        GlbRhsDisp(2*j2  ) = GlbRhsDisp(2*j2  ) + NeumannBC(ig,2,2);
      end
      % For edges 
      NeumannBC = zeros(NumNeumannEgsElas,1);
      j1 = QuadriMesh.edge(NeumannEdgeElas,1);
      j2 = QuadriMesh.edge(NeumannEdgeElas,2);
      x1 = QuadriMesh.node(j1,1);  y1 = QuadriMesh.node(j1,2);
      x2 = QuadriMesh.node(j2,1);  y2 = QuadriMesh.node(j2,2);
      LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
      nml = QuadriMesh.EgNml(NeumannEdgeElas,:);
      NumQuadPts = size(GAUSSQUAD.LINE,1);
      for k=1:NumQuadPts
        qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
        tNval = EqnBIC.fxntN(qp,tnew);
        psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
        NeumannBC = NeumannBC ...
          + GAUSSQUAD.LINE(k,3) * sum(tNval.*nml,2) * psi .* LenEgNeumann;
      end
      GlbRhsDisp(2*NumNds+NeumannEdgeElas) = ...
      GlbRhsDisp(2*NumNds+NeumannEdgeElas) + NeumannBC;
    end

    % For solution 
    NumerDisp_GS(:,m+1) = zeros(DOFsDisp,1);
  
    % For Elas.: Enforcing Dirichlet boundary conditions as essential 
    if (NumDirichletEgsElas>0)  % Assuming NumDirichletEgsElas>0 
      % For nodes: 
      DirichletNdElas = QuadriMesh.node(DirichletNodeElas,:);
      uDval = EqnBIC.fxnuD(DirichletNdElas,tnew);
      NumerDisp_GS(2*DirichletNodeElas-1,m+1) = uDval(:,1);
      NumerDisp_GS(2*DirichletNodeElas,m+1  ) = uDval(:,2);
    end

    if (NumDirichletEgsElas>0)
      intgrl1 = zeros(NumDirichletEgsElas,1);
      intrpl = zeros(NumDirichletEgsElas,1);
      intgrl2 = zeros(NumDirichletEgsElas,1);
      x1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeElas,1),1);
      y1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeElas,1),2);
      x2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeElas,2),1);
      y2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeElas,2),2);
      LenDirichletEdgeElas = sqrt((x2-x1).^2+(y2-y1).^2);
      nml = QuadriMesh.EgNml(DirichletEdgeElas,:);
      NumQuadPts = size(GAUSSQUAD.LINE,1);
      for k=1:NumQuadPts
        qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
        uD = EqnBIC.fxnu(qp,tnew);
        psi = GAUSSQUAD.LINE(k,1)*GAUSSQUAD.LINE(k,2);
        intgrl1 = intgrl1 + GAUSSQUAD.LINE(k,3)*psi*LenDirichletEdgeElas;
        intrpl = GAUSSQUAD.LINE(k,1)*EqnBIC.fxnu(QuadriMesh.node(QuadriMesh.edge(DirichletEdgeElas,1),:),tnew)...
               + GAUSSQUAD.LINE(k,2)*EqnBIC.fxnu(QuadriMesh.node(QuadriMesh.edge(DirichletEdgeElas,2),:),tnew); 
        res = uD - intrpl;
        intgrl2 = intgrl2 + GAUSSQUAD.LINE(k,3)*sum(res.*nml,2).*LenDirichletEdgeElas;
      end
      NumerDisp_GS(2*NumNds+DirichletEdgeElas,m+1) = intgrl2./intgrl1;
    end
 
    % Reducing...
    GlbRhsDisp = GlbRhsDisp - GlbMatElas * NumerDisp_GS(:,m+1);

    % For flag 
    flagElas = ones(DOFsDisp,1);
    flagElas(2*DirichletNodeElas-1) = 0;
    flagElas(2*DirichletNodeElas  ) = 0;
    flagElas(2*NumNds+DirichletEdgeElas) = 0;
    FreeNdEmFreeEgElas = find(flagElas);
 
    % Solving the reduced linear system 
    NumerDisp_GS(FreeNdEmFreeEgElas,m+1) = ... 
      GlbMatElas(FreeNdEmFreeEgElas,FreeNdEmFreeEgElas) ...
        \ GlbRhsDisp(FreeNdEmFreeEgElas);

    % Now for the fluid part
    % Assembling GlbRHS: For fluid 
    GlbRhsPres = zeros(DOFsPres,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];
      sval = EqnBIC.fxns(qp,tnew);
      GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) ...
        + GAUSSQUAD.RECT(k,3) * sval .* jac * Deltat;
    end
    
    % Assembling GlbRHSPres: Pressure at old time
    GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) ...
      + c0 * area .* squeeze(NumerPres(1:NumEms,n));
  
    % Assembling GlbRHSPres: Solid displacement divergence at old time 
    tmp = zeros(NumEms,12);
    EltDispDivOld = zeros(NumEms,1);
    for j=1:4  % Treating vertices and edges in one loop 
      k = QuadriMesh.elem(:,j);
      tmp(:,2*j-1) = NumerDisp(2*k-1,n);
      tmp(:,2*j  ) = NumerDisp(2*k,  n);
      k = QuadriMesh.elem2edge(:,j);
      tmp(:,8+j) = NumerDisp(2*NumNds+k,n);
    end
    EltDispDivOld(:,1) = sum(EltMatP0DivBR1.*tmp,2);
    GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) + alpha * EltDispDivOld(:,1);

    % Assembling GlbRHSPres: Solid displacement divergence at old GS step 
    tmp = zeros(NumEms,12);
    EltDispDivGS = zeros(NumEms,1);
    for j=1:4  % Treating vertices and edges in one loop 
      k = QuadriMesh.elem(:,j);
      tmp(:,2*j-1) = NumerDisp_GS(2*k-1,m+1);
      tmp(:,2*j  ) = NumerDisp_GS(2*k,m+1);
      k = QuadriMesh.elem2edge(:,j);
      tmp(:,8+j) = NumerDisp_GS(2*NumNds+k,m+1);
    end
    EltDispDivGS(:,1) = sum(EltMatP0DivBR1.*tmp,2);
    GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) - alpha * EltDispDivGS(:,1);
  
    % For fluid/Darcy: Incorporating Neumann (flux) boundary conditions 
    if (NumNeumannEgsDarcy>0)
      NeumannBC = zeros(NumNeumannEgsDarcy,1);
      x1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdgeDarcy,1),1);
      y1 = QuadriMesh.node(QuadriMesh.edge(NeumannEdgeDarcy,1),2);
      x2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdgeDarcy,2),1);
      y2 = QuadriMesh.node(QuadriMesh.edge(NeumannEdgeDarcy,2),2);
      LenNeumannEg = sqrt((x2-x1).^2+(y2-y1).^2);
      NumQuadPts = size(GAUSSQUAD.LINE,1);
      for k=1:NumQuadPts
        qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
        uNval = EqnBIC.fxnuN(qp,tnew);
        NeumannBC = NeumannBC + GAUSSQUAD.LINE(k,3) * uNval .* LenNeumannEg; 
      end
      GlbRhsPres(NumEms+NeumannEdgeDarcy) = ...
      GlbRhsPres(NumEms+NeumannEdgeDarcy) - NeumannBC; 
    end
    
    % For solution
    NumerPres_GS(:,m+1) = zeros(DOFsPres,1);
  
    % For Darcy: Enforcing Dirichlet boundary conditions as essential
    if (NumDirichletEgsDarcy>0)  % Assuming NumDirichletEgsDarcy>0 
      pDavg = zeros(NumDirichletEgsDarcy,1);
      x1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeDarcy,1),1);
      y1 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeDarcy,1),2);
      x2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeDarcy,2),1);
      y2 = QuadriMesh.node(QuadriMesh.edge(DirichletEdgeDarcy,2),2);
      NumQuadPts = size(GAUSSQUAD.LINE,1);
      for k=1:NumQuadPts
        qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
        pDval = EqnBIC.fxnpD(qp,tnew);
        pDavg = pDavg + GAUSSQUAD.LINE(k,3) * pDval;
      end
      NumerPres_GS(NumEms+DirichletEdgeDarcy,m+1) = pDavg;
    end
  
    % Reducing...
    GlbRhsPres = GlbRhsPres - GlbMatDarcy*NumerPres_GS(:,m+1);
  
    % For flag 
    flagPres = ones(DOFsPres,1);
    % Darcy D.B.C edges
    flagPres(NumEms+DirichletEdgeDarcy) = 0;
    FreeNdEmFreeEgPres = find(flagPres);
  
    GlbMatDarcy(FreeNdEmFreeEgPres,FreeNdEmFreeEgPres);
  
    % Solving the reduced linear system 
    NumerPres_GS(FreeNdEmFreeEgPres,m+1) ... 
      = GlbMatDarcy(FreeNdEmFreeEgPres,FreeNdEmFreeEgPres) ...
        \ GlbRhsPres(FreeNdEmFreeEgPres);

    % Now calculating errors in GS iteration 
    [L2ErrDisp_GS, L2ErrPres_GS] ... 
      = LinPoroElas_OS_QuadriIE_CGEQ1_WGP0P0AC0_ErrGS( ...
        EqnBIC,QuadriMesh,tyme, NumerDisp_GS(:,m+1), NumerPres_GS(:,m+1),...
        NumerDisp_GS(:,m), NumerPres_GS(:,m), GAUSSQUAD,tnew);
    
    if (L2ErrDisp_GS>1e-12 || L2ErrPres_GS>1e-12) 
      % m
      % [L2ErrDspl_GS, L2ErrPres_GS] 
      old_NumerDspl_GS(:,1) = NumerDisp_GS(:,m+1);
      old_NumerPres_GS(:,1) = NumerPres_GS(:,m+1);
    end 

    if (L2ErrDisp_GS<1e-12 && L2ErrPres_GS<1e-12) 
      % m 
      % [L2ErrDspl_GS, L2ErrPres_GS] 
      NumerDisp(:,n+1) = NumerDisp_GS(:,m+1);
      NumerPres(:,n+1) = NumerPres_GS(:,m+1);
      break; 
    end
  end

  %% Auxiliary quantities 
  % Elementwise 5 coeffs. in WG bas.fxns. 
  cof = zeros(NumEms,5);
  cof(:,1) = NumerPres(1:NumEms,n+1);
  for j=1:4
    cof(:,1+j) = NumerPres(NumEms+QuadriMesh.elem2edge(:,j),n+1);
  end
  % Elementwise cofs. in AC0 nmlz.Piola.bas.: Disc.wk.grad. of WG0 bas.fxns.
  % CDWGB = zeros(NumEms,5,4) 
  CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
  CDWGS = zeros(NumEms,4);
  for j=1:4
    CDWGS(:,j) = sum(cof.*squeeze(CDWGB(:,:,j)),2);
  end
  % GM = zeros(NumEms,4,4);
  % GMK = zeros(NumEms,4,4);
  GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);
  GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh,PermK,GAUSSQUAD);
  % Projection coeffs. in AC0 nmlz.Piola.bas. for permeability * bas.fxn.
  ProjCof = zeros(NumEms,4,4);  % j-th fxn. projected onto k-th fxn. 
  for ie=1:NumEms  % TO BE REVISED FOR EFFECIENCY! 
    EltGM = squeeze(GM(ie,:,:));
    for j=1:4
      EltRHS = squeeze(GMK(ie,:,j));
      ProjCof(ie,j,:) = EltGM\EltRHS';
    end
  end

  %% Numerical velocity: Coeffs. in AC0 elementwise 
  for k=1:4 
    NumerVelCof(:,k,n) = -sum(CDWGS.*squeeze(ProjCof(:,:,k)),2);
  end
end
fprintf('\n');

return;