Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24854)
+++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24855)
@@ -84,7 +84,11 @@
 
 		end % }}}
-		function mesh_connectivity(self,mesh) % {{{
-
-			plotr=1;
+		function mesh_connectivity(self,mesh,varargin) % {{{
+			
+			%retrieve options: 
+			options=pairoptions(varargin{:}); 
+
+			subsetregions=getfieldvalue(options,'regions',1:length(self.boxes));
+			errornotfound=getfieldvalue(options,'errornotfound',1);
 
 			%The way we run this is through the O2 zones defined in boxes. We go through 
@@ -119,5 +123,5 @@
 
 			%Go through O2 regions: 
-			for i=1:length(self.boxes),
+			for i=subsetregions,
 				string=self.boxes(i).RGI_CODE; 
 				disp(['progressing with region ' num2str(i) ' ' string]);
@@ -138,15 +142,15 @@
 					radius=sqrt( (lat0-minlat)^2+(long0-maxlong)^2);
 
-		%some radius adjustment:  {{{
-		switch i,
-			case 4, radius=40;
-			case 12, radius=25;
-			case 33, radius=5;
-			case 41, radius=75;
-			case 42, radius=45;
-			case 68, radius=10;
-			case 82, radius=30;
-			otherwise,
-		end % }}}
+					%some radius adjustment:  {{{
+					switch i,
+						case 4, radius=40;
+						case 12, radius=25;
+						case 33, radius=5;
+						case 41, radius=75;
+						case 42, radius=45;
+						case 68, radius=10;
+						case 82, radius=30;
+						otherwise,
+					end % }}}
 
 					[lids,rids]=self.gidtolid(glaciers);
@@ -197,8 +201,10 @@
 							end
 						end
-						self.glacier_connectivity(glaciers(j))=el;
-						if ~found, 
-							error(sprintf('could not find element for glacier %i with lid %i',j,lids(j)));
+						if ~found,
+							if errornotfound,
+								error(sprintf('could not find element for glacier %i with lid %i',j,lids(j)));
+							end
 						end
+						if(found)self.glacier_connectivity(glaciers(j))=el; end;
 					end
 				end
@@ -208,4 +214,58 @@
 			for j=1:length(self.glacier_connectivity),
 				el=self.glacier_connectivity(j);
+				if ~el,continue; end;
+				count=self.element_connectivity(el,ny);
+				if count>ny,
+					error('need to enlarge connectivity table to at least');
+				end
+				self.element_connectivity(el,count+1)=j;
+				self.element_connectivity(el,ny)=count+1;
+			end
+
+			%Reduce the number of columns (we did not initially, so we chose an arbitrary ny:
+			nmax=max(self.element_connectivity(:,end));
+			self.element_connectivity=self.element_connectivity(:,[1:nmax,ny]);
+
+
+		end % }}}
+		function mesh_connectivity2d(self,md,proj,varargin) % {{{
+			
+			%retrieve options: 
+			options=pairoptions(varargin{:}); 
+
+			subsetregions=getfieldvalue(options,'regions',1:self.nregions());
+
+			%initialize glacier connectivity: 
+			self.glacier_connectivity=zeros(self.nglaciers,1);
+
+			%assume maximum width for connectivity table and initialize 
+			%as sparse matrix: 
+			ny=self.nglaciers(); self.element_connectivity=sparse(md.mesh.numberofelements,ny);
+				
+			disp('Building element and glacier connectivity table: ');
+			[lat0,long0]=projlatlong(proj);
+		
+			[mpartition,npartition]=self.partition();
+			for i=subsetregions,
+				region=self.regions(i);
+				disp(sprintf(' progress for region: %s',region.name));
+
+				%project lat,long: 
+				[xlid,ylid]=gdaltransform(region.CenLon,region.CenLat,'EPSG:4326',proj);
+
+				%go through lids: 
+				x0=xlid; y0=ylid;
+				x1=md.mesh.x(md.mesh.elements(:,1)); y1=md.mesh.y(md.mesh.elements(:,1));
+				x2=md.mesh.x(md.mesh.elements(:,2)); y2=md.mesh.y(md.mesh.elements(:,2));
+				x3=md.mesh.x(md.mesh.elements(:,3)); y3=md.mesh.y(md.mesh.elements(:,3));
+				in=isintrianglearraytotal(x0,x1,x2,x3,y0,y1,y2,y3);
+				[els,glacs]=find(in);
+				self.glacier_connectivity(mpartition(i)-1+glacs)=els;
+			end
+
+			%build element connectivity table: 
+			for j=1:length(self.glacier_connectivity),
+				el=self.glacier_connectivity(j);
+				if ~el,continue; end;
 				count=self.element_connectivity(el,ny);
 				if count>ny,
@@ -332,4 +392,13 @@
 		end
 		%}}}
+		function listboxes(self) % {{{
+
+			for i=1:length(self.boxes),
+				b=self.boxes(i);
+				disp(sprintf('Region #%i: %s (%s)',i,b.FULL_NAME,b.RGI_CODE));
+			end
+			
+		end
+		%}}}
 	end
 end
Index: /issm/trunk-jpl/src/m/contrib/larour/isintrianglearray.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/isintrianglearray.m	(revision 24855)
+++ /issm/trunk-jpl/src/m/contrib/larour/isintrianglearray.m	(revision 24855)
@@ -0,0 +1,24 @@
+function isin=isintrianglearray(x0,x1,x2,x3,y0,y1,y2,y3)
+
+	isin=ones(length(x1),1);
+
+	dx = x0 - x3; dy = y0 - y3;
+	y23 = y2 - y3; x32 = x3 - x2; y31 = y3 - y1; x13 = x1 - x3; 
+	det = y23 .* x13 - x32 .* y31;
+	minD = min(det, 0); maxD = max(det, 0);
+
+
+	a = y23 .* dx + x32 .* dy;
+
+	pos=find( a < minD | a > maxD);
+	isin(pos)=0;
+					
+	b = y31 .* dx + x13 .* dy;
+	pos=find(b < minD | b > maxD);
+	isin(pos)=0;
+						
+	c = det - a - b;
+	pos=find(c < minD | c > maxD);
+	isin(pos)=0;
+
+end
Index: /issm/trunk-jpl/src/m/contrib/larour/isintrianglearraytotal.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/isintrianglearraytotal.m	(revision 24855)
+++ /issm/trunk-jpl/src/m/contrib/larour/isintrianglearraytotal.m	(revision 24855)
@@ -0,0 +1,28 @@
+function isin=isintrianglearraytotal(x0,x1,x2,x3,y0,y1,y2,y3)
+
+	isin=ones(length(x1),length(x0));
+
+	x0=repmat(x0',length(x1),1);
+	y0=repmat(y0',length(x1),1);
+
+	dx = x0 - x3; 
+	dy = y0 - y3;
+
+	y23 = y2 - y3; x32 = x3 - x2; y31 = y3 - y1; x13 = x1 - x3; 
+	det = y23 .* x13 - x32 .* y31;
+	minD = min(det, 0); maxD = max(det, 0);
+
+	a = y23 .* dx + x32 .* dy;
+
+	pos=find( a < minD | a > maxD);
+	isin(pos)=0;
+					
+	b = y31 .* dx + x13 .* dy;
+	pos=find(b < minD | b > maxD);
+	isin(pos)=0;
+						
+	c = det - a - b;
+	pos=find(c < minD | c > maxD);
+	isin(pos)=0;
+
+end