Changeset 13984


Ignore:
Timestamp:
11/19/12 14:11:00 (12 years ago)
Author:
jschierm
Message:

CHG: Change properties in python md.mesh and md.mask to int or bool.

Location:
issm/trunk-jpl/src/m
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py

    r13970 r13984  
    2525                if not os.path.exists(icefrontfile):
    2626                        raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
    27                 [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
    28                 nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
     27                [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
     28                nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
    2929        else:
    30                 nodeonicefront=numpy.zeros((md.mesh.numberofvertices))
     30                nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool)
    3131
    3232#       pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
    … …  
    5656        #segment on Neumann (Ice Front)
    5757#       pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
    58         pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
     58        pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0]-1],nodeonicefront[md.mesh.segments[:,1]-1]))[0]
    5959        if   md.mesh.dimension==2:
    6060                pressureload=md.mesh.segments[pos,:]
    … …  
    6262#               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)];
    6363                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]))
    64                 pressureload=numpy.zeros((0,5))
     64                pressureload=numpy.zeros((0,5),int)
    6565                for i in xrange(1,md.mesh.numberoflayers):
    6666#                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
    … …  
    6969        #Add water or air enum depending on the element
    7070#       pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
    71         pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
     71        pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)))
    7272
    7373        #plug onto model

  • issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py

    r13857 r13984  
    2727                if not os.path.exists(icefrontfile):
    2828                        raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
    29                 [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
    30                 vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
     29                [nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
     30                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
    3131        else:
    3232                #Guess where the ice front is
    33                 vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices))
    34                 vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:].astype(int)-1]=1
    35                 vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice).astype(float)
     33                vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices),bool)
     34                vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:]-1]=True
     35                vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice)
    3636
    3737#       pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
    … …  
    6262        #segment on Neumann (Ice Front)
    6363#       pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2)));
    64         pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0].astype(int)-1],vertexonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
     64        pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0]-1],vertexonicefront[md.mesh.segments[:,1]-1]))[0]
    6565        if   md.mesh.dimension==2:
    6666                pressureload=md.mesh.segments[pos,:]
    … …  
    6868#               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)];
    6969                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]))
    70                 pressureload=numpy.zeros((0,5))
     70                pressureload=numpy.zeros((0,5),int)
    7171                for i in xrange(1,md.mesh.numberoflayers):
    7272#                       pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
    … …  
    7575        #Add water or air enum depending on the element
    7676#       pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))];
    77         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)))
     77        pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)+0*md.mask.elementongroundedice[pressureload[:,-1]-1].reshape(-1,1)))
    7878
    7979        #plug onto model

  • issm/trunk-jpl/src/m/classes/dependent.py

    r13865 r13984  
    4545                        #process the file and retrieve segments
    4646                        mesh=options.getfieldvalue('mesh')
    47                         self.segments=MeshProfileIntersection(mesh.elements,mesh.x,mesh.y,self.exp)
     47                        self.segments=MeshProfileIntersection(mesh.elements.astype(float),mesh.x,mesh.y,self.exp)
    4848        # }}}
    4949

  • issm/trunk-jpl/src/m/classes/initialization.py

    r13626 r13984  
    7070                        md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
    7171                        #Triangle with zero velocity
    72                         if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements.astype(int)-1]),axis=1)==0,\
    73                                                        numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements.astype(int)-1]),axis=1)==0)):
     72                        if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\
     73                                                       numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)):
    7474                                md.checkmessage("at least one triangle has all its vertices with a zero velocity")
    7575                if ThermalAnalysisEnum() in analyses:

  • issm/trunk-jpl/src/m/classes/mask.py

    r12958 r13984  
    3333                string="%s\n%s"%(string,fielddisplay(self,"elementonfloatingice","element on floating ice flags list"))
    3434                string="%s\n%s"%(string,fielddisplay(self,"vertexonfloatingice","vertex on floating ice flags list"))
    35                 string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice  list"))
     35                string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list"))
    3636                string="%s\n%s"%(string,fielddisplay(self,"vertexongroundedice","vertex on grounded ice flags list"))
    3737                string="%s\n%s"%(string,fielddisplay(self,"elementonwater","element on water flags list"))

  • issm/trunk-jpl/src/m/classes/mesh.py

    r13951 r13984  
    9595                string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (NaN for element on the upper layer)"))
    9696                string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (NaN for vertex on the lower surface)"))
    97                 string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (NaN for element on the lower layer"))
     97                string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (NaN for element on the lower layer)"))
    9898                string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list"))
    9999                string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)"))

  • issm/trunk-jpl/src/m/classes/model/model.py

    r13966 r13984  
    207207                #kick out all elements with 3 dirichlets
    208208                spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
    209                 spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1
     209                spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1
    210210                flag=numpy.ones(md1.mesh.numberofvertices)
    211211                flag[spc_node]=0
    212                 pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements.astype(int)-1],axis=1)))[0]
     212                pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0]
    213213                flag_elem[pos]=0
    214214
    215215                #extracted elements and nodes lists
    216216                pos_elem=numpy.nonzero(flag_elem)[0]
    217                 pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1
     217                pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1
    218218
    219219                #keep track of some fields
    … …  
    234234                elements_1=copy.deepcopy(md1.mesh.elements)
    235235                elements_2=elements_1[pos_elem,:]
    236                 elements_2[:,0]=Pnode[elements_2[:,0].astype(int)-1]
    237                 elements_2[:,1]=Pnode[elements_2[:,1].astype(int)-1]
    238                 elements_2[:,2]=Pnode[elements_2[:,2].astype(int)-1]
     236                elements_2[:,0]=Pnode[elements_2[:,0]-1]
     237                elements_2[:,1]=Pnode[elements_2[:,1]-1]
     238                elements_2[:,2]=Pnode[elements_2[:,2]-1]
    239239                if md1.mesh.dimension==3:
    240                         elements_2[:,3]=Pnode[elements_2[:,3].astype(int)-1]
    241                         elements_2[:,4]=Pnode[elements_2[:,4].astype(int)-1]
    242                         elements_2[:,5]=Pnode[elements_2[:,5].astype(int)-1]
     240                        elements_2[:,3]=Pnode[elements_2[:,3]-1]
     241                        elements_2[:,4]=Pnode[elements_2[:,4]-1]
     242                        elements_2[:,5]=Pnode[elements_2[:,5]-1]
    243243
    244244                #OK, now create the new model!
    … …  
    316316                        md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
    317317                        md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:]
    318                         md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0].astype(int)-1]
    319                         md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1].astype(int)-1]
    320                         md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2].astype(int)-1]
     318                        md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1]
     319                        md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1]-1]
     320                        md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2]-1]
    321321
    322322                        md2.mesh.x2d=md1.mesh.x[pos_node_2d]
    … …  
    324324
    325325                #Edges
    326                 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...
     326                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...
    327327                        #renumber first two columns
    328328                        pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
    329                         md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0].astype(int)-1]
    330                         md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1].astype(int)-1]
    331                         md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2].astype(int)-1]
    332                         md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3].astype(int)-1]
     329                        md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0]-1]
     330                        md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1]-1]
     331                        md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2]-1]
     332                        md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
    333333                        #remove edges when the 2 vertices are not in the domain.
    334334                        md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
    … …  
    361361                #recreate segments
    362362                if md1.mesh.dimension==2:
    363                         [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
    364                         [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
     363                        [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements.astype(float),md2.mesh.numberofvertices)
     364                        md2.mesh.vertexconnectivity=md2.mesh.vertexconnectivity.astype(int)
     365                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements.astype(float),md2.mesh.vertexconnectivity.astype(float))
     366                        md2.mesh.elementconnectivity=md2.mesh.elementconnectivity.astype(int)
    365367                        md2.mesh.segments=contourenvelope(md2)
    366                         md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2)
    367                         md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
     368                        md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool)
     369                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
    368370                else:
    369371                        #First do the connectivity for the contourenvelope in 2d
    370                         [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d)
    371                         [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)
     372                        [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d.astype(float),md2.mesh.numberofvertices2d)
     373                        md2.mesh.vertexconnectivity=md2.mesh.vertexconnectivity.astype(int)
     374                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d.astype(float),md2.mesh.vertexconnectivity.astype(float))
     375                        md2.mesh.elementconnectivity=md2.mesh.elementconnectivity.astype(int)
    372376                        md2.mesh.segments=contourenvelope(md2)
    373                         md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers)
    374                         md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
     377                        md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
     378                        md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
    375379                        md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
    376380                        #Then do it for 3d as usual
    377                         [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
    378                         [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
     381                        [md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements.astype(float),md2.mesh.numberofvertices)
     382                        md2.mesh.vertexconnectivity=md2.mesh.vertexconnectivity.astype(int)
     383                        [md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements.astype(float),md2.mesh.vertexconnectivity.astype(float))
     384                        md2.mesh.elementconnectivity=md2.mesh.elementconnectivity.astype(int)
    379385
    380386                #Boundary conditions: Dirichlets on new boundary
    381387                #Catch the elements that have not been extracted
    382388                orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
    383                 orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1
     389                orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1
    384390                #Figure out which node are on the boundary between md2 and md1
    385391                nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
    … …  
    443449
    444450                #Keep track of pos_node and pos_elem
    445                 md2.mesh.extractedvertices=pos_node.astype(float)+1
    446                 md2.mesh.extractedelements=pos_elem.astype(float)+1
     451                md2.mesh.extractedvertices=pos_node+1
     452                md2.mesh.extractedelements=pos_elem+1
    447453
    448454                return md2
    … …  
    526532
    527533                #Extrude elements
    528                 elements3d=numpy.empty((0,6))
     534                elements3d=numpy.empty((0,6),int)
    529535                for i in xrange(numlayers-1):
    530536                        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
    … …  
    603609
    604610                #bedinfo and surface info
    605                 md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',1)
    606                 md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',md.mesh.numberoflayers-1)
    607                 md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',1)
    608                 md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',md.mesh.numberoflayers)
     611                md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',1)
     612                md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',md.mesh.numberoflayers-1)
     613                md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
     614                md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
    609615
    610616                #elementstype
    … …  
    635641                #in 3d, pressureload: [node1 node2 node3 node4 element]
    636642                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
    637                 pressureload=numpy.empty((0,6))
     643                pressureload=numpy.empty((0,6),int)
    638644                for i in xrange(numlayers-1):
    639645                        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]))))
    … …  
    642648                #connectivity
    643649                md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
     650                md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(float)
    644651                md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=float('NaN')
    645652                for i in xrange(1,numlayers-1):
    … …  
    647654                                =md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
    648655                md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
     656                md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
    649657
    650658                #materials

  • issm/trunk-jpl/src/m/consistency/checkfield.py

    r13043 r13984  
    143143        if options.getfieldvalue('forcing',0):
    144144                if   numpy.size(field,0)==md.mesh.numberofvertices:
    145                         if len(numpy.shape(field))>1 and not numpy.size(field,1)==1:
     145                        if numpy.ndim(field)>1 and not numpy.size(field,1)==1:
    146146                                md = md.checkmessage(options.getfieldvalue('message',\
    147147                                        "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))

  • issm/trunk-jpl/src/m/extrusion/project3d.py

    r13151 r13984  
    3838
    3939        vector1d=False
    40         if isinstance(vector2d,numpy.ndarray) and len(vector2d.shape)==1:
     40        if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1:
    4141                vector1d=True
    4242                vector2d=vector2d.reshape(numpy.size(vector2d),1)
    … …  
    4949                #Initialize 3d vector
    5050                if   numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d:
    51                         projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices,  numpy.size(vector2d,axis=1)))
     51                        projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices,  numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
    5252                elif numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d+1:
    53                         projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))
     53                        projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
    5454                        projected_vector[-1,:]=vector2d[-1,:]
    5555                        vector2d=vector2d[:-1,:]
    … …  
    6868                #Initialize 3d vector
    6969                if   numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d:
    70                         projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements,  numpy.size(vector2d,axis=1)))
     70                        projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements,  numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
    7171                elif numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d+1:
    72                         projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))
     72                        projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
    7373                        projected_vector[-1,:]=vector2d[-1,:]
    7474                        vector2d=vector2d[:-1,:]

  • issm/trunk-jpl/src/m/geometry/FlagElements.py

    r13741 r13984  
    2424        if   isinstance(region,(str,unicode)):
    2525                if   not region:
    26                         flag=numpy.zeros(md.mesh.numberofelements,'bool')
     26                        flag=numpy.zeros(md.mesh.numberofelements,bool)
    2727                        invert=0
    2828                elif strcmpi(region,'all'):
    29                         flag=numpy.ones(md.mesh.numberofelements,'bool')
     29                        flag=numpy.ones(md.mesh.numberofelements,bool)
    3030                        invert=0
    3131                else:
    … …  
    4444                                xlim,ylim=basinzoom('basin',region)
    4545                                flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0])).astype(float)
    46                                 flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1)
     46                                flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool)
    4747                        else:
    4848                                #ok, flag elements
    49                                 [flag,dum]=ContourToMesh(md.mesh.elements[:,0:3].copy(),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1)
     49                                [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)
     50                                flag=flag.astype(bool)
    5051
    5152                if invert:

  • issm/trunk-jpl/src/m/geometry/GetAreas.py

    r13472 r13984  
    3333        #initialization
    3434        areas=numpy.zeros(nels)
    35         x1=x[index[:,0].astype(int)-1]
    36         x2=x[index[:,1].astype(int)-1]
    37         x3=x[index[:,2].astype(int)-1]
    38         y1=y[index[:,0].astype(int)-1]
    39         y2=y[index[:,1].astype(int)-1]
    40         y3=y[index[:,2].astype(int)-1]
     35        x1=x[index[:,0]-1]
     36        x2=x[index[:,1]-1]
     37        x3=x[index[:,2]-1]
     38        y1=y[index[:,0]-1]
     39        y2=y[index[:,1]-1]
     40        y3=y[index[:,2]-1]
    4141
    4242        #compute the volume of each element
    … …  
    4646        else:
    4747                #V=area(triangle)*1/3(z1+z2+z3)
    48                 thickness=numpy.mean(z[index[:,3:6].astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1])
     48                thickness=numpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1])
    4949                areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness
    5050

  • issm/trunk-jpl/src/m/mesh/ComputeHessian.py

    r13482 r13984  
    4545
    4646        #Compute gradient for each element
    47         grad_elx=numpy.sum(field[index.astype(int)-1,0]*alpha,axis=1)
    48         grad_ely=numpy.sum(field[index.astype(int)-1,0]*beta,axis=1)
     47        grad_elx=numpy.sum(field[index-1,0]*alpha,axis=1)
     48        grad_ely=numpy.sum(field[index-1,0]*beta,axis=1)
    4949
    5050        #Compute gradient for each node (average of the elements around)
    … …  
    5555
    5656        #Compute hessian for each element
    57         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)))
     57        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)))
    5858
    5959        if strcmpi(type,'node'):

  • issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py

    r13482 r13984  
    3939
    4040        #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
    41         x1=x[index[:,0].astype(int)-1]
    42         x2=x[index[:,1].astype(int)-1]
    43         x3=x[index[:,2].astype(int)-1]
    44         y1=y[index[:,0].astype(int)-1]
    45         y2=y[index[:,1].astype(int)-1]
    46         y3=y[index[:,2].astype(int)-1]
     41        x1=x[index[:,0]-1]
     42        x2=x[index[:,1]-1]
     43        x3=x[index[:,2]-1]
     44        y1=y[index[:,0]-1]
     45        y2=y[index[:,1]-1]
     46        y3=y[index[:,2]-1]
    4747        invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2))
    4848

  • issm/trunk-jpl/src/m/mesh/bamg.py

    r13917 r13984  
    276276                else:
    277277                        bamg_mesh.Vertices=numpy.hstack((md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),numpy.ones((md.mesh.numberofvertices,1))))
    278                         bamg_mesh.Triangles=numpy.hstack((md.mesh.elements,numpy.ones((md.mesh.numberofelements,1))))
     278                        bamg_mesh.Triangles=numpy.hstack((md.mesh.elements.astype(float),numpy.ones((md.mesh.numberofelements,1))))
    279279
    280280                if isinstance(md.rifts.riftstruct,dict):
    … …  
    320320        md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
    321321        md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
    322         md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy()
    323         md.mesh.edges=bamgmesh_out['IssmEdges'].copy()
    324         md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy()
    325         md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy()
     322        md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
     323        md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
     324        md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
     325        md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
    326326
    327327        #Fill in rest of fields:
    … …  
    331331        md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
    332332        md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
    333         md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
    334         md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices)
    335         md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
    336         md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
    337         md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
    338         md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
    339         md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
     333        md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
     334        md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
     335        md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
     336        md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
     337        md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
     338        md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
     339        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
    340340        md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
    341341        md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
     342        md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
    342343
    343344        #Check for orphan

  • issm/trunk-jpl/src/m/mesh/meshconvert.py

    r13436 r13984  
    2727
    2828        #call Bamg
    29         bamgmesh_out,bamggeom_out=BamgConvertMesh(index,x,y)
     29        bamgmesh_out,bamggeom_out=BamgConvertMesh(index.astype(float),x,y)
    3030
    3131        # plug results onto model
    … …  
    3535        md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
    3636        md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
    37         md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy()
    38         md.mesh.edges=bamgmesh_out['IssmEdges'].copy()
    39         md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy()
    40         md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy()
     37        md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
     38        md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
     39        md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
     40        md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
    4141
    4242        #Fill in rest of fields:
    … …  
    4646        md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
    4747        md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
    48         md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
    49         md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices)
    50         md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
    51         md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
    52         md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
    53         md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
    54         md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
     48        md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
     49        md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
     50        md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
     51        md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
     52        md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
     53        md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
     54        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
    5555
    5656        return md

  • issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py

    r13524 r13984  
    1616       
    1717                #first, flag nodes that belong to the domain outline
    18                 flags=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)
     18                flags=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)
    1919
    2020                tips=rift.tips
    … …  
    4141                        B=node_connected_to_tip
    4242
    43                         elements=numpy.empty(0)
     43                        elements=numpy.empty(0,int)
    4444
    4545                        while flags(B):    #as long as B does not belong to the domain outline, keep looking.
    … …  
    6262               
    6363                        #replace tip in elements
    64                         newelements=md.mesh.elements[elements.astype(int)-1,:]
     64                        newelements=md.mesh.elements[elements-1,:]
    6565                        pos=numpy.nonzero(newelements==tip)
    6666                        newelements[pos]=num
    67                         md.mesh.elements[elements.astype(int)-1,:]=newelements
     67                        md.mesh.elements[elements-1,:]=newelements
    6868                        rift.tips=numpy.concatenate((rift.tips,num))
    6969
    … …  
    8181        md.mesh.numberofvertices=numpy.size(md.mesh.x)
    8282        md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
    83         md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
    84         md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
     83        md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
     84        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
    8585        md.rifts.numrifts=length(md.rifts.riftstruct)
    8686        md.flowequation.element_equation=3.*numpy.ones(md.mesh.numberofelements)
    87         md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
    88         md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
    89         md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
    90         md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
     87        md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
     88        md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
     89        md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
     90        md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
    9191
    9292        return md

  • issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py

    r13710 r13984  
    2222
    2323        #Call MEX file
    24         [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)
     24        [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))
     25        md.mesh.elements=md.mesh.elements.astype(int)
    2526        md.mesh.x=md.mesh.x.reshape(-1)
    2627        md.mesh.y=md.mesh.y.reshape(-1)
     28        md.mesh.segments=md.mesh.segments.astype(int)
     29        md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
    2730        if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
    2831                raise RuntimeError("TriMeshProcessRifts did not find any rift")
    … …  
    3336        md.mesh.numberofvertices=numpy.size(md.mesh.x)
    3437        md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
    35         md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
    36         md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
    37         md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
    38         md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
    39         md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
    40         md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
     38        md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
     39        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
     40        md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
     41        md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
     42        md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
     43        md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
    4144
    4245        #get coordinates of rift tips
    … …  
    4649
    4750        #In case we have rifts that open up the domain outline, we need to open them:
    48         [flags,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)
     51        [flags,dum]=ContourToMesh(md.mesh.elements.astype(float),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0)
    4952        found=0
    5053        for rift in md.rifts.riftstruct:

  • issm/trunk-jpl/src/m/mesh/triangle.py

    r13520 r13984  
    4444        #Mesh using TriMesh
    4545        [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area)
     46        md.mesh.elements=md.mesh.elements.astype(int)
     47        md.mesh.segments=md.mesh.segments.astype(int)
     48        md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
    4649
    4750        #Fill in rest of fields:
    … …  
    4952        md.mesh.numberofvertices = numpy.size(md.mesh.x)
    5053        md.mesh.z = numpy.zeros(md.mesh.numberofvertices)
    51         md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices)
    52         md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1.
    53         md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices)
    54         md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices)
    55         md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements)
    56         md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements)
     54        md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool)
     55        md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True
     56        md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices,bool)
     57        md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices,bool)
     58        md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements,bool)
     59        md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements,bool)
    5760
    5861        #Now, build the connectivity tables for this mesh.
    59         [md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
    60         [md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
     62        [md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements.astype(float), md.mesh.numberofvertices)
     63        md.mesh.vertexconnectivity=md.mesh.vertexconnectivity.astype(int)
     64        [md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements.astype(float), md.mesh.vertexconnectivity.astype(float))
     65        md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
    6166
    6267        #type of model

  • issm/trunk-jpl/src/m/parameterization/contourenvelope.py

    r13857 r13984  
    4040        #Computing connectivity
    4141        if numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices and numpy.size(md.mesh.vertexconnectivity,axis=0)!=md.mesh.numberofvertices2d:
    42                 [md.mesh.vertexconnectivity]=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices)
     42                [md.mesh.vertexconnectivity]=NodeConnectivity(md.mesh.elements.astype(float),md.mesh.numberofvertices)
     43                md.mesh.vertexconnectivity=md.mesh.vertexconnectivity.astype(int)
    4344        if numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements and numpy.size(md.mesh.elementconnectivity,axis=0)!=md.mesh.numberofelements2d:
    44                 [md.mesh.elementconnectivity]=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity)
     45                [md.mesh.elementconnectivity]=ElementConnectivity(md.mesh.elements.astype(float),md.mesh.vertexconnectivity.astype(float))
     46                md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
    4547
    4648        #get nodes inside profile
    … …  
    6365                if isfile:
    6466                        #get flag list of elements and nodes inside the contour
    65                         nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1)
     67                        nodein=ContourToMesh(mesh.elements.astype(float),mesh.x,mesh.y,file,'node',1)
    6668                        elemin=(numpy.sum(nodein(mesh.elements),axis=1)==numpy.size(mesh.elements,axis=1))
    6769                        #modify element connectivity
    … …  
    7678                        pos=numpy.nonzero(flags)
    7779                        elemin[pos]=1
    78                         nodein[mesh.elements[pos,:].astype(int)-1]=1
     80                        nodein[mesh.elements[pos,:]-1]=1
    7981
    8082                        #modify element connectivity
    … …  
    9395        pos=numpy.nonzero(elementonboundary)[0]
    9496        num_segments=numpy.size(pos)
    95         segments=numpy.zeros((num_segments*3,3))
     97        segments=numpy.zeros((num_segments*3,3),int)
    9698        count=0
    9799
    98100        for el1 in pos:
    99                 els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]].astype(int)-1
     101                els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]]-1
    100102                if numpy.size(els2)>1:
    101103                        flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:])

  • issm/trunk-jpl/src/m/parameterization/setflowequation.py

    r13566 r13984  
    7878        #Initialize node fields
    7979        nodeonhutter=numpy.zeros(md.mesh.numberofvertices)
    80         nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1
     80        nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=1
    8181        nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices)
    82         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     82        nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1
    8383        nodeonl1l2=numpy.zeros(md.mesh.numberofvertices)
    84         nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:].astype(int)-1]=1
     84        nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=1
    8585        nodeonpattyn=numpy.zeros(md.mesh.numberofvertices)
    86         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
     86        nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1
    8787        nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
    8888        noneflag=numpy.zeros(md.mesh.numberofelements)
    … …  
    9696                                              numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
    9797#               fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
    98                 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
     98                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
    9999                stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0
    100                 nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     100                nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1
    101101
    102102        #Then complete with NoneApproximation or the other model used if there is no stokes
    … …  
    104104                if   any(pattynflag):    #fill with pattyn
    105105                        pattynflag[numpy.logical_not(stokesflag)]=1
    106                         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
     106                        nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1
    107107                elif any(macayealflag):    #fill with macayeal
    108108                        macayealflag[numpy.logical_not(stokesflag)]=1
    109                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     109                        nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1
    110110                else:    #fill with none
    111111                        noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1
    … …  
    137137                        nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1
    138138                        #Macayeal elements in contact with this layer become MacAyealPattyn elements
    139                         matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealpattyn)[0])
     139                        matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0])
    140140                        commonelements=numpy.sum(matrixelements,axis=1)!=0
    141141                        commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
    … …  
    143143                        macayealpattynflag[numpy.nonzero(commonelements)]=1
    144144                        nodeonmacayeal[:]=0
    145                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     145                        nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1
    146146
    147147                        #rule out elements that don't touch the 2 boundaries
    148148                        pos=numpy.nonzero(macayealpattynflag)[0]
    149149                        elist=numpy.zeros(numpy.size(pos),dtype=int)
    150                         elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
    151                         elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1).astype(bool)
     150                        elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
     151                        elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
    152152                        pos1=numpy.nonzero(elist==1)[0]
    153153                        macayealflag[pos[pos1]]=1
    … …  
    159159                        #Recompute nodes associated to these elements
    160160                        nodeonmacayeal[:]=0
    161                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     161                        nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1
    162162                        nodeonpattyn[:]=0
    163                         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
     163                        nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1
    164164                        nodeonmacayealpattyn[:]=0
    165                         nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1
     165                        nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=1
    166166
    167167                elif any(pattynflag) and any(stokesflag):    #coupling pattyn stokes
    … …  
    169169                        nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1
    170170                        #Stokes elements in contact with this layer become PattynStokes elements
    171                         matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonpattynstokes)[0])
     171                        matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0])
    172172                        commonelements=numpy.sum(matrixelements,axis=1)!=0
    173173                        commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
    … …  
    175175                        pattynstokesflag[numpy.nonzero(commonelements)]=1
    176176                        nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
    177                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     177                        nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1
    178178
    179179                        #rule out elements that don't touch the 2 boundaries
    180180                        pos=numpy.nonzero(pattynstokesflag)[0]
    181181                        elist=numpy.zeros(numpy.size(pos),dtype=int)
    182                         elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
    183                         elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
     182                        elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
     183                        elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
    184184                        pos1=numpy.nonzero(elist==1)[0]
    185185                        stokesflag[pos[pos1]]=1
    … …  
    191191                        #Recompute nodes associated to these elements
    192192                        nodeonstokes[:]=0
    193                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     193                        nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1
    194194                        nodeonpattyn[:]=0
    195                         nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
     195                        nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=1
    196196                        nodeonpattynstokes[:]=0
    197                         nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1
     197                        nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=1
    198198
    199199                elif any(stokesflag) and any(macayealflag):
    … …  
    201201                        nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1
    202202                        #Stokes elements in contact with this layer become MacAyealStokes elements
    203                         matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealstokes)[0])
     203                        matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0])
    204204                        commonelements=numpy.sum(matrixelements,axis=1)!=0
    205205                        commonelements[numpy.nonzero(macayealflag)]=0    #only one layer: the elements previously in macayeal
    … …  
    207207                        macayealstokesflag[numpy.nonzero(commonelements)]=1
    208208                        nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
    209                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     209                        nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1
    210210
    211211                        #rule out elements that don't touch the 2 boundaries
    212212                        pos=numpy.nonzero(macayealstokesflag)[0]
    213213                        elist=numpy.zeros(numpy.size(pos),dtype=int)
    214                         elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
    215                         elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1).astype(bool)
     214                        elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
     215                        elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
    216216                        pos1=numpy.nonzero(elist==1)[0]
    217217                        macayealflag[pos[pos1]]=1
    … …  
    223223                        #Recompute nodes associated to these elements
    224224                        nodeonmacayeal[:]=0
    225                         nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
     225                        nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=1
    226226                        nodeonstokes[:]=0
    227                         nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
     227                        nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=1
    228228                        nodeonmacayealstokes[:]=0
    229                         nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1
     229                        nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=1
    230230
    231231                elif any(stokesflag) and any(hutterflag):

  • issm/trunk-jpl/src/m/parameterization/setmask.py

    r13449 r13984  
    3737        vertexonfloatingice = numpy.zeros(md.mesh.numberofvertices,'bool')
    3838        vertexongroundedice = numpy.zeros(md.mesh.numberofvertices,'bool')
    39         vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:].astype(int)-1]=True
     39        vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=True
    4040        vertexonfloatingice[numpy.nonzero(numpy.logical_not(vertexongroundedice))]=True
    4141        #}}}
    4242
    4343        #Return:
    44         md.mask.elementonfloatingice = elementonfloatingice.astype(float)
    45         md.mask.vertexonfloatingice = vertexonfloatingice.astype(float)
    46         md.mask.elementongroundedice = elementongroundedice.astype(float)
    47         md.mask.vertexongroundedice = vertexongroundedice.astype(float)
    48         md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices)
    49         md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements)
     44        md.mask.elementonfloatingice = elementonfloatingice
     45        md.mask.vertexonfloatingice = vertexonfloatingice
     46        md.mask.elementongroundedice = elementongroundedice
     47        md.mask.vertexongroundedice = vertexongroundedice
     48        md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices,bool)
     49        md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements,bool)
    5050
    5151        return md

  • issm/trunk-jpl/src/m/solve/WriteData.py

    r13899 r13984  
    105105                elif isinstance(data,(list,tuple)):
    106106                        data=numpy.array(data).reshape(-1,1)
    107                 if len(data.shape) == 1:
     107                if numpy.ndim(data) == 1:
    108108                        if numpy.size(data):
    109109                                data=data.reshape(numpy.size(data),1)
    … …  
    138138                elif isinstance(data,(list,tuple)):
    139139                        data=numpy.array(data).reshape(-1,1)
    140                 if len(data.shape) == 1:
     140                if numpy.ndim(data) == 1:
    141141                        if numpy.size(data):
    142142                                data=data.reshape(numpy.size(data),1)
    … …  
    171171                elif isinstance(data,(list,tuple)):
    172172                        data=numpy.array(data).reshape(-1,1)
    173                 if len(data.shape) == 1:
     173                if numpy.ndim(data) == 1:
    174174                        if numpy.size(data):
    175175                                data=data.reshape(numpy.size(data),1)
    … …  
    207207                        elif isinstance(matrix,(list,tuple)):
    208208                                matrix=numpy.array(matrix).reshape(-1,1)
    209                         if len(matrix.shape) == 1:
     209                        if numpy.ndim(matrix) == 1:
    210210                                if numpy.size(matrix):
    211211                                        matrix=matrix.reshape(numpy.size(matrix),1)
    … …  
    231231                        elif isinstance(matrix,(list,tuple)):
    232232                                matrix=numpy.array(matrix).reshape(-1,1)
    233                         if len(matrix.shape) == 1:
     233                        if numpy.ndim(matrix) == 1:
    234234                                matrix=matrix.reshape(numpy.size(matrix),1)
    235235
Note: See TracChangeset for help on using the changeset viewer.