function QuadriMesh = TZdom_QuadriMesh_GenLogiRect(...
  xa,ya, xb,yb, xc,yc, xd,yd, nx,ny, delta, status) 
%% Generating a logically-rectangular quadrilateral mesh 
% based on prescribed perturbation pattern of a uniform rectangular mesh 
% for a trapezoidal domain with vertices (xa,ya), (xb,yb), (xc,yc), (xd,yd) 
% (in counterclockwise orientation) 
% e.g., Cook's membrane in planar linear elasticity 
% Input:  
%   xa, ya, xb, yb, xc, yc, xd, yd:  Coordinates of 4 vertices 
%       nx: Number of partitions in the x-direction 
%       ny: Number of partitions in the y-direction 
%    delta: Pertubation magnitude 
%   status: 1 or 2 
% Output: 
%   QuadriMesh: a structure of arrays for primary mesh info 
%% JL20160810: TO BE REVISED FOR EFFECIENCY
% James Liu, ColoState; 2012/07--2016/09

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

%% Partitions for [0,1]^2 (the unit square) 
x = linspace(0,1,nx+1);
y = linspace(0,1,ny+1)';
[X,Y] = meshgrid(x,y);  % Labelling from bottom to top, from left to right 

%% JL20160612: TO BE REVISED FOR EFFECIENCY 
%% Perturbation and prescribed perturbation pattern 
k = 1;
for i=1:(nx+1)
  for j=1:(ny+1)
    X(k) = X(k) + delta * sin((i-1)/nx*pi*2) * sin((j-1)/ny*pi*2);
    Y(k) = Y(k) + delta * sin((i-1)/nx*pi*2) * sin((j-1)/ny*pi*2);
    k = k + 1;
  end
end

%% Transformation from the unit square to the trapezoid 
% Coefficients 
A1 = xa;  A2 = xb - xa;  A3 = xd - xa;  A4 = (xa+xc) - (xb+xd);
B1 = ya;  B2 = yb - ya;  B3 = yd - ya;  B4 = (ya+yc) - (yb+yd);
% transformation
x = A1 + A2*X + A3*Y + A4*(X.*Y);
y = B1 + B2*X + B3*Y + B4*(X.*Y);

%% Setting the lexicongraphical order for all nodes 
QuadriMesh.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) 
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 on 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);
xa = QuadriMesh.node(k1,1);  ya = QuadriMesh.node(k1,2);
xb = QuadriMesh.node(k2,1);  yb = QuadriMesh.node(k2,2);
xc = QuadriMesh.node(k3,1);  yc = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
QuadriMesh.area = 0.5*( (xb-xa).*(yc-ya) - (xc-xa).*(yb-ya)...
                      + (xc-xa).*(y4-ya) - (x4-xa).*(yc-ya));
QuadriMesh.EmCntr = [0.25*(xa+xb+xc+x4), 0.25*(ya+yb+yc+y4)];

%% Coefficients for the bilinear mapping 
QuadriMesh.CofA = zeros(QuadriMesh.NumEms,4);       
QuadriMesh.CofA(:,1) = xa;
QuadriMesh.CofA(:,2) = xb-xa;
QuadriMesh.CofA(:,3) = x4-xa;
QuadriMesh.CofA(:,4) = (xa+xc)-(xb+x4);
QuadriMesh.CofB = zeros(QuadriMesh.NumEms,4);
QuadriMesh.CofB(:,1) = ya;
QuadriMesh.CofB(:,2) = yb-ya;
QuadriMesh.CofB(:,3) = y4-ya;
QuadriMesh.CofB(:,4) = (ya+yc)-(yb+y4);

%% Computing elementwise diameter
diag13 = sqrt((xa-xc).^2+(ya-yc).^2);
diag24 = sqrt((xb-x4).^2+(yb-y4).^2);
QuadriMesh.diam = max(diag13,diag24);

%% 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;