Changeset 21254


Ignore:
Timestamp:
10/11/16 01:24:07 (8 years ago)
Author:
bdef
Message:

CHG: replacing npy by np for numpy abreviation

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

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py

    r17567 r21254  
    11import os
    2 import numpy as npy
     2import numpy as np
    33def PattynSMB(md,Tf):
    44        """
     
    2525
    2626        #Double check lat and long exist:
    27         if npy.any(npy.isnan(md.mesh.lat)):
     27        if np.any(np.isnan(md.mesh.lat)):
    2828                raise IOError('PattynSMB error message: md.mesh.lat field required')
    2929
     
    3636        # Calculate summer temperature, Eqn (12)
    3737        # No melting at PIG in mean conditions - need about 6 degress Tf to start having a negative yearly SMB
    38         Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*npy.abs(md.mesh.lat) + Tf
     38        Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*np.abs(md.mesh.lat) + Tf
    3939        Tms= Tms[0]
    4040
     
    4343
    4444        # Calculate Ablation, Eqn (10) (use for both Antarctica & Greenland), max melt is 10m/a
    45         ABLdot=0.*npy.ones(md.mesh.numberofvertices)
    46         pos=npy.nonzero(Tms>=0)
    47         ABLdot[pos]=npy.minimum(1.4*Tms[pos],10)
     45        ABLdot=0.*np.ones(md.mesh.numberofvertices)
     46        pos=np.nonzero(Tms>=0)
     47        ABLdot[pos]=np.minimum(1.4*Tms[pos],10)
    4848
    4949        smb=ACCdot-ABLdot
  • issm/trunk-jpl/src/m/classes/outputdefinition.py

    r21049 r21254  
    22from checkfield import checkfield
    33from WriteData import WriteData
    4 import numpy as npy
     4import numpy as np
    55
    66class outputdefinition(object):
     
    4040                        data.append(classdefinition)
    4141
    42                 data=npy.unique(data);
     42                data=np.unique(data);
    4343                WriteData(fid,prefix,'data',data,'name','md.outputdefinition.list','format','StringArray');
    4444        # }}}
  • issm/trunk-jpl/src/m/classes/qmu.py

    r21049 r21254  
    55from checkfield import checkfield
    66from WriteData import WriteData
    7 import MatlabFuncs as m
    87
    98class qmu(object):
     
    4140
    4241                s+="%s\n" % fielddisplay(self,'isdakota','is qmu analysis activated?')
    43                 for i,variable in enumerate(self.variables.iteritems()):
    44                         s+="         variables%s:  (arrays of each variable class)\n" % \
    45                                         string_dim(self.variables,i)
     42                for variable in self.variables:
     43                        s+="         variables%s:  (arrays of each variable class)\n" % len(variable)
    4644                        fnames=vars(variable)
    4745                        maxlen=0
     
    5351                                                (maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname)))
    5452
    55                 for i,response in enumerate(self.responses.iteritems()):
    56                         s+="         responses%s:  (arrays of each response class)\n" % \
    57                                         string_dim(self.responses,i)
     53                for response in self.responses:
     54                        s+="         responses%s:  (arrays of each response class)\n" % len(responses)
    5855                        fnames=vars(response)
    5956                        maxlen=0
     
    6764                s+="%s\n" % fielddisplay(self,'numberofresponses','number of responses')
    6865
    69                 for i,method in enumerate(self.method.iteritems()):
     66                for method in self.method:
    7067                        if isinstance(method,'dakota_method'):
    7168                                s+="            method%s :    '%s'\n" % \
    72                                                 (string_dim(method,i),method.method)
     69                                                (len(method),method.method)
    7370
    74                 for i,param in enumerate(self.params.iteritems()):
    75                         s+="         params%s:  (array of method-independent parameters)\n" % \
    76                                         string_dim(self.params,i)
     71                for param in self.params:
     72                        s+="         params%s:  (array of method-independent parameters)\n" % len(param)
    7773                        fnames=vars(param)
    7874                        maxlen=0
     
    8480                                                (maxlen+1,fname,any2str(getattr(param,fname)))
    8581
    86                 for i,result in enumerate(self.results.iteritems()):
    87                         s+="         results%s:  (information from dakota files)\n" % \
    88                                         string_dim(self.results,i)
     82                for result in self.results:
     83                        s+="         results%s:  (information from dakota files)\n" % len(self.result)
    8984                        fnames=vars(result)
    9085                        maxlen=0
     
    132127                                md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas")
    133128
    134                 if not m.strcmpi(md.cluster.name,'none'):
     129                if md.cluster.name!='none':
    135130                        if not md.settings.waitonlock:
    136131                                md.checkmessage("waitonlock should be activated when running qmu in parallel mode!")
  • issm/trunk-jpl/src/m/classes/settings.py

    r21049 r21254  
    1616                self.lowmem              = 0
    1717                self.output_frequency    = 0
    18                 self.recording_frequency    = 0
     18                self.recording_frequency = 0
    1919                self.waitonlock          = 0
    2020
  • issm/trunk-jpl/src/m/coordsystems/ll2xy.py

    r17712 r21254  
    1 import numpy as npy
     1import numpy as np
    22
    33def ll2xy(lat,lon,sgn=-1,central_meridian=0,standard_parallel=71):
     
    3636        ex2 = .006693883
    3737        # Eccentricity of the Hughes ellipsoid
    38         ex = npy.sqrt(ex2)
     38        ex = np.sqrt(ex2)
    3939       
    40         latitude = npy.abs(lat) * npy.pi/180.
    41         longitude = (lon + delta) * npy.pi/180.
     40        latitude = np.abs(lat) * np.pi/180.
     41        longitude = (lon + delta) * np.pi/180.
    4242       
    4343        # compute X and Y in grid coordinates.
    44         T = npy.tan(npy.pi/4-latitude/2) / ((1-ex*npy.sin(latitude))/(1+ex*npy.sin(latitude)))**(ex/2)
     44        T = np.tan(np.pi/4-latitude/2) / ((1-ex*np.sin(latitude))/(1+ex*np.sin(latitude)))**(ex/2)
    4545       
    4646        if (90 - slat) <  1.e-5:
    47                 rho = 2.*re*T/npy.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
     47                rho = 2.*re*T/np.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
    4848        else:
    49                 sl  = slat*npy.pi/180.
    50                 tc  = npy.tan(npy.pi/4.-sl/2.)/((1.-ex*npy.sin(sl))/(1.+ex*npy.sin(sl)))**(ex/2.)
    51                 mc  = npy.cos(sl)/npy.sqrt(1.0-ex2*(npy.sin(sl)**2))
     49                sl  = slat*np.pi/180.
     50                tc  = np.tan(np.pi/4.-sl/2.)/((1.-ex*np.sin(sl))/(1.+ex*np.sin(sl)))**(ex/2.)
     51                mc  = np.cos(sl)/np.sqrt(1.0-ex2*(np.sin(sl)**2))
    5252                rho = re*mc*T/tc
    5353       
    54         y = -rho * sgn * npy.cos(sgn*longitude)
    55         x =  rho * sgn * npy.sin(sgn*longitude)
     54        y = -rho * sgn * np.cos(sgn*longitude)
     55        x =  rho * sgn * np.sin(sgn*longitude)
    5656
    57         cnt1=npy.nonzero(latitude>= npy.pi/2.)[0]
     57        cnt1=np.nonzero(latitude>= np.pi/2.)[0]
    5858       
    5959        if cnt1:
  • issm/trunk-jpl/src/m/coordsystems/xy2ll.py

    r17712 r21254  
    1 import numpy as npy
     1import numpy as np
    22from math import pi
    33
     
    3939        # if x,y passed as lists, convert to numpy arrays
    4040        if type(x) != "numpy.ndarray":
    41                 x=npy.array(x)
     41                x=np.array(x)
    4242        if type(y) != "numpy.ndarray":
    43                 y=npy.array(y)
     43                y=np.array(y)
    4444
    4545        ## Conversion constant from degrees to radians
     
    5050        ex2 = .006693883
    5151        ## Eccentricity of the Hughes ellipsoid
    52         ex = npy.sqrt(ex2)
     52        ex = np.sqrt(ex2)
    5353       
    5454        sl = slat*pi/180.
    55         rho = npy.sqrt(x**2 + y**2)
    56         cm = npy.cos(sl) / npy.sqrt(1.0 - ex2 * (npy.sin(sl)**2))
    57         T = npy.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*npy.sin(sl)) / (1.0 + ex*npy.sin(sl)))**(ex / 2.0)
     55        rho = np.sqrt(x**2 + y**2)
     56        cm = np.cos(sl) / np.sqrt(1.0 - ex2 * (np.sin(sl)**2))
     57        T = np.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*np.sin(sl)) / (1.0 + ex*np.sin(sl)))**(ex / 2.0)
    5858       
    5959        if abs(slat-90.) < 1.e-5:
    60                 T = rho*npy.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
     60                T = rho*np.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
    6161        else:
    6262                T = rho * T / (re * cm)
    6363       
    64         chi = (pi / 2.0) - 2.0 * npy.arctan(T)
     64        chi = (pi / 2.0) - 2.0 * np.arctan(T)
    6565        lat = chi + ((ex2 / 2.0) + (5.0 * ex2**2.0 / 24.0) + (ex2**3.0 / 12.0)) * \
    66                 npy.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
    67                 npy.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * npy.sin(6.0 * chi)
     66                np.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
     67                np.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * np.sin(6.0 * chi)
    6868       
    6969        lat = sgn * lat
    70         lon = npy.arctan2(sgn * x,-sgn * y)
     70        lon = np.arctan2(sgn * x,-sgn * y)
    7171        lon = sgn * lon
    7272       
    73         res1 = npy.nonzero(rho <= 0.1)[0]
     73        res1 = np.nonzero(rho <= 0.1)[0]
    7474        if len(res1) > 0:
    7575                lat[res1] = 90. * sgn
  • issm/trunk-jpl/src/m/exp/expcoarsen.py

    r17480 r21254  
    11import os.path
    2 import numpy as npy
     2import numpy as np
    33from collections import OrderedDict
    44from expread import expread
     
    3434        for contour in  contours:
    3535               
    36                 numpoints=npy.size(contour['x'])
     36                numpoints=np.size(contour['x'])
    3737
    3838                j=0
     
    4444
    4545                        #see whether we keep this point or not
    46                         distance=npy.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
     46                        distance=np.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
    4747                        if distance<resolution and j<numpoints-2:   #do not remove last point
    48                                 x=npy.delete(x,j+1,0)
    49                                 y=npy.delete(y,j+1,0)
     48                                x=np.delete(x,j+1,0)
     49                                y=np.delete(y,j+1,0)
    5050                                numpoints=numpoints-1
    5151                        else:
    52                                 division=int(npy.floor(distance/resolution)+1)
     52                                division=int(np.floor(distance/resolution)+1)
    5353                                if division>=2:
    54                                         xi=npy.linspace(x[j],x[j+1],division)
    55                                         yi=npy.linspace(y[j],y[j+1],division)
     54                                        xi=np.linspace(x[j],x[j+1],division)
     55                                        yi=np.linspace(y[j],y[j+1],division)
    5656                                       
    57                                         x=npy.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
    58                                         y=npy.hstack((y[0:j+1],yi[1:-1],y[j+1:]))
     57                                        x=np.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
     58                                        y=np.hstack((y[0:j+1],yi[1:-1],y[j+1:]))
    5959
    6060                                        #update current point
     
    6565                                        j=j+1
    6666               
    67                 if npy.size(x)>1:
     67                if np.size(x)>1:
    6868                        #keep the (x,y) contour arond
    6969                        newcontour=OrderedDict()
    70                         newcontour['nods']=npy.size(x)
     70                        newcontour['nods']=np.size(x)
    7171                        newcontour['density']=contour['density']
    7272                        newcontour['x']=x
  • issm/trunk-jpl/src/m/exp/expdisp.py

    r19421 r21254  
    11from expread import expread
    2 import numpy as npy
     2import numpy as np
    33
    44def expdisp(domainoutline,ax,linestyle='--k',linewidth=1,unitmultiplier=1.):
  • issm/trunk-jpl/src/m/extrusion/DepthAverage.py

    r17810 r21254  
    1 import numpy as npy
     1import numpy as np
    22from project2d import project2d
    33
     
    1919
    2020        # coerce to array in case float is passed
    21         if type(vector)!=npy.ndarray:
     21        if type(vector)!=np.ndarray:
    2222                print 'coercing array'
    23                 vector=npy.array(value)
     23                vector=np.array(value)
    2424
    2525        vec2d=False
     
    3030        #nods data
    3131        if vector.shape[0]==md.mesh.numberofvertices:
    32                 vector_average=npy.zeros(md.mesh.numberofvertices2d)
     32                vector_average=np.zeros(md.mesh.numberofvertices2d)
    3333                for i in xrange(1,md.mesh.numberoflayers):
    3434                        vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i))
     
    3737        #element data
    3838        elif vector.shape[0]==md.mesh.numberofelements:
    39                 vector_average=npy.zeros(md.mesh.numberofelements2d)
     39                vector_average=np.zeros(md.mesh.numberofelements2d)
    4040                for i in xrange(1,md.mesh.numberoflayers):
    4141                        vector_average=vector_average+project2d(md,vector,i)*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i))
  • issm/trunk-jpl/src/m/extrusion/project2d.py

    r18216 r21254  
    1 import numpy as npy
     1import numpy as np
    22
    33def project2d(md3d,value,layer):
     
    2525       
    2626        # coerce to array in case float is passed
    27         if type(value)!=npy.ndarray:
     27        if type(value)!=np.ndarray:
    2828                print 'coercing array'
    29                 value=npy.array(value)
     29                value=np.array(value)
    3030
    3131        vec2d=False
  • issm/trunk-jpl/src/m/geometry/slope.py

    r17686 r21254  
    1 import numpy as npy
     1import numpy as np
    22from GetNodalFunctionsCoeff import  GetNodalFunctionsCoeff
    33
     
    3333        alpha,beta=GetNodalFunctionsCoeff(index,x,y)[0:2]
    3434
    35         summation=npy.array([[1],[1],[1]])
    36         sx=npy.dot(surf[index-1]*alpha,summation).reshape(-1,)
    37         sy=npy.dot(surf[index-1]*beta,summation).reshape(-1,)
     35        summation=np.array([[1],[1],[1]])
     36        sx=np.dot(surf[index-1]*alpha,summation).reshape(-1,)
     37        sy=np.dot(surf[index-1]*beta,summation).reshape(-1,)
    3838
    39         s=npy.sqrt(sx**2+sy**2)
     39        s=np.sqrt(sx**2+sy**2)
    4040
    4141        if md.mesh.dimension()==3:
    4242                sx=project3d(md,'vector',sx,'type','element')
    4343                sy=project3d(md,'vector',sy,'type','element')
    44                 s=npy.sqrt(sx**2+sy**2)
     44                s=np.sqrt(sx**2+sy**2)
    4545
    4646        return (sx,sy,s)
  • issm/trunk-jpl/src/m/interp/SectionValues.py

    r20910 r21254  
    11import os
    22from expread import expread
    3 import numpy as npy
     3import numpy as np
    44from project2d import project2d
    55#from InterpFromMesh2d import InterpFromMesh2d
     
    4141
    4242        #initialization
    43         X=npy.array([]) #X-coordinate
    44         Y=npy.array([]) #Y-coordinate
    45         S=npy.array([0.])  #curvilinear coordinate
     43        X=np.array([]) #X-coordinate
     44        Y=np.array([]) #Y-coordinate
     45        S=np.array([0.])  #curvilinear coordinate
    4646       
    4747        for i in xrange(nods-1):
     
    5353                s_start=S[-1]
    5454       
    55                 length_segment=npy.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
    56                 portion=npy.ceil(length_segment/res_h)
     55                length_segment=np.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
     56                portion=np.ceil(length_segment/res_h)
    5757       
    58                 x_segment=npy.zeros(portion)
    59                 y_segment=npy.zeros(portion)
    60                 s_segment=npy.zeros(portion)
     58                x_segment=np.zeros(portion)
     59                y_segment=np.zeros(portion)
     60                s_segment=np.zeros(portion)
    6161
    6262                for j in xrange(int(portion)):
     
    6666       
    6767                #plug into X and Y
    68                 X=npy.append(X,x_segment)
    69                 Y=npy.append(Y,y_segment)
    70                 S=npy.append(S,s_segment)
     68                X=np.append(X,x_segment)
     69                Y=np.append(Y,y_segment)
     70                S=np.append(S,s_segment)
    7171
    72         X=npy.append(X,x[nods-1])
    73         Y=npy.append(Y,y[nods-1])
     72        X=np.append(X,x[nods-1])
     73        Y=np.append(Y,y[nods-1])
    7474       
    7575        #Number of nodes:
     
    7777       
    7878        #Compute Z
    79         Z=npy.zeros(numberofnodes)
     79        Z=np.zeros(numberofnodes)
    8080       
    8181        #New mesh and Data interpolation
     
    8888       
    8989                #Compute index
    90                 index=npy.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T
     90                index=np.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T
    9191       
    9292        else:
     
    102102       
    103103                #Some useful parameters
    104                 layers=int(npy.ceil(npy.mean(md.geometry.thickness)/res_v))
     104                layers=int(np.ceil(np.mean(md.geometry.thickness)/res_v))
    105105                nodesperlayer=int(numberofnodes)
    106106                nodestot=int(nodesperlayer*layers)
     
    109109       
    110110                #initialization
    111                 X3=npy.zeros(nodesperlayer*layers)
    112                 Y3=npy.zeros(nodesperlayer*layers)
    113                 Z3=npy.zeros(nodesperlayer*layers)
    114                 S3=npy.zeros(nodesperlayer*layers)
    115                 index3=npy.zeros((elementstot,4))
     111                X3=np.zeros(nodesperlayer*layers)
     112                Y3=np.zeros(nodesperlayer*layers)
     113                Z3=np.zeros(nodesperlayer*layers)
     114                S3=np.zeros(nodesperlayer*layers)
     115                index3=np.zeros((elementstot,4))
    116116       
    117117                #Get new coordinates in 3d
     
    123123       
    124124                        if i<layers-1:  #Build index3 with quads
    125                                 ids=npy.vstack((npy.arange(i,nodestot-layers,layers),npy.arange(i+1,nodestot-layers,layers),npy.arange(i+layers+1,nodestot,layers),npy.arange(i+layers,nodestot,layers))).T
     125                                ids=np.vstack((np.arange(i,nodestot-layers,layers),np.arange(i+1,nodestot-layers,layers),np.arange(i+layers+1,nodestot,layers),np.arange(i+layers,nodestot,layers))).T
    126126                                index3[(i-1)*elementsperlayer:i*elementsperlayer,:]=ids
    127127
    128128                #Interpolation of data on specified points
    129                 data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,npy.nan)[0]
     129                data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,np.nan)[0]
    130130       
    131131                #build outputs
  • issm/trunk-jpl/src/m/interp/averaging.py

    r20144 r21254  
    1 import numpy as npy
     1import numpy as np
    22from GetAreas import GetAreas
    33import MatlabFuncs as m
     
    3939        #initialization
    4040        if layer==0:
    41                 weights=npy.zeros(md.mesh.numberofvertices,)
     41                weights=np.zeros(md.mesh.numberofvertices,)
    4242                data=data.flatten(1)
    4343        else:
    44                 weights=npy.zeros(md.mesh.numberofvertices2d,)
     44                weights=np.zeros(md.mesh.numberofvertices2d,)
    4545                data=data[(layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:]
    4646       
     
    6969        index=index-1 # since python indexes starting from zero
    7070        line=index.flatten(1)
    71         areas=npy.vstack(areas).reshape(-1,)
    72         summation=1./rep*npy.ones(rep,)
     71        areas=np.vstack(areas).reshape(-1,)
     72        summation=1./rep*np.ones(rep,)
    7373        linesize=rep*numberofelements
    7474       
    7575        #update weights that holds the volume of all the element holding the node i
    76         weights=csc_matrix( (npy.tile(areas,(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
     76        weights=csc_matrix( (np.tile(areas,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
    7777       
    7878        #initialization
    7979        if len(data)==numberofelements:
    80                 average_node=csc_matrix( (npy.tile(areas*data,(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
     80                average_node=csc_matrix( (np.tile(areas*data,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
    8181                average_node=average_node/weights
    8282                average_node = csc_matrix(average_node)
     
    8585
    8686        #loop over iteration
    87         for i in npy.arange(1,iterations+1):
    88                 average_el=npy.asarray(npy.dot(average_node.todense()[index].reshape(numberofelements,rep),npy.vstack(summation))).reshape(-1,)
    89                 average_node=csc_matrix( (npy.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
     87        for i in np.arange(1,iterations+1):
     88                average_el=np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements,rep),np.vstack(summation))).reshape(-1,)
     89                average_node=csc_matrix( (np.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
    9090                average_node=average_node/weights
    9191                average_node=csc_matrix(average_node)
    9292       
    9393        #return output as a full matrix (C code do not like sparse matrices)
    94         average=npy.asarray(average_node.todense()).reshape(-1,)
     94        average=np.asarray(average_node.todense()).reshape(-1,)
    9595
    9696        return average
  • issm/trunk-jpl/src/m/interp/holefiller.py

    r17867 r21254  
    1 import numpy as npy
     1import numpy as np
    22from scipy.spatial import cKDTree
    33
     
    3030        filled=data
    3131       
    32         XYGood=npy.dstack([x[goodids],y[goodids]])[0]
    33         XYBad=npy.dstack([x[badids],y[badids]])[0]
     32        XYGood=np.dstack([x[goodids],y[goodids]])[0]
     33        XYBad=np.dstack([x[badids],y[badids]])[0]
    3434        tree=cKDTree(XYGood)
    3535        nearest=tree.query(XYBad,k=knn)[1]
     
    4242                        for j in range(knn):
    4343                                neardat.append(filled[goodids][nearest[i][j]])
    44                                 filled[badids[i]]=npy.mean(neardat)
     44                                filled[badids[i]]=np.mean(neardat)
    4545                               
    4646        return filled
  • issm/trunk-jpl/src/m/interp/interp.py

    r18109 r21254  
    11# module for inperpolating/smoothing data
    2 import numpy as npy
     2import numpy as np
    33from scipy.interpolate import CloughTocher2DInterpolator, Rbf
    44from scipy.spatial import cKDTree
     
    5151        xlim=[min(xi)-dx,max(xi)+dx]
    5252        ylim=[min(yi)-dy,max(yi)+dy]
    53         xflag=npy.logical_and(x>xlim[0],x<xlim[1])
    54         yflag=npy.logical_and(y>ylim[0],y<ylim[1])
    55         bothind=npy.nonzero(npy.logical_and(xflag,yflag))
     53        xflag=np.logical_and(x>xlim[0],x<xlim[1])
     54        yflag=np.logical_and(y>ylim[0],y<ylim[1])
     55        bothind=np.nonzero(np.logical_and(xflag,yflag))
    5656        subdata=data[bothind]
    5757        subx=x[bothind]
    5858        suby=y[bothind]
    59         points=npy.array([subx,suby]).T
     59        points=np.array([subx,suby]).T
    6060
    6161        # mask out any nan's in the data and corresponding coordinate points
    62         mask=npy.isnan(subdata)
    63         ind=npy.nonzero(mask)[0]
     62        mask=np.isnan(subdata)
     63        ind=np.nonzero(mask)[0]
    6464        if len(ind) and fill_nans:
    6565                print "         WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully."
    66         subdata=npy.delete(subdata,ind)
    67         points=npy.delete(points,ind,axis=0)
     66        subdata=np.delete(subdata,ind)
     67        points=np.delete(points,ind,axis=0)
    6868
    6969        if maxiter:
     
    7676        if not fill_nans:
    7777                # identify nan's in xi,yi using nearest neighbors
    78                 xyinterp=npy.dstack([xi,yi])[0]
    79                 xg,yg=npy.meshgrid(subx,suby)
    80                 xydata=npy.dstack([subx,suby])[0]
     78                xyinterp=np.dstack([xi,yi])[0]
     79                xg,yg=np.meshgrid(subx,suby)
     80                xydata=np.dstack([subx,suby])[0]
    8181                tree=cKDTree(xydata)
    8282                nearest=tree.query(xyinterp)[1]
    83                 pos=npy.nonzero(npy.isnan(subdata[nearest]))
     83                pos=np.nonzero(np.isnan(subdata[nearest]))
    8484                interpdata[pos]=subdata[nearest][pos]
    8585
    8686        return interpdata
    8787#}}}
    88 def GridSplineToMesh2d(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False):#{{{
     88def GridSplineToMesh2d(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False):#{{{
    8989        '''
    9090        python analog to InterpFromGridToMesh.  This routine uses
     
    108108
    109109        Usage:
    110                 interpdata=GridToMesh(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False)
     110                interpdata=GridToMesh(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False)
    111111
    112112        Examples:
     
    114114        '''
    115115
    116         if npy.ndim(x)==2:
     116        if np.ndim(x)==2:
    117117                x=x.reshape(-1,)
    118         if npy.ndim(y)==2:
     118        if np.ndim(y)==2:
    119119                y=y.reshape(-1,)
    120120        if len(x) != data.shape[1]+1 and len(x) != data.shape[1]:
     
    134134        if len(x)==data.shape[1] and len(y)==data.shape[0]:
    135135                print '         x,y taken to define the center of data grid cells'
    136                 xind=npy.nonzero(npy.logical_and(x>xlim[0],x<xlim[1]))[0]
    137                 yind=npy.nonzero(npy.logical_and(y>ylim[0],y<ylim[1]))[0]
    138                 xg,yg=npy.meshgrid(x[xind],y[yind])
     136                xind=np.nonzero(np.logical_and(x>xlim[0],x<xlim[1]))[0]
     137                yind=np.nonzero(np.logical_and(y>ylim[0],y<ylim[1]))[0]
     138                xg,yg=np.meshgrid(x[xind],y[yind])
    139139                subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
    140140        elif len(x)==data.shape[1]+1 and len(y)==data.shape[0]+1:
    141141                print '         x,y taken to define the corners of data grid cells'
    142                 xcenter=npy.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),npy.float)
    143                 ycenter=npy.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),npy.float)
    144                 xind=npy.nonzero(npy.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
    145                 yind=npy.nonzero(npy.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
    146                 xg,yg=npy.meshgrid(xcenter[xind],ycenter[yind])
     142                xcenter=np.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),np.float)
     143                ycenter=np.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),np.float)
     144                xind=np.nonzero(np.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
     145                yind=np.nonzero(np.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
     146                xg,yg=np.meshgrid(xcenter[xind],ycenter[yind])
    147147                subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
    148148        else:
    149149                raise ValueError('x and y have inconsistent sizes: both should have length ncols(data)/nrows(data) or ncols(data)+1/nrows(data)+1')
    150150
    151         points=npy.array([xg.ravel(),yg.ravel()]).T
     151        points=np.array([xg.ravel(),yg.ravel()]).T
    152152        flatsubdata=subdata.ravel()
    153153
    154154        if plotonly:
    155                 plt.imshow(npy.flipud(subdata),origin='upper')
     155                plt.imshow(np.flipud(subdata),origin='upper')
    156156                plt.show()
    157157                return
    158158
    159159        # mask out any nan's in the data and corresponding coordinate points
    160         mask=npy.isnan(flatsubdata)
    161         ind=npy.nonzero(mask)[0]
     160        mask=np.isnan(flatsubdata)
     161        ind=np.nonzero(mask)[0]
    162162        if len(ind) and fill_nans:
    163163                print "         WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully."
    164         goodsubdata=npy.delete(flatsubdata,ind)
    165         goodpoints=npy.delete(points,ind,axis=0)
     164        goodsubdata=np.delete(flatsubdata,ind)
     165        goodpoints=np.delete(points,ind,axis=0)
    166166
    167167        # create spline and index spline at mesh points
     
    171171        if not fill_nans:
    172172                # identify nan's in xi,yi using nearest neighbors
    173                 xyinterp=npy.dstack([xi,yi])[0]
    174                 xydata=npy.dstack([xg.ravel(),yg.ravel()])[0]
     173                xyinterp=np.dstack([xi,yi])[0]
     174                xydata=np.dstack([xg.ravel(),yg.ravel()])[0]
    175175                tree=cKDTree(xydata)
    176176                nearest=tree.query(xyinterp)[1]
    177                 pos=npy.nonzero(npy.isnan(flatsubdata[nearest]))
     177                pos=np.nonzero(np.isnan(flatsubdata[nearest]))
    178178                interpdata[pos]=flatsubdata[nearest][pos]
    179179
  • issm/trunk-jpl/src/m/inversions/parametercontroldrag.py

    r16108 r21254  
    2626
    2727        #weights
    28         weights=options.getfieldvalue('weights',npy.ones(md.mesh.numberofvertices))
    29         if npy.size(weights)!=md.mesh.numberofvertices:
     28        weights=options.getfieldvalue('weights',np.ones(md.mesh.numberofvertices))
     29        if np.size(weights)!=md.mesh.numberofvertices:
    3030                md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices)
    3131        else:
     
    3434        #nsteps
    3535        nsteps=options.getfieldvalue('nsteps',100);
    36         if (npy.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):
     36        if (np.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):
    3737                md.inversion.nsteps=100
    3838        else:
     
    4141        #cm_min
    4242        cm_min=options.getfieldvalue('cm_min',ones(md.mesh.numberofvertices))
    43         if (npy.size(cm_min)==1):
     43        if (np.size(cm_min)==1):
    4444                md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices)
    45         elif (npy.size(cm_min)==md.mesh.numberofvertices):
     45        elif (np.size(cm_min)==md.mesh.numberofvertices):
    4646                md.inversion.min_parameters=cm_min
    4747        else:
     
    5050        #cm_max
    5151        cm_max=options.getfieldvalue('cm_max',250*ones(md.mesh.numberofvertices))
    52         if (npy.size(cm_max)==1):
     52        if (np.size(cm_max)==1):
    5353                md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices)
    54         elif (npy.size(cm_max)==md.mesh.numberofvertices):
     54        elif (np.size(cm_max)==md.mesh.numberofvertices):
    5555                md.inversion.max_parameters=cm_max
    5656        else:
     
    5959        #eps_cm
    6060        eps_cm=optoins.getfieldvalue('eps_cm',float('nan'))
    61         if (npy.size(eps_cm)~=1 | eps_cm<0 ):
     61        if (np.size(eps_cm)~=1 | eps_cm<0 ):
    6262                md.inversion.cost_function_threshold=float('nan')
    6363        else:
     
    6666        #maxiter
    6767        maxiter=options.getfieldvalue('maxiter',10*ones(md.inversion.nsteps))
    68         if (npy.any(maxiter<0) | npy.any(floor(maxiter)~=maxiter)):
     68        if (np.any(maxiter<0) | np.any(floor(maxiter)~=maxiter)):
    6969                md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps)
    7070        else:
     
    7575        #cm_jump
    7676        cm_jump=options.getfieldvalue('cm_jump',0.8*ones(md.inversion.nsteps))
    77         if !npy.isreal(cm_jump):
     77        if !np.isreal(cm_jump):
    7878                md.inversion.step_threshold=0.8*ones(md.inversion.nsteps)
    7979        else:
  • issm/trunk-jpl/src/m/materials/TMeltingPoint.py

    r17196 r21254  
    1 import numpy as npy
     1import numpy as np
    22
    33def TMeltingPoint(reftemp,pressure):
     
    1717
    1818        #ensure ref is same dimension as pressure
    19         ref=reftemp*npy.ones_like(pressure)
     19        ref=reftemp*np.ones_like(pressure)
    2020
    2121        return reftemp-beta*pressure
  • issm/trunk-jpl/src/m/mech/analyticaldamage.py

    r17965 r21254  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33#from plotmodel import plotmodel
     
    5454
    5555        if isinstance(sigmab,(int,float)):
    56                 sigmab=sigmab*npy.ones((md.mesh.numberofvertices,))
     56                sigmab=sigmab*np.ones((md.mesh.numberofvertices,))
    5757
    5858        # check inputs
     
    6161        if not '2d' in md.mesh.__doc__:
    6262                raise StandardError('only 2d (planview) model supported currently')
    63         if npy.any(md.flowequation.element_equation!=2):
     63        if np.any(md.flowequation.element_equation!=2):
    6464                print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
    6565
     
    7676        n=averaging(md,md.materials.rheology_n,0)
    7777       
    78         D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/npy.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/npy.sign(ex)
     78        D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/np.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/np.sign(ex)
    7979       
    8080        # D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
    81         pos=npy.nonzero(D>1)
     81        pos=np.nonzero(D>1)
    8282        D[pos]=0
    8383       
    84         backstress=npy.zeros((md.mesh.numberofvertices,))
     84        backstress=np.zeros((md.mesh.numberofvertices,))
    8585
    8686        # backstress to bring D down to one
    87         backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
     87        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np.sign(ex[pos])*(2.+a[pos])*np.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
    8888       
    89         pos=npy.nonzero(D<0)
     89        pos=np.nonzero(D<0)
    9090        #mask=ismember(1:md.mesh.numberofvertices,pos);
    9191        D[pos]=0
    9292       
    9393        # backstress to bring negative damage to zero
    94         backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
     94        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np.sign(ex[pos])*(2.+a[pos])*np.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
    9595       
    96         pos=npy.nonzero(backstress<0)
     96        pos=np.nonzero(backstress<0)
    9797        backstress[pos]=0
    9898       
    9999        # rigidity from Thomas relation for D=0 and backstress=0
    100         B=npy.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(npy.abs(ex)**(1./n))
    101         pos=npy.nonzero(B<0)
     100        B=np.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(np.abs(ex)**(1./n))
     101        pos=np.nonzero(B<0)
    102102        B[pos]=md.materials.rheology_B[pos]
    103103       
  • issm/trunk-jpl/src/m/mech/backstressfrominversion.py

    r18011 r21254  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33from thomasparams import thomasparams
     
    6565        if tempmask:
    6666                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
    67                 pos=npy.nonzero(Bi>md.materials.rheology_B)
     67                pos=np.nonzero(Bi>md.materials.rheology_B)
    6868                Bi[pos]=md.materials.rheology_B[pos]
    6969       
    7070        # analytical backstress solution
    71         backstress=T-Bi*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
    72         backstress[npy.nonzero(backstress<0)]=0
     71        backstress=T-Bi*np.sign(ex0)*(2+a0)*np.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
     72        backstress[np.nonzero(backstress<0)]=0
    7373
    7474        return backstress
  • issm/trunk-jpl/src/m/mech/calcbackstress.py

    r18011 r21254  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33from thomasparams import thomasparams
     
    6161       
    6262        # analytical backstress solution
    63         backstress=T-(1.-D)*B*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
    64         backstress[npy.nonzero(backstress<0)]=0
     63        backstress=T-(1.-D)*B*np.sign(ex0)*(2+a0)*np.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
     64        backstress[np.nonzero(backstress<0)]=0
    6565
    6666        return backstress
  • issm/trunk-jpl/src/m/mech/damagefrominversion.py

    r17624 r21254  
    1 import numpy as npy
     1import numpy as np
    22
    33def damagefrominversion(md):
     
    2424        if any(md.flowequation.element_equation!=2):
    2525                raise StandardError('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
    26         if npy.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
     26        if np.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
    2727                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1,)
    2828        else:
    2929                Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
    30         if npy.ndim(md.materials.rheology_B)==2:
     30        if np.ndim(md.materials.rheology_B)==2:
    3131                BT=md.materials.rheology_B.reshape(-1,)
    3232        else:
    3333                BT=md.materials.rheology_B
    3434
    35         damage=npy.zeros_like(Bi)
     35        damage=np.zeros_like(Bi)
    3636
    3737        # Damage where Bi softer than B(T)
    38         pos=npy.nonzero(Bi<BT)[0]
     38        pos=np.nonzero(Bi<BT)[0]
    3939        damage[pos]=1.-Bi[pos]/BT[pos]
    4040       
    41         pos=npy.nonzero(damage<0)
     41        pos=np.nonzero(damage<0)
    4242        damage[pos]=0
    4343
  • issm/trunk-jpl/src/m/mech/mechanicalproperties.py

    r19466 r21254  
    1 import numpy as npy
     1import numpy as np
    22from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff
    33from results import results
     
    2727        #       raise StandardError('only 2D model supported currently')
    2828
    29         if npy.any(md.flowequation.element_equation!=2):
     29        if np.any(md.flowequation.element_equation!=2):
    3030                print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
    3131
     
    3535            if len(damage)!=md.mesh.numberofvertices:
    3636                raise ValueError('if damage is supplied it should be of size ' + md.mesh.numberofvertices)
    37             if npy.ndim(damage)==2:
     37            if np.ndim(damage)==2:
    3838                damage=damage.reshape(-1,)
    3939        else: damage=None
    4040
    41         if npy.ndim(vx)==2:
     41        if np.ndim(vx)==2:
    4242                vx=vx.reshape(-1,)
    43         if npy.ndim(vy)==2:
     43        if np.ndim(vy)==2:
    4444                vy=vy.reshape(-1,)
    4545       
     
    4848        numberofvertices=md.mesh.numberofvertices
    4949        index=md.mesh.elements
    50         summation=npy.array([[1],[1],[1]])
    51         directionsstress=npy.zeros((numberofelements,4))
    52         directionsstrain=npy.zeros((numberofelements,4))
    53         valuesstress=npy.zeros((numberofelements,2))
    54         valuesstrain=npy.zeros((numberofelements,2))
     50        summation=np.array([[1],[1],[1]])
     51        directionsstress=np.zeros((numberofelements,4))
     52        directionsstrain=np.zeros((numberofelements,4))
     53        valuesstress=np.zeros((numberofelements,2))
     54        valuesstrain=np.zeros((numberofelements,2))
    5555       
    5656        #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
     
    6060        vxlist=vx[index-1]/md.constants.yts
    6161        vylist=vy[index-1]/md.constants.yts
    62         ux=npy.dot((vxlist*alpha),summation).reshape(-1,)
    63         uy=npy.dot((vxlist*beta),summation).reshape(-1,)
    64         vx=npy.dot((vylist*alpha),summation).reshape(-1,)
    65         vy=npy.dot((vylist*beta),summation).reshape(-1,)
     62        ux=np.dot((vxlist*alpha),summation).reshape(-1,)
     63        uy=np.dot((vxlist*beta),summation).reshape(-1,)
     64        vx=np.dot((vylist*alpha),summation).reshape(-1,)
     65        vy=np.dot((vylist*beta),summation).reshape(-1,)
    6666        uyvx=(vx+uy)/2.
    6767        #clear vxlist vylist
    6868       
    6969        #compute viscosity
    70         nu=npy.zeros((numberofelements,))
    71         B_bar=npy.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
     70        nu=np.zeros((numberofelements,))
     71        B_bar=np.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
    7272        power=((md.materials.rheology_n-1.)/(2.*md.materials.rheology_n)).reshape(-1,)
    7373        second_inv=(ux**2.+vy**2.+((uy+vx)**2.)/4.+ux*vy).reshape(-1,)
    7474       
    7575        #some corrections
    76         location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
     76        location=np.nonzero(np.logical_and(second_inv==0,power!=0))
    7777        nu[location]=10^18      #arbitrary maximum viscosity to apply where there is no effective shear
    7878       
    7979        if 'matice' in md.materials.__module__:
    80                 location=npy.nonzero(second_inv)
     80                location=np.nonzero(second_inv)
    8181                nu[location]=B_bar[location]/(second_inv[location]**power[location])
    82                 location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
     82                location=np.nonzero(np.logical_and(second_inv==0,power==0))
    8383                nu[location]=B_bar[location]
    84                 location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
     84                location=np.nonzero(np.logical_and(second_inv==0,power!=0))
    8585                nu[location]=10^18
    8686        elif 'matdamageice' in md.materials.__module__ and damage is not None:
    8787                print 'computing damage-dependent properties!'
    88                 Zinv=npy.dot(1-damage[index-1],summation/3.).reshape(-1,)
    89                 location=npy.nonzero(second_inv)
    90                 nu[location]=Zinv[location]*B_bar[location]/npy.power(second_inv[location],power[location])
    91                 location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
     88                Zinv=np.dot(1-damage[index-1],summation/3.).reshape(-1,)
     89                location=np.nonzero(second_inv)
     90                nu[location]=Zinv[location]*B_bar[location]/np.power(second_inv[location],power[location])
     91                location=np.nonzero(np.logical_and(second_inv==0,power==0))
    9292                nu[location]=Zinv[location]*B_bar[location]
    9393                #clear Zinv
     
    101101       
    102102        #compute principal properties of stress
    103         for i in npy.arange(numberofelements):
     103        for i in np.arange(numberofelements):
    104104       
    105105                #compute stress and strainrate matrices
    106                 stress=npy.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
    107                 strain=npy.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
     106                stress=np.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
     107                strain=np.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
    108108       
    109109                #eigenvalues and vectors for stress
    110                 value,directions=npy.linalg.eig(stress);
     110                value,directions=np.linalg.eig(stress);
    111111                idx=abs(value).argsort()[::-1] # sort in descending order
    112112                value=value[idx]
     
    116116
    117117                #eigenvalues and vectors for strain
    118                 value,directions=npy.linalg.eig(strain);
     118                value,directions=np.linalg.eig(strain);
    119119                idx=abs(value).argsort()[::-1] # sort in descending order
    120120                value=value[idx]
     
    133133        stress.principalvalue2=valuesstress[:,1]
    134134        stress.principalaxis2=directionsstress[:,2:4]
    135         stress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
     135        stress.effectivevalue=1./np.sqrt(2.)*np.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
    136136        md.results.stress=stress
    137137       
     
    144144        strainrate.principalvalue2=valuesstrain[:,1]*md.constants.yts
    145145        strainrate.principalaxis2=directionsstrain[:,2:4]
    146         strainrate.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
     146        strainrate.effectivevalue=1./np.sqrt(2.)*np.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
    147147        md.results.strainrate=strainrate
    148148       
     
    155155        deviatoricstress.principalvalue2=valuesstress[:,1]
    156156        deviatoricstress.principalaxis2=directionsstress[:,2:4]
    157         deviatoricstress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
     157        deviatoricstress.effectivevalue=1./np.sqrt(2.)*np.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
    158158        md.results.deviatoricstress=deviatoricstress
    159159
  • issm/trunk-jpl/src/m/mech/robintemperature.py

    r17836 r21254  
    1 import numpy as npy
     1import numpy as np
    22from scipy.special import erf
    33
     
    3434       
    3535        #create vertical coordinate variable
    36         zstar=npy.sqrt(2.*alphaT*thickness/accumrate)
     36        zstar=np.sqrt(2.*alphaT*thickness/accumrate)
    3737       
    38         tprofile=surftemp+npy.sqrt(2.*thickness*npy.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
     38        tprofile=surftemp+np.sqrt(2.*thickness*np.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
    3939
    4040        return tprofile
  • issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.py

    r18068 r21254  
    1 import numpy as npy
     1import numpy as np
    22
    33def steadystateiceshelftemp(md,surfacetemp,basaltemp):
     
    4242        temperature=(Ts+Tb)/2  # where wi~=0
    4343       
    44         pos=npy.nonzero(abs(wi)>=1e-4) # to avoid division by zero
     44        pos=np.nonzero(abs(wi)>=1e-4) # to avoid division by zero
    4545
    46         npy.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
     46        np.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
    4747        #calculate depth-averaged temperature (in Celsius)
    4848        try:
    49                 temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*npy.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(npy.exp(Hi[pos]*wi[pos]/ki)-1))
     49                temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*np.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(np.exp(Hi[pos]*wi[pos]/ki)-1))
    5050        except FloatingPointError:
    5151                print 'steadystateiceshelf warning: overflow encountered in multipy/divide/exp, trying another formulation.'
    52                 temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/npy.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-npy.exp(-Hi[pos]*wi[pos]/ki)))
     52                temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/np.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-np.exp(-Hi[pos]*wi[pos]/ki)))
    5353       
    5454        #temperature should not be less than surface temp
    55         pos=npy.nonzero(temperature<Ts)
     55        pos=np.nonzero(temperature<Ts)
    5656        temperature[pos]=Ts[pos]
    5757       
    5858        # NaN where melt rates are too high (infinity/infinity in exponential)
    59         pos=npy.nonzero(npy.isnan(temperature))
     59        pos=np.nonzero(np.isnan(temperature))
    6060        temperature[pos]=Ts[pos]
    6161       
  • issm/trunk-jpl/src/m/mech/thomasparams.py

    r17903 r21254  
    1 import numpy as npy
     1import numpy as np
    22from averaging import averaging
    33
     
    7575       
    7676        # checks: any of e1 or e2 equal to zero?
    77         pos=npy.nonzero(e1==0)
    78         if npy.any(pos==1):
     77        pos=np.nonzero(e1==0)
     78        if np.any(pos==1):
    7979                print 'WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1'
    8080                e1[pos]=1.e-13
    81         pos=npy.nonzero(e2==0)
    82         if npy.any(pos==1):
     81        pos=np.nonzero(e2==0)
     82        if np.any(pos==1):
    8383                print 'WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1'
    8484                e2[pos]=1.e-13
     
    8989       
    9090        if coordsys=='principal':
    91                 b=npy.zeros((md.mesh.numberofvertices,))
     91                b=np.zeros((md.mesh.numberofvertices,))
    9292                ex=e1
    9393                a=e2/e1
    94                 pos=npy.nonzero(npy.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
     94                pos=np.nonzero(np.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
    9595                a[pos]=e1[pos]/e2[pos]
    9696                ex[pos]=e2[pos]
    97                 pos2=npy.nonzero(e1<0 & e2<0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal compression
     97                pos2=np.nonzero(e1<0 & e2<0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal compression
    9898                a[pos2]=e1[pos2]/e2[pos2]
    9999                ex[pos2]=e2[pos2]
    100                 pos3=npy.nonzero(e1>0 & e2>0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal tension
     100                pos3=np.nonzero(e1>0 & e2>0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal tension
    101101                a[pos3]=e1[pos3]/e2[pos3]
    102102                ex[pos3]=e2[pos3]
    103                 ind=npy.nonzero(e1<0 & e2<0)
     103                ind=np.nonzero(e1<0 & e2<0)
    104104                a[ind]=-a[ind] # where both strain rates are compressive, enforce negative alpha
    105                 sigxx=(npy.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
     105                sigxx=(np.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
    106106        elif coordsys=='xy':
    107107                ex=exx
     
    110110        elif coordsys=='longitudinal':
    111111                # using longitudinal strain rates defined by observed velocity vector
    112                 velangle=npy.arctan(md.initialization.vy/md.initialization.vx)
    113                 pos=npy.nonzero(md.initialization.vx==0)
    114                 velangle[pos]=npy.pi/2
    115                 ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
     112                velangle=np.arctan(md.initialization.vy/md.initialization.vx)
     113                pos=np.nonzero(md.initialization.vx==0)
     114                velangle[pos]=np.pi/2
     115                ex=0.5*(exx+eyy)+0.5*(exx-eyy)*np.cos(2.*velangle)+exy*np.sin(2.*velangle)
    116116                ey=exx+eyy-ex # trace of strain rate tensor is invariant
    117                 exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
     117                exy=-0.5*(exx-eyy)*np.sin(2.*velangle)+exy*np.cos(2.*velangle)
    118118                a=ey/ex
    119119                b=exy/ex
     
    124124        # a < -1 in areas of strong lateral compression or longitudinal compression and
    125125        # theta flips sign at a = -2
    126         pos=npy.nonzero(npy.abs((npy.abs(a)-2.))<1.e-3)
     126        pos=np.nonzero(np.abs((np.abs(a)-2.))<1.e-3)
    127127        if len(pos)>0:
    128128                print 'Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2'
     
    131131        if eq=='Weertman1D':
    132132                theta=1./8
    133                 a=npy.zeros((md.mesh.numberofvertices,))
     133                a=np.zeros((md.mesh.numberofvertices,))
    134134        elif eq=='Weertman2D':
    135135                theta=1./9
    136                 a=npy.ones((md.mesh.numberofvertices,))
     136                a=np.ones((md.mesh.numberofvertices,))
    137137        elif eq=='Thomas':
    138                 theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(npy.abs(2.+a)**n)
     138                theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(np.abs(2.+a)**n)
    139139        else:
    140140                raise ValueError('argument passed to "eq" not valid')
  • issm/trunk-jpl/src/m/plot/applyoptions.py

    r21253 r21254  
    1 import numpy as npy
     1import numpy as np
    22from cmaptools import truncate_colormap
    33from plot_contour import plot_contour
     
    223223                        cb.locator=MaxNLocator(nbins=5) # default 5 ticks
    224224                if options.exist('colorbartickspacing'):
    225                         locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))
     225                        locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))
    226226                        cb.set_ticks(locs)
    227227                if options.exist('colorbarlines'):
    228                         locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))
    229                         cb.add_lines(locs,['k' for i in range(len(locs))],npy.ones_like(locs))
     228                        locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))
     229                        cb.add_lines(locs,['k' for i in range(len(locs))],np.ones_like(locs))
    230230                if options.exist('colorbarlineatvalue'):
    231231                        locs=options.getfieldvalue('colorbarlineatvalue')
    232232                        colors=options.getfieldvalue('colorbarlineatvaluecolor',['k' for i in range (len(locs))])
    233                         widths=options.getfieldvalue('colorbarlineatvaluewidth',npy.ones_like(locs))
     233                        widths=options.getfieldvalue('colorbarlineatvaluewidth',np.ones_like(locs))
    234234                        cb.add_lines(locs,colors,widths)
    235235                if options.exist('colorbartitle'):
  • issm/trunk-jpl/src/m/plot/colormaps/cmaptools.py

    r17649 r21254  
    1 import numpy as npy
     1import numpy as np
    22
    33try:
     
    2121
    2222        new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name,
    23                 a=minval, b=maxval), cmap(npy.linspace(minval, maxval, n)))
     23                a=minval, b=maxval), cmap(np.linspace(minval, maxval, n)))
    2424       
    2525        return new_cmap
  • issm/trunk-jpl/src/m/plot/plot_contour.py

    r20921 r21254  
    2626                pass
    2727        elif datatype==3: # quiver (vector) data
    28                 data=npy.sqrt(datain**2)
     28                data=np.sqrt(datain**2)
    2929        else:
    3030                raise ValueError('datatype not supported in call to plot_contour')
  • issm/trunk-jpl/src/m/plot/plot_overlay.py

    r19432 r21254  
    1 import numpy as npy
     1import numpy as np
    22from processmesh import processmesh
    33from processdata import processdata
     
    2828        if data=='none' or data==None:
    2929                imageonly=1
    30                 data=npy.float('nan')*npy.ones((md.mesh.numberofvertices,))
     30                data=np.float('nan')*np.ones((md.mesh.numberofvertices,))
    3131                datatype=1
    3232        else:
     
    7070
    7171        # normalize array
    72         arr=arr/npy.float(npy.max(arr.ravel()))
     72        arr=arr/np.float(np.max(arr.ravel()))
    7373        arr=1.-arr # somehow the values got flipped
    7474
     
    9696        dy=trans[5]     
    9797       
    98         xarr=npy.arange(xmin,xmax,dx)
    99         yarr=npy.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
    100         xg,yg=npy.meshgrid(xarr,yarr)
     98        xarr=np.arange(xmin,xmax,dx)
     99        yarr=np.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
     100        xg,yg=np.meshgrid(xarr,yarr)
    101101        overlaylims=options.getfieldvalue('overlaylims',[min(arr.ravel()),max(arr.ravel())])
    102102        norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1])
     
    118118                #test
    119119                #m.ax=ax
    120                 meridians=npy.arange(-180.,181.,1.)
    121                 parallels=npy.arange(-80.,80.,1.)
     120                meridians=np.arange(-180.,181.,1.)
     121                parallels=np.arange(-80.,80.,1.)
    122122                m.drawparallels(parallels,labels=[0,0,1,1]) # labels=[left,right,top,bottom]
    123123                m.drawmeridians(meridians,labels=[1,1,0,0])
    124124                m.drawcoastlines()
    125                 pc=m.pcolormesh(xg, yg, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax)
     125                pc=m.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax)
    126126
    127127        else:
    128                 pc=ax.pcolormesh(xg, yg, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
     128                pc=ax.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
    129129       
    130130        #rasterization?
  • issm/trunk-jpl/src/m/plot/plot_streamlines.py

    r20144 r21254  
    1 import numpy as npy
     1import numpy as np
    22from processmesh import processmesh
    33from processdata import processdata
     
    4545
    4646    # format data for matplotlib streamplot function
    47     yg,xg=npy.mgrid[min(md.mesh.y):max(md.mesh.y):100j,min(md.mesh.x):max(md.mesh.x):100j]
     47    yg,xg=np.mgrid[min(md.mesh.y):max(md.mesh.y):100j,min(md.mesh.x):max(md.mesh.x):100j]
    4848    ug=griddata((x,y),u,(xg,yg),method='linear')
    4949    vg=griddata((x,y),v,(xg,yg),method='linear')
     
    6060    if linewidth=='vel':
    6161        scale=options.getfieldvalue('streamlineswidthscale',3)
    62         vel=npy.sqrt(ug**2+vg**2)
    63         linewidth=scale*vel/npy.amax(vel)
     62        vel=np.sqrt(ug**2+vg**2)
     63        linewidth=scale*vel/np.amax(vel)
    6464        linewidth[linewidth<0.5]=0.5
    6565
  • issm/trunk-jpl/src/m/plot/plotmodel.py

    r21253 r21254  
    1 import numpy as npy
     1import numpy as np
    22from plotoptions import plotoptions
    33
     
    3434                nr=True
    3535        else:
    36                 nrows=npy.ceil(numberofplots/subplotwidth)
     36                nrows=np.ceil(numberofplots/subplotwidth)
    3737                nr=False
    3838       
  • issm/trunk-jpl/src/m/plot/processmesh.py

    r21253 r21254  
    11from math import isnan
    2 import numpy as npy
     2import numpy as np
    33
    44def processmesh(md,data,options):
     
    3535                        z=md.mesh.z
    3636                except AttributeError:
    37                         z=npy.zeros_like(md.mesh.x)
     37                        z=np.zeros_like(md.mesh.x)
    3838               
    3939                try:
  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

    r20900 r21254  
    33from collections import OrderedDict
    44import results as resultsclass
    5 import MatlabFuncs as m
    65
    76def parseresultsfromdisk(md,filename,iosplit):
    87        if iosplit:
    9                 results=parseresultsfromdiskiosplit(md,filename)
     8                saveres=parseresultsfromdiskiosplit(md,filename)
    109        else:
    11                 results=parseresultsfromdiskioserial(md,filename)
    12 
    13         return results
     10                saveres=parseresultsfromdiskioserial(md,filename)
     11
     12        return saveres
    1413
    1514def parseresultsfromdiskioserial(md,filename):    # {{{
    16 
    1715        #Open file
    1816        try:
     
    2220
    2321        #initialize results:
    24         results=[]
    25         results.append(None)
     22        saveres=[]
    2623
    2724        #Read fields until the end of the file.
    28         result=ReadData(fid,md)
     25        loadres=ReadData(fid,md)
    2926
    3027        counter=0
    3128        check_nomoresteps=0
    32         step=result['step']
    33 
    34         while result:
    35 
     29        step=loadres['step']
     30        print ('step is now {}'.format(step))
     31
     32        while loadres:
     33                #check that the new result does not add a step, which would be an error:
    3634                if check_nomoresteps:
    37                         #check that the new result does not add a step, which would be an error:
    38                         if result['step']>=1:
     35                        if loadres['step']>=1:
    3936                                raise TypeError("parsing results for a steady-state core, which incorporates transient results!")
    4037
    4138                #Check step, increase counter if this is a new step
    42                 if(step!=result['step'] and result['step']>1):
     39                if(step!=loadres['step'] and loadres['step']>1):
    4340                        counter = counter + 1
    44                         step    = result['step']
    45 
    46                 #Add result
    47                 if result['step']==0:
     41                        step    = loadres['step']
     42
     43                #Add result
     44                if loadres['step']==0:
    4845                        #if we have a step = 0, this is a steady state solution, don't expect more steps.
    4946                        index = 0;
    5047                        check_nomoresteps=1
    51        
    52                 elif result['step']==1:
     48                elif loadres['step']==1:
    5349                        index = 0
    5450                else:
    5551                        index = counter;
    56        
    57                 if index > len(results)-1:
    58                         for i in xrange(len(results)-1,index-1):
    59                                 results.append(None)
    60                         results.append(resultsclass.results())
    6152               
    62                 elif results[index] is None:
    63                         results[index]=resultsclass.results()
    64 
     53                if index > len(saveres)-1:
     54                        for i in xrange(len(saveres)-1,index-1):
     55                                saveres.append(None)
     56                        saveres.append(resultsclass.results())
     57                elif saveres[index] is None:
     58                        saveres[index]=resultsclass.results()
    6559                       
    6660                #Get time and step
    67                 if result['step'] != -9999.:
    68                         setattr(results[index],'step',result['step'])
    69                 if result['time'] != -9999.:
    70                         setattr(results[index],'time',result['time'])
    71        
    72                 #Add result
    73                 if hasattr(results[index],result['fieldname']) and not m.strcmp(result['fieldname'],'SolutionType'):
    74                         setattr(results[index],result['fieldname'],numpy.vstack((getattr(results[index],result['fieldname']),result['field'])))
    75                 else:
    76                         setattr(results[index],result['fieldname'],result['field'])
    77 
    78                 #read next result
    79                 result=ReadData(fid,md)
     61                if loadres['step'] != -9999.:
     62                        saveres[index].__dict__['step']=loadres['step']
     63                if loadres['time'] != -9999.:
     64                        saveres[index].__dict__['time']=loadres['time']
     65
     66                #Add result
     67                saveres[index].__dict__[loadres['fieldname']]=loadres['field']
     68
     69                #read next result
     70                loadres=ReadData(fid,md)
    8071
    8172        fid.close()
    8273
    83         return results
     74        return saveres
    8475        # }}}
    8576def parseresultsfromdiskiosplit(md,filename):    # {{{
     
    9182                raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename)
    9283
    93         results=[]
     84        saveres=[]
    9485
    9586        #if we have done split I/O, ie, we have results that are fragmented across patches,
    9687        #do a first pass, and figure out the structure of results
    97         result=ReadDataDimensions(fid)
    98         while result:
     88        loadres=ReadDataDimensions(fid)
     89        while loadres:
    9990
    10091                #Get time and step
    101                 if result['step'] > len(results):
    102                         for i in xrange(len(results),result['step']-1):
    103                                 results.append(None)
    104                         results.append(resultsclass.results())
    105                 setattr(results[result['step']-1],'step',result['step'])
    106                 setattr(results[result['step']-1],'time',result['time'])
    107 
    108                 #Add result
    109                 setattr(results[result['step']-1],result['fieldname'],float('NaN'))
    110 
    111                 #read next result
    112                 result=ReadDataDimensions(fid)
     92                if loadres['step'] > len(saveres):
     93                        for i in xrange(len(saveres),loadres['step']-1):
     94                                saveres.append(None)
     95                        saveres.append(resultsclass.results())
     96                setattr(saveres[loadres['step']-1],'step',loadres['step'])
     97                setattr(saveres[loadres['step']-1],'time',loadres['time'])
     98
     99                #Add result
     100                setattr(saveres[loadres['step']-1],loadres['fieldname'],float('NaN'))
     101
     102                #read next result
     103                loadres=ReadDataDimensions(fid)
    113104
    114105        #do a second pass, and figure out the size of the patches
    115106        fid.seek(0)    #rewind
    116         result=ReadDataDimensions(fid)
    117         while result:
    118 
    119                 #read next result
    120                 result=ReadDataDimensions(fid)
     107        loadres=ReadDataDimensions(fid)
     108        while loadres:
     109
     110                #read next result
     111                loadres=ReadDataDimensions(fid)
    121112
    122113        #third pass, this time to read the real information
    123114        fid.seek(0)    #rewind
    124         result=ReadData(fid,md)
    125         while result:
     115        loadres=ReadData(fid,md)
     116        while loadres:
    126117
    127118                #Get time and step
    128                 if result['step']> len(results):
    129                         for i in xrange(len(results),result['step']-1):
    130                                 results.append(None)
    131                         results.append(resultsclass.results())
    132                 setattr(results[result['step']-1],'step',result['step'])
    133                 setattr(results[result['step']-1],'time',result['time'])
    134 
    135                 #Add result
    136                 setattr(results[result['step']-1],result['fieldname'],result['field'])
    137 
    138                 #read next result
    139                 result=ReadData(fid,md)
     119                if loadres['step']> len(saveres):
     120                        for i in xrange(len(saveres),loadres['step']-1):
     121                                saveres.append(None)
     122                        saveres.append(saveresclass.saveres())
     123                setattr(saveres[loadres['step']-1],'step',loadres['step'])
     124                setattr(saveres[loadres['step']-1],'time',loadres['time'])
     125
     126                #Add result
     127                setattr(saveres[loadres['step']-1],loadres['fieldname'],loadres['field'])
     128
     129                #read next result
     130                loadres=ReadData(fid,md)
    140131
    141132        #close file
    142133        fid.close()
    143134
    144         return results
     135        return saveres
    145136        # }}}
    146137def ReadData(fid,md):    # {{{
     
    177168                #Process units here FIXME: this should not be done here!
    178169                yts=md.constants.yts
    179                 if m.strcmp(fieldname,'BalancethicknessThickeningRate'):
    180                         field = field*yts
    181                 elif m.strcmp(fieldname,'Time'):
     170                if fieldname=='BalancethicknessThickeningRate':
     171                        field = field*yts
     172                elif fieldname=='Time':
    182173                        field = field/yts
    183                 elif m.strcmp(fieldname,'HydrologyWaterVx'):
    184                         field = field*yts
    185                 elif m.strcmp(fieldname,'HydrologyWaterVy'):
    186                         field = field*yts
    187                 elif m.strcmp(fieldname,'Vx'):
    188                         field = field*yts
    189                 elif m.strcmp(fieldname,'Vy'):
    190                         field = field*yts
    191                 elif m.strcmp(fieldname,'Vz'):
    192                         field = field*yts
    193                 elif m.strcmp(fieldname,'Vel'):
    194                         field = field*yts
    195                 elif m.strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'):
    196                         field = field*yts
    197                 elif m.strcmp(fieldname,'TotalFloatingBmb'):
     174                elif fieldname=='HydrologyWaterVx':
     175                        field = field*yts
     176                elif fieldname=='HydrologyWaterVy':
     177                        field = field*yts
     178                elif fieldname=='Vx':
     179                        field = field*yts
     180                elif fieldname=='Vy':
     181                        field = field*yts
     182                elif fieldname=='Vz':
     183                        field = field*yts
     184                elif fieldname=='Vel':
     185                        field = field*yts
     186                elif fieldname=='BasalforcingsGroundediceMeltingRate':
     187                        field = field*yts
     188                elif fieldname=='TotalFloatingBmb':
    198189                        field = field/10.**12.*yts #(GigaTon/year)
    199                 elif m.strcmp(fieldname,'TotalGroundedBmb'):
     190                elif fieldname=='TotalGroundedBmb':
    200191                        field = field/10.**12.*yts #(GigaTon/year)
    201                 elif m.strcmp(fieldname,'TotalSmb'):
     192                elif fieldname=='TotalSmb':
    202193                        field = field/10.**12.*yts #(GigaTon/year)
    203                 elif m.strcmp(fieldname,'SmbMassBalance'):
    204                         field = field*yts
    205                 elif m.strcmp(fieldname,'CalvingCalvingrate'):
    206                         field = field*yts
    207 
    208                 result=OrderedDict()
    209                 result['fieldname']=fieldname
    210                 result['time']=time
    211                 result['step']=step
    212                 result['field']=field
     194                elif fieldname=='SmbMassBalance':
     195                        field = field*yts
     196                elif fieldname=='CalvingCalvingrate':
     197                        field = field*yts
     198
     199                saveres=OrderedDict()
     200                saveres['fieldname']=fieldname
     201                saveres['time']=time
     202                saveres['step']=step
     203                saveres['field']=field
    213204
    214205        except struct.error as e:
    215                 result=None
    216 
    217         return result
     206                saveres=None
     207
     208        return saveres
    218209        # }}}
    219210def ReadDataDimensions(fid):    # {{{
     
    246237                        raise TypeError("cannot read data of type %d" % type)
    247238
    248                 result=OrderedDict()
    249                 result['fieldname']=fieldname
    250                 result['time']=time
    251                 result['step']=step
    252                 result['M']=M
    253                 result['N']=N
     239                saveres=OrderedDict()
     240                saveres['fieldname']=fieldname
     241                saveres['time']=time
     242                saveres['step']=step
     243                saveres['M']=M
     244                saveres['N']=N
    254245
    255246        except struct.error as e:
    256                 result=None
    257 
    258         return result
    259         # }}}
     247                saveres=None
     248
     249        return saveres
     250        # }}}
Note: See TracChangeset for help on using the changeset viewer.