function HybridMesh = RectDom_HybridMesh_GenLogiRectYlowPtb(...
  xa,xb,nx,yc,yd,ny,my,delta)
%% For a given rectangular domain, 
% generating a hybrid mesh (tri.+quadri.) that is logically rectangular 
% Y-dir. lower part for quadrilateralss; with prescribed perturbation pattern 
% 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 
%   my: for the lower my partitons, elements are quadrilaterals 
%       for the upper, elements are triangles 
%   delta: magnitude for perterbation for interior nodes 
% Output: 
%   HybridMesh: 
% JL20170111: This function is mainly for simplicity not efficiency 
% James Liu, Zhuoran Wang, ColoState; 2007/01--2017/01

%% 
HybridMesh.NumNds  = (nx+1)*(ny+1);  % Number of nodes 
HybridMesh.NumEmsQ = my * ny;
HybridMesh.NumEmsT = 2 * nx * (ny-my);
HybridMesh.NumEms  = HybridMesh.NumEmsQ + HybridMesh.NumEmsT;

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

%% JL20170111: TO BE REVISED FOR EFFECIENCY 
%% Prescribed perturbation pattern in WheeXueYotov,_NumerMath_2012, p.193
for j=2:ny
  for i=2:nx
    X(i,j) = X(i,j) + 0.06 * sin((i-1)/nx*pi*2) * sin((j-1)/ny*pi*2);
    Y(i,j) = Y(i,j) - 0.05 * sin((i-1)/nx*pi*2) * sin((j-1)/ny*pi*2);
  end
end

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

%% 
HybridMesh.labelemQ = (1:nx*my)';
HybridMesh.labelemT = nx*my + (1:2*nx*(ny-my))';

%% element-vs-nodes 
% Lower part: quadrilaterals, counterclockwise orientation 
k = 1;
for j=1:my
  for i=1:nx 
    HybridMesh.elemQ(k,1) = (i-1) * (ny+1) + j;
    HybridMesh.elemQ(k,2) = HybridMesh.elemQ(k,1) + (ny+1);
    HybridMesh.elemQ(k,3) = HybridMesh.elemQ(k,1) + (ny+2);
    HybridMesh.elemQ(k,4) = HybridMesh.elemQ(k,1) + 1;
    k = k + 1;
  end  
end
% Upper part: triangles, counterclockwise orientation 
k = 1;
for j=(my+1):ny
  for i=1:nx 
    % Left triangles: Starting at the right-angle vertex 
    HybridMesh.elemT(k,1) = (i-1) * (ny+1) + (j+1);
    HybridMesh.elemT(k,2) = HybridMesh.elemT(k,1) - 1; 
    HybridMesh.elemT(k,3) = HybridMesh.elemT(k,1) + (ny+1); 
    k = k + 1;
    % Right triangles: Starting at the right-angle vertex 
    HybridMesh.elemT(k,1) = i * (ny+1) + j;
    HybridMesh.elemT(k,2) = HybridMesh.elemT(k,1) + 1;
    HybridMesh.elemT(k,3) = HybridMesh.elemT(k,1) - (ny+1);
    k = k + 1;
  end
end

%% Finishing primary mesh info
HybridMesh.flag = 1;

%% Edges 
NumHrznEgs = nx*(ny+1);
NumVertEgs = (nx+1)*ny;
NumDiagEgs = nx*(ny-my);
NumEgs = NumHrznEgs + NumVertEgs + NumDiagEgs;
HybridMesh.NumEgs = NumEgs;
HybridMesh.edge = zeros(NumEgs,2);
% For the horizontal edges 
k = 1;
for j=1:ny+1
  for i=1:nx
    HybridMesh.edge(k,1) = (i-1)*(ny+1) + j;
    HybridMesh.edge(k,2) = HybridMesh.edge(k,1) + (ny+1); 
    k = k + 1;
  end
end
% For the vertical edges 
k = NumHrznEgs + 1;
for j=1:ny
  for i=1:(nx+1)
    HybridMesh.edge(k,1) = (i-1)*(ny+1) + j;
    HybridMesh.edge(k,2) = HybridMesh.edge(k,1) + 1;
    k = k + 1;
  end
