Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 13565)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 13566)
@@ -137,5 +137,5 @@
 			nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1
 			#Macayeal elements in contact with this layer become MacAyealPattyn elements
-			matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealpattyn))
+			matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealpattyn)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
 			commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
@@ -146,12 +146,12 @@
 
 			#rule out elements that don't touch the 2 boundaries
-			pos=numpy.nonzero(macayealpattynflag)
-			elist=numpy.zeros(len(pos))
-			elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1)
-			elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1),axis=1)
-			pos1=numpy.nonzero(elist==1)
+			pos=numpy.nonzero(macayealpattynflag)[0]
+			elist=numpy.zeros(numpy.size(pos),dtype=int)
+			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1).astype(bool)
+			pos1=numpy.nonzero(elist==1)[0]
 			macayealflag[pos[pos1]]=1
 			macayealpattynflag[pos[pos1]]=0
-			pos2=numpy.nonzero(elist==-1)
+			pos2=numpy.nonzero(elist==-1)[0]
 			pattynflag[pos[pos2]]=1
 			macayealpattynflag[pos[pos2]]=0
@@ -169,5 +169,5 @@
 			nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1
 			#Stokes elements in contact with this layer become PattynStokes elements
-			matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonpattynstokes))
+			matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonpattynstokes)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
 			commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
@@ -178,12 +178,12 @@
 
 			#rule out elements that don't touch the 2 boundaries
-			pos=numpy.nonzero(pattynstokesflag)
-			elist=numpy.zeros(len(pos))
-			elist = elist + numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1)
-			elist = elist - numpy.any(numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1)
-			pos1=numpy.nonzero(elist==1)
+			pos=numpy.nonzero(pattynstokesflag)[0]
+			elist=numpy.zeros(numpy.size(pos),dtype=int)
+			elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
+			pos1=numpy.nonzero(elist==1)[0]
 			stokesflag[pos[pos1]]=1
 			pattynstokesflag[pos[pos1]]=0
-			pos2=numpy.nonzero(elist==-1)
+			pos2=numpy.nonzero(elist==-1)[0]
 			pattynflag[pos[pos2]]=1
 			pattynstokesflag[pos[pos2]]=0
@@ -201,5 +201,5 @@
 			nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1
 			#Stokes elements in contact with this layer become MacAyealStokes elements
-			matrixelements=ismember(md.mesh.elements,numpy.nonzero(nodeonmacayealstokes))
+			matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealstokes)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
 			commonelements[numpy.nonzero(macayealflag)]=0    #only one layer: the elements previously in macayeal
@@ -210,12 +210,12 @@
 
 			#rule out elements that don't touch the 2 boundaries
-			pos=numpy.nonzero(macayealstokesflag)
-			elist=numpy.zeros(len(pos))
-			elist = elist + numpy.any(numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1),axis=1)
-			elist = elist - numpy.any(numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1),axis=1)
-			pos1=numpy.nonzero(elist==1)
+			pos=numpy.nonzero(macayealstokesflag)[0]
+			elist=numpy.zeros(numpy.size(pos),dtype=int)
+			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1).astype(bool)
+			pos1=numpy.nonzero(elist==1)[0]
 			macayealflag[pos[pos1]]=1
 			macayealstokesflag[pos[pos1]]=0
-			pos2=numpy.nonzero(elist==-1)
+			pos2=numpy.nonzero(elist==-1)[0]
 			stokesflag[pos[pos2]]=1
 			macayealstokesflag[pos[pos2]]=0