function msh=kmg2du(cp,msize,xc)
%KMG2DU Uniform symmetric mesh of Omega = (ax,bx)*(ay,by)
% or, if xc (the center coordinates) is given, Omega=(theta1,theta2)*(r1,r2)
% in polar coordinates. With polar coordinates, KMG2DU generates a 
% circular mesh (ring, disk). If r1=0 or theta2=2*pi, the duplicated nodes
% and boundaries are removed from the mesh.
%
% Boundary sides are numbered as follows :
% s1 line y=ay (lower side); s2 line x=bx (right side); 
% s3 line y=by (upper side); s4 line x=ax (left side)
%
%-------------------------------------------------------------------------
%      msh=kmg2du(cp,msize)
%      msh=kmg2du(cp,msize,xc)
% 
% Inputs:
%   cp      : Corner points cp=[ax,ay
%                               bx,by]
%           if xc is given, then cp must be in polar coordinates
%           cp=[theta1, r1
%               theta2, r2]
%   msize   : Number of points in each direction msize=[nx
%                                                       ny]  
%   xc      : center (cartesian) coordinates xc=[x0;y0]
%             in this case, cp is given in polar coordinates
%
% Ouputs:
%   msh structure with fields  
%      p    : Nodes coordinates, array np*2
%      t    : Triangle vertices, array nt*3
%     be    : Boundary edges 4-by-1 cell-array with be{i} edges of side si
%             3-by-1 or 2-by-1 or 1-by-1 cell-array for circular domains
%     bn    : Boundary nodes 4-by-1 cell-array with bn{i} nodes of side si
%             3-by-1 or 2-by-1 or 1-by-1 cell-array for circular domains
%-------------------------------------------------------------------------
% (c) 2013-2021, KOKO J., LIMOS/ISIMA, Jonas.Koko@isima.fr
%-------------------------------------------------------------------------
%
eps0=10^-10;
ax = cp(1,1); ay = cp(1,2);   
bx = cp(2,1); by = cp(2,2);  
nx = msize(1); ny = msize(2);  

nnx =round((nx-1)/2) +1; nny=round((ny-1)/2)+1;
np=nnx*nny; nq=(nnx-1)*(nny-1); nt=2*nq;               

% vertce coordinates
hx=(bx-ax)/(nnx-1); hy=(by-ay)/(nny-1);
[x,y]=ndgrid(ax:hx:bx,ay:hy:by);
p=[x(:),y(:)];


% quadrangles
ip=[1:nnx*nny]';
ib1=[1:nnx]'; ib2=nnx*[1:nny]';
ib3=[nnx*(nny-1)+1:nnx*nny]'; ib4=[1:nnx:nnx*nny]';

ib23=union(ib2,ib3); ib34=union(ib3,ib4);
ib14=union(ib1,ib4); ib12=union(ib1,ib2);
iq1=setdiff(ip,ib23); iq2=setdiff(ip,ib34);
iq3=setdiff(ip,ib14); iq4=setdiff(ip,ib12);

% Extract all edges
eh0=[iq1 iq2
     iq2 iq3 
     iq3 iq4
     iq4 iq1];
[eh,ih]=sort(eh0,2);

% Remove duplicates
ee=int64(eh(:,1)*(np+1)+eh(:,2));
[~,iee,je]=unique(ee); 
e=eh(iee,:);
clear eh0 ih iee ee

% Form new nodes (mid-point of edges)
ne=size(e,1);
p0=p;
p=[p0 
   (p0(e(:,1),:)+p0(e(:,2),:))/2
   (p0(iq1,:)+p0(iq2,:)+p0(iq3,:)+p0(iq4,:))/4];

% Form new edges
pm=[(np+1):(np+ne)]'; lmp=pm(je);
imq1=lmp(1:nq); imq2=lmp(nq+1:2*nq); 
imq3=lmp(2*nq+1:3*nq); imq4=lmp(3*nq+1:4*nq);
imqc=[np+ne+1:np+ne+nq]';

% Triangulation
t=[imqc iq1 imq1
   imqc imq1 iq2
   imqc iq2 imq2
   imqc imq2 iq3
   imqc iq3 imq3
   imqc imq3 iq4
   imqc iq4 imq4
   imqc imq4 iq1];

msh.p = p;
msh.t = t;

clear iq1 iq2 iq3 iq4 pm lmp
clear imq1 imq2 imq3 imq4 

% Extract boundary edges and nodes
% Extract all edges
np=size(p,1);
eh0=[t(:,1) t(:,2)
     t(:,2) t(:,3)
     t(:,3) t(:,1)];
[eh,ih]=sort(eh0,2);

