Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.py	(revision 25678)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py	(revision 25679)
@@ -147,4 +147,8 @@
     #}}}
 
+    def ncaps(self): #{{{
+        return len(self.icecaps)
+    #}}}
+
     def continents(self):  # {{{
         list = []
@@ -166,4 +170,29 @@
             raise RuntimeError('addbasin method only takes a \'basin\' class object as input')
         self.basins.append(bas)
+    #}}}
+
+    def intersections2d(self, *args): #{{{
+        options = pairoptions(*args)
+        force = options.getfieldvalue('force', 0)
+
+        # Initialize, to avoid issues of having more transitions than meshes
+        self.transitions = []
+        self.eltransitions = []
+
+        # For elements
+        xe = np.mean(self.earth.mesh.x[self.earth.mesh.elements - 1], axis=1)
+        ye = np.mean(self.earth.mesh.y[self.earth.mesh.elements - 1], axis=1)
+
+        for i in range(len(self.icecaps)):
+            mdi = self.icecaps[i]
+
+            # For elements
+            xei = np.mean(mdi.mesh.x[mdi.mesh.elements - 1], axis=1)
+            yei = np.mean(mdi.mesh.y[mdi.mesh.elements - 1], axis=1)
+
+            print('Computing vertex intersections for basin {}'.format(self.basins[i].name))
+
+            self.transitions.append(meshintersect2d(self.earth.mesh.x, self.earth.mesh.y, mdi.mesh.x, mdi.mesh.y, 'force', force))
+            self.eltransitions.append(meshintersect2d(xe, ye, xei, yei, 'force', force))
     #}}}
 
@@ -315,9 +344,9 @@
         if loneedgesdetect:
             edges = loneedges(md)
+            # TODO: Reconfigure the following in the process of bringing plotting online
             plotmodel(md, 'data', md.mask.land_levelset)
             for i in range(len(edges)):
                 ind1 = edges(i, 1)
                 ind2 = edges(i, 2)
-                # TODO: Reconfigure the following in the process of bringing plotting online
                 plot3([md.mesh.x[ind1], md.mesh.x[ind2]], [md.mesh.y[ind1], md.mesh.y[ind2]], [md.mesh.z[ind1], md.mesh.z[ind2]], 'g*-')
 
@@ -327,4 +356,30 @@
         # Create mesh radius
         self.earth.mesh.r = planetradius('earth') * np.ones((md.mesh.numberofvertices, ))
+    #}}}
+
+    def caticecaps2d(self, *args):  #{{{
+        # Recover options
+        options = pairoptions(*args)
+        tolerance = options.getfieldvalue('tolerance', 1e-5)
+        loneedgesdetect = options.getfieldvalue('loneedgesdetect', 0)
+        models = self.icecaps
+
+        # Plug all models together
+        md = models[0]
+        for i in range(1, len(models)):
+            md = modelmerge2d(md, models[i], 'tolerance', tolerance)
+
+        # Look for lone edges if asked for it
+        if loneedgesdetect:
+            edges = loneedges(md)
+            # TODO: Reconfigure the following in the process of bringing plotting online
+            plotmodel(md, 'data', md.mask.land_levelset)
+            for i in range(len(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*-')
+
+        # Plug into earth
+        self.earth = md
     #}}}