function QuadriMesh = CircSecDom_QuadriMesh_GenUnfm(...
  ra,rb,nr,thetac,thetad,ntheta,status)
%% Generating a quadrilateral mesh for a (non-close) circular sector domain
% Polar coodinate description: 
%   0<ra<=r<rb, 0<=thetac<=theta<=thetad<2*pi
% Logically rectangular partitions: 
%       nr: number of partitions in r-direction
%   ntheta: number of partitions in theta direction 
% status=1: basic mesh info 
% status=2: more mesh info 
% status=3: boundary edge marks
% James Liu, ColoState; 2012/07--2016/10

%% For status=1
QuadriMesh.NumNds = (nr+1) * (ntheta+1);
QuadriMesh.NumEms = nr * ntheta;

%% Generating mesh info: nodes 
r = linspace(ra,rb,nr+1); 
theta = linspace(thetac,thetad,ntheta+1);
[R,S] = meshgrid(r,theta);
QuadriMesh.node = [R(:).*cos(S(:)), R(:).*sin(S(:))];

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

%% Generating mesh info: element-vs-nodes (counterclockwise orientation) 
QuadriMesh.elem = zeros(QuadriMesh.NumEms,4);
lbl = LblNd(1:ntheta,1:nr);  % Right to left, inner to outer 
QuadriMesh.elem(:,1) = lbl(:);                             % right-inner 
QuadriMesh.elem(:,2) = QuadriMesh.elem(:,1) + (ntheta+1);  % right-outer 
QuadriMesh.elem(:,3) = QuadriMesh.elem(:,1) + (ntheta+2);  % left-outer 
QuadriMesh.elem(:,4) = QuadriMesh.elem(:,1) + 1;           % left-inner 

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

%% Further work for status=2
% For status=2 (Secondary mesh info)
NumEgsVert = (nr+1) * ntheta;
NumEgsHori = nr * (ntheta+1);
QuadriMesh.NumEgs = NumEgsVert + NumEgsHori;

%% Generating secondary mesh info: edges 
QuadriMesh.edge = zeros(QuadriMesh.NumEgs,2);
% Angular edges 
lbl = LblNd(1:ntheta, 1:nr+1);
QuadriMesh.edge(1:NumEgsVert,1) = lbl(:);
QuadriMesh.edge(1:NumEgsVert,2) = QuadriMesh.edge(1:NumEgsVert,1) + 1;
% Radial edges 
lbl = LblNd(1:ntheta+1, 1:nr);
QuadriMesh.edge(NumEgsVert+1:end,1) = lbl(:); 
QuadriMesh.edge(NumEgsVert+1:end,2) = ...
  QuadriMesh.edge(NumEgsVert+1:end,1) + (ntheta+1); 

%% Generating secondary mesh info: element-vs-edges 
QuadriMesh.elem2edge = zeros(QuadriMesh.NumEms,4);
lbl = LblEgVert(1:ntheta, 1:nr);
QuadriMesh.elem2edge(:,4) = lbl(:);
QuadriMesh.elem2edge(:,2) = QuadriMesh.elem2edge(:,4) + ntheta;
lbl = LblEgHori(1:ntheta, 1:nr);
QuadriMesh.elem2edge(:,1) = NumEgsVert + lbl(:);
QuadriMesh.elem2edge(:,3) = QuadriMesh.elem2edge(:,1) + 1;

%% Generating secondary mesh info: edge-vs-elements based on elem2edge 
QuadriMesh.edge2elem = zeros(QuadriMesh.NumEgs,2);
CntEmsEg = zeros(QuadriMesh.NumEgs,1);
for ie=1:QuadriMesh.NumEms
  LblEg = QuadriMesh.elem2edge(ie,1:4);
  CntEmsEg(LblEg) = CntEmsEg(LblEg) + 1;
  for k=1:4
    QuadriMesh.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(QuadriMesh.edge2elem(:,1)>QuadriMesh.edge2elem(:,2));
tmp = QuadriMesh.edge2elem(ig,1);
QuadriMesh.edge2elem(ig,1) = QuadriMesh.edge2elem(ig,2);
QuadriMesh.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(QuadriMesh.edge2elem(:,1)==0);
QuadriMesh.edge2elem(ig,1) = QuadriMesh.edge2elem(ig,2);
QuadriMesh.edge2elem(ig,2) = 0;

%% Generating secondary mesh info: area and center 
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
QuadriMesh.area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
                      + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));
QuadriMesh.EmCntr = [0.25*(x1+x2+x3+x4), 0.25*(y1+y2+y3+y4)];

%% Coefficients for the bilinear mapping 
QuadriMesh.CofA = zeros(QuadriMesh.NumEms,4);       
QuadriMesh.CofA(:,1) = x1;
QuadriMesh.CofA(:,2) = x2-x1;
QuadriMesh.CofA(:,3) = x4-x1;
QuadriMesh.CofA(:,4) = (x1+x3)-(x2+x4);
QuadriMesh.CofB = zeros(QuadriMesh.NumEms,4);
QuadriMesh.CofB(:,1) = y1;
QuadriMesh.CofB(:,2) = y2-y1;
QuadriMesh.CofB(:,3) = y4-y1;
QuadriMesh.CofB(:,4) = (y1+y3)-(y2+y4);

%% Finishing secondary mesh info
QuadriMesh.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 

%% Further work for status=2
% Setting boundary edge marks 
QuadriMesh.BndryEdge = zeros(QuadriMesh.NumEgs,1);
QuadriMesh.BndryEdge(1:ntheta) = 4;                             % inner 
QuadriMesh.BndryEdge(nr*ntheta+(1:ntheta)) = 2;                 % outer 
QuadriMesh.BndryEdge((nr+1)*ntheta+1+(ntheta+1)*(0:nr-1)) = 1;  % right/start 
QuadriMesh.BndryEdge((nr+2)*ntheta+1+(ntheta+1)*(0:nr-1)) = 3;  % left/end 

return;