function [NumerDisp, NumerPres] = ...
  LinPoroElas_MS_QuadriIE_WGP02P02AC02P0_WGP0P0AC0_AsmSlv( ...
  EqnBIC,BndryCondTypeElas,BndryCondTypeDarcy,QuadriMesh,tyme,GAUSSQUAD) 
%% Solving linear poroelasticity by implicit Euler & monolithic system 
% James Liu, ColoState; Zhuoran Wang, SYSU; 2020/09--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;

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

%% Sorting out boundary edges 
% For Elas.: Dirichlet, Neumann 
DirichletEdgeElas = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==1);
NeumannEdgeElas   = find(BndryCondTypeElas(QuadriMesh.BndryEdge+1)==2);
NumDirichletEgsElas = size(DirichletEdgeElas,1);
NumNeumannEgsElas = size(NeumannEdgeElas,1);
% 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*NumEms + 2*NumEgs;
DOFsPres = NumEms + NumEgs;
NumerDisp = zeros(DOFsDisp, NT+1);
NumerPres = zeros(DOFsPres, NT+1);

%% Elas.: Assembling GlbMatElas 
% Computing elementwise matrices: strain-div formulation 
% GMAC02,GMST both size NumEms*8*8;  
% CDWGB size NumEms*10*8;  CDWDB size NumEms*10*1; 
GMAC02 = Hdiv_QuadriAC02_NmlzPiolaBas_GramMat(QuadriMesh,GAUSSQUAD);
GMST = Hdiv_QuadriAC02_NmlzPiolaBas_GramMatST(QuadriMesh,GAUSSQUAD);
CDWGB = WG_QuadriP02P02AC02_CofNmlzPiolaBas_DiscWkGradBasFxn(QuadriMesh,GAUSSQUAD);
EltMatStrnStrn = zeros(NumEms,10,10);
for i=1:10
  for j=1:10
    EltMatStrnStrn(:,i,j) = 0;
    for k=1:8 
      for l=1:8 
        EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j) ...        
          + CDWGB(:,i,k) .* (0.5*(GMAC02(:,k,l)+GMST(:,k,l))) .* CDWGB(:,j,l);
      end
    end
  end
end
CDWDB = WG_QuadriP02P02P0_CofNmlzBas_DiscWkDivBasFxn(QuadriMesh);
EltMatDivDiv = zeros(NumEms,10,10);
for i=1:10
  for j=1:10
    EltMatDivDiv(:,i,j) = CDWDB(:,i) .* CDWDB(:,j) .* area;
  end
end
% Elementwise matrices altogether: EltMat size NumEms*10*10; 
% JL20170923: FOR CONVENIENCE NOW, BUT USE TOO MUCH MEMORY 
EltMat = (2*mu) * EltMatStrnStrn + lambda * EltMatDivDiv;
% Positions for assembly for elasticity 
pos = zeros(NumEms,10);
% For element interiors 
pos(:,1) = 2*(0:NumEms-1) + 1;
pos(:,2) = 2*(0:NumEms-1) + 2;
% For elementwise 4 edges 
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  pos(:,2+2*j-1) = 2*NumEms + 2*k-1;
  pos(:,2+2*j  ) = 2*NumEms + 2*k;
end
% Now assembly 
GlbMatElas = sparse(DOFsDisp,DOFsDisp);
for i=1:10
  II = pos(:,i); 
  for j=1:10
    JJ = pos(:,j);
    GlbMatElas = GlbMatElas ... 
      + sparse(II,JJ, EltMat(:,i,j), DOFsDisp,DOFsDisp);
  end
end

%% Darcy: Assembling GlbMatDarcy 
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 EFFICIENCY 
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
% GlbMatDarcyMass = sparse(DOFsPres,DOFsPres);
II = (1:NumEms);  JJ = II;
GlbMatDarcyMass = sparse(II,JJ, area, DOFsPres,DOFsPres);
% Positions for assembly for Darcy "stiffness" 
pos = zeros(NumEms,5);
pos(:,1) = (1:NumEms);
for j=1:4
  pos(:,1+j) = NumEms + QuadriMesh.elem2edge(:,j);
end
% Now assembly 
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;

%% Darcy-Elas.: Assembling GlbMatDarcyElas 
% Coupling: Darcy P0 test fxns. with Elas. P0^2 trial fxns. 
EltMat = zeros(NumEms,10);  % Emws. 1 P0 (intr.)trial vs. 10 P0^2 test 
for i=1:10 
  EltMat(:,i) = CDWDB(:,i) .* area;
