Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 13983) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 13984) @@ -25,8 +25,8 @@ if not os.path.exists(icefrontfile): raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile) - [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) - nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float) + [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) + nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) else: - nodeonicefront=numpy.zeros((md.mesh.numberofvertices)) + nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool) # pos=find(md.mesh.vertexonboundary & ~nodeonicefront); @@ -56,5 +56,5 @@ #segment on Neumann (Ice Front) # pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); - pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0] + pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0]-1],nodeonicefront[md.mesh.segments[:,1]-1]))[0] if md.mesh.dimension==2: pressureload=md.mesh.segments[pos,:] @@ -62,5 +62,5 @@ # pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2])) - pressureload=numpy.zeros((0,5)) + pressureload=numpy.zeros((0,5),int) for i in xrange(1,md.mesh.numberoflayers): # pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; @@ -69,5 +69,5 @@ #Add water or air enum depending on the element # pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))]; - pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1))) + pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1))) #plug onto model Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 13983) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 13984) @@ -27,11 +27,11 @@ if not os.path.exists(icefrontfile): raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile) - [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) - vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float) + [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2) + vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)) else: #Guess where the ice front is - vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices)) - vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:].astype(int)-1]=1 - vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice).astype(float) + vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices),bool) + vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:]-1]=True + vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice) # pos=find(md.mesh.vertexonboundary & ~vertexonicefront); @@ -62,5 +62,5 @@ #segment on Neumann (Ice Front) # pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); - pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0].astype(int)-1],vertexonicefront[md.mesh.segments[:,1].astype(int)-1]))[0] + pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0]-1],vertexonicefront[md.mesh.segments[:,1]-1]))[0] if md.mesh.dimension==2: pressureload=md.mesh.segments[pos,:] @@ -68,5 +68,5 @@ # pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2])) - pressureload=numpy.zeros((0,5)) + pressureload=numpy.zeros((0,5),int) for i in xrange(1,md.mesh.numberoflayers): # pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; @@ -75,5 +75,5 @@ #Add water or air enum depending on the element # pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))]; - pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)+0.*md.mask.elementongroundedice[pressureload[:,-1].astype('int')-1].reshape(-1,1))) + pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)+0*md.mask.elementongroundedice[pressureload[:,-1]-1].reshape(-1,1))) #plug onto model Index: /issm/trunk-jpl/src/m/classes/dependent.py =================================================================== --- /issm/trunk-jpl/src/m/classes/dependent.py (revision 13983) +++ /issm/trunk-jpl/src/m/classes/dependent.py (revision 13984) @@ -45,5 +45,5 @@ #process the file and retrieve segments mesh=options.getfieldvalue('mesh') - self.segments=MeshProfileIntersection(mesh.elements,mesh.x,mesh.y,self.exp) + self.segments=MeshProfileIntersection(mesh.elements.astype(float),mesh.x,mesh.y,self.exp) # }}} Index: /issm/trunk-jpl/src/m/classes/initialization.py =================================================================== --- /issm/trunk-jpl/src/m/classes/initialization.py (revision 13983) +++ /issm/trunk-jpl/src/m/classes/initialization.py (revision 13984) @@ -70,6 +70,6 @@ md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices]) #Triangle with zero velocity - if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements.astype(int)-1]),axis=1)==0,\ - numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements.astype(int)-1]),axis=1)==0)): + if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\ + numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)): md.checkmessage("at least one triangle has all its vertices with a zero velocity") if ThermalAnalysisEnum() in analyses: Index: /issm/trunk-jpl/src/m/classes/mask.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.py (revision 13983) +++ /issm/trunk-jpl/src/m/classes/mask.py (revision 13984) @@ -33,5 +33,5 @@ string="%s\n%s"%(string,fielddisplay(self,"elementonfloatingice","element on floating ice flags list")) string="%s\n%s"%(string,fielddisplay(self,"vertexonfloatingice","vertex on floating ice flags list")) - string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list")) + string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list")) string="%s\n%s"%(string,fielddisplay(self,"vertexongroundedice","vertex on grounded ice flags list")) string="%s\n%s"%(string,fielddisplay(self,"elementonwater","element on water flags list")) Index: /issm/trunk-jpl/src/m/classes/mesh.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh.py (revision 13983) +++ /issm/trunk-jpl/src/m/classes/mesh.py (revision 13984) @@ -95,5 +95,5 @@ string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (NaN for element on the upper layer)")) string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (NaN for vertex on the lower surface)")) - string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (NaN for element on the lower layer")) + string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (NaN for element on the lower layer)")) string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list")) string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)")) Index: /issm/trunk-jpl/src/m/classes/model/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.py (revision 13983) +++ /issm/trunk-jpl/src/m/classes/model/model.py (revision 13984) @@ -207,13 +207,13 @@ #kick out all elements with 3 dirichlets spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] - spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1 + spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1 flag=numpy.ones(md1.mesh.numberofvertices) flag[spc_node]=0 - pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements.astype(int)-1],axis=1)))[0] + pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0] flag_elem[pos]=0 #extracted elements and nodes lists pos_elem=numpy.nonzero(flag_elem)[0] - pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1 + pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1 #keep track of some fields @@ -234,11 +234,11 @@ elements_1=copy.deepcopy(md1.mesh.elements) elements_2=elements_1[pos_elem,:] - elements_2[:,0]=Pnode[elements_2[:,0].astype(int)-1] - elements_2[:,1]=Pnode[elements_2[:,1].astype(int)-1] - elements_2[:,2]=Pnode[elements_2[:,2].astype(int)-1] + elements_2[:,0]=Pnode[elements_2[:,0]-1] + elements_2[:,1]=Pnode[elements_2[:,1]-1] + elements_2[:,2]=Pnode[elements_2[:,2]-1] if md1.mesh.dimension==3: - elements_2[:,3]=Pnode[elements_2[:,3].astype(int)-1] - elements_2[:,4]=Pnode[elements_2[:,4].astype(int)-1] - elements_2[:,5]=Pnode[elements_2[:,5].astype(int)-1] + elements_2[:,3]=Pnode[elements_2[:,3]-1] + elements_2[:,4]=Pnode[elements_2[:,4]-1] + elements_2[:,5]=Pnode[elements_2[:,5]-1] #OK, now create the new model! @@ -316,7 +316,7 @@ md2.mesh.numberofvertices2d=numpy.size(pos_node_2d) md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:] - md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0].astype(int)-1] - md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1].astype(int)-1] - md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2].astype(int)-1] + md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1] + md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1]-1] + md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2]-1] md2.mesh.x2d=md1.mesh.x[pos_node_2d] @@ -324,11 +324,11 @@ #Edges - if len(numpy.shape(md2.mesh.edges))>1 and numpy.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some NaNs... + 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].astype(int)-1] - md2.mesh.edges[: ,1]=Pnode[md2.mesh.edges[:,1].astype(int)-1] - md2.mesh.edges[: ,2]=Pelem[md2.mesh.edges[:,2].astype(int)-1] - md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3].astype(int)-1] + 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],:] @@ -361,25 +361,31 @@ #recreate segments if md1.mesh.dimension==2: - [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices) - [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity) + [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements.astype(float),md2.mesh.numberofvertices) + md2.mesh.vertexconnectivity=md2.mesh.vertexconnectivity.astype(int) + [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements.astype(float),md2.mesh.vertexconnectivity.astype(float)) + md2.mesh.elementconnectivity=md2.mesh.elementconnectivity.astype(int) md2.mesh.segments=contourenvelope(md2) - md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2) - md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1 + md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool) + md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True else: #First do the connectivity for the contourenvelope in 2d - [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d) - [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity) + [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d.astype(float),md2.mesh.numberofvertices2d) + md2.mesh.vertexconnectivity=md2.mesh.vertexconnectivity.astype(int) + [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d.astype(float),md2.mesh.vertexconnectivity.astype(float)) + md2.mesh.elementconnectivity=md2.mesh.elementconnectivity.astype(int) md2.mesh.segments=contourenvelope(md2) - md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers) - md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1 + md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool) + md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers) #Then do it for 3d as usual - [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices) - [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity) + [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements.astype(float),md2.mesh.numberofvertices) + md2.mesh.vertexconnectivity=md2.mesh.vertexconnectivity.astype(int) + [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements.astype(float),md2.mesh.vertexconnectivity.astype(float)) + md2.mesh.elementconnectivity=md2.mesh.elementconnectivity.astype(int) #Boundary conditions: Dirichlets on new boundary #Catch the elements that have not been extracted orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0] - orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1 + orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1 #Figure out which node are on the boundary between md2 and md1 nodestoflag1=numpy.intersect1d(orphans_node,pos_node) @@ -443,6 +449,6 @@ #Keep track of pos_node and pos_elem - md2.mesh.extractedvertices=pos_node.astype(float)+1 - md2.mesh.extractedelements=pos_elem.astype(float)+1 + md2.mesh.extractedvertices=pos_node+1 + md2.mesh.extractedelements=pos_elem+1 return md2 @@ -526,5 +532,5 @@ #Extrude elements - elements3d=numpy.empty((0,6)) + elements3d=numpy.empty((0,6),int) for i in xrange(numlayers-1): elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices)))) #Create the elements of the 3d mesh for the non extruded part @@ -603,8 +609,8 @@ #bedinfo and surface info - md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',1) - md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',md.mesh.numberoflayers-1) - md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',1) - md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',md.mesh.numberoflayers) + md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',1) + md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',md.mesh.numberoflayers-1) + md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1) + md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers) #elementstype @@ -635,5 +641,5 @@ #in 3d, pressureload: [node1 node2 node3 node4 element] pressureload_layer1=numpy.hstack((md.diagnostic.icefront[:,0:2],md.diagnostic.icefront[:,1:2]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,0:1]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,2:4])) #Add two columns on the first layer - pressureload=numpy.empty((0,6)) + pressureload=numpy.empty((0,6),int) for i in xrange(numlayers-1): pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+i*md.mesh.numberofvertices2d,pressureload_layer1[:,4:5]+i*md.mesh.numberofelements2d,pressureload_layer1[:,5:6])))) @@ -642,4 +648,5 @@ #connectivity md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1)) + md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(float) md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=float('NaN') for i in xrange(1,numlayers-1): @@ -647,4 +654,5 @@ =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0 + md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) #materials Index: /issm/trunk-jpl/src/m/consistency/checkfield.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 13983) +++ /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 13984) @@ -143,5 +143,5 @@ if options.getfieldvalue('forcing',0): if numpy.size(field,0)==md.mesh.numberofvertices: - if len(numpy.shape(field))>1 and not numpy.size(field,1)==1: + if numpy.ndim(field)>1 and not numpy.size(field,1)==1: md = md.checkmessage(options.getfieldvalue('message',\ "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname)) Index: /issm/trunk-jpl/src/m/extrusion/project3d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 13983) +++ /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 13984) @@ -38,5 +38,5 @@ vector1d=False - if isinstance(vector2d,numpy.ndarray) and len(vector2d.shape)==1: + if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1: vector1d=True vector2d=vector2d.reshape(numpy.size(vector2d),1) @@ -49,7 +49,7 @@ #Initialize 3d vector if numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d: - projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1))) + projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) elif numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d+1: - projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1))) + projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) projected_vector[-1,:]=vector2d[-1,:] vector2d=vector2d[:-1,:] @@ -68,7 +68,7 @@ #Initialize 3d vector if numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d: - projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1))) + projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) elif numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d+1: - projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1))) + projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) projected_vector[-1,:]=vector2d[-1,:] vector2d=vector2d[:-1,:] Index: /issm/trunk-jpl/src/m/geometry/FlagElements.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/FlagElements.py (revision 13983) +++ /issm/trunk-jpl/src/m/geometry/FlagElements.py (revision 13984) @@ -24,8 +24,8 @@ if isinstance(region,(str,unicode)): if not region: - flag=numpy.zeros(md.mesh.numberofelements,'bool') + flag=numpy.zeros(md.mesh.numberofelements,bool) invert=0 elif strcmpi(region,'all'): - flag=numpy.ones(md.mesh.numberofelements,'bool') + flag=numpy.ones(md.mesh.numberofelements,bool) invert=0 else: @@ -44,8 +44,9 @@ xlim,ylim=basinzoom('basin',region) flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.xxlim[0]),numpy.logical_and(md.mesh.yylim[0])).astype(float) - flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1) + flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool) else: #ok, flag elements - [flag,dum]=ContourToMesh(md.mesh.elements[:,0:3].copy(),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1) + [flag,dum]=ContourToMesh(md.mesh.elements[:,0:3].astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1) + flag=flag.astype(bool) if invert: Index: /issm/trunk-jpl/src/m/geometry/GetAreas.py =================================================================== --- /issm/trunk-jpl/src/m/geometry/GetAreas.py (revision 13983) +++ /issm/trunk-jpl/src/m/geometry/GetAreas.py (revision 13984) @@ -33,10 +33,10 @@ #initialization areas=numpy.zeros(nels) - x1=x[index[:,0].astype(int)-1] - x2=x[index[:,1].astype(int)-1] - x3=x[index[:,2].astype(int)-1] - y1=y[index[:,0].astype(int)-1] - y2=y[index[:,1].astype(int)-1] - y3=y[index[:,2].astype(int)-1] + x1=x[index[:,0]-1] + x2=x[index[:,1]-1] + x3=x[index[:,2]-1] + y1=y[index[:,0]-1] + y2=y[index[:,1]-1] + y3=y[index[:,2]-1] #compute the volume of each element @@ -46,5 +46,5 @@ else: #V=area(triangle)*1/3(z1+z2+z3) - thickness=numpy.mean(z[index[:,3:6].astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1]) + thickness=numpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1]) areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness Index: /issm/trunk-jpl/src/m/mesh/ComputeHessian.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/ComputeHessian.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/ComputeHessian.py (revision 13984) @@ -45,6 +45,6 @@ #Compute gradient for each element - grad_elx=numpy.sum(field[index.astype(int)-1,0]*alpha,axis=1) - grad_ely=numpy.sum(field[index.astype(int)-1,0]*beta,axis=1) + grad_elx=numpy.sum(field[index-1,0]*alpha,axis=1) + grad_ely=numpy.sum(field[index-1,0]*beta,axis=1) #Compute gradient for each node (average of the elements around) @@ -55,5 +55,5 @@ #Compute hessian for each element - hessian=numpy.hstack((numpy.sum(gradx[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*beta,axis=1).reshape(-1,1))) + hessian=numpy.hstack((numpy.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*beta,axis=1).reshape(-1,1))) if strcmpi(type,'node'): Index: /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py (revision 13984) @@ -39,10 +39,10 @@ #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma - x1=x[index[:,0].astype(int)-1] - x2=x[index[:,1].astype(int)-1] - x3=x[index[:,2].astype(int)-1] - y1=y[index[:,0].astype(int)-1] - y2=y[index[:,1].astype(int)-1] - y3=y[index[:,2].astype(int)-1] + x1=x[index[:,0]-1] + x2=x[index[:,1]-1] + x3=x[index[:,2]-1] + y1=y[index[:,0]-1] + y2=y[index[:,1]-1] + y3=y[index[:,2]-1] invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2)) Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13984) @@ -276,5 +276,5 @@ else: bamg_mesh.Vertices=numpy.hstack((md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),numpy.ones((md.mesh.numberofvertices,1)))) - bamg_mesh.Triangles=numpy.hstack((md.mesh.elements,numpy.ones((md.mesh.numberofelements,1)))) + bamg_mesh.Triangles=numpy.hstack((md.mesh.elements.astype(float),numpy.ones((md.mesh.numberofelements,1)))) if isinstance(md.rifts.riftstruct,dict): @@ -320,8 +320,8 @@ md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() - md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy() - md.mesh.edges=bamgmesh_out['IssmEdges'].copy() - md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy() - md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy() + md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) + md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) + md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) + md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) #Fill in rest of fields: @@ -331,13 +331,14 @@ md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0) md.mesh.z=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) - md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) - md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) - md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) - md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) + md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) + md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0 + md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) #Check for orphan Index: /issm/trunk-jpl/src/m/mesh/meshconvert.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/meshconvert.py (revision 13984) @@ -27,5 +27,5 @@ #call Bamg - bamgmesh_out,bamggeom_out=BamgConvertMesh(index,x,y) + bamgmesh_out,bamggeom_out=BamgConvertMesh(index.astype(float),x,y) # plug results onto model @@ -35,8 +35,8 @@ md.mesh.x=bamgmesh_out['Vertices'][:,0].copy() md.mesh.y=bamgmesh_out['Vertices'][:,1].copy() - md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy() - md.mesh.edges=bamgmesh_out['IssmEdges'].copy() - md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy() - md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy() + md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int) + md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int) + md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int) + md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int) #Fill in rest of fields: @@ -46,11 +46,11 @@ md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0) md.mesh.z=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) - md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) - md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) - md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) - md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) + md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) + md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True return md Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 13984) @@ -16,5 +16,5 @@ #first, flag nodes that belong to the domain outline - flags=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0) + flags=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0) tips=rift.tips @@ -41,5 +41,5 @@ B=node_connected_to_tip - elements=numpy.empty(0) + elements=numpy.empty(0,int) while flags(B): #as long as B does not belong to the domain outline, keep looking. @@ -62,8 +62,8 @@ #replace tip in elements - newelements=md.mesh.elements[elements.astype(int)-1,:] + newelements=md.mesh.elements[elements-1,:] pos=numpy.nonzero(newelements==tip) newelements[pos]=num - md.mesh.elements[elements.astype(int)-1,:]=newelements + md.mesh.elements[elements-1,:]=newelements rift.tips=numpy.concatenate((rift.tips,num)) @@ -81,12 +81,12 @@ md.mesh.numberofvertices=numpy.size(md.mesh.x) md.mesh.z=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x)) - md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 + md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True md.rifts.numrifts=length(md.rifts.riftstruct) md.flowequation.element_equation=3.*numpy.ones(md.mesh.numberofelements) - md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) - md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) - md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) - md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) return md Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 13984) @@ -22,7 +22,10 @@ #Call MEX file - [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers) + [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements.astype(float),md.mesh.x,md.mesh.y,md.mesh.segments.astype(float),md.mesh.segmentmarkers.astype(float)) + md.mesh.elements=md.mesh.elements.astype(int) md.mesh.x=md.mesh.x.reshape(-1) md.mesh.y=md.mesh.y.reshape(-1) + md.mesh.segments=md.mesh.segments.astype(int) + md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int) if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct: raise RuntimeError("TriMeshProcessRifts did not find any rift") @@ -33,10 +36,10 @@ md.mesh.numberofvertices=numpy.size(md.mesh.x) md.mesh.z=numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x)) - md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 - md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) - md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) - md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) - md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) + md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool) #get coordinates of rift tips @@ -46,5 +49,5 @@ #In case we have rifts that open up the domain outline, we need to open them: - [flags,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0) + [flags,dum]=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0) found=0 for rift in md.rifts.riftstruct: Index: /issm/trunk-jpl/src/m/mesh/triangle.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/triangle.py (revision 13983) +++ /issm/trunk-jpl/src/m/mesh/triangle.py (revision 13984) @@ -44,4 +44,7 @@ #Mesh using TriMesh [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(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) #Fill in rest of fields: @@ -49,14 +52,16 @@ md.mesh.numberofvertices = numpy.size(md.mesh.x) md.mesh.z = numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices) - md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1. - md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices) - md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices) - md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements) - md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements) + md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True + md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices,bool) + md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements,bool) + md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements,bool) #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) + [md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements.astype(float), md.mesh.numberofvertices) + md.mesh.vertexconnectivity=md.mesh.vertexconnectivity.astype(int) + [md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements.astype(float), md.mesh.vertexconnectivity.astype(float)) + md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) #type of model Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/contourenvelope.py (revision 13983) +++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.py (revision 13984) @@ -40,7 +40,9 @@ #Computing connectivity if numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d: - [md.mesh.vertexconnectivity]=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices) + [md.mesh.vertexconnectivity]=NodeConnectivity(md.mesh.elements.astype(float),md.mesh.numberofvertices) + md.mesh.vertexconnectivity=md.mesh.vertexconnectivity.astype(int) if numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d: - [md.mesh.elementconnectivity]=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity) + [md.mesh.elementconnectivity]=ElementConnectivity(md.mesh.elements.astype(float),md.mesh.vertexconnectivity.astype(float)) + md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int) #get nodes inside profile @@ -63,5 +65,5 @@ if isfile: #get flag list of elements and nodes inside the contour - nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1) + nodein=ContourToMesh(mesh.elements.astype(float),mesh.x,mesh.y,file,'node',1) elemin=(numpy.sum(nodein(mesh.elements),axis=1)==numpy.size(mesh.elements,axis=1)) #modify element connectivity @@ -76,5 +78,5 @@ pos=numpy.nonzero(flags) elemin[pos]=1 - nodein[mesh.elements[pos,:].astype(int)-1]=1 + nodein[mesh.elements[pos,:]-1]=1 #modify element connectivity @@ -93,9 +95,9 @@ pos=numpy.nonzero(elementonboundary)[0] num_segments=numpy.size(pos) - segments=numpy.zeros((num_segments*3,3)) + segments=numpy.zeros((num_segments*3,3),int) count=0 for el1 in pos: - els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]].astype(int)-1 + els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]]-1 if numpy.size(els2)>1: flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:]) Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/setflowequation.py (revision 13983) +++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py (revision 13984) @@ -78,11 +78,11 @@ #Initialize node fields nodeonhutter=numpy.zeros(md.mesh.numberofvertices) - nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1 + nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=1 nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices) - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 nodeonl1l2=numpy.zeros(md.mesh.numberofvertices) - nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:].astype(int)-1]=1 + nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=1 nodeonpattyn=numpy.zeros(md.mesh.numberofvertices) - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 nodeonstokes=numpy.zeros(md.mesh.numberofvertices) noneflag=numpy.zeros(md.mesh.numberofelements) @@ -96,7 +96,7 @@ numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int) #find all the nodes on the boundary of the domain without icefront # fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront - fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements.astype(int)-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront + fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int) #find all the nodes on the boundary of the domain without icefront stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 #Then complete with NoneApproximation or the other model used if there is no stokes @@ -104,8 +104,8 @@ if any(pattynflag): #fill with pattyn pattynflag[numpy.logical_not(stokesflag)]=1 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 elif any(macayealflag): #fill with macayeal macayealflag[numpy.logical_not(stokesflag)]=1 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 else: #fill with none noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1 @@ -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.astype(int)-1,numpy.nonzero(nodeonmacayealpattyn)[0]) + matrixelements=ismember(md.mesh.elements-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 @@ -143,11 +143,11 @@ macayealpattynflag[numpy.nonzero(commonelements)]=1 nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 #rule out elements that don't touch the 2 boundaries 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) + elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) pos1=numpy.nonzero(elist==1)[0] macayealflag[pos[pos1]]=1 @@ -159,9 +159,9 @@ #Recompute nodes associated to these elements nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 nodeonpattyn[:]=0 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 nodeonmacayealpattyn[:]=0 - nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1 + nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=1 elif any(pattynflag) and any(stokesflag): #coupling pattyn stokes @@ -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.astype(int)-1,numpy.nonzero(nodeonpattynstokes)[0]) + matrixelements=ismember(md.mesh.elements-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 @@ -175,11 +175,11 @@ pattynstokesflag[numpy.nonzero(commonelements)]=1 nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 #rule out elements that don't touch the 2 boundaries 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) + elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool) pos1=numpy.nonzero(elist==1)[0] stokesflag[pos[pos1]]=1 @@ -191,9 +191,9 @@ #Recompute nodes associated to these elements nodeonstokes[:]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 nodeonpattyn[:]=0 - nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1 + nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1 nodeonpattynstokes[:]=0 - nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1 + nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=1 elif any(stokesflag) and any(macayealflag): @@ -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.astype(int)-1,numpy.nonzero(nodeonmacayealstokes)[0]) + matrixelements=ismember(md.mesh.elements-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 @@ -207,11 +207,11 @@ macayealstokesflag[numpy.nonzero(commonelements)]=1 nodeonstokes=numpy.zeros(md.mesh.numberofvertices) - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 #rule out elements that don't touch the 2 boundaries 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) + elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool) + elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1] ,axis=1).astype(bool) pos1=numpy.nonzero(elist==1)[0] macayealflag[pos[pos1]]=1 @@ -223,9 +223,9 @@ #Recompute nodes associated to these elements nodeonmacayeal[:]=0 - nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1 + nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1 nodeonstokes[:]=0 - nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1 + nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1 nodeonmacayealstokes[:]=0 - nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1 + nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=1 elif any(stokesflag) and any(hutterflag): Index: /issm/trunk-jpl/src/m/parameterization/setmask.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/setmask.py (revision 13983) +++ /issm/trunk-jpl/src/m/parameterization/setmask.py (revision 13984) @@ -37,15 +37,15 @@ vertexonfloatingice = numpy.zeros(md.mesh.numberofvertices,'bool') vertexongroundedice = numpy.zeros(md.mesh.numberofvertices,'bool') - vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:].astype(int)-1]=True + vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=True vertexonfloatingice[numpy.nonzero(numpy.logical_not(vertexongroundedice))]=True #}}} #Return: - md.mask.elementonfloatingice = elementonfloatingice.astype(float) - md.mask.vertexonfloatingice = vertexonfloatingice.astype(float) - md.mask.elementongroundedice = elementongroundedice.astype(float) - md.mask.vertexongroundedice = vertexongroundedice.astype(float) - md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices) - md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements) + md.mask.elementonfloatingice = elementonfloatingice + md.mask.vertexonfloatingice = vertexonfloatingice + md.mask.elementongroundedice = elementongroundedice + md.mask.vertexongroundedice = vertexongroundedice + md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices,bool) + md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements,bool) return md Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 13983) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 13984) @@ -105,5 +105,5 @@ elif isinstance(data,(list,tuple)): data=numpy.array(data).reshape(-1,1) - if len(data.shape) == 1: + if numpy.ndim(data) == 1: if numpy.size(data): data=data.reshape(numpy.size(data),1) @@ -138,5 +138,5 @@ elif isinstance(data,(list,tuple)): data=numpy.array(data).reshape(-1,1) - if len(data.shape) == 1: + if numpy.ndim(data) == 1: if numpy.size(data): data=data.reshape(numpy.size(data),1) @@ -171,5 +171,5 @@ elif isinstance(data,(list,tuple)): data=numpy.array(data).reshape(-1,1) - if len(data.shape) == 1: + if numpy.ndim(data) == 1: if numpy.size(data): data=data.reshape(numpy.size(data),1) @@ -207,5 +207,5 @@ elif isinstance(matrix,(list,tuple)): matrix=numpy.array(matrix).reshape(-1,1) - if len(matrix.shape) == 1: + if numpy.ndim(matrix) == 1: if numpy.size(matrix): matrix=matrix.reshape(numpy.size(matrix),1) @@ -231,5 +231,5 @@ elif isinstance(matrix,(list,tuple)): matrix=numpy.array(matrix).reshape(-1,1) - if len(matrix.shape) == 1: + if numpy.ndim(matrix) == 1: matrix=matrix.reshape(numpy.size(matrix),1)