Index: /issm/trunk-jpl/src/m/exp/exp_to_levelset.m
===================================================================
--- /issm/trunk-jpl/src/m/exp/exp_to_levelset.m	(revision 18773)
+++ /issm/trunk-jpl/src/m/exp/exp_to_levelset.m	(revision 18773)
@@ -0,0 +1,143 @@
+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(length(segments),1);
+		for j=1:length(segments),
+			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(length(segments),1);
+		for j=1:length(segments),
+			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;