function segments=contourenvelope(md,varargin) %CONTOURENVELOPE - build a set of segments enveloping a contour .exp % % Usage: % segments=contourenvelope(md,varargin) % % Example: % segments=contourenvelope(md,'Stream.exp'); % segments=contourenvelope(md,md.mask.elementonfloatingice) % segments=contourenvelope(md); %some checks if nargin>2, help contourenvelope error('contourenvelope error message: bad usage'); end if nargin==2, flags=varargin{1}; if ischar(flags), file=flags; file=varargin{1}; if ~exist(file), error(['contourenvelope error message: file ' file ' not found']); end isfile=1; elseif isnumeric(flags), %do nothing for now isfile=0; else error('contourenvelope error message: second argument should a file or an elements flag'); end end %Now, build the connectivity tables for this mesh. %Computing connectivity if size(md.nodeconnectivity,1)~=md.numberofnodes, md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes); end if size(md.elementconnectivity,1)~=md.numberofelements, md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity); end %get nodes inside profile elementconnectivity=md.elementconnectivity; if nargin==2, if isfile, %get flag list of elements and nodes inside the contour nodein=ContourToMesh(md.elements,md.x,md.y,file,'node',1); elemin=(sum(nodein(md.elements),2)==size(md.elements,2)); %modify element connectivity elemout=find(~elemin); elementconnectivity(elemout,:)=0; elementconnectivity(find(ismember(elementconnectivity,elemout)))=0; else %get flag list of elements and nodes inside the contour nodein=zeros(md.numberofnodes,1); elemin=zeros(md.numberofelements,1); pos=find(flags); elemin(pos)=1; nodein(md.elements(pos,:))=1; %modify element connectivity elemout=find(~elemin); elementconnectivity(elemout,:)=0; elementconnectivity(find(ismember(elementconnectivity,elemout)))=0; end end %Find element on boundary %First: find elements on the boundary of the domain flag=elementconnectivity; if nargin==2, flag(find(flag))=elemin(flag(find(flag))); end elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0); %Find segments on boundary pos=find(elementonboundary); num_segments=length(pos); segments=zeros(num_segments,3); count=1; for i=1:num_segments, el1=pos(i); els2=elementconnectivity(el1,find(elementconnectivity(el1,:))); if length(els2)>1, flag=intersect(md.elements(els2(1),:),md.elements(els2(2),:)); nods1=md.elements(el1,:); nods1(find(nods1==flag))=[]; segments(count,:)=[nods1 el1]; ord1=find(nods1(1)==md.elements(el1,:)); ord2=find(nods1(2)==md.elements(el1,:)); %swap segment nodes if necessary if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), temp=segments(count,1); segments(count,1)=segments(count,2); segments(count,2)=temp; end segments(count,1:2)=fliplr(segments(count,1:2)); count=count+1; else nods1=md.elements(el1,:); flag=setdiff(nods1,md.elements(els2,:)); for j=1:3, nods=nods1; nods(j)=[]; if any(ismember(flag,nods)), segments(count,:)=[nods el1]; ord1=find(nods(1)==md.elements(el1,:)); ord2=find(nods(2)==md.elements(el1,:)); if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), temp=segments(count,1); segments(count,1)=segments(count,2); segments(count,2)=temp; end segments(count,1:2)=fliplr(segments(count,1:2)); count=count+1; end end end end