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