Index: /issm/trunk/src/m/classes/public/geography2.m
===================================================================
--- /issm/trunk/src/m/classes/public/geography2.m	(revision 3115)
+++ /issm/trunk/src/m/classes/public/geography2.m	(revision 3116)
@@ -25,4 +25,12 @@
 else
 	error('Invalid area option option');
+end
+
+%Now, build the connectivity tables for this mesh.
+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
 
@@ -87,5 +95,4 @@
 elementonwater=double(~elementonland);
 
-
 %correct for islands:
 gridoniceshelf=double(gridoniceshelf & ~gridonicesheet);
@@ -96,43 +103,22 @@
 elementonicesheet=double(~elementoniceshelf & ~elementonwater);
 
+%Deal with segments on neumann:
 
+%Get current connectivity
+elementconnectivity=md.elementconnectivity;
 
-%Deal with segments on neumann: 
-elementconnectivity=md.elementconnectivity;
+%put 0 for elements on water
 pos=find(elementconnectivity);
-elementconnectivity(pos)=elementonwater(elementconnectivity(pos));
+elementconnectivity(pos)=elementconnectivity(pos).*(~elementonwater(elementconnectivity(pos)));
 
+%put line of ones for elements on water
 pos=find(elementonwater);
-elementconnectivity(pos,:)=0;
+elementconnectivity(pos,:)=1;% line of ones for elements on water so they won't be considered
 
-num_segments=sum(sum(elementconnectivity));
-segments=zeros(num_segments,3);
+%resort lines (zeros must be at the last column for findsegments)
+elementconnectivity=sort(elementconnectivity,2,'descend');
 
-icefrontelementsonland=find(sum(elementconnectivity,2));
-
-count=1;
-for i=1:numel(icefrontelementsonland),
-	el1=icefrontelementsonland(i);
-	els2=md.elementconnectivity(el1,find(elementconnectivity(el1,:)));
-	for j=1:numel(els2),
-		el2=els2(j);
-		%we have a segment between land el1 and water el2:
-		segments(count,:)=[intersect(md.elements(el1,:),md.elements(el2,:)) el1];
-		
-		ord1=find(segments(count,1)==md.elements(el1,:));
-		ord2=find(segments(count,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;
-	end
-end
-
+%call findsegments to build segment using THIS conectivity
+md.segments=findsegments(md,'elementconnectivity',elementconnectivity);
 
 %some final checks: