function [sln,CofQhg] = StokesDarcy_QuadriCGBR1WGAC0_1meshBT_AsmSlv( ... 
  EqnBC, QuadriMesh, BndryCondType, GAUSSQUAD)
% The 1-(quadri.)mesh approach with Darcy at bottom and Stokes at top 
% -- Labels for mesh components resp. for Stokes and Darcy; 
% -- Assembly for Stokes, Darcy, Interface; 
% -- Enforcing bndry.conds. for Stokes and Darcy; 
% -- Solving the reduced monolithic discrete system; 
% For Darcy: Assuming PermK is an elementwise constant 2x2 SPD matrix 
% Graham Harper, Sandia, James Liu, ColoState; 2019/07--2020/08 

%% Equation info 
   nu = EqnBC.nu;
alpha = EqnBC.alpha;

%% Mesh info 
% Assuming QuadriMesh.flag>=2 
CofA = QuadriMesh.CofA;
CofB = QuadriMesh.CofB;

%% More mesh info 
LblEmStokes = QuadriMesh.StokesEms;
LblEmDarcy  = QuadriMesh.DarcyEms;
%
% Assuming vertical concatenation of Stokes and Darcy 
NumNdsStokes = nnz(QuadriMesh.posStokesNd);
NumEgsStokes = nnz(QuadriMesh.posStokesEg);
NumEmsStokes = length(LblEmStokes);
NumDOFsStokes = 2*NumNdsStokes + NumEgsStokes + NumEmsStokes;
%
NumEgsDarcy = nnz(QuadriMesh.posDarcyEg);
NumEmsDarcy = length(LblEmDarcy);
NumDOFsDarcy = NumEmsDarcy + NumEgsDarcy;
%
NumNdsInterface = length(QuadriMesh.InterfaceNds);
NumEgsInterface = length(QuadriMesh.InterfaceEgs);
%
DOFsDarcy  = (1:NumDOFsDarcy);
DOFsStokes = DOFsDarcy + (1:NumDOFsStokes)';

%% Sorting out boundary edges: Dirichlet, Neumann 
StokesDirichletEdge = find(QuadriMesh.BndryEdge>=3 & QuadriMesh.BndryEdge<=5 & BndryCondType(QuadriMesh.BndryEdge+1)==1);
StokesNeumannEdge   = find(QuadriMesh.BndryEdge>=3 & QuadriMesh.BndryEdge<=5 & BndryCondType(QuadriMesh.BndryEdge+1)==2);
NumDirichletEgsStokes = length(StokesDirichletEdge);
NumNeumannEgsStokes   = length(StokesNeumannEdge);
% StokesDirichletEdge'
% StokesNeumannEdge'
DarcyDirichletEdge = find((QuadriMesh.BndryEdge<=2 | QuadriMesh.BndryEdge>=6) & BndryCondType(QuadriMesh.BndryEdge+1)==1);
DarcyNeumannEdge   = find((QuadriMesh.BndryEdge<=2 | QuadriMesh.BndryEdge>=6) & BndryCondType(QuadriMesh.BndryEdge+1)==2);
NumDirichletEgsDarcy = length(DarcyDirichletEdge);
NumNeumannEgsDarcy   = length(DarcyNeumannEdge);
% DarcyDirichletEdge'
% DarcyNeumannEdge'

%% JL20160808: TO BE REVISED FOR EFFICIENCY !
StokesDirichletNodeFlag = zeros(QuadriMesh.NumNds,1);
for ig=1:NumDirichletEgsStokes
  j1 = QuadriMesh.edge(StokesDirichletEdge(ig),1);
  j2 = QuadriMesh.edge(StokesDirichletEdge(ig),2);
  StokesDirichletNodeFlag(j1) = 1;
  StokesDirichletNodeFlag(j2) = 1;
end
StokesDirichletNode = find(StokesDirichletNodeFlag);

%% Setup for the gigantic monolithic discrete linear system 
DOFs = NumDOFsDarcy + NumDOFsStokes;
GlbMat = sparse(DOFs,DOFs);
GlbRHS = zeros(DOFs,1);

%% JL20191107: UNDER THE ASSUMPTION "Darcy elements go first" 
%% Sampling the permeability (2x2 SPD matrix)
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBC.fxnK_Darcy, QuadriMesh, GAUSSQUAD);

%% Assembly of GlbMat: For Darcy
% Darcy: First computing element stiffness matrices 
% GMK = zeros(NumEms,4,4);
GM  = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat( QuadriMesh, GAUSSQUAD);
GMK = Hdiv_QuadriAC0_NmlzPiolaBas_GramMatK(QuadriMesh, PermK, GAUSSQUAD);
% CDWGB = zeros(NumEms,5,4);
CDWGB = WG_QuadriP0P0AC0_CofAC0NmlzPiolaBas_DiscWkGradBasFxn(... 
  QuadriMesh, GAUSSQUAD);
ArrayEE = zeros(NumEmsDarcy,1);
ArrayEG = zeros(NumEmsDarcy,4);
ArrayGG = zeros(NumEmsDarcy,4,4);
%% JL20160517: TO BE REVISED FOR EFFECIENCY !!
for ie=1:NumEmsDarcy
  EGM = squeeze(GMK(LblEmDarcy(ie),:,:));
  ArrayEE(ie) = squeeze(CDWGB(LblEmDarcy(ie),1,:))' * EGM * squeeze(CDWGB(LblEmDarcy(ie),1,:));
  for i=1:4
    ArrayEG(ie,i) = squeeze(CDWGB(LblEmDarcy(ie),1,:))' * EGM * squeeze(CDWGB(LblEmDarcy(ie),1+i,:));
    for j=1:4
      ArrayGG(ie,i,j) = squeeze(CDWGB(LblEmDarcy(ie),1+i,:))' * EGM * squeeze(CDWGB(LblEmDarcy(ie),1+j,:));
    end
  end
end
% Darcy: Secondly assembling element stiffness matrices 
% Part 1: element-element interaction 
II = (1:NumEmsDarcy);  JJ = II;
GlbMat = GlbMat + sparse(II,JJ,ArrayEE,DOFs,DOFs);
% Part 2: element-edge interaction 
II = (1:NumEmsDarcy);
for j=1:4  % Looping over edges per element 
  JJ = NumEmsDarcy + QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(LblEmDarcy,j));
  GlbMat = GlbMat + sparse(II,JJ,ArrayEG(:,j),DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II,ArrayEG(:,j),DOFs,DOFs);
end
% Part 3: edge-edge interaction 
for i=1:4
  II = NumEmsDarcy + QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(LblEmDarcy,i));
  for j=i:4  % Utilizing symmetry 
    JJ = NumEmsDarcy + QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(LblEmDarcy,j));
    if (j==i) 
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
    else
      GlbMat = GlbMat + sparse(II,JJ,ArrayGG(:,i,j),DOFs,DOFs);
      GlbMat = GlbMat + sparse(JJ,II,ArrayGG(:,i,j),DOFs,DOFs);
    end
  end
end

%% Assembly of GlbMat: For Stokes 
% Stokes: First computing element stiffness matrices 
% For velocity strain-strain interaction 
EltMatStrnStrn = zeros(NumEmsStokes,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(LblEmStokes,3)+CofB(LblEmStokes,4)*xhat) .* (CofA(LblEmStokes,2)+CofA(LblEmStokes,4)*yhat)...
      - (CofA(LblEmStokes,3)+CofA(LblEmStokes,4)*xhat) .* (CofB(LblEmStokes,2)+CofB(LblEmStokes,4)*yhat);
  [BFGrad,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  StrnPhi = zeros(NumEmsStokes,12,2,2);
  StrnPhi(:,:,1,1) = BFGrad(LblEmStokes,:,1,1);
  StrnPhi(:,:,2,2) = BFGrad(LblEmStokes,:,2,2);
  StrnPhi(:,:,1,2) = 0.5*(BFGrad(LblEmStokes,:,1,2)+BFGrad(LblEmStokes,:,2,1));
  StrnPhi(:,:,2,1) = StrnPhi(:,:,1,2);
  for i=1:12
    for j=1:12
      EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j) ...
        + (2*nu) *( 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
% Modification of signs for edge-bubble functions 
for i=1:12
  for j=9:12    
    EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j) .* QuadriMesh.EmEgSign(LblEmStokes,j-8);
  end
end
for i=9:12
  for j=1:12
    EltMatStrnStrn(:,i,j) = EltMatStrnStrn(:,i,j) .* QuadriMesh.EmEgSign(LblEmStokes,i-8);
  end
