function [NumerDisp,NumerPres]...
  = LinPoroElas_MS_QuadriIE_CGEQ1_WGP0P0AC0_AsmSlv(...
    EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,QuadriMesh,tyme,GAUSSQUAD)
%% Solving linear poroelasticity by implicit Euler & monolithic system 
% Special ed. for Mandel's problem: Partial Dirichlet boundary conditions 
% James Liu, ColoState; Zhuoran Wang, SYSU; 2017/07--2021/04 

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

%% Time-marching info 
NT = tyme.NT;  MT = tyme.MT;  Deltat = tyme.Deltat;

%% Boundary conditions/data 
% 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)==3);
SymmetryEdgeY     = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==4);
SymmetryEdgeY2    = find(BndryCondTypeElas(QuadriMesh.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 = 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 partial D.B.C., X-direction.
Dirichlet_SymmetryEgsX_NodeElasFlag = zeros(NumNds,1);
for ig=1:NumSymmetryEgsX 
  k1 = QuadriMesh.edge(SymmetryEdgeX(ig),1);
  k2 = QuadriMesh.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 = QuadriMesh.edge(SymmetryEdgeY(ig),1);
  k2 = QuadriMesh.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 = QuadriMesh.edge(SymmetryEdgeY2(ig),1);
  k2 = QuadriMesh.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(QuadriMesh.BndryEdge+1)==1);
NeumannEdgeDarcy   = find(BndryCondTypeDarcy(QuadriMesh.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,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);
% Sampling permeability 
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBIC.fxnK, QuadriMesh ,GAUSSQUAD);
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
% Altogether 
GlbMatDarcy = c0*GlbMatDarcyMass + Deltat*GlbMatDarcyStf;
clear EltMat;

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

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

%% Computing projections/interpolations 
% [ProjPres,IntpltDspl] = LinPoroElas_QuadriMesh_ProjPresP0P0_IntpltDsplBR1(...
%   EqnBIC,QuadriMesh,tyme,GAUSSQUAD);
% Initial conditions: 
% NumerPres(:,1) = ProjPres(:,1);
% NumerDspl(:,1) = IntpltDspl(:,1);
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 solid 
  % Computing elementwise contributions of 12 basis functions 
  EltRHS = 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
      EltRHS(:,j) = EltRHS(:,j) + GAUSSQUAD.RECT(k,3)...
        * jac .*( BFVal(:,j,1).*fval(:,1) + BFVal(:,j,2).*fval(:,2) );
    end  
  end
  % Modification for 4 edge-based basis functions 
  for j=1:4
    EltRHS(:,8+j) = EltRHS(:,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) + EltRHS(ie,2*j-1);
      GlbRhsDisp(2*k  ) = GlbRhsDisp(2*k  ) + EltRHS(ie,2*j  );
      k = QuadriMesh.elem2edge(ie,j);
      GlbRhsDisp(2*NumNds+k) = GlbRhsDisp(2*NumNds+k) + EltRHS(ie,8+j);
    end
  end
  % 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 GlbRHS: Pressure at old time 
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) ...
    + area .* squeeze(NumerPres(1:NumEms,n)) * c0;
  % Assembling GlbRHS: Solid displacement divergence term at old time 
  tmp = zeros(NumEms,12);
  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 = sum(EltMatP0DivBR1.*tmp,2);
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) + EltDispDivOld * alpha;

  % For solid/elas.: Incorporating Neumann (traction) boundary conditions 
  if NumNeumannEgsElas>0
    % For 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 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 * Deltat; 
  end

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

  % For solution 
  sln = zeros(DOFsDisp+DOFsPres,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);
    sln(2*DirichletNodeElas-1) = uDval(:,1);
    sln(2*DirichletNodeElas  ) = uDval(:,2);
    % For edges 
    intgrl1 = 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,:);
    uD1 = EqnBIC.fxnu([x1,y1],tnew);
    uD2 = EqnBIC.fxnu([x2,y2],tnew);
    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);
      intrpl = GAUSSQUAD.LINE(k,1) * uD1 + GAUSSQUAD.LINE(k,2) * uD2;
      res = uD - intrpl;
      psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
      intgrl1 = intgrl1 + GAUSSQUAD.LINE(k,3) * sum(res.*nml,2) .* LenDirichletEdgeElas;
      intgrl2 = intgrl2 + GAUSSQUAD.LINE(k,3) * psi * LenDirichletEdgeElas;
    end
    sln(2*NumNds+DirichletEdgeElas) = intgrl1./intgrl2;
  end

  % 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
    sln(DOFsDisp+NumEms+DirichletEdgeDarcy) = pDavg;
  end
  
  % Reducing...
  GlbRHS = GlbRHS - GlbMat*sln;

  % For flag 
  flag = ones(DOFsDisp+DOFsPres,1);
  % Elas. Dirichlet nodes
  flag(2*DirichletNodeElas-1) = 0;
  flag(2*DirichletNodeElas  ) = 0;
  % Elas. Dirichlet edges 
  flag(2*NumNds+DirichletEdgeElas) = 0;
  % Darcy Dirichlet 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;