Changeset 21254
- Timestamp:
- 10/11/16 01:24:07 (8 years ago)
- Location:
- issm/trunk-jpl/src/m
- Files:
-
- 33 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py
r17567 r21254 1 1 import os 2 import numpy as np y2 import numpy as np 3 3 def PattynSMB(md,Tf): 4 4 """ … … 25 25 26 26 #Double check lat and long exist: 27 if np y.any(npy.isnan(md.mesh.lat)):27 if np.any(np.isnan(md.mesh.lat)): 28 28 raise IOError('PattynSMB error message: md.mesh.lat field required') 29 29 … … 36 36 # Calculate summer temperature, Eqn (12) 37 37 # 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*np y.abs(md.mesh.lat) + Tf38 Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*np.abs(md.mesh.lat) + Tf 39 39 Tms= Tms[0] 40 40 … … 43 43 44 44 # Calculate Ablation, Eqn (10) (use for both Antarctica & Greenland), max melt is 10m/a 45 ABLdot=0.*np y.ones(md.mesh.numberofvertices)46 pos=np y.nonzero(Tms>=0)47 ABLdot[pos]=np y.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) 48 48 49 49 smb=ACCdot-ABLdot -
issm/trunk-jpl/src/m/classes/outputdefinition.py
r21049 r21254 2 2 from checkfield import checkfield 3 3 from WriteData import WriteData 4 import numpy as np y4 import numpy as np 5 5 6 6 class outputdefinition(object): … … 40 40 data.append(classdefinition) 41 41 42 data=np y.unique(data);42 data=np.unique(data); 43 43 WriteData(fid,prefix,'data',data,'name','md.outputdefinition.list','format','StringArray'); 44 44 # }}} -
issm/trunk-jpl/src/m/classes/qmu.py
r21049 r21254 5 5 from checkfield import checkfield 6 6 from WriteData import WriteData 7 import MatlabFuncs as m8 7 9 8 class qmu(object): … … 41 40 42 41 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) 46 44 fnames=vars(variable) 47 45 maxlen=0 … … 53 51 (maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname))) 54 52 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) 58 55 fnames=vars(response) 59 56 maxlen=0 … … 67 64 s+="%s\n" % fielddisplay(self,'numberofresponses','number of responses') 68 65 69 for i,method in enumerate(self.method.iteritems()):66 for method in self.method: 70 67 if isinstance(method,'dakota_method'): 71 68 s+=" method%s : '%s'\n" % \ 72 ( string_dim(method,i),method.method)69 (len(method),method.method) 73 70 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) 77 73 fnames=vars(param) 78 74 maxlen=0 … … 84 80 (maxlen+1,fname,any2str(getattr(param,fname))) 85 81 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) 89 84 fnames=vars(result) 90 85 maxlen=0 … … 132 127 md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas") 133 128 134 if not m.strcmpi(md.cluster.name,'none'):129 if md.cluster.name!='none': 135 130 if not md.settings.waitonlock: 136 131 md.checkmessage("waitonlock should be activated when running qmu in parallel mode!") -
issm/trunk-jpl/src/m/classes/settings.py
r21049 r21254 16 16 self.lowmem = 0 17 17 self.output_frequency = 0 18 self.recording_frequency 18 self.recording_frequency = 0 19 19 self.waitonlock = 0 20 20 -
issm/trunk-jpl/src/m/coordsystems/ll2xy.py
r17712 r21254 1 import numpy as np y1 import numpy as np 2 2 3 3 def ll2xy(lat,lon,sgn=-1,central_meridian=0,standard_parallel=71): … … 36 36 ex2 = .006693883 37 37 # Eccentricity of the Hughes ellipsoid 38 ex = np y.sqrt(ex2)38 ex = np.sqrt(ex2) 39 39 40 latitude = np y.abs(lat) * npy.pi/180.41 longitude = (lon + delta) * np y.pi/180.40 latitude = np.abs(lat) * np.pi/180. 41 longitude = (lon + delta) * np.pi/180. 42 42 43 43 # compute X and Y in grid coordinates. 44 T = np y.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) 45 45 46 46 if (90 - slat) < 1.e-5: 47 rho = 2.*re*T/np y.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))47 rho = 2.*re*T/np.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex)) 48 48 else: 49 sl = slat*np y.pi/180.50 tc = np y.tan(npy.pi/4.-sl/2.)/((1.-ex*npy.sin(sl))/(1.+ex*npy.sin(sl)))**(ex/2.)51 mc = np y.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)) 52 52 rho = re*mc*T/tc 53 53 54 y = -rho * sgn * np y.cos(sgn*longitude)55 x = rho * sgn * np y.sin(sgn*longitude)54 y = -rho * sgn * np.cos(sgn*longitude) 55 x = rho * sgn * np.sin(sgn*longitude) 56 56 57 cnt1=np y.nonzero(latitude>= npy.pi/2.)[0]57 cnt1=np.nonzero(latitude>= np.pi/2.)[0] 58 58 59 59 if cnt1: -
issm/trunk-jpl/src/m/coordsystems/xy2ll.py
r17712 r21254 1 import numpy as np y1 import numpy as np 2 2 from math import pi 3 3 … … 39 39 # if x,y passed as lists, convert to numpy arrays 40 40 if type(x) != "numpy.ndarray": 41 x=np y.array(x)41 x=np.array(x) 42 42 if type(y) != "numpy.ndarray": 43 y=np y.array(y)43 y=np.array(y) 44 44 45 45 ## Conversion constant from degrees to radians … … 50 50 ex2 = .006693883 51 51 ## Eccentricity of the Hughes ellipsoid 52 ex = np y.sqrt(ex2)52 ex = np.sqrt(ex2) 53 53 54 54 sl = slat*pi/180. 55 rho = np y.sqrt(x**2 + y**2)56 cm = np y.cos(sl) / npy.sqrt(1.0 - ex2 * (npy.sin(sl)**2))57 T = np y.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) 58 58 59 59 if abs(slat-90.) < 1.e-5: 60 T = rho*np y.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re60 T = rho*np.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re 61 61 else: 62 62 T = rho * T / (re * cm) 63 63 64 chi = (pi / 2.0) - 2.0 * np y.arctan(T)64 chi = (pi / 2.0) - 2.0 * np.arctan(T) 65 65 lat = chi + ((ex2 / 2.0) + (5.0 * ex2**2.0 / 24.0) + (ex2**3.0 / 12.0)) * \ 66 np y.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \67 np y.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) 68 68 69 69 lat = sgn * lat 70 lon = np y.arctan2(sgn * x,-sgn * y)70 lon = np.arctan2(sgn * x,-sgn * y) 71 71 lon = sgn * lon 72 72 73 res1 = np y.nonzero(rho <= 0.1)[0]73 res1 = np.nonzero(rho <= 0.1)[0] 74 74 if len(res1) > 0: 75 75 lat[res1] = 90. * sgn -
issm/trunk-jpl/src/m/exp/expcoarsen.py
r17480 r21254 1 1 import os.path 2 import numpy as np y2 import numpy as np 3 3 from collections import OrderedDict 4 4 from expread import expread … … 34 34 for contour in contours: 35 35 36 numpoints=np y.size(contour['x'])36 numpoints=np.size(contour['x']) 37 37 38 38 j=0 … … 44 44 45 45 #see whether we keep this point or not 46 distance=np y.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) 47 47 if distance<resolution and j<numpoints-2: #do not remove last point 48 x=np y.delete(x,j+1,0)49 y=np y.delete(y,j+1,0)48 x=np.delete(x,j+1,0) 49 y=np.delete(y,j+1,0) 50 50 numpoints=numpoints-1 51 51 else: 52 division=int(np y.floor(distance/resolution)+1)52 division=int(np.floor(distance/resolution)+1) 53 53 if division>=2: 54 xi=np y.linspace(x[j],x[j+1],division)55 yi=np y.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) 56 56 57 x=np y.hstack((x[0:j+1],xi[1:-1],x[j+1:]))58 y=np y.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:])) 59 59 60 60 #update current point … … 65 65 j=j+1 66 66 67 if np y.size(x)>1:67 if np.size(x)>1: 68 68 #keep the (x,y) contour arond 69 69 newcontour=OrderedDict() 70 newcontour['nods']=np y.size(x)70 newcontour['nods']=np.size(x) 71 71 newcontour['density']=contour['density'] 72 72 newcontour['x']=x -
issm/trunk-jpl/src/m/exp/expdisp.py
r19421 r21254 1 1 from expread import expread 2 import numpy as np y2 import numpy as np 3 3 4 4 def expdisp(domainoutline,ax,linestyle='--k',linewidth=1,unitmultiplier=1.): -
issm/trunk-jpl/src/m/extrusion/DepthAverage.py
r17810 r21254 1 import numpy as np y1 import numpy as np 2 2 from project2d import project2d 3 3 … … 19 19 20 20 # coerce to array in case float is passed 21 if type(vector)!=np y.ndarray:21 if type(vector)!=np.ndarray: 22 22 print 'coercing array' 23 vector=np y.array(value)23 vector=np.array(value) 24 24 25 25 vec2d=False … … 30 30 #nods data 31 31 if vector.shape[0]==md.mesh.numberofvertices: 32 vector_average=np y.zeros(md.mesh.numberofvertices2d)32 vector_average=np.zeros(md.mesh.numberofvertices2d) 33 33 for i in xrange(1,md.mesh.numberoflayers): 34 34 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)) … … 37 37 #element data 38 38 elif vector.shape[0]==md.mesh.numberofelements: 39 vector_average=np y.zeros(md.mesh.numberofelements2d)39 vector_average=np.zeros(md.mesh.numberofelements2d) 40 40 for i in xrange(1,md.mesh.numberoflayers): 41 41 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 np y1 import numpy as np 2 2 3 3 def project2d(md3d,value,layer): … … 25 25 26 26 # coerce to array in case float is passed 27 if type(value)!=np y.ndarray:27 if type(value)!=np.ndarray: 28 28 print 'coercing array' 29 value=np y.array(value)29 value=np.array(value) 30 30 31 31 vec2d=False -
issm/trunk-jpl/src/m/geometry/slope.py
r17686 r21254 1 import numpy as np y1 import numpy as np 2 2 from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff 3 3 … … 33 33 alpha,beta=GetNodalFunctionsCoeff(index,x,y)[0:2] 34 34 35 summation=np y.array([[1],[1],[1]])36 sx=np y.dot(surf[index-1]*alpha,summation).reshape(-1,)37 sy=np y.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,) 38 38 39 s=np y.sqrt(sx**2+sy**2)39 s=np.sqrt(sx**2+sy**2) 40 40 41 41 if md.mesh.dimension()==3: 42 42 sx=project3d(md,'vector',sx,'type','element') 43 43 sy=project3d(md,'vector',sy,'type','element') 44 s=np y.sqrt(sx**2+sy**2)44 s=np.sqrt(sx**2+sy**2) 45 45 46 46 return (sx,sy,s) -
issm/trunk-jpl/src/m/interp/SectionValues.py
r20910 r21254 1 1 import os 2 2 from expread import expread 3 import numpy as np y3 import numpy as np 4 4 from project2d import project2d 5 5 #from InterpFromMesh2d import InterpFromMesh2d … … 41 41 42 42 #initialization 43 X=np y.array([]) #X-coordinate44 Y=np y.array([]) #Y-coordinate45 S=np y.array([0.]) #curvilinear coordinate43 X=np.array([]) #X-coordinate 44 Y=np.array([]) #Y-coordinate 45 S=np.array([0.]) #curvilinear coordinate 46 46 47 47 for i in xrange(nods-1): … … 53 53 s_start=S[-1] 54 54 55 length_segment=np y.sqrt((x_end-x_start)**2+(y_end-y_start)**2)56 portion=np y.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) 57 57 58 x_segment=np y.zeros(portion)59 y_segment=np y.zeros(portion)60 s_segment=np y.zeros(portion)58 x_segment=np.zeros(portion) 59 y_segment=np.zeros(portion) 60 s_segment=np.zeros(portion) 61 61 62 62 for j in xrange(int(portion)): … … 66 66 67 67 #plug into X and Y 68 X=np y.append(X,x_segment)69 Y=np y.append(Y,y_segment)70 S=np y.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) 71 71 72 X=np y.append(X,x[nods-1])73 Y=np y.append(Y,y[nods-1])72 X=np.append(X,x[nods-1]) 73 Y=np.append(Y,y[nods-1]) 74 74 75 75 #Number of nodes: … … 77 77 78 78 #Compute Z 79 Z=np y.zeros(numberofnodes)79 Z=np.zeros(numberofnodes) 80 80 81 81 #New mesh and Data interpolation … … 88 88 89 89 #Compute index 90 index=np y.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T90 index=np.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T 91 91 92 92 else: … … 102 102 103 103 #Some useful parameters 104 layers=int(np y.ceil(npy.mean(md.geometry.thickness)/res_v))104 layers=int(np.ceil(np.mean(md.geometry.thickness)/res_v)) 105 105 nodesperlayer=int(numberofnodes) 106 106 nodestot=int(nodesperlayer*layers) … … 109 109 110 110 #initialization 111 X3=np y.zeros(nodesperlayer*layers)112 Y3=np y.zeros(nodesperlayer*layers)113 Z3=np y.zeros(nodesperlayer*layers)114 S3=np y.zeros(nodesperlayer*layers)115 index3=np y.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)) 116 116 117 117 #Get new coordinates in 3d … … 123 123 124 124 if i<layers-1: #Build index3 with quads 125 ids=np y.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))).T125 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 126 126 index3[(i-1)*elementsperlayer:i*elementsperlayer,:]=ids 127 127 128 128 #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,np y.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] 130 130 131 131 #build outputs -
issm/trunk-jpl/src/m/interp/averaging.py
r20144 r21254 1 import numpy as np y1 import numpy as np 2 2 from GetAreas import GetAreas 3 3 import MatlabFuncs as m … … 39 39 #initialization 40 40 if layer==0: 41 weights=np y.zeros(md.mesh.numberofvertices,)41 weights=np.zeros(md.mesh.numberofvertices,) 42 42 data=data.flatten(1) 43 43 else: 44 weights=np y.zeros(md.mesh.numberofvertices2d,)44 weights=np.zeros(md.mesh.numberofvertices2d,) 45 45 data=data[(layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:] 46 46 … … 69 69 index=index-1 # since python indexes starting from zero 70 70 line=index.flatten(1) 71 areas=np y.vstack(areas).reshape(-1,)72 summation=1./rep*np y.ones(rep,)71 areas=np.vstack(areas).reshape(-1,) 72 summation=1./rep*np.ones(rep,) 73 73 linesize=rep*numberofelements 74 74 75 75 #update weights that holds the volume of all the element holding the node i 76 weights=csc_matrix( (np y.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)) 77 77 78 78 #initialization 79 79 if len(data)==numberofelements: 80 average_node=csc_matrix( (np y.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)) 81 81 average_node=average_node/weights 82 82 average_node = csc_matrix(average_node) … … 85 85 86 86 #loop over iteration 87 for i in np y.arange(1,iterations+1):88 average_el=np y.asarray(npy.dot(average_node.todense()[index].reshape(numberofelements,rep),npy.vstack(summation))).reshape(-1,)89 average_node=csc_matrix( (np y.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)) 90 90 average_node=average_node/weights 91 91 average_node=csc_matrix(average_node) 92 92 93 93 #return output as a full matrix (C code do not like sparse matrices) 94 average=np y.asarray(average_node.todense()).reshape(-1,)94 average=np.asarray(average_node.todense()).reshape(-1,) 95 95 96 96 return average -
issm/trunk-jpl/src/m/interp/holefiller.py
r17867 r21254 1 import numpy as np y1 import numpy as np 2 2 from scipy.spatial import cKDTree 3 3 … … 30 30 filled=data 31 31 32 XYGood=np y.dstack([x[goodids],y[goodids]])[0]33 XYBad=np y.dstack([x[badids],y[badids]])[0]32 XYGood=np.dstack([x[goodids],y[goodids]])[0] 33 XYBad=np.dstack([x[badids],y[badids]])[0] 34 34 tree=cKDTree(XYGood) 35 35 nearest=tree.query(XYBad,k=knn)[1] … … 42 42 for j in range(knn): 43 43 neardat.append(filled[goodids][nearest[i][j]]) 44 filled[badids[i]]=np y.mean(neardat)44 filled[badids[i]]=np.mean(neardat) 45 45 46 46 return filled -
issm/trunk-jpl/src/m/interp/interp.py
r18109 r21254 1 1 # module for inperpolating/smoothing data 2 import numpy as np y2 import numpy as np 3 3 from scipy.interpolate import CloughTocher2DInterpolator, Rbf 4 4 from scipy.spatial import cKDTree … … 51 51 xlim=[min(xi)-dx,max(xi)+dx] 52 52 ylim=[min(yi)-dy,max(yi)+dy] 53 xflag=np y.logical_and(x>xlim[0],x<xlim[1])54 yflag=np y.logical_and(y>ylim[0],y<ylim[1])55 bothind=np y.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)) 56 56 subdata=data[bothind] 57 57 subx=x[bothind] 58 58 suby=y[bothind] 59 points=np y.array([subx,suby]).T59 points=np.array([subx,suby]).T 60 60 61 61 # mask out any nan's in the data and corresponding coordinate points 62 mask=np y.isnan(subdata)63 ind=np y.nonzero(mask)[0]62 mask=np.isnan(subdata) 63 ind=np.nonzero(mask)[0] 64 64 if len(ind) and fill_nans: 65 65 print " WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully." 66 subdata=np y.delete(subdata,ind)67 points=np y.delete(points,ind,axis=0)66 subdata=np.delete(subdata,ind) 67 points=np.delete(points,ind,axis=0) 68 68 69 69 if maxiter: … … 76 76 if not fill_nans: 77 77 # identify nan's in xi,yi using nearest neighbors 78 xyinterp=np y.dstack([xi,yi])[0]79 xg,yg=np y.meshgrid(subx,suby)80 xydata=np y.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] 81 81 tree=cKDTree(xydata) 82 82 nearest=tree.query(xyinterp)[1] 83 pos=np y.nonzero(npy.isnan(subdata[nearest]))83 pos=np.nonzero(np.isnan(subdata[nearest])) 84 84 interpdata[pos]=subdata[nearest][pos] 85 85 86 86 return interpdata 87 87 #}}} 88 def GridSplineToMesh2d(x,y,data,xi,yi,default_value=np y.nan,plotonly=False,fill_nans=False):#{{{88 def GridSplineToMesh2d(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False):#{{{ 89 89 ''' 90 90 python analog to InterpFromGridToMesh. This routine uses … … 108 108 109 109 Usage: 110 interpdata=GridToMesh(x,y,data,xi,yi,default_value=np y.nan,plotonly=False,fill_nans=False)110 interpdata=GridToMesh(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False) 111 111 112 112 Examples: … … 114 114 ''' 115 115 116 if np y.ndim(x)==2:116 if np.ndim(x)==2: 117 117 x=x.reshape(-1,) 118 if np y.ndim(y)==2:118 if np.ndim(y)==2: 119 119 y=y.reshape(-1,) 120 120 if len(x) != data.shape[1]+1 and len(x) != data.shape[1]: … … 134 134 if len(x)==data.shape[1] and len(y)==data.shape[0]: 135 135 print ' x,y taken to define the center of data grid cells' 136 xind=np y.nonzero(npy.logical_and(x>xlim[0],x<xlim[1]))[0]137 yind=np y.nonzero(npy.logical_and(y>ylim[0],y<ylim[1]))[0]138 xg,yg=np y.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]) 139 139 subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1] 140 140 elif len(x)==data.shape[1]+1 and len(y)==data.shape[0]+1: 141 141 print ' x,y taken to define the corners of data grid cells' 142 xcenter=np y.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),npy.float)143 ycenter=np y.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),npy.float)144 xind=np y.nonzero(npy.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]145 yind=np y.nonzero(npy.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]146 xg,yg=np y.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]) 147 147 subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1] 148 148 else: 149 149 raise ValueError('x and y have inconsistent sizes: both should have length ncols(data)/nrows(data) or ncols(data)+1/nrows(data)+1') 150 150 151 points=np y.array([xg.ravel(),yg.ravel()]).T151 points=np.array([xg.ravel(),yg.ravel()]).T 152 152 flatsubdata=subdata.ravel() 153 153 154 154 if plotonly: 155 plt.imshow(np y.flipud(subdata),origin='upper')155 plt.imshow(np.flipud(subdata),origin='upper') 156 156 plt.show() 157 157 return 158 158 159 159 # mask out any nan's in the data and corresponding coordinate points 160 mask=np y.isnan(flatsubdata)161 ind=np y.nonzero(mask)[0]160 mask=np.isnan(flatsubdata) 161 ind=np.nonzero(mask)[0] 162 162 if len(ind) and fill_nans: 163 163 print " WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully." 164 goodsubdata=np y.delete(flatsubdata,ind)165 goodpoints=np y.delete(points,ind,axis=0)164 goodsubdata=np.delete(flatsubdata,ind) 165 goodpoints=np.delete(points,ind,axis=0) 166 166 167 167 # create spline and index spline at mesh points … … 171 171 if not fill_nans: 172 172 # identify nan's in xi,yi using nearest neighbors 173 xyinterp=np y.dstack([xi,yi])[0]174 xydata=np y.dstack([xg.ravel(),yg.ravel()])[0]173 xyinterp=np.dstack([xi,yi])[0] 174 xydata=np.dstack([xg.ravel(),yg.ravel()])[0] 175 175 tree=cKDTree(xydata) 176 176 nearest=tree.query(xyinterp)[1] 177 pos=np y.nonzero(npy.isnan(flatsubdata[nearest]))177 pos=np.nonzero(np.isnan(flatsubdata[nearest])) 178 178 interpdata[pos]=flatsubdata[nearest][pos] 179 179 -
issm/trunk-jpl/src/m/inversions/parametercontroldrag.py
r16108 r21254 26 26 27 27 #weights 28 weights=options.getfieldvalue('weights',np y.ones(md.mesh.numberofvertices))29 if np y.size(weights)!=md.mesh.numberofvertices:28 weights=options.getfieldvalue('weights',np.ones(md.mesh.numberofvertices)) 29 if np.size(weights)!=md.mesh.numberofvertices: 30 30 md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices) 31 31 else: … … 34 34 #nsteps 35 35 nsteps=options.getfieldvalue('nsteps',100); 36 if (np y.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):36 if (np.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps): 37 37 md.inversion.nsteps=100 38 38 else: … … 41 41 #cm_min 42 42 cm_min=options.getfieldvalue('cm_min',ones(md.mesh.numberofvertices)) 43 if (np y.size(cm_min)==1):43 if (np.size(cm_min)==1): 44 44 md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices) 45 elif (np y.size(cm_min)==md.mesh.numberofvertices):45 elif (np.size(cm_min)==md.mesh.numberofvertices): 46 46 md.inversion.min_parameters=cm_min 47 47 else: … … 50 50 #cm_max 51 51 cm_max=options.getfieldvalue('cm_max',250*ones(md.mesh.numberofvertices)) 52 if (np y.size(cm_max)==1):52 if (np.size(cm_max)==1): 53 53 md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices) 54 elif (np y.size(cm_max)==md.mesh.numberofvertices):54 elif (np.size(cm_max)==md.mesh.numberofvertices): 55 55 md.inversion.max_parameters=cm_max 56 56 else: … … 59 59 #eps_cm 60 60 eps_cm=optoins.getfieldvalue('eps_cm',float('nan')) 61 if (np y.size(eps_cm)~=1 | eps_cm<0 ):61 if (np.size(eps_cm)~=1 | eps_cm<0 ): 62 62 md.inversion.cost_function_threshold=float('nan') 63 63 else: … … 66 66 #maxiter 67 67 maxiter=options.getfieldvalue('maxiter',10*ones(md.inversion.nsteps)) 68 if (np y.any(maxiter<0) | npy.any(floor(maxiter)~=maxiter)):68 if (np.any(maxiter<0) | np.any(floor(maxiter)~=maxiter)): 69 69 md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps) 70 70 else: … … 75 75 #cm_jump 76 76 cm_jump=options.getfieldvalue('cm_jump',0.8*ones(md.inversion.nsteps)) 77 if !np y.isreal(cm_jump):77 if !np.isreal(cm_jump): 78 78 md.inversion.step_threshold=0.8*ones(md.inversion.nsteps) 79 79 else: -
issm/trunk-jpl/src/m/materials/TMeltingPoint.py
r17196 r21254 1 import numpy as np y1 import numpy as np 2 2 3 3 def TMeltingPoint(reftemp,pressure): … … 17 17 18 18 #ensure ref is same dimension as pressure 19 ref=reftemp*np y.ones_like(pressure)19 ref=reftemp*np.ones_like(pressure) 20 20 21 21 return reftemp-beta*pressure -
issm/trunk-jpl/src/m/mech/analyticaldamage.py
r17965 r21254 1 import numpy as np y1 import numpy as np 2 2 from averaging import averaging 3 3 #from plotmodel import plotmodel … … 54 54 55 55 if isinstance(sigmab,(int,float)): 56 sigmab=sigmab*np y.ones((md.mesh.numberofvertices,))56 sigmab=sigmab*np.ones((md.mesh.numberofvertices,)) 57 57 58 58 # check inputs … … 61 61 if not '2d' in md.mesh.__doc__: 62 62 raise StandardError('only 2d (planview) model supported currently') 63 if np y.any(md.flowequation.element_equation!=2):63 if np.any(md.flowequation.element_equation!=2): 64 64 print 'Warning: the model has some non SSA elements. These will be treated like SSA elements' 65 65 … … 76 76 n=averaging(md,md.materials.rheology_n,0) 77 77 78 D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/np y.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) 79 79 80 80 # D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed 81 pos=np y.nonzero(D>1)81 pos=np.nonzero(D>1) 82 82 D[pos]=0 83 83 84 backstress=np y.zeros((md.mesh.numberofvertices,))84 backstress=np.zeros((md.mesh.numberofvertices,)) 85 85 86 86 # backstress to bring D down to one 87 backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np y.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]) 88 88 89 pos=np y.nonzero(D<0)89 pos=np.nonzero(D<0) 90 90 #mask=ismember(1:md.mesh.numberofvertices,pos); 91 91 D[pos]=0 92 92 93 93 # backstress to bring negative damage to zero 94 backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np y.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]) 95 95 96 pos=np y.nonzero(backstress<0)96 pos=np.nonzero(backstress<0) 97 97 backstress[pos]=0 98 98 99 99 # rigidity from Thomas relation for D=0 and backstress=0 100 B=np y.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(npy.abs(ex)**(1./n))101 pos=np y.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) 102 102 B[pos]=md.materials.rheology_B[pos] 103 103 -
issm/trunk-jpl/src/m/mech/backstressfrominversion.py
r18011 r21254 1 import numpy as np y1 import numpy as np 2 2 from averaging import averaging 3 3 from thomasparams import thomasparams … … 65 65 if tempmask: 66 66 Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar 67 pos=np y.nonzero(Bi>md.materials.rheology_B)67 pos=np.nonzero(Bi>md.materials.rheology_B) 68 68 Bi[pos]=md.materials.rheology_B[pos] 69 69 70 70 # analytical backstress solution 71 backstress=T-Bi*np y.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))72 backstress[np y.nonzero(backstress<0)]=071 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 73 73 74 74 return backstress -
issm/trunk-jpl/src/m/mech/calcbackstress.py
r18011 r21254 1 import numpy as np y1 import numpy as np 2 2 from averaging import averaging 3 3 from thomasparams import thomasparams … … 61 61 62 62 # analytical backstress solution 63 backstress=T-(1.-D)*B*np y.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))64 backstress[np y.nonzero(backstress<0)]=063 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 65 65 66 66 return backstress -
issm/trunk-jpl/src/m/mech/damagefrominversion.py
r17624 r21254 1 import numpy as np y1 import numpy as np 2 2 3 3 def damagefrominversion(md): … … 24 24 if any(md.flowequation.element_equation!=2): 25 25 raise StandardError('Warning: the model has some non-SSA elements. These will be treated like SSA elements') 26 if np y.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:26 if np.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2: 27 27 Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1,) 28 28 else: 29 29 Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar 30 if np y.ndim(md.materials.rheology_B)==2:30 if np.ndim(md.materials.rheology_B)==2: 31 31 BT=md.materials.rheology_B.reshape(-1,) 32 32 else: 33 33 BT=md.materials.rheology_B 34 34 35 damage=np y.zeros_like(Bi)35 damage=np.zeros_like(Bi) 36 36 37 37 # Damage where Bi softer than B(T) 38 pos=np y.nonzero(Bi<BT)[0]38 pos=np.nonzero(Bi<BT)[0] 39 39 damage[pos]=1.-Bi[pos]/BT[pos] 40 40 41 pos=np y.nonzero(damage<0)41 pos=np.nonzero(damage<0) 42 42 damage[pos]=0 43 43 -
issm/trunk-jpl/src/m/mech/mechanicalproperties.py
r19466 r21254 1 import numpy as np y1 import numpy as np 2 2 from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff 3 3 from results import results … … 27 27 # raise StandardError('only 2D model supported currently') 28 28 29 if np y.any(md.flowequation.element_equation!=2):29 if np.any(md.flowequation.element_equation!=2): 30 30 print 'Warning: the model has some non SSA elements. These will be treated like SSA elements' 31 31 … … 35 35 if len(damage)!=md.mesh.numberofvertices: 36 36 raise ValueError('if damage is supplied it should be of size ' + md.mesh.numberofvertices) 37 if np y.ndim(damage)==2:37 if np.ndim(damage)==2: 38 38 damage=damage.reshape(-1,) 39 39 else: damage=None 40 40 41 if np y.ndim(vx)==2:41 if np.ndim(vx)==2: 42 42 vx=vx.reshape(-1,) 43 if np y.ndim(vy)==2:43 if np.ndim(vy)==2: 44 44 vy=vy.reshape(-1,) 45 45 … … 48 48 numberofvertices=md.mesh.numberofvertices 49 49 index=md.mesh.elements 50 summation=np y.array([[1],[1],[1]])51 directionsstress=np y.zeros((numberofelements,4))52 directionsstrain=np y.zeros((numberofelements,4))53 valuesstress=np y.zeros((numberofelements,2))54 valuesstrain=np y.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)) 55 55 56 56 #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma … … 60 60 vxlist=vx[index-1]/md.constants.yts 61 61 vylist=vy[index-1]/md.constants.yts 62 ux=np y.dot((vxlist*alpha),summation).reshape(-1,)63 uy=np y.dot((vxlist*beta),summation).reshape(-1,)64 vx=np y.dot((vylist*alpha),summation).reshape(-1,)65 vy=np y.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,) 66 66 uyvx=(vx+uy)/2. 67 67 #clear vxlist vylist 68 68 69 69 #compute viscosity 70 nu=np y.zeros((numberofelements,))71 B_bar=np y.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,) 72 72 power=((md.materials.rheology_n-1.)/(2.*md.materials.rheology_n)).reshape(-1,) 73 73 second_inv=(ux**2.+vy**2.+((uy+vx)**2.)/4.+ux*vy).reshape(-1,) 74 74 75 75 #some corrections 76 location=np y.nonzero(npy.logical_and(second_inv==0,power!=0))76 location=np.nonzero(np.logical_and(second_inv==0,power!=0)) 77 77 nu[location]=10^18 #arbitrary maximum viscosity to apply where there is no effective shear 78 78 79 79 if 'matice' in md.materials.__module__: 80 location=np y.nonzero(second_inv)80 location=np.nonzero(second_inv) 81 81 nu[location]=B_bar[location]/(second_inv[location]**power[location]) 82 location=np y.nonzero(npy.logical_and(second_inv==0,power==0))82 location=np.nonzero(np.logical_and(second_inv==0,power==0)) 83 83 nu[location]=B_bar[location] 84 location=np y.nonzero(npy.logical_and(second_inv==0,power!=0))84 location=np.nonzero(np.logical_and(second_inv==0,power!=0)) 85 85 nu[location]=10^18 86 86 elif 'matdamageice' in md.materials.__module__ and damage is not None: 87 87 print 'computing damage-dependent properties!' 88 Zinv=np y.dot(1-damage[index-1],summation/3.).reshape(-1,)89 location=np y.nonzero(second_inv)90 nu[location]=Zinv[location]*B_bar[location]/np y.power(second_inv[location],power[location])91 location=np y.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)) 92 92 nu[location]=Zinv[location]*B_bar[location] 93 93 #clear Zinv … … 101 101 102 102 #compute principal properties of stress 103 for i in np y.arange(numberofelements):103 for i in np.arange(numberofelements): 104 104 105 105 #compute stress and strainrate matrices 106 stress=np y.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])107 strain=np y.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]] ]) 108 108 109 109 #eigenvalues and vectors for stress 110 value,directions=np y.linalg.eig(stress);110 value,directions=np.linalg.eig(stress); 111 111 idx=abs(value).argsort()[::-1] # sort in descending order 112 112 value=value[idx] … … 116 116 117 117 #eigenvalues and vectors for strain 118 value,directions=np y.linalg.eig(strain);118 value,directions=np.linalg.eig(strain); 119 119 idx=abs(value).argsort()[::-1] # sort in descending order 120 120 value=value[idx] … … 133 133 stress.principalvalue2=valuesstress[:,1] 134 134 stress.principalaxis2=directionsstress[:,2:4] 135 stress.effectivevalue=1./np y.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) 136 136 md.results.stress=stress 137 137 … … 144 144 strainrate.principalvalue2=valuesstrain[:,1]*md.constants.yts 145 145 strainrate.principalaxis2=directionsstrain[:,2:4] 146 strainrate.effectivevalue=1./np y.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) 147 147 md.results.strainrate=strainrate 148 148 … … 155 155 deviatoricstress.principalvalue2=valuesstress[:,1] 156 156 deviatoricstress.principalaxis2=directionsstress[:,2:4] 157 deviatoricstress.effectivevalue=1./np y.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) 158 158 md.results.deviatoricstress=deviatoricstress 159 159 -
issm/trunk-jpl/src/m/mech/robintemperature.py
r17836 r21254 1 import numpy as np y1 import numpy as np 2 2 from scipy.special import erf 3 3 … … 34 34 35 35 #create vertical coordinate variable 36 zstar=np y.sqrt(2.*alphaT*thickness/accumrate)36 zstar=np.sqrt(2.*alphaT*thickness/accumrate) 37 37 38 tprofile=surftemp+np y.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)) 39 39 40 40 return tprofile -
issm/trunk-jpl/src/m/mech/steadystateiceshelftemp.py
r18068 r21254 1 import numpy as np y1 import numpy as np 2 2 3 3 def steadystateiceshelftemp(md,surfacetemp,basaltemp): … … 42 42 temperature=(Ts+Tb)/2 # where wi~=0 43 43 44 pos=np y.nonzero(abs(wi)>=1e-4) # to avoid division by zero44 pos=np.nonzero(abs(wi)>=1e-4) # to avoid division by zero 45 45 46 np y.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation46 np.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation 47 47 #calculate depth-averaged temperature (in Celsius) 48 48 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])*np y.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)) 50 50 except FloatingPointError: 51 51 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])/np y.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))) 53 53 54 54 #temperature should not be less than surface temp 55 pos=np y.nonzero(temperature<Ts)55 pos=np.nonzero(temperature<Ts) 56 56 temperature[pos]=Ts[pos] 57 57 58 58 # NaN where melt rates are too high (infinity/infinity in exponential) 59 pos=np y.nonzero(npy.isnan(temperature))59 pos=np.nonzero(np.isnan(temperature)) 60 60 temperature[pos]=Ts[pos] 61 61 -
issm/trunk-jpl/src/m/mech/thomasparams.py
r17903 r21254 1 import numpy as np y1 import numpy as np 2 2 from averaging import averaging 3 3 … … 75 75 76 76 # checks: any of e1 or e2 equal to zero? 77 pos=np y.nonzero(e1==0)78 if np y.any(pos==1):77 pos=np.nonzero(e1==0) 78 if np.any(pos==1): 79 79 print 'WARNING: first principal strain rate equal to zero. Value set to 1e-13 s^-1' 80 80 e1[pos]=1.e-13 81 pos=np y.nonzero(e2==0)82 if np y.any(pos==1):81 pos=np.nonzero(e2==0) 82 if np.any(pos==1): 83 83 print 'WARNING: second principal strain rate equal to zero. Value set to 1e-13 s^-1' 84 84 e2[pos]=1.e-13 … … 89 89 90 90 if coordsys=='principal': 91 b=np y.zeros((md.mesh.numberofvertices,))91 b=np.zeros((md.mesh.numberofvertices,)) 92 92 ex=e1 93 93 a=e2/e1 94 pos=np y.nonzero(npy.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension94 pos=np.nonzero(np.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension 95 95 a[pos]=e1[pos]/e2[pos] 96 96 ex[pos]=e2[pos] 97 pos2=np y.nonzero(e1<0 & e2<0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal compression97 pos2=np.nonzero(e1<0 & e2<0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal compression 98 98 a[pos2]=e1[pos2]/e2[pos2] 99 99 ex[pos2]=e2[pos2] 100 pos3=np y.nonzero(e1>0 & e2>0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal tension100 pos3=np.nonzero(e1>0 & e2>0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal tension 101 101 a[pos3]=e1[pos3]/e2[pos3] 102 102 ex[pos3]=e2[pos3] 103 ind=np y.nonzero(e1<0 & e2<0)103 ind=np.nonzero(e1<0 & e2<0) 104 104 a[ind]=-a[ind] # where both strain rates are compressive, enforce negative alpha 105 sigxx=(np y.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B105 sigxx=(np.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B 106 106 elif coordsys=='xy': 107 107 ex=exx … … 110 110 elif coordsys=='longitudinal': 111 111 # using longitudinal strain rates defined by observed velocity vector 112 velangle=np y.arctan(md.initialization.vy/md.initialization.vx)113 pos=np y.nonzero(md.initialization.vx==0)114 velangle[pos]=np y.pi/2115 ex=0.5*(exx+eyy)+0.5*(exx-eyy)*np y.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) 116 116 ey=exx+eyy-ex # trace of strain rate tensor is invariant 117 exy=-0.5*(exx-eyy)*np y.sin(2.*velangle)+exy*npy.cos(2.*velangle)117 exy=-0.5*(exx-eyy)*np.sin(2.*velangle)+exy*np.cos(2.*velangle) 118 118 a=ey/ex 119 119 b=exy/ex … … 124 124 # a < -1 in areas of strong lateral compression or longitudinal compression and 125 125 # theta flips sign at a = -2 126 pos=np y.nonzero(npy.abs((npy.abs(a)-2.))<1.e-3)126 pos=np.nonzero(np.abs((np.abs(a)-2.))<1.e-3) 127 127 if len(pos)>0: 128 128 print 'Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2' … … 131 131 if eq=='Weertman1D': 132 132 theta=1./8 133 a=np y.zeros((md.mesh.numberofvertices,))133 a=np.zeros((md.mesh.numberofvertices,)) 134 134 elif eq=='Weertman2D': 135 135 theta=1./9 136 a=np y.ones((md.mesh.numberofvertices,))136 a=np.ones((md.mesh.numberofvertices,)) 137 137 elif eq=='Thomas': 138 theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(np y.abs(2.+a)**n)138 theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(np.abs(2.+a)**n) 139 139 else: 140 140 raise ValueError('argument passed to "eq" not valid') -
issm/trunk-jpl/src/m/plot/applyoptions.py
r21253 r21254 1 import numpy as np y1 import numpy as np 2 2 from cmaptools import truncate_colormap 3 3 from plot_contour import plot_contour … … 223 223 cb.locator=MaxNLocator(nbins=5) # default 5 ticks 224 224 if options.exist('colorbartickspacing'): 225 locs=np y.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))225 locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing')) 226 226 cb.set_ticks(locs) 227 227 if options.exist('colorbarlines'): 228 locs=np y.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))229 cb.add_lines(locs,['k' for i in range(len(locs))],np y.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)) 230 230 if options.exist('colorbarlineatvalue'): 231 231 locs=options.getfieldvalue('colorbarlineatvalue') 232 232 colors=options.getfieldvalue('colorbarlineatvaluecolor',['k' for i in range (len(locs))]) 233 widths=options.getfieldvalue('colorbarlineatvaluewidth',np y.ones_like(locs))233 widths=options.getfieldvalue('colorbarlineatvaluewidth',np.ones_like(locs)) 234 234 cb.add_lines(locs,colors,widths) 235 235 if options.exist('colorbartitle'): -
issm/trunk-jpl/src/m/plot/colormaps/cmaptools.py
r17649 r21254 1 import numpy as np y1 import numpy as np 2 2 3 3 try: … … 21 21 22 22 new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, 23 a=minval, b=maxval), cmap(np y.linspace(minval, maxval, n)))23 a=minval, b=maxval), cmap(np.linspace(minval, maxval, n))) 24 24 25 25 return new_cmap -
issm/trunk-jpl/src/m/plot/plot_contour.py
r20921 r21254 26 26 pass 27 27 elif datatype==3: # quiver (vector) data 28 data=np y.sqrt(datain**2)28 data=np.sqrt(datain**2) 29 29 else: 30 30 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 np y1 import numpy as np 2 2 from processmesh import processmesh 3 3 from processdata import processdata … … 28 28 if data=='none' or data==None: 29 29 imageonly=1 30 data=np y.float('nan')*npy.ones((md.mesh.numberofvertices,))30 data=np.float('nan')*np.ones((md.mesh.numberofvertices,)) 31 31 datatype=1 32 32 else: … … 70 70 71 71 # normalize array 72 arr=arr/np y.float(npy.max(arr.ravel()))72 arr=arr/np.float(np.max(arr.ravel())) 73 73 arr=1.-arr # somehow the values got flipped 74 74 … … 96 96 dy=trans[5] 97 97 98 xarr=np y.arange(xmin,xmax,dx)99 yarr=np y.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)100 xg,yg=np y.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) 101 101 overlaylims=options.getfieldvalue('overlaylims',[min(arr.ravel()),max(arr.ravel())]) 102 102 norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1]) … … 118 118 #test 119 119 #m.ax=ax 120 meridians=np y.arange(-180.,181.,1.)121 parallels=np y.arange(-80.,80.,1.)120 meridians=np.arange(-180.,181.,1.) 121 parallels=np.arange(-80.,80.,1.) 122 122 m.drawparallels(parallels,labels=[0,0,1,1]) # labels=[left,right,top,bottom] 123 123 m.drawmeridians(meridians,labels=[1,1,0,0]) 124 124 m.drawcoastlines() 125 pc=m.pcolormesh(xg, yg, np y.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) 126 126 127 127 else: 128 pc=ax.pcolormesh(xg, yg, np y.flipud(arr), cmap=mpl.cm.Greys, norm=norm)128 pc=ax.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm) 129 129 130 130 #rasterization? -
issm/trunk-jpl/src/m/plot/plot_streamlines.py
r20144 r21254 1 import numpy as np y1 import numpy as np 2 2 from processmesh import processmesh 3 3 from processdata import processdata … … 45 45 46 46 # format data for matplotlib streamplot function 47 yg,xg=np y.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] 48 48 ug=griddata((x,y),u,(xg,yg),method='linear') 49 49 vg=griddata((x,y),v,(xg,yg),method='linear') … … 60 60 if linewidth=='vel': 61 61 scale=options.getfieldvalue('streamlineswidthscale',3) 62 vel=np y.sqrt(ug**2+vg**2)63 linewidth=scale*vel/np y.amax(vel)62 vel=np.sqrt(ug**2+vg**2) 63 linewidth=scale*vel/np.amax(vel) 64 64 linewidth[linewidth<0.5]=0.5 65 65 -
issm/trunk-jpl/src/m/plot/plotmodel.py
r21253 r21254 1 import numpy as np y1 import numpy as np 2 2 from plotoptions import plotoptions 3 3 … … 34 34 nr=True 35 35 else: 36 nrows=np y.ceil(numberofplots/subplotwidth)36 nrows=np.ceil(numberofplots/subplotwidth) 37 37 nr=False 38 38 -
issm/trunk-jpl/src/m/plot/processmesh.py
r21253 r21254 1 1 from math import isnan 2 import numpy as np y2 import numpy as np 3 3 4 4 def processmesh(md,data,options): … … 35 35 z=md.mesh.z 36 36 except AttributeError: 37 z=np y.zeros_like(md.mesh.x)37 z=np.zeros_like(md.mesh.x) 38 38 39 39 try: -
issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
r20900 r21254 3 3 from collections import OrderedDict 4 4 import results as resultsclass 5 import MatlabFuncs as m6 5 7 6 def parseresultsfromdisk(md,filename,iosplit): 8 7 if iosplit: 9 results=parseresultsfromdiskiosplit(md,filename)8 saveres=parseresultsfromdiskiosplit(md,filename) 10 9 else: 11 results=parseresultsfromdiskioserial(md,filename)12 13 return results10 saveres=parseresultsfromdiskioserial(md,filename) 11 12 return saveres 14 13 15 14 def parseresultsfromdiskioserial(md,filename): # {{{ 16 17 15 #Open file 18 16 try: … … 22 20 23 21 #initialize results: 24 results=[] 25 results.append(None) 22 saveres=[] 26 23 27 24 #Read fields until the end of the file. 28 result=ReadData(fid,md)25 loadres=ReadData(fid,md) 29 26 30 27 counter=0 31 28 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: 36 34 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: 39 36 raise TypeError("parsing results for a steady-state core, which incorporates transient results!") 40 37 41 38 #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): 43 40 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: 48 45 #if we have a step = 0, this is a steady state solution, don't expect more steps. 49 46 index = 0; 50 47 check_nomoresteps=1 51 52 elif result['step']==1: 48 elif loadres['step']==1: 53 49 index = 0 54 50 else: 55 51 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())61 52 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() 65 59 66 60 #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) 80 71 81 72 fid.close() 82 73 83 return results74 return saveres 84 75 # }}} 85 76 def parseresultsfromdiskiosplit(md,filename): # {{{ … … 91 82 raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename) 92 83 93 results=[]84 saveres=[] 94 85 95 86 #if we have done split I/O, ie, we have results that are fragmented across patches, 96 87 #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: 99 90 100 91 #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) 113 104 114 105 #do a second pass, and figure out the size of the patches 115 106 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) 121 112 122 113 #third pass, this time to read the real information 123 114 fid.seek(0) #rewind 124 result=ReadData(fid,md)125 while result:115 loadres=ReadData(fid,md) 116 while loadres: 126 117 127 118 #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) 140 131 141 132 #close file 142 133 fid.close() 143 134 144 return results135 return saveres 145 136 # }}} 146 137 def ReadData(fid,md): # {{{ … … 177 168 #Process units here FIXME: this should not be done here! 178 169 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': 182 173 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': 198 189 field = field/10.**12.*yts #(GigaTon/year) 199 elif m.strcmp(fieldname,'TotalGroundedBmb'):190 elif fieldname=='TotalGroundedBmb': 200 191 field = field/10.**12.*yts #(GigaTon/year) 201 elif m.strcmp(fieldname,'TotalSmb'):192 elif fieldname=='TotalSmb': 202 193 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']=fieldname210 result['time']=time211 result['step']=step212 result['field']=field194 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 213 204 214 205 except struct.error as e: 215 result=None216 217 return result206 saveres=None 207 208 return saveres 218 209 # }}} 219 210 def ReadDataDimensions(fid): # {{{ … … 246 237 raise TypeError("cannot read data of type %d" % type) 247 238 248 result=OrderedDict()249 result['fieldname']=fieldname250 result['time']=time251 result['step']=step252 result['M']=M253 result['N']=N239 saveres=OrderedDict() 240 saveres['fieldname']=fieldname 241 saveres['time']=time 242 saveres['step']=step 243 saveres['M']=M 244 saveres['N']=N 254 245 255 246 except struct.error as e: 256 result=None257 258 return result259 # }}} 247 saveres=None 248 249 return saveres 250 # }}}
Note:
See TracChangeset
for help on using the changeset viewer.