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