Index: /issm/trunk/src/m/classes/public/contourenvelope.m
===================================================================
--- /issm/trunk/src/m/classes/public/contourenvelope.m	(revision 2517)
+++ /issm/trunk/src/m/classes/public/contourenvelope.m	(revision 2517)
@@ -0,0 +1,82 @@
+function segments=contourenvelope(md,file)
+%CONTOURENVELOPE - build a set of segments enveloping a contour .exp
+%
+%   Usage:
+%      segments=contourenvelope(md,file)
+
+%some checks
+if ~exist(file),
+	error(['thicknessevolution error message: file ' file ' not found']);
+end
+
+%Now, build the connectivity tables for this mesh.
+%Computing connectivity
+if size(md.nodeconnectivity,1)~=md.numberofgrids,
+	md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+end
+if size(md.elementconnectivity,1)~=md.numberofelements,
+	md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+end
+
+%get nodes inside profile
+%get flag list of elements and nodes inside the contour
+nodein=ContourToMesh(md.elements,md.x,md.y,expread(file,1),'node',1);
+elemin=(sum(nodein(md.elements),2)==size(md.elements,2));
+%modify element connectivity
+elemout=find(~elemin);
+elementconnectivity=md.elementconnectivity;
+elementconnectivity(elemout,:)=0;
+elementconnectivity(find(ismember(elementconnectivity,elemout)))=0;
+
+%Find element on boundary
+%First: find elements on the boundary of the domain
+flag=elementconnectivity;
+flag(find(flag))=elemin(flag(find(flag)));
+elementonboundary=double(sum(flag,2)<3 & 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 grids 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