end
% For pressure & velocity div 
EltMatPresDiv = zeros(NumEmsStokes,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(LblEmStokes,3)+CofB(LblEmStokes,4)*xhat) .* (CofA(LblEmStokes,2)+CofA(LblEmStokes,4)*yhat)...
      - (CofA(LblEmStokes,3)+CofA(LblEmStokes,4)*xhat) .* (CofB(LblEmStokes,2)+CofB(LblEmStokes,4)*yhat);
  [~,BFDiv] = CG_QuadriBR1_BasFxnGradDiv(QuadriMesh,xhat,yhat);
  for i=1:12
    EltMatPresDiv(:,i) = EltMatPresDiv(:,i) ...
      + BFDiv(LblEmStokes,i) .* jac * GAUSSQUAD.RECT(k,3);
  end
end
% Modification of signs for edge-bubble functions 
for i=9:12
  EltMatPresDiv(:,i) = EltMatPresDiv(:,i) .* QuadriMesh.EmEgSign(LblEmStokes,i-8);
end
% 
% Stokes: Secondly assembling element stiffness matrices 
%% JL20180224: THIS IS RISKY 
% For velocity grad-grad 
II = zeros(NumEmsStokes,12);
for i=1:4
  II(:,2*i-1) = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(QuadriMesh.elem(LblEmStokes,i)) - 1;
  II(:,2*i  ) = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(QuadriMesh.elem(LblEmStokes,i));
  II(:,8+i) =   NumDOFsDarcy + 2*NumNdsStokes ... 
    + QuadriMesh.posStokesEg(QuadriMesh.elem2edge(LblEmStokes,i));
end
JJ = II;
for i=1:12
  for j=1:12 
    GlbMat = GlbMat + sparse(II(:,i),JJ(:,j), EltMatStrnStrn(:,i,j), DOFs,DOFs);
  end
end
% For pressure & velocity div
JJ = NumDOFsDarcy + 2*NumNdsStokes + NumEgsStokes + (1:NumEmsStokes);
for i=1:12
  GlbMat = GlbMat + sparse(II(:,i),JJ, -EltMatPresDiv(:,i), DOFs,DOFs);
  GlbMat = GlbMat + sparse(JJ,II(:,i),  EltMatPresDiv(:,i), DOFs,DOFs);
end

% QuadriMesh.InterfaceEgs'

%% Assembly of GlbMat: For interface interaction 
% Computing and assembling the two matrices simultaneously
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig=(QuadriMesh.InterfaceEgs')
  x1 = QuadriMesh.node(QuadriMesh.edge(ig,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(ig,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(ig,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(ig,2),2);
  II = NumEmsDarcy + QuadriMesh.posDarcyEg(ig);
  JJ = zeros(1,5);
  for j=1:2  % Looping over vertices on the edge 
    JJ(2*j-1) = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(QuadriMesh.edge(ig,j)) - 1;
    JJ(2*j  ) = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(QuadriMesh.edge(ig,j));
  end
  JJ(5) = NumDOFsDarcy + 2*NumNdsStokes + QuadriMesh.posStokesEg(ig);
  tangent = [x2-x1,y2-y1]/QuadriMesh.EgLen(ig);  % Unit tangent vector 
  normal = [tangent(2), -tangent(1)]; 
  LocalMat = zeros(1,5);
  for k=1:NumQuadPts
    shat = GAUSSQUAD.LINE(k,2);
    StokesVelShapeVal = [normal * ((1-shat)*[1;0]), ...  % Bas.fxn. 1
                         normal * ((1-shat)*[0;1]), ...  % Bas.fxn. 2
                         normal * (    shat*[1;0]), ...  % Bas.fxn. 3
                         normal * (    shat*[0;1]), ...  % Bas.fxn. 4
                                  (1-shat) * shat];      % Bas.fxn. 5 
    LocalMat = LocalMat + StokesVelShapeVal * GAUSSQUAD.LINE(k,3) * QuadriMesh.EgLen(ig);
  end
  for j=1:5  % Only 4 basis functions are used, Why? (JL20200718) 
    GlbMat = GlbMat + sparse(II,JJ(j), -LocalMat(j), DOFs,DOFs);
    GlbMat = GlbMat + sparse(JJ(j),II,  LocalMat(j), DOFs,DOFs);
  end
end

%% Assembly of GlbMat: For BJS interface condition 
for ig=(QuadriMesh.InterfaceEgs')
  % Note: Edge bubble functions do not contribute to BJS condition! 
  % instead formulate this as a local 4x4 problem  
  x1 = QuadriMesh.node(QuadriMesh.edge(ig,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(ig,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(ig,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(ig,2),2);
  tangent = [x2-x1,y2-y1]/QuadriMesh.EgLen(ig);  % Unit tangent vector 
  II = zeros(1,4);
  for j=1:2  % Two vertices
    II(2*j-1) = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(QuadriMesh.edge(ig,j)) - 1;
    II(2*j  ) = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(QuadriMesh.edge(ig,j));
  end
  % II = repmat(II,4,1);
  JJ = II;  
  LocalMat = zeros(4,4);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    shat = GAUSSQUAD.LINE(k,2);
    StokesVelShapeVal = [tangent*((1-shat)*[1;0]), ... % shape fxn 1
                         tangent*((1-shat)*[0;1]), ... % shape fxn 2
                         tangent*(   shat *[1;0]), ... % shape fxn 3
                         tangent*(   shat *[0;1])];    % shape fxn 4
    % Grahma's Notes: No need for shape fxn. #5 since tangent*normal = 0 
    % JL20200731: sqrt(1) TO BE REVISED !! 
    LocalMat = LocalMat + StokesVelShapeVal' * StokesVelShapeVal * ...
      GAUSSQUAD.LINE(k,3) * QuadriMesh.EgLen(ig) * alpha*nu/sqrt(1);

  end
  for i=1:4
    for j=1:4
      GlbMat = GlbMat + sparse(II(i),JJ(j), LocalMat(i,j), DOFs,DOFs);
    end
  end
end
% size(II)

%% Auxiliary but needed stuff 
%% Darcy: Computing local L2-projection of g=-(K/nu)f into elementwise AC0 
CofQhg = zeros(NumEmsDarcy,4);
EltRHSg = zeros(NumEmsDarcy,4);
MatJ = zeros(NumEmsDarcy,2,2);
MatP = zeros(NumEmsDarcy,2,2);
qp = zeros(NumEmsDarcy,2);
for k=1:NumQuadPts
  xhat = GAUSSQUADRECT1(k,1);
  yhat = GAUSSQUADRECT1(k,2);
  J11 = CofA(LblEmDarcy,2) + CofA(LblEmDarcy,4)*yhat;  
  J12 = CofA(LblEmDarcy,3) + CofA(LblEmDarcy,4)*xhat;
  J21 = CofB(LblEmDarcy,2) + CofB(LblEmDarcy,4)*yhat;  
  J22 = CofB(LblEmDarcy,3) + CofB(LblEmDarcy,4)*xhat;
  jac = J11.*J22 - J12.*J21;
  MatJ(:,1,1) = J11;  MatJ(:,1,2) = J12;
  MatJ(:,2,1) = J21;  MatJ(:,2,2) = J22;
  for i=1:2
    for j=1:2
      MatP(:,i,j) = MatJ(:,i,j)./jac;
    end
  end
  qp(:,1) = CofA(LblEmDarcy,1) + CofA(LblEmDarcy,2)*xhat + CofA(LblEmDarcy,3)*yhat + CofA(LblEmDarcy,4)*xhat*yhat;
  qp(:,2) = CofB(LblEmDarcy,1) + CofB(LblEmDarcy,2)*xhat + CofB(LblEmDarcy,3)*yhat + CofB(LblEmDarcy,4)*xhat*yhat;
  X = qp(:,1) - QuadriMesh.EmCntr(LblEmDarcy,1);
  Y = qp(:,2) - QuadriMesh.EmCntr(LblEmDarcy,2);
  U = MatP(:,1,1) .* xhat - MatP(:,1,2) .* yhat;
  V = MatP(:,2,1) .* xhat - MatP(:,2,2) .* yhat;
  gval = EqnBC.fxng_Darcy(qp);
  EltRHSg(:,1) = EltRHSg(:,1) +  gval(:,1) .* jac *  GAUSSQUAD.RECT(k,3);
  EltRHSg(:,2) = EltRHSg(:,2) +  gval(:,2) .* jac *  GAUSSQUAD.RECT(k,3);
  EltRHSg(:,3) = EltRHSg(:,3) + (gval(:,1).*X + gval(:,2).*Y) .* jac *  GAUSSQUAD.RECT(k,3);
  EltRHSg(:,4) = EltRHSg(:,4) + (gval(:,1).*U + gval(:,2).*V) .* jac *  GAUSSQUAD.RECT(k,3);
end
for ie=1:NumEmsDarcy
  CofQhg(ie,:) = ... 
    squeeze(GM(LblEmDarcy(ie),:,:)) \ squeeze(EltRHSg(ie,:))';
end

%% Assembly of GlbRHS: 
GlbRHS = zeros(DOFs,1);
% For Darcy: Source (or sink)  
qp = zeros(NumEmsDarcy,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);
  qp(1:NumEmsDarcy,1) = CofA(LblEmDarcy,1) + CofA(LblEmDarcy,2)*xhat + CofA(LblEmDarcy,3)*yhat + CofA(LblEmDarcy,4)*xhat*yhat;
  qp(1:NumEmsDarcy,2) = CofB(LblEmDarcy,1) + CofB(LblEmDarcy,2)*xhat + CofB(LblEmDarcy,3)*yhat + CofB(LblEmDarcy,4)*xhat*yhat;
  jac = (CofB(LblEmDarcy,3)+CofB(LblEmDarcy,4)*xhat) .* (CofA(LblEmDarcy,2)+CofA(LblEmDarcy,4)*yhat) ...
      - (CofA(LblEmDarcy,3)+CofA(LblEmDarcy,4)*xhat) .* (CofB(LblEmDarcy,2)+CofB(LblEmDarcy,4)*yhat);
  GlbRHS(1:NumEmsDarcy) = GlbRHS(1:NumEmsDarcy) ...
    + EqnBC.fxns_Darcy(qp) .* jac * GAUSSQUAD.RECT(k,3);
end
% For Darcy: Additional/gravity term 
% Computing local L2-projection coeff. fxng_Darcy into AC0 
GM = Hdiv_QuadriAC0_NmlzPiolaBas_GramMat(QuadriMesh, GAUSSQUAD);
EltRHS = zeros(NumEmsDarcy, 5);
for i=1:5
  for j=1:4
    for k=1:4
      EltRHS(:,i) = EltRHS(:,i) + CDWGB(LblEmDarcy,i,j) .* GM(LblEmDarcy,j,k) .* CofQhg(:,k);
    end
  end
end
% Now assembly for the additional/gravity term to GlbRHS 
pos = zeros(NumEmsDarcy, 5);
pos(:,1) = LblEmDarcy;
for i=2:5
  pos(:,i) = NumEmsDarcy + QuadriMesh.posDarcyEg(QuadriMesh.elem2edge(LblEmDarcy,i-1));
end
for i=1:5
  GlbRHS(pos(:,i)) = GlbRHS(pos(:,i)) - EltRHS(:,i);  % JL20200726: "-"? 
end
% For Stokes: Source (body force) 
EltRHS = zeros(NumEmsStokes,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);
  qp = [CofA(LblEmStokes,1) + CofA(LblEmStokes,2)*xhat + CofA(LblEmStokes,3)*yhat + CofA(LblEmStokes,4)*xhat*yhat, ...
        CofB(LblEmStokes,1) + CofB(LblEmStokes,2)*xhat + CofB(LblEmStokes,3)*yhat + CofB(LblEmStokes,4)*xhat*yhat];
  jac = (CofB(LblEmStokes,3)+CofB(LblEmStokes,4)*xhat) .* (CofA(LblEmStokes,2)+CofA(LblEmStokes,4)*yhat) ...
      - (CofA(LblEmStokes,3)+CofA(LblEmStokes,4)*xhat) .* (CofB(LblEmStokes,2)+CofB(LblEmStokes,4)*yhat);
  fval = EqnBC.fxnf_Stokes(qp);
  BFVal = CG_QuadriBR1_BasFxnVal(QuadriMesh,xhat,yhat);
  for j=1:12
    EltRHS(:,j) = EltRHS(:,j) + GAUSSQUAD.RECT(k,3) * jac .* ...
      (BFVal(LblEmStokes,j,1).*fval(:,1) + BFVal(LblEmStokes,j,2).*fval(:,2));
  end 
end
% Modification for the 4 edge-based basis functions 
for j=1:4
  EltRHS(:,8+j) = EltRHS(:,8+j) .* QuadriMesh.EmEgSign(LblEmStokes,j);
end
% Now assembly of contributions from elementwise 12 basis functions 
for ie=1:NumEmsStokes
  for j=1:4  % Treating vertices and edges in one loop 
    k = QuadriMesh.elem(LblEmStokes(ie),j);
    loc1 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(k) - 1;
    loc2 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(k);
    GlbRHS(loc1) = GlbRHS(loc1) + EltRHS(ie, 2*j-1);  % For vertex
    GlbRHS(loc2) = GlbRHS(loc2) + EltRHS(ie, 2*j  );  % For vertex
    k = QuadriMesh.elem2edge(LblEmStokes(ie),j);
    loc = NumDOFsDarcy + 2*NumNdsStokes + QuadriMesh.posStokesEg(k);
    GlbRHS(loc) = GlbRHS(loc) + EltRHS(ie, 8+j);  % For edge
  end
end

%% Assembly/Modification of GlbRHS: Incorporating natural boundary conditions 
% For Stokes: Neumann (traction) = natural 
if NumNeumannEgsStokes>0 
  % For end-nodes 
  NeumannBC = zeros(NumNeumannEgsStokes,2,2);
  j1 = QuadriMesh.edge(StokesNeumannEdge,1);
  j2 = QuadriMesh.edge(StokesNeumannEdge,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 = EqnBC.fxntN_Stokes(qp);
    NeumannBC(:,1,1) = NeumannBC(:,1,1) ... 
      + GAUSSQUAD.LINE(k,1) * tNval(:,1) .* LenEgNeumann * GAUSSQUAD.LINE(k,3);
    NeumannBC(:,1,2) = NeumannBC(:,1,2) ...
      + GAUSSQUAD.LINE(k,1) * tNval(:,2) .* LenEgNeumann * GAUSSQUAD.LINE(k,3);
    NeumannBC(:,2,1) = NeumannBC(:,2,1) ...
      + GAUSSQUAD.LINE(k,2) * tNval(:,1) .* LenEgNeumann * GAUSSQUAD.LINE(k,3);
    NeumannBC(:,2,2) = NeumannBC(:,2,2) ... 
      + GAUSSQUAD.LINE(k,2) * tNval(:,2) .* LenEgNeumann * GAUSSQUAD.LINE(k,3);
  end
  loc11 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(j1) - 1;
  loc12 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(j1);
  loc21 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(j2) - 1;
  loc22 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(j2);  
  GlbRHS(loc11) = GlbRHS(loc11) + NeumannBC(:,1,1);  % But "+" here 
  GlbRHS(loc12) = GlbRHS(loc12) + NeumannBC(:,1,2);  % But "+" here 
  GlbRHS(loc21) = GlbRHS(loc21) + NeumannBC(:,2,1);  % But "+" here 
  GlbRHS(loc22) = GlbRHS(loc22) + NeumannBC(:,2,2);  % But "+" here 
  % For edges 
  NeumannBC = zeros(NumNeumannEgsStokes,1);
  j1 = QuadriMesh.edge(StokesNeumannEdge,1);
  j2 = QuadriMesh.edge(StokesNeumannEdge,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(StokesNeumannEdge,:);
  NumQuadPts = size(GAUSSQUAD.LINE,1);
  for k=1:NumQuadPts
    qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
    tNval = EqnBC.fxntN_Stokes(qp);
    psi = GAUSSQUAD.LINE(k,1) * GAUSSQUAD.LINE(k,2);
    NeumannBC = NeumannBC ... 
      + sum(tNval.*nml,2) * psi .* LenEgNeumann * GAUSSQUAD.LINE(k,3);
  end
  loc = NumDOFsDarcy + 2*NumNdsStokes + QuadriMesh.posStokesEg(StokesNeumannEdge);
  GlbRHS(loc) = GlbRHS(loc) + NeumannBC;
end
% For Darcy: Neumann = natural 
if NumNeumannEgsDarcy>0 
  NeumannBC = zeros(NumNeumannEgsDarcy,1);
  x1 = QuadriMesh.node(QuadriMesh.edge(DarcyNeumannEdge,1),1);
  y1 = QuadriMesh.node(QuadriMesh.edge(DarcyNeumannEdge,1),2);
  x2 = QuadriMesh.node(QuadriMesh.edge(DarcyNeumannEdge,2),1);
  y2 = QuadriMesh.node(QuadriMesh.edge(DarcyNeumannEdge,2),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];
    uNval = EqnBC.fxnuN_Darcy(qp);
    NeumannBC = NeumannBC + uNval .* LenEgNeumann * GAUSSQUAD.LINE(k,3);
  end
  loc = NumEmsDarcy + QuadriMesh.posDarcyEg(DarcyNeumannEdge);  
  GlbRHS(loc) = GlbRHS(loc) - NeumannBC;  % Note: "-" negative sign 
end

%% Enforcing essential boundary conditions: Darcy & Stokes  
sln = zeros(DOFs,1);
CnstrDOFs = zeros(DOFs,1);
% For Stokes Dirichlet bndry. cond.: Assuming NumDirichletEgsStokes>0 
% For Stokes: Dirichlet nodes: Interpolation 
uDval = EqnBC.fxnuD_Stokes(QuadriMesh.node(StokesDirichletNode,:));
sln(NumDOFsDarcy + 2*QuadriMesh.posStokesNd(StokesDirichletNode)-1) = uDval(:,1);
sln(NumDOFsDarcy + 2*QuadriMesh.posStokesNd(StokesDirichletNode)  ) = uDval(:,2);
CnstrDOFs(NumDOFsDarcy + 2*QuadriMesh.posStokesNd(StokesDirichletNode)-1) = 1;
CnstrDOFs(NumDOFsDarcy + 2*QuadriMesh.posStokesNd(StokesDirichletNode)  ) = 1;
% For Stokes: Dirichlet edges: Based on projection of interpolation residual 
intgrl = zeros(NumDirichletEgsStokes,1);
nml = QuadriMesh.EgNml(StokesDirichletEdge,:);
j1 = QuadriMesh.edge(StokesDirichletEdge,1);
j2 = QuadriMesh.edge(StokesDirichletEdge,2);
loc1 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(j1);
loc2 = NumDOFsDarcy + 2*QuadriMesh.posStokesNd(j2);
tmp1 = sln([loc1-1, loc1]);
tmp2 = sln([loc2-1, loc2]);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts
  s = GAUSSQUAD.LINE(k,2);
  qp = GAUSSQUAD.LINE(k,1) * QuadriMesh.node(j1,:) ... 
     + GAUSSQUAD.LINE(k,2) * QuadriMesh.node(j2,:);
  val = EqnBC.fxnuD_Stokes(qp);
  intplt = (1-s) * tmp1 + s * tmp2;
  res = val - intplt;
  intgrl = intgrl + dot(res,nml,2) * GAUSSQUAD.LINE(k,3);
end
loc = NumDOFsDarcy + 2*NumNdsStokes + QuadriMesh.posStokesEg(StokesDirichletEdge);
sln(loc) = intgrl * 6;  % JL20200726: Double-check "6" 
CnstrDOFs(loc) = 1;
% For Darcy: Dirichlet edges: Assuming NumDirichletEgsDarcy>0 
pDavg = zeros(NumDirichletEgsDarcy,1);
x1 = QuadriMesh.node(QuadriMesh.edge(DarcyDirichletEdge,1),1);
y1 = QuadriMesh.node(QuadriMesh.edge(DarcyDirichletEdge,1),2);
x2 = QuadriMesh.node(QuadriMesh.edge(DarcyDirichletEdge,2),1);
y2 = QuadriMesh.node(QuadriMesh.edge(DarcyDirichletEdge,2),2);
NumQuadPts = size(GAUSSQUAD.LINE,1);
for k=1:NumQuadPts
  qp = GAUSSQUAD.LINE(k,1)*[x1,y1] + GAUSSQUAD.LINE(k,2)*[x2,y2];
  pDavg = pDavg + GAUSSQUAD.LINE(k,3) * EqnBC.fxnpD_Darcy(qp);
end
loc = NumEmsDarcy + QuadriMesh.posDarcyEg(DarcyDirichletEdge);
sln(loc) = pDavg; 
CnstrDOFs(loc) = 1;

% %% For uniqueness of Stokes pressure 
% loc = NumDOFsDarcy + 2*NumNdsStokes + NumEgsStokes + 1;
% sln(loc) = 1;
% CnstrDOFs(loc) = 1;

%% Reduction 
GlbRHS = GlbRHS - GlbMat*sln;

%% Solving the reduced global linear system directly 
FreeDOFs = find(~CnstrDOFs);
sln(FreeDOFs) = GlbMat(FreeDOFs,FreeDOFs) \ GlbRHS(FreeDOFs);

return;