end
% For the diagonal edges 
k = NumHrznEgs + NumVertEgs + 1;
for j=(my+1):ny
  for i=1:nx
    HybridMesh.edge(k,1) = (i-1)*(ny+1) + j;
    HybridMesh.edge(k,2) = HybridMesh.edge(k,1) + (ny+2);
    k = k + 1;
  end
end

%% Generating elemQedge, elemTedge 
% For the quadrilateral elements 
HybridMesh.elemQedge = zeros(HybridMesh.NumEmsQ,4);
k = 1;
for j=1:my
  for i=1:nx
    HybridMesh.elemQedge(k,1) = (j-1)*nx + i;
    HybridMesh.elemQedge(k,2) = NumHrznEgs + (j-1)*(nx+1) + (i+1);
    HybridMesh.elemQedge(k,3) = HybridMesh.elemQedge(k,1) + nx;
    HybridMesh.elemQedge(k,4) = HybridMesh.elemQedge(k,2) - 1;
    k = k + 1;
  end
end
% For the triangular elements 
HybridMesh.elemTedge = zeros(HybridMesh.NumEmsT,3);
k = 1;
for j=(my+1):ny
  for i=1:nx
    % For the left triangles 
    HybridMesh.elemTedge(k,1) = NumHrznEgs + NumVertEgs + (j-my-1)*nx + i;
    HybridMesh.elemTedge(k,2) = j*nx + i;
    HybridMesh.elemTedge(k,3) = NumHrznEgs + (j-1)*(nx+1) + i;
    k = k + 1;
    % For the right triangles 
    HybridMesh.elemTedge(k,1) = NumHrznEgs + NumVertEgs + (j-my-1)*nx + i;
    HybridMesh.elemTedge(k,2) = (j-1)*nx + i;
    HybridMesh.elemTedge(k,3) = NumHrznEgs + (j-1)*(nx+1) + (i+1);
    k = k + 1;
  end
end

%% JL20160716: TO BE REVISED FOR EFFECIENCY 
%% Generating edge2elem based on elemQedge, elemTedge, labelemQ, labelemT  
HybridMesh.edge2elem = zeros(HybridMesh.NumEgs,2);
cntEmsEg = zeros(HybridMesh.NumEgs,1);
% For quadrilateral elements 
for ie=1:HybridMesh.NumEmsQ
  for j=1:4
    k = HybridMesh.elemQedge(ie,j);
    cntEmsEg(k) = cntEmsEg(k) + 1;
    HybridMesh.edge2elem(k,cntEmsEg(k)) = HybridMesh.labelemQ(ie);
  end
end
% For triangular elements 
for ie=1:HybridMesh.NumEmsT
  for j=1:3
    k = HybridMesh.elemTedge(ie,j);
    cntEmsEg(k) = cntEmsEg(k) + 1;
    HybridMesh.edge2elem(k,cntEmsEg(k)) = HybridMesh.labelemT(ie);
  end
end

%% 
HybridMesh.EmCntr = zeros(HybridMesh.NumEms,2);
% For quadrilateral elements
k1 = HybridMesh.elemQ(:,1);  k2 = HybridMesh.elemQ(:,2);  
k3 = HybridMesh.elemQ(:,3);  k4 = HybridMesh.elemQ(:,4);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2);
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
x4 = HybridMesh.node(k4,1);  y4 = HybridMesh.node(k4,2);
xc = 0.25*(x1+x2+x3+x4);
yc = 0.25*(y1+y2+y3+y4);
HybridMesh.EmCntr(HybridMesh.labelemQ,:) = [xc,yc];
% For triangular elements
k1 = HybridMesh.elemT(:,1);  
k2 = HybridMesh.elemT(:,2);  
k3 = HybridMesh.elemT(:,3);
x1 = HybridMesh.node(k1,1);  y1 = HybridMesh.node(k1,2); 
x2 = HybridMesh.node(k2,1);  y2 = HybridMesh.node(k2,2);
x3 = HybridMesh.node(k3,1);  y3 = HybridMesh.node(k3,2);
xc = (1.0/3)*(x1+x2+x3);
yc = (1.0/3)*(y1+y2+y3);
HybridMesh.EmCntr(HybridMesh.labelemT,:) = [xc,yc];

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