source: issm/trunk/src/m/model/contourenvelope.m@ 9725

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.mesh.numberofvertices,
        md.nodeconnectivity=NodeConnectivity(md.elements,md.mesh.numberofvertices);
end
if size(md.elementconnectivity,1)~=md.mesh.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.mesh.numberofvertices,1); 
                elemin=zeros(md.mesh.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