Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 17562)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 17563)
@@ -421,28 +421,30 @@
 
 			%Edges
-			if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
-				%renumber first two columns
-				pos=find(md2.mesh.edges(:,4)~=-1);
-				md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1));
-				md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2));
-				md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
-				md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
-				%remove edges when the 2 vertices are not in the domain.
-				md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
-				%Replace all zeros by -1 in the last two columns
-				pos=find(md2.mesh.edges(:,3)==0);
-				md2.mesh.edges(pos,3)=-1;
-				pos=find(md2.mesh.edges(:,4)==0);
-				md2.mesh.edges(pos,4)=-1;
-				%Invert -1 on the third column with last column (Also invert first two columns!!)
-				pos=find(md2.mesh.edges(:,3)==-1);
-				md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
-				md2.mesh.edges(pos,4)=-1;
-				values=md2.mesh.edges(pos,2);
-				md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
-				md2.mesh.edges(pos,1)=values;
-				%Finally remove edges that do not belong to any element
-				pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
-				md2.mesh.edges(pos,:)=[];
+			if(strcmp(meshtype(md.mesh),'2Dhorizontal')),
+				if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
+					%renumber first two columns
+					pos=find(md2.mesh.edges(:,4)~=-1);
+					md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1));
+					md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2));
+					md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
+					md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
+					%remove edges when the 2 vertices are not in the domain.
+					md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
+					%Replace all zeros by -1 in the last two columns
+					pos=find(md2.mesh.edges(:,3)==0);
+					md2.mesh.edges(pos,3)=-1;
+					pos=find(md2.mesh.edges(:,4)==0);
+					md2.mesh.edges(pos,4)=-1;
+					%Invert -1 on the third column with last column (Also invert first two columns!!)
+					pos=find(md2.mesh.edges(:,3)==-1);
+					md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
+					md2.mesh.edges(pos,4)=-1;
+					values=md2.mesh.edges(pos,2);
+					md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
+					md2.mesh.edges(pos,1)=values;
+					%Finally remove edges that do not belong to any element
+					pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
+					md2.mesh.edges(pos,:)=[];
+				end
 			end
 
Index: /issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.py	(revision 17562)
+++ /issm/trunk-jpl/src/m/classes/model.py	(revision 17563)
@@ -327,28 +327,29 @@
 
 		#Edges
-		if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some NaNs...
-			#renumber first two columns
-			pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
-			md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0]-1]
-			md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1]-1]
-			md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2]-1]
-			md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
-			#remove edges when the 2 vertices are not in the domain.
-			md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
-			#Replace all zeros by -1 in the last two columns
-			pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0]
-			md2.mesh.edges[pos,2]=-1
-			pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0]
-			md2.mesh.edges[pos,3]=-1
-			#Invert -1 on the third column with last column (Also invert first two columns!!)
-			pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0]
-			md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3]
-			md2.mesh.edges[pos,3]=-1
-			values=md2.mesh.edges[pos,1]
-			md2.mesh.edges[pos,1]=md2.mesh.edges[pos,0]
-			md2.mesh.edges[pos,0]=values
-			#Finally remove edges that do not belong to any element
-			pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]
-			md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0)
+		if m.strcmp(md.mesh.meshtype(),'2Dhorizontal'):
+			if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some NaNs...
+				#renumber first two columns
+				pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
+				md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0]-1]
+				md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1]-1]
+				md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2]-1]
+				md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
+				#remove edges when the 2 vertices are not in the domain.
+				md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
+				#Replace all zeros by -1 in the last two columns
+				pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0]
+				md2.mesh.edges[pos,2]=-1
+				pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0]
+				md2.mesh.edges[pos,3]=-1
+				#Invert -1 on the third column with last column (Also invert first two columns!!)
+				pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0]
+				md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3]
+				md2.mesh.edges[pos,3]=-1
+				values=md2.mesh.edges[pos,1]
+				md2.mesh.edges[pos,1]=md2.mesh.edges[pos,0]
+				md2.mesh.edges[pos,0]=values
+				#Finally remove edges that do not belong to any element
+				pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]
+				md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0)
 
 		#Penalties
Index: /issm/trunk-jpl/src/m/classes/rifts.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/rifts.py	(revision 17562)
+++ /issm/trunk-jpl/src/m/classes/rifts.py	(revision 17563)
@@ -5,4 +5,5 @@
 from WriteData import WriteData
 from isnans import isnans
+import MatlabFuncs as m
 
 class rifts(object):
Index: /issm/trunk-jpl/src/m/classes/thermal.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/thermal.py	(revision 17562)
+++ /issm/trunk-jpl/src/m/classes/thermal.py	(revision 17563)
@@ -4,4 +4,5 @@
 from checkfield import checkfield
 from WriteData import WriteData
+import MatlabFuncs as m
 
 class thermal(object):