end
% Positions for assembly for GlbMatDarcyElas 
pos = zeros(NumEms,10);
pos(:,1) = 2*(0:NumEms-1) + 1;
pos(:,2) = 2*(0:NumEms-1) + 2;
for j=1:4
  k = QuadriMesh.elem2edge(:,j);
  pos(:,2+2*j-1) = 2*NumEms + 2*k-1;
  pos(:,2+2*j  ) = 2*NumEms + 2*k;
end
% Now assembly for GlbMatDarcyElas 
GlbMatDarcyElas = sparse(DOFsPres,DOFsDisp);
II = (1:NumEms); 
for j=1:10
  JJ = pos(:,j);
  GlbMatDarcyElas = GlbMatDarcyElas ... 
    + sparse(II,JJ, EltMat(:,j), DOFsPres,DOFsDisp);
end

%% Elas.-Darcy: Assembling/Producing GlbMatElasDarcy 
% Coupling: Elas. P0^2 test fxns. with Darcy P0 trial fxns. 
GlbMatElasDarcy = -GlbMatDarcyElas'; 

%% MS: A big monolithic system 
GlbMat = [GlbMatElas,  alpha * GlbMatElasDarcy; ...
          alpha * GlbMatDarcyElas, GlbMatDarcy];

%% JL20201225: TEMPORARILY all set to 0 
%% Setting up initial conditions 
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);

  % Assembling GlbRHS: For Elas. body force 
  GlbRhsDisp = zeros(DOFsDisp,1);
  II1 = 2*(0:NumEms-1) + 1;
  II2 = 2*(0:NumEms-1) + 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);
    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);
    GlbRhsDisp(II1) = GlbRhsDisp(II1) + GAUSSQUAD.RECT(k,3) * (jac.*fval(:,1));
    GlbRhsDisp(II2) = GlbRhsDisp(II2) + GAUSSQUAD.RECT(k,3) * (jac.*fval(:,2));
  end
  
  % Assembling GlbRHS: For Darcy source term 
  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) * (jac.*sval) * Deltat;
  end

  % Assembling GlbRHS: Fluid pressure at old time 
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) ...
    + area.*squeeze(NumerPres(1:NumEms,n)) * c0;

  % Assembling GlbRHS: Solid dilation at old time 
  tmp = zeros(NumEms,10);
  tmp(:,1) = NumerDisp(2*(0:NumEms-1)+1, n);
  tmp(:,2) = NumerDisp(2*(0:NumEms-1)+2, n);
  for j=1:4
    k = QuadriMesh.elem2edge(:,j);
    tmp(:,2+2*j-1) = NumerDisp(2*NumEms+2*k-1,n);
    tmp(:,2+2*j  ) = NumerDisp(2*NumEms+2*k,  n);
  end
  EltDispDivOld = sum(CDWDB.*tmp,2) .* area;
  GlbRhsPres(1:NumEms) = GlbRhsPres(1:NumEms) + EltDispDivOld * alpha;

  % For Elas.: Incorporating Neumann (traction) boundary conditions 
  if (NumNeumannEgsElas>0)
    NeumannBC = zeros(NumNeumannEgsElas,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 = NeumannBC + GAUSSQUAD.LINE(k,3) * tNval;
    end
    NeumannBC = NeumannBC .* [LenEgNeumann, LenEgNeumann];
    GlbRhsDisp(2*NumEms+NeumannEdgeElas) = ...
      GlbRhsDisp(2*NumEms+NeumannEdgeElas) + NeumannBC;
  end

  % For 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  
    uDavg = zeros(NumDirichletEgsElas,2);
    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);
    NumQuadPts = size(GAUSSQUAD.LINE,1);
    for k=1:NumQuadPts
      qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
      uDval = EqnBIC.fxnuD(qp,tnew);
      uDavg = uDavg + GAUSSQUAD.LINE(k,3) * uDval;
    end
    sln(2*NumEms+2*DirichletEdgeElas-1) = uDavg(:,1);
    sln(2*NumEms+2*DirichletEdgeElas  ) = uDavg(:,2);
  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);
  % For Elas. Dirichlet boundary conditions 
  flag(2*NumEms+2*DirichletEdgeElas-1) = 0;
  flag(2*NumEms+2*DirichletEdgeElas  ) = 0;
  % For Darcy Dirichlet boundary conditions: 
  flag(DOFsDisp+NumEms+DirichletEdgeDarcy) = 0;

  % Solving the reduced linear system 
  FreeEmEg = find(flag);
  sln(FreeEmEg) = GlbMat(FreeEmEg,FreeEmEg) \ GlbRHS(FreeEmEg);

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

return;