Index: /issm/trunk-jpl/src/m/mesh/triangle.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.m	(revision 27798)
+++ /issm/trunk-jpl/src/m/mesh/triangle.m	(revision 27799)
@@ -48,17 +48,18 @@
 removeorphans=1;
 if removeorphans,
-	orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
-	for i=1:length(orphan),
+	uniqueelements=sort(unique(elements(:)));
+	orphans=find(~ismember([1:length(x)],uniqueelements));
+	for i=1:length(orphans),
 		disp('WARNING: removing orphans');
 		%get rid of the orphan node i
 		%update x and y
-		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)];
+		x=[x(1:orphans(i)-(i-1)-1); x(orphans(i)-(i-1)+1:end)];
+		y=[y(1:orphans(i)-(i-1)-1); y(orphans(i)-(i-1)+1:end)];
 		%update elements
-		pos=find(elements>orphan(i)-(i-1));
+		pos=find(elements>orphans(i)-(i-1));
 		elements(pos)=elements(pos)-1;
 		%update segments
-		pos1=find(segments(:,1)>orphan(i)-(i-1));
-		pos2=find(segments(:,2)>orphan(i)-(i-1));
+		pos1=find(segments(:,1)>orphans(i)-(i-1));
+		pos2=find(segments(:,2)>orphans(i)-(i-1));
 		segments(pos1,1)=segments(pos1,1)-1;
 		segments(pos2,2)=segments(pos2,2)-1;
@@ -77,5 +78,6 @@
 md.mesh.numberofelements=size(md.mesh.elements,1);
 md.mesh.numberofvertices=length(md.mesh.x);
-md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 
 %Now, build the connectivity tables for this mesh.
Index: /issm/trunk-jpl/src/m/mesh/triangle.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 27798)
+++ /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 27799)
@@ -10,5 +10,5 @@
 
 def triangle(md, domainname, *args):
-    """TRIANGLE - create model mesh using the triangle package
+    """triangle - create model mesh using the triangle package
 
     This routine creates a model mesh using Triangle and a domain outline, to within a certain resolution
@@ -27,6 +27,6 @@
     """
 
-    #Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m  resolution node would
-    #be made of 1000 * 1000 area squares).
+    # Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m  resolution node would
+    # be made of 1000 * 1000 area squares).
 
     if len(args) == 1:
@@ -37,9 +37,9 @@
         resolution = args[1]
 
-    #Check that mesh was not already run, and warn user:
+    # Check that mesh was not already run, and warn user:
     if md.mesh.numberofelements:
         choice = input('This model already has a mesh. Are you sure you want to go ahead? (y / n)')
         if choice not in ['y', 'n']:
-            print("bad answer try you should use 'y' or 'n' ... exiting")
+            print('bad answer try you should use \'y\' or \'n\' ... exiting')
             return None
         if choice == 'n':
@@ -49,16 +49,42 @@
     area = resolution ** 2
 
-    #Check that file exist (this is a very very common mistake)
+    # Check that file exist (this is a very very common mistake)
     if not os.path.exists(domainname):
-        raise IOError("file '%s' not found" % domainname)
+        raise IOError('file {} not found'.format(domainname))
 
-    #Mesh using Triangle
+    # Mesh using Triangle
+    elements, x, y, segments, segmentmarkers = Triangle_python(domainname, riftname, area)
+
+    # Check that all the created nodes belong to at least one element
+    removeorphans = 1
+    if removeorphans:
+        uniqueelements = np.sort(np.unique(elements))
+        orphans = np.nonzero((~np.isin(range(1, len(x)), uniqueelements)).astype(int))[0]
+        for i in range(0, len(orphans)):
+            print('WARNING: removing orphans')
+            # Get rid of the orphan node i
+            # Update x and y
+            x = np.concatenate((x[0:(orphans[i] - i)], x[(orphans[i] - i + 1):]))
+            y = np.concatenate((y[0:(orphans[i] - i)], y[(orphans[i] - i + 1):]))
+            # Update elements
+            pos = np.nonzero((elements > (orphans[i] - i)).flatten(order='F'))[0]
+            elementstmp = elements.flatten(order='F')
+            elementstmp[pos] -= 1
+            elements = elementstmp.reshape(np.shape(elements), order='F')
+            # Update segments
+            pos1 = np.nonzero(segments[:,0] > (orphans[i] - i))[0]
+            pos2 = np.nonzero(segments[:,1] > (orphans[i] - i))[0]
+            segments[pos1, 0] -= 1
+            segments[pos2, 1] -= 1
+
+    # Plug into md
     md.mesh = mesh2d()
-    md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.segments, md.mesh.segmentmarkers = Triangle_python(domainname, riftname, area)
-    md.mesh.elements = md.mesh.elements.astype(int)
-    md.mesh.segments = md.mesh.segments.astype(int)
-    md.mesh.segmentmarkers = md.mesh.segmentmarkers.astype(int)
+    md.mesh.x = x
+    md.mesh.y = y
+    md.mesh.elements = elements.astype(int)
+    md.mesh.segments = segments.astype(int)
+    md.mesh.segmentmarkers = segmentmarkers.astype(int)
 
-    #Fill in rest of fields:
+    # Fill in rest of fields:
     md.mesh.numberofelements = np.size(md.mesh.elements, axis=0)
     md.mesh.numberofvertices = np.size(md.mesh.x)
@@ -66,5 +92,5 @@
     md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
 
-    #Now, build the connectivity tables for this mesh.
+    # Now, build the connectivity tables for this mesh.
     md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
     md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)