Index: /issm/trunk-jpl/src/m/consistency/checkfield.py
===================================================================
--- /issm/trunk-jpl/src/m/consistency/checkfield.py	(revision 22864)
+++ /issm/trunk-jpl/src/m/consistency/checkfield.py	(revision 22865)
@@ -9,6 +9,6 @@
 
 	   Used to check model consistency.,
-	   Requires: 
-	   'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname)) 
+	   Requires:
+	   'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname))
              If 'field' is provided, it will assume the argument following 'field' is a numeric array.
 
@@ -40,5 +40,5 @@
 		fieldname=options.getfieldvalue('fieldname','no fieldname')
 	else:
-		fieldname=options.getfieldvalue('fieldname') 
+		fieldname=options.getfieldvalue('fieldname')
 		exec("field=md.{}".format(fieldname))
 
@@ -80,5 +80,5 @@
 				if (not np.size(field,0)==fieldsize[0]) or (not np.size(field,1)==fieldsize[1]):
 					md = md.checkmessage(options.getfieldvalue('message',"field '%s' size should be %d x %d" % (fieldname,fieldsize[0],fieldsize[1])))
-	
+
 	#Check numel
 	if options.exist('numel'):
@@ -164,5 +164,5 @@
 
 		if options.getfieldvalue('singletimeseries',0):
-			field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1)
+			field2 = np.reshape(field[0],int(np.prod(np.shape(field[0]))),1)
 
 		if np.any(field2>upperbound):
@@ -176,5 +176,5 @@
 
 		if options.getfieldvalue('singletimeseries',0):
-			field2 = np.reshape(field[0],np.prod(np.shape(field[0])),1)
+			field2 = np.reshape(field[0],int(np.prod(np.shape(field[0]))),1)
 
 		if np.any(field2>=upperbound):
@@ -228,4 +228,2 @@
 
 	return md
-
-
Index: /issm/trunk-jpl/src/m/mesh/ComputeHessian.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/ComputeHessian.py	(revision 22864)
+++ /issm/trunk-jpl/src/m/mesh/ComputeHessian.py	(revision 22865)
@@ -33,34 +33,33 @@
 	linesize=3*numberofelements
 
-	#get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma 
+	#get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma
 	[alpha,beta,dum]=GetNodalFunctionsCoeff(index,x,y)
 	areas=GetAreas(index,x,y)
 
 	#compute weights that hold the volume of all the element holding the node i
-	weights=m.sparse(line,np.ones((linesize,1)),np.tile(areas.reshape(-1,),(3,1)),numberofnodes,1)
+	weights=m.sparse(line,np.ones((linesize,1),dtype=int),np.tile(areas,(3,1)),numberofnodes,1)
 
 	#compute field on nodes if on elements
 	if np.size(field,axis=0)==numberofelements:
-		field=m.sparse(line,np.ones((linesize,1)),np.tile(areas*field,(3,1)),numberofnodes,1)/weights
+		field=m.sparse(line,np.ones((linesize,1),dtype=int),np.tile(areas*field,(3,1)),numberofnodes,1)/weights
 
 	#Compute gradient for each element
-	grad_elx=np.sum(field[index-1,0]*alpha,axis=1) 
-	grad_ely=np.sum(field[index-1,0]*beta,axis=1)
+	grad_elx=np.sum(field[index-1]*alpha,axis=1)
+	grad_ely=np.sum(field[index-1]*beta,axis=1)
 
 	#Compute gradient for each node (average of the elements around)
-	gradx=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_elx).reshape(-1,),(3,1)),numberofnodes,1)
-	grady=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_ely).reshape(-1,),(3,1)),numberofnodes,1)
+	gradx=m.sparse(line,np.ones((linesize,1),dtype=int),np.tile((areas*grad_elx),(3,1)),numberofnodes,1)
+	grady=m.sparse(line,np.ones((linesize,1),dtype=int),np.tile((areas*grad_ely),(3,1)),numberofnodes,1)
 	gradx=gradx/weights
 	grady=grady/weights
 
 	#Compute hessian for each element
-	hessian=np.vstack((np.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,),np.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,),np.sum(grady[index-1,0]*beta,axis=1).reshape(-1,))).T
+	hessian=np.vstack((np.sum(gradx[index-1,0]*alpha,axis=1),np.sum(grady[index-1,0]*alpha,axis=1),np.sum(grady[index-1,0]*beta,axis=1))).T
 
 	if m.strcmpi(type,'node'):
 		#Compute Hessian on the nodes (average of the elements around)
-		hessian=np.hstack((m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,0]).reshape(-1,),(3,1)),numberofnodes,1)/weights,
-											 m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,1]).reshape(-1,),(3,1)),numberofnodes,1)/weights,
-											 m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,2]).reshape(-1,),(3,1)),numberofnodes,1)/weights ))
+		hessian=np.hstack((m.sparse(line,np.ones((linesize,1),dtype=int),np.tile((areas*hessian[:,0]),(3,1)),numberofnodes,1)/weights,
+											 m.sparse(line,np.ones((linesize,1),dtype=int),np.tile((areas*hessian[:,1]),(3,1)),numberofnodes,1)/weights,
+											 m.sparse(line,np.ones((linesize,1),dtype=int),np.tile((areas*hessian[:,2]),(3,1)),numberofnodes,1)/weights ))
 
 	return hessian
-
Index: /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py	(revision 22864)
+++ /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py	(revision 22865)
@@ -105,3 +105,2 @@
 
 	return y
-