% Remove duplicates
ee=int64(eh(:,1)*(np+1)+eh(:,2));
[~,iee,je]=unique(ee); 
e=eh(iee,:);

ne=size(e,1); 
he=accumarray(je,1,[ne 1]); ip=find(he==1);
e=e(ip,:); i1=e(:,1); i2=e(:,2);


ref=zeros(np,1);

 
iy1 = find(abs(p(:,2)-ay)<eps0); ref(iy1) = 1;
il1 = find(ref(i1) == 1 & ref(i2) == 1 );
ref(iy1)=0;

ix2 = find(abs(p(:,1)-bx)<eps0); ref(ix2) = 2;
il2 = find(ref(i1) == 2 & ref(i2) == 2);
ref(ix2)=0;

iy2 = find(abs(p(:,2)-by)<eps0); ref(iy2) = 3;
il3 = find(ref(i1) == 3 & ref(i2) == 3);
ref(iy2)=0;

ix1 = find(abs(p(:,1)-ax)<eps0); ref(ix1) = 4;
il4 = find(ref(i1) == 4 & ref(i2) == 4);

be = {e(il1,:)
      e(il2,:)
      e(il3,:)
      e(il4,:)};
bn = {union(e(il1,1),e(il1,2))
      union(e(il2,1),e(il2,2))
      union(e(il3,1),e(il3,2))
      union(e(il4,1),e(il4,2))};

msh.be = be;
msh.bn = bn;

if (nargin == 2)
    return
end

clear ref iy1 il1 ix2 il2 iy2 il3 ix1 il4

% convert polar -> cartesians
th=msh.p(:,1); r=msh.p(:,2);
x=xc(1)+r.*cos(th); y=xc(2)+r.*sin(th);
msh.p(:,1)=x; msh.p(:,2)=y;
np0=size(msh.p,1);

clear x y

% search for duplicated nodes & remove internal boundaries
period=0;
i=[1:np0]';
if (abs(cp(2,1)-2*pi)<10^-10)
    % i=[1:np0]';
    ii=msh.bn{2}; jj=msh.bn{4};
    i(ii) = jj; [~, I]=histc(i,unique(i));
    msh.t=[renum(I,msh.t(:,1)) renum(I,msh.t(:,2)) renum(I,msh.t(:,3))];
    msh.p(ii,:)=[];
    
    be_ext = [renum(I,msh.be{1}(:,1)) renum(I,msh.be{1}(:,2))];
    be_int = [renum(I,msh.be{3}(:,1)) renum(I,msh.be{3}(:,2))];
    msh=rmfield(msh,'be');
    msh.be{1} = be_int;
    msh.be{2} = be_ext;
    
    msh=rmfield(msh,'bn');
    msh.bn{1} = union(be_int(:,1),be_int(:,2));
    msh.bn{2} = union(be_ext(:,1),be_ext(:,2));
    period = 1;
    np0=size(msh.p,1);
    i=[1:np0]';
end

if (abs(cp(1,2)) < 10^-10)
    ic=find(sqrt((msh.p(:,1)-xc(1)).^2+(msh.p(:,2)-xc(2)).^2) < 10^-10);
    i(ic)=ic(1);
    [~, I]=histc(i,unique(i));
    msh.t=[renum(I,msh.t(:,1)) renum(I,msh.t(:,2)) renum(I,msh.t(:,3))]; 
    msh.p(ic(2:end),:)=[];
 
    % remove internal boundary
    if (period == 1)
        msh.be{1} = [renum(I,msh.be{1}(:,1)) renum(I,msh.be{1}(:,2))];
        msh.bn{1} = union(msh.be{1}(:,1),msh.bn{1}(:,2));
        msh.be(2) = [];
        msh.bn(2) = [];
    else
        msh.be(1) = []; msh.bn(1) = [];
        for k = 1:3
           msh.be{k} = [renum(I,msh.be{k}(:,1)) renum(I,msh.be{k}(:,2))];
           msh.bn{k} = union(msh.be{k}(:,1), msh.be{k}(:,2));
        end
        be = msh.be{3}; bn = msh.bn{3};
        msh.be{3} = msh.be{1}; msh.bn{3} = msh.bn{1};
        msh.be{1} = be; msh.bn{1} = bn;
    end
   
end

% delete elements/faces with duplicate nodes
it = min(diff(sort(msh.t,2),1,2),[],2); msh.t(it==0,:) = [];
clear it
nt = size(msh.be,1);
for k = 1:nt
    it = min(diff(sort(msh.be{k},2),1,2),[],2); msh.be{k}(it==0,:) = [];
end
 

    



