function msh=kmg3ddu(cp,msize,xc)
%KMG3DU Uniform symmetric mesh of Omega = (ax,bx)*(ay,by)*(az,bz)
% or, if xc (the center coordinates) is given, 
% Omega=(theta1,theta2)*(phi1,phi2)*(r1,r2) in spherical coordinates.
% With spherical coordinates, KMG3DU generates a tetrahedral mesh of 
% a spherical domain (ark, dome, 3D ring, sphere). If r1=0 or theta1=0 or 
% theta2=pi or phi2=2*pi, the duplicated nodes and boundaries triangles are 
% removed from the mesh.

% Boundary sides are numbered as follows :
% s1 side z=az (lower side); s2 side z=bz (upper side) 
% s3 side y=ay (front side); s4 side x=bx (right side)
% s5 side y=by (back side);  s6 side x=ax (left side)
%--------------------------------------------------------------------
%      msh=kmg3du(cp,msize)
%      msh=kmg3du(cp,msize,xc)
% 
% Input:
%   cp      : Corner points cp=[xmin,ymin,zmin
%                               xmax,ymax,zmax]
%           if xc is given, then cp must be in spherical coordinates
%           cp=[theta1, phi1, r1
%               theta2, phi2, r2]
%           0 <= theta1 < theta2 <= pi
%           0 <= phi1 < phi2 <= 2*pi
%           0 <= r1 < r2
%   msize   : Number of points in each direction msize=[nx
%                                                       ny
%                                                       nz]        
%
% Ouputs:
%   msh structure with fields  
%      p    : Nodes coordinates, array np*3
%      t    : Tetrahedra vertices, array nt*4
%     bt    : Boundary triangles 6-by-1 cell-array with 
%             bt{i}(:,1:3) triangle vertices of side si (cartesian domain)
%             or 5-by-1, 4-by-1, 3-by-1, 2-by-1, or 1 (sphere)
%     bn    : Boundary nodes 6-by-1 cell-array with bn{i} nodes of side si
%             or 5-by-1, 4-by-1, 3-by-1, 2-by-1, or 1 (sphere)
%-------------------------------------------------------------------------
% (c) 2013-2021, KOKO J., LIMOS/ISIMA, Jonas.Koko@isima.fr
%-------------------------------------------------------------------------
eps0=10^-10;
ax = cp(1,1); ay = cp(1,2); az = cp(1,3);  
bx = cp(2,1); by = cp(2,2); bz = cp(2,3);
nx = msize(1); ny = msize(2); nz = msize(3); 

nnx =round((nx-1)/2)+1; nny=round((ny-1)/2)+1; nnz=round((nz-1)/2)+1;
np=nnx*nny*nnz; nq=(nnx-1)*(nny-1)*(nnz-1); nt=6*nq; 
 
% nodes coordinates
hx=(bx-ax)/(nnx-1); hy=(by-ay)/(nny-1); hz=(bz-az)/(nnz-1);
[x,y,z]=ndgrid(ax:hx:bx,ay:hy:by,az:hz:bz);
p=[x(:),y(:),z(:)];

% quadrangles
ip=[1:nnx*nny*nnz]'; ijk=reshape(ip,[nnx,nny,nnz]);
b1=squeeze(ijk(:,1,:));  ib1=b1(:);
b2=squeeze(ijk(nnx,:,:)); ib2=b2(:);
b3=squeeze(ijk(:,nny,:)); ib3=b3(:);
b4=squeeze(ijk(1,:,:));  ib4=b4(:);
b5=squeeze(ijk(:,:,1));  ib5=b5(:);
b6=squeeze(ijk(:,:,nnz)); ib6=b6(:);

ib236=union(ib2,union(ib3,ib6)); iq1=setdiff(ip,ib236);
ib346=union(ib3,union(ib4,ib6)); iq2=setdiff(ip,ib346);
ib345=union(ib3,union(ib4,ib5)); iq3=setdiff(ip,ib345);
ib235=union(ib2,union(ib3,ib5)); iq4=setdiff(ip,ib235);
ib126=union(ib1,union(ib2,ib6)); iq5=setdiff(ip,ib126);
ib146=union(ib1,union(ib4,ib6)); iq6=setdiff(ip,ib146);
ib145=union(ib1,union(ib4,ib5)); iq7=setdiff(ip,ib145);
ib125=union(ib1,union(ib2,ib5)); iq8=setdiff(ip,ib125);

% extract edges & faces
eh0=[iq1 iq2
    iq2 iq3
    iq3 iq4
    iq4 iq1
    iq5 iq6
    iq6 iq7
    iq7 iq8
    iq8 iq5
    iq1 iq5
    iq2 iq6
    iq3 iq7
    iq4 iq8];

fh0=[iq1 iq2 iq3 iq4
     iq5 iq6 iq7 iq8
     iq1 iq2 iq6 iq5
     iq2 iq3 iq7 iq6
     iq3 iq4 iq8 iq7
     iq1 iq4 iq8 iq5];

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

[fh,ih]=sort(fh0,2);
%ff=int64(fh(:,1)*(np+1)+fh(:,2)+fh(:,3)+fh(:,4));
[~,iff,jf]=unique(fh,'rows'); 
f=fh(iff,:);

clear eh0 fh0 ih ee ff iee iff 

% Form new nodes (mid-point of edges)
ne=size(e,1); nf=size(f,1);
nef=ne+nf;
p0=p;

msh.p=[p0 
       (p0(e(:,1),:)+p0(e(:,2),:))/2
       (p0(f(:,1),:)+p0(f(:,2),:)+p0(f(:,3),:)+p0(f(:,4),:))/4
       (p0(iq1,:)+p0(iq2,:)+p0(iq3,:)+p0(iq4,:)+p0(iq5,:)+p0(iq6,:)...
                                           +p0(iq7,:)+p0(iq8,:))/8];
                                       
% Form new edges
pme=[(np+1):(np+ne)]'; lmpe=pme(je);
pmf=[np+ne+1:np+ne+nf]'; lmpf=pmf(jf);
ime1=lmpe(1:nq);        ime2=lmpe(nq+1:2*nq); 
ime3=lmpe(2*nq+1:3*nq); ime4=lmpe(3*nq+1:4*nq);
ime5=lmpe(4*nq+1:5*nq); ime6=lmpe(5*nq+1:6*nq);
ime7=lmpe(6*nq+1:7*nq); ime8=lmpe(7*nq+1:8*nq);
ime9=lmpe(8*nq+1:9*nq); ime10=lmpe(9*nq+1:10*nq);
ime11=lmpe(10*nq+1:11*nq); ime12=lmpe(11*nq+1:12*nq);

imf1=lmpf(1:nq);        imf2=lmpf(nq+1:2*nq);
imf3=lmpf(2*nq+1:3*nq); imf4=lmpf(3*nq+1:4*nq);
imf5=lmpf(4*nq+1:5*nq); imf6=lmpf(5*nq+1:6*nq);
imqc=[np+ne+nf+1:np+ne+nf+nq]';

% triangulation
t1=kuhntriang(imqc,imf6,ime9,imf3,imf1,ime4,iq1,ime1);
t2=kuhntriang(imqc,imf3,ime10,imf4,imf1,ime1,iq2,ime2);
t3=kuhntriang(imqc,imf4,ime11,imf5,imf1,ime2,iq3,ime3);
t4=kuhntriang(imqc,imf5,ime12,imf6,imf1,ime3,iq4,ime4);
t5=kuhntriang(imqc,imf6,ime9,imf3,imf2,ime8,iq5,ime5);
t6=kuhntriang(imqc,imf3,ime10,imf4,imf2,ime5,iq6,ime6);
t7=kuhntriang(imqc,imf4,ime11,imf5,imf2,ime6,iq7,ime7);
t8=kuhntriang(imqc,imf5,ime12,imf6,imf2,ime7,iq8,ime8);

msh.t=[t1
       t2
       t3
       t4
       t5
       t6
       t7
       t8];

clear iq1 iq2 iq3 iq4 iq5 iq6 iq7 iq8
clear p0 ime1 ime2 ime3 ime4 ime5 ime6 ime7 ime8 ime8 ime10 ime11 ime12
clear ifm1 ifm2 ifm3 ifm4 imf5 imqc

% Extract boundary faces
f0=[msh.t(:,1) msh.t(:,2) msh.t(:,3)
    msh.t(:,1) msh.t(:,2) msh.t(:,4)
    msh.t(:,1) msh.t(:,4) msh.t(:,3)
    msh.t(:,2) msh.t(:,3) msh.t(:,4)];

fs=sort(f0,2);
[~,iff,jf]=unique(fs,'rows'); 
f=fs(iff,:);
 
nf=size(f,1); 
he=accumarray(jf,1,[nf 1]); ip=find(he==1);
f = f(ip,:); i1=f(:,1); i2=f(:,2); i3=f(:,3);


np=size(msh.p,1);
ref=zeros(np,1);

iz1 = find(abs(msh.p(:,3)-az)<eps0); ref(iz1) = 1; %size(iz1)
if1 = find(ref(i1) == 1 & ref(i2) == 1 & ref(i3) == 1); %size(if1)
ref(iz1)=0;

iz2 = find(abs(msh.p(:,3)-bz)<eps0); ref(iz2) = 2; %size(iz2)
if2 = find(ref(i1) == 2 & ref(i2) == 2 & ref(i3) == 2);
ref(iz2)=0;

iy1 = find(abs(msh.p(:,2)-ay)<eps0); ref(iy1) = 3;
if3 = find(ref(i1) == 3 & ref(i2) == 3 & ref(i3) == 3);
ref(iy1)=0;

ix2 = find(abs(msh.p(:,1)-bx)<eps0); ref(ix2) = 4;
if4 = find(ref(i1) == 4 & ref(i2) == 4 & ref(i3) == 4);
ref(ix2)=0;

iy2 = find(abs(msh.p(:,2)-by)<eps0); ref(iy2) = 5;
if5 = find(ref(i1) == 5 & ref(i2) == 5 & ref(i3) == 5);
ref(iy2)=0;

ix1 = find(abs(msh.p(:,1)-ax)<eps0); ref(ix1) = 6;
if6 = find(ref(i1) == 6 & ref(i2) == 6 & ref(i3) == 6);

% Form the boundary triangles/side
msh.bt = {f(if1,:)
          f(if2,:)
          f(if3,:)
          f(if4,:)
          f(if5,:)
          f(if6,:)};

% Boundary nodes/side
msh.bn = {union(f(if1,1),union(f(if1,2),f(if1,3)))
          union(f(if2,1),union(f(if2,2),f(if2,3)))
          union(f(if3,1),union(f(if3,2),f(if3,3)))
          union(f(if4,1),union(f(if4,2),f(if4,3)))
          union(f(if5,1),union(f(if5,2),f(if5,3)))
          union(f(if6,1),union(f(if6,2),f(if6,3)))};
      
if (nargin == 2)
    return
end

clear if1 if2 if3 if4 if5 if6

% convert spherical -> cartesians
th = msh.p(:,1); phi = msh.p(:,2); r = msh.p(:,3);
x = xc(1)+r.*sin(th).*cos(phi); 
y = xc(2)+r.*sin(th).*sin(phi);
z = xc(3)+r.*cos(th);
msh.p(:,1) = x; msh.p(:,2) = y; msh.p(:,3) = z;

np = size(msh.p,1); 
i = [1:np]';

%return
% north pole (bt{6})
northp = 0;
if (abs(cp(1,1))<10^-10)
    ip = msh.bn{6};
    [~,ia,ja]=uniquetol(msh.p(ip,:),10^-10,'ByRows',true);
    ip1 = ip(ia); ip2 = setdiff(ip,ip1);
    i(ip) = ip1(ja);
    [~, I]=histc(i,unique(i));
    msh.t=[renum(I,msh.t(:,1)) renum(I,msh.t(:,2)) renum(I,msh.t(:,3)) renum(I,msh.t(:,4))]; 
    msh.p(ip2,:)=[];
    
    for k = 1:5
        msh.bt{k} = [renum(I,msh.bt{k}(:,1)) renum(I,msh.bt{k}(:,2)) renum(I,msh.bt{k}(:,3))];
        msh.bn{k} = union(msh.bt{k}(:,3), union(msh.bt{k}(:,2),msh.bt{k}(:,1)));
    end
    msh.bt(6) = []; msh.bn(6) = [];
    np = size(msh.p,1); i = [1:np]';
    fprintf('north-pole np=%d  max-num=%d %d %d %d \n',np,max(msh.t));
    northp = 1;
end

% south pole (bt{4})
southp = 0;
if (abs(cp(2,1)-pi)<10^-10)
    ip = msh.bn{4};
    [~,ia,ja]=uniquetol(msh.p(ip,:),10^-10,'ByRows',true);
    ip1 = ip(ia); ip2 = setdiff(ip,ip1);
    i(ip) = ip1(ja); 
    %i(ip) = ip(1);
    [~, I]=histc(i,unique(i));
    msh.t=[renum(I,msh.t(:,1)) renum(I,msh.t(:,2)) renum(I,msh.t(:,3)) renum(I,msh.t(:,4))]; 
    msh.p(ip2,:)=[];
    
    msh.bt(4) = []; msh.bn(4) = [];
    if (northp == 1)
        nt = 4; 
    else 
        nt = 5;
    end
    for k = 1:nt
        msh.bt{k} = [renum(I,msh.bt{k}(:,1)) renum(I,msh.bt{k}(:,2)) renum(I,msh.bt{k}(:,3))];
        msh.bn{k} = union(msh.bt{k}(:,3), union(msh.bt{k}(:,2),msh.bt{k}(:,1)));
    end
    np = size(msh.p,1); i = [1:np]';
    southp = 1;
    fprintf('south-pole np=%d  max-num=%d %d %d %d \n',np,max(msh.t));
end

% circular form (bt{3} = bt{5})
if (abs(cp(1,2))<10^-10 & abs(cp(2,2)-2*pi)<10^-10)
    
    fprintf('circular-dom0 np=%d  max-num=%d %d %d %d \n',size(msh.p,1),max(msh.t));
    
    ii=msh.bn{3}; 
    if (southp == 0)
        jj=msh.bn{5};
    else
        jj=msh.bn{4};
    end
 
    ip = union(ii,jj);
    [~,ia,ja]=uniquetol(msh.p(ip,:),10^-10,'ByRows',true);
    ip1 = ip(ia); ip2 = setdiff(ip,ip1);
    i(ip) = ip1(ja);
    %[length(ip) length(ii) length(jj) length(ip1)]
    
    [~, I]=histc(i,unique(i)); 
    %[max(i) max(I)]
    msh.t=[renum(I,msh.t(:,1)) renum(I,msh.t(:,2)) renum(I,msh.t(:,3)) renum(I,msh.t(:,4))];
    msh.p(ip2,:)=[];
    
    if (southp == 0)
        msh.bt(5) = []; msh.bn(5) = [];
    else
        msh.bt(4) = []; msh.bn(4) = [];
    end
    msh.bt(3) = []; msh.bn(3) = [];
    nt = size(msh.bt,1);
    for k = 1:nt
        msh.bt{k} = [renum(I,msh.bt{k}(:,1)) renum(I,msh.bt{k}(:,2)) renum(I,msh.bt{k}(:,3))];
        msh.bn{k} = union(msh.bt{k}(:,3), union(msh.bt{k}(:,2),msh.bt{k}(:,1)));
    end
     
    np=size(msh.p,1); i=[1:np]';
    %fprintf('circular-dom np=%d  max-num=%d %d %d %d \n',np,max(msh.t));
    
end

% center (delete bt{1})
if (abs(cp(1,3))< 10^-10)
    ip = msh.bn{1};
    max(ip)%
    %[~,ia,ja]=uniquetol(msh.p(ip,:),10^-10,'ByRows',true);
    %ip1 = ip(ia); ip2 = setdiff(ip,ip1);
    %i(ip) = ip1(ja);
    i(ip) = ip(1); 
    [~, I]=histc(i,unique(i));
    %max(I)
    %max(msh.t(:,1))
    msh.t=[renum(I,msh.t(:,1)) renum(I,msh.t(:,2)) renum(I,msh.t(:,3)) renum(I,msh.t(:,4))];
    msh.p(ip(2:end),:)=[];
    
    msh.bt(1) = []; msh.bn(1) = [];
    nt = size(msh.bt,1);
    for k = 1:nt
        msh.bt{k} = [renum(I,msh.bt{k}(:,1)) renum(I,msh.bt{k}(:,2)) renum(I,msh.bt{k}(:,3))];
        msh.bn{k} = union(msh.bt{k}(:,3), union(msh.bt{k}(:,2),msh.bt{k}(:,1)));
    end
    
    np=size(msh.p,1);
    fprintf('center np=%d  max-num=%d %d %d %d \n',np,max(msh.t));
    
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.bt,1);
for k = 1:nt
    it = min(diff(sort(msh.bt{k},2),1,2),[],2); msh.bt{k}(it==0,:) = [];
end
 
end 

% Kuhn triangulation of a cube into 6 tetrahedra
function t=kuhntriang(i1,i2,i3,i4,i5,i6,i7,i8)
  t=[i1 i2 i6 i7
     i1 i6 i5 i7
     i1 i2 i7 i3
     i1 i4 i3 i7
     i1 i5 i8 i7
     i1 i4 i7 i8];
end