function [xyz,skip,dx2,dy2,tri]=genmesh(bounds,nx,ny) xyz=zeros(nx*ny,2); skip=zeros(nx*ny,1); count=0; for i=1:ny for j=1:nx p=nx*(i-1)+j; xyz(p,1)=(j-1)/(nx-1)*(bounds(2)-bounds(1))+bounds(1); xyz(p,2)=(i-1)/(ny-1)*(bounds(4)-bounds(3))+bounds(3); if norm([xyz(p,1) xyz(p,2)])<1 skip(p)=-1; count=count+1; elseif norm([xyz(p,1) xyz(p,2)])==1 skip(p)=0; else skip(p)=0; end end end skip=[skip; -1*ones(count,1)]; xyz=[xyz; zeros(count,2)]; %Initialize a Map dx2=containers.Map('KeyType','double', 'ValueType','double'); dy2=containers.Map('KeyType','double', 'ValueType','double'); xb= containers.Map('KeyType','double', 'ValueType','double'); yb= containers.Map('KeyType','double', 'ValueType','double'); node_count=nx*ny+1; index=nx*ny+1; for i=2:ny-1 for j=2:nx-1 p=nx*(i-1)+j; if skip(p)==-1 continue end if skip(p+1)==-1 skip(p)=skip(p)+1; x=-sqrt(1-xyz(p,2)^2); xb (p)=node_count; node_count=node_count+1; dx2(p)=x-xyz(p,1); xyz(index,:)=[x xyz(p,2)]; index=index+1; end if skip(p+nx)==-1 skip(p)=skip(p)+2; y=-sqrt(1-xyz(p,1)^2); yb (p)=node_count; node_count=node_count+1; dy2(p)=y-xyz(p,2); xyz(index,:)=[xyz(p,1) y]; index=index+1; end if skip(p-1)==-1 skip(p)=skip(p)+4; x=sqrt(1-xyz(p,2)^2); xb (p)=node_count; node_count=node_count+1; dx2(p)=x-xyz(p,1); xyz(index,:)=[x xyz(p,2)]; index=index+1; end if skip(p-nx)==-1 skip(p)=skip(p)+8; y=sqrt(1-xyz(p,1)^2); yb(p)=node_count; node_count=node_count+1; dy2(p)=y-xyz(p,2); xyz(index,:)=[xyz(p,1) y]; index=index+1; end end end tri=zeros(2*(nx-1)*(ny-1),3); %we don't know exactly how many there will be, but it is at least this index=1; for i=1:ny-1 for j=1:nx-1 p=nx*(i-1)+j; selem=[skip(p) skip(p+1) skip(p+nx+1) skip(p+nx)]; for k=1:4 if selem(k)>0 selem(k)=0; end end elem=[p p+1 p+nx+1 p+nx -7 -7]; if sum(selem)==0 || sum(selem)==-4 tri(index+0,:)=[elem(1),elem(2),elem(3)]; tri(index+1,:)=[elem(1),elem(3),elem(4)]; index=index+2; elseif sum(selem)==-1 || sum(selem)==-3 for g=1:4 if selem(g)==(-3-sum(selem))/2 spot=g; end end if spot==1 elseif spot==2 elem=[elem(2:4) elem(1)]; elseif spot==3 elem=[elem(3:4) elem(1:2)]; elseif spot==4 elem=[elem(4) elem(1:3)]; end if sum(selem)==-3 if spot==1 || spot==3 elem(5:6)= [xb(elem(1)) yb(elem(1))]; elseif spot==2 || spot==4 elem(5:6)=[yb(elem(1)) xb(elem(1))]; else disp('error') end end if sum(selem)==-1 if spot==1 || spot==3 elem(5:6)=[xb(elem(2)) yb(elem(4))]; elseif spot==2 || spot==4 elem(5:6)=[yb(elem(2)) xb(elem(4))]; else disp('error') end end tri(index+0,:)=[elem(1),elem(5),elem(6)]; tri(index+1,:)=[elem(5),elem(2),elem(3)]; tri(index+2,:)=[elem(5),elem(3),elem(6)]; tri(index+3,:)=[elem(6),elem(3),elem(4)]; index=index+4; elseif sum(selem)==-2 if isequal(selem(1:2),[0 -1]) elem(5:6)=[xb(elem(1)) xb(elem(4))]; elseif isequal(selem(1:2),[0 0]) elem=[elem(2:4) elem(1)]; elem(5:6)=[yb(elem(1)) yb(elem(4))]; elseif isequal(selem(1:2),[-1 0]) elem=[elem(3:4) elem(1:2)]; elem(5:6)=[xb(elem(1)) xb(elem(4))]; elseif isequal(selem(1:2),[-1 -1]) elem=[elem(4) elem(1:3)]; elem(5:6)=[yb(elem(1)) yb(elem(4))]; end tri(index+0,:)=[elem(1),elem(5),elem(6)]; tri(index+1,:)=[elem(1),elem(6),elem(4)]; tri(index+2,:)=[elem(3),elem(5),elem(2)]; tri(index+3,:)=[elem(6),elem(5),elem(3)]; index=index+4; else disp('error') disp(elem) end end end end