Index: sm/trunk-jpl/src/m/exp/exp_to_levelset.m
===================================================================
--- /issm/trunk-jpl/src/m/exp/exp_to_levelset.m	(revision 20179)
+++ 	(revision )
@@ -1,143 +1,0 @@
-function  levelset=exp_to_levelset(md,contourname)
-%EXP_TO_LEVELSET - from an exp contour file, build a levelset which has positive sign inside the contour, and negative outside.
-%
-%   Usage:
-%      levelset=exp_to_levelset(md,contourname)
-%
-%   Where: 
-%      - 'md'     : model specifying the mesh
-%      - 'contourname'  : exp file which defines the closed contour
-%      - 'levelset'  : levelset vector on each vertex of the mesh, size md.mesh.numberofvertices
-%
-%   Example:
-%      level=exp_to_levelset(md,'DomainOutline.exp');
-%
-%   See also MESHPROFILEINTERSECTION, SHPREAD, EXPREAD
-	
-	tolerance=1e-5; %tolerance level used for detection of edges, points coinciding, etc ...
-		
-	%from exp contour, build list of segments throughout our mesh:
-	segments=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,contourname);
-
-	%cleanup: remove 0 length segments: 
-	cleanup=1;
-	if cleanup,
-		flags=zeros(size(segments,1),1);
-		for j=1:size(segments,1),
-			segment=segments(j,:);
-			x1=segment(1); x2=segment(3); 
-			y1=segment(2); y2=segment(4); 
-			Q1=[x1;y1];Q2=[x2;y2];
-			if norm(Q2-Q1)==0, 
-				flags(j)=1;
-			end
-		end
-		pos=find(flags);
-		segments(pos,:)=[];
-	end
-
-	%simplify segments. levelsets are P1 through our triangular elements, which means we can only have one segment per element. 
-	%We simplify by merging all segments into one for each element.
-	simplify=1;
-	if simplify,
-		ids=unique(segments(:,5));
-		newsegs=zeros(length(ids),5);
-		newsegs(:,5)=ids;
-		for j=1:length(ids),
-			pos=find(segments(:,5)==ids(j));
-			segs=segments(pos,1:4);
-			count=0;
-			while(size(segs,1)>1),
-				count=count+1;
-				if count>10,
-					%keep first segment cause something funky in the configuration is happening.
-					segs=segs(1,:);
-					break;
-				end
-				if (abs(segs(1,3)-segs(2,1))<tolerance & abs(segs(1,4)-segs(2,2))<tolerance),
-					%the second segment disappears, extending the first segment:
-					if size(segs,1)==2,
-						segs=[segs(1,1) segs(1,2) segs(2,3) segs(2,4)];
-					else
-						segs=[segs(1,1), segs(1,2), segs(2,3), segs(2,4); segs(3:end,:)];
-					end
-				else
-					%put the first segment at the end
-					segs=[segs(2:end,:); segs(1,:)];
-				end
-			end
-			newsegs(j,1:4)=segs;
-		end
-		segments=newsegs;
-	end
-
-	%cleanup: remove 0 length segments: 
-	if cleanup,
-		flags=zeros(size(segments,1),1);
-		for j=1:size(segments,1),
-			segment=segments(j,:);
-			x1=segment(1); x2=segment(3); 
-			y1=segment(2); y2=segment(4); 
-			Q1=[x1;y1];Q2=[x2;y2];
-			if norm(Q2-Q1)<tolerance,
-				flags(j)=1;
-			end
-		end
-		pos=find(flags);
-		segments(pos,:)=[];
-	end
-
-	%now, using the segments, create level set flags on each element that touches the segment.
-	distance=zeros(md.mesh.numberofvertices,1);
-	flags=zeros(md.mesh.numberofvertices,1);
-	for j=1:size(segments,1),
-		segment=segments(j,:);
-		id=segment(5);
-		x1=segment(1); x2=segment(3); 
-		y1=segment(2); y2=segment(4); 
-		for k=1:3,
-			node=md.mesh.elements(id,k);
-			xk=md.mesh.x(node); yk=md.mesh.y(node);
-			
-			%on which side of the segment is this node k? 
-			P=[xk;yk]; Q1=[x1;y1];Q2=[x2;y2];
-			d=(det([Q2-Q1,P-Q1]))/norm(Q2-Q1);
-			if(~flags(node)),
-				distance(node)=d;
-				flags(node)=1;
-			end
-		end
-	end
-		
-	%now go through all the other nodes and complete (similar to phase unwrapping: we march through the domain, and each time we 
-	%are connected to a node which has a levelset value, we adopt the same one.
-	complete=1;
-	if complete,
-		while(length(find(flags==0))),
-			sumf=sum(flags(md.mesh.elements),2);
-			pos=find(sumf==2);
-			if ~isempty(pos),
-				for j=1:length(pos),
-					el=pos(j);
-					f=flags(md.mesh.elements(el,:));
-					ids=find(f==1); 
-					d=distance(md.mesh.elements(el,ids));
-					if d(1)*d(2)<0,
-						%we have hit a snag here, can't decide!
-						continue;
-					else
-						d=d(1);
-					end
-					if abs(d)<tolerance,
-						continue;
-					end
-					ids=find(f==0);
-					if ~isempty(ids),
-						flags(md.mesh.elements(el,ids))=1;
-						distance(md.mesh.elements(el,ids))=d;
-					end
-				end
-			end
-		end
-	end
-	levelset=distance;