function RectMesh = RectDom_RectMesh_GenUnfm(xa,xb,nx,yc,yd,ny,status)
%% Generating a uniform rectangular mesh on a given rectangular domain
% (assuming it is parallel to the x-, y-axes)
% Input:  
%   xa: left  xb: right  nx: number of partitions in the x-direction 
%   yc: bottom  yd: top  ny: number of partitions in the y-direction 
% Output: 
%   RectMesh: a structure of arrays for primary mesh info 
% JL20160106: This function is mainly for simplicity not efficiency 
% Originally authored by Rachel Cali, ColoState, Spring 2007 
% James Liu, ColoState; 2007/01--2016/09

%% For status=1 (Primary mesh info) 
RectMesh.NumNds = (nx+1)*(ny+1);  % Number of nodes 
RectMesh.NumEms = nx*ny;  % Number of rectangular elements 

x = linspace(xa,xb,nx+1);
y = linspace(yc,yd,ny+1)';
[X,Y] = meshgrid(x,y);

RectMesh.xa = xa;  RectMesh.xb = xb;  RectMesh.yc = yc;  RectMesh.yd = yd;
RectMesh.nx = nx;  RectMesh.hx = (xb-xa)/nx;  RectMesh.x = x;
RectMesh.ny = ny;  RectMesh.hy = (yd-yc)/ny;  RectMesh.y = y;

% Setting the lexicongraphical order for all nodes 
RectMesh.node = [X(:),Y(:)];

% Auxiliary arrays 
LblNd = reshape((1:(nx+1)*(ny+1)),ny+1,nx+1);
LblEgVert = reshape(1:(nx+1)*ny,ny,nx+1);
LblEgHori = reshape(1:nx*(ny+1),ny+1,nx);

% Setting up element-vs-nodes (counterclockwise orientation) 
RectMesh.elem = zeros(RectMesh.NumEms,4);
lbl = LblNd(1:ny,1:nx);  % Bottom to top, left to right 
RectMesh.elem(:,1) = lbl(:);                       % lower-left 
RectMesh.elem(:,2) = RectMesh.elem(:,1) + (ny+1);  % lower-right 
RectMesh.elem(:,3) = RectMesh.elem(:,1) + (ny+2);  % upper-right 
RectMesh.elem(:,4) = RectMesh.elem(:,1) + 1;       % upper-left 

if (status==1)
  RectMesh.flag = 1;
  return;
end

%% For status=2 (Secondary mesh info)
NumEgsVert = (nx+1)*ny;
NumEgsHori = nx*(ny+1);
RectMesh.NumEgs = NumEgsVert + NumEgsHori;

% Setting up edges 
% disp('Setting up edges...');
RectMesh.edge = zeros(RectMesh.NumEgs,2);
% Vertical edges 
lbl = LblNd(1:ny,1:nx+1);
RectMesh.edge(1:NumEgsVert,1) = lbl(:);
RectMesh.edge(1:NumEgsVert,2) = RectMesh.edge(1:NumEgsVert,1) + 1;
% Horizontal edges
lbl = LblNd(1:ny+1,1:nx);
RectMesh.edge(NumEgsVert+1:end,1) = lbl(:); 
RectMesh.edge(NumEgsVert+1:end,2) = RectMesh.edge(NumEgsVert+1:end,1) + (ny+1); 

% Generating secondary mesh info on element-vs-edges 
% disp('Setting up element-vs-edges...'); 
RectMesh.elem2edge = zeros(RectMesh.NumEms,4);
lbl = LblEgVert(1:ny,1:nx);
RectMesh.elem2edge(:,4) = lbl(:);
RectMesh.elem2edge(:,2) = RectMesh.elem2edge(:,4) + ny;
lbl = LblEgHori(1:ny,1:nx);
RectMesh.elem2edge(:,1) = NumEgsVert + lbl(:);
RectMesh.elem2edge(:,3) = RectMesh.elem2edge(:,1) + 1;

% Generating secondary mesh info on edge-vs-elements based on elem2edge 
RectMesh.edge2elem = zeros(RectMesh.NumEgs,2);
CntEmsEg = zeros(RectMesh.NumEgs,1);
for ie=1:RectMesh.NumEms
  LblEg = RectMesh.elem2edge(ie,1:4);
  CntEmsEg(LblEg) = CntEmsEg(LblEg) + 1;
  for k=1:4
    RectMesh.edge2elem(LblEg(k),CntEmsEg(LblEg(k))) = ie;
  end
end

% Adjusting 
% for ig=1:RectMesh.NumEgs
%   if RectMesh.edge2elem(ig,1)>RectMesh.edge2elem(ig,2)
%     tmp = RectMesh.edge2elem(ig,1);
%     RectMesh.edge2elem(ig,1) = RectMesh.edge2elem(ig,2);
%     RectMesh.edge2elem(ig,2) = tmp;
%   end
% end
ig = find(RectMesh.edge2elem(:,1)>RectMesh.edge2elem(:,2));
tmp = RectMesh.edge2elem(ig,1);
RectMesh.edge2elem(ig,1) = RectMesh.edge2elem(ig,2);
RectMesh.edge2elem(ig,2) = tmp;
% for ig=1:RectMesh.NumEgs
%   if RectMesh.edge2elem(ig,1)==0
%     RectMesh.edge2elem(ig,1) = RectMesh.edge2elem(ig,2);
%     RectMesh.edge2elem(ig,2) = 0;
%   end
% end
ig = find(RectMesh.edge2elem(:,1)==0);
RectMesh.edge2elem(ig,1) = RectMesh.edge2elem(ig,2);
RectMesh.edge2elem(ig,2) = 0;

% area 
x1 = RectMesh.node(RectMesh.elem(:,1),1);
y1 = RectMesh.node(RectMesh.elem(:,1),2);
x2 = RectMesh.node(RectMesh.elem(:,3),1);
y2 = RectMesh.node(RectMesh.elem(:,3),2);
RectMesh.area = (x2-x1) .* (y2-y1);

% edge length 
x1 = RectMesh.node(RectMesh.edge(:,1),1);
y1 = RectMesh.node(RectMesh.edge(:,1),2);
x2 = RectMesh.node(RectMesh.edge(:,2),1);
y2 = RectMesh.node(RectMesh.edge(:,2),2);
RectMesh.LenEg = sqrt((x2-x1).^2+(y2-y1.^2));

%% Finishing secondary mesh info
RectMesh.flag = 2;

%% Afternotes: Mesh info is organized in the following format 
% status for the structure: 
%   1 - primary mesh info only; 
%   2 - secondary mesh info completed 
% NumNds, NumEms, NumEgs 
% node, elem 
% edge, elem2edge, edge2elem, area, LenEg 

return;