Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 25619)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 25620)
@@ -168,4 +168,31 @@
 			self.basins{end+1}=bas;
 		end % }}}
+		function intersections2d(self,varargin) % {{{
+
+			options=pairoptions(varargin{:}); 
+			force=getfieldvalue(options,'force',0);
+			
+			%initialize, to avoid issues of having more transitions than meshes. 
+			self.transitions={};
+			self.eltransitions={};
+
+			%for elements:
+			xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
+			ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
+			
+			for i=1:length(self.icecaps),
+				mdi=self.icecaps{i};
+		
+				%for elements: 
+				xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
+				yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
+		
+				disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+			
+				self.transitions{end+1}=meshintersect2d(self.earth.mesh.x,self.earth.mesh.y,mdi.mesh.x,mdi.mesh.y,'force',force);
+
+				self.eltransitions{end+1}=meshintersect2d(xe,ye,xei,yei,'force',force);
+			end
+		end % }}}
 		function intersections(self,varargin) % {{{
 
@@ -360,14 +387,43 @@
 
 		end % }}}
-		function viscousiterations(self) % {{{
-			for  i=1:length(self.icecaps),
-				ic=self.icecaps{i}; 
-				mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
-				for j=2:length(ic.results.TransientSolution)-1,
-					mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
-				end
-				disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
-			end
-		end % }}}
+		function caticecaps2d(self,varargin) % {{{
+
+			%recover options: 
+			options=pairoptions(varargin{:}); 
+			tolerance=getfieldvalue(options,'tolerance',1e-5);
+			loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
+			models=self.icecaps;
+
+			%Plug all models together:
+			md=models{1}; 
+			for i=2:length(models),
+				md=modelmerge2d(md,models{i},'tolerance',tolerance);
+			end
+
+			%Look for lone edges if asked for it: {{{
+			if loneedgesdetect,
+				edges=loneedges(md);
+				hold on;
+				for i=1:length(edges),
+					ind1=edges(i,1);
+					ind2=edges(i,2);
+					plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'g*-');
+				end
+			end %}}}
+
+			%Plug into earth: 
+			self.earth=md;
+
+		end % }}}
+function viscousiterations(self) % {{{
+	for  i=1:length(self.icecaps),
+		ic=self.icecaps{i}; 
+		mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
+		for j=2:length(ic.results.TransientSolution)-1,
+			mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
+		end
+		disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
+	end
+end % }}}
 		function maxtimestep(self) % {{{
 			for  i=1:length(self.icecaps),