Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/augment2dmesh.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/augment2dmesh.m	(revision 21362)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/augment2dmesh.m	(revision 21363)
@@ -19,5 +19,5 @@
 %The innner segments of mhband and the outer segments of mh are identical. Go into  the elements of 
 %mhband and set them to their md1 equivalent: 
-tol=1; %1 meter 
+tol=getfieldvalue(options,'tolerance',1); %1 meter 
 for i=1:length(mhband.segments),
 	node2=mhband.segments(i,1);
@@ -30,4 +30,6 @@
 			pos=find(mhband.elements==node2); mhband.elements(pos)=node1;
 			segs=mhband.segments(:,1:2); pos=find(segs==node2); segs(pos)=node1; mhband.segments(:,1:2)=segs;
+			%flag the segment for future erasal: 
+			pos=find(mhband.segments(:,1)==node1 | mhband.segments(:,2)==node1); mhband.segments(pos,3)=NaN;
 			break;
 		end
@@ -39,7 +41,10 @@
 mh.x=[mh.x;mhband.x];
 mh.y=[mh.y;mhband.y];
-mh.lat=[mh.lat;mhband.lat];
-mh.long=[mh.long;mhband.long];
-mh.segments=[mh.segments;mhband.segments];
+if ~isnan(mh.lat), mh.lat=[mh.lat;mhband.lat]; end
+if ~isnan(mh.long), mh.long=[mh.long;mhband.long];  end
+mh.segments=[mhband.segments]; 
+%get rid of segments that are internal: 
+pos=find(~isnan(mh.segments(:,3)));
+mh.segments=mh.segments(pos,:); 
 
 %Remove orphans:
@@ -53,6 +58,6 @@
 	x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)];
 	y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)];
-	lat=[lat(1:orphan(i)-(i-1)-1); lat(orphan(i)-(i-1)+1:end)];
-	long=[long(1:orphan(i)-(i-1)-1); long(orphan(i)-(i-1)+1:end)];
+	if ~isnan(lat), lat=[lat(1:orphan(i)-(i-1)-1); lat(orphan(i)-(i-1)+1:end)]; end
+		if ~isnan(long), long=[long(1:orphan(i)-(i-1)-1); long(orphan(i)-(i-1)+1:end)]; end
 	%update elements
 	pos=find(elements>orphan(i)-(i-1));
@@ -73,2 +78,10 @@
 mh.numberofelements=length(mh.elements);
 mh.numberofvertices=length(mh.x);
+
+%vertex on boundaries: 
+mh.vertexonboundary=zeros(mh.numberofvertices,1);
+mh.vertexonboundary(mh.segments(:,1:2))=1;
+
+%connectivities: 
+mh.vertexconnectivity=NodeConnectivity(mh.elements,mh.numberofvertices);
+mh.elementconnectivity=ElementConnectivity(mh.elements,mh.vertexconnectivity);