function QuadriMesh = QuadriMesh_Enrich_StokesDarcy(QuadriMesh)
%% Generating a list of geometry-to-DOF information for QuadriMesh
% Assumes QuadriMesh has two lists, DarcyEms and StokesEms, which list the
% Darcy elements and the Stokes elements.  This does not assume number of
% DoFs per edge or number of DoFs per vertex.  Those are handled in the
% AsmSlv routine.  This eases the burden of indexing the global system for
% the AsmSlv, VelPresFlux, and ShowResults routines.
% Graham Harper, James Liu, ColoState; 2014/07--2019/11 

% if (~isfield(QuadriMesh,'StokesEms') && ~isfield(QuadriMesh,'DarcyEms'))  % If no labels are given
%   error(['Error! Please create a vector QuadriMesh.StokesEms or ' ...
%     'QuadriMesh.DarcyEms that defines which cells are Stokes cells ' ...
%     'or which cells are Darcy cells. Aborting...']);
% end
% 
% if(~isfield(QuadriMesh, 'StokesEms'))  % If only Darcy elements are given
%   StokesEms = zeros(QuadriMesh.NumEms - length(QuadriMesh.DarcyEms),1);
%   % TODO: Fill in StokesEms as elements that are not Darcy
% end
% 
% if(~isfield(QuadriMesh, 'DarcyEms'))
%   DarcyEms = zeros(QuadriMesh.NumEms - length(QuadriMesh.StokesEms),1);
%   % TODO: Fill in DarcyEms as elements that are not Stokes
% end

%% For this special case 
nx = QuadriMesh.nx;
ny = QuadriMesh.ny;

%% For this special case 
QuadriMesh.StokesEms = zeros(1, nx*(ny/2));
QuadriMesh.DarcyEms  = zeros(1, nx*(ny/2));
for i=1:nx
  loc = (i-1)*(ny/2) + (1:(ny/2));
  QuadriMesh.StokesEms(loc) = (i-1)*ny + (1:(ny/2)) + (ny/2);
  QuadriMesh.DarcyEms(loc)  = (i-1)*ny + (1:(ny/2));
end

%% For this special case 
QuadriMesh.BndryEdge(nx*ny+(1:(ny/2))) = 2;     % Right-lower for Darcy 
QuadriMesh.BndryEdge(nx*ny+((ny/2)+1:ny)) = 3;  % Right-upper for Stokes 
QuadriMesh.BndryEdge((ny/2)+1:ny) = 5;          % Left-upper  for Stokes 
QuadriMesh.BndryEdge(1:(ny/2)) = 6;             % Left-lower  for Darcy 
% Right-lower boundary normal 
QuadriMesh.EgNml(nx*ny+(1:(ny/2)), 1) = 1;     
QuadriMesh.EgNml(nx*ny+(1:(ny/2)), 2) = 0;
% Right-upper boundary normal 
QuadriMesh.EgNml(nx*ny+((ny/2)+1:ny), 1) = 1;
QuadriMesh.EgNml(nx*ny+((ny/2)+1:ny), 2) = 0;
% Left-upper boundary normal 
QuadriMesh.EgNml((ny/2)+1:ny, 1) = -1;
QuadriMesh.EgNml((ny/2)+1:ny, 2) =  0;
% Left-lower boundary normal 
QuadriMesh.EgNml(1:(ny/2), 1) = -1;
QuadriMesh.EgNml(1:(ny/2), 2) =  0;

%% Computing posDarcyNd, posDarcyEg, posStokesNd, posStokesEg 
% posStokesEg maps a global mesh edge index to a DoF index for all edges
% adjacent to a Stokes cell.  If the edge is not adjacent to a Stokes cell,
% it maps to 0.
% posStokesNd does the same except for nodes.
% posDarcyEg and posDarcyNd do the same except for the Darcy cells.
posStokesNd = zeros(QuadriMesh.NumNds,1);
posStokesEg = zeros(QuadriMesh.NumEgs,1);
posDarcyNd = zeros(QuadriMesh.NumNds,1);
posDarcyEg = zeros(QuadriMesh.NumEgs,1);

%% JL20200727: SPECIAL TREATMENT FOR NOW 
% Assuming nx, ny each is an even number 
% For nodes 
kDarcy = 0;  kStokes = 0;  
for i=0:nx 
  for j=0:(ny/2)
    k = i*(ny+1) + j+1;
    kDarcy = kDarcy + 1;
    posDarcyNd(k) = kDarcy;
    % posStokesNd(k) = 0;
  end
  for j=(ny/2):ny 
    k = i*(ny+1) + j+1;
    kStokes = kStokes + 1;
    % posDarcyNd(k) = 0;
    posStokesNd(k) = kStokes;
  end
end
% For vertical edges 
kDarcy = 0;  kStokes = 0;  
for i=0:nx
  for j=1:(ny/2)
    k = i*ny + j;
    kDarcy = kDarcy + 1;
    posDarcyEg(k) = kDarcy;
  end
  for j=(ny/2)+1:ny
    k = i*ny + j;
    kStokes = kStokes + 1;
    posStokesEg(k) = kStokes;
  end
end
% For horizontal edges 
for i=1:nx
  for j=0:(ny/2)-1
    k = (nx+1)*ny + (i-1)*(ny+1) + j+1;
    kDarcy = kDarcy + 1;
    posDarcyEg(k) = kDarcy;    
  end
  j = ny/2;
  k = (nx+1)*ny + (i-1)*(ny+1) + j+1;
  kDarcy = kDarcy + 1;
  posDarcyEg(k) = kDarcy;    
  kStokes = kStokes + 1;
  posStokesEg(k) = kStokes;
  for j=(ny/2)+1:ny
    k = (nx+1)*ny + (i-1)*(ny+1) + j+1;
    kStokes = kStokes + 1;
    posStokesEg(k) = kStokes;
  end
end

%% 
QuadriMesh.posDarcyNd = posDarcyNd;
QuadriMesh.posDarcyEg = posDarcyEg;
QuadriMesh.posStokesNd = posStokesNd;
QuadriMesh.posStokesEg = posStokesEg;
QuadriMesh.InterfaceNds = find(posDarcyNd & posStokesNd);
QuadriMesh.InterfaceEgs = find(posDarcyEg & posStokesEg);

return;