0001 function ind=contour_maker(index,x,y,node_on_boundary,element_on_boundary);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 pos=find(node_on_boundary);
0012 boundary_chaos=[x(pos) y(pos)];
0013 ind=convhull(boundary_chaos(:,1),boundary_chaos(:,2));
0014
0015 xa=boundary_chaos(ind(1),1);
0016 xb=boundary_chaos(ind(2),1);
0017 ya=boundary_chaos(ind(1),2);
0018 yb=boundary_chaos(ind(2),2);
0019
0020
0021
0022 n00=find(x==xa & y==ya);
0023 n10=find(x==xb & y==yb);
0024 n1=n10;
0025 n0=n00;
0026 cont=[n0 n1];
0027
0028
0029 while (n1~=n00),
0030
0031 els=find(index(:,1)==n1 | index(:,2)==n1 | index(:,3)==n1);
0032 possible_nodes=zeros(2*length(els),1);
0033 for n=1:length(els),
0034 other_pos=find(index(els(n),:)~=n1);
0035 possible_nodes(2*n-1,:)= index(els(n),other_pos(1));
0036 possible_nodes(2*n,:)=index(els(n),other_pos(2));
0037 end
0038 potential_node=possible_nodes(1);
0039
0040 n0n1=[x(n1)-x(n0)
0041 y(n1)-y(n0)];
0042 n0n1_unit=n0n1/norm(n0n1);
0043 ny=flipud(n0n1_unit);
0044 ny(1)=-ny(1);
0045 possible_angles=zeros(length(possible_nodes),1);
0046
0047 for n=1:length(possible_nodes),
0048 if (possible_nodes(n)~=n0),
0049 u=[x(possible_nodes(n))-x(n1)
0050 y(possible_nodes(n))-y(n1)];
0051 ux=sum(u.*n0n1_unit);
0052 uy=sum((u-ux*n0n1_unit).*ny);
0053 possible_angles(n)=atan2(uy,ux);
0054 else
0055 possible_angles(n)=NaN;
0056 end
0057 end
0058 pos=find(possible_angles==min(possible_angles));
0059 potential_node=possible_nodes(pos(1));
0060
0061 cont=[cont
0062 [n1 potential_node]];
0063 n0=n1;
0064 n1=potential_node;
0065 end
0066
0067 ind=cont(2:length(cont),1);
0068
0069