function [sln] ... 
  = StokesDarcy_QuadriCGBR1WGAC0_3domH_AsmSlv( ... 
    EqnBC1, BndryCondType1, QuadriMesh1, ... 
    EqnBC2, BndryCondType2, QuadriMesh2, ... 
    EqnBC3, BndryCondType3, QuadriMesh3, GAUSSQUAD)
%% Stokes-Darcy-Stokes 3-domain horizontal coupling 
% James Liu, ColoState; Graham Harper, Sandia; 2019/07--2020/09 

%% Equation info 
%    nu = EqnBC.nu;
% alpha = EqnBC.alpha;
nu = 1;  alpha = 1;
kappa = EqnBC2.kappa;

%% Domain #1/2/3: Stokes/Darcy/Stokes on QuadriMesh1/2/3 
[GlbMat1, GlbRHS1, GlbVecEss1, FlagDOFs1] = Stokes_CG_QuadriBR1Q0_Asm( ... 
  EqnBC1, BndryCondType1, QuadriMesh1, GAUSSQUAD); 
[GlbMat2, GlbRHS2, GlbVecEss2, FlagDOFs2] = DarcyGen_WG_QuadriP0P0AC0_Asm( ... 
  EqnBC2, BndryCondType2, QuadriMesh2, GAUSSQUAD);
[GlbMat3, GlbRHS3, GlbVecEss3, FlagDOFs3] = Stokes_CG_QuadriBR1Q0_Asm( ... 
  EqnBC3, BndryCondType3, QuadriMesh3, GAUSSQUAD);

%% Setup for DOFs 
% DOFs1Stokes = size(GlbMat1,1);
% DOFs2Darcy  = size(GlbMat2,1);
% DOFs3Stokes = size(GlbMat3,1);
DOFs1Stokes = 2*QuadriMesh1.NumNds + QuadriMesh1.NumEgs + QuadriMesh1.NumEms;
DOFs3Stokes = 2*QuadriMesh3.NumNds + QuadriMesh3.NumEgs + QuadriMesh3.NumEms;
DOFs2Darcy  = QuadriMesh2.NumEms + QuadriMesh2.NumEgs;
DOFs = DOFs1Stokes + DOFs2Darcy + DOFs3Stokes;

%% Revising GlbMat1, GlbMat3 for the BJS conditions 
% For interface #1: Using QuadriMesh1 info: nml=[1;0], tan=[0;1] 
nx = QuadriMesh1.nx;  ny = QuadriMesh1.ny;
tangent =[0,1];
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig=ny*nx + (1:ny)  % Edge labels in QuadriMesh1 
  % 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 
  LocalMat = zeros(4,4);
  for k=1:NumQuadPts
    r = GAUSSQUAD.LINE(k,2);
    StokesVelShapeVal = [tangent*((1-r)*[1;0]); ... % bas.fxn.#1
                         tangent*((1-r)*[0;1]); ... % bas.fxn.#2
                         tangent*(   r *[1;0]); ... % bas.fxn.#3
                         tangent*(   r *[0;1])];    % bas.fxn.#4
    % Grahma's Notes: No need for bas.fxn.#5 since tangent*normal = 0 
    % JL20200731: sqrt(1) TO BE REVISED !! 
    LocalMat = LocalMat + StokesVelShapeVal * StokesVelShapeVal' * ...
      GAUSSQUAD.LINE(k,3) * QuadriMesh1.EgLen(ig) * alpha*nu/sqrt(kappa);
  end
  II = zeros(4,1);
  for j=1:2  % Two vertices for each edge, two components 
    II(2*j-1) = 2*QuadriMesh1.edge(ig,j) - 1;
    II(2*j  ) = 2*QuadriMesh1.edge(ig,j);
  end
  JJ = II;  
  for i=1:4
    for j=1:4
      GlbMat1 = GlbMat1 + sparse(II(i),JJ(j), LocalMat(i,j), DOFs1Stokes,DOFs1Stokes);
    end
  end
end
% For interface #2: Using QuadriMesh3 info: nml=[-1;0], tan=[0;-1] 
nx = QuadriMesh3.nx;  ny = QuadriMesh3.ny;
tangent =[0,-1];
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig = (1:ny)  % Edge labels in QuadriMesh3 
  % 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 
  LocalMat = zeros(4,4);
  for k=1:NumQuadPts
    r = GAUSSQUAD.LINE(k,2);
    StokesVelShapeVal = [tangent*((1-r)*[1;0]); ... % bas.fxn.#1
                         tangent*((1-r)*[0;1]); ... % bas.fxn.#2
                         tangent*(   r *[1;0]); ... % bas.fxn.#3
                         tangent*(   r *[0;1])];    % bas.fxn.#4
    % Grahma's Notes: No need for bas.fxn.#5 since tangent*normal = 0 
    % JL20200731: sqrt(1) TO BE REVISED !! 
    LocalMat = LocalMat + StokesVelShapeVal * StokesVelShapeVal' * ...
      GAUSSQUAD.LINE(k,3) * QuadriMesh3.EgLen(ig) * alpha*nu/sqrt(kappa);
  end
  II = zeros(4,1);
  for j=1:2  % Two vertices for each edge, two components 
    II(2*j-1) = 2*QuadriMesh3.edge(ig,j) - 1;
    II(2*j  ) = 2*QuadriMesh3.edge(ig,j);
  end
  JJ = II;  
  for i=1:4
    for j=1:4
      GlbMat3 = GlbMat3 + sparse(II(i),JJ(j), LocalMat(i,j), DOFs3Stokes,DOFs3Stokes);
    end
  end
end

%% Assembling GlbMat1/2/3 into GlbMat 
% GlbMat = sparse(DOFs,DOFs);
GlbMat = [ ... 
  GlbMat1, sparse(DOFs1Stokes,DOFs2Darcy), sparse(DOFs1Stokes,DOFs3Stokes); ...
  sparse(DOFs2Darcy,DOFs1Stokes), GlbMat2, sparse(DOFs2Darcy, DOFs3Stokes); ...
  sparse(DOFs3Stokes,DOFs1Stokes), sparse(DOFs3Stokes,DOFs2Darcy), GlbMat3];

%% Computing and Assembing the INTERFACE conditions 
% For interface #1: Using QuadriMesh1 info: nml=[1;0], tan=[0;1] 
nx = QuadriMesh1.nx;  ny = QuadriMesh1.ny;
nml =[1,0];
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig = ny*nx + (1:ny)  % Edge labels in QuadriMesh1 
  % x1 = QuadriMesh1.node(QuadriMesh1.edge(ig,1),1);
  % y1 = QuadriMesh1.node(QuadriMesh1.edge(ig,1),2);
  % x2 = QuadriMesh1.node(QuadriMesh1.edge(ig,2),1);
  % y2 = QuadriMesh1.node(QuadriMesh1.edge(ig,2),2);
  LocalMat = zeros(5,1);
  for k=1:NumQuadPts
    r = GAUSSQUAD.LINE(k,2);
    StokesVelShapeVal = [nml * ((1-r)*[1;0]); ...  % Bas.fxn.#1
                         nml * ((1-r)*[0;1]); ...  % Bas.fxn.#2
                         nml * (    r*[1;0]); ...  % Bas.fxn.#3
                         nml * (    r*[0;1]); ...  % Bas.fxn.#4
                                    (1-r)*r ];     % Bas.fxn.#5 
    LocalMat = LocalMat + StokesVelShapeVal * GAUSSQUAD.LINE(k,3) * QuadriMesh1.EgLen(ig);
  end 
  II = zeros(5,1);
  for j=1:2  % For each edge, its two vertices/nodes and 2 components 
    II(2*j-1) = 2*QuadriMesh1.edge(ig,j) - 1;
    II(2*j  ) = 2*QuadriMesh1.edge(ig,j);
  end
  II(5) = 2*QuadriMesh1.NumNds + ig;
  JJ = DOFs1Stokes + QuadriMesh2.NumEms + ig - ny*nx;
  for i=1:5  % All 5 basis functions are used 
    GlbMat = GlbMat + sparse(II(i),JJ,  LocalMat(i), DOFs,DOFs);
    GlbMat = GlbMat + sparse(JJ,II(i), -LocalMat(i), DOFs,DOFs);
  end 
end 
% For interface #2: Using QuadriMesh3 info: nml=[-1;0], tan=[0;1] 
nx = QuadriMesh3.nx;  ny = QuadriMesh3.ny;
nml =[-1,0];
NumQuadPts = size(GAUSSQUAD.LINE,1);
for ig=(1:ny)  % Edge labels in QuadriMesh3 
  LocalMat = zeros(5,1);
  for k=1:NumQuadPts
    r = GAUSSQUAD.LINE(k,2);
    StokesVelShapeVal = [nml * ((1-r)*[1;0]); ...  % Bas.fxn.#1
                         nml * ((1-r)*[0;1]); ...  % Bas.fxn.#2
                         nml * (    r*[1;0]); ...  % Bas.fxn.#3
                         nml * (    r*[0;1]); ...  % Bas.fxn.#4
                                   (1-r)*r ];      % Bas.fxn.#5 
    LocalMat = LocalMat + StokesVelShapeVal * GAUSSQUAD.LINE(k,3) * QuadriMesh3.EgLen(ig);
  end 
  II = zeros(5,1);
  for j=1:2  % For each edge, its two vertices/nodes and 2 components 
    II(2*j-1) = DOFs1Stokes + DOFs2Darcy + 2*QuadriMesh3.edge(ig,j) - 1;
    II(2*j  ) = DOFs1Stokes + DOFs2Darcy + 2*QuadriMesh3.edge(ig,j);
  end
  II(5) = DOFs1Stokes + DOFs2Darcy + 2*QuadriMesh3.NumNds + ig;
  JJ = DOFs1Stokes + QuadriMesh2.NumEms + ny*nx + ig ;
  for i=1:5  % All 5 basis functions are used 
    GlbMat = GlbMat + sparse(II(i),JJ,  LocalMat(i), DOFs,DOFs);
    GlbMat = GlbMat + sparse(JJ,II(i), -LocalMat(i), DOFs,DOFs);
  end 
end 

%% Assembling GlbRHS1/2/3 into GlbRHS 
GlbRHS = [GlbRHS1; GlbRHS2; GlbRHS3];

%% Enforcing the essential conditions 
sln = zeros(DOFs,1);
sln((1:DOFs1Stokes)') = GlbVecEss1; 
sln(DOFs1Stokes+(1:DOFs2Darcy)') = GlbVecEss2; 
sln(DOFs1Stokes+DOFs2Darcy+(1:DOFs3Stokes)') = GlbVecEss3;

%% Reduction 
GlbRHS = GlbRHS - GlbMat*sln;

%% Solving the reduced global linear system directly 
FreeDOFs = find(~[FlagDOFs1; FlagDOFs2; FlagDOFs3]);
sln(FreeDOFs) = GlbMat(FreeDOFs,FreeDOFs) \ GlbRHS(FreeDOFs);

return;