contour_icefront

PURPOSE ^

CONTOUR_ICEFRONT - ???

SYNOPSIS ^

function ind=contour_icefront(index,x,y,node_on_icefront,element_on_boundary);

DESCRIPTION ^

CONTOUR_ICEFRONT - ???

   This program takes element_on_boundary and node_on_boundary 
   x and  y and index and makes an ordered boundary with it,
   turning counter clockwise.

   Usage:
      ind=contour_icefront(index,x,y,node_on_icefront,element_on_boundary);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function ind=contour_icefront(index,x,y,node_on_icefront,element_on_boundary);
0002 %CONTOUR_ICEFRONT - ???
0003 %
0004 %   This program takes element_on_boundary and node_on_boundary
0005 %   x and  y and index and makes an ordered boundary with it,
0006 %   turning counter clockwise.
0007 %
0008 %   Usage:
0009 %      ind=contour_icefront(index,x,y,node_on_icefront,element_on_boundary);
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 pos=find(temp);
0019 pos2=find(x(pos)==max(x(pos)));
0020 xb=x(pos(pos2));
0021 yb=y(pos(pos2));
0022 
0023 
0024 
0025    
0026 n00=find(x==xa & y==ya)
0027 n10=find(x==xb & y==yb)
0028 n1=n10;
0029 n0=n00;
0030 cont=[n0 n1];
0031 
0032    
0033 while (n1~=n00),
0034    
0035    els=find(index(:,1)==n1 | index(:,2)==n1 | index(:,3)==n1);
0036    possible_nodes=zeros(2*length(els),1);
0037    for n=1:length(els),
0038    other_pos=find(index(els(n),:)~=n1);
0039    possible_nodes(2*n-1,:)= index(els(n),other_pos(1));
0040    possible_nodes(2*n,:)=index(els(n),other_pos(2));
0041     end
0042 potential_node=possible_nodes(1);
0043 
0044 n0n1=[x(n1)-x(n0)
0045    y(n1)-y(n0)];
0046 n0n1_unit=n0n1/norm(n0n1);
0047 ny=flipud(n0n1_unit);
0048 ny(1)=-ny(1);
0049 possible_angles=zeros(length(possible_nodes),1);
0050 
0051 for n=1:length(possible_nodes),
0052    if (possible_nodes(n)~=n0),
0053         u=[x(possible_nodes(n))-x(n1)
0054           y(possible_nodes(n))-y(n1)];
0055           ux=sum(u.*n0n1_unit);
0056           uy=sum((u-ux*n0n1_unit).*ny);
0057          possible_angles(n)=atan2(uy,ux);
0058       else
0059          possible_angles(n)=NaN;
0060    end
0061 end
0062 pos=find(possible_angles==min(possible_angles));
0063 potential_node=possible_nodes(pos(1));
0064 
0065 cont=[cont
0066    [n1 potential_node]];
0067 n0=n1;
0068 n1=potential_node;
0069 end
0070 
0071 ind=cont(2:length(cont),1);
0072 
0073

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003