Index: /issm/trunk/src/m/classes/public/geography2.m
===================================================================
--- /issm/trunk/src/m/classes/public/geography2.m	(revision 1125)
+++ /issm/trunk/src/m/classes/public/geography2.m	(revision 1126)
@@ -1,4 +1,3 @@
 function md=geography2(md,landname,iceshelfname,icesheetname)
-
 %Get assigned fields
 x=md.x;
@@ -7,10 +6,10 @@
 
 %recover elements and grids on land.
-[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',1);
+[gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
 
 %use segments to make sure gridonland is correct at the boundary:
-gridonland(md.segments(:,1))=1;
-gridonland(md.segments(:,2))=1;
-
+pos=find(md.segmentmarkers==2);
+gridonland(md.segments(pos,1))=1;
+gridonland(md.segments(pos,2))=1;
 %now, any elements with 3 grids on land should be on land:
 elementsonwater=find(~elementonland);
@@ -42,5 +41,5 @@
 	elementoniceshelf=ones(md.numberofelements,1);
 else
-	[gridoniceshelf,elementoniceshelf]=ContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',1);
+	[gridoniceshelf,elementoniceshelf]=ContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',2);
 end
 
@@ -49,5 +48,5 @@
 	elementonicesheet=zeros(md.numberofelements,1);
 else
-	[gridonicesheet,elementonicesheet]=ContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',1);
+	[gridonicesheet,elementonicesheet]=ContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',2);
 end
 
@@ -65,4 +64,5 @@
 gridonwater=double(~gridonland);
 elementonwater=double(~elementonland);
+
 
 %correct for islands:
@@ -89,8 +89,9 @@
 
 %swap segments so they are all oriented the same way
+maxcon=size(md.nodeconnectivity,2);
 for i=1:length(segments),
-	pair=find(  (md.elements(:,1)==segments(i,1)  | md.elements(:,2)==segments(i,1)  | md.elements(:,3)==segments(i,1) ) & ...
-	            (md.elements(:,1)==segments(i,2)  | md.elements(:,2)==segments(i,2)  | md.elements(:,3)==segments(i,2) )   ...
-			 );
+
+	pair=intersect(md.nodeconnectivity(segments(i,1),1:md.nodeconnectivity(segments(i,1),maxcon)),md.nodeconnectivity(segments(i,2),1:md.nodeconnectivity(segments(i,2),maxcon)));
+
     if md.elementonwater(pair(1)),
 		segments(i,3)=pair(2);
@@ -98,4 +99,5 @@
 		segments(i,3)=pair(1);
 	end
+	
 	ord1=find(segments(i,1)==md.elements(segments(i,3),:));
 	ord2=find(segments(i,2)==md.elements(segments(i,3),:));
@@ -108,8 +110,7 @@
 	end
 	segments(i,1:2)=fliplr(segments(i,1:2));
+	
 end
 md.segmentonneumann_diag=segments;
 md.counter=2;
 md.segmentmarkers(:)=1;
-
-