function [NumerDisp,NumerPres] = ...
  LinPoroElas_MS_TriIE_CGEP1_WGP0P0RT0_AsmSlv_Mandel( ... 
  EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,TriMesh,tyme,GAUSSQUAD)
%% Solving linear poroelasticity by implicit Euler & monolithic system 
% Zhuoran Wang, SYSU; James Liu, ColoState; 2017/07--2021/01 

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

%% 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;

% For Elas: Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdgeElas = find(BndryCondTypeElas((TriMesh.BndryEdge+1),1)==1);
NeumannEdgeElas   = find(BndryCondTypeElas((TriMesh.BndryEdge+1),1)==2);
SymmetryEdgeX     = find(BndryCondTypeElas((TriMesh.BndryEdge+1),2)==3);
SymmetryEdgeY     = find(BndryCondTypeElas((TriMesh.BndryEdge+1),2)==4);
SymmetryEdgeY2    = find(BndryCondTypeElas((TriMesh.BndryEdge+1),2)==5);
% DirichletEdgeElas = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==1);
% NeumannEdgeElas   = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==2);
% SymmetryEdgeX     = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==3);
% SymmetryEdgeY     = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==4);
% SymmetryEdgeY2    = find(BndryCondTypeElas(TriMesh.BndryEdge+1)==5);
NumDirichletEgsElas = size(DirichletEdgeElas,1);
NumNeumannEgsElas   = size(NeumannEdgeElas,1);
NumSymmetryEgsX     = size(SymmetryEdgeX,1);
NumSymmetryEgsY     = size(SymmetryEdgeY,1);
NumSymmetryEgsY2    = size(SymmetryEdgeY2,1);

% JL20160808: TO BE REVISED FOR EFFICIENCY !
DirichletNodeElasFlag = zeros(NumNds,1);
for ig=1:NumDirichletEgsElas 
  k1 = TriMesh.edge(DirichletEdgeElas(ig),1);
  k2 = TriMesh.edge(DirichletEdgeElas(ig),2);
  DirichletNodeElasFlag(k1) = 1;
  DirichletNodeElasFlag(k2) = 1;
end
DirichletNodeElas = find(DirichletNodeElasFlag);
NumDirichletNdsElas = sum(DirichletNodeElasFlag);
% For partial D.B.C., X-direction.
Dirichlet_SymmetryEgsX_NodeElasFlag = zeros(NumNds,1);
for ig=1:NumSymmetryEgsX 
  k1 = TriMesh.edge(SymmetryEdgeX(ig),1);
  k2 = TriMesh.edge(SymmetryEdgeX(ig),2);
  Dirichlet_SymmetryEgsX_NodeElasFlag(k1) = 1;
  Dirichlet_SymmetryEgsX_NodeElasFlag(k2) = 1;
end
Dirichlet_SymmetryEgsX_NodeElas = find(Dirichlet_SymmetryEgsX_NodeElasFlag);
NumDirichlet_SymmetryEgsX_NodeElas = sum(Dirichlet_SymmetryEgsX_NodeElasFlag);
% For partial D.B.C., Y-direction.
Dirichlet_SymmetryEgsY_NodeElasFlag = zeros(NumNds,1);
for ig=1:NumSymmetryEgsY
  k1 = TriMesh.edge(SymmetryEdgeY(ig),1);
  k2 = TriMesh.edge(SymmetryEdgeY(ig),2);
  Dirichlet_SymmetryEgsY_NodeElasFlag(k1) = 1;
  Dirichlet_SymmetryEgsY_NodeElasFlag(k2) = 1;
end
Dirichlet_SymmetryEgsY_NodeElas = find(Dirichlet_SymmetryEgsY_NodeElasFlag);
NumDirichlet_SymmetryEgsY_NodeElas = sum(Dirichlet_SymmetryEgsY_NodeElasFlag);
% For partial D.B.C., Y-direction-2.
Dirichlet_SymmetryEgsY2_NodeElasFlag = zeros(NumNds,1);
for ig=1:NumSymmetryEgsY2 
  k1 = TriMesh.edge(SymmetryEdgeY2(ig),1);
  k2 = TriMesh.edge(SymmetryEdgeY2(ig),2);
  Dirichlet_SymmetryEgsY2_NodeElasFlag(k1) = 1;
  Dirichlet_SymmetryEgsY2_NodeElasFlag(k2) = 1;
end
Dirichlet_SymmetryEgsY2_NodeElas = find(Dirichlet_SymmetryEgsY2_NodeElasFlag);
NumDirichlet_SymmetryEgsY2_NodeElas = sum(Dirichlet_SymmetryEgsY2_NodeElasFlag);
% For Darcy: Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdgeDarcy = find(BndryCondTypeDarcy(TriMesh.BndryEdge+1)==1);
NeumannEdgeDarcy   = find(BndryCondTypeDarcy(TriMesh.BndryEdge+1)==2);
NumDirichletEgsDarcy = size(DirichletEdgeDarcy,1);
NumNeumannEgsDarcy = size(NeumannEdgeDarcy,1);

%% Setting up the discrete linear system 
DOFsDisp = 2*NumNds + NumEgs;
DOFsPres = NumEms + NumEgs;
NumerDisp = zeros(DOFsDisp,NT+1);
NumerPres = zeros(DOFsPres,NT+1);

%% Assembling GlbMatElas 
% Element matrix: strain-strain 
EltMatStrnStrn = zeros(NumEms,9,9);
NumQuadPts = size(GAUSSQUAD.TRIG,1);
for k=1:NumQuadPts
  [BFGrad,BFDiv] = CG_TriBR1_BasFxnGradDiv(TriMesh,...
    GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
  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);
  for i=1:9
    for j=1:9
      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)) .* area * GAUSSQUAD.TRIG(k,4);
    end
  end
end
% Element matrix: div(avg)-div(avg) 
EltMatDivDiv = zeros(NumEms,9,9);
EltDivAvg = zeros(NumEms,9);
CP1BasGrad = CG_TriP1_BasFxnGrad(TriMesh);
for j=1:3  % Treating vertices and edges in one loop 
  EltDivAvg(:,2*j-1) = CP1BasGrad(:,j,1);
  EltDivAvg(:,2*j  ) = CP1BasGrad(:,j,2);
  EltDivAvg(:,6+j) = TriMesh.LenEg(TriMesh.elem2edge(:,j))./(6*area);
end
for i=1:9
  for j=1:9
    EltMatDivDiv(:,i,j) = EltDivAvg(:,i) .* EltDivAvg(:,j) .* area;
  end
end
% Altogether 
% EltMat = zeros(NumEms,9,9);
EltMat = (2*mu)*EltMatStrnStrn + lambda*EltMatDivDiv;
% Modification for signs of edge normals 
for i=1:3
  for j=1:9
    EltMat(:,6+i,j) = EltMat(:,6+i,j) .* TriMesh.EmEgSign(:,i);
  end
end
for i=1:9
  for j=1:3
    EltMat(:,i,6+j) = EltMat(:,i,6+j) .* TriMesh.EmEgSign(:,j);
  end
end
% Now assembling GlbMatElas 
GlbMatElas = sparse(DOFsDisp,DOFsDisp);
pos = zeros(NumEms,9);
for j=1:3  % Treating vertices and edges in one loop 
  k = TriMesh.elem(:,j);
  pos(:,2*j-1) = 2*k - 1;
  pos(:,2*j  ) = 2*k;
  k = TriMesh.elem2edge(:,j);
  pos(:,6+j) = 2*NumNds + k;
end
for i=1:9 
  II = pos(:,i);
  for j=1:9 
    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);
% Sampling permeability 
PermK = Darcy_SmplnPerm_TriMesh(EqnBIC.fxnK,TriMesh,GAUSSQUAD);
GMK = Hdiv_TriRT0_NmlzBas_GramMatK(TriMesh,PermK);
CDWGB = WG_TriP0P0RT0_CofRT0NmlzBas_DiscWkGradBasFxn(TriMesh);
EltMat = zeros(NumEms,4,4);
% JL20171231: TO BE REVISED FOR EFFECIENCY 
for ie=1:NumEms
  EltGM = squeeze(GMK(ie,:,:));
  for i=1:4
    for j=1:4
      EltMat(ie,i,j) = squeeze(CDWGB(ie,i,:))' *EltGM* squeeze(CDWGB(ie,j,:));
    end
  end
end
pos = zeros(NumEms,4);
pos(:,1) = (1:NumEms);
for j=1:3
  pos(:,1+j) = NumEms + TriMesh.elem2edge(:,j);
end
GlbMatDarcyStf = sparse(DOFsPres,DOFsPres);
for i=1:4
  II = pos(:,i);
  for j=1:4
    JJ = pos(:,j);
    GlbMatDarcyStf = GlbMatDarcyStf ...
      + sparse(II,JJ, EltMat(:,i,j), DOFsPres,DOFsPres);
  end
end
% Altogether 
GlbMatDarcy = c0*GlbMatDarcyMass + Deltat*GlbMatDarcyStf;
clear EltMat;

%% JL20171231: TO BE REVISED 
%% Assembling GlbMatDarcyElas 
% Coupling: P0 (test fxns.) vs. averages of BR1 (trial fxns.) 
EltMatP0DivBR1 = zeros(NumEms,9);
for j=1:3  % For nodes and edges (treated simultaneously) 
  EltMatP0DivBR1(:,2*j-1) = CP1BasGrad(:,j,1).*area;
  EltMatP0DivBR1(:,2*j  ) = CP1BasGrad(:,j,2).*area;
  EltMatP0DivBR1(:,6+j) = TriMesh.EgLen(TriMesh.elem2edge(j))/6 ...
    .* TriMesh.EmEgSign(:,j);
end
GlbMatDarcyElas = sparse(DOFsPres,DOFsDisp);
pos = zeros(NumEms,9);
for j=1:3  % For nodes and edges (treated simultaneously) 
 pos(:,2*j-1) = 2*TriMesh.elem(:,j) - 1;
 pos(:,2*j)   = 2*TriMesh.elem(:,j);
 pos(:,6+j) = 2*NumNds + TriMesh.elem2edge(:,j);
end
II = (1:NumEms);
for j=1:9
  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';

%% MS: The monolithic system 
GlbMat = [GlbMatElas,  GlbMatElasDarcy; ...
          GlbMatDarcyElas, GlbMatDarcy];

%% JL20210116: TEMPORARILY all 0 
%% Computing projections/interpolations 
NumerPres(:,1) = 0;
NumerDisp(:,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);

  % JL20171221: TO BE REVISED 
  % Assembling GlbRHS: For displacement 
  % Computing elementwise contributions of 12 basis functions 
  EltRHS = 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),:);
    fval = EqnBIC.fxnf(qp,tnew);
    BFVal = CG_TriBR1_BasFxnVal(TriMesh,...
      GAUSSQUAD.TRIG(k,1), GAUSSQUAD.TRIG(k,2), GAUSSQUAD.TRIG(k,3));
    for j=1:9
      EltRHS(:,j) = EltRHS(:,j) + GAUSSQUAD.TRIG(k,4)...
        *( BFVal(:,j,1).*fval(:,1) + BFVal(:,j,2).*fval(:,2) ).*area;
    end
  end
  % Modification for 4 edge-based basis functions 
  for j=1:3
    EltRHS(:,6+j) = EltRHS(:,6+j) .* TriMesh.EmEgSign(:,j);
  end
  % Now assembly 
  GlbRhsDisp = zeros(DOFsDisp,1);
  for ie=1:NumEms
    for j=1:3  % Treating vertices and edges in one loop 
      k = TriMesh.elem(ie,j);
      GlbRhsDisp(2*k-1) = GlbRhsDisp(2*k-1) + EltRHS(ie,2*j-1);
      GlbRhsDisp(2*k  ) = GlbRhsDisp(2*k  ) + EltRHS(ie,2*j  );
      k = TriMesh.elem2edge(ie,j);
      GlbRhsDisp(2*NumNds+k) = GlbRhsDisp(2*NumNds+k) + EltRHS(ie,6+j);
    end
  end
  
  % Assembling GlbRHS: For pressure 
  GlbRhsPres = zeros(DOFsPres,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),:);
    sval = EqnBIC.fxns(qp,tnew);
    GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms)...
      + GAUSSQUAD.TRIG(k,4) * sval.* area * Deltat;
  end
  
  % Assembling GlbRHS: Pressure at old time 
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) ...
    + c0 * area.*squeeze(NumerPres(1:NumEms,n));

  % Assembling GlbRHS: Solid displacement divergence term at old time 
  tmp = zeros(NumEms,9);
  for j=1:3  % Treating vertices and edges in one loop 
    k = TriMesh.elem(:,j);
    tmp(:,2*j-1) = NumerDisp(2*k-1,n);
    tmp(:,2*j  ) = NumerDisp(2*k,  n);
    k = TriMesh.elem2edge(:,j);
    tmp(:,6+j) = NumerDisp(2*NumNds+k,n);
  end
  EltDispDivOld = sum(squeeze(EltMatP0DivBR1).*tmp,2);
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) + alpha * EltDispDivOld;

  % For Elas.: Incorporating Neumann (traction) boundary conditions 
  if NumNeumannEgsElas>0
    % For end-nodes 
    NeumannBC = zeros(NumNeumannEgsElas,2,2);
    j1 = TriMesh.edge(NeumannEdgeElas,1);  
    j2 = TriMesh.edge(NeumannEdgeElas,2);
    x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2);
    x2 = TriMesh.node(j2,1);  y2 = TriMesh.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 = TriMesh.edge(NeumannEdgeElas(ig),1);
      j2 = TriMesh.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 = TriMesh.edge(NeumannEdgeElas,1);
    j2 = TriMesh.edge(NeumannEdgeElas,2);
    x1 = TriMesh.node(j1,1);  y1 = TriMesh.node(j1,2);
    x2 = TriMesh.node(j2,1);  y2 = TriMesh.node(j2,2);
    LenEgNeumann = sqrt((x2-x1).^2+(y2-y1).^2);
    nml = TriMesh.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 Darcy: Incorporating Neumann (flux) boundary conditions 
  if (NumNeumannEgsDarcy>0)
    NeumannBC = zeros(NumNeumannEgsDarcy,1);
    x1 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,1),1);
    y1 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,1),2);
    x2 = TriMesh.node(TriMesh.edge(NeumannEdgeDarcy,2),1);
    y2 = TriMesh.node(TriMesh.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

  % Assembling GlbRHS: Combination 
  GlbRHS = [GlbRhsDisp; GlbRhsPres];

  % For solution 
  sln = zeros(DOFsDisp+DOFsPres,1);

  % Enforcing partial Dirichlet boundary conditions (x-direction) 
  if (NumSymmetryEgsX>0) 
    % For nodes: 
    Dirichlet_SymmetryEgsX_NdElas = TriMesh.node(Dirichlet_SymmetryEgsX_NodeElas,:);
    % upDval = EqnBIC.fxnuDpartial(Dirichlet_SymmetryEgsX_NdElas,tnew,3);
    sln(2*Dirichlet_SymmetryEgsX_NodeElas-1) = 0;
  end  
  if (NumSymmetryEgsX>0)
    sln(2*NumNds+SymmetryEdgeX) = 0;
  end
  
  % Enforcing partial Dirichlet boundary conditions (y-direction)   
  if (NumSymmetryEgsY>0) 
    % For nodes: 
    Dirichlet_SymmetryEgsY_NdElas = TriMesh.node(Dirichlet_SymmetryEgsY_NodeElas,:);
    upDval = EqnBIC.fxnuDpartial(Dirichlet_SymmetryEgsY_NdElas,tnew,4);
    %  sln(2*Dirichlet_SymmetryEgsY_NodeElas) = upDval(:,2);
    sln(2*Dirichlet_SymmetryEgsY_NodeElas) = 0;
  end
  if (NumSymmetryEgsY>0)
    sln(2*NumNds+SymmetryEdgeY) = 0;
  end
  
  % Enforcing partial Dirichlet boundary conditions (y-direction nonzero)   
  if (NumSymmetryEgsY2>0) 
    % For nodes: 
    Dirichlet_SymmetryEgsY2_NdElas = TriMesh.node(Dirichlet_SymmetryEgsY2_NodeElas,:);
    upDval = EqnBIC.fxnuDpartial(Dirichlet_SymmetryEgsY2_NdElas,tnew,5);
    sln(2*Dirichlet_SymmetryEgsY2_NodeElas) = upDval(:,2);
    % sln(2*Dirichlet_SymmetryEgsY2_NodeElas) = 1;
  end
  if (NumSymmetryEgsY2>0)
    sln(2*NumNds+SymmetryEdgeY2) = 0;
  end

  % For Elas.: Enforcing Dirichlet boundary conditions as essential 
  if (NumDirichletEgsElas>0)  % Assuming NumDirichletEgsElas>0 
    % For nodes: 
    DirichletNdElas = TriMesh.node(DirichletNodeElas,:);
    uDval = EqnBIC.fxnuD(DirichletNdElas,tnew);
    sln(2*DirichletNodeElas-1) = uDval(:,1);
    sln(2*DirichletNodeElas  ) = uDval(:,2);
  end
  if (NumDirichletEgsElas>0)
    intgrl1 = zeros(NumDirichletEgsElas,1);
    intrpl = zeros(NumDirichletEgsElas,1);
    intgrl2 = zeros(NumDirichletEgsElas,1);
    x1 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,1),1);
    y1 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,1),2);
    x2 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,2),1);
    y2 = TriMesh.node(TriMesh.edge(DirichletEdgeElas,2),2);
    LenDirichletEdgeElas = sqrt((x2-x1).^2+(y2-y1).^2);
    nml = TriMesh.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(TriMesh.node(TriMesh.edge(DirichletEdgeElas,1),:),tnew)...
              + GAUSSQUAD.LINE(k,2)*EqnBIC.fxnu(TriMesh.node(TriMesh.edge(DirichletEdgeElas,2),:),tnew); 
      res = uD - intrpl;
      intgrl2 = intgrl2 + GAUSSQUAD.LINE(k,3)*sum(res.*nml,2).*LenDirichletEdgeElas;
    end
    sln(2*NumNds+DirichletEdgeElas) = intgrl2./intgrl1;
  end
    
  % For Darcy: Enforcing Dirichlet boundary conditions as essential 
  if (NumDirichletEgsDarcy>0)  % Assuming NumDirichletEgsDarcy>0 
    pDavg = zeros(NumDirichletEgsDarcy,1);
    x1 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,1),1);
    y1 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,1),2);
    x2 = TriMesh.node(TriMesh.edge(DirichletEdgeDarcy,2),1);
    y2 = TriMesh.node(TriMesh.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
    sln(DOFsDisp+NumEms+DirichletEdgeDarcy) = pDavg;
  end
  
  % Reducing...
  GlbRHS = GlbRHS - GlbMat*sln;

  % For flag 
  flag = ones(DOFsDisp+DOFsPres,1);
  % Elas D.B.C.nodes
  flag(2*DirichletNodeElas-1) = 0;
  flag(2*DirichletNodeElas  ) = 0;
  % Elas D.B.C.edges
  flag(2*NumNds+DirichletEdgeElas) = 0;
  % Elas partial D.B.C. edges
  flag(2*NumNds+SymmetryEdgeX) = 0;
  flag(2*NumNds+SymmetryEdgeY) = 0;
  flag(2*NumNds+SymmetryEdgeY2) = 0;  % For y = 1
  % Elas partial D.B.C. nodes
  flag(2*Dirichlet_SymmetryEgsX_NodeElas-1) = 0;
  flag(2*Dirichlet_SymmetryEgsY_NodeElas) = 0;
  flag(2*Dirichlet_SymmetryEgsY2_NodeElas) = 0;  % For y = 1 
  % Darcy D.B.C edges, 
  flag(DOFsDisp+NumEms+DirichletEdgeDarcy) = 0;
  FreeNdEmFreeEg = find(flag);
  
  % Solving the reduced linear system 
  sln(FreeNdEmFreeEg) = GlbMat(FreeNdEmFreeEg,FreeNdEmFreeEg) \ GlbRHS(FreeNdEmFreeEg);

  % Extracting and storing 
  NumerDisp(:,n+1) = sln(1:DOFsDisp);
  NumerPres(:,n+1) = sln(DOFsDisp+(1:DOFsPres));
end
fprintf('\n');

return;