Index: /issm/trunk/src/m/classes/public/geography2.m
===================================================================
--- /issm/trunk/src/m/classes/public/geography2.m	(revision 1222)
+++ /issm/trunk/src/m/classes/public/geography2.m	(revision 1223)
@@ -1,3 +1,4 @@
 function md=geography2(md,landname,iceshelfname,icesheetname)
+
 %Get assigned fields
 x=md.x;
@@ -8,9 +9,5 @@
 [gridonland,elementonland]=ContourToMesh(elements,x,y,expread(landname,1),'element and node',2);
 
-%use segments to make sure gridonland is correct at the boundary:
-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:
+%any element with 3 grids on land should be on land:
 elementsonwater=find(~elementonland);
 wrongelements=elementsonwater(find(( gridonland(md.elements(elementsonwater,1)) + gridonland(md.elements(elementsonwater,2)) + gridonland(md.elements(elementsonwater,3)) ...
@@ -18,8 +15,13 @@
 elementonland(wrongelements)=1;
 
-%finally, any element with its barycentre in the water should be in the water!
-xelem=x(md.elements)*[1;1;1]/3; yelem=y(md.elements)*[1;1;1]/3;
+%any element with its barycentre on land should be on land:
+weights={[1;1;1],[2;1;1],[1;2;1],[1;1;2]};
+for i=1:length(weights),
+	xelem=x(md.elements)*weights{i}/sum(weights{i});
+	yelem=y(md.elements)*weights{i}/sum(weights{i});
+end
 baryonland=ContourToNodes(xelem,yelem,expread(landname,1),1);
 pos=find(~baryonland); elementonland(pos)=0;
+pos=find(baryonland); elementonland(pos)=1;
 
 %figure out which elements on land are actually in the middle of the ocean!
@@ -29,7 +31,13 @@
 connectedtolandsum=sum(connectedtoland,2);
 waterelements=pos1(find(connectedtolandsum==3));
+elementonland(waterelements)=0;
 
-%now use waterelements to correct elementonland: 
-elementonland(waterelements)=0;
+%figure out which elements on water  are actually in the middle of the land!
+pos1=find(~elementonland); 
+connectedtowater=md.elementconnectivity(pos1,:);
+pos=find(connectedtowater); connectedtowater(pos)=elementonland(connectedtowater(pos));
+connectedtowatersum=sum(connectedtowater,2);
+landelements=pos1(find(connectedtowatersum==3));
+elementonland(landelements)=1;
 
 %recover arrays of ice shelf grids and elements, and ice sheet grids and elements.
@@ -74,4 +82,43 @@
 elementonicesheet=double(~elementoniceshelf & ~elementonwater);
 
+
+
+%Deal with segments on neumann: 
+elementconnectivity=md.elementconnectivity;
+pos=find(elementconnectivity);
+elementconnectivity(pos)=elementonwater(elementconnectivity(pos));
+
+pos=find(elementonwater);
+elementconnectivity(pos,:)=0;
+
+num_segments=sum(sum(elementconnectivity));
+segments=zeros(num_segments,3);
+
+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
+
 %Return: 
 md.gridoniceshelf=gridoniceshelf;
@@ -84,32 +131,4 @@
 md.elementonicesheet=elementonicesheet;
 
-%Deal with segments on neumann: 
-pos=find(md.segmentmarkers==2); %internal segments.
-segments=md.segments(pos,:);
-
-%swap segments so they are all oriented the same way
-maxcon=size(md.nodeconnectivity,2);
-for i=1:length(segments),
-
-	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);
-	else
-		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),:));
-	
-	%swap segment grids if necessary
-	if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
-		temp=segments(i,1);
-		segments(i,1)=segments(i,2);
-		segments(i,2)=temp;
-	end
-	segments(i,1:2)=fliplr(segments(i,1:2));
-	
-end
 md.segmentonneumann_diag=segments;
 md.counter=2;