function QuadriMesh = RectDom_QuadriMesh_GenTZasymp(xa,xb,yc,yd,m,status)
%% Generating a quadrilateral mesh that is trapezoidal & asympotically refined 
% Ref: Arnold,Boffi,Falk, Math.Comput., 2002, p.918 Fig.2 row3 
% Input: 
%   xa: left  xb: right  yc: bottom  yd: top (of the domain, usually [0,1]^2) 
%   m>=1 a positive integer (2^m: # of partitions for x & y-directions) 
%   status: 1 or 2 or 3 for primary/secondary/tertiary mesh info 
% Output: 
%   QuadriMesh: a structure of arrays for mesh info 
% JL20160106: This function is mainly for simplicity not efficiency 
% James Liu, ColoState; 2012/07--2016/09

%% Generating auxiliary quantities 
A = zeros(3,3);
A(:,1) = [0; 1/3; 1];
A(:,2) = [0; 2/3; 1];
A(:,3) = [0; 1/3; 1];
lenthy = (yd-yc);

%% Generating y-coordinates 
if (m==1)
  n = 2^m+1;  % n=3 
  Y = A*lenthy + ones(n)*yc;
else  % Expanding for m>=2
  for k=2:m
    n = 2^k + 1;
    Y = zeros(n,n);
    for j=1:2:n
      for i=1:2:n
        Y(i,j) = A((i+1)/2,(j+1)/2);
      end
    end
    for j=1:2:n
      for i=2:2:(n-1)
        Y(i,j) = 0.5*(Y(i-1,j)+Y(i+1,j));
      end
    end
    for j=2:2:(n-1)
      Y(:,j) = 0.5*(Y(:,j-1)+Y(:,j+1));
    end
    Y = Y*lenthy + ones(n)*yc;
    A = Y;
  end
end

%% Generating x-ccordinates 
x = linspace(xa,xb,n);
X = repmat(x,n,1);

%% For status=1 (Primary mesh info) 
QuadriMesh.NumNds = n^2;  % Number of nodes 
QuadriMesh.NumEms = (n-1)^2;  % Number of rectangular elements 

%% Nodes 
QuadriMesh.node = [X(:),Y(:)];

%% Auxiliary arrays 
nx = n-1;
ny = n-1;
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) 
QuadriMesh.elem = zeros(QuadriMesh.NumEms,4);
lbl = LblNd(1:ny,1:nx);  % Bottom to top, left to right 
QuadriMesh.elem(:,1) = lbl(:);                         % lower-left 
QuadriMesh.elem(:,2) = QuadriMesh.elem(:,1) + (ny+1);  % lower-right 
QuadriMesh.elem(:,3) = QuadriMesh.elem(:,1) + (ny+2);  % upper-right 
QuadriMesh.elem(:,4) = QuadriMesh.elem(:,1) + 1;       % upper-left 

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

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

%% Setting up vertical and horizontal edges respectively 
% Vertical edges 
lbl = LblNd(1:ny,1:nx+1);
QuadriMesh.edge(1:NumEgsVert,1) = lbl(:);
QuadriMesh.edge(1:NumEgsVert,2) = QuadriMesh.edge(1:NumEgsVert,1) + 1;
% Horizonatl edges 
lbl = LblNd(1:ny+1,1:nx);
QuadriMesh.edge(NumEgsVert+1:end,1) = lbl(:); 
QuadriMesh.edge(NumEgsVert+1:end,2) = QuadriMesh.edge(NumEgsVert+1:end,1) + (ny+1); 

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

%% Generating secondary mesh info on 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 
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;
ig = find(QuadriMesh.edge2elem(:,1)==0);
QuadriMesh.edge2elem(ig,1) = QuadriMesh.edge2elem(ig,2);
QuadriMesh.edge2elem(ig,2) = 0;

%% Elementwise 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 
% node, elem 
% NumEgs, edge, elem2edge, edge2elem, 
% area, EmCntr 